栅格图斑聚类分析及C++实现
算法作用
在进行非监督分类之后，往往需要进一步处理才能获得理想的分类结果。栅格图斑聚类分析就是其中一种方式。
通过栅格图斑聚类，我们最终会输出两个图像：一是像元点Data值为像元点所对应的图斑序号；二是像元点Data值为图斑面积值的图像。
算法原理
图斑的聚类需要判断像元A与左侧、上方像元的关系：
| 与左侧关系 |
与上方关系 |
处理方式 |
| 不同 |
不同 |
产生一个新的图斑 |
| 相同 |
不同 |
合并入左侧图斑 |
| 不同 |
相同 |
合并入上方图斑 |
| 相同 |
相同 |
判断两图斑是否属于同一图斑序号
若相同则合并进入该图斑
若不同则需要合并两个图斑 |
分析过程
读取数据
老一套，首先在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");
|
创建数据集，打开经过K-Means聚类处理后的图像，并获得条带数据。
1
2
3
4
5
| GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("Gtiff");
GDALDataset* poData = (GDALDataset*)GDALOpen(inPath, GA_ReadOnly);
int xSize = poData->GetRasterXSize();
int ySize = poData->GetRasterYSize();
GDALRasterBand* poBand = poData->GetRasterBand(1);
|
同理创建用于输出的数据集，以及对应的条带数据备用。
1
2
3
4
| GDALDataset* powData = poDriver->Create(outPath, xSize, ySize, 1, GDT_Float64, NULL);
GDALDataset* posData = poDriver->Create(statPath, xSize, ySize, 1, GDT_Float64, NULL);
GDALRasterBand* powBand = powData->GetRasterBand(1);
GDALRasterBand* posBand = posData->GetRasterBand(1);
|
数据处理
通过按行读取数据的方式，获得上下两行数据：lastLine和nextLine，便于对比像元与上方像元的数值。
图像处理时需要注意第一列以及第一行的处理方式。这里在处理第一行时将lastLine初始化为-1进行相等判断，即第一行与“上方像元”不同；处理第一列时进行特判，即第一列与“左侧像元”不同。
合并图斑时注意左侧图斑与上方图斑序号不同时的情况，此时需要将左侧图斑全部合并入上方像元所在的图斑，因此需要遍历左侧像元上方全部的行列。
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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
| int* pDstData = new int[xSize * ySize];
for (int i = 0; i < xSize * ySize; i++) pDstData[i] = -1;
int flag = 0;
int* lastLine = new int[xSize];
int* nextLine = new int[xSize];
for (int i = 0; i < xSize; i++) {
lastLine[i] = -1;
nextLine[i] = -1;
}
for (int i = 0; i < ySize; i++) {
for (int j = 0; j < xSize; j++) lastLine[j] = nextLine[j];
poBand->RasterIO(GF_Read, 0, i, xSize, 1, nextLine, xSize, 1, GDT_UInt32, 0, 0);
for (int j = 0; j < xSize; j++) {
if (j == 0 || nextLine[j] != nextLine[j - 1]) {
if (nextLine[j] == lastLine[j]) pDstData[i * xSize + j] = pDstData[(i - 1) * xSize + j];
else {
pDstData[i * xSize + j] = flag;
flag++;
}
}
else {
if (nextLine[j] == lastLine[j]) {
if (pDstData[(i - 1) * xSize + j] == pDstData[i * xSize + j - 1]) {
pDstData[i * xSize + j] = pDstData[i * xSize + j - 1];
}
else {
pDstData[i * xSize + j] = pDstData[(i - 1) * xSize + j];
for (int k = 0; k <= i; k++) {
for (int l = 0; l < xSize; l++) {
if (pDstData[k * xSize + l] == pDstData[i * xSize + j - 1]) pDstData[k * xSize + l] = pDstDatz[(i - 1) * xSize + j];
}
}
}
}
else pDstData[i * xSize + j] = pDstData[i * xSize + j - 1];
}
}
}
int* psData = new int[xSize * ySize];
int* flagCount = new int[flag + 1];
for (int i = 0; i < flag + 1; i++) flagCount[i] = 0;
for (int i = 0; i < ySize; i++) {
for (int j = 0; j < xSize; j++) {
for (int k = 0; k < flag + 1; k++) {
if (pDstData[i * xSize + j] == k) flagCount[k]++;
}
}
}
for (int i = 0; i < ySize; i++) {
for (int j = 0; j < xSize; j++) {
int k = pDstData[i * xSize + j];
psData[i * xSize + j] = flagCount[k];
}
}
powBand->RasterIO(GF_Write, 0, 0, xSize, ySize, pDstData, xSize, ySize, GDT_UInt32, 0, 0);
posBand->RasterIO(GF_Write, 0, 0, xSize, ySize, psData, xSize, ySize, GDT_UInt32, 0, 0);
|
样例
截取自ArcMap。
输入图像
结果图像
统计图像
