テクスチャの特徴:Greey Level Co-occurrence Matrix


以下はGLACMを勉強している時に発生した問題といい資料です.後で調べるのに便利です.
opencvのlegacyライブラリにはGLACM(Gree Level Co-occurrence Matrix'es)の方法が含まれています.texture.cppファイルにあります.
しかし、直接opencvの中のGLACMを使って、createの時、メモリエラーの問題を申し込みます.http://answers.opencv.org/question/41/allocation-error-with-cvcreateglcm-method/
解決案は次の文章を参考にしてこの問題を解決しました.
http://www.cnblogs.com/xiangshancuizhu/archive/2011/07/15/2107656.html
完全なGCM例:
http://download.csdn.net/detail/hwb506/7466757
いくつかの参照文書:
http://download.csdn.net/detail/hwb506/7466817
ホームページの参考:
http://blog.csdn.net/cxf7394373/article/details/6988229
http://blog.sina.com.cn/s/blog_66 e 58 dcd 0100 lqgn.
http://blog.163.com/wujiaxing009@126/blog/static/71988399202344846222/
http://www.fp.ucalgary.ca/mhallbey/examples.htm
http://matlab.izmiran.ru/help/toolbox/images/enhanc15.html
http://www.mathworks.co.jp/jp/help/images/ref/graycomatrix.html
opencvはコードを持っています
/*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 [[email protected]]

\****************************************************************************************/

#include "precomp.hpp"

#include 
#include 

#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 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; stepLoopforwardLookupTable[srcImageData[lineOffset + colLoop]];

            for( stepLoop=0; stepLoop=0 && row2>=0 && col2forwardLookupTable[ 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; sideLoop1matrices = 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
" "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; sideLoop1reverseLookupTable[ sideLoop1 ]; for( sideLoop2=0; sideLoop2reverseLookupTable[ 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[ actualSideLoop1 ] > 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; sideLoop1reverseLookupTable[ sideLoop1 ]; for( sideLoop2=0; sideLoop2reverseLookupTable[ 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; }