GDAL图像分割算法实现

一、图像分割在arcgis中叫做分割栅格，在工具箱里面打开，可以看到基本界面如下；arcgis提供了两种分割方法：按数量和按尺寸分割。

二、按数量分割是输入列和行，输出列*行个被分割的图像，同时还要保存原始图像的地理信息和投影信息。原始图像宽除以列得到输出图像的宽，原始图像高除以行得到输出图像的高，如果除不尽则向上取整，所以在边缘部分会有黑边。arcgis分幅方法是：从下往上、从左往右裁剪。本文从上往下、从左往右，能够正确处理黑边。

三、编译环境：VS2012+GDAL，需要自己配置gdal和修改文件路径，可以参考https://blog.youkuaiyun.com/HB_Programmer/article/details/81808963，分割算法是我自己所想，里面重点是分割算法的具体实现，针对不同格式的数据还需自己加入判断。

四、源代码

#include <iostream>
#include "gdal_priv.h"//包含头文件

using namespace std;
//声明函数
void CutByNumber(const char *pszSrcFile,//输入文件
                 const char *pszDstFile,//输出文件夹
                 int iStartX,int iStartY,//起始行列号row,col
                 int iSizeX,int iSizeY,//裁剪列数和行数imgWidth,imgHeight
                 const char *pszFormat,//输出文件格式
                 int x,int y);
//声明全局变量
int a=0;

int main()
{
	GDALAllRegister();//注册驱动
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8","NO");//支持中文路径

	GDALDataset  *poDataset;
	//请输入你的正确的路径
	string pszSrcFile = "D:/Desktop/1.png";
	poDataset = (GDALDataset *) GDALOpen( pszSrcFile.data(), GA_ReadOnly );
    if( poDataset == NULL )
    {
        cout<<"poDataset is NULL"<<endl;
        return 0;
    }
	//获取图像宽、高、波段数
	int width = poDataset->GetRasterXSize();
	int height = poDataset->GetRasterYSize();

	//起始行列号
    int row = 0;
    int col = 0;
	//裁剪个数
	int x,y;
	cout<<"输入裁剪列和行(以空格隔开)：";
	cin>>x>>y;
    //裁剪大小
    int imgWidth = (ceil)(width/(float)x);
    int imgHeight = (ceil)(height/(float)y);
	
	//输出图像格式
    string type = ".png";
    //执行裁剪
    GDALRasterBand *poBand;
    string pszDstFile = "D:/Desktop/0";// D:/Desktop是输出文件夹 /0是输出图像的前缀
    string name2;
    int b = -1;
    for(int i=0;i<x*y;i++)
    {
        name2 = QString::number(i);
        pszDstFile = pszDstFile + name2.toLatin1() + type.toLatin1();

        row = imgWidth*(i%x);

        if(i%x == 0)
            b++;
        col = imgHeight*b;
        CutByNumber(pszSrcFile.data(),pszDstFile.data(),row,col,imgWidth,imgHeight,"GTiff",x,y);

        //输出统计信息
        GDALDataset *pDstDS = (GDALDataset*)GDALOpen(pszDstFile.data(),GA_ReadOnly);
        //生成均值和标准差
        double min,max,mean,stdDev;
        for(int i=1;i<=pDstDS->GetRasterCount();i++)
        {
            poBand = pDstDS->GetRasterBand(i);
            poBand->ComputeStatistics(FALSE,&min,&max,&mean,&stdDev,NULL,NULL);
        }
        //统计直方图
        unsigned long long panHistogram[256] = {0};
        for(int i=1;i<=pDstDS->GetRasterCount();i++)
        {
            poBand = pDstDS->GetRasterBand(i);
            poBand->GetHistogram(-0.5,255.5,256,panHistogram,TRUE,FALSE,NULL,NULL);
        }
        GDALClose((GDALDatasetH)pDstDS);
        pszDstFile = "D:/Desktop/0";
    }
}

//实现裁剪函数
void CutByNumber(const char *pszSrcFile,//输入文件
                 const char *pszDstFile,//输出文件夹
                 int iStartX,int iStartY,//起始行列号row,col
                 int iSizeX,int iSizeY,//裁剪列数和行数imgWidth,imgHeight
                 const char *pszFormat,//输出文件格式
                 int x,int y)
{
	GDALAllRegister();
    CPLSetConfigOption( "GDAL_FILENAME_IS_UTF8", "NO" );

    GDALDataset *pSrcDS = (GDALDataset*)GDALOpen(pszSrcFile,GA_ReadOnly);

    if(pSrcDS == NULL)
    {
        cout<<"Open Image failed! "<<endl;
        return ;
    }

    a++;

    //获取数据类型和波段信息
    GDALDataType dataType = pSrcDS->GetRasterBand(1)->GetRasterDataType();
    int iBandCount = pSrcDS->GetRasterCount();
    int width = pSrcDS->GetRasterXSize();
    int height = pSrcDS->GetRasterYSize();
    //裁剪后的图像宽高
    int iDstWidth = iSizeX;
    int iDstHeight = iSizeY;

    //计算裁剪后图像左上角坐标
    double adfGeoTransform[6] = {0};
    pSrcDS->GetGeoTransform(adfGeoTransform);
    adfGeoTransform[0] = adfGeoTransform[0] + iStartX*adfGeoTransform[1] + iStartY*adfGeoTransform[2];
    adfGeoTransform[3] = adfGeoTransform[3] + iStartX*adfGeoTransform[4] + iStartY*adfGeoTransform[5];

    //创建输出文件并设置空间参考和坐标信息
    GDALDriver *poDriver = (GDALDriver*)GDALGetDriverByName(pszFormat);
    GDALDataset *pDstDS = poDriver->Create(pszDstFile,iDstWidth,iDstHeight,iBandCount,dataType,NULL);
    pDstDS->SetGeoTransform(adfGeoTransform);
    pDstDS->SetProjection(pSrcDS->GetProjectionRef());

    int *pBandMap = new int[iBandCount];
    for(int i = 0;i < iBandCount;i++)
        pBandMap[i] = i + 1;

    //根据类型判断，申请不同类型的缓存
    if(dataType == GDT_Byte)//如果是8bit图像
    {
        unsigned char *pDataBuff = new unsigned char[iDstWidth*iDstHeight*iBandCount];

        //防止越界
        int val1 = 0;
        int val2 = 0;

        if(a%x == 0 && iDstWidth/(float)x != 0)
        {
            val1 = x*iDstWidth - width;
            iDstWidth = iDstWidth - val1;
        }

        if(a > x*(y-1) && iSizeY/(float)y != 0)
        {
            val2 = y*iDstHeight - height;
            iDstHeight -= val2;
        }
        //上部的黑边只能加在下面，改坐标系或者从左下角读取数据
        pSrcDS->RasterIO(GF_Read,iStartX,iStartY,iDstWidth,iDstHeight,pDataBuff,
                         iDstWidth,iDstHeight,dataType,iBandCount,pBandMap,0,0,0);

        pDstDS->RasterIO(GF_Write,0,0,iDstWidth,iDstHeight,pDataBuff,
                         iDstWidth,iDstHeight,dataType,iBandCount,pBandMap,0,0,0);

    }
    else
    {
        //其他类型的图像，与8bit类型，就是申请的缓存类型不同
    }

    delete[] pBandMap;
    GDALClose((GDALDatasetH)pSrcDS);
    GDALClose((GDALDatasetH)pDstDS);
}