c++实现影像重采样

2024-07-14 04:43| 来源: 网络整理| 查看: 265

c++实现影像重采样

概要

本片文章使用C++语言GDAL库对遥感影像数据进行重采样，分别按照重采样到输出影像的宽、高和按照分辨率的倍数两种方式进行重采样。

脚本介绍

引入其他.h和.cpp文件如下：这里PathLib.h是我个人工具类，大家可以按照这个来组织管理自己的工具方法，process.h这个是进度条工具，网上可以直接下载使用，这里我在下面把代码拷贝进来，大家也可以直接使用；最主要的是resample.h和resample.cpp两个重采样脚本以及main主函数；在这里插入图片描述

技术细节

这里话不多说，直接上代码！ 1、process.h代码如下：

#pragma once #pragma once #include #include using namespace std; #ifndef MAX #define MAX(a, b) (((a) > (b)) ? (a) : (b)) #endif // !1 #ifndef MIN #define MIN(a, b) (((a) > (b)) ? (b) : (a)) #endif // !1 class CProcessBase { public: /** *@brief构造函数 */ CProcessBase() { m_dPosition = 0.0; m_iStepCount = 100; m_iCurStep = 0; m_bIsContinue = true; } /** *@brief析构函数 */ virtual~CProcessBase() {} /** *@brief设置进度信息 *@parampszMsg 进度信息 */ virtual void SetMessage(const char* pszMsg) = 0; /** *@brief设置进度值 *@paramdPosition 进度值 *@return 返回是否取消的状态，true为不取消，false为取消 */ virtual bool SetPosition(double dPosition) = 0; /** *@brief进度条前进一步，返回true表示继续，false表示取消 *@return 返回是否取消的状态，true为不取消，false为取消 */ virtual bool StepIt() = 0; /** *@brief设置进度个数 *@paramiStepCount 进度个数 */ virtual void SetStepCount(int iStepCount) { ReSetProcess(); m_iStepCount = iStepCount; } /** *@brief获取进度信息 *@return 返回当前进度信息 */ string GetMessage() { return m_strMessage; } /** *@brief获取进度值 *@return 返回当前进度值 */ double GetPosition() { return m_dPosition; } /** *@brief重置进度条 */ void ReSetProcess() { m_dPosition = 0.0; m_iStepCount = 100; m_iCurStep = 0; m_bIsContinue = true; } /*!进度信息*/ string m_strMessage; /*!进度值*/ double m_dPosition; /*!进度个数*/ int m_iStepCount; /*!进度当前个数*/ int m_iCurStep; /*!是否取消，值为false时表示计算取消*/ bool m_bIsContinue; }; /** *@brief控制台进度条类 * *提供控制台程序的进度条类接口，来反映当前算法的进度值 */ class CConsoleProcess : public CProcessBase { public: /** *@brief构造函数 */ CConsoleProcess() { m_dPosition = 0.0; m_iStepCount = 100; m_iCurStep = 0; }; /** *@brief析构函数 */ ~CConsoleProcess() { //remove(m_pszFile); }; /** *@brief设置进度信息 *@parampszMsg 进度信息 */ void SetMessage(const char* pszMsg) { m_strMessage = pszMsg; printf("%s\n", pszMsg); } /** *@brief设置进度值 *@paramdPosition 进度值 *@return返回是否取消的状态，true为不取消，false为取消 */ bool SetPosition(double dPosition) { m_dPosition = dPosition; TermProgress(m_dPosition); m_bIsContinue = true; return true; } /** *@brief进度条前进一步 *@return返回是否取消的状态，true为不取消，false为取消 */ bool StepIt() { m_iCurStep++; m_dPosition = m_iCurStep * 1.0 / m_iStepCount; TermProgress(m_dPosition); m_bIsContinue = true; return true; } private: void TermProgress(double dfComplete) { static int nLastTick = -1; int nThisTick = (int)(dfComplete * 40.0); nThisTick = MIN(40, MAX(0, nThisTick)); if (nThisTick = 39) nLastTick = -1; if (nThisTick nLastTick) { nLastTick++; if (nLastTick % 4 == 0) fprintf(stdout, "%d", (nLastTick / 4) * 10); else fprintf(stdout, "."); } if (nThisTick == 40) fprintf(stdout, "-done.\n"); else fflush(stdout); } }; /** *\brief调用GDAL进度条接口 * *该函数用于将GDAL算法中的进度信息导出到CProcessBase基类中，供给界面显示 * *@paramdfComplete 完成进度值，其取值为0.0到1.0之间 *@parampszMessage 进度信息 *@parampProgressArgCProcessBase的指针 * *@return返回TRUE表示继续计算，否则为取消 */ inline int ALGTermProgress(double dfComplete, const char* pszMessage, void* pProgressArg) { if (pProgressArg != NULL) { CProcessBase* pProcess = (CProcessBase*)pProgressArg; pProcess->m_bIsContinue = pProcess->SetPosition(dfComplete); if (pProcess->m_bIsContinue) return true; else return false; } else return true; }

2、resample.h代码如下：

#pragma once #include "gdal_priv.h" #include "gdalwarper.h" #include "ogr_spatialref.h" #include "iostream" //#include "PathLib.h" #include "process.h" #include "gdal_priv.h" #include "malloc.h" #include "fstream" #include "string" using namespace std; bool ImageResample2(const char* pszSrcFile, const char* pszDstFile, double X, double Y, GDALResampleAlg eResampleMethod, const char* pszFormat, CProcessBase* pProcess); bool ImageResample3(const char* pszSrcFile, const char* pszDstFile, double dResX, double dResY, GDALResampleAlg eResampleMethod, const char* pszFormat, CProcessBase* pProcess);

3、resample.cpp代码如下：

#pragma once #include "resample.h" 重采样类型选择 //GDALResampleAlg eResampleMethod; //eResampleMethod = GRA_NearestNeighbour; //最邻近内插法 eResampleMethod = GRA_Bilinear; //双线性内插法 eResampleMethod = GRA_Cubic; //三次卷积法内插法 /*按照输出影像宽高重采样*/ bool ImageResample2(const char* pszSrcFile, const char* pszDstFile, double X, double Y, GDALResampleAlg eResampleMethod, const char* pszFormat, CProcessBase* pProcess) { // 前面需要对输入输出的文件路径进行判断，这里不进行判断直接处理 if (pProcess != NULL) { pProcess->ReSetProcess(); pProcess->SetMessage("开始重采样..."); } GDALAllRegister(); GDALDataset* pSrcDS = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly); GDALDataType eDT = pSrcDS->GetRasterBand(1)->GetRasterDataType(); int iBandCount = pSrcDS->GetRasterCount(); int iSrcWidth = pSrcDS->GetRasterXSize(); int iSrcHeight = pSrcDS->GetRasterYSize(); // 根据采样比例计算重采样后的图像宽高 /*int iDstWidth = iSrcWidth; int iDstHeight = iSrcHeight;*/ double dResX = 0.0; double dResY = 0.0; int iDstWidth = 0; int iDstHeight = 0; if (X == 0.0 || Y == 0.0) { dResX = 1.0; dResY = 1.0; iDstWidth = iSrcWidth; iDstHeight = iSrcHeight; } else { dResX = (X * 1.0) / (iSrcWidth * 1.0); dResY = (Y * 1.0) / (iSrcHeight * 1.0); iDstWidth = X; iDstHeight = Y; } /*double dResX = (X * 1.0) / (iSrcWidth * 1.0); double dResY = (Y * 1.0) / (iSrcHeight * 1.0);*/ /*int iDstWidth = static_cast(iSrcWidth * dResX + 0.5); int iDstHeight = static_cast(iSrcHeight * dResY + 0.5);*/ double adfGeoTransform[6] = { 0.0 }; pSrcDS->GetGeoTransform(adfGeoTransform); // 计算采样后的图像的分辨率 adfGeoTransform[1] = adfGeoTransform[1] / dResX; adfGeoTransform[5] = adfGeoTransform[5] / dResY; // 创建输出文件并设置空间参考和坐标信息 GDALDriver* poDriver = (GDALDriver*)GDALGetDriverByName(pszFormat); GDALDataset* pDstDS = poDriver->Create(pszDstFile, iDstWidth, iDstHeight, iBandCount, eDT, NULL); pDstDS->SetGeoTransform(adfGeoTransform); pDstDS->SetProjection(pSrcDS->GetProjectionRef()); // 构造坐标转换关系 void* hTransformArg = NULL; hTransformArg = GDALCreateGenImgProjTransformer2((GDALDatasetH)pSrcDS, (GDALDatasetH)pDstDS, NULL); GDALTransformerFunc pfnTransformer = GDALGenImgProjTransform; // 构造GDALWarp的变换选项 GDALWarpOptions* psWO = GDALCreateWarpOptions(); psWO->papszWarpOptions = CSLDuplicate(NULL); psWO->eWorkingDataType = eDT; psWO->eResampleAlg = eResampleMethod; psWO->hSrcDS = (GDALDatasetH)pSrcDS; psWO->hDstDS = (GDALDatasetH)pDstDS; psWO->pfnTransformer = pfnTransformer; psWO->pTransformerArg = hTransformArg; psWO->pfnProgress = ALGTermProgress; psWO->pProgressArg = pProcess; psWO->nBandCount = iBandCount; psWO->panSrcBands = (int*)CPLMalloc(psWO->nBandCount * sizeof(int)); psWO->panDstBands = (int*)CPLMalloc(psWO->nBandCount * sizeof(int)); for (int i = 0; i panSrcBands[i] = i + 1; psWO->panDstBands[i] = i + 1; } // 创建GDALWarp执行对象，并使用GDALWarpOptions来进行初始化 GDALWarpOperation oWO; oWO.Initialize(psWO); // 执行处理 oWO.ChunkAndWarpImage(0, 0, iDstWidth, iDstHeight); // 释放资源和关闭文件 GDALDestroyGenImgProjTransformer(psWO->pTransformerArg); GDALDestroyWarpOptions(psWO); //GDALDriver* pDriverPNG = GetGDALDriverManager()->GetDriverByName("PNG"); //pDriverPNG->CreateCopy("C:\\Users\\17733634114\\Desktop\\out\\data\\GF1_PMS1_E117.0_N39.4_20210128_L1A0005436461-PAN1-rc-ac-rpc_resample.png", pDstDS, TRUE, 0, 0, 0); //创建jpg文件 //cout ReSetProcess(); pProcess->SetMessage("开始重采样..."); } GDALAllRegister(); GDALDataset* pSrcDS = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly); GDALDataType eDT = pSrcDS->GetRasterBand(1)->GetRasterDataType(); int iBandCount = pSrcDS->GetRasterCount(); int iSrcWidth = pSrcDS->GetRasterXSize(); int iSrcHeight = pSrcDS->GetRasterYSize(); // 根据采样比例计算重采样后的图像宽高 /*int iDstWidth = iSrcWidth; int iDstHeight = iSrcHeight;*/ int iDstWidth = static_cast(iSrcWidth * dResX + 0.5); int iDstHeight = static_cast(iSrcHeight * dResY + 0.5); double adfGeoTransform[6] = { 0.0 }; pSrcDS->GetGeoTransform(adfGeoTransform); // 计算采样后的图像的分辨率 adfGeoTransform[1] = adfGeoTransform[1] / dResX; adfGeoTransform[5] = adfGeoTransform[5] / dResY; // 创建输出文件并设置空间参考和坐标信息 GDALDriver* poDriver = (GDALDriver*)GDALGetDriverByName(pszFormat); GDALDataset* pDstDS = poDriver->Create(pszDstFile, iDstWidth, iDstHeight, iBandCount, eDT, NULL); pDstDS->SetGeoTransform(adfGeoTransform); pDstDS->SetProjection(pSrcDS->GetProjectionRef()); // 构造坐标转换关系 void* hTransformArg = NULL; hTransformArg = GDALCreateGenImgProjTransformer2((GDALDatasetH)pSrcDS, (GDALDatasetH)pDstDS, NULL); GDALTransformerFunc pfnTransformer = GDALGenImgProjTransform; // 构造GDALWarp的变换选项 GDALWarpOptions* psWO = GDALCreateWarpOptions(); psWO->papszWarpOptions = CSLDuplicate(NULL); psWO->eWorkingDataType = eDT; psWO->eResampleAlg = eResampleMethod; psWO->hSrcDS = (GDALDatasetH)pSrcDS; psWO->hDstDS = (GDALDatasetH)pDstDS; psWO->pfnTransformer = pfnTransformer; psWO->pTransformerArg = hTransformArg; psWO->pfnProgress = ALGTermProgress; psWO->pProgressArg = pProcess; psWO->nBandCount = iBandCount; psWO->panSrcBands = (int*)CPLMalloc(psWO->nBandCount * sizeof(int)); psWO->panDstBands = (int*)CPLMalloc(psWO->nBandCount * sizeof(int)); for (int i = 0; i panSrcBands[i] = i + 1; psWO->panDstBands[i] = i + 1; } // 创建GDALWarp执行对象，并使用GDALWarpOptions来进行初始化 GDALWarpOperation oWO; oWO.Initialize(psWO); // 执行处理 oWO.ChunkAndWarpImage(0, 0, iDstWidth, iDstHeight); // 释放资源和关闭文件 GDALDestroyGenImgProjTransformer(psWO->pTransformerArg); GDALDestroyWarpOptions(psWO); GDALClose((GDALDatasetH)pSrcDS); GDALClose((GDALDatasetH)pDstDS); if (pProcess != NULL) pProcess->SetMessage("重采样完成！"); return true; }

4、main.cpp代码如下：

#include "resample.h" 重采样类型选择 //GDALResampleAlg eResampleMethod; //eResampleMethod = GRA_NearestNeighbour; //最邻近内插法 eResampleMethod = GRA_Bilinear; //双线性内插法 eResampleMethod = GRA_Cubic; //三次卷积法内插法 using namespace std; int main(int argc, const char* argv[]) { #ifdef _DEBUG printf("debug model\n"); const char* input = "C:\\Users\\Desktop\\PNG\\ground.tif"; const char* output = "C:\\Users\\Desktop\\PNG\\resample.tif"; double xx = stold("0.5"); double yy = stold("0.5"); #else if (argc != 5) { cout

【本文地址】

c++实现影像重采样

c++实现影像重采样

今日新闻

推荐新闻