GDAL图像分割算法实现

阳光穿透心脏的1/2处 2022-05-14 09:40 316阅读 0赞

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

![70][]

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

三、编译环境：VS2012+GDAL，需要自己配置gdal和修改文件路径，可以参考[https://blog.csdn.net/HB\_Programmer/article/details/81808963][https_blog.csdn.net_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);
    }

[70]: /images/20220514/b5380fab4b054553be81f6196acb5473.png
[https_blog.csdn.net_HB_Programmer_article_details_81808963]: https://blog.csdn.net/HB_Programmer/article/details/81808963