栅格图斑聚类分析及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。
输入图像
结果图像
统计图像