算法作用

在进行非监督分类之后，往往需要进一步处理才能获得理想的分类结果。栅格图斑聚类分析就是其中一种方式。

通过栅格图斑聚类，我们最终会输出两个图像：一是像元点Data值为像元点所对应的图斑序号；二是像元点Data值为图斑面积值的图像。

算法原理

图斑的聚类需要判断像元A与左侧、上方像元的关系：

与左侧关系	与上方关系	处理方式
不同	不同	产生一个新的图斑
相同	不同	合并入左侧图斑
不同	相同	合并入上方图斑
相同	相同	判断两图斑是否属于同一图斑序号 若相同则合并进入该图斑 若不同则需要合并两个图斑

分析过程

读取数据

老一套，首先在MFC工程中添加GDAL库，并在程序中引用并注册。

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

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

cpp
复制代码
1GDALAllRegister();
2CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");

创建数据集，打开经过K-Means聚类处理后的图像，并获得条带数据。

cpp
复制代码
1GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("Gtiff");//获得Gtiff驱动
2GDALDataset* poData = (GDALDataset*)GDALOpen(inPath, GA_ReadOnly);//创建数据集获取图像数据
3int xSize = poData->GetRasterXSize();
4int ySize = poData->GetRasterYSize();
5GDALRasterBand* poBand = poData->GetRasterBand(1);//获取条带数据

同理创建用于输出的数据集，以及对应的条带数据备用。

cpp
复制代码
1GDALDataset* powData = poDriver->Create(outPath, xSize, ySize, 1, GDT_Float64, NULL);//结果图像
2GDALDataset* posData = poDriver->Create(statPath, xSize, ySize, 1, GDT_Float64, NULL);//统计图像
3GDALRasterBand* powBand = powData->GetRasterBand(1);
4GDALRasterBand* posBand = posData->GetRasterBand(1);

数据处理

通过按行读取数据的方式，获得上下两行数据：lastLine和nextLine，便于对比像元与上方像元的数值。

图像处理时需要注意第一列以及第一行的处理方式。这里在处理第一行时将lastLine初始化为-1进行相等判断，即第一行与“上方像元”不同；处理第一列时进行特判，即第一列与“左侧像元”不同。

合并图斑时注意左侧图斑与上方图斑序号不同时的情况，此时需要将左侧图斑全部合并入上方像元所在的图斑，因此需要遍历左侧像元上方全部的行列。

cpp
复制代码
1int* pDstData = new int[xSize * ySize];//用于存储像元点对应图斑序号
2for (int i = 0; i < xSize * ySize; i++) pDstData[i] = -1;
3int flag = 0;//图斑序号值
4int* lastLine = new int[xSize];
5int* nextLine = new int[xSize];
6for (int i = 0; i < xSize; i++) {
7	lastLine[i] = -1;
8	nextLine[i] = -1;
9}
10//处理结果图像
11for (int i = 0; i < ySize; i++) {
12	for (int j = 0; j < xSize; j++) lastLine[j] = nextLine[j];
13	poBand->RasterIO(GF_Read, 0, i, xSize, 1, nextLine, xSize, 1, GDT_UInt32, 0, 0);
14		
15	for (int j = 0; j < xSize; j++) {
16		//判断是否与左侧同图斑
17		if (j == 0 || nextLine[j] != nextLine[j - 1]) {//左不等
18			//判断是否与上侧同图斑
19			if (nextLine[j] == lastLine[j]) pDstData[i * xSize + j] = pDstData[(i - 1) * xSize + j];
20			else {
21				pDstData[i * xSize + j] = flag;
22				flag++;
23			}
24		}
25		else {//左相等
26			if (nextLine[j] == lastLine[j]) {//上相等
27				if (pDstData[(i - 1) * xSize + j] == pDstData[i * xSize + j - 1]) {//同图斑
28					pDstData[i * xSize + j] = pDstData[i * xSize + j - 1];
29				}
30				else {//不同图斑
31					pDstData[i * xSize + j] = pDstData[(i - 1) * xSize + j];
32					for (int k = 0; k <= i; k++) {
33						for (int l = 0; l < xSize; l++) {
34							if (pDstData[k * xSize + l] == pDstData[i * xSize + j - 1]) pDstData[k * xSize + l] = pDstDatz[(i - 1) * xSize + j];
35						}
36					}
37				}
38			}
39			else pDstData[i * xSize + j] = pDstData[i * xSize + j - 1];
40		}
41	}	
42}
43//处理统计图像
44int* psData = new int[xSize * ySize];
45int* flagCount = new int[flag + 1];
46for (int i = 0; i < flag + 1; i++) flagCount[i] = 0;
47for (int i = 0; i < ySize; i++) {
48	for (int j = 0; j < xSize; j++) {
49		for (int k = 0; k < flag + 1; k++) {
50    	if (pDstData[i * xSize + j] == k) flagCount[k]++;
51		}
52	}
53}
54for (int i = 0; i < ySize; i++) {
55	for (int j = 0; j < xSize; j++) {
56		int k = pDstData[i * xSize + j];
57		psData[i * xSize + j] = flagCount[k];
58	}
59}
60//输出图像
61powBand->RasterIO(GF_Write, 0, 0, xSize, ySize, pDstData, xSize, ySize, GDT_UInt32, 0, 0);
62posBand->RasterIO(GF_Write, 0, 0, xSize, ySize, psData, xSize, ySize, GDT_UInt32, 0, 0);

样例

截取自ArcMap。

输入图像

输入图像

结果图像

结果图像

统计图像

统计图像