基于GDAL的多波段图像压缩及二值分割

内容概述

将多波段各个像素灰度值加和，并且在保留仿射参数和投影信息情况下生成单波段图像，用于经典算法（otsu，em等）本文将以otsu为例，生成简单二值图像
原基于pie二次开发，摘取部分函数算法成文，记录学习过程，有需求的话后续会分享pie二次开发学习过程

环境

语言：C++【100%】
vs版本：【2015】
GDAL版本：【2.3.X】
OPENCV版本：【3.x】

多波段图像压缩

将多波段各个像素灰度值加和，并且在保留仿射参数和投影信息情况下生成单波段图像，用于经典算法（otsu，em等）本文将以otsu为例

代码部分

注释标注的还算全面

头文件

#pragma once
#define BYTE unsigned char

class Arr
{
public:
	int imgX = 0;
	int imgY = 0;
	GDALDataType datatype;
	
	double* EndArr = new double[imgX * imgY];
	BYTE *EndImg = (BYTE *)CPLMalloc(sizeof(double)*imgX*imgY);

};


void GdalAlgo(const char* filepath1, const char* filepath2, Arr &arr);

源文件

#include "gdal_priv.h"
#include"Algo.h"
#include<iostream>  

using namespace std;


void toGeoCoord(int x, int y, double* coords, double* transform)
{
	coords[0] = transform[0] + x * transform[1] + y * transform[2];
	coords[1] = transform[3] + x * transform[4] + y * transform[5];
}

void GdalAlgo(const char* filepath1, const char* filepath2, Arr &arr)
{
	GDALAllRegister(); //注册驱动

	GDALDataset* PD1;//建立数据集读第一张图片
	PD1 = (GDALDataset*)GDALOpen(filepath1, GA_ReadOnly);//指针指向文件路径
	if (PD1 == NULL)
	{
		//openfail
		return;
	}
	int iXSize1 = PD1->GetRasterXSize();//获取图像的宽
	int iYSize1 = PD1->GetRasterYSize();//获取图像的高
	int iBandCount1 = PD1->GetRasterCount();//获取图像的波段数
	GDALDataType datatype1 = PD1->GetRasterBand(1)->GetRasterDataType();//获取影像数据类型
	//######################################################
	double adfGeoTransform[6];//获取仿射变换六参数信息
	PD1->GetGeoTransform(adfGeoTransform);
	const char* pszProj = PD1->GetProjectionRef(); // 获得WKT形式的投影信息
	//通过仿射变换将影像左上角和右下角行列号坐标转化为地理坐标
	double adfULCoord[2];
	double adfLRCoord[2];
	toGeoCoord(0, 0, adfULCoord, adfGeoTransform);
	toGeoCoord(iYSize1 - 1, iYSize1 - 1, adfLRCoord, adfGeoTransform);
	//########################################################

	GDALDataset* PD2;//建立数据集读第二张图片 
	PD2 = (GDALDataset*)GDALOpen(filepath2, GA_ReadOnly);//指针指向文件路径
	if (PD2 == NULL)
	{
		//openfail
		return;
	}
	int iXSize2 = PD2->GetRasterXSize();//获取图像的宽
	int iYSize2 = PD2->GetRasterYSize();//获取图像的高
	int iBandCount2 = PD2->GetRasterCount();//获取图像的波段数
	GDALDataType datatype2 = PD2->GetRasterBand(1)->GetRasterDataType();//获取影像数据类型


	if (iXSize1 == iXSize2 && iYSize1 == iYSize2 && iBandCount1 == iBandCount2 && datatype1 == datatype2)//鲁棒
	{
		int imgX = iXSize1; int imgY = iYSize1; GDALDataType datatype = datatype1;
		arr.imgX = imgX; arr.imgY = imgY; arr.datatype = datatype;//对象获得长宽及数据类型
		
		int size_m1 = sizeof(double)*imgX*imgY;
		BYTE *pData1 = (BYTE *)CPLMalloc(size_m1); // 分配内存1

		int size_m2 = sizeof(double)*imgX*imgY;
		BYTE *pData2 = (BYTE *)CPLMalloc(size_m2); // 分配内存2

		int size_m3 = sizeof(double)*imgX*imgY;
		BYTE *pData3 = (BYTE *)CPLMalloc(size_m3); // 分配内存3

		double *tempp1 = new double[imgX * imgY];// 分配临时内存
		double *tempp2 = new double[imgX * imgY];// 分配临时内存
		double *tempp3 = new double[imgX * imgY];// 分配临时内存
		for (int i = 0; i < imgX; i++)//初始化（也可置于下方循环中）
		{
			for (int j = 0; j < imgY; j++)
			{
				tempp1[i * imgX + j] = 0;
				tempp2[i * imgX + j] = 0;
				tempp3[i * imgX + j] = 0;
			}
		}

		for (int ibands = 1; ibands <= iBandCount1; ibands++)//将多波段图像压缩为单波段图像（一维数组），为传入otsu算法做准备
		{
			GDALRasterBand *poBand1 = PD1->GetRasterBand(ibands);
			GDALRasterBand *poBand2 = PD2->GetRasterBand(ibands);

			poBand1->RasterIO(GF_Read, 0, 0, imgX, imgY, pData1, imgX, imgY, GDT_Byte, 0, 0);
			poBand2->RasterIO(GF_Read, 0, 0, imgX, imgY, pData2, imgX, imgY, GDT_Byte, 0, 0);

			for (int i = 0; i < imgY; i++)//循环获得整图所有波段累计灰度值
			{
				for (int j = 0; j < imgX; j++)
				{
					tempp1[i * imgY + j] += pData1[i * imgY + j];
					tempp2[i * imgY + j] += pData2[i * imgY + j];

					//pData3[i * imgX + j] += pData4[i * imgX + j];
					//pData4[i * imgX + j] = pData1[i * imgX + j] - pData2[i * imgX + j];
				}
			}
		}
		for (int i = 0; i < imgY; i++)//做差获得2时相的一维灰度差
		{
			for (int j = 0; j < imgX; j++)
			{
			
				pData3[i * imgY + j] = fabs(tempp1[i * imgY + j] - tempp2[i * imgY + j]);
				tempp3[i * imgY + j] = fabs(tempp1[i * imgY + j] - tempp2[i * imgY + j]);
		
				//此段数据转换留于自行试写
				//BYTE tempa;
				//tempa = static_cast<BYTE>(a3);
			}
		}
		arr.EndArr = tempp3;//对象获得最终 double 类型一维数组
		arr.EndImg = pData3;//对象获得最终 无数据 类型一维数组
	}


	GDALDriver* tempDriver;//写图像驱动
	GDALDataset* tempData;
	const char* pszFormat = "GTiff";
	const char* pszDstFile = "NewPic.tif";
	tempDriver = GetGDALDriverManager()->GetDriverByName(pszFormat);
	tempData = tempDriver->Create(pszDstFile, arr.imgX, arr.imgY, 1, GDT_Byte, NULL);//创建文件用来写入计算后的数据

	GDALRasterBand* temptif = tempData->GetRasterBand(1);
	tempData->SetGeoTransform(adfGeoTransform);
	tempData->SetProjection(pszProj);
	temptif->RasterIO(GF_Write, 0, 0, arr.imgX, arr.imgY, arr.EndImg, arr.imgX, arr.imgY, GDT_Byte, 0, 0);//写temp图像完成
	
    //release space
    GDALClose(tempData);
	delete PD1;
	delete PD2;

	
}

原代码有调用qt部分，已做更改，打包即用。
突然懒得写了，后续代码心情好了再补上XD