Update to 2.0.0 tree from current Fremantle build
[opencv] / src / cvaux / cvtexture.cpp
diff --git a/src/cvaux/cvtexture.cpp b/src/cvaux/cvtexture.cpp
new file mode 100644 (file)
index 0000000..434baaf
--- /dev/null
@@ -0,0 +1,648 @@
+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+/****************************************************************************************\
+
+      Calculation of a texture descriptors from GLCM (Grey Level Co-occurrence Matrix'es)
+      The code was submitted by Daniel Eaton [danieljameseaton@yahoo.com]
+
+\****************************************************************************************/
+
+#include "_cvaux.h"
+
+#include <math.h>
+#include <assert.h>
+
+#define CV_MAX_NUM_GREY_LEVELS_8U  256
+
+struct CvGLCM
+{
+    int matrixSideLength;
+    int numMatrices;
+    double*** matrices;
+
+    int  numLookupTableElements;
+    int  forwardLookupTable[CV_MAX_NUM_GREY_LEVELS_8U];
+    int  reverseLookupTable[CV_MAX_NUM_GREY_LEVELS_8U];
+
+    double** descriptors;
+    int numDescriptors;
+    int descriptorOptimizationType;
+    int optimizationType;
+};
+
+
+static void icvCreateGLCM_LookupTable_8u_C1R( const uchar* srcImageData, int srcImageStep,
+                                             CvSize srcImageSize, CvGLCM* destGLCM,
+                                             int* steps, int numSteps, int* memorySteps );
+
+static void
+icvCreateGLCMDescriptors_AllowDoubleNest( CvGLCM* destGLCM, int matrixIndex );
+
+
+CV_IMPL CvGLCM*
+cvCreateGLCM( const IplImage* srcImage,
+              int stepMagnitude,
+              const int* srcStepDirections,/* should be static array..
+                                          or if not the user should handle de-allocation */
+              int numStepDirections,
+              int optimizationType )
+{
+    static const int defaultStepDirections[] = { 0,1, -1,1, -1,0, -1,-1 };
+
+    int* memorySteps = 0;
+    CvGLCM* newGLCM = 0;
+    int* stepDirections = 0;
+
+    CV_FUNCNAME( "cvCreateGLCM" );
+
+    __BEGIN__;
+
+    uchar* srcImageData = 0;
+    CvSize srcImageSize;
+    int srcImageStep;
+    int stepLoop;
+    const int maxNumGreyLevels8u = CV_MAX_NUM_GREY_LEVELS_8U;
+
+    if( !srcImage )
+        CV_ERROR( CV_StsNullPtr, "" );
+
+    if( srcImage->nChannels != 1 )
+        CV_ERROR( CV_BadNumChannels, "Number of channels must be 1");
+
+    if( srcImage->depth != IPL_DEPTH_8U )
+        CV_ERROR( CV_BadDepth, "Depth must be equal IPL_DEPTH_8U");
+
+    // no Directions provided, use the default ones - 0 deg, 45, 90, 135
+    if( !srcStepDirections )
+    {
+        srcStepDirections = defaultStepDirections;
+    }
+
+    CV_CALL( stepDirections = (int*)cvAlloc( numStepDirections*2*sizeof(stepDirections[0])));
+    memcpy( stepDirections, srcStepDirections, numStepDirections*2*sizeof(stepDirections[0]));
+
+    cvGetImageRawData( srcImage, &srcImageData, &srcImageStep, &srcImageSize );
+
+    // roll together Directions and magnitudes together with knowledge of image (step)
+    CV_CALL( memorySteps = (int*)cvAlloc( numStepDirections*sizeof(memorySteps[0])));
+
+    for( stepLoop = 0; stepLoop < numStepDirections; stepLoop++ )
+    {
+        stepDirections[stepLoop*2 + 0] *= stepMagnitude;
+        stepDirections[stepLoop*2 + 1] *= stepMagnitude;
+
+        memorySteps[stepLoop] = stepDirections[stepLoop*2 + 0]*srcImageStep +
+                                stepDirections[stepLoop*2 + 1];
+    }
+
+    CV_CALL( newGLCM = (CvGLCM*)cvAlloc(sizeof(newGLCM)));
+    memset( newGLCM, 0, sizeof(newGLCM) );
+
+    newGLCM->matrices = 0;
+    newGLCM->numMatrices = numStepDirections;
+    newGLCM->optimizationType = optimizationType;
+
+    if( optimizationType <= CV_GLCM_OPTIMIZATION_LUT )
+    {
+        int lookupTableLoop, imageColLoop, imageRowLoop, lineOffset = 0;
+
+        // if optimization type is set to lut, then make one for the image
+        if( optimizationType == CV_GLCM_OPTIMIZATION_LUT )
+        {
+            for( imageRowLoop = 0; imageRowLoop < srcImageSize.height;
+                                   imageRowLoop++, lineOffset += srcImageStep )
+            {
+                for( imageColLoop = 0; imageColLoop < srcImageSize.width; imageColLoop++ )
+                {
+                    newGLCM->forwardLookupTable[srcImageData[lineOffset+imageColLoop]]=1;
+                }
+            }
+
+            newGLCM->numLookupTableElements = 0;
+
+            for( lookupTableLoop = 0; lookupTableLoop < maxNumGreyLevels8u; lookupTableLoop++ )
+            {
+                if( newGLCM->forwardLookupTable[ lookupTableLoop ] != 0 )
+                {
+                    newGLCM->forwardLookupTable[ lookupTableLoop ] =
+                        newGLCM->numLookupTableElements;
+                    newGLCM->reverseLookupTable[ newGLCM->numLookupTableElements ] =
+                        lookupTableLoop;
+
+                    newGLCM->numLookupTableElements++;
+                }
+            }
+        }
+        // otherwise make a "LUT" which contains all the gray-levels (for code-reuse)
+        else if( optimizationType == CV_GLCM_OPTIMIZATION_NONE )
+        {
+            for( lookupTableLoop = 0; lookupTableLoop <maxNumGreyLevels8u; lookupTableLoop++ )
+            {
+                newGLCM->forwardLookupTable[ lookupTableLoop ] = lookupTableLoop;
+                newGLCM->reverseLookupTable[ lookupTableLoop ] = lookupTableLoop;
+            }
+            newGLCM->numLookupTableElements = maxNumGreyLevels8u;
+        }
+
+        newGLCM->matrixSideLength = newGLCM->numLookupTableElements;
+        icvCreateGLCM_LookupTable_8u_C1R( srcImageData, srcImageStep, srcImageSize,
+                                          newGLCM, stepDirections,
+                                          numStepDirections, memorySteps );
+    }
+    else if( optimizationType == CV_GLCM_OPTIMIZATION_HISTOGRAM )
+    {
+        CV_ERROR( CV_StsBadFlag, "Histogram-based method is not implemented" );
+
+    /*  newGLCM->numMatrices *= 2;
+        newGLCM->matrixSideLength = maxNumGreyLevels8u*2;
+
+        icvCreateGLCM_Histogram_8uC1R( srcImageStep, srcImageSize, srcImageData,
+                                       newGLCM, numStepDirections,
+                                       stepDirections, memorySteps );
+    */
+    }
+
+    __END__;
+
+    cvFree( &memorySteps );
+    cvFree( &stepDirections );
+
+    if( cvGetErrStatus() < 0 )
+    {
+        cvFree( &newGLCM );
+    }
+
+    return newGLCM;
+}
+
+
+CV_IMPL void
+cvReleaseGLCM( CvGLCM** GLCM, int flag )
+{
+    CV_FUNCNAME( "cvReleaseGLCM" );
+
+    __BEGIN__;
+
+    int matrixLoop;
+
+    if( !GLCM )
+        CV_ERROR( CV_StsNullPtr, "" );
+
+    if( *GLCM )
+        EXIT; // repeated deallocation: just skip it.
+
+    if( (flag == CV_GLCM_GLCM || flag == CV_GLCM_ALL) && (*GLCM)->matrices )
+    {
+        for( matrixLoop = 0; matrixLoop < (*GLCM)->numMatrices; matrixLoop++ )
+        {
+            if( (*GLCM)->matrices[ matrixLoop ] )
+            {
+                cvFree( (*GLCM)->matrices[matrixLoop] );
+                cvFree( (*GLCM)->matrices + matrixLoop );
+            }
+        }
+
+        cvFree( &((*GLCM)->matrices) );
+    }
+
+    if( (flag == CV_GLCM_DESC || flag == CV_GLCM_ALL) && (*GLCM)->descriptors )
+    {
+        for( matrixLoop = 0; matrixLoop < (*GLCM)->numMatrices; matrixLoop++ )
+        {
+            cvFree( (*GLCM)->descriptors + matrixLoop );
+        }
+        cvFree( &((*GLCM)->descriptors) );
+    }
+
+    if( flag == CV_GLCM_ALL )
+    {
+        cvFree( GLCM );
+    }
+
+    __END__;
+}
+
+
+static void
+icvCreateGLCM_LookupTable_8u_C1R( const uchar* srcImageData,
+                                  int srcImageStep,
+                                  CvSize srcImageSize,
+                                  CvGLCM* destGLCM,
+                                  int* steps,
+                                  int numSteps,
+                                  int* memorySteps )
+{
+    int* stepIncrementsCounter = 0;
+
+    CV_FUNCNAME( "icvCreateGLCM_LookupTable_8u_C1R" );
+
+    __BEGIN__;
+
+    int matrixSideLength = destGLCM->matrixSideLength;
+    int stepLoop, sideLoop1, sideLoop2;
+    int colLoop, rowLoop, lineOffset = 0;
+    double*** matrices = 0;
+
+    // allocate memory to the matrices
+    CV_CALL( destGLCM->matrices = (double***)cvAlloc( sizeof(matrices[0])*numSteps ));
+    matrices = destGLCM->matrices;
+
+    for( stepLoop=0; stepLoop<numSteps; stepLoop++ )
+    {
+        CV_CALL( matrices[stepLoop] = (double**)cvAlloc( sizeof(matrices[0])*matrixSideLength ));
+        CV_CALL( matrices[stepLoop][0] = (double*)cvAlloc( sizeof(matrices[0][0])*
+                                                  matrixSideLength*matrixSideLength ));
+
+        memset( matrices[stepLoop][0], 0, matrixSideLength*matrixSideLength*
+                                          sizeof(matrices[0][0]) );
+
+        for( sideLoop1 = 1; sideLoop1 < matrixSideLength; sideLoop1++ )
+        {
+            matrices[stepLoop][sideLoop1] = matrices[stepLoop][sideLoop1-1] + matrixSideLength;
+        }
+    }
+
+    CV_CALL( stepIncrementsCounter = (int*)cvAlloc( numSteps*sizeof(stepIncrementsCounter[0])));
+    memset( stepIncrementsCounter, 0, numSteps*sizeof(stepIncrementsCounter[0]) );
+
+    // generate GLCM for each step
+    for( rowLoop=0; rowLoop<srcImageSize.height; rowLoop++, lineOffset+=srcImageStep )
+    {
+        for( colLoop=0; colLoop<srcImageSize.width; colLoop++ )
+        {
+            int pixelValue1 = destGLCM->forwardLookupTable[srcImageData[lineOffset + colLoop]];
+
+            for( stepLoop=0; stepLoop<numSteps; stepLoop++ )
+            {
+                int col2, row2;
+                row2 = rowLoop + steps[stepLoop*2 + 0];
+                col2 = colLoop + steps[stepLoop*2 + 1];
+
+                if( col2>=0 && row2>=0 && col2<srcImageSize.width && row2<srcImageSize.height )
+                {
+                    int memoryStep = memorySteps[ stepLoop ];
+                    int pixelValue2 = destGLCM->forwardLookupTable[ srcImageData[ lineOffset + colLoop + memoryStep ] ];
+
+                    // maintain symmetry
+                    matrices[stepLoop][pixelValue1][pixelValue2] ++;
+                    matrices[stepLoop][pixelValue2][pixelValue1] ++;
+
+                    // incremenet counter of total number of increments
+                    stepIncrementsCounter[stepLoop] += 2;
+                }
+            }
+        }
+    }
+
+    // normalize matrices. each element is a probability of gray value i,j adjacency in direction/magnitude k
+    for( sideLoop1=0; sideLoop1<matrixSideLength; sideLoop1++ )
+    {
+        for( sideLoop2=0; sideLoop2<matrixSideLength; sideLoop2++ )
+        {
+            for( stepLoop=0; stepLoop<numSteps; stepLoop++ )
+            {
+                matrices[stepLoop][sideLoop1][sideLoop2] /= double(stepIncrementsCounter[stepLoop]);
+            }
+        }
+    }
+
+    destGLCM->matrices = matrices;
+
+    __END__;
+
+    cvFree( &stepIncrementsCounter );
+
+    if( cvGetErrStatus() < 0 )
+        cvReleaseGLCM( &destGLCM, CV_GLCM_GLCM );
+}
+
+
+CV_IMPL void
+cvCreateGLCMDescriptors( CvGLCM* destGLCM, int descriptorOptimizationType )
+{
+    CV_FUNCNAME( "cvCreateGLCMDescriptors" );
+
+    __BEGIN__;
+
+    int matrixLoop;
+
+    if( !destGLCM )
+        CV_ERROR( CV_StsNullPtr, "" );
+
+    if( !(destGLCM->matrices) )
+        CV_ERROR( CV_StsNullPtr, "Matrices are not allocated" );
+
+    CV_CALL( cvReleaseGLCM( &destGLCM, CV_GLCM_DESC ));
+
+    if( destGLCM->optimizationType != CV_GLCM_OPTIMIZATION_HISTOGRAM )
+    {
+        destGLCM->descriptorOptimizationType = destGLCM->numDescriptors = descriptorOptimizationType;
+    }
+    else
+    {
+        CV_ERROR( CV_StsBadFlag, "Histogram-based method is not implemented" );
+//      destGLCM->descriptorOptimizationType = destGLCM->numDescriptors = CV_GLCMDESC_OPTIMIZATION_HISTOGRAM;
+    }
+
+    CV_CALL( destGLCM->descriptors = (double**)
+            cvAlloc( destGLCM->numMatrices*sizeof(destGLCM->descriptors[0])));
+
+    for( matrixLoop = 0; matrixLoop < destGLCM->numMatrices; matrixLoop ++ )
+    {
+        CV_CALL( destGLCM->descriptors[ matrixLoop ] =
+                (double*)cvAlloc( destGLCM->numDescriptors*sizeof(destGLCM->descriptors[0][0])));
+        memset( destGLCM->descriptors[matrixLoop], 0, destGLCM->numDescriptors*sizeof(double) );
+
+        switch( destGLCM->descriptorOptimizationType )
+        {
+            case CV_GLCMDESC_OPTIMIZATION_ALLOWDOUBLENEST:
+                icvCreateGLCMDescriptors_AllowDoubleNest( destGLCM, matrixLoop );
+                break;
+            default:
+                CV_ERROR( CV_StsBadFlag,
+                "descriptorOptimizationType different from CV_GLCMDESC_OPTIMIZATION_ALLOWDOUBLENEST\n"
+                "is not supported" );
+            /*
+            case CV_GLCMDESC_OPTIMIZATION_ALLOWTRIPLENEST:
+                icvCreateGLCMDescriptors_AllowTripleNest( destGLCM, matrixLoop );
+                break;
+            case CV_GLCMDESC_OPTIMIZATION_HISTOGRAM:
+                if(matrixLoop < destGLCM->numMatrices>>1)
+                    icvCreateGLCMDescriptors_Histogram( destGLCM, matrixLoop);
+                    break;
+            */
+        }
+    }
+
+    __END__;
+
+    if( cvGetErrStatus() < 0 )
+        cvReleaseGLCM( &destGLCM, CV_GLCM_DESC );
+}
+
+
+static void
+icvCreateGLCMDescriptors_AllowDoubleNest( CvGLCM* destGLCM, int matrixIndex )
+{
+    int sideLoop1, sideLoop2;
+    int matrixSideLength = destGLCM->matrixSideLength;
+
+    double** matrix = destGLCM->matrices[ matrixIndex ];
+    double* descriptors = destGLCM->descriptors[ matrixIndex ];
+
+    double* marginalProbability =
+        (double*)cvAlloc( matrixSideLength * sizeof(marginalProbability[0]));
+    memset( marginalProbability, 0, matrixSideLength * sizeof(double) );
+
+    double maximumProbability = 0;
+    double marginalProbabilityEntropy = 0;
+    double correlationMean = 0, correlationStdDeviation = 0, correlationProductTerm = 0;
+
+    for( sideLoop1=0; sideLoop1<matrixSideLength; sideLoop1++ )
+    {
+        int actualSideLoop1 = destGLCM->reverseLookupTable[ sideLoop1 ];
+
+        for( sideLoop2=0; sideLoop2<matrixSideLength; sideLoop2++ )
+        {
+            double entryValue = matrix[ sideLoop1 ][ sideLoop2 ];
+
+            int actualSideLoop2 = destGLCM->reverseLookupTable[ sideLoop2 ];
+            int sideLoopDifference = actualSideLoop1 - actualSideLoop2;
+            int sideLoopDifferenceSquared = sideLoopDifference*sideLoopDifference;
+
+            marginalProbability[ sideLoop1 ] += entryValue;
+            correlationMean += actualSideLoop1*entryValue;
+
+            maximumProbability = MAX( maximumProbability, entryValue );
+
+            if( actualSideLoop2 > actualSideLoop1 )
+            {
+                descriptors[ CV_GLCMDESC_CONTRAST ] += sideLoopDifferenceSquared * entryValue;
+            }
+
+            descriptors[ CV_GLCMDESC_HOMOGENITY ] += entryValue / ( 1.0 + sideLoopDifferenceSquared );
+
+            if( entryValue > 0 )
+            {
+                descriptors[ CV_GLCMDESC_ENTROPY ] += entryValue * log( entryValue );
+            }
+
+            descriptors[ CV_GLCMDESC_ENERGY ] += entryValue*entryValue;
+        }
+
+        if( marginalProbability>0 )
+            marginalProbabilityEntropy += marginalProbability[ actualSideLoop1 ]*log(marginalProbability[ actualSideLoop1 ]);
+    }
+
+    marginalProbabilityEntropy = -marginalProbabilityEntropy;
+
+    descriptors[ CV_GLCMDESC_CONTRAST ] += descriptors[ CV_GLCMDESC_CONTRAST ];
+    descriptors[ CV_GLCMDESC_ENTROPY ] = -descriptors[ CV_GLCMDESC_ENTROPY ];
+    descriptors[ CV_GLCMDESC_MAXIMUMPROBABILITY ] = maximumProbability;
+
+    double HXY = 0, HXY1 = 0, HXY2 = 0;
+
+    HXY = descriptors[ CV_GLCMDESC_ENTROPY ];
+
+    for( sideLoop1=0; sideLoop1<matrixSideLength; sideLoop1++ )
+    {
+        double sideEntryValueSum = 0;
+        int actualSideLoop1 = destGLCM->reverseLookupTable[ sideLoop1 ];
+
+        for( sideLoop2=0; sideLoop2<matrixSideLength; sideLoop2++ )
+        {
+            double entryValue = matrix[ sideLoop1 ][ sideLoop2 ];
+
+            sideEntryValueSum += entryValue;
+
+            int actualSideLoop2 = destGLCM->reverseLookupTable[ sideLoop2 ];
+
+            correlationProductTerm += (actualSideLoop1 - correlationMean) * (actualSideLoop2 - correlationMean) * entryValue;
+
+            double clusterTerm = actualSideLoop1 + actualSideLoop2 - correlationMean - correlationMean;
+
+            descriptors[ CV_GLCMDESC_CLUSTERTENDENCY ] += clusterTerm * clusterTerm * entryValue;
+            descriptors[ CV_GLCMDESC_CLUSTERSHADE ] += clusterTerm * clusterTerm * clusterTerm * entryValue;
+
+            double HXYValue = marginalProbability[ actualSideLoop1 ] * marginalProbability[ actualSideLoop2 ];
+            if( HXYValue>0 )
+            {
+                double HXYValueLog = log( HXYValue );
+                HXY1 += entryValue * HXYValueLog;
+                HXY2 += HXYValue * HXYValueLog;
+            }
+        }
+
+        correlationStdDeviation += (actualSideLoop1-correlationMean) * (actualSideLoop1-correlationMean) * sideEntryValueSum;
+    }
+
+    HXY1 =- HXY1;
+    HXY2 =- HXY2;
+
+    descriptors[ CV_GLCMDESC_CORRELATIONINFO1 ] = ( HXY - HXY1 ) / ( correlationMean );
+    descriptors[ CV_GLCMDESC_CORRELATIONINFO2 ] = sqrt( 1.0 - exp( -2.0 * (HXY2 - HXY ) ) );
+
+    correlationStdDeviation = sqrt( correlationStdDeviation );
+
+    descriptors[ CV_GLCMDESC_CORRELATION ] = correlationProductTerm / (correlationStdDeviation*correlationStdDeviation );
+
+    delete [] marginalProbability;
+}
+
+
+CV_IMPL double cvGetGLCMDescriptor( CvGLCM* GLCM, int step, int descriptor )
+{
+    double value = DBL_MAX;
+
+    CV_FUNCNAME( "cvGetGLCMDescriptor" );
+
+    __BEGIN__;
+
+    if( !GLCM )
+        CV_ERROR( CV_StsNullPtr, "" );
+
+    if( !(GLCM->descriptors) )
+        CV_ERROR( CV_StsNullPtr, "" );
+
+    if( (unsigned)step >= (unsigned)(GLCM->numMatrices))
+        CV_ERROR( CV_StsOutOfRange, "step is not in 0 .. GLCM->numMatrices - 1" );
+
+    if( (unsigned)descriptor >= (unsigned)(GLCM->numDescriptors))
+        CV_ERROR( CV_StsOutOfRange, "descriptor is not in 0 .. GLCM->numDescriptors - 1" );
+
+    value = GLCM->descriptors[step][descriptor];
+
+    __END__;
+
+    return value;
+}
+
+
+CV_IMPL void
+cvGetGLCMDescriptorStatistics( CvGLCM* GLCM, int descriptor,
+                               double* _average, double* _standardDeviation )
+{
+    CV_FUNCNAME( "cvGetGLCMDescriptorStatistics" );
+
+    if( _average )
+        *_average = DBL_MAX;
+
+    if( _standardDeviation )
+        *_standardDeviation = DBL_MAX;
+
+    __BEGIN__;
+
+    int matrixLoop, numMatrices;
+    double average = 0, squareSum = 0;
+
+    if( !GLCM )
+        CV_ERROR( CV_StsNullPtr, "" );
+
+    if( !(GLCM->descriptors))
+        CV_ERROR( CV_StsNullPtr, "Descriptors are not calculated" );
+
+    if( (unsigned)descriptor >= (unsigned)(GLCM->numDescriptors) )
+        CV_ERROR( CV_StsOutOfRange, "Descriptor index is out of range" );
+
+    numMatrices = GLCM->numMatrices;
+
+    for( matrixLoop = 0; matrixLoop < numMatrices; matrixLoop++ )
+    {
+        double temp = GLCM->descriptors[ matrixLoop ][ descriptor ];
+        average += temp;
+        squareSum += temp*temp;
+    }
+
+    average /= numMatrices;
+
+    if( _average )
+        *_average = average;
+
+    if( _standardDeviation )
+        *_standardDeviation = sqrt( (squareSum - average*average*numMatrices)/(numMatrices-1));
+
+    __END__;
+}
+
+
+CV_IMPL IplImage*
+cvCreateGLCMImage( CvGLCM* GLCM, int step )
+{
+    IplImage* dest = 0;
+
+    CV_FUNCNAME( "cvCreateGLCMImage" );
+
+    __BEGIN__;
+
+    float* destData;
+    int sideLoop1, sideLoop2;
+
+    if( !GLCM )
+        CV_ERROR( CV_StsNullPtr, "" );
+
+    if( !(GLCM->matrices) )
+        CV_ERROR( CV_StsNullPtr, "Matrices are not allocated" );
+
+    if( (unsigned)step >= (unsigned)(GLCM->numMatrices) )
+        CV_ERROR( CV_StsOutOfRange, "The step index is out of range" );
+
+    dest = cvCreateImage( cvSize( GLCM->matrixSideLength, GLCM->matrixSideLength ), IPL_DEPTH_32F, 1 );
+    destData = (float*)(dest->imageData);
+
+    for( sideLoop1 = 0; sideLoop1 < GLCM->matrixSideLength;
+                        sideLoop1++, (float*&)destData += dest->widthStep )
+    {
+        for( sideLoop2=0; sideLoop2 < GLCM->matrixSideLength; sideLoop2++ )
+        {
+            double matrixValue = GLCM->matrices[step][sideLoop1][sideLoop2];
+            destData[ sideLoop2 ] = (float)matrixValue;
+        }
+    }
+
+    __END__;
+
+    if( cvGetErrStatus() < 0 )
+        cvReleaseImage( &dest );
+
+    return dest;
+}
+