灰度共生矩阵

灰度直方图是对整张图像所有像素点的灰度值进行统计的结果。

灰度共生矩阵是对图像上保持某种距离的两像素点的灰度状况进行统计的结果。

取图像 (N×N) 中任意一点 （x，y）及偏离它的另一点 （x+a，y+b），设该点对的灰度值为 （g1，g2）。令点（x，y） 在整个画面上移动，则会得到各种 （g1，g2）值，设灰度值的级数为 k，则（g1，g2） 的组合共有 k 的平方种。对于整个画面，统计出每一种 （g1，g2）值出现的次数，然后排列成一个方阵，再用（g1，g2） 出现的总次数将它们归一化为出现的概率P（g1，g2）这样的方阵称为灰度共生矩阵。距离差分值（a，b） 取不同的数值组合，可以得到不同情况下的联合概率矩阵。（a，b） 取值要根据纹理周期分布的特性来选择，对于较细的纹理，选取（1，0）、（1，1）、（2，0）等小的差分值。

当 a=1，b=0时，像素对是水平的，即0度扫描；

当 a=0，b=1 时，像素对是垂直的，即90度扫描；

当 a=1，b=1时，像素对是右对角线的，即45度扫描；

当 a=-1，b=1时，像素对是左对角线，即135度扫描。

#include "cbir.h"

#define GLCM_DISTEN 3  //偏离距离 
#define GLCM_CLASS 8 //灰度降级  

/*
功能：
    求图像4个方向的灰度共生矩阵，并计算4个特征值，保存于featherVector[]
参数：
    Grey_Image：灰度图像
    angleDirection：方向
    featherVector：保存特征值
*/

int calGLCM(Mat Grey_Image, int angleDirection, double* featureVector)
{
    int i, j;
    int width, height;

    width = Grey_Image.cols;
    height = Grey_Image.rows;

    int * glcm = new int[GLCM_CLASS * GLCM_CLASS];//灰度共生矩 8 * 8 一维数组
    int * histImage = new int[width * height];//降级后的图像 300 * 200 一维数组

    if (NULL == glcm || NULL == histImage)
        return 2;

    //降低图像灰度等级，把图像256个灰度级按顺序分为8个区间  
    uchar *data = (uchar*)Grey_Image.data;
    for (i = 0; i < height; i++) {
        for (j = 0; j < width; j++) {
            histImage[i * width + j] = (int)(data[width * i + j] * GLCM_CLASS / 256);// 除32取整 = 降级后的灰度
        }
    }

    //初始化共生矩阵  
    for (i = 0; i < GLCM_CLASS; i++)
        for (j = 0; j < GLCM_CLASS; j++)
            glcm[i * GLCM_CLASS + j] = 0;// 8 x 8 置零

    //计算灰度共生矩阵  
    int k, l;

    // 0°水平
    if (angleDirection == 0)
    {
        // 遍历降级后的图像 300 * 200
        for (i = 0; i < height; i++)
        {
            for (j = 0; j < width; j++)
            {
                l = histImage[i * width + j];//取图像点（i,j）像素值

                if (j + GLCM_DISTEN < width)
                {
                    k = histImage[i * width + j + GLCM_DISTEN];//取图像点(i,j+d)像素值
                    glcm[l * GLCM_CLASS + k]++;// 将灰度共生矩中（value(i,j),value(i,j+d)）值加1
                }
            }
        }
    }

    // 90°垂直
    else if (angleDirection == 90)
    {
        for (i = 0; i < height; i++)
        {
            for (j = 0; j < width; j++)
            {
                l = histImage[i * width + j];//取图像点（i,j）像素值

                if (i + GLCM_DISTEN < height)
                {
                    k = histImage[(i + GLCM_DISTEN) * width + j];//取图像点（i+d,j）像素值
                    glcm[l * GLCM_CLASS + k]++;// 将灰度共生矩中（value(i,j),value(i+d,j)）值加1
                }
            }
        }
    }

    // 45°
    else if (angleDirection == 45)
    {
        for (i = 0; i < height; i++)
        {
            for (j = 0; j < width; j++)
            {
                l = histImage[i * width + j];//取图像点（i,j）像素值

                if (j + GLCM_DISTEN < width && i + GLCM_DISTEN < height)
                {
                    k = histImage[(i + GLCM_DISTEN) * width + j + GLCM_DISTEN];//取图像点（i+d,j+d）像素值
                    glcm[l * GLCM_CLASS + k]++;// 将灰度共生矩中（value(i,j),value(i+d,j+d)）值加1
                }
            }
        }
    }

    // 135°
    else if (angleDirection == 135)
    {
        for (i = 0; i < height; i++)
        {
            for (j = 0; j < width; j++)
            {
                l = histImage[i * width + j];//取图像点（i,j）像素值

                if (j - GLCM_DISTEN >= 0 && i - GLCM_DISTEN >= 0)
                {
                    k = histImage[(i - GLCM_DISTEN) * width + j - GLCM_DISTEN];//取图像点（i+d,j-d）像素值
                    glcm[l * GLCM_CLASS + k]++;// 将灰度共生矩中（value(i,j),value(i+d,j-d)）值加1
                }
            }
        }
    }

    //计算特征值  
    double entropy = 0, energy = 0, contrast = 0, homogenity = 0;

    for (i = 0; i < GLCM_CLASS; i++)
    {
        for (j = 0; j < GLCM_CLASS; j++)
        {
            //熵  
            if (glcm[i * GLCM_CLASS + j] > 0)
                entropy -= glcm[i * GLCM_CLASS + j] * log10(double(glcm[i * GLCM_CLASS + j]));

            //能量  
            energy += glcm[i * GLCM_CLASS + j] * glcm[i * GLCM_CLASS + j];

            //对比度  
            contrast += (i - j) * (i - j) * glcm[i * GLCM_CLASS + j];

            //一致性  
            homogenity += 1.0 / (1 + (i - j) * (i - j)) * glcm[i * GLCM_CLASS + j];
        }
    }

    //返回特征值  
    i = 0;
    featureVector[i++] = entropy;
    featureVector[i++] = energy;
    featureVector[i++] = contrast;
    featureVector[i++] = homogenity;

    delete[] glcm;
    delete[] histImage;

    return 0;
}


/*
功能：
    纹理特征：灰度共生矩
参数：
    src:传入图像
*/

double AveDev[8];// 保存 8 个特征值

double getEigenValueByHDGSJ(Mat src)
{
    Mat dst = rgb2grey(src);//RNG -> GRAY
    imshow("1", dst);
    waitKey(1000);

    // 垂直 0°
    double GLCM_0[4] = { 0 };
    calGLCM(dst, 0, GLCM_0);
    cout << setprecision(20) << " 0°熵：" << GLCM_0[0] << endl;
    cout << setprecision(20) << " 0°能量：" << GLCM_0[1] << endl;
    cout << setprecision(20) << " 0°对比度：" << GLCM_0[2] << endl;
    cout << setprecision(20) << " 0°一致性：" << GLCM_0[3] << endl;

    double Average0 = 0.25 * (GLCM_0[0] + GLCM_0[1] + GLCM_0[2] + GLCM_0[3]);
    cout << setprecision(20) << " 0°均值：" << Average0 << endl;

    double temp0 = 1;
    for (int i = 0; i < 4; i++)
    {
        temp0 += GLCM_0[i] * GLCM_0[i];
    }
    double StaDev0 = sqrt(0.25 * temp0);
    cout << setprecision(20) << " 0°标准差：" << StaDev0 << endl;

    cout << endl;

    // 水平 90°
    double GLCM_90[4] = { 0 };
    calGLCM(dst, 90, GLCM_90);
    cout << setprecision(20) << " 90°熵：" << GLCM_90[0] << endl;
    cout << setprecision(20) << " 90°能量：" << GLCM_90[1] << endl;
    cout << setprecision(20) << " 90°对比度：" << GLCM_90[2] << endl;
    cout << setprecision(20) << " 90°一致性：" << GLCM_90[3] << endl;

    double Average90 = 0.25 * (GLCM_90[0] + GLCM_90[1] + GLCM_90[2] + GLCM_90[3]);
    cout << setprecision(20) << " 90°均值：" << Average90 << endl;

    double temp90 = 1;
    for (int i = 0; i < 4; i++)
    {
        temp90 += GLCM_90[i] * GLCM_90[i];
    }
    double StaDev90 = sqrt(0.25 * temp90);
    cout << setprecision(20) << " 90°标准差：" << StaDev90 << endl;

    cout << endl;

    // 45°
    double GLCM_45[4] = { 0 };
    calGLCM(dst, 45, GLCM_45);
    cout << setprecision(20) << " 45°熵：" << GLCM_45[0] << endl;
    cout << setprecision(20) << " 45°能量：" << GLCM_45[1] << endl;
    cout << setprecision(20) << " 45°对比度：" << GLCM_45[2] << endl;
    cout << setprecision(20) << " 45°一致性：" << GLCM_45[3] << endl;

    double Average45 = 0.25 * (GLCM_45[0] + GLCM_45[1] + GLCM_45[2] + GLCM_45[3]);
    cout << setprecision(20) << " 45°均值：" << Average45 << endl;

    double temp45 = 1;
    for (int i = 0; i < 4; i++)
    {
        temp45 += GLCM_45[i] * GLCM_45[i];
    }
    double StaDev45 = sqrt(0.25 * temp45);
    cout << setprecision(20) << " 45°标准差：" << StaDev45 << endl;

    cout << endl;

    // 135°
    double GLCM_135[4] = { 0 };
    calGLCM(dst, 135, GLCM_135);
    cout << setprecision(20) << " 135°熵：" << GLCM_135[0] << endl;
    cout << setprecision(20) << " 135°能量：" << GLCM_135[1] << endl;
    cout << setprecision(20) << " 135°对比度：" << GLCM_135[2] << endl;
    cout << setprecision(20) << " 135°一致性：" << GLCM_135[3] << endl;

    double Average135 = 0.25 * (GLCM_135[0] + GLCM_135[1] + GLCM_135[2] + GLCM_135[3]);
    cout << setprecision(20) << " 135°均值：" << Average135 << endl;

    double temp135 = 1;
    for (int i = 0; i < 4; i++)
    {
        temp135 += GLCM_135[i] * GLCM_135[i];
    }
    double StaDev135 = sqrt(0.25 * temp135);
    cout << setprecision(20) << " 135°标准差：" << StaDev135 << endl;

    cout << endl;

    AveDev[0] = Average0;
    AveDev[1] = Average90;
    AveDev[2] = Average45;
    AveDev[3] = Average135;
    AveDev[4] = StaDev0;
    AveDev[5] = StaDev90;
    AveDev[6] = StaDev45;
    AveDev[7] = StaDev135;

    return 0;

}