利用灰度共生矩阵提取图像纹理特征

1. 灰度共生矩阵概念

灰度共生矩阵定义为像素对的联合概率分布，是一个对称矩阵，它不仅反映图像灰度在相邻的方向、相邻间隔、变化幅度的综合信息，也反映了相同的灰度级像素之间的位置分布特征，是计算纹理特征的基础。

在图像中任意取一点(x,y)及偏离它的一点(x+a,y+b)(其中，a、b为整数，人为定义)构成点对。设该点对的灰度值为(f1,f2)，再令点(x,y)在整幅图像上移动，则会得到不同的(f1,f2)值。

设图像的最大灰度级为L，则f1与f2的组合共有L*L种。对于整幅图像，统计出每一种(f1,f2)值出现的次数，然后排列成一个方阵，再用(f1,f2)出现的总次数将他们归一化为出现的概率P(f1,f2)，由此产生的矩阵为灰度共生矩阵。θ方向上的间隔为d的灰度共生矩阵实际上是θ方向间隔为d的灰度变化量的联合概率分布。
2. 共生矩阵的计算

公式(1)中d表示像素间隔，(k,l), (m,n)分别为原像素和偏移后的像素坐标，其中k,m为纵坐标，D为图像范围[Image Processing, Analysis, and Machine Vision (Sonka 3rd Edition2007)]

举例说明，假设原图像如图1.a所示

对1.b中蓝色字表示原像素灰度值，红字为偏移后像素灰度值。则对矩阵元素P0°,1 (0,0)表示1.a中在0°方向上（包括正和负方向）相距为1的(0,0)点对有两对，考虑正负方向的加倍效果，P0°,1 (0,0)=4。同样由于公式(1)对距离d定义的双向性，使得灰度共生矩阵为对称矩阵。

为了减小计算量，可将d定义为沿θ正方向。则(1)式变为

由1.a得到的新的灰度共生矩阵为

3. 共生矩阵计算纹理特征

能量(Energy)：是灰度共生矩阵各元素值的平方和，是对图像纹理的灰度变化稳定程度的度量，反应了图像灰度分布均匀程度和纹理粗细度。能量值大表明当前纹理是一种规则变化较为稳定的纹理。

熵(Entropy)：是图像包含信息量的随机性度量。当共生矩阵中所有值均相等或者像素值表现出最大的随机性时，熵最大；因此熵值表明了图像灰度分布的复杂程度，熵值越大，图像越复杂。

最大概率(Maximum probability)：表示图像中出现次数最多的纹理特征。

对比度(Contrast)：度量矩阵的值是如何分布和图像中局部变化的多少，反应了图像的清晰度和纹理的沟纹深浅。纹理的沟纹越深，反差越大，效果清晰；反之，对比值小，则沟纹浅，效果模糊。对公式(6)，典型的有κ=2，λ=1。

倒数差分矩(Inverse difference moment)：反映图像纹理的同质性，度量图像纹理局部变化的多少。其值大则说明图像纹理的不同区域间缺少变化，局部非常均匀。

相关性(Correlation)：自相关反应了图像纹理的一致性。如果图像中有水平方向纹理，则水平方向共生矩阵Correlation值大于其余方向共生矩阵Correlation的值。它度量空间灰度共生矩阵元素在行或列方向上的相似程度，因此，相关值大小反映了图像中局部灰度相关性。当矩阵元素值均匀相等时，相关值就大;相反，如果矩阵像元值相差很大则相关值小。

其中μx, μy为均值，σx, σy为标准差，计算公式如下
4. 算法实现

#pragma once
#include<iostream>
#include <cassert>
#include <vector>
#include <iterator>
#include <functional>
#include <algorithm>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

typedef vector<vector<int> > VecGLCM;

typedef struct _GLCMFeatures
{
    _GLCMFeatures()
        : energy(0.0)
        , entropy(0.0)
        , contrast(0.0)
        , idMoment(0.0)
    {

    }

    double energy;      // 能量 
    double entropy;     // 熵
    double contrast;    // 对比度
    double idMoment;    // 逆差分矩， inverse difference moment

} GLCMFeatures;

class GLCM
{
public:
    GLCM();
    ~GLCM();

public:
    // 枚举灰度共生矩阵的方向
    enum
    {
        GLCM_HORIZATION = 0,        // 水平
        GLCM_VERTICAL = 1,          // 垂直
        GLCM_ANGLE45 = 2,           // 45度角
        GLCM_ANGLE135 = 3           // 135度角
    };

public:
    // 计算灰度共生矩阵
//  void calGLCM(IplImage* inputImg, VecGLCM& vecGLCM, int angle);
    void calGLCM(cv::Mat &inputImg, VecGLCM &vecGLCM, int angle);

    // 计算特征值
    void getGLCMFeatures(VecGLCM& vecGLCM, GLCMFeatures& features);
public:
    // 初始化灰度共生矩阵
    void initGLCM(VecGLCM& vecGLCM, int size = 16);
    // 设置灰度划分等级，默认值为 16
    void setGrayLevel(int grayLevel) { m_grayLevel = grayLevel; }
    // 获取灰度等级
    int getGrayLevel() const { return m_grayLevel; }
private:
    // 计算水平灰度共生矩阵
    void getHorisonGLCM(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight);
    // 计算垂直灰度共生矩阵
    void getVertialGLCM(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight);
    // 计算 45 度灰度共生矩阵
    void getGLCM45(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight);
    // 计算 135 度灰度共生矩阵
    void getGLCM135(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight);

private:
    int m_grayLevel;        // 将灰度共生矩阵划分为 grayLevel 个等级

};

#include "GLCM.h"

GLCM::GLCM() : m_grayLevel(16)
{

}

GLCM::~GLCM()
{

}

//==============================================================================
// 函数名称: initGLCM
// 参数说明: vecGLCM,要进行初始化的共生矩阵,为二维方阵
//          size, 二维矩阵的大小,必须与图像划分的灰度等级相等
// 函数功能: 初始化二维矩阵
//==============================================================================

void GLCM::initGLCM(VecGLCM& vecGLCM, int size)
{
    assert(size == m_grayLevel);
    vecGLCM.resize(size);
    for (int i = 0; i < size; ++i)
    {
        vecGLCM[i].resize(size);
    }

    for (int i = 0; i < size; ++i)
    {
        for (int j = 0; j < size; ++j)
        {
            vecGLCM[i][j] = 0;
        }
    }
}

//==============================================================================
// 函数名称: getHorisonGLCM
// 参数说明: src,要进行处理的矩阵,源数据
//          dst,输出矩阵,计算后的矩阵，即要求的灰度共生矩阵
//          imgWidth, 图像宽度
//          imgHeight, 图像高度
// 函数功能: 计算水平方向的灰度共生矩阵
//==============================================================================

void GLCM::getHorisonGLCM(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight)
{
    int height = imgHeight;
    int width = imgWidth;

    for (int i = 0; i < height; ++i)
    {
        for (int j = 0; j < width - 1; ++j)
        {
            int rows = src[i][j];
            int cols = src[i][j + 1];
            dst[rows][cols]++;
        }
    }


}

//==============================================================================
// 函数名称: getVertialGLCM
// 参数说明: src,要进行处理的矩阵,源数据
//          dst,输出矩阵,计算后的矩阵，即要求的灰度共生矩阵
//          imgWidth, 图像宽度
//          imgHeight, 图像高度
// 函数功能: 计算垂直方向的灰度共生矩阵
//==============================================================================

void GLCM::getVertialGLCM(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight)
{
    int height = imgHeight;
    int width = imgWidth;
    for (int i = 0; i < height - 1; ++i)
    {
        for (int j = 0; j < width; ++j)
        {
            int rows = src[i][j];
            int cols = src[i + 1][j];
            dst[rows][cols]++;
        }
    }
}

//==============================================================================
// 函数名称: getGLCM45
// 参数说明: src,要进行处理的矩阵,源数据
//          dst,输出矩阵,计算后的矩阵，即要求的灰度共生矩阵
//          imgWidth, 图像宽度
//          imgHeight, 图像高度
// 函数功能: 计算45度的灰度共生矩阵
//==============================================================================

void GLCM::getGLCM45(VecGLCM &src, VecGLCM &dst, int imgWidth, int imgHeight)
{
    int height = imgHeight;
    int width = imgWidth;
    for (int i = 0; i < height - 1; ++i)
    {
        for (int j = 0; j < width - 1; ++j)
        {
            int rows = src[i][j];
            int cols = src[i + 1][j + 1];
            dst[rows][cols]++;
        }
    }
}


//==============================================================================
// 函数名称: getGLCM135
// 参数说明: src,要进行处理的矩阵,源数据
//          dst,输出矩阵,计算后的矩阵，即要求的灰度共生矩阵
//          imgWidth, 图像宽度
//          imgHeight, 图像高度
// 函数功能: 计算 135 度的灰度共生矩阵
//==============================================================================

void GLCM::getGLCM135(VecGLCM& src, VecGLCM& dst, int imgWidth, int imgHeight)
{
    int height = imgHeight;
    int width = imgWidth;
    for (int i = 0; i < height - 1; ++i)
    {
        for (int j = 1; j < width; ++j)
        {
            int rows = src[i][j];
            int cols = src[i + 1][j - 1];
            dst[rows][cols]++;
        }
    }
}

//==============================================================================
// 函数名称: calGLCM
// 参数说明: inputImg,要进行纹理特征计算的图像,为灰度图像
//          vecGLCM, 输出矩阵,根据灰度图像计算出的灰度共生阵
//          angle,灰度共生矩阵的方向,有水平、垂直、45度、135度四个方向
// 函数功能: 计算灰度共生矩阵
//==============================================================================

// void GLCM::calGLCM(IplImage* inputImg, VecGLCM& vecGLCM, int angle)
// {
//  assert(inputImg->nChannels == 1);
//  IplImage* src = NULL;
//  src = cvCreateImage(cvGetSize(inputImg), IPL_DEPTH_32S, inputImg->nChannels);
//  cvConvert(inputImg, src);
// 
//  int height = src->height;
//  int width = src->width;
//  int maxGrayLevel = 0;
//  // 寻找最大像素灰度最大值
//  for (int i = 0; i < height; ++i)
//  {
//      for (int j = 0; j < width; ++j)
//      {
//          int grayVal = cvGetReal2D(src, i, j);
//          if (grayVal > maxGrayLevel)
//          {
//              maxGrayLevel = grayVal;
//          }
// 
//      }
//  }// end for i
// 
//  ++maxGrayLevel;
//  VecGLCM tempVec;
//  // 初始化动态数组
//  tempVec.resize(height);
//  for (int i = 0; i < height; ++i)
//  {
//      tempVec[i].resize(width);
//  }
// 
//  if (maxGrayLevel > 16)//若灰度级数大于16，则将图像的灰度级缩小至16级，减小灰度共生矩阵的大小。
//  {
//      for (int i = 0; i < height; ++i)
//      {
//          for (int j = 0; j < width; ++j)
//          {
//              int tmpVal = cvGetReal2D(src, i, j);
//              tmpVal /= m_grayLevel;
//              tempVec[i][j] = tmpVal;
//          }
//      }
// 
//      if (angle == GLCM_HORIZATION)  // 水平方向
//          getHorisonGLCM(tempVec, vecGLCM, width, height);
//      if (angle == GLCM_VERTICAL)    // 垂直方向
//          getVertialGLCM(tempVec, vecGLCM, width, height);
//      if (angle == GLCM_ANGLE45)     // 45 度灰度共生阵
//          getGLCM45(tempVec, vecGLCM, width, height);
//      if (angle == GLCM_ANGLE135)    // 135 度灰度共生阵
//          getGLCM135(tempVec, vecGLCM, width, height);
//  }
//  else//若灰度级数小于16，则生成相应的灰度共生矩阵
//  {
//      for (int i = 0; i < height; ++i)
//      {
//          for (int j = 1; j < width; ++j)
//          {
//              int tmpVal = cvGetReal2D(src, i, j);
//              tempVec[i][j] = tmpVal;
//          }
//      }
// 
//      if (angle == GLCM_HORIZATION)  // 水平方向
//          getHorisonGLCM(tempVec, vecGLCM, width, height);
//      if (angle == GLCM_VERTICAL)    // 垂直方向
//          getVertialGLCM(tempVec, vecGLCM, width, height);
//      if (angle == GLCM_ANGLE45)     // 45 度灰度共生阵
//          getGLCM45(tempVec, vecGLCM, width, height);
//      if (angle == GLCM_ANGLE135)    // 135 度灰度共生阵
//          getGLCM135(tempVec, vecGLCM, width, height);
//  }
// 
//  cvReleaseImage(&src);
// }

void GLCM::calGLCM(cv::Mat & inputImg, VecGLCM & vecGLCM, int angle)
{
    assert(inputImg.channels() == 1);
    int height = inputImg.rows;
    int width = inputImg.cols;
    int maxGrayLevel = 0;

    // 寻找最大像素灰度最大值
    for (int i = 0 ; i < height ; ++i)
    {
        for (int j = 0 ; j < width ; ++j)
        {
            int grayVal = inputImg.at<uchar>(i, j);
            if (grayVal > maxGrayLevel)
            {
                maxGrayLevel = grayVal;
            }
        }
    }

    ++maxGrayLevel;
    VecGLCM tempVec;
    tempVec.resize(height);
    for (int i = 0 ; i < height ; ++i)
    {
        tempVec[i].resize(width);
    }

    if (maxGrayLevel > 16)
    {
        for (int i = 0 ; i < height ; ++i)
        {
            for (int j = 0 ; j < width ; ++j)
            {
                int tmpVal = inputImg.at<uchar>(i, j);
                tmpVal /= m_grayLevel;
                tempVec[i][j] = tmpVal;
            }
        }

        if (angle == GLCM_HORIZATION)
        {
            getHorisonGLCM(tempVec, vecGLCM, width, height);
        }
        else if (angle == GLCM_VERTICAL)
        {
            getVertialGLCM(tempVec, vecGLCM, width, height);
        }
        else if (angle == GLCM_ANGLE45)
        {
            getGLCM45(tempVec, vecGLCM, width, height);
        }
        else if (angle == GLCM_ANGLE135)
        {
            getGLCM135(tempVec, vecGLCM, width, height);
        }
    }
    else
    {
        for (int i = 0; i < height; ++i)
        {
            for (int j = 0; j < width; ++j)
            {
                int tmpVal = inputImg.at<uchar>(i, j);
                tempVec[i][j] = tmpVal;
            }
        }
        if (angle == GLCM_HORIZATION)
        {
            getHorisonGLCM(tempVec, vecGLCM, width, height);
        }
        else if (angle == GLCM_VERTICAL)
        {
            getVertialGLCM(tempVec, vecGLCM, width, height);
        }
        else if (angle == GLCM_ANGLE45)
        {
            getGLCM45(tempVec, vecGLCM, width, height);
        }
        else if (angle == GLCM_ANGLE135)
        {
            getGLCM135(tempVec, vecGLCM, width, height);
        }
    }
}

//==============================================================================
// 函数名称: getGLCMFeatures
// 参数说明: vecGLCM, 输入矩阵,灰度共生阵
//          features,灰度共生矩阵计算的特征值,主要包含了能量、熵、对比度、逆差分矩
// 函数功能: 根据灰度共生矩阵计算的特征值
//==============================================================================

void GLCM::getGLCMFeatures(VecGLCM& vecGLCM, GLCMFeatures& features)
{
    int total = 0;

    for (int i = 0; i < m_grayLevel; ++i)
    {
        for (int j = 0; j < m_grayLevel; ++j)
        {
            total += vecGLCM[i][j];     // 求所有图像的灰度值的和
        }
    }

    vector<vector<double> > temp;
    temp.resize(m_grayLevel);
    for (int i = 0; i < m_grayLevel; ++i)
    {
        temp[i].resize(m_grayLevel);
    }

    // 归一化
    for (int i = 0; i < m_grayLevel; ++i)
    {
        for (int j = 0; j < m_grayLevel; ++j)
        {
            temp[i][j] = (double)vecGLCM[i][j] / (double)total;
        }
    }

    for (int i = 0; i < m_grayLevel; ++i)
    {
        for (int j = 0; j < m_grayLevel; ++j)
        {
            features.energy += temp[i][j] * temp[i][j];

            if (temp[i][j] > 0)
                features.entropy -= temp[i][j] * log(temp[i][j]);               //熵     

            features.contrast += (double)(i - j)*(double)(i - j)*temp[i][j];        //对比度
            features.idMoment += temp[i][j] / (1 + (double)(i - j)*(double)(i - j));//逆差矩
        }
    }
}

#include "GLCM.h"

int main()
{

    IplImage* img = cvLoadImage("1.jpg", 0);
//  cv::Mat src = cv::imread("1.jpg", 0);
    GLCM glcm;
    VecGLCM vec;
    GLCMFeatures features;
    glcm.initGLCM(vec);
    // 水平
    glcm.calGLCM(img, vec, GLCM::GLCM_HORIZATION);
    glcm.getGLCMFeatures(vec, features);
    // 垂直
    glcm.calGLCM(img, vec, GLCM::GLCM_VERTICAL);
    glcm.getGLCMFeatures(vec, features);
    // 45 度
    glcm.calGLCM(img, vec, GLCM::GLCM_ANGLE45);
    glcm.getGLCMFeatures(vec, features);
    // 135 度
    glcm.calGLCM(img, vec, GLCM::GLCM_ANGLE135);
    glcm.getGLCMFeatures(vec, features);

    cout << "asm = " << features.energy << endl;
    cout << "eng = " << features.entropy << endl;
    cout << "Con = " << features.contrast << endl;
    cout << "Idm = " << features.idMoment << endl;
    system("pause");
    return 0;
}