基于自然断点法实现对多/单波段Tiff图像分类

由于在一维数据中自然断点法和K-means算法效果一致，因此笔者使用C++语言基于Alglib（用于实现K-mans）和GDAL（处理Tiff格式数据）两个外部库实现Tiff图像的自然断点分类，并按升序分类以GByte格式输出结果。

处理后效果展示：

上代码：

#include "gdal_priv.h"
#include "ogrsf_frmts.h"
#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "dataanalysis.h"

using namespace alglib;
using namespace std;

//基于Float/Double类型的多/单波段自然间断法分类
//initialTiff：原始数据存放路径  newTiff：自然断点法分类后数据存放路径 nInterval：分类数
void NaturalbreakpointTiff(std::string initialTiff, std::string newTiff, int nInterval)
{
	GDALAllRegister();
	//设置支持中文路径
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");

	const char * pszFile = initialTiff.c_str();
	GDALDataset *poDataset = (GDALDataset*)GDALOpen(pszFile, GA_ReadOnly);//使用只读方式打开图像 
	if (!poDataset)
	{
		printf("File: %s不能打开！\n", pszFile);
	}


	//仿射参数
	double padfTransform0[6];
	if (poDataset->GetGeoTransform(padfTransform0) == CE_Failure)
	{
		printf("获取仿射变换参数失败");
	}
	int iImgSizeX0 = poDataset->GetRasterXSize();
	int iImgSizeY0 = poDataset->GetRasterYSize();
	int nCount = poDataset->GetRasterCount(); 	//影像的波段数
	GDALDataType gdal_data_type = poDataset->GetRasterBand(1)->GetRasterDataType();//获取栅格类型
	double dnodata = poDataset->GetRasterBand(1)->GetNoDataValue();//获取空值对应大小
	int *panBandMap = new int[nCount];
	for (int i = 0; i < nCount; i++)
	{
		panBandMap[i] = i + 1;
	}
	
  if(gdal_data_type==GDT_Float32)
  {
	  float *pF_afScanblock1 = new float[iImgSizeX0*iImgSizeY0 * nCount];
	  poDataset->RasterIO(GF_Read, 0, 0, iImgSizeX0, iImgSizeY0, pF_afScanblock1, iImgSizeX0, iImgSizeY0, gdal_data_type, nCount, panBandMap, 0, 0, 0);
	  unsigned char *pUC_afScanblock1 = new unsigned char[iImgSizeX0*iImgSizeY0 * nCount];
	  for (int nBand = 0; nBand < nCount; nBand++)
	  {
		  int nEffectNum = 0;
		  for (int j = 0; j < iImgSizeX0*iImgSizeY0; j++)
		  {
			  if (abs((double)pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + j] - dnodata) > 0.000001)
			  {
				  nEffectNum++;
			  }
		  }
		  clusterizerstate s;
		  kmeansreport rep;
		  real_2d_array xy;
		  xy.setlength(nEffectNum, 1);
		  int nIndex = 0;
		  for (int j = 0; j < iImgSizeX0*iImgSizeY0; j++)
		  {
			  if (abs((double)pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + j] - dnodata) > 0.000001)
			  {
				  xy[nIndex][0] = (double)pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + j];
				  nIndex++;
			  }
		  }
		  clusterizercreate(s);
		  clusterizersetpoints(s, xy, 2);//2    Euclidean distance
		  clusterizersetkmeanslimits(s, 5, 0);
		  clusterizerrunkmeans(s, nInterval, rep);

		  double *pdValue = new double[nInterval];
		  double *pdValue_sort = new double[nInterval];
		  nIndex = 0;
		  for (int i = 0; i < iImgSizeX0*iImgSizeY0; i++)
		  {
			  if (abs((double)pF_afScanblock1[i] - dnodata) > 0.000001)
			  {
				  int a = rep.cidx[nIndex];
				  pdValue[rep.cidx[nIndex]] = pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i];
				  pdValue_sort[rep.cidx[nIndex]] = pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i];
				  nIndex++;
			  }

		  }
		  sort(pdValue_sort, pdValue_sort + nInterval);
		  vector<int>vNum;
		  for (int i = 0; i < nInterval; i++)
		  {
			  double _dValue = pdValue[i];
			  nIndex = 0;
			  while (abs(pdValue_sort[nIndex] - _dValue) > 0.00001)
			  {
				  nIndex++;
			  }
			  vNum.push_back(nIndex + 1);
		  }
		  nIndex = 0;
		  for (int i = 0; i < iImgSizeX0*iImgSizeY0; i++)
		  {
			  if (abs((double)pF_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] - dnodata) > 0.000001)
			  {
				  pUC_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] = vNum[rep.cidx[nIndex]];
				  nIndex++;
			  }
			  else
			  {
				  pUC_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] = 100;
			  }
		  }
		  delete[]pdValue_sort;
		  delete[]pdValue;
	  }

	  //写入新的tif
	  const char* pszDstFilename = newTiff.c_str();
	  GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	  GDALDataset* poDstDS = poDriver->Create(pszDstFilename, iImgSizeX0, iImgSizeY0, nCount, gdal_data_type, NULL);
	  poDstDS->SetProjection(poDataset->GetProjectionRef());//给它设置投影
	  poDstDS->SetGeoTransform(padfTransform0);//给设置空间转换的六参数
	  for (int i = 0; i < nCount; i++)
	  {
		  poDstDS->GetRasterBand(i + 1)->SetNoDataValue(100);//将空值设置为“无数据值”
	  }
	  //保存全波段影像
	  poDstDS->RasterIO(GF_Write, 0, 0, iImgSizeX0, iImgSizeY0, pUC_afScanblock1, iImgSizeX0, iImgSizeY0, GDT_Byte, nCount, panBandMap, 0, 0, 0);

	  //释放内存
	  GDALClose(poDstDS);
	  GDALClose(poDataset);
	  delete[]pF_afScanblock1;
	  delete[]pUC_afScanblock1;

  }
  else if (gdal_data_type == GDT_Float64)
  {
	  double *pD_afScanblock1 = new double[iImgSizeX0*iImgSizeY0 *nCount];
	  poDataset->RasterIO(GF_Read, 0, 0, iImgSizeX0, iImgSizeY0, pD_afScanblock1, iImgSizeX0, iImgSizeY0, gdal_data_type, nCount, panBandMap, 0, 0, 0);
	  unsigned char *pUC_afScanblock1 = new unsigned char[iImgSizeX0*iImgSizeY0 * nCount];
	  for (int nBand = 0; nBand < nCount; nBand++)
	  {
		  int nEffectNum = 0;
		  for (int j = 0; j < iImgSizeX0*iImgSizeY0; j++)
		  {
			  if (abs((double)pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 +j] - dnodata) > 0.000001)
			  {
				  nEffectNum++;
			  }
		  }
		  clusterizerstate s;
		  kmeansreport rep;
		  real_2d_array xy;
		  xy.setlength(nEffectNum, 1);
		  int nIndex = 0;
		  for (int j = 0; j < iImgSizeX0*iImgSizeY0; j++)
		  {
			  if (abs((double)pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + j] - dnodata) > 0.000001)
			  {
				  xy[nIndex][0] = (double)pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + j];
				  nIndex++;
			  }
		  }
		  clusterizercreate(s);
		  clusterizersetpoints(s, xy, 2);//2    Euclidean distance
		  clusterizersetkmeanslimits(s, 5, 0);
		  clusterizerrunkmeans(s, nInterval, rep);

		  double *pdValue = new double[nInterval];
		  double *pdValue_sort = new double[nInterval];
		  nIndex = 0;
		  for (int i = 0; i < iImgSizeX0*iImgSizeY0; i++)
		  {
			  if (abs((double)pD_afScanblock1[i] - dnodata) > 0.000001)
			  {
				  int a = rep.cidx[nIndex];
				  pdValue[rep.cidx[nIndex]] = pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i];
				  pdValue_sort[rep.cidx[nIndex]] = pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i];
				  nIndex++;
			  }

		  }
		  sort(pdValue_sort, pdValue_sort + nInterval);
		  vector<int>vNum;
		  for (int i = 0; i < nInterval; i++)
		  {
			  double _dValue = pdValue[i];
			  nIndex = 0;
			  while (abs(pdValue_sort[nIndex] - _dValue) > 0.00001)
			  {
				  nIndex++;
			  }
			  vNum.push_back(nIndex + 1);
		  }
		  nIndex = 0;
		  for (int i = 0; i < iImgSizeX0*iImgSizeY0; i++)
		  {
			  if (abs((double)pD_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] - dnodata) > 0.000001)
			  {
				  pUC_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] = vNum[rep.cidx[nIndex]];
				  nIndex++;
			  }
			  else
			  {
				  pUC_afScanblock1[nBand*iImgSizeX0*iImgSizeY0 + i] = 100;
			  }
		  }
		  delete[]pdValue_sort;
		  delete[]pdValue;
	  }
	 
	  //写入新的tif
	  const char* pszDstFilename = newTiff.c_str();
	  GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	  GDALDataset* poDstDS = poDriver->Create(pszDstFilename, iImgSizeX0, iImgSizeY0, nCount, gdal_data_type, NULL);
	  poDstDS->SetProjection(poDataset->GetProjectionRef());//给它设置投影
	  poDstDS->SetGeoTransform(padfTransform0);//给设置空间转换的六参数
	  for (int i = 0; i < nCount; i++)
	  {
		  poDstDS->GetRasterBand(i + 1)->SetNoDataValue(100);//将空值设置为“无数据值”
	  }
	  //保存全波段影像
	  poDstDS->RasterIO(GF_Write, 0, 0, iImgSizeX0, iImgSizeY0, pUC_afScanblock1, iImgSizeX0, iImgSizeY0, GDT_Byte, nCount, panBandMap, 0, 0, 0);

	  //释放内存
	  GDALClose(poDstDS);
	  GDALClose(poDataset);
	  delete[]pD_afScanblock1;
	  delete[]pUC_afScanblock1;
  }
 
	delete[]panBandMap;
}