算法作用
在进行非监督分类之后，往往需要进一步处理才能获得理想的分类结果。栅格图斑聚类分析就是其中一种方式。
通过栅格图斑聚类，我们最终会输出两个图像：一是像元点Data值为像元点所对应的图斑序号；二是像元点Data值为图斑面积值的图像。
算法原理
图斑的聚类需要判断像元A与左侧、上方像元的关系：
| 与左侧关系 | 与上方关系 | 处理方式 |
|---|---|---|
| 不同 | 不同 | 产生一个新的图斑 |
| 相同 | 不同 | 合并入左侧图斑 |
| 不同 | 相同 | 合并入上方图斑 |
| 相同 | 相同 | 判断两图斑是否属于同一图斑序号 若相同则合并进入该图斑 若不同则需要合并两个图斑 |
分析过程
读取数据
老一套，首先在MFC工程中添加GDAL库，并在程序中引用并注册。
cpp
1# "gdal/include/gdal_priv.h"
2# "gdal/include/gdalwarper.h"
3# "gdal/include/ogrsf_frmts.h"
4# 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。
输入图像
结果图像
统计图像
