任务

用C++和GDAL实现K均值分类

原理

关于K均值算法的帖子很多，原理也不复杂：KMeans 算法

代码

程序结构：

在这里插入图片描述

主函数 main.cpp

/* 
 *实验8 K均值分类
 *时间 20210715
 *这里使用的影像是单波段的 uint8
 */

#include <iostream>
#include <gdal.h>
#include "gdal_priv.h"
#include "func.h"

using namespace std;

int main()
{
	// 注册
	GDALAllRegister();
	// 路径
	char ImgPath[] = "D:\\Practice\\org\\exp8_Img.tif";
	char SavePath[] = "D:\\Practice\\res\\exp8_classification.tif";
	// 打开数据
	GDALDataset* ImgSet = (GDALDataset*)GDALOpen(ImgPath, GA_ReadOnly);
	// 获取大小
	int mX = ImgSet->GetRasterXSize();
	int mY = ImgSet->GetRasterYSize();
	GDALDataType mDataType = ImgSet->GetRasterBand(1)->GetRasterDataType();

	// 调用K均值分类的函数进行分类
	unsigned char* Reslut = new unsigned char[mX*mY];		//保存分类结果
	int Num = 3;				//类别数
	mKmeans(ImgSet, Reslut, Num);

	// 创建输出的数据集
	GDALDriver* hDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset* mSaveSet = hDriver->Create(SavePath, mX, mY, 1, GDT_Byte, NULL);

	// 获取影像的地理信息和投影信息
	double geoInformation[6];
	ImgSet->GetGeoTransform(geoInformation);
	const char* gdalProjection = ImgSet->GetProjectionRef();

	// 将影像的地理信息写到纠正后影像数据集中
	mSaveSet->SetGeoTransform(geoInformation);
	mSaveSet->SetProjection(gdalProjection);

	// 写入数据集
	mSaveSet->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, mX, mY, Reslut, mX, mY, GDT_Byte, 0, 0);


	// 释放内存,关闭数据集
	delete[] Reslut;
	GDALClose(ImgSet);
	GDALClose(mSaveSet);
	GDALDestroyDriverManager();

	cout << "程序运行完毕！！！" << endl;
	getchar();
	return 0;
}

头文件 func.h

#pragma once
#include "gdal.h"
#include "gdal_priv.h"
#include <iostream>

using namespace std;

// 进行K均值分类计算的函数,三个输入分别为待分类影像，分类结果和类别个数
bool mKmeans(GDALDataset* mImgSet, unsigned char* Reslut, int Num);

函数文件func.cpp

#include "func.h"


bool mKmeans(GDALDataset* mImgSet, unsigned char* Reslut, int Num)
{
	// 获取大小
	int mX = mImgSet->GetRasterXSize();
	int mY = mImgSet->GetRasterYSize();
	GDALDataType mDataType = mImgSet->GetRasterBand(1)->GetRasterDataType();
	// 取波段
	GDALRasterBand* mBand = mImgSet->GetRasterBand(1);
	// 取数据
	unsigned char* mBand_ = new unsigned char[mX*mY*mDataType];
	mBand->RasterIO(GF_Read, 0, 0, mX, mY, mBand_, mX, mY, mDataType, 0, 0);

	// 用一个数组来表示类别，类别有0 1 2 3 …… Num-1
	unsigned char* label = new unsigned char[mX*mY];

	// 保存聚类中心,并赋予初始值
	double* Center = new double[Num];
	for (int i = 0; i < Num; i++)
	{
		Center[i] =  i * 255 / (Num - 1);
	}
	// 用来临时保存新的聚类中心
	double* NewCenter = new double[Num];
	memset(NewCenter, 0, sizeof(double));

	// 循环迭代，指导分类结束
	int iter = 0;		//记录迭代次数
	int max_iter = 50;	//最大迭代次数，超过这个次数还没有分类成功说明不收敛
	while (true)
	{
		cout << "Center[0]= " << Center[0] << "  Center[1]= " << Center[1] << "  Center[2]= " << Center[2] << endl;
		// 记录每种类别的像素数量
		double* cnt = new double[Num];
		memset(cnt, 0, sizeof(double));

		// 记录每种类别像素值的和
		double* sum = new double[Num];
		memset(sum, 0, sizeof(double));

		// 到每个聚类中心的距离
		double* d = new double[Num];
		memset(d, 0, sizeof(double));

		// 循环遍历所有像素
		for (int i = 0; i < mX*mY; i++)
		{
			// 计算到每个聚类中心的距离
			double minDist = 255;	//用来保存每个像素的最小类中心距离
			int tag = -1;			// 最小距离类中心的类别

			for (int j = 0; j < Num; j++)
			{
				d[j] = abs(double(mBand_[i]) - Center[j]);
				if (d[j] <= minDist) { minDist = d[j], tag = j; }
			}
			
			// 计算每种类别的像素值的和
			sum[tag] += double(mBand_[i]);
			cnt[tag] ++;
			Reslut[i] = 255.0 * tag / (Num - 1);
	
		}
	
		// 计算新的聚类中心
		for (int j = 0; j < Num; j++)
		{
			if (cnt[j] == 0) { cout << "错误！" << endl; }
			else { NewCenter[j] = sum[j] / cnt[j]; }
		}

		// 判断是否跳出循环
		bool flag = 1;
		// 若新的聚类中心与之前的聚类中心相同，跳出循环，分类结束

		for (int j = 0; j < Num; j++)
		{
			
			if (NewCenter[j] != Center[j]) { flag = 0; }

			Center[j] = NewCenter[j];
		}
		cout << endl;
		if (flag == 1) { cout << "分类成功，迭代次数为 " << iter << endl; break; }

		// 达到最大循环次数，也跳出循环
		if (iter == max_iter) { cout << "达到最大循环次数，看来是不收敛了！！！" << endl; break;}


		iter++;

	}

	return true;
}