#include "stdafx.h"
#include "ImageOperatorBase.h"
#include "BaseTool.h"
#include "GeoOperator.h"
#include <Eigen/Core>
#include <Eigen/Dense>
#include <omp.h>
#include <io.h>
#include <stdio.h>
#include <stdlib.h>
#include <gdal.h>
#include <gdal_utils.h>
#include <gdal_priv.h>
#include <gdalwarper.h>
#include <proj.h>
#include <string.h>
#include <memory>
#include <iostream>
#include "FileOperator.h"
#include <opencv2/opencv.hpp>
#include <QMessageBox>
#include <QDir>
#include <fcntl.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <QProgressDialog>
#include <gdal_priv.h>
#include <ogr_spatialref.h>  // OGRSpatialReference 用于空间参考转换
#include <gdal_alg.h>        // 用于 GDALWarp 操作



/**
 *  输入数据是ENVI格式数据
 */

std::shared_ptr<GDALDataset> OpenDataset(const QString& in_path, GDALAccess rwmode)
{
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	GDALDataset* dataset_ptr = (GDALDataset*)(GDALOpen(in_path.toUtf8().constData(), rwmode));
	std::shared_ptr<GDALDataset> rasterDataset(dataset_ptr, CloseDataset);
	return rasterDataset;
}

void CloseDataset(GDALDataset* ptr)
{
	GDALClose(ptr);
	ptr = NULL;
}

int TIFF2ENVI(QString in_tiff_path, QString out_envi_path)
{
	std::shared_ptr<GDALDataset> ds		   = OpenDataset(in_tiff_path);
	const char*					 args[]	   = { "-of", "ENVI", NULL };
	GDALTranslateOptions*		 psOptions = GDALTranslateOptionsNew((char**)args, NULL);
	GDALClose(GDALTranslate(out_envi_path.toUtf8().constData(), ds.get(), psOptions, NULL));
	GDALTranslateOptionsFree(psOptions);
	return 0;
}

int ENVI2TIFF(QString in_envi_path, QString out_tiff_path)
{
	std::shared_ptr<GDALDataset> ds		   = OpenDataset(in_envi_path);
	const char*					 args[]	   = { "-of", "Gtiff", NULL };
	GDALTranslateOptions*		 psOptions = GDALTranslateOptionsNew((char**)args, NULL);
	GDALClose(GDALTranslate(out_tiff_path.toUtf8().constData(), ds.get(), psOptions, NULL));
	GDALTranslateOptionsFree(psOptions);
	return 0;
}

int CreateDataset(QString new_file_path, int height, int width, int band_num, double* gt,
				  QString projection, GDALDataType gdal_dtype, bool need_gt)
{
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	GDALDriver*					 poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
	std::shared_ptr<GDALDataset> poDstDS(poDriver->Create(new_file_path.toUtf8().constData(), width,
														  height, band_num, gdal_dtype, NULL));
	if(need_gt) {
		poDstDS->SetProjection(projection.toUtf8().constData());
		poDstDS->SetGeoTransform(gt);
	} else {
	}
	GDALFlushCache((GDALDatasetH)poDstDS.get());
	return 0;
}

int saveDataset(QString new_file_path, int start_line, int start_cols, int band_ids, int datacols,
				int datarows, void* databuffer)
{
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	std::shared_ptr<GDALDataset> poDstDS	   = OpenDataset(new_file_path, GA_Update);
	GDALDataType				 gdal_datatype = poDstDS->GetRasterBand(1)->GetRasterDataType();
	poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_cols, start_line, datacols, datarows,
											   databuffer, datacols, datarows, gdal_datatype, 0, 0);
	GDALFlushCache(poDstDS.get());
	return 0;
}

int block_num_pre_memory(int block_width, int height, GDALDataType gdal_datatype, double memey_size)
{
	// 计算大小
	int size_meta = 0;

	if(gdal_datatype == GDT_Byte) {
		size_meta = 1;
	} else if(gdal_datatype == GDT_UInt16) {
		size_meta = 2; // 只有双通道才能构建 复数矩阵
	} else if(gdal_datatype == GDT_UInt16) {
		size_meta = 2;
	} else if(gdal_datatype == GDT_Int16) {
		size_meta = 2;
	} else if(gdal_datatype == GDT_UInt32) {
		size_meta = 4;
	} else if(gdal_datatype == GDT_Int32) {
		size_meta = 4;
	}
	// else if (gdal_datatype == GDT_UInt64) {
	//	size_meta = 8;
	// }
	// else if (gdal_datatype == GDT_Int64) {
	//	size_meta = 8;
	// }
	else if(gdal_datatype == GDT_Float32) {
		size_meta = 4;
	} else if(gdal_datatype == GDT_Float64) {
		size_meta = 4;
	} else if(gdal_datatype == GDT_CInt16) {
		size_meta = 2;
	} else if(gdal_datatype == GDT_CInt32) {
		size_meta = 2;
	} else if(gdal_datatype == GDT_CFloat32) {
		size_meta = 4;
	} else if(gdal_datatype == GDT_CFloat64) {
		size_meta = 8;
	} else {
	}
	int block_num = int(memey_size / (size_meta * block_width));
	block_num	  = block_num > height ? height : block_num; // 行数
	block_num	  = block_num < 1 ? 1 : block_num;
	return block_num;
}

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
ReadComplexMatrixData(int start_line, int width, int line_num,
					  std::shared_ptr<GDALDataset> rasterDataset, GDALDataType gdal_datatype)
{
	int band_num = rasterDataset->GetRasterCount();
	if(gdal_datatype == 0) {
		return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(0, 0);
	} else if(gdal_datatype < 8) {
		if(band_num != 2) {
			return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(0, 0);
		}
	} else if(gdal_datatype < 12) {
		if(band_num != 1) {
			return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(0, 0);
		}

	} else {
	}
	bool																   _flag = false;
	Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> data_mat(
		line_num * width, 2); // 必须强制行优先
	if(gdal_datatype == GDT_Byte) {
		Eigen::MatrixX<char> real_mat(line_num * width, 1);
		Eigen::MatrixX<char> imag_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  real_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // real
		rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  imag_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // imag
		data_mat.col(0) = (real_mat.array().cast<double>()).array();
		data_mat.col(1) = (imag_mat.array().cast<double>()).array();
		_flag			= true;
	} else if(gdal_datatype == GDT_UInt16) {
		Eigen::MatrixX<unsigned short> real_mat(line_num * width, 1);
		Eigen::MatrixX<unsigned short> imag_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  real_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // real
		rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  imag_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // imag
		data_mat.col(0) = (real_mat.array().cast<double>()).array();
		data_mat.col(1) = (imag_mat.array().cast<double>()).array();
		_flag			= true;
	} else if(gdal_datatype == GDT_Int16) {
		Eigen::MatrixX<short> real_mat(line_num * width, 1);
		Eigen::MatrixX<short> imag_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  real_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // real
		rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  imag_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // imag
		data_mat.col(0) = (real_mat.array().cast<double>()).array();
		data_mat.col(1) = (imag_mat.array().cast<double>()).array();
		_flag			= true;
	} else if(gdal_datatype == GDT_UInt32) {
		Eigen::MatrixX<unsigned int> real_mat(line_num * width, 1);
		Eigen::MatrixX<unsigned int> imag_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  real_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // real
		rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  imag_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // imag
		data_mat.col(0) = (real_mat.array().cast<double>()).array();
		data_mat.col(1) = (imag_mat.array().cast<double>()).array();
		_flag			= true;
	} else if(gdal_datatype == GDT_Int32) {
		Eigen::MatrixX<int> real_mat(line_num * width, 1);
		Eigen::MatrixX<int> imag_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  real_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // real
		rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  imag_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // imag
		data_mat.col(0) = (real_mat.array().cast<double>()).array();
		data_mat.col(1) = (imag_mat.array().cast<double>()).array();
		_flag			= true;
	}
	// else if (gdal_datatype == GDT_UInt64) {
	//	Eigen::MatrixX<unsigned long> real_mat(line_num * width, 1);
	//	Eigen::MatrixX<unsigned long> imag_mat(line_num * width, 1);
	//	rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
	//real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
	//	rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
	//imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag 	data_mat.col(0) =
	//(real_mat.array().cast<double>()).array(); 	data_mat.col(1) =
	//(imag_mat.array().cast<double>()).array(); 	_flag = true;
	// }
	// else if (gdal_datatype == GDT_Int64) {
	//	Eigen::MatrixX<long> real_mat(line_num * width, 1);
	//	Eigen::MatrixX<long> imag_mat(line_num * width, 1);
	//	rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
	//real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real
	//	rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
	//imag_mat.data(), width, line_num, gdal_datatype, 0, 0); // imag 	data_mat.col(0) =
	//(real_mat.array().cast<double>()).array(); 	data_mat.col(1) =
	//(imag_mat.array().cast<double>()).array(); 	_flag = true;
	// }
	else if(gdal_datatype == GDT_Float32) {
		Eigen::MatrixX<float> real_mat(line_num * width, 1);
		Eigen::MatrixX<float> imag_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  real_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // real
		rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  imag_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // imag
		data_mat.col(0) = (real_mat.array().cast<double>()).array();
		data_mat.col(1) = (imag_mat.array().cast<double>()).array();
		_flag			= true;
	} else if(gdal_datatype == GDT_Float64) {
		Eigen::MatrixX<double> real_mat(line_num * width, 1);
		Eigen::MatrixX<double> imag_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  real_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // real
		rasterDataset->GetRasterBand(2)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  imag_mat.data(), width, line_num, gdal_datatype,
												  0, 0); // imag
		data_mat.col(0) = (real_mat.array().cast<double>()).array();
		data_mat.col(1) = (imag_mat.array().cast<double>()).array();
		_flag			= true;
	} else if(gdal_datatype == GDT_CInt16) {
		Eigen::MatrixX<std::complex<short>> complex_short_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  complex_short_mat.data(), width, line_num,
												  gdal_datatype, 0, 0); // real
		data_mat.col(0) = (complex_short_mat.real().array().cast<double>()).array();
		data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
		_flag			= true;
	} else if(gdal_datatype == GDT_CInt32) {
		Eigen::MatrixX<std::complex<int>> complex_short_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  complex_short_mat.data(), width, line_num,
												  gdal_datatype, 0, 0); // real
		data_mat.col(0) = (complex_short_mat.real().array().cast<double>()).array();
		data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
		_flag			= true;
	} else if(gdal_datatype == GDT_CFloat32) {
		Eigen::MatrixX<std::complex<double>> complex_short_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  complex_short_mat.data(), width, line_num,
												  gdal_datatype, 0, 0); // real
		data_mat.col(0) = (complex_short_mat.real().array().cast<double>()).array();
		data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
		_flag			= true;
	} else if(gdal_datatype == GDT_CFloat64) {
		Eigen::MatrixX<std::complex<double>> complex_short_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
												  complex_short_mat.data(), width, line_num,
												  gdal_datatype, 0, 0); // real
		data_mat.col(0) = (complex_short_mat.real().array().cast<double>()).array();
		data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
		_flag			= true;
	} else {
	}
	// 保存数据

	if(_flag) {
		return data_mat;
	} else {
		return Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>(
			0, 0); // 必须强制行优先;
	}
}

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
ReadMatrixDoubleData(int start_line, int width, int line_num,
					 std::shared_ptr<GDALDataset> rasterDataset, GDALDataType gdal_datatype,
					 int band_idx)
{
	// 构建矩阵块，使用eigen 进行矩阵计算，加速计算
	bool																   _flag = false;
	Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> data_mat(
		line_num * width, 1); // 必须强制行优先
	if(gdal_datatype == GDT_Byte) {
		Eigen::MatrixX<char> real_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
														 real_mat.data(), width, line_num,
														 gdal_datatype, 0, 0); // real
		data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
		_flag			= true;
	} else if(gdal_datatype == GDT_UInt16) {
		Eigen::MatrixX<unsigned short> real_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
														 real_mat.data(), width, line_num,
														 gdal_datatype, 0, 0); // real
		data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
		_flag			= true;
	} else if(gdal_datatype == GDT_Int16) {
		Eigen::MatrixX<short> real_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
														 real_mat.data(), width, line_num,
														 gdal_datatype, 0, 0); // real
		data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
		_flag			= true;
	} else if(gdal_datatype == GDT_UInt32) {
		Eigen::MatrixX<unsigned int> real_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
														 real_mat.data(), width, line_num,
														 gdal_datatype, 0, 0); // real
		data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
		_flag			= true;
	} else if(gdal_datatype == GDT_Int32) {
		Eigen::MatrixX<int> real_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
														 real_mat.data(), width, line_num,
														 gdal_datatype, 0, 0); // real
		data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
		_flag			= true;
	}
	// else if (gdal_datatype == GDT_UInt64) {
	//	Eigen::MatrixX<unsigned long> real_mat(line_num * width, 1);
	//	rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
	//real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real 	data_mat.col(0) =
	//((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0; 	_flag = true;
	// }
	// else if (gdal_datatype == GDT_Int64) {
	//	Eigen::MatrixX<long> real_mat(line_num * width, 1);
	//	rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
	//real_mat.data(), width, line_num, gdal_datatype, 0, 0); // real 	data_mat.col(0) =
	//((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0; 	_flag = true;
	// }
	else if(gdal_datatype == GDT_Float32) {
		Eigen::MatrixX<float> real_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
														 real_mat.data(), width, line_num,
														 gdal_datatype, 0, 0); // real
		data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
		_flag			= true;
	} else if(gdal_datatype == GDT_Float64) {
		Eigen::MatrixX<double> real_mat(line_num * width, 1);
		rasterDataset->GetRasterBand(band_idx)->RasterIO(GF_Read, 0, start_line, width, line_num,
														 real_mat.data(), width, line_num,
														 gdal_datatype, 0, 0); // real
		data_mat.col(0) = ((real_mat.array().cast<double>()).array().pow(2)).log10() * 10.0;
		_flag			= true;
	} else {
	}

	return data_mat;
}

Eigen::MatrixXd getGeoTranslationArray(QString in_path)
{
	return Eigen::MatrixXd();
}

ImageGEOINFO getImageINFO(QString in_path)
{
	std::shared_ptr<GDALDataset> df		  = OpenDataset(in_path);
	int							 width	  = df->GetRasterXSize();
	int							 heigh	  = df->GetRasterYSize();
	int							 band_num = df->GetRasterCount();
	ImageGEOINFO				 result;
	result.width   = width;
	result.height  = heigh;
	result.bandnum = band_num;

	return result;
}

GDALDataType getGDALDataType(QString fileptah)
{
	omp_lock_t lock;
	omp_init_lock(&lock);
	omp_set_lock(&lock);
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(
		fileptah.toUtf8().constData(), GA_ReadOnly)); // 锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷

	GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();

	GDALClose((GDALDatasetH)rasterDataset);
	omp_unset_lock(&lock);	 // 锟酵放伙拷斤拷
	omp_destroy_lock(&lock); // 劫伙拷斤拷

	return gdal_datatype;
}


int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width,
				 int new_height, GDALResampleAlg eResample)
{
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
	GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile, GA_ReadOnly);
	if(pDSrc == NULL) {
		return -1;
	}

	GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	if(pDriver == NULL) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		return -2;
	}
	int			 width		= pDSrc->GetRasterXSize();
	int			 height		= pDSrc->GetRasterYSize();
	int			 nBandCount = pDSrc->GetRasterCount();
	GDALDataType dataType	= pDSrc->GetRasterBand(1)->GetRasterDataType();

	char*		 pszSrcWKT	= NULL;
	pszSrcWKT				= const_cast<char*>(pDSrc->GetProjectionRef());

	// 锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7
	if(strlen(pszSrcWKT) <= 0) {
		OGRSpatialReference oSRS;
		oSRS.importFromEPSG(4326);
		// oSRS.SetUTM(50, true);   //锟斤拷锟斤拷锟斤拷  锟斤拷锟斤拷120锟斤拷
		// oSRS.SetWellKnownGeogCS("WGS84");
		oSRS.exportToWkt(&pszSrcWKT);
	}
	qDebug() << "GDALCreateGenImgProjTransformer " << Qt::endl;
	void* hTransformArg;
	hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT,
													FALSE, 0.0, 1);
	qDebug() << "no proj " << Qt::endl;
	//(没锟斤拷投影锟斤拷影锟斤拷锟斤拷锟斤拷锟斤拷卟锟酵?拷锟?1锟?7)
	if(hTransformArg == NULL) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		return -3;
	}
	qDebug() << "has proj " << Qt::endl;
	double dGeoTrans[6] = { 0 };
	int	   nNewWidth = 0, nNewHeight = 0;
	if(GDALSuggestedWarpOutput((GDALDatasetH)pDSrc, GDALGenImgProjTransform, hTransformArg,
							   dGeoTrans, &nNewWidth, &nNewHeight)
	   != CE_None) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		return -3;
	}

	// GDALDestroyGenImgProjTransformer(hTransformArg);
	qDebug() << "create init GDALDataset " ;
	GDALDataset* pDDst =
		pDriver->Create(pszOutFile, new_width, new_height, nBandCount, dataType, NULL);
	if(pDDst == NULL) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		return -2;
	}

	pDDst->SetProjection(pszSrcWKT);
	pDDst->SetGeoTransform(gt);

	GDALWarpOptions* psWo  = GDALCreateWarpOptions();
	CPLSetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS"); // 使用所有可用的CPU核心
	CPLSetConfigOption("GDAL_CACHEMAX", "16000");         // 设置缓存大小为500MB
	// psWo->papszWarpOptions = CSLDuplicate(NULL);
	psWo->eWorkingDataType = dataType;
	psWo->eResampleAlg	   = eResample;

	psWo->hSrcDS		   = (GDALDatasetH)pDSrc;
	psWo->hDstDS		   = (GDALDatasetH)pDDst;
	qDebug() << "GDALCreateGenImgProjTransformer" << Qt::endl;
	psWo->pfnTransformer  = GDALGenImgProjTransform;
	psWo->pTransformerArg = GDALCreateGenImgProjTransformer(
		(GDALDatasetH)pDSrc, pszSrcWKT, (GDALDatasetH)pDDst, pszSrcWKT, FALSE, 0.0, 1);
	;

	qDebug() << "GDALCreateGenImgProjTransformer has created" << Qt::endl;
	psWo->nBandCount  = nBandCount;
	psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int));
	psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int));
	for(int i = 0; i < nBandCount; i++) {
		psWo->panSrcBands[i] = i + 1;
		psWo->panDstBands[i] = i + 1;
	}

	GDALWarpOperation oWo;
	if(oWo.Initialize(psWo) != CE_None) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
		return -3;
	}
	qDebug() << "ChunkAndWarpImage:" << new_width << "," << new_height << Qt::endl;
	oWo.ChunkAndWarpMulti(0, 0, new_width, new_height);
	GDALFlushCache(pDDst);
	qDebug() << "ChunkAndWarpImage over" << Qt::endl;
	// GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
	// GDALDestroyWarpOptions(psWo);
	GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
	GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
	GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	return 0;
}

int ResampleGDALs(const char* pszSrcFile, int band_ids, GDALRIOResampleAlg eResample)
{
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
	GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile, GA_Update);
	if(pDSrc == NULL) {
		return -1;
	}

	GDALDataType		 gdal_datatype = pDSrc->GetRasterBand(1)->GetRasterDataType();

	GDALRasterBand*		 demBand	   = pDSrc->GetRasterBand(band_ids);

	int					 width		   = pDSrc->GetRasterXSize();
	int					 height		   = pDSrc->GetRasterYSize();
	int					 start_col = 0, start_row = 0, rows_count = 0, cols_count;

	int					 row_delta = int(120000000 / width);

	GDALRasterIOExtraArg psExtraArg;
	INIT_RASTERIO_EXTRA_ARG(psExtraArg);
	psExtraArg.eResampleAlg = eResample;

	do {
		rows_count = start_row + row_delta < height ? row_delta : height - start_row;
		cols_count = width;

		if(gdal_datatype == GDALDataType::GDT_UInt16) {
			unsigned short* temp = new unsigned short[rows_count * cols_count];
			demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp,
							  cols_count, rows_count, gdal_datatype, 0, 0);
			demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp,
							  cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg);
			delete[] temp;

		} else if(gdal_datatype == GDALDataType::GDT_Int16) {
			short* temp = new short[rows_count * cols_count];
			demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp,
							  cols_count, rows_count, gdal_datatype, 0, 0);
			demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp,
							  cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg);
			delete[] temp;
		} else if(gdal_datatype == GDALDataType::GDT_Float32) {
			float* temp = new float[rows_count * cols_count];
			demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp,
							  cols_count, rows_count, gdal_datatype, 0, 0);
			demBand->RasterIO(GF_Write, start_col, start_row, cols_count, rows_count, temp,
							  cols_count, rows_count, gdal_datatype, 0, 0, &psExtraArg);
			delete[] temp;
		}
		start_row = start_row + rows_count;
	} while(start_row < height);
	GDALClose((GDALDatasetH)pDSrc);

	return 0;
}

int alignRaster(QString inputPath, QString referencePath, QString outputPath, GDALResampleAlg eResample)
{

	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
 
	GDALDataset* pDSrc = (GDALDataset*)GDALOpen(inputPath.toUtf8().constData(), GA_ReadOnly);
	if (pDSrc == NULL) {
		return -1;
	}

	GDALDataset* tDSrc = (GDALDataset*)GDALOpen(referencePath.toUtf8().constData(), GA_ReadOnly);
	if (tDSrc == NULL) {
		return -1;
	}

	long new_width = tDSrc->GetRasterXSize();
	long new_height = tDSrc->GetRasterYSize();

	GDALDriver* pDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (pDriver == NULL) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		return -2;
	}
	int			 width = pDSrc->GetRasterXSize();
	int			 height = pDSrc->GetRasterYSize();
	int			 nBandCount = pDSrc->GetRasterCount();
	GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType(); 

	char* pszSrcWKT = NULL;
	pszSrcWKT = const_cast<char*>(pDSrc->GetProjectionRef());

	if (strlen(pszSrcWKT) <= 0) {
		OGRSpatialReference oSRS;
		oSRS.importFromEPSG(4326);
		oSRS.exportToWkt(&pszSrcWKT);
	}
	qDebug() << "GDALCreateGenImgProjTransformer " << Qt::endl;
	void* hTransformArg;
	hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT,FALSE, 0.0, 1);
	qDebug() << "no proj " << Qt::endl;
	if (hTransformArg == NULL) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		return -3;
	}
	qDebug() << "has proj " << Qt::endl;
	std::shared_ptr<double> dGeoTrans(new double[6], delArrPtr);
	int	   nNewWidth = 0,
		   nNewHeight = 0;

	if (GDALSuggestedWarpOutput((GDALDatasetH)pDSrc, GDALGenImgProjTransform, hTransformArg,dGeoTrans.get(), &nNewWidth, &nNewHeight) != CE_None) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		return -3;
	}

	GDALDataset* pDDst =	pDriver->Create(outputPath.toUtf8().constData(), new_width, new_height, pDSrc->GetRasterCount(), dataType, NULL);
	if (pDDst == NULL) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		return -2;
	}

	std::shared_ptr<double> gt(new double[6], delArrPtr);
	tDSrc->GetGeoTransform(gt.get());
	pDDst->SetProjection(pszSrcWKT);
	pDDst->SetGeoTransform(gt.get());

	GDALWarpOptions* psWo = GDALCreateWarpOptions();

	psWo->eWorkingDataType = dataType;
	psWo->eResampleAlg = eResample;

	psWo->hSrcDS = (GDALDatasetH)pDSrc;
	psWo->hDstDS = (GDALDatasetH)pDDst;

 	psWo->pfnTransformer = GDALGenImgProjTransform;
	psWo->pTransformerArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, (GDALDatasetH)pDDst, pszSrcWKT, FALSE, 0.0, 1);
	

	qDebug() << "GDALCreateGenImgProjTransformer has created" << Qt::endl;
	psWo->nBandCount = nBandCount;
	psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int));
	psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int));
	for (int i = 0; i < nBandCount; i++) {
		psWo->panSrcBands[i] = i + 1;
		psWo->panDstBands[i] = i + 1;
	}

	GDALWarpOperation oWo;
	if (oWo.Initialize(psWo) != CE_None) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
		GDALClose((GDALDatasetH)(GDALDatasetH)tDSrc);
		return -3;
	}
	qDebug() << "ChunkAndWarpImage:" << new_width << "," << new_height << Qt::endl;
	oWo.ChunkAndWarpMulti(0, 0, new_width, new_height);
	GDALFlushCache(pDDst);
	qDebug() << "ChunkAndWarpImage over" << Qt::endl;
	GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
	GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
	GDALClose((GDALDatasetH)(GDALDatasetH)tDSrc);
	return 0;


}


int saveMatrixXcd2TiFF(Eigen::MatrixXcd data, QString out_tiff_path)
{
	int				rows = data.rows();
	int				cols = data.cols();

	Eigen::MatrixXd gt	 = Eigen::MatrixXd::Zero(2, 3);

	gdalImage		image_tiff =
		CreategdalImage(out_tiff_path, rows, cols, 2, gt, "", false, true); // 注意这里保留仿真结果
	// 保存二进制文件
	Eigen::MatrixXd real_img = data.array().real();
	Eigen::MatrixXd imag_img = data.array().imag();
	image_tiff.saveImage(real_img, 0, 0, 1);
	image_tiff.saveImage(imag_img, 0, 0, 2);
	return -1;
}

void clipRaster(QString inRasterPath, QString outRasterPath, long minRow, long maxRow, long minCol, long maxCol)
{
	long rownum = maxRow - minRow + 1;
	long colnum = maxCol - minCol + 1;

	gdalImage inimg(inRasterPath);
	Eigen::MatrixXd gt = inimg.gt;
	
	Landpoint lp = inimg.getLandPoint(minRow, minCol, 0);

	gt(0, 0) = lp.lon;
	gt(1, 0) = lp.lat;

	gdalImage outimg= CreategdalImageDouble(outRasterPath, rownum, colnum, inimg.band_num, gt, inimg.projection, true, true, true);

	for (long bi = 1; bi < inimg.band_num + 1; bi++) {
		Eigen::MatrixXd brasterData = inimg.getData(minRow, minCol, rownum, colnum, bi);
		outimg.saveImage(brasterData, 0, 0, bi);
		qDebug() << "writer raster band : " << bi;
	}

	qDebug() << "writer raster overring";

}



 
long getProjectEPSGCodeByLon_Lat(double long, double lat, ProjectStripDelta stripState)
{
	long EPSGCode = 0;
	switch (stripState) {
	case ProjectStripDelta::Strip_3: {
		break;
	};
	case ProjectStripDelta::Strip_6: {
		break;
	}
	default: {
		EPSGCode = -1;
		break;
	}
	}
	qDebug() << QString(" EPSG code : %1").arg(EPSGCode);
	return EPSGCode;
}
long getProjectEPSGCodeByLon_Lat_inStrip3(double lon, double lat)
{
	// EPSG 4534 ~ 4554 3 度带
	// 首先判断是否是在 中国带宽范围
	// 中心经度范围 ：75E ~ 135E 实际范围 73.5E ~ 136.5E,
	// 纬度范围 3N ~ 54N，放宽到 0N~ 60N
	if (73.5 <= lon && lon <= 136.5 && 0 <= lat && lat <= 60) { // 中国境内
		long code = trunc((lon - 73.5) / 3) + 4534;
		return code;
	}
	else { // 非中国境内 使用 高斯克吕格
		bool isSouth = lat < 0; // 简单判断南北半球，这里仅为示例，实际应用可能需要更细致的逻辑
		long prefix = isSouth ? 327000 : 326000;
		// std::string prefix = isSouth ? "327" : "326";
		lon = fmod(lon + 360.0, 360.0);
		long zone = std::floor((lon + 180.0) / 3.0);
		prefix = prefix + zone;
		return prefix;
	}
	return 0;
}
long getProjectEPSGCodeByLon_Lat_inStrip6(double lon, double lat)
{
	// EPSG 4502 ~ 4512   6度带
	// 首先判断是否是在 中国带宽范围
	// 中心经度范围 ：75E ~ 135E 实际范围 72.0E ~ 138E,
	// 纬度范围 3N ~ 54N，放宽到 0N~ 60N
	if (73.5 <= lon && lon <= 136.5 && 0 <= lat && lat <= 60) { // 中国境内
		long code = trunc((lon - 72.0) / 6) + 4502;
		return code;
	}
	else { // 非中国境内 使用 UTM// 确定带号，6度带从1开始到60，每6度一个带
		int	   zone = static_cast<int>((lon + 180.0) / 6.0) + 1;
		bool   isSouth = lon < 0;				   // 判断是否在南半球
		long epsgCodeBase = isSouth ? 32700 : 32600; // 计算EPSG代码
		long prefix = static_cast<int>(epsgCodeBase + zone);
		return prefix;
	}
	return 0;
}

QString GetProjectionNameFromEPSG(long epsgCode)
{
	qDebug() << "============= GetProjectionNameFromEPSG ======================";
	OGRSpatialReference oSRS;

	// 设置EPSG代码
	if (oSRS.importFromEPSG(epsgCode) != OGRERR_NONE) {
		qDebug() << "epsgcode not  recognition";
		return "";
	}

	// 获取并输出坐标系的描述（名称）
	const char* pszName = oSRS.GetAttrValue("GEOGCS");
	if (pszName) {
		qDebug() << "Coordinate system name for EPSG " + QString::number(epsgCode)
			<< " is: " + QString::fromStdString(pszName);
		return QString::fromStdString(pszName);
	}
	else {
		qDebug() << "Unable to retrieve the name for EPSG " + QString::number(epsgCode);
		return "";
	}

	//		char*				wkt = NULL;
	//		// 转换为WKT格式
	//		oSRS.exportToWkt(&wkt);
	//
	//		qDebug() << wkt;
	//
	//		// 从WKT中解析投影名称，这里简化处理，实际可能需要更复杂的逻辑来准确提取名称
	//		std::string wktStr(wkt);
	//		long start = wktStr.find("PROJCS[\"") + 8; // 找到"PROJCS["后的第一个双引号位置
	//		// 从start位置开始找下一个双引号，这之间的内容即为投影名称
	//		int	end		 = wktStr.find('\"', start);
	//		QString projName = QString::fromStdString(wktStr.substr(start, end - start));
	//
	//		// 释放WKT字符串内存
	//		CPLFree(wkt);

	// return projName;
}



long GetEPSGFromRasterFile(QString filepath)
{
	qDebug() << "============= GetEPSGFromRasterFile ======================";
	//		QTextCodec* codec = QTextCodec::codecForLocale(); // 获取系统默认编码的文本编解码器
	//		QByteArray	byteArray = codec->fromUnicode(filepath); // 将QString转换为QByteArray
	//,这个应该会自动释放 		const char* charArray = byteArray.constData(); //
	// 获取QByteArray的const char*指针

	{
		if (QFile(filepath).exists()) {
			qDebug() << "info: the image found.\n";
		}
		else {
			return -1;
		}

		GDALAllRegister();
		CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注册GDAL驱动
		//			qDebug()<<filepath.toLocal8Bit().constData()<<std::endl;
		// 打开影像文件
		GDALDataset* poDataset;
		poDataset = (GDALDataset*)GDALOpen(filepath.toUtf8().data(), GA_ReadOnly);
		if (NULL == poDataset) {
			qDebug() << "Error: Unable to open the image.\n";
			return -1;
		}

		// 获取影像的投影信息
		const char* pszProjection = poDataset->GetProjectionRef();

		qDebug() << QString::fromUtf8(pszProjection);

		// 创建SpatialReference对象
		OGRSpatialReference oSRS;
		if (oSRS.importFromWkt((char**)&pszProjection) != OGRERR_NONE) {
			qDebug() << ("Error: Unable to import projection information.\n");
			GDALClose(poDataset);
			return -1;
		}

		qDebug() << pszProjection ;

		const char* epscodestr = oSRS.GetAuthorityCode(nullptr);
		if (NULL == epscodestr || nullptr == epscodestr) {
			qDebug() << "EPSG code  string could not be determined from the spatial reference.";
			GDALClose(poDataset);

			return -1;
		}
		long epsgCode = atoi(epscodestr); // 获取EPSG代码

		if (epsgCode != 0) {
			GDALClose(poDataset);
			qDebug() << QString("file %1 :epsg Code %2").arg(filepath).arg(epsgCode);
			return epsgCode;
		}
		else {
			qDebug() << "EPSG code could not be determined from the spatial reference.";
			GDALClose(poDataset);

			return -1;
		}
	}
}

std::shared_ptr<PointRaster> GetCenterPointInRaster(QString filepath)
{
	qDebug() << "============= GetCenterPointInRaster ======================";
	//		QTextCodec* codec = QTextCodec::codecForLocale(); // 获取系统默认编码的文本编解码器
	//		QByteArray	byteArray = codec->fromUnicode(filepath); // 将QString转换为QByteArray
	//,这个应该会自动释放 		const char* charArray = byteArray.constData(); //
	// 获取QByteArray的const char*指针
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	// qDebug()<<filepath.toLocal8Bit().constData()<<std::endl;
	GDALDataset* poDataset = (GDALDataset*)GDALOpen(filepath.toUtf8().data(), GA_ReadOnly);
	if (nullptr == poDataset || NULL == poDataset) {
		qDebug() << "Could not open dataset";
		return nullptr; // 表示读取失败
	}

	double x, y, z;

	bool   flag = false;

	{
		double adfGeoTransform[6];
		if (poDataset->GetGeoTransform(adfGeoTransform) != CE_None) {
			qDebug() << "Failed to get GeoTransform";
			return nullptr;
		}

		double dfWidth = poDataset->GetRasterXSize();
		double dfHeight = poDataset->GetRasterYSize();

		// 计算中心点坐标（像素坐标）
		double dfCenterX = adfGeoTransform[0] + dfWidth * adfGeoTransform[1] / 2.0
			+ dfHeight * adfGeoTransform[2] / 2.0;
		double dfCenterY = adfGeoTransform[3] + dfWidth * adfGeoTransform[4] / 2.0
			+ dfHeight * adfGeoTransform[5] / 2.0;

		OGRSpatialReference oSRS;
		oSRS.importFromWkt(poDataset->GetProjectionRef());

		if (oSRS.IsGeographic()) {
			qDebug() << "Center coords (already in geographic): (" + QString::number(dfCenterX)
				+ ", " + QString::number(dfCenterY) + ")";
			flag = true;
			x = dfCenterX;
			y = dfCenterY;
		}
		else {
			// 如果不是地理坐标系，转换到WGS84
			OGRSpatialReference oSRS_WGS84;
			oSRS_WGS84.SetWellKnownGeogCS("WGS84");

			OGRCoordinateTransformation* poCT =
				OGRCreateCoordinateTransformation(&oSRS, &oSRS_WGS84);
			if (poCT == nullptr) {
				qDebug() << "Failed to create coordinate transformation";
				return nullptr;
			}

			// double dfLon, dfLat;
			if (poCT->Transform(1, &dfCenterX, &dfCenterY)) {
				qDebug() << "Center coords (transformed to WGS84): ("
					+ QString::number(dfCenterX) + ", " + QString::number(dfCenterY)
					<< ")";
				flag = true;
				x = dfCenterX;
				y = dfCenterY;
			}
			else {
				qDebug() << "Transformation failed.";
			}
			OGRCoordinateTransformation::DestroyCT(poCT);
		}
	}
	if (nullptr == poDataset || NULL == poDataset) {}
	else {
		GDALClose(poDataset);
	}

	if (flag) {
		std::shared_ptr<PointRaster> RasterCenterPoint = std::make_shared<PointRaster>();
		RasterCenterPoint->x = x;
		RasterCenterPoint->y = y;
		RasterCenterPoint->z = 0;
		return RasterCenterPoint;
	}
	else {
		return nullptr;
	}
}
CoordinateSystemType getCoordinateSystemTypeByEPSGCode(long epsg_code)
{
	OGRSpatialReference oSRS;
	if (oSRS.importFromEPSG(epsg_code) == OGRERR_NONE) {
		if (oSRS.IsGeographic()) {
			return CoordinateSystemType::GeoCoordinateSystem;
		}
		else if (oSRS.IsProjected()) {
			return CoordinateSystemType::ProjectCoordinateSystem;
		}
		else {
			return CoordinateSystemType::UNKNOW;
		}
	}
	else {
		return CoordinateSystemType::UNKNOW;
	}
}


void resampleRaster(const char* inputRaster, const char* outputRaster, double targetPixelSizeX, double targetPixelSizeY) {
	// 初始化GDAL
	GDALAllRegister();
	// 打开输入栅格文件
	GDALDataset* poDataset = (GDALDataset*)GDALOpen(inputRaster, GA_ReadOnly);
	if (poDataset == nullptr) {
		std::cerr << "Failed to open raster file." << std::endl;
		return;
	}

	// 获取原始栅格的空间参考
	double adfGeoTransform[6];
	if (poDataset->GetGeoTransform(adfGeoTransform) != CE_None) {
		std::cerr << "Failed to get geotransform." << std::endl;
		GDALClose(poDataset);
		return;
	}

	// 获取原始栅格的尺寸
	int nXSize = poDataset->GetRasterXSize();
	int nYSize = poDataset->GetRasterYSize();

	// 计算目标栅格的尺寸
	double targetXSize = (adfGeoTransform[1] * nXSize) / targetPixelSizeX;
	double targetYSize = (adfGeoTransform[5] * nYSize) / targetPixelSizeY;

	// 创建目标数据集（输出栅格）
	GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (poDriver == nullptr) {
		std::cerr << "Failed to get GTiff driver." << std::endl;
		GDALClose(poDataset);
		return;
	}

	// 创建输出数据集
	GDALDataset* poOutDataset = poDriver->Create(outputRaster, (int)targetXSize, (int)targetYSize, poDataset->GetRasterCount(), GDT_Float32, nullptr);
	if (poOutDataset == nullptr) {
		std::cerr << "Failed to create output raster." << std::endl;
		GDALClose(poDataset);
		return;
	}

	// 设置输出数据集的地理变换（坐标系）
	double targetGeoTransform[6] = {
		adfGeoTransform[0], targetPixelSizeX, 0, adfGeoTransform[3], 0, -targetPixelSizeY
	};
	poOutDataset->SetGeoTransform(targetGeoTransform);

	// 设置输出数据集的投影信息
	poOutDataset->SetProjection(poDataset->GetProjectionRef());

	// 进行重采样
	for (int i = 0; i < poDataset->GetRasterCount(); i++) {
		GDALRasterBand* poBand = poDataset->GetRasterBand(i + 1);
		GDALRasterBand* poOutBand = poOutDataset->GetRasterBand(i + 1);

		// 使用GDAL的重采样方法，选择一个适当的重采样方式
		poOutBand->RasterIO(GF_Write, 0, 0, (int)targetXSize, (int)targetYSize,
			nullptr, (int)targetXSize, (int)targetYSize,
			poBand->GetRasterDataType(), 0, 0, nullptr);
	}

	// 关闭数据集
	GDALClose(poDataset);
	GDALClose(poOutDataset);

	qDebug() << "Resampling completed." ;
}

void transformRaster(const char* inputFile, const char* outputFile, int sourceEPSG, int targetEPSG) {
	// 初始化 GDAL 库
	GDALAllRegister();

	// 打开源栅格文件
	GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(inputFile, GA_ReadOnly);
	if (poSrcDS == nullptr) {
		qDebug() << "Failed to open input file:" << inputFile;
		return;
	}

	// 获取源栅格的基本信息
	int nXSize = poSrcDS->GetRasterXSize();
	int nYSize = poSrcDS->GetRasterYSize();
	int nBands = poSrcDS->GetRasterCount();
	GDALDataType eDT = poSrcDS->GetRasterBand(1)->GetRasterDataType();

	// 创建目标栅格文件
	GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (poDriver == nullptr) {
		qDebug() << "GTiff driver not available.";
		GDALClose(poSrcDS);
		return;
	}

	GDALDataset* poDstDS = poDriver->Create(outputFile, nXSize, nYSize, nBands, eDT, nullptr);
	if (poDstDS == nullptr) {
		qDebug() << "Failed to create output file:" << outputFile;
		GDALClose(poSrcDS);
		return;
	}

	// 设置目标栅格的空间参考系统
	OGRSpatialReference oSRS;
	oSRS.importFromEPSG(targetEPSG);

	char* pszWKT = nullptr;
	oSRS.exportToWkt(&pszWKT);
	poDstDS->SetProjection(pszWKT);
	CPLFree(pszWKT);

	// 复制元数据
	poDstDS->SetMetadata(poSrcDS->GetMetadata());

	// 复制每个波段的数据
	for (int i = 1; i <= nBands; ++i) {
		GDALRasterBand* poSrcBand = poSrcDS->GetRasterBand(i);
		GDALRasterBand* poDstBand = poDstDS->GetRasterBand(i);

		float* pafScanline = (float*)CPLMalloc(sizeof(float) * nXSize);

		for (int j = 0; j < nYSize; ++j) {
			poSrcBand->RasterIO(GF_Read, 0, j, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);
			poDstBand->RasterIO(GF_Write, 0, j, nXSize, 1, pafScanline, nXSize, 1, GDT_Float32, 0, 0);
		}

		CPLFree(pafScanline);
	}

	// 关闭数据集
	GDALClose(poSrcDS);
	GDALClose(poDstDS);

	qDebug() << "Raster transformation completed successfully.";
}



ErrorCode transformCoordinate(double x, double y, int sourceEPSG, int targetEPSG, Point2 &p) {
	// 创建源坐标系（原始坐标系）
	OGRSpatialReference sourceSRS;
	sourceSRS.importFromEPSG(sourceEPSG);  // 使用 EPSG 编码来指定坐标系

	// 创建目标坐标系（目标坐标系）
	OGRSpatialReference targetSRS;
	targetSRS.importFromEPSG(targetEPSG);  // WGS84 坐标系 EPSG:4326

	// 创建坐标变换对象
	OGRCoordinateTransformation* transform = OGRCreateCoordinateTransformation(&sourceSRS, &targetSRS);
	if (transform == nullptr) {
		qDebug() << "Failed to create coordinate transformation." ;
		return ErrorCode::FAIL;
	}

	// 转换坐标
	double transformedX = x;
	double transformedY = y;
	if (transform->Transform(1, &transformedX, &transformedY)) {
		qDebug() << "Original Coordinates: (" << x << ", " << y << ")" ;
		qDebug() << "Transformed Coordinates (EPSG:" << targetEPSG << "): (" << transformedX << ", " << transformedY << ")" ;
	}
	else {
		qDebug() << "Coordinate transformation failed." ;
	}

	// 清理
	delete transform;
	p.x = transformedX;
	p.y = transformedY;
	return ErrorCode::SUCCESS;

}

void cropRasterByLatLon(const char* inputFile, const char* outputFile,double minLon, double maxLon, double minLat, double maxLat) {
	// 初始化 GDAL 库
	GDALAllRegister();

	// 打开栅格数据集
	GDALDataset* poDataset = (GDALDataset*)GDALOpen(inputFile, GA_ReadOnly);
	if (poDataset == nullptr) {
		std::cerr << "Failed to open input raster." << std::endl;
		return;
	}

	// 获取栅格数据的地理参考信息
	double adfGeoTransform[6];
	if (poDataset->GetGeoTransform(adfGeoTransform) != CE_None) {
		std::cerr << "Failed to get geotransform." << std::endl;
		GDALClose(poDataset);
		return;
	}

	// 获取输入影像的投影信息
	const char* projection = poDataset->GetProjectionRef();

	// 根据经纬度计算出裁剪区域对应的栅格像素坐标
	int xMin = (int)((minLon - adfGeoTransform[0]) / adfGeoTransform[1]);
	int xMax = (int)((maxLon - adfGeoTransform[0]) / adfGeoTransform[1]);
	int yMin = (int)((maxLat - adfGeoTransform[3]) / adfGeoTransform[5]);
	int yMax = (int)((minLat - adfGeoTransform[3]) / adfGeoTransform[5]);

	// 创建裁剪区域的目标栅格数据集
	GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (poDriver == nullptr) {
		std::cerr << "Failed to get GTiff driver." << std::endl;
		GDALClose(poDataset);
		return;
	}

	// 创建输出栅格数据集，指定尺寸
	int width = xMax - xMin;
	int height = yMax - yMin;
	GDALDataset* poOutDataset = poDriver->Create(outputFile, width, height, poDataset->GetRasterCount(), GDT_Float32, nullptr);
	if (poOutDataset == nullptr) {
		std::cerr << "Failed to create output raster." << std::endl;
		GDALClose(poDataset);
		return;
	}

	// 设置输出栅格的投影信息和地理变换
	poOutDataset->SetProjection(projection);
	double newGeoTransform[6] = { adfGeoTransform[0] + xMin * adfGeoTransform[1], adfGeoTransform[1], 0.0, adfGeoTransform[3] + yMin * adfGeoTransform[5], 0.0, adfGeoTransform[5] };
	poOutDataset->SetGeoTransform(newGeoTransform);

	// 循环读取源数据并写入目标数据集
	for (int i = 0; i < poDataset->GetRasterCount(); ++i) {
		GDALRasterBand* poBand = poDataset->GetRasterBand(i + 1);
		GDALRasterBand* poOutBand = poOutDataset->GetRasterBand(i + 1);

		// 读取源数据
		int* pData = new int[width * height];
		poBand->RasterIO(GF_Read, xMin, yMin, width, height, pData, width, height, GDT_Int32, 0, 0);

		// 写入目标数据
		poOutBand->RasterIO(GF_Write, 0, 0, width, height, pData, width, height, GDT_Int32, 0, 0);

		delete[] pData;
	}

	qDebug() << "Raster cropped and saved to: " << outputFile ;

	// 清理
	GDALClose(poDataset);
	GDALClose(poOutDataset);
}


ErrorCode  DEM2XYZRasterAndSlopRaster(QString dempath, QString demxyzpath, QString demsloperPath)
{
	gdalImage demds(dempath);
	gdalImage demxyz = CreategdalImage(demxyzpath, demds.height, demds.width, 3, demds.gt, demds.projection, true, true);// X,Y,Z

	// 分块计算并转换为XYZ

	Eigen::MatrixXd demArr = demds.getData(0, 0, demds.height, demds.width, 1);
	Eigen::MatrixXd demR = demArr;
	Landpoint LandP{ 0,0,0 };
	Point3 GERpoint{ 0,0,0 };
	double R = 0;
	double dem_row = 0, dem_col = 0, dem_alt = 0;

	long line_invert = 1000;

	double rowidx = 0;
	double colidx = 0;
	for (int max_rows_ids = 0; max_rows_ids < demds.height; max_rows_ids = max_rows_ids + line_invert) {
		Eigen::MatrixXd demdata = demds.getData(max_rows_ids, 0, line_invert, demds.width, 1);
		Eigen::MatrixXd xyzdata_x = demdata.array() * 0;
		Eigen::MatrixXd xyzdata_y = demdata.array() * 0;
		Eigen::MatrixXd xyzdata_z = demdata.array() * 0;

		int datarows = demdata.rows();
		int datacols = demdata.cols();

		for (int i = 0; i < datarows; i++) {
			for (int j = 0; j < datacols; j++) {
				rowidx = i + max_rows_ids;
				colidx = j;
				demds.getLandPoint(rowidx, colidx, demdata(i, j), LandP); // 获取地理坐标
				LLA2XYZ(LandP, GERpoint); // 经纬度转换为地心坐标系
				xyzdata_x(i, j) = GERpoint.x;
				xyzdata_y(i, j) = GERpoint.y;
				xyzdata_z(i, j) = GERpoint.z;
			}
		}
		demxyz.saveImage(xyzdata_x, max_rows_ids, 0, 1);
		demxyz.saveImage(xyzdata_y, max_rows_ids, 0, 2);
		demxyz.saveImage(xyzdata_z, max_rows_ids, 0, 3);
	}


	// 计算坡向角
	gdalImage demsloperxyz = CreategdalImage(demsloperPath, demds.height, demds.width, 4, demds.gt, demds.projection, true, true);// X,Y,Z,cosangle

	line_invert = 1000;
	long start_ids = 0;
	long dem_rows = 0, dem_cols = 0;

	for (start_ids = 1; start_ids < demds.height; start_ids = start_ids + line_invert) {
		Eigen::MatrixXd demdata = demds.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 1);
		long startlineid = start_ids;
		Eigen::MatrixXd demsloper_x = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 1);
		Eigen::MatrixXd demsloper_y = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 2);
		Eigen::MatrixXd demsloper_z = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 3);
		Eigen::MatrixXd demsloper_angle = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 4);

		Landpoint p0, p1, p2, p3, p4, pslopeVector, pp;
		Vector3D  slopeVector;

		dem_rows = demsloper_y.rows();
		dem_cols = demsloper_y.cols();
		double sloperAngle = 0;
		Vector3D Zaxis = { 0,0,1 };

		double rowidx = 0, colidx = 0;

		for (long i = 1; i < dem_rows - 1; i++) {
			for (long j = 1; j < dem_cols - 1; j++) {
				rowidx = i + startlineid;
				colidx = j;
				demds.getLandPoint(rowidx, colidx, demdata(i, j), p0);
				demds.getLandPoint(rowidx - 1, colidx, demdata(i - 1, j), p1);
				demds.getLandPoint(rowidx, colidx - 1, demdata(i, j - 1), p2);
				demds.getLandPoint(rowidx + 1, colidx, demdata(i + 1, j), p3);
				demds.getLandPoint(rowidx, colidx + 1, demdata(i, j + 1), p4);

				pslopeVector = getSlopeVector(p0, p1, p2, p3, p4); // 地面坡向矢量
				slopeVector = { pslopeVector.lon,pslopeVector.lat,pslopeVector.ati };
				pp = LLA2XYZ(p0);
				Zaxis.x = pp.lon;
				Zaxis.y = pp.lat;
				Zaxis.z = pp.ati;
				sloperAngle = getCosAngle(slopeVector, Zaxis); // 地面坡向角

				demsloper_x(i, j) = slopeVector.x;
				demsloper_y(i, j) = slopeVector.y;
				demsloper_z(i, j) = slopeVector.z;
				demsloper_angle(i, j) = sloperAngle;
			}
		}
		demsloperxyz.saveImage(demsloper_x, start_ids - 1, 0, 1);
		demsloperxyz.saveImage(demsloper_y, start_ids - 1, 0, 2);
		demsloperxyz.saveImage(demsloper_z, start_ids - 1, 0, 3);
		demsloperxyz.saveImage(demsloper_angle, start_ids - 1, 0, 4);
	}


	return ErrorCode::SUCCESS;
}


void ConvertCoordinateSystem(QString inRasterPath, QString outRasterPath, long outepsgcode) {
	// 注册所有GDAL驱动
	GDALAllRegister();

	// 打开输入栅格文件
	GDALDataset* srcDataset = (GDALDataset*)GDALOpen(inRasterPath.toUtf8().constData(), GA_ReadOnly);
	if (!srcDataset) {
		// 错误处理：输出文件打开失败
		// qDebug() << "无法打开输入文件：" << inRasterPath;
		return;
	}

	// 创建目标坐标系
	OGRSpatialReference targetSRS;
	if (targetSRS.importFromEPSG(outepsgcode) != OGRERR_NONE) {
		GDALClose(srcDataset);
		// qDebug() << "无效的EPSG代码：" << outepsgcode;
		return;
	}
	GDALDataType datetype = srcDataset->GetRasterBand(1)->GetRasterDataType();
	// 获取目标坐标系的WKT表示
	char* targetSRSWkt = nullptr;
	targetSRS.exportToWkt(&targetSRSWkt);

	bool flag = (datetype == GDT_Byte ||	datetype == GDT_Int8 ||	datetype == GDT_Int16 ||datetype == GDT_UInt16 ||		datetype == GDT_Int32 ||		datetype == GDT_UInt32 ||		datetype == GDT_Int64 ||		datetype == GDT_UInt64);

	// 创建重投影后的虚拟数据集（Warped VRT）
	GDALDataset* warpedVRT = flag? (GDALDataset*)GDALAutoCreateWarpedVRT(
		srcDataset,
		nullptr,          // 输入坐标系（默认使用源数据）
		targetSRSWkt,     // 目标坐标系
		GRA_NearestNeighbour,     // 重采样方法：双线性插值
		0.0,              // 最大误差（0表示自动计算）
		nullptr           // 其他选项
	) :(GDALDataset*)GDALAutoCreateWarpedVRT(
		srcDataset,
		nullptr,          // 输入坐标系（默认使用源数据）
		targetSRSWkt,     // 目标坐标系
		GRA_Bilinear,     // 重采样方法：双线性插值
		0.0,              // 最大误差（0表示自动计算）
		nullptr           // 其他选项
	);
	CPLFree(targetSRSWkt); // 释放WKT内存

	if (!warpedVRT) {
		GDALClose(srcDataset);
		qDebug() << u8"创建投影转换VRT失败";
		return;
	}

	// 获取输出驱动（GeoTIFF格式）
	QString filesuffer = getFileExtension(outRasterPath).toLower();
	bool isTiff = filesuffer.contains("tif");
	GDALDriver* driver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");
	if (!driver) {
		GDALClose(warpedVRT);
		GDALClose(srcDataset);
		// qDebug() << "无法获取GeoTIFF驱动";
		return;
	}

	// 创建输出栅格文件
	GDALDataset* dstDataset = driver->CreateCopy(
		outRasterPath.toUtf8().constData(), // 输出文件路径
		warpedVRT,                          // 输入数据集（VRT）
		false,                              // 是否严格复制
		nullptr,                            // 创建选项
		nullptr,                            // 进度回调
		nullptr                             // 回调参数
	);

	if (!dstDataset) {
		// qDebug() << "创建输出文件失败：" << outRasterPath;
		GDALClose(warpedVRT);
		GDALClose(srcDataset);
		return;
	}

	// 释放资源
	GDALClose(dstDataset);
	GDALClose(warpedVRT);
	GDALClose(srcDataset);
}


void ResampleByReferenceRasterB(QString pszSrcFile, QString RefrasterBPath, QString pszOutFile, GDALResampleAlg eResample) {
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
	GDALDataset* pDSrc = (GDALDataset*)GDALOpen(pszSrcFile.toLocal8Bit().constData(), GA_ReadOnly);
	if (pDSrc == NULL) {
		qDebug() << u8"do not  open In Raster file: " << pszSrcFile;
		return  ;
	}

	GDALDataset* pDRef = (GDALDataset*)GDALOpen(RefrasterBPath.toLocal8Bit().constData(), GA_ReadOnly);
	if (pDRef == NULL) {
		qDebug() << u8"do not open Ref Raster file: " << RefrasterBPath;
		return;
	}

	QString filesuffer = getFileExtension(pszOutFile).toLower();
	bool isTiff = filesuffer.contains("tif");
	GDALDriver* pDriver = isTiff ? GetGDALDriverManager()->GetDriverByName("GTiff") : GetGDALDriverManager()->GetDriverByName("ENVI");

	if (pDriver == NULL) {
		qDebug() << "not open driver";
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		return ;
	}
	int			 width = pDSrc->GetRasterXSize();
	int			 height = pDSrc->GetRasterYSize();
	int			 nBandCount = pDSrc->GetRasterCount();
	GDALDataType dataType = pDSrc->GetRasterBand(1)->GetRasterDataType();

	char* pszSrcWKT = NULL;
	pszSrcWKT = const_cast<char*>(pDSrc->GetProjectionRef());

	// 锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7
	if (strlen(pszSrcWKT) <= 0) {
		OGRSpatialReference oSRS;
		oSRS.importFromEPSG(4326);
		// oSRS.SetUTM(50, true);   //锟斤拷锟斤拷锟斤拷  锟斤拷锟斤拷120锟斤拷
		// oSRS.SetWellKnownGeogCS("WGS84");
		oSRS.exportToWkt(&pszSrcWKT);
	}

	char* pdstSrcWKT = NULL;
	pdstSrcWKT = const_cast<char*>(pDRef->GetProjectionRef());

	// 锟斤拷锟矫伙拷锟酵队帮拷锟斤拷锟轿?拷锟斤拷锟揭伙拷锟?1锟?7
	if (strlen(pdstSrcWKT) <= 0)
	{
		OGRSpatialReference oSRS;
		oSRS.importFromEPSG(4326);
		// oSRS.SetUTM(50, true);   //锟斤拷锟斤拷锟斤拷  锟斤拷锟斤拷120锟斤拷
		// oSRS.SetWellKnownGeogCS("WGS84");
		oSRS.exportToWkt(&pdstSrcWKT);
	}


	int			 new_width = pDRef->GetRasterXSize();
	int			 new_height = pDRef->GetRasterYSize();
	double gt[6]  ;
	pDRef->GetGeoTransform(gt);

	// GDALDestroyGenImgProjTransformer(hTransformArg);
	qDebug() << "create init GDALDataset ";
	GDALDataset* pDDst =
		pDriver->Create(pszOutFile.toLocal8Bit().constData(), new_width, new_height, nBandCount, dataType, NULL);
	if (pDDst == NULL) {
		qDebug() << "not create init GDALDataset ";
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
		GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
		return  ;
	}

	pDDst->SetProjection(pdstSrcWKT);
	pDDst->SetGeoTransform(gt);


	qDebug() << "GDALCreateGenImgProjTransformer " << Qt::endl;
	void* hTransformArg;
	hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT,
		FALSE, 0.0, 1);
	qDebug() << "no proj ";
	//(没锟斤拷投影锟斤拷影锟斤拷锟斤拷锟斤拷锟斤拷卟锟酵?拷锟?1锟?7)
	if (hTransformArg == NULL) {
		qDebug() << "hTransformArg create failure";
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
		GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
		return;
	}
	qDebug() << "has proj ";
	double dGeoTrans[6] = { 0 };
	int	   nNewWidth = 0, nNewHeight = 0;
	if (GDALSuggestedWarpOutput((GDALDatasetH)pDSrc, GDALGenImgProjTransform, hTransformArg,
		dGeoTrans, &nNewWidth, &nNewHeight)
		!= CE_None) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
		GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
		return  ;
	}



	GDALWarpOptions* psWo = GDALCreateWarpOptions();
	CPLSetConfigOption("GDAL_NUM_THREADS", "ALL_CPUS"); // 使用所有可用的CPU核心
	CPLSetConfigOption("GDAL_CACHEMAX", "16000");         // 设置缓存大小为500MB
	// psWo->papszWarpOptions = CSLDuplicate(NULL);
	psWo->eWorkingDataType = dataType;



	psWo->eResampleAlg = eResample;

	psWo->hSrcDS = (GDALDatasetH)pDSrc;
	psWo->hDstDS = (GDALDatasetH)pDDst;
	qDebug() << "GDALCreateGenImgProjTransformer"  ;
	psWo->pfnTransformer = GDALGenImgProjTransform;
	psWo->pTransformerArg = GDALCreateGenImgProjTransformer(
		(GDALDatasetH)pDSrc, pszSrcWKT, (GDALDatasetH)pDDst, pszSrcWKT, FALSE, 0.0, 1);
	;

	qDebug() << "GDALCreateGenImgProjTransformer has created" << Qt::endl;
	psWo->nBandCount = nBandCount;
	psWo->panSrcBands = (int*)CPLMalloc(nBandCount * sizeof(int));
	psWo->panDstBands = (int*)CPLMalloc(nBandCount * sizeof(int));
	for (int i = 0; i < nBandCount; i++) {
		psWo->panSrcBands[i] = i + 1;
		psWo->panDstBands[i] = i + 1;
	}

	GDALWarpOperation oWo;
	if (oWo.Initialize(psWo) != CE_None) {
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
		GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
		return ;
	}
	qDebug() << "ChunkAndWarpImage:" << new_width << "," << new_height << Qt::endl;
	oWo.ChunkAndWarpMulti(0, 0, new_width, new_height);
	GDALFlushCache(pDDst);
	qDebug() << "ChunkAndWarpImage over" << Qt::endl;
	// GDALDestroyGenImgProjTransformer(psWo->pTransformerArg);
	// GDALDestroyWarpOptions(psWo);
	GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
	GDALClose((GDALDatasetH)(GDALDatasetH)pDDst);
	GDALClose((GDALDatasetH)(GDALDatasetH)pDRef);
	GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	return ;
}



void ResampleByReferenceRasterB(QString InrasterAPath, QString RefrasterBPath, QString OutrasterCPath) {
    // 注册所有GDAL驱动
    GDALAllRegister();

    // 打开参考栅格B
    GDALDataset* refDS = (GDALDataset*)GDALOpen(RefrasterBPath.toUtf8().constData(), GA_ReadOnly);
    if (!refDS) {
        qDebug() << "无法打开参考栅格B：" << RefrasterBPath;
        return;
    }

    // 获取参考栅格的地理变换、投影和尺寸
    double geotransform[6];
    if (refDS->GetGeoTransform(geotransform) != CE_None) {
        qDebug() << "获取参考栅格的地理变换失败。";
        GDALClose(refDS);
        return;
    }

    const char* proj = refDS->GetProjectionRef();
    int cols = refDS->GetRasterXSize();
    int rows = refDS->GetRasterYSize();
    GDALClose(refDS); // 获取信息后关闭参考栅格

    // 打开输入栅格A
    GDALDataset* srcDS = (GDALDataset*)GDALOpen(InrasterAPath.toUtf8().constData(), GA_ReadOnly);
    if (!srcDS) {
        qDebug() << "无法打开输入栅格A：" << InrasterAPath;
        return;
    }

    // 获取输入栅格的波段数和数据类型
    int nBands = srcDS->GetRasterCount();
    if (nBands == 0) {
        qDebug() << "输入栅格没有波段数据。";
        GDALClose(srcDS);
        return;
    }
    GDALDataType dataType = srcDS->GetRasterBand(1)->GetRasterDataType();

    // 创建输出栅格C
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
    if (!driver) {
        qDebug() << "无法获取GeoTIFF驱动。";
        GDALClose(srcDS);
        return;
    }

    GDALDataset* dstDS = driver->Create(
        OutrasterCPath.toUtf8().constData(),
        cols,
        rows,
        nBands,
        dataType,
        nullptr
    );
    if (!dstDS) {
        qDebug() << "无法创建输出栅格：" << OutrasterCPath;
        GDALClose(srcDS);
        return;
    }

    // 设置输出栅格的地理变换和投影
    dstDS->SetGeoTransform(geotransform);
    dstDS->SetProjection(proj);

    // 配置GDAL Warp选项
    GDALWarpOptions* psWO = GDALCreateWarpOptions();
    psWO->hSrcDS = srcDS;
    psWO->hDstDS = dstDS;
    psWO->nBandCount = nBands;
    psWO->panSrcBands = (int*)CPLMalloc(nBands * sizeof(int));
    psWO->panDstBands = (int*)CPLMalloc(nBands * sizeof(int));
    for (int i = 0; i < nBands; ++i) {
        psWO->panSrcBands[i] = i + 1;
        psWO->panDstBands[i] = i + 1;
    }
    psWO->eResampleAlg = GRA_NearestNeighbour; // 使用最近邻重采样

    // 初始化坐标转换器
    psWO->pfnTransformer = GDALGenImgProjTransform;
    psWO->pTransformerArg = GDALCreateGenImgProjTransformer(
        srcDS, GDALGetProjectionRef(srcDS),
        dstDS, GDALGetProjectionRef(dstDS),
        FALSE, 0.0, 1
    );
    if (!psWO->pTransformerArg) {
        qDebug() << "创建坐标转换器失败。";
        GDALDestroyWarpOptions(psWO);
        GDALClose(srcDS);
        GDALClose(dstDS);
        return;
    }

    // 执行Warp操作
    GDALWarpOperation oWarp;
    if (oWarp.Initialize(psWO) != CE_None) {
        qDebug() << "初始化Warp操作失败。";
        GDALDestroyGenImgProjTransformer(psWO->pTransformerArg);
        GDALDestroyWarpOptions(psWO);
        GDALClose(srcDS);
        GDALClose(dstDS);
        return;
    }

    CPLErr eErr = oWarp.ChunkAndWarpImage(0, 0, cols, rows);
    if (eErr != CE_None) {
        qDebug() << "执行Warp操作失败。";
    }

    // 清理资源
    GDALDestroyGenImgProjTransformer(psWO->pTransformerArg);
    GDALDestroyWarpOptions(psWO);
    GDALClose(srcDS);
    GDALClose(dstDS);

    qDebug() << "重采样完成，结果已保存至：" << OutrasterCPath;
}


void CreateSARIntensityByLookTable(QString IntensityRasterPath,
	QString LookTableRasterPath,
	QString SARIntensityPath,
	long min_rid, long max_rid,
	long min_cid, long max_cid,
	std::function<void(long, long)> processBarShow 
)
{
	gdalImage looktableds(LookTableRasterPath);
	gdalImage geoIntensity(IntensityRasterPath);
	gdalImage SARIntensity= CreategdalImageDouble(SARIntensityPath,max_rid-min_rid,max_cid-min_cid,1);

	long blockYSize = Memory1GB / looktableds.width / 8 * 2;
	
	Eigen::MatrixXd SARData = SARIntensity.getData(0, 0, SARIntensity.height, SARIntensity.width, 1);
	SARData = SARData.array() * 0;
	// 分块处理
	for (int yOff = 0; yOff < looktableds.height; yOff += blockYSize)
	{
		processBarShow(yOff, looktableds.height);
		qDebug() << "Process : [" << yOff * 100.0 / looktableds.height << " % ]";
		Eigen::MatrixXd rowData = looktableds.getData(yOff, 0, blockYSize, looktableds.width, 1);
		Eigen::MatrixXd colData = looktableds.getData(yOff, 0, blockYSize, looktableds.width, 2);
		Eigen::MatrixXd geoData = geoIntensity.getData(yOff, 0, blockYSize, looktableds.width, 1);

		for (long i = 0; i < rowData.rows(); i++) {
			for (long j = 0; j < rowData.cols(); j++) {
				long r =round( rowData(i,j))-min_rid;
				long c = round(colData(i, j))-min_cid;

				if (r >= 0 && r < SARIntensity.height && c >= 0 && c < SARIntensity.width) {
					SARData(r, c) = SARData(r, c) + geoData(i, j);
				}
			}
		}
	}
	SARIntensity.saveImage(SARData, 0, 0, 1);
	qDebug() << "Process : [ 100 % ]";
	processBarShow(1000,1000);
}


bool saveEigenMatrixXd2Bin(Eigen::MatrixXd data, QString dataStrPath)
{

	Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);

	gdalImage img = CreategdalImageDouble(dataStrPath, data.rows(), data.cols(), 1, gt, "", false, true, true);

	img.saveImage(data, 0, 0, 1);

	return true;
}