图像特征算子--之灰度共生矩阵原理分析与实现（八）
[email protected]
灰度共生矩阵最早由Robert M.在提出，早期称为灰度空间依赖矩阵（Gray-Tone Spatial-Dependence Matrices），在其论文里，根据灰度空间依赖矩阵可以计算28种纹理特征，详细内容可以参考：Textural Features for Image Classification。由于纹理是由灰度分布在空间位置上反复出现而形成的，因而在图像空间中相隔某距离的两象素之间会存在一定的灰度关系，即图像中灰度的空间相关特性。灰度共生矩阵就是一种通过研究灰度的空间相关特性来描述纹理的常用方法。
灰度共生矩阵定义为像素对的联合分布概率，是一个对称矩阵，它不仅反映图像灰度在相邻的方向、相邻间隔、变化幅度的综合信息，但也反映了相同的灰度级像素之间的位置分布特征，是计算纹理特征的基础。 设f(x,y)为一幅数字图像，其大小为M×N，灰度级别为Ng，则满足一定空间关系的灰度共生矩阵为：
其中#(x)表示集合x中的元素个数，显然P为Ng×Ng的矩阵，若(x1,y1)与(x2,y2)间距离为d,两者与坐标横轴的夹角为θ，则可以得到各种间距及角度的灰度共生矩阵(i,j,d,θ)。其中元素(i,j)的值表示一个灰度为i，另一个灰度为j的两个相距为d的像素对在角的方向上出现的次数。经常，我们在使用中求0°、45°、90°和135°四个方向的共生矩阵，公式如下：
为了更进一步的明白灰度共生矩阵的计算方式，下面通过一个例子来说明，如下图是一个4×4的图像块，将灰度量化到0~3及位置关系：
图像块 位置表示
其灰度矩阵为：
依次我们可以求得四个方向的矩阵为（这里距离d=1）：
这里简单已0°处的#（0,0）的计算做下说明：对于上图的图像块，存在(0,0)关系的有第一行的(1,1)、(1,2)和(1,2)、(1,1)，第二行有(2,1)、(2,2)和(2,2)、(2,1)，则其值为4；同理，对于#(1,0)有(1,3)、(1,2)和（2,,3)、(2,2)，则其值为2。
通过上面的介绍，我们计算得到共生矩阵之后，往往不是直接应用计算的灰度共生矩阵，而是在此基础上计算纹理特征量，我们经常用反差、能量、熵、相关性等特征量来表示纹理特征。
1）ASM （angular second moment)特征（或称能量特征）
其计算公式如下：
Asm是灰度共生矩阵元素值的平方和，所以也称能量，反映了图像灰度分布均匀程度和纹理粗细度。如果共生矩阵的所有值均相等，则ASM值小；相反，如果其中一些值大而其它值小，则ASM值大。当共生矩阵中元素集中分布时，此时ASM值大。ASM值大表明一种较均一和规则变化的纹理模式。
2）对比度（Contrast）
其计算公式如下：
对比度度量矩阵的值是如何分布和图像中局部变化的多少，反应了图像的清晰度和纹理的沟纹深浅。纹理的沟纹越深，反差越大，效果清晰；反之，对比值小，则沟纹浅，效果模糊。
3）熵（entropy）
其计算公式如下：
熵是图像包含信息量的随机性度量，纹理信息也属于图像的信息。当共生矩阵中所有元素有最大的随机性、空间共生矩阵中所有值几乎相等时，共生矩阵中元素分散分布时，熵较大。因此熵值表明了图像灰度分布的复杂程度，熵值越大，图像越复杂。
4）自相关（correlation）
其计算公式为：
相关性也称为同质性，用来度量图像的灰度级在行或列方向上的相似程度，因此值的大小反应了局部灰度相关性，值越大，相关性也越大。自相关反应了图像纹理的一致性。如果图像中有水平方向纹理，则水平方向矩阵的COR大于其余矩阵的COR值。它度量空间灰度共生矩阵元素在行或列方向上的相似程度，因此，相关值大小反映了图像中局部灰度相关性。当矩阵元素值均匀相等时，相关值就大;相反，如果矩阵像元值相差很大则相关值小。
5）逆差矩（IDM：Inverse Difference Moment）
其计算公式为：
逆差矩反映图像纹理的同质性，度量图像纹理局部变化的多少。其值大则说明图像纹理的不同区域间缺少变化，局部非常均匀。如果灰度共生矩阵对角元素有较大值，IDM就会取较大的值。因此连续灰度的图像会有较大IDM值。
有关灰度共生矩阵更多的特征值计算可以参考问候给出的参考资料，里面介绍多达20多种与灰度矩阵相关特征、如方差、方差和、熵差等。
如果以上特征用来作为图像分类、识别等，可以将上述计算的特征组合起来形成一个向量，如，当距离差分值（a,b）取四种值的时候，可以综合得到向量：
Feature=[ASM1, CON1, IDM1, ENT1, COR1, ..., ASM4, CON4, IDM4, ENT4, COR4]
综合后的向量就可以看做是对图像纹理的一种描述，可以进一步用来分类、识别、检索等。
接着，来进一步介绍灰度共生矩阵的实现（不一一描述了，不懂的差异参考资料）：
enum GLCM_DIRECTION
{
GLCM_ANGLE_HORIZONTAL,
GLCM_ANGLE_VERTICAL,
GLCM_ANGLE_DIGONAL45,
GLCM_ANGLE_DIGONAL135
};
// binImage：灰度图
// iDir:方向, 1--0°, 2--45°,3--90°,4--135°
// 归一化系数：0°（2Ny(Nx-D）, 45°(2(Ny-D)(Nx-D), 90°(2NX(Ny-D), 135°(2(Ny-D)(Nx-D))
// iDist: 1,2,3,4....
// GRAY_LEVEL:灰度层级，建议取值8、16、32、64
void CalcGLCM(IplImage *binImage,vector<double>& features, int iDist/* =1 */, GLCM_DIRECTION glcmDir/* =GLCM_ANGLE_HORIZONTAL */, int GRAY_LEVEL /* = 16 */)
{
int i, j;
int height, width, widthStep;
double factor = 0; // 归一化因子
if ( binImage == NULL)
return;
height = binImage->height;
width = binImage->width;
widthStep = binImage->widthStep;
int *glcm = new int[GRAY_LEVEL*GRAY_LEVEL];
double *glcmNorm = new double[GRAY_LEVEL*GRAY_LEVEL]; // 归一化的值
int *histImage = new int[width*height];
ZeroMemory(glcm, sizeof(int)*GRAY_LEVEL*GRAY_LEVEL);
ZeroMemory(glcmNorm, sizeof(double)*GRAY_LEVEL*GRAY_LEVEL);
ZeroMemory(histImage, sizeof(int)*width*height);
// 灰度量化到0-GRAY_LEVEL
uchar* data = (uchar*)binImage->imageData;
for ( i=0; i<height; i++)
{
for ( j=0; j<width; j++)
{
histImage[i*width+j] = (int)(data[i*widthStep+j]*GRAY_LEVEL/256);
}
}
// 计算灰度共生矩阵
int k, l;
if ( glcmDir == GLCM_ANGLE_HORIZONTAL) // 水平方向（0°）
{
factor = 2*height*(width-iDist);
for ( i=0; i<height; i++)
{
for ( j=0; j<width; j++)
{
l = histImage[i*width+j];
if ( j+iDist>=0 && j+iDist<width)
{
k = histImage[i*width+j+iDist];
glcm[l*GRAY_LEVEL+k]++;
}
if ( j-iDist>=0 && j-iDist<width)
{
k = histImage[i*width+j-iDist];
glcm[l*GRAY_LEVEL+k]++;
}
}// for j
} // for i
}
else if ( glcmDir == GLCM_ANGLE_DIGONAL45) // 对角方向（45°）
{
factor = 2*(height-iDist)*(width-iDist);
for ( i=0; i<height; i++)
{
for ( j=0; j<width; j++)
{
l = histImage[i*width+j];
if ( j+iDist>=0 && j+iDist<width && i-iDist>=0 && i-iDist<height)
{
k = histImage[(i+iDist)*width+j+iDist];
glcm[l*GRAY_LEVEL+k]++;
}
if ( j-iDist>=0 && j-iDist<width && i+iDist>=0 && i+iDist<height )
{
k = histImage[(i-iDist)*width+j-iDist];
glcm[l*GRAY_LEVEL+k]++;
}
}// for j
} // for i
}
else if ( glcmDir == GLCM_ANGLE_VERTICAL) // 垂直方向 90°
{
factor = 2*width*(height-iDist);
for ( i=0; i<height; i++)
{
for ( j=0; j<width; j++)
{
l = histImage[i*width+j];
if ( i+iDist>=0 && i+iDist<height)
{
k = histImage[(i+iDist)*width+j];
glcm[l*GRAY_LEVEL+k]++;
}
if ( i-iDist>=0 && i-iDist<height)
{
k = histImage[(i-iDist)*width+j];
glcm[l*GRAY_LEVEL+k]++;
}
}// for j
} // for i
}
else if ( glcmDir == GLCM_ANGLE_DIGONAL135) // 对角方向（135°）
{
factor = 2*(width-iDist)*(height-iDist);
for ( i=0; i<height; i++)
{
for ( j=0; j<width; j++)
{
l = histImage[i*width+j];
if ( j+iDist>=0 && j+iDist<width && i+iDist>=0 && i+iDist<height)
{
k = histImage[(i+iDist)*width+j+iDist];
glcm[l*GRAY_LEVEL+k]++;
}
if ( j-iDist>=0 && j-iDist<width && i-iDist>=0 && i-iDist<height )
{
k = histImage[(i-iDist)*width+j-iDist];
glcm[l*GRAY_LEVEL+k]++;
}
}// for j
} // for i
}
// 归一化
for ( i=0; i<GRAY_LEVEL; i++)
{
for ( j=0; j<GRAY_LEVEL;j++)
{
glcmNorm[i*GRAY_LEVEL+j] = (double)glcm[i*GRAY_LEVEL+j]/factor;
}
}
// 计算特征值
double entropy = 0, energy = 0, contrast = 0, homogenity = 0;
for ( i=0; i<GRAY_LEVEL; i++)
{
for ( j=0; j<GRAY_LEVEL; j++)
{
// 熵
if ( glcmNorm[i*GRAY_LEVEL+j]>0)
entropy -= glcmNorm[i*GRAY_LEVEL+j]* log10(double(glcmNorm[i*GRAY_LEVEL+j]));
// 能量
energy += glcmNorm[i*GRAY_LEVEL+j]*glcmNorm[i*GRAY_LEVEL+j];
// 对比度
contrast += (i-j)*(i-j)*glcmNorm[i*GRAY_LEVEL+j];
// 一致性
homogenity += 1.0/(1+(i-j)*(i-j))*glcmNorm[i*GRAY_LEVEL+j];
}
}
// 组成特征向量
features.clear();
features.push_back(entropy);
features.push_back(energy);
features.push_back(contrast);
features.push_back(homogenity);
delete[]glcm;
delete[]glcmNorm;
delete[]histImage;
}
在上面的算法实现中，对灰度共生矩阵进行了归一化（主要是考虑到不归一化时的计算结果值会比较大），简单测试了下，觉得计算得到的数据值应该是正常的，如图：
参考资料：
1、灰度共生矩阵
2、Textural Features for Image Classification.