git.maemo.org Git - opencv/blob - cvaux/src/cvtexture.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   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.
   8 //
   9 //
  10 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  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.
  25 //
  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.
  28 //
  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.
  39 //
  40 //M*/
  41
  42 /****************************************************************************************\
  43
  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]
  46
  47 \****************************************************************************************/
  48
  49 #include "_cvaux.h"
  50
  51 #include <math.h>
  52 #include <assert.h>
  53
  54 #define CV_MAX_NUM_GREY_LEVELS_8U  256
  55
  56 struct CvGLCM
  57 {
  58     int matrixSideLength;
  59     int numMatrices;
  60     double*** matrices;
  61
  62     int  numLookupTableElements;
  63     int  forwardLookupTable[CV_MAX_NUM_GREY_LEVELS_8U];
  64     int  reverseLookupTable[CV_MAX_NUM_GREY_LEVELS_8U];
  65
  66     double** descriptors;
  67     int numDescriptors;
  68     int descriptorOptimizationType;
  69     int optimizationType;
  70 };
  71
  72
  73 static void icvCreateGLCM_LookupTable_8u_C1R( const uchar* srcImageData, int srcImageStep,
  74                                              CvSize srcImageSize, CvGLCM* destGLCM,
  75                                              int* steps, int numSteps, int* memorySteps );
  76
  77 static void
  78 icvCreateGLCMDescriptors_AllowDoubleNest( CvGLCM* destGLCM, int matrixIndex );
  79
  80
  81 CV_IMPL CvGLCM*
  82 cvCreateGLCM( const IplImage* srcImage,
  83               int stepMagnitude,
  84               const int* srcStepDirections,/* should be static array..
  85                                           or if not the user should handle de-allocation */
  86               int numStepDirections,
  87               int optimizationType )
  88 {
  89     static const int defaultStepDirections[] = { 0,1, -1,1, -1,0, -1,-1 };
  90
  91     int* memorySteps = 0;
  92     CvGLCM* newGLCM = 0;
  93     int* stepDirections = 0;
  94
  95     CV_FUNCNAME( "cvCreateGLCM" );
  96
  97     __BEGIN__;
  98
  99     uchar* srcImageData = 0;
 100     CvSize srcImageSize;
 101     int srcImageStep;
 102     int stepLoop;
 103     const int maxNumGreyLevels8u = CV_MAX_NUM_GREY_LEVELS_8U;
 104
 105     if( !srcImage )
 106         CV_ERROR( CV_StsNullPtr, "" );
 107
 108     if( srcImage->nChannels != 1 )
 109         CV_ERROR( CV_BadNumChannels, "Number of channels must be 1");
 110
 111     if( srcImage->depth != IPL_DEPTH_8U )
 112         CV_ERROR( CV_BadDepth, "Depth must be equal IPL_DEPTH_8U");
 113
 114     // no Directions provided, use the default ones - 0 deg, 45, 90, 135
 115     if( !srcStepDirections )
 116     {
 117         srcStepDirections = defaultStepDirections;
 118     }
 119
 120     CV_CALL( stepDirections = (int*)cvAlloc( numStepDirections*2*sizeof(stepDirections[0])));
 121     memcpy( stepDirections, srcStepDirections, numStepDirections*2*sizeof(stepDirections[0]));
 122
 123     cvGetImageRawData( srcImage, &srcImageData, &srcImageStep, &srcImageSize );
 124
 125     // roll together Directions and magnitudes together with knowledge of image (step)
 126     CV_CALL( memorySteps = (int*)cvAlloc( numStepDirections*sizeof(memorySteps[0])));
 127
 128     for( stepLoop = 0; stepLoop < numStepDirections; stepLoop++ )
 129     {
 130         stepDirections[stepLoop*2 + 0] *= stepMagnitude;
 131         stepDirections[stepLoop*2 + 1] *= stepMagnitude;
 132
 133         memorySteps[stepLoop] = stepDirections[stepLoop*2 + 0]*srcImageStep +
 134                                 stepDirections[stepLoop*2 + 1];
 135     }
 136
 137     CV_CALL( newGLCM = (CvGLCM*)cvAlloc(sizeof(newGLCM)));
 138     memset( newGLCM, 0, sizeof(newGLCM) );
 139
 140     newGLCM->matrices = 0;
 141     newGLCM->numMatrices = numStepDirections;
 142     newGLCM->optimizationType = optimizationType;
 143
 144     if( optimizationType <= CV_GLCM_OPTIMIZATION_LUT )
 145     {
 146         int lookupTableLoop, imageColLoop, imageRowLoop, lineOffset = 0;
 147
 148         // if optimization type is set to lut, then make one for the image
 149         if( optimizationType == CV_GLCM_OPTIMIZATION_LUT )
 150         {
 151             for( imageRowLoop = 0; imageRowLoop < srcImageSize.height;
 152                                    imageRowLoop++, lineOffset += srcImageStep )
 153             {
 154                 for( imageColLoop = 0; imageColLoop < srcImageSize.width; imageColLoop++ )
 155                 {
 156                     newGLCM->forwardLookupTable[srcImageData[lineOffset+imageColLoop]]=1;
 157                 }
 158             }
 159
 160             newGLCM->numLookupTableElements = 0;
 161
 162             for( lookupTableLoop = 0; lookupTableLoop < maxNumGreyLevels8u; lookupTableLoop++ )
 163             {
 164                 if( newGLCM->forwardLookupTable[ lookupTableLoop ] != 0 )
 165                 {
 166                     newGLCM->forwardLookupTable[ lookupTableLoop ] =
 167                         newGLCM->numLookupTableElements;
 168                     newGLCM->reverseLookupTable[ newGLCM->numLookupTableElements ] =
 169                         lookupTableLoop;
 170
 171                     newGLCM->numLookupTableElements++;
 172                 }
 173             }
 174         }
 175         // otherwise make a "LUT" which contains all the gray-levels (for code-reuse)
 176         else if( optimizationType == CV_GLCM_OPTIMIZATION_NONE )
 177         {
 178             for( lookupTableLoop = 0; lookupTableLoop <maxNumGreyLevels8u; lookupTableLoop++ )
 179             {
 180                 newGLCM->forwardLookupTable[ lookupTableLoop ] = lookupTableLoop;
 181                 newGLCM->reverseLookupTable[ lookupTableLoop ] = lookupTableLoop;
 182             }
 183             newGLCM->numLookupTableElements = maxNumGreyLevels8u;
 184         }
 185
 186         newGLCM->matrixSideLength = newGLCM->numLookupTableElements;
 187         icvCreateGLCM_LookupTable_8u_C1R( srcImageData, srcImageStep, srcImageSize,
 188                                           newGLCM, stepDirections,
 189                                           numStepDirections, memorySteps );
 190     }
 191     else if( optimizationType == CV_GLCM_OPTIMIZATION_HISTOGRAM )
 192     {
 193         CV_ERROR( CV_StsBadFlag, "Histogram-based method is not implemented" );
 194
 195     /*  newGLCM->numMatrices *= 2;
 196         newGLCM->matrixSideLength = maxNumGreyLevels8u*2;
 197
 198         icvCreateGLCM_Histogram_8uC1R( srcImageStep, srcImageSize, srcImageData,
 199                                        newGLCM, numStepDirections,
 200                                        stepDirections, memorySteps );
 201     */
 202     }
 203
 204     __END__;
 205
 206     cvFree( &memorySteps );
 207     cvFree( &stepDirections );
 208
 209     if( cvGetErrStatus() < 0 )
 210     {
 211         cvFree( &newGLCM );
 212     }
 213
 214     return newGLCM;
 215 }
 216
 217
 218 CV_IMPL void
 219 cvReleaseGLCM( CvGLCM** GLCM, int flag )
 220 {
 221     CV_FUNCNAME( "cvReleaseGLCM" );
 222
 223     __BEGIN__;
 224
 225     int matrixLoop;
 226
 227     if( !GLCM )
 228         CV_ERROR( CV_StsNullPtr, "" );
 229
 230     if( *GLCM )
 231         EXIT; // repeated deallocation: just skip it.
 232
 233     if( (flag == CV_GLCM_GLCM || flag == CV_GLCM_ALL) && (*GLCM)->matrices )
 234     {
 235         for( matrixLoop = 0; matrixLoop < (*GLCM)->numMatrices; matrixLoop++ )
 236         {
 237             if( (*GLCM)->matrices[ matrixLoop ] )
 238             {
 239                 cvFree( (*GLCM)->matrices[matrixLoop] );
 240                 cvFree( (*GLCM)->matrices + matrixLoop );
 241             }
 242         }
 243
 244         cvFree( &((*GLCM)->matrices) );
 245     }
 246
 247     if( (flag == CV_GLCM_DESC || flag == CV_GLCM_ALL) && (*GLCM)->descriptors )
 248     {
 249         for( matrixLoop = 0; matrixLoop < (*GLCM)->numMatrices; matrixLoop++ )
 250         {
 251             cvFree( (*GLCM)->descriptors + matrixLoop );
 252         }
 253         cvFree( &((*GLCM)->descriptors) );
 254     }
 255
 256     if( flag == CV_GLCM_ALL )
 257     {
 258         cvFree( GLCM );
 259     }
 260
 261     __END__;
 262 }
 263
 264
 265 static void
 266 icvCreateGLCM_LookupTable_8u_C1R( const uchar* srcImageData,
 267                                   int srcImageStep,
 268                                   CvSize srcImageSize,
 269                                   CvGLCM* destGLCM,
 270                                   int* steps,
 271                                   int numSteps,
 272                                   int* memorySteps )
 273 {
 274     int* stepIncrementsCounter = 0;
 275
 276     CV_FUNCNAME( "icvCreateGLCM_LookupTable_8u_C1R" );
 277
 278     __BEGIN__;
 279
 280     int matrixSideLength = destGLCM->matrixSideLength;
 281     int stepLoop, sideLoop1, sideLoop2;
 282     int colLoop, rowLoop, lineOffset = 0;
 283     double*** matrices = 0;
 284
 285     // allocate memory to the matrices
 286     CV_CALL( destGLCM->matrices = (double***)cvAlloc( sizeof(matrices[0])*numSteps ));
 287     matrices = destGLCM->matrices;
 288
 289     for( stepLoop=0; stepLoop<numSteps; stepLoop++ )
 290     {
 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 ));
 294
 295         memset( matrices[stepLoop][0], 0, matrixSideLength*matrixSideLength*
 296                                           sizeof(matrices[0][0]) );
 297
 298         for( sideLoop1 = 1; sideLoop1 < matrixSideLength; sideLoop1++ )
 299         {
 300             matrices[stepLoop][sideLoop1] = matrices[stepLoop][sideLoop1-1] + matrixSideLength;
 301         }
 302     }
 303
 304     CV_CALL( stepIncrementsCounter = (int*)cvAlloc( numSteps*sizeof(stepIncrementsCounter[0])));
 305     memset( stepIncrementsCounter, 0, numSteps*sizeof(stepIncrementsCounter[0]) );
 306
 307     // generate GLCM for each step
 308     for( rowLoop=0; rowLoop<srcImageSize.height; rowLoop++, lineOffset+=srcImageStep )
 309     {
 310         for( colLoop=0; colLoop<srcImageSize.width; colLoop++ )
 311         {
 312             int pixelValue1 = destGLCM->forwardLookupTable[srcImageData[lineOffset + colLoop]];
 313
 314             for( stepLoop=0; stepLoop<numSteps; stepLoop++ )
 315             {
 316                 int col2, row2;
 317                 row2 = rowLoop + steps[stepLoop*2 + 0];
 318                 col2 = colLoop + steps[stepLoop*2 + 1];
 319
 320                 if( col2>=0 && row2>=0 && col2<srcImageSize.width && row2<srcImageSize.height )
 321                 {
 322                     int memoryStep = memorySteps[ stepLoop ];
 323                     int pixelValue2 = destGLCM->forwardLookupTable[ srcImageData[ lineOffset + colLoop + memoryStep ] ];
 324
 325                     // maintain symmetry
 326                     matrices[stepLoop][pixelValue1][pixelValue2] ++;
 327                     matrices[stepLoop][pixelValue2][pixelValue1] ++;
 328
 329                     // incremenet counter of total number of increments
 330                     stepIncrementsCounter[stepLoop] += 2;
 331                 }
 332             }
 333         }
 334     }
 335
 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++ )
 338     {
 339         for( sideLoop2=0; sideLoop2<matrixSideLength; sideLoop2++ )
 340         {
 341             for( stepLoop=0; stepLoop<numSteps; stepLoop++ )
 342             {
 343                 matrices[stepLoop][sideLoop1][sideLoop2] /= double(stepIncrementsCounter[stepLoop]);
 344             }
 345         }
 346     }
 347
 348     destGLCM->matrices = matrices;
 349
 350     __END__;
 351
 352     cvFree( &stepIncrementsCounter );
 353
 354     if( cvGetErrStatus() < 0 )
 355         cvReleaseGLCM( &destGLCM, CV_GLCM_GLCM );
 356 }
 357
 358
 359 CV_IMPL void
 360 cvCreateGLCMDescriptors( CvGLCM* destGLCM, int descriptorOptimizationType )
 361 {
 362     CV_FUNCNAME( "cvCreateGLCMDescriptors" );
 363
 364     __BEGIN__;
 365
 366     int matrixLoop;
 367
 368     if( !destGLCM )
 369         CV_ERROR( CV_StsNullPtr, "" );
 370
 371     if( !(destGLCM->matrices) )
 372         CV_ERROR( CV_StsNullPtr, "Matrices are not allocated" );
 373
 374     CV_CALL( cvReleaseGLCM( &destGLCM, CV_GLCM_DESC ));
 375
 376     if( destGLCM->optimizationType != CV_GLCM_OPTIMIZATION_HISTOGRAM )
 377     {
 378         destGLCM->descriptorOptimizationType = destGLCM->numDescriptors = descriptorOptimizationType;
 379     }
 380     else
 381     {
 382         CV_ERROR( CV_StsBadFlag, "Histogram-based method is not implemented" );
 383 //      destGLCM->descriptorOptimizationType = destGLCM->numDescriptors = CV_GLCMDESC_OPTIMIZATION_HISTOGRAM;
 384     }
 385
 386     CV_CALL( destGLCM->descriptors = (double**)
 387             cvAlloc( destGLCM->numMatrices*sizeof(destGLCM->descriptors[0])));
 388
 389     for( matrixLoop = 0; matrixLoop < destGLCM->numMatrices; matrixLoop ++ )
 390     {
 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) );
 394
 395         switch( destGLCM->descriptorOptimizationType )
 396         {
 397             case CV_GLCMDESC_OPTIMIZATION_ALLOWDOUBLENEST:
 398                 icvCreateGLCMDescriptors_AllowDoubleNest( destGLCM, matrixLoop );
 399                 break;
 400             default:
 401                 CV_ERROR( CV_StsBadFlag,
 402                 "descriptorOptimizationType different from CV_GLCMDESC_OPTIMIZATION_ALLOWDOUBLENEST\n"
 403                 "is not supported" );
 404             /*
 405             case CV_GLCMDESC_OPTIMIZATION_ALLOWTRIPLENEST:
 406                 icvCreateGLCMDescriptors_AllowTripleNest( destGLCM, matrixLoop );
 407                 break;
 408             case CV_GLCMDESC_OPTIMIZATION_HISTOGRAM:
 409                 if(matrixLoop < destGLCM->numMatrices>>1)
 410                     icvCreateGLCMDescriptors_Histogram( destGLCM, matrixLoop);
 411                     break;
 412             */
 413         }
 414     }
 415
 416     __END__;
 417
 418     if( cvGetErrStatus() < 0 )
 419         cvReleaseGLCM( &destGLCM, CV_GLCM_DESC );
 420 }
 421
 422
 423 static void
 424 icvCreateGLCMDescriptors_AllowDoubleNest( CvGLCM* destGLCM, int matrixIndex )
 425 {
 426     int sideLoop1, sideLoop2;
 427     int matrixSideLength = destGLCM->matrixSideLength;
 428
 429     double** matrix = destGLCM->matrices[ matrixIndex ];
 430     double* descriptors = destGLCM->descriptors[ matrixIndex ];
 431
 432     double* marginalProbability =
 433         (double*)cvAlloc( matrixSideLength * sizeof(marginalProbability[0]));
 434     memset( marginalProbability, 0, matrixSideLength * sizeof(double) );
 435
 436     double maximumProbability = 0;
 437     double marginalProbabilityEntropy = 0;
 438     double correlationMean = 0, correlationStdDeviation = 0, correlationProductTerm = 0;
 439
 440     for( sideLoop1=0; sideLoop1<matrixSideLength; sideLoop1++ )
 441     {
 442         int actualSideLoop1 = destGLCM->reverseLookupTable[ sideLoop1 ];
 443
 444         for( sideLoop2=0; sideLoop2<matrixSideLength; sideLoop2++ )
 445         {
 446             double entryValue = matrix[ sideLoop1 ][ sideLoop2 ];
 447
 448             int actualSideLoop2 = destGLCM->reverseLookupTable[ sideLoop2 ];
 449             int sideLoopDifference = actualSideLoop1 - actualSideLoop2;
 450             int sideLoopDifferenceSquared = sideLoopDifference*sideLoopDifference;
 451
 452             marginalProbability[ sideLoop1 ] += entryValue;
 453             correlationMean += actualSideLoop1*entryValue;
 454
 455             maximumProbability = MAX( maximumProbability, entryValue );
 456
 457             if( actualSideLoop2 > actualSideLoop1 )
 458             {
 459                 descriptors[ CV_GLCMDESC_CONTRAST ] += sideLoopDifferenceSquared * entryValue;
 460             }
 461
 462             descriptors[ CV_GLCMDESC_HOMOGENITY ] += entryValue / ( 1.0 + sideLoopDifferenceSquared );
 463
 464             if( entryValue > 0 )
 465             {
 466                 descriptors[ CV_GLCMDESC_ENTROPY ] += entryValue * log( entryValue );
 467             }
 468
 469             descriptors[ CV_GLCMDESC_ENERGY ] += entryValue*entryValue;
 470         }
 471
 472         if( marginalProbability>0 )
 473             marginalProbabilityEntropy += marginalProbability[ actualSideLoop1 ]*log(marginalProbability[ actualSideLoop1 ]);
 474     }
 475
 476     marginalProbabilityEntropy = -marginalProbabilityEntropy;
 477
 478     descriptors[ CV_GLCMDESC_CONTRAST ] += descriptors[ CV_GLCMDESC_CONTRAST ];
 479     descriptors[ CV_GLCMDESC_ENTROPY ] = -descriptors[ CV_GLCMDESC_ENTROPY ];
 480     descriptors[ CV_GLCMDESC_MAXIMUMPROBABILITY ] = maximumProbability;
 481
 482     double HXY = 0, HXY1 = 0, HXY2 = 0;
 483
 484     HXY = descriptors[ CV_GLCMDESC_ENTROPY ];
 485
 486     for( sideLoop1=0; sideLoop1<matrixSideLength; sideLoop1++ )
 487     {
 488         double sideEntryValueSum = 0;
 489         int actualSideLoop1 = destGLCM->reverseLookupTable[ sideLoop1 ];
 490
 491         for( sideLoop2=0; sideLoop2<matrixSideLength; sideLoop2++ )
 492         {
 493             double entryValue = matrix[ sideLoop1 ][ sideLoop2 ];
 494
 495             sideEntryValueSum += entryValue;
 496
 497             int actualSideLoop2 = destGLCM->reverseLookupTable[ sideLoop2 ];
 498
 499             correlationProductTerm += (actualSideLoop1 - correlationMean) * (actualSideLoop2 - correlationMean) * entryValue;
 500
 501             double clusterTerm = actualSideLoop1 + actualSideLoop2 - correlationMean - correlationMean;
 502
 503             descriptors[ CV_GLCMDESC_CLUSTERTENDENCY ] += clusterTerm * clusterTerm * entryValue;
 504             descriptors[ CV_GLCMDESC_CLUSTERSHADE ] += clusterTerm * clusterTerm * clusterTerm * entryValue;
 505
 506             double HXYValue = marginalProbability[ actualSideLoop1 ] * marginalProbability[ actualSideLoop2 ];
 507             if( HXYValue>0 )
 508             {
 509                 double HXYValueLog = log( HXYValue );
 510                 HXY1 += entryValue * HXYValueLog;
 511                 HXY2 += HXYValue * HXYValueLog;
 512             }
 513         }
 514
 515         correlationStdDeviation += (actualSideLoop1-correlationMean) * (actualSideLoop1-correlationMean) * sideEntryValueSum;
 516     }
 517
 518     HXY1 =- HXY1;
 519     HXY2 =- HXY2;
 520
 521     descriptors[ CV_GLCMDESC_CORRELATIONINFO1 ] = ( HXY - HXY1 ) / ( correlationMean );
 522     descriptors[ CV_GLCMDESC_CORRELATIONINFO2 ] = sqrt( 1.0 - exp( -2.0 * (HXY2 - HXY ) ) );
 523
 524     correlationStdDeviation = sqrt( correlationStdDeviation );
 525
 526     descriptors[ CV_GLCMDESC_CORRELATION ] = correlationProductTerm / (correlationStdDeviation*correlationStdDeviation );
 527
 528     delete [] marginalProbability;
 529 }
 530
 531
 532 CV_IMPL double cvGetGLCMDescriptor( CvGLCM* GLCM, int step, int descriptor )
 533 {
 534     double value = DBL_MAX;
 535
 536     CV_FUNCNAME( "cvGetGLCMDescriptor" );
 537
 538     __BEGIN__;
 539
 540     if( !GLCM )
 541         CV_ERROR( CV_StsNullPtr, "" );
 542
 543     if( !(GLCM->descriptors) )
 544         CV_ERROR( CV_StsNullPtr, "" );
 545
 546     if( (unsigned)step >= (unsigned)(GLCM->numMatrices))
 547         CV_ERROR( CV_StsOutOfRange, "step is not in 0 .. GLCM->numMatrices - 1" );
 548
 549     if( (unsigned)descriptor >= (unsigned)(GLCM->numDescriptors))
 550         CV_ERROR( CV_StsOutOfRange, "descriptor is not in 0 .. GLCM->numDescriptors - 1" );
 551
 552     value = GLCM->descriptors[step][descriptor];
 553
 554     __END__;
 555
 556     return value;
 557 }
 558
 559
 560 CV_IMPL void
 561 cvGetGLCMDescriptorStatistics( CvGLCM* GLCM, int descriptor,
 562                                double* _average, double* _standardDeviation )
 563 {
 564     CV_FUNCNAME( "cvGetGLCMDescriptorStatistics" );
 565
 566     if( _average )
 567         *_average = DBL_MAX;
 568
 569     if( _standardDeviation )
 570         *_standardDeviation = DBL_MAX;
 571
 572     __BEGIN__;
 573
 574     int matrixLoop, numMatrices;
 575     double average = 0, squareSum = 0;
 576
 577     if( !GLCM )
 578         CV_ERROR( CV_StsNullPtr, "" );
 579
 580     if( !(GLCM->descriptors))
 581         CV_ERROR( CV_StsNullPtr, "Descriptors are not calculated" );
 582
 583     if( (unsigned)descriptor >= (unsigned)(GLCM->numDescriptors) )
 584         CV_ERROR( CV_StsOutOfRange, "Descriptor index is out of range" );
 585
 586     numMatrices = GLCM->numMatrices;
 587
 588     for( matrixLoop = 0; matrixLoop < numMatrices; matrixLoop++ )
 589     {
 590         double temp = GLCM->descriptors[ matrixLoop ][ descriptor ];
 591         average += temp;
 592         squareSum += temp*temp;
 593     }
 594
 595     average /= numMatrices;
 596
 597     if( _average )
 598         *_average = average;
 599
 600     if( _standardDeviation )
 601         *_standardDeviation = sqrt( (squareSum - average*average*numMatrices)/(numMatrices-1));
 602
 603     __END__;
 604 }
 605
 606
 607 CV_IMPL IplImage*
 608 cvCreateGLCMImage( CvGLCM* GLCM, int step )
 609 {
 610     IplImage* dest = 0;
 611
 612     CV_FUNCNAME( "cvCreateGLCMImage" );
 613
 614     __BEGIN__;
 615
 616     float* destData;
 617     int sideLoop1, sideLoop2;
 618
 619     if( !GLCM )
 620         CV_ERROR( CV_StsNullPtr, "" );
 621
 622     if( !(GLCM->matrices) )
 623         CV_ERROR( CV_StsNullPtr, "Matrices are not allocated" );
 624
 625     if( (unsigned)step >= (unsigned)(GLCM->numMatrices) )
 626         CV_ERROR( CV_StsOutOfRange, "The step index is out of range" );
 627
 628     dest = cvCreateImage( cvSize( GLCM->matrixSideLength, GLCM->matrixSideLength ), IPL_DEPTH_32F, 1 );
 629     destData = (float*)(dest->imageData);
 630
 631     for( sideLoop1 = 0; sideLoop1 < GLCM->matrixSideLength;
 632                         sideLoop1++, (float*&)destData += dest->widthStep )
 633     {
 634         for( sideLoop2=0; sideLoop2 < GLCM->matrixSideLength; sideLoop2++ )
 635         {
 636             double matrixValue = GLCM->matrices[step][sideLoop1][sideLoop2];
 637             destData[ sideLoop2 ] = (float)matrixValue;
 638         }
 639     }
 640
 641     __END__;
 642
 643     if( cvGetErrStatus() < 0 )
 644         cvReleaseImage( &dest );
 645
 646     return dest;
 647 }
 648