1
/*M///////////////////////////////////////////////////////////////////////////////////////
3
// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
5
// By downloading, copying, installing or using the software you agree to this license.
6
// If you do not agree to this license, do not download, install,
7
// copy or use the software.
10
// Intel License Agreement
11
// For Open Source Computer Vision Library
13
// Copyright (C) 2000, Intel Corporation, all rights reserved.
14
// Third party copyrights are property of their respective owners.
16
// Redistribution and use in source and binary forms, with or without modification,
17
// are permitted provided that the following conditions are met:
19
// * Redistribution's of source code must retain the above copyright notice,
20
// this list of conditions and the following disclaimer.
22
// * Redistribution's in binary form must reproduce the above copyright notice,
23
// this list of conditions and the following disclaimer in the documentation
24
// and/or other materials provided with the distribution.
26
// * The name of Intel Corporation may not be used to endorse or promote products
27
// derived from this software without specific prior written permission.
29
// This software is provided by the copyright holders and contributors "as is" and
30
// any express or implied warranties, including, but not limited to, the implied
31
// warranties of merchantability and fitness for a particular purpose are disclaimed.
32
// In no event shall the Intel Corporation or contributors be liable for any direct,
33
// indirect, incidental, special, exemplary, or consequential damages
34
// (including, but not limited to, procurement of substitute goods or services;
35
// loss of use, data, or profits; or business interruption) however caused
36
// and on any theory of liability, whether in contract, strict liability,
37
// or tort (including negligence or otherwise) arising in any way out of
38
// the use of this software, even if advised of the possibility of such damage.
42
/****************************************************************************************\
44
Calculation of a texture descriptors from GLCM (Grey Level Co-occurrence Matrix'es)
45
The code was submitted by Daniel Eaton [danieljameseaton@yahoo.com]
47
\****************************************************************************************/
49
#include "precomp.hpp"
54
#define CV_MAX_NUM_GREY_LEVELS_8U 256
62
int numLookupTableElements;
63
int forwardLookupTable[CV_MAX_NUM_GREY_LEVELS_8U];
64
int reverseLookupTable[CV_MAX_NUM_GREY_LEVELS_8U];
68
int descriptorOptimizationType;
73
static void icvCreateGLCM_LookupTable_8u_C1R( const uchar* srcImageData, int srcImageStep,
74
CvSize srcImageSize, CvGLCM* destGLCM,
75
int* steps, int numSteps, int* memorySteps );
78
icvCreateGLCMDescriptors_AllowDoubleNest( CvGLCM* destGLCM, int matrixIndex );
82
cvCreateGLCM( const IplImage* srcImage,
84
const int* srcStepDirections,/* should be static array..
85
or if not the user should handle de-allocation */
86
int numStepDirections,
87
int optimizationType )
89
static const int defaultStepDirections[] = { 0,1, -1,1, -1,0, -1,-1 };
93
int* stepDirections = 0;
95
CV_FUNCNAME( "cvCreateGLCM" );
99
uchar* srcImageData = 0;
103
const int maxNumGreyLevels8u = CV_MAX_NUM_GREY_LEVELS_8U;
106
CV_ERROR( CV_StsNullPtr, "" );
108
if( srcImage->nChannels != 1 )
109
CV_ERROR( CV_BadNumChannels, "Number of channels must be 1");
111
if( srcImage->depth != IPL_DEPTH_8U )
112
CV_ERROR( CV_BadDepth, "Depth must be equal IPL_DEPTH_8U");
114
// no Directions provided, use the default ones - 0 deg, 45, 90, 135
115
if( !srcStepDirections )
117
srcStepDirections = defaultStepDirections;
120
CV_CALL( stepDirections = (int*)cvAlloc( numStepDirections*2*sizeof(stepDirections[0])));
121
memcpy( stepDirections, srcStepDirections, numStepDirections*2*sizeof(stepDirections[0]));
123
cvGetImageRawData( srcImage, &srcImageData, &srcImageStep, &srcImageSize );
125
// roll together Directions and magnitudes together with knowledge of image (step)
126
CV_CALL( memorySteps = (int*)cvAlloc( numStepDirections*sizeof(memorySteps[0])));
128
for( stepLoop = 0; stepLoop < numStepDirections; stepLoop++ )
130
stepDirections[stepLoop*2 + 0] *= stepMagnitude;
131
stepDirections[stepLoop*2 + 1] *= stepMagnitude;
133
memorySteps[stepLoop] = stepDirections[stepLoop*2 + 0]*srcImageStep +
134
stepDirections[stepLoop*2 + 1];
137
CV_CALL( newGLCM = (CvGLCM*)cvAlloc(sizeof(newGLCM)));
138
memset( newGLCM, 0, sizeof(*newGLCM) );
140
newGLCM->matrices = 0;
141
newGLCM->numMatrices = numStepDirections;
142
newGLCM->optimizationType = optimizationType;
144
if( optimizationType <= CV_GLCM_OPTIMIZATION_LUT )
146
int lookupTableLoop, imageColLoop, imageRowLoop, lineOffset = 0;
148
// if optimization type is set to lut, then make one for the image
149
if( optimizationType == CV_GLCM_OPTIMIZATION_LUT )
151
for( imageRowLoop = 0; imageRowLoop < srcImageSize.height;
152
imageRowLoop++, lineOffset += srcImageStep )
154
for( imageColLoop = 0; imageColLoop < srcImageSize.width; imageColLoop++ )
156
newGLCM->forwardLookupTable[srcImageData[lineOffset+imageColLoop]]=1;
160
newGLCM->numLookupTableElements = 0;
162
for( lookupTableLoop = 0; lookupTableLoop < maxNumGreyLevels8u; lookupTableLoop++ )
164
if( newGLCM->forwardLookupTable[ lookupTableLoop ] != 0 )
166
newGLCM->forwardLookupTable[ lookupTableLoop ] =
167
newGLCM->numLookupTableElements;
168
newGLCM->reverseLookupTable[ newGLCM->numLookupTableElements ] =
171
newGLCM->numLookupTableElements++;
175
// otherwise make a "LUT" which contains all the gray-levels (for code-reuse)
176
else if( optimizationType == CV_GLCM_OPTIMIZATION_NONE )
178
for( lookupTableLoop = 0; lookupTableLoop <maxNumGreyLevels8u; lookupTableLoop++ )
180
newGLCM->forwardLookupTable[ lookupTableLoop ] = lookupTableLoop;
181
newGLCM->reverseLookupTable[ lookupTableLoop ] = lookupTableLoop;
183
newGLCM->numLookupTableElements = maxNumGreyLevels8u;
186
newGLCM->matrixSideLength = newGLCM->numLookupTableElements;
187
icvCreateGLCM_LookupTable_8u_C1R( srcImageData, srcImageStep, srcImageSize,
188
newGLCM, stepDirections,
189
numStepDirections, memorySteps );
191
else if( optimizationType == CV_GLCM_OPTIMIZATION_HISTOGRAM )
193
CV_ERROR( CV_StsBadFlag, "Histogram-based method is not implemented" );
195
/* newGLCM->numMatrices *= 2;
196
newGLCM->matrixSideLength = maxNumGreyLevels8u*2;
198
icvCreateGLCM_Histogram_8uC1R( srcImageStep, srcImageSize, srcImageData,
199
newGLCM, numStepDirections,
200
stepDirections, memorySteps );
206
cvFree( &memorySteps );
207
cvFree( &stepDirections );
209
if( cvGetErrStatus() < 0 )
219
cvReleaseGLCM( CvGLCM** GLCM, int flag )
221
CV_FUNCNAME( "cvReleaseGLCM" );
228
CV_ERROR( CV_StsNullPtr, "" );
231
EXIT; // repeated deallocation: just skip it.
233
if( (flag == CV_GLCM_GLCM || flag == CV_GLCM_ALL) && (*GLCM)->matrices )
235
for( matrixLoop = 0; matrixLoop < (*GLCM)->numMatrices; matrixLoop++ )
237
if( (*GLCM)->matrices[ matrixLoop ] )
239
cvFree( (*GLCM)->matrices[matrixLoop] );
240
cvFree( (*GLCM)->matrices + matrixLoop );
244
cvFree( &((*GLCM)->matrices) );
247
if( (flag == CV_GLCM_DESC || flag == CV_GLCM_ALL) && (*GLCM)->descriptors )
249
for( matrixLoop = 0; matrixLoop < (*GLCM)->numMatrices; matrixLoop++ )
251
cvFree( (*GLCM)->descriptors + matrixLoop );
253
cvFree( &((*GLCM)->descriptors) );
256
if( flag == CV_GLCM_ALL )
266
icvCreateGLCM_LookupTable_8u_C1R( const uchar* srcImageData,
274
int* stepIncrementsCounter = 0;
276
CV_FUNCNAME( "icvCreateGLCM_LookupTable_8u_C1R" );
280
int matrixSideLength = destGLCM->matrixSideLength;
281
int stepLoop, sideLoop1, sideLoop2;
282
int colLoop, rowLoop, lineOffset = 0;
283
double*** matrices = 0;
285
// allocate memory to the matrices
286
CV_CALL( destGLCM->matrices = (double***)cvAlloc( sizeof(matrices[0])*numSteps ));
287
matrices = destGLCM->matrices;
289
for( stepLoop=0; stepLoop<numSteps; stepLoop++ )
291
CV_CALL( matrices[stepLoop] = (double**)cvAlloc( sizeof(matrices[0])*matrixSideLength ));
292
CV_CALL( matrices[stepLoop][0] = (double*)cvAlloc( sizeof(matrices[0][0])*
293
matrixSideLength*matrixSideLength ));
295
memset( matrices[stepLoop][0], 0, matrixSideLength*matrixSideLength*
296
sizeof(matrices[0][0]) );
298
for( sideLoop1 = 1; sideLoop1 < matrixSideLength; sideLoop1++ )
300
matrices[stepLoop][sideLoop1] = matrices[stepLoop][sideLoop1-1] + matrixSideLength;
304
CV_CALL( stepIncrementsCounter = (int*)cvAlloc( numSteps*sizeof(stepIncrementsCounter[0])));
305
memset( stepIncrementsCounter, 0, numSteps*sizeof(stepIncrementsCounter[0]) );
307
// generate GLCM for each step
308
for( rowLoop=0; rowLoop<srcImageSize.height; rowLoop++, lineOffset+=srcImageStep )
310
for( colLoop=0; colLoop<srcImageSize.width; colLoop++ )
312
int pixelValue1 = destGLCM->forwardLookupTable[srcImageData[lineOffset + colLoop]];
314
for( stepLoop=0; stepLoop<numSteps; stepLoop++ )
317
row2 = rowLoop + steps[stepLoop*2 + 0];
318
col2 = colLoop + steps[stepLoop*2 + 1];
320
if( col2>=0 && row2>=0 && col2<srcImageSize.width && row2<srcImageSize.height )
322
int memoryStep = memorySteps[ stepLoop ];
323
int pixelValue2 = destGLCM->forwardLookupTable[ srcImageData[ lineOffset + colLoop + memoryStep ] ];
326
matrices[stepLoop][pixelValue1][pixelValue2] ++;
327
matrices[stepLoop][pixelValue2][pixelValue1] ++;
329
// incremenet counter of total number of increments
330
stepIncrementsCounter[stepLoop] += 2;
336
// normalize matrices. each element is a probability of gray value i,j adjacency in direction/magnitude k
337
for( sideLoop1=0; sideLoop1<matrixSideLength; sideLoop1++ )
339
for( sideLoop2=0; sideLoop2<matrixSideLength; sideLoop2++ )
341
for( stepLoop=0; stepLoop<numSteps; stepLoop++ )
343
matrices[stepLoop][sideLoop1][sideLoop2] /= double(stepIncrementsCounter[stepLoop]);
348
destGLCM->matrices = matrices;
352
cvFree( &stepIncrementsCounter );
354
if( cvGetErrStatus() < 0 )
355
cvReleaseGLCM( &destGLCM, CV_GLCM_GLCM );
360
cvCreateGLCMDescriptors( CvGLCM* destGLCM, int descriptorOptimizationType )
362
CV_FUNCNAME( "cvCreateGLCMDescriptors" );
369
CV_ERROR( CV_StsNullPtr, "" );
371
if( !(destGLCM->matrices) )
372
CV_ERROR( CV_StsNullPtr, "Matrices are not allocated" );
374
CV_CALL( cvReleaseGLCM( &destGLCM, CV_GLCM_DESC ));
376
if( destGLCM->optimizationType != CV_GLCM_OPTIMIZATION_HISTOGRAM )
378
destGLCM->descriptorOptimizationType = destGLCM->numDescriptors = descriptorOptimizationType;
382
CV_ERROR( CV_StsBadFlag, "Histogram-based method is not implemented" );
383
// destGLCM->descriptorOptimizationType = destGLCM->numDescriptors = CV_GLCMDESC_OPTIMIZATION_HISTOGRAM;
386
CV_CALL( destGLCM->descriptors = (double**)
387
cvAlloc( destGLCM->numMatrices*sizeof(destGLCM->descriptors[0])));
389
for( matrixLoop = 0; matrixLoop < destGLCM->numMatrices; matrixLoop ++ )
391
CV_CALL( destGLCM->descriptors[ matrixLoop ] =
392
(double*)cvAlloc( destGLCM->numDescriptors*sizeof(destGLCM->descriptors[0][0])));
393
memset( destGLCM->descriptors[matrixLoop], 0, destGLCM->numDescriptors*sizeof(double) );
395
switch( destGLCM->descriptorOptimizationType )
397
case CV_GLCMDESC_OPTIMIZATION_ALLOWDOUBLENEST:
398
icvCreateGLCMDescriptors_AllowDoubleNest( destGLCM, matrixLoop );
401
CV_ERROR( CV_StsBadFlag,
402
"descriptorOptimizationType different from CV_GLCMDESC_OPTIMIZATION_ALLOWDOUBLENEST\n"
403
"is not supported" );
405
case CV_GLCMDESC_OPTIMIZATION_ALLOWTRIPLENEST:
406
icvCreateGLCMDescriptors_AllowTripleNest( destGLCM, matrixLoop );
408
case CV_GLCMDESC_OPTIMIZATION_HISTOGRAM:
409
if(matrixLoop < destGLCM->numMatrices>>1)
410
icvCreateGLCMDescriptors_Histogram( destGLCM, matrixLoop);
418
if( cvGetErrStatus() < 0 )
419
cvReleaseGLCM( &destGLCM, CV_GLCM_DESC );
424
icvCreateGLCMDescriptors_AllowDoubleNest( CvGLCM* destGLCM, int matrixIndex )
426
int sideLoop1, sideLoop2;
427
int matrixSideLength = destGLCM->matrixSideLength;
429
double** matrix = destGLCM->matrices[ matrixIndex ];
430
double* descriptors = destGLCM->descriptors[ matrixIndex ];
432
double* marginalProbability =
433
(double*)cvAlloc( matrixSideLength * sizeof(marginalProbability[0]));
434
memset( marginalProbability, 0, matrixSideLength * sizeof(double) );
436
double maximumProbability = 0;
437
double marginalProbabilityEntropy = 0;
438
double correlationMean = 0, correlationStdDeviation = 0, correlationProductTerm = 0;
440
for( sideLoop1=0; sideLoop1<matrixSideLength; sideLoop1++ )
442
int actualSideLoop1 = destGLCM->reverseLookupTable[ sideLoop1 ];
444
for( sideLoop2=0; sideLoop2<matrixSideLength; sideLoop2++ )
446
double entryValue = matrix[ sideLoop1 ][ sideLoop2 ];
448
int actualSideLoop2 = destGLCM->reverseLookupTable[ sideLoop2 ];
449
int sideLoopDifference = actualSideLoop1 - actualSideLoop2;
450
int sideLoopDifferenceSquared = sideLoopDifference*sideLoopDifference;
452
marginalProbability[ sideLoop1 ] += entryValue;
453
correlationMean += actualSideLoop1*entryValue;
455
maximumProbability = MAX( maximumProbability, entryValue );
457
if( actualSideLoop2 > actualSideLoop1 )
459
descriptors[ CV_GLCMDESC_CONTRAST ] += sideLoopDifferenceSquared * entryValue;
462
descriptors[ CV_GLCMDESC_HOMOGENITY ] += entryValue / ( 1.0 + sideLoopDifferenceSquared );
466
descriptors[ CV_GLCMDESC_ENTROPY ] += entryValue * log( entryValue );
469
descriptors[ CV_GLCMDESC_ENERGY ] += entryValue*entryValue;
472
if( marginalProbability[ actualSideLoop1 ] > 0 )
473
marginalProbabilityEntropy += marginalProbability[ actualSideLoop1 ]*log(marginalProbability[ actualSideLoop1 ]);
476
marginalProbabilityEntropy = -marginalProbabilityEntropy;
478
descriptors[ CV_GLCMDESC_CONTRAST ] += descriptors[ CV_GLCMDESC_CONTRAST ];
479
descriptors[ CV_GLCMDESC_ENTROPY ] = -descriptors[ CV_GLCMDESC_ENTROPY ];
480
descriptors[ CV_GLCMDESC_MAXIMUMPROBABILITY ] = maximumProbability;
482
double HXY = 0, HXY1 = 0, HXY2 = 0;
484
HXY = descriptors[ CV_GLCMDESC_ENTROPY ];
486
for( sideLoop1=0; sideLoop1<matrixSideLength; sideLoop1++ )
488
double sideEntryValueSum = 0;
489
int actualSideLoop1 = destGLCM->reverseLookupTable[ sideLoop1 ];
491
for( sideLoop2=0; sideLoop2<matrixSideLength; sideLoop2++ )
493
double entryValue = matrix[ sideLoop1 ][ sideLoop2 ];
495
sideEntryValueSum += entryValue;
497
int actualSideLoop2 = destGLCM->reverseLookupTable[ sideLoop2 ];
499
correlationProductTerm += (actualSideLoop1 - correlationMean) * (actualSideLoop2 - correlationMean) * entryValue;
501
double clusterTerm = actualSideLoop1 + actualSideLoop2 - correlationMean - correlationMean;
503
descriptors[ CV_GLCMDESC_CLUSTERTENDENCY ] += clusterTerm * clusterTerm * entryValue;
504
descriptors[ CV_GLCMDESC_CLUSTERSHADE ] += clusterTerm * clusterTerm * clusterTerm * entryValue;
506
double HXYValue = marginalProbability[ actualSideLoop1 ] * marginalProbability[ actualSideLoop2 ];
509
double HXYValueLog = log( HXYValue );
510
HXY1 += entryValue * HXYValueLog;
511
HXY2 += HXYValue * HXYValueLog;
515
correlationStdDeviation += (actualSideLoop1-correlationMean) * (actualSideLoop1-correlationMean) * sideEntryValueSum;
521
descriptors[ CV_GLCMDESC_CORRELATIONINFO1 ] = ( HXY - HXY1 ) / ( correlationMean );
522
descriptors[ CV_GLCMDESC_CORRELATIONINFO2 ] = sqrt( 1.0 - exp( -2.0 * (HXY2 - HXY ) ) );
524
correlationStdDeviation = sqrt( correlationStdDeviation );
526
descriptors[ CV_GLCMDESC_CORRELATION ] = correlationProductTerm / (correlationStdDeviation*correlationStdDeviation );
528
delete [] marginalProbability;
532
CV_IMPL double cvGetGLCMDescriptor( CvGLCM* GLCM, int step, int descriptor )
534
double value = DBL_MAX;
536
CV_FUNCNAME( "cvGetGLCMDescriptor" );
541
CV_ERROR( CV_StsNullPtr, "" );
543
if( !(GLCM->descriptors) )
544
CV_ERROR( CV_StsNullPtr, "" );
546
if( (unsigned)step >= (unsigned)(GLCM->numMatrices))
547
CV_ERROR( CV_StsOutOfRange, "step is not in 0 .. GLCM->numMatrices - 1" );
549
if( (unsigned)descriptor >= (unsigned)(GLCM->numDescriptors))
550
CV_ERROR( CV_StsOutOfRange, "descriptor is not in 0 .. GLCM->numDescriptors - 1" );
552
value = GLCM->descriptors[step][descriptor];
561
cvGetGLCMDescriptorStatistics( CvGLCM* GLCM, int descriptor,
562
double* _average, double* _standardDeviation )
564
CV_FUNCNAME( "cvGetGLCMDescriptorStatistics" );
569
if( _standardDeviation )
570
*_standardDeviation = DBL_MAX;
574
int matrixLoop, numMatrices;
575
double average = 0, squareSum = 0;
578
CV_ERROR( CV_StsNullPtr, "" );
580
if( !(GLCM->descriptors))
581
CV_ERROR( CV_StsNullPtr, "Descriptors are not calculated" );
583
if( (unsigned)descriptor >= (unsigned)(GLCM->numDescriptors) )
584
CV_ERROR( CV_StsOutOfRange, "Descriptor index is out of range" );
586
numMatrices = GLCM->numMatrices;
588
for( matrixLoop = 0; matrixLoop < numMatrices; matrixLoop++ )
590
double temp = GLCM->descriptors[ matrixLoop ][ descriptor ];
592
squareSum += temp*temp;
595
average /= numMatrices;
600
if( _standardDeviation )
601
*_standardDeviation = sqrt( (squareSum - average*average*numMatrices)/(numMatrices-1));
608
cvCreateGLCMImage( CvGLCM* GLCM, int step )
612
CV_FUNCNAME( "cvCreateGLCMImage" );
617
int sideLoop1, sideLoop2;
620
CV_ERROR( CV_StsNullPtr, "" );
622
if( !(GLCM->matrices) )
623
CV_ERROR( CV_StsNullPtr, "Matrices are not allocated" );
625
if( (unsigned)step >= (unsigned)(GLCM->numMatrices) )
626
CV_ERROR( CV_StsOutOfRange, "The step index is out of range" );
628
dest = cvCreateImage( cvSize( GLCM->matrixSideLength, GLCM->matrixSideLength ), IPL_DEPTH_32F, 1 );
629
destData = (float*)(dest->imageData);
631
for( sideLoop1 = 0; sideLoop1 < GLCM->matrixSideLength;
632
sideLoop1++, (float*&)destData += dest->widthStep )
634
for( sideLoop2=0; sideLoop2 < GLCM->matrixSideLength; sideLoop2++ )
636
double matrixValue = GLCM->matrices[step][sideLoop1][sideLoop2];
637
destData[ sideLoop2 ] = (float)matrixValue;
643
if( cvGetErrStatus() < 0 )
644
cvReleaseImage( &dest );
added
removed
modules/legacy/src/texture.cpp
