灰度共生矩阵代码实现-CSDN博客

本文链接：https://blog.csdn.net/breeze5428/article/details/34086425

注：理论部分从百度百科复制。若引用理论部分，请标明出自百度百科。

概念

由于纹理是由灰度分布在空间位置上反复出现而形成的，因而在图像空间中相隔某距离的两象素之间会存在一定的灰度关系，即图像中灰度的空间相关特性。灰度共生矩阵就是一种通过研究灰度的空间相关特性来描述纹理的常用方法。

灰度共生矩阵的生成

灰度直方图是对图像上单个象素具有某个灰度进行统计的结果，而灰度共生矩阵是对图像上保持某距离的两象素分别具有某灰度的状况进行统计得到的。

取图像(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度扫描。

这样，两个象素灰度级同时发生的概率，就将（x，y）的空间坐标转化为“灰度对” （g1，g2）的描述，形成了灰度共生矩阵。

实验中对灰度共生矩阵进行了如下的归一化：

灰度共生矩阵的特征

直觉上来说，如果图像的是由具有相似灰度值的像素块构成，则灰度共生矩阵的对角元素会有比较大的值；如果图像像素灰度值在局部有变化，那么偏离对角线的元素会有比较大的值。

通常可以用一些标量来表征灰度共生矩阵的特征，令 G表示灰度共生矩阵常用的特征有

ASM 能量（angular second moment)

也即每个矩阵元素的平方和。

如果灰度共生矩阵中的值集中在某一块（比如对连续灰度值图像，值集中在对角线；对结构化的图像，值集中在偏离对角线的位置），则ASM有较大值，若 G中的值分布较均匀（如噪声严重的图像），则ASM有较小的值。

能量是灰度共生矩阵元素值的平方和，所以也称能量，反映了图像灰度分布均匀程度和纹理粗细度。如果共生矩阵的所有值均相等，则ASM值小；相反，如果其中一些值大而其它值小，则ASM值大。当共生矩阵中元素集中分布时，此时ASM值大。ASM值大表明一种较均一和规则变化的纹理模式。

对比度（contrast）

如果偏离对角线的元素有较大值，即图像亮度值变化很快，则CON会有较大取值，这也符合对比度的定义。其中。反映了图像的清晰度和纹理沟纹深浅的程度。纹理沟纹越深，其对比度越大，视觉效果越清晰；反之，对比度小，则沟纹浅，效果模糊。灰度差即对比度大的象素对越多，这个值越大。灰度公生矩阵中远离对角线的元素值越大，CON越大。

IDM（inverse different moment）

如果灰度共生矩阵对角元素有较大值，IDM就会取较大的值。因此连续灰度的图像会有较大IDM值。

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

熵（entropy）

若灰度共生矩阵值分布均匀，也即图像近于随机或噪声很大，熵会有较大值。

熵是图像所具有的信息量的度量，纹理信息也属于图像的信息，是一个随机性的度量，当共生矩阵中所有元素有最大的随机性、空间共生矩阵中所有值几乎相等时，共生矩阵中元素分散分布时，熵较大。它表示了图像中纹理的非均匀程度或复杂程度。

自相关（correlation）

其中

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

最后，可以用一个向量将以上特征综合在一起。例如，当距离差分值（a,b）取四种值的时候，可以综合得到向量：

h=[ASM1, CON1, IDM1, ENT1, COR1, ..., ASM4, CON4, IDM4, ENT4, COR4]

综合后的向量就可以看做是对图像纹理的一种描述，可以进一步用来分类、识别、检索等。

代码实现

以下提取共生矩阵特征的C++代码，需要opencv才可以运行。

</pre><pre name="code" class="cpp">#define GLCM_DIS 1  //灰度共生矩阵的统计距离
#define GLCM_CLASS 64 //计算灰度共生矩阵的图像灰度值等级化
#define GLCM_ANGLE_HORIZATION 0  //水平
#define GLCM_ANGLE_VERTICAL   1	 //垂直
#define GLCM_ANGLE_DIGONAL_45    2  //对角
#define GLCM_ANGLE_DIGONAL_135    3//对角
#include <opencv/cv.h>
#include <opencv/highgui.h>
static int calGLCM(IplImage* bWavelet,int angleDirection,std::vector<float>&featureVector)
{
	int i,j;
	int width,height;
	if(NULL == bWavelet)
		return 1;

	width = bWavelet->width;
	height = bWavelet->height;
	float * glcm = new float[GLCM_CLASS * GLCM_CLASS];
	int * histImage = new int[width * height];
	int  normalized_factor=0;
	if(NULL == glcm || NULL == histImage)
		return 2;

	//灰度等级化---分GLCM_CLASS个等级
	uchar *data =(uchar*) bWavelet->imageData;
	for(i = 0;i < height;i++){
		for(j = 0;j < width;j++){
			histImage[i * width + j] = (int)(data[bWavelet->widthStep * i + j] * GLCM_CLASS / 256);
		}
	}

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

	//计算灰度共生矩阵
	int k,l;//（l，k）是点对
	//水平方向
	if(angleDirection == GLCM_ANGLE_HORIZATION)
	{
		for (i = 0;i < height;i++)
		{
			for (j = 0;j < width;j++)
			{
				if(j + GLCM_DIS < width)
				{
					l = histImage[i * width + j];//当前象元
					k = histImage[i * width + j + GLCM_DIS];
					glcm[l * GLCM_CLASS + k]++;
					++normalized_factor;
				}
			}
		}
	}
	//垂直方向
	else if(angleDirection == GLCM_ANGLE_VERTICAL)
	{
		for (i = 0;i < height;i++)
		{
			for (j = 0;j < width;j++)
			{
				if(i + GLCM_DIS < height)
				{
					l = histImage[i * width + j];
					k = histImage[(i + GLCM_DIS) * width + j];
					glcm[l * GLCM_CLASS + k]++;
					++normalized_factor;
				}
			}
		}
	}
	//45度对角方向
		else if(angleDirection == GLCM_ANGLE_DIGONAL_45)
			{
				for (i = 0;i < height;i++)
				{
					for (j = 0;j < width;j++)
					{
						if(j - GLCM_DIS >= 0  && i + GLCM_DIS < height)
						{
							l = histImage[i * width + j];
							k = histImage[(i + GLCM_DIS) * width + j - GLCM_DIS];
							glcm[l * GLCM_CLASS + k]++;
							++normalized_factor;
						}
					}
				}
			}

	//135度对角方向
	else if(angleDirection == GLCM_ANGLE_DIGONAL_135)
	{
		for (i = 0;i < height;i++)
		{
			for (j = 0;j < width;j++)
			{
				if(j + GLCM_DIS < width && i + GLCM_DIS < height)
				{
					l = histImage[i * width + j];
					k = histImage[(i + GLCM_DIS) * width + j + GLCM_DIS];
					glcm[l * GLCM_CLASS + k]++;
					++normalized_factor;
				}
			}
		}
	}

	for (i = 0;i < GLCM_CLASS;i++)
		for (j = 0;j < GLCM_CLASS;j++)
			glcm[i * GLCM_CLASS + j] = glcm[i * GLCM_CLASS + j]/normalized_factor;

	/************************计算特征值*********************/
	float u_i=0,u_j=0,s_i_2=0,s_j_2=0;
	for (i = 0;i < GLCM_CLASS;i++)
	{
		for (j = 0;j < GLCM_CLASS;j++)
		{
			u_i+=i*glcm[i * GLCM_CLASS + j];
			u_j+=j*glcm[i * GLCM_CLASS + j];
		}
	}

	for (i = 0;i < GLCM_CLASS;i++)
	{
		for (j = 0;j < GLCM_CLASS;j++)
		{
			s_i_2+=(i-u_i)*(i-u_i)*glcm[i * GLCM_CLASS + j];//s_i的平方
			s_j_2+=(j-u_j)*(j-u_j)*glcm[i * GLCM_CLASS + j];
		}
	}

	float entropy = 0,energy = 0,contrast = 0,IDM = 0,correlation=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];
			//一致性
			IDM += 1.0 / (1 + (i - j) * (i - j)) * glcm[i * GLCM_CLASS + j];
			//自相关性
			correlation+=(i*j*glcm[i * GLCM_CLASS + j]-u_i*u_j)/(sqrt(s_i_2)*sqrt(s_j_2)+0.000001);
		}
	}


	//返回特征值
	featureVector.push_back(entropy);
	featureVector.push_back(energy);
	featureVector.push_back(contrast);
	featureVector.push_back(IDM);
	featureVector.push_back(correlation);
	delete[] glcm;
	delete[] histImage;
	return 0;
}

参考：

http://baike.baidu.com/view/4293543.htm?fr=aladdin

http://blog.csdn.net/cxf7394373/article/details/6988229