灰度共生矩阵生成分析及C++实现

算法目的

通过创建灰度共生矩阵GLCM来表示纹理的空间分布情况。

算法原理

  1. 选择不同方向(0° , 45°, 90°, 135°)

  2. 对原图灰度值进行量化，如8级、16级等

  3. 遍历所有像元，统计特定步长下不同灰度级出现的次数。

    image.png

  4. 输出获得的灰度共生矩阵。

分析过程

在MFC工程中添加GDAL库，并在程序中引用并注册。

1
2
3
4
#include "gdal/include/gdal_priv.h"
#include "gdal/include/gdalwarper.h"
#include "gdal/include/ogrsf_frmts.h"
#pragma comment(lib, "gdal/lib/gdal.lib")

初始化对话框时，注册驱动器并禁用UTF8编码。

1
2
GDALAllRegister();
CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");

获取图像的灰度信息，这里创建了一维数组pData用于存储数据。

1
2
3
4
5
6
7
GDALDriver* poDriver=GetGDALDriverManager()->GetDriverByName("GTIFF");
GDALDataset* poDataset = (GDALDataset*)GDALOpen(inPath, GA_ReadOnly);
GDALRasterBand* poBand = poDataset->GetRasterBand(band+1);
int xSize = poDataset->GetRasterXSize();
int ySize = poDataset->GetRasterYSize();
int* pData = new int[xSize * ySize];
poBand->RasterIO(GF_Read, 0, 0, xSize, ySize, pData, xSize, ySize, GDT_Int32, 0, 0);

对原图进行量化，灰度级数以gray表示；同时根据灰度级数创建最终输出的灰度共生矩阵。

1
2
3
4
5
6
7
8
9
10
11
12
int max = pData[0];
for (int i = 0; i < xSize * ySize; i++) {	
if (pData[i] > max) max = pData[i];
}
int* pDataClass = new int[xSize * ySize];
for (int i = 0; i < xSize * ySize; i++) {
	pDataClass[i] = (int)(pData[i] / (double)(max) * gray);
	if (max == pData[i]) pDataClass[i] = gray - 1;
}

int* glcm = new int[gray * gray];
for (int i = 0; i < gray * gray; i++) glcm[i] = 0;

根据不同的方向，进行灰度共生矩阵的统计。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
int p1 = 0, p2 = 0;
if (strcmp("0°", ori) == 0) {
	for (int i = 0; i < ySize; i++) {
		for (int j = 0; j < xSize - step; j++) {
			p1 = pDataClass[i * xSize + j];
			p2 = pDataClass[i * xSize + j + step];
			glcm[p1 * gray + p2]++;
		}
	}
}
else if (strcmp("45°", ori) == 0) {
	for (int i = step; i < ySize; i++) {
		for (int j = 0; j < xSize - step; j++) {
			p1 = pDataClass[i * xSize + j];
			p2 = pDataClass[(i - step) * xSize + j + step];
			glcm[p1 * gray + p2]++;
		}
	}
}
else if (strcmp("90°", ori) == 0) {
	for (int i = step; i < ySize; i++) {
		for (int j = 0; j < xSize; j++) {
			p1 = pDataClass[i * xSize + j];
			p2 = pDataClass[(i - step) * xSize + j];
			glcm[p1 * gray + p2]++;
		}
	}
}
else if (strcmp("135°", ori) == 0) {
	for (int i = step; i < ySize; i++) {
		for (int j = step; j < xSize; j++) {
			p1 = pDataClass[i * xSize + j];
			p2 = pDataClass[(i - step) * xSize + j - step];
			glcm[p1 * gray + p2]++;
		}
	}
}

处理完毕，写进文本文档里输出，完成任务。

1
2
3
4
5
6
7
8
std::ofstream outfile;
outfile.open(outPath);
for (int i = 0; i < gray; i++) {
	for (int j = 0; j < gray; j++) {
		outfile << glcm[i * gray + j] << '\t';
	}
  outfile << '\n';
}