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

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;
}