#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 操作

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


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 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;
}




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 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;
}




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;
}



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

	// 打开输入栅格文件
	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");
	 
	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 cropRasterByLatLon(const char* inputFile, const char* outputFile, double minLon, double maxLon, double minLat, double maxLat) {
	// 初始化 GDAL 库
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");

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

	// 获取栅格数据的地理参考信息
	double adfGeoTransform[6];
	if (poDataset->GetGeoTransform(adfGeoTransform) != CE_None) {
		qDebug() << "Failed to get geotransform." ;
		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) {
		qDebug() << "Failed to get GTiff driver." ;
		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) {
		qDebug() << "Failed to create output raster." ;
		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((GDALDatasetH)(GDALDatasetH)poDataset);
	GDALClose((GDALDatasetH)(GDALDatasetH)poOutDataset);
	GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
}

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 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.";
}


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

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

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

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

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

	// 设置输出数据集的地理变换（坐标系）
	double targetGeoTransform[6] = {
		adfGeoTransform[0], targetPixelSizeX, 0, adfGeoTransform[3], 0, -targetPixelSizeY
	};
	GDALClose((GDALDatasetH)(GDALDatasetH)poDataset);
	GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	ResampleGDAL(inputRaster, outputRaster, targetGeoTransform, targetXSize, targetYSize,GRA_Bilinear);
	qDebug() << "Resampling completed.";
}

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;
	}
}



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;
		}
	}
}



long getProjectEPSGCodeByLon_Lat(double lon, double lat, ProjectStripDelta stripState)
{
	long EPSGCode = 0;
	switch (stripState) {
	case ProjectStripDelta::Strip_3: {
		EPSGCode = getProjectEPSGCodeByLon_Lat_inStrip3(lon, lat);
		break;
	};
	case ProjectStripDelta::Strip_6: {
		EPSGCode = getProjectEPSGCodeByLon_Lat_inStrip6(lon, lat);
		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;
}



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;
	}
}



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", "4000");         // 设置缓存大小为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;


}



/** 2025.03.25 下面的函数存在 Eigen使用****************************/



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;
}


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


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";

}



ErrorCode  DEM2XYZRasterAndSlopRaster(QString dempath, QString demxyzpath, QString demsloperPath)
{

	omp_lock_t lock;
	omp_init_lock(&lock);


	QString DEMPath = dempath;
	QString XYZPath = demxyzpath;
	QString SLOPERPath = demsloperPath;
	gdalImage demds(DEMPath);
	gdalImage demxyz = CreategdalImageDouble(XYZPath, demds.height, demds.width, 3, demds.gt, demds.projection, true, true, false);// X,Y,Z

	int64_t linecount = 0;
	qDebug() << u8"start dem (lon,lat,ati) -> dem [X,Y,Z]";
	int64_t line_invert= Memory1MB / 8.0 / demds.width * 200;

#pragma omp parallel for
	for (int64_t max_rows_ids = 0; max_rows_ids < demds.height; max_rows_ids = max_rows_ids + line_invert) {

		long temp_line_invert = max_rows_ids + line_invert < demds.height ? line_invert : demds.height - max_rows_ids;
		double R = 0;
		double dem_row = 0, dem_col = 0, dem_alt = 0;
		Eigen::MatrixXd demdata = demds.getData(max_rows_ids, 0, temp_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++) {
			Landpoint LandP{ 0,0,0 };
			Point3 GERpoint{ 0,0,0 };
			double rowidx = 0;
			double colidx = 0;
			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;
			}
		}
		omp_set_lock(&lock);
		
		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);
		linecount = linecount + temp_line_invert;
		qDebug() << "dem -> XYZ [" << linecount * 100.0 / demds.height << "] %";
		omp_unset_lock(&lock);	 // 锟酵放伙拷斤拷
	}

	// 计算坡向角
	qDebug() << u8"start dem (lon,lat,ati) -> dem Sloper[X,Y,Z]";
	gdalImage demsloperxyz = CreategdalImageDouble(SLOPERPath, demds.height, demds.width, 4, demds.gt, demds.projection, true, true, false);// X,Y,Z,cosangle
	
	linecount = 0;
	line_invert = Memory1MB / 8.0 / demds.width * 200;;
#pragma omp parallel for
	for (int64_t start_ids = 1; start_ids < demds.height; start_ids = start_ids + line_invert) {

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

		dem_rows = demsloper_x.rows();
		dem_cols = demsloper_x.cols();


		for (long i = 0; i < dem_rows ; i++) {
			Landpoint p0, p1, p2, p3, p4, pslopeVector, pp;
			Vector3D  slopeVector;
			double sloperAngle = 0;
			Vector3D Zaxis = { 0,0,1 };

			double rowidx = 0, colidx = 0;
			for (long j = 1; j < dem_cols - 1; j++) {
				rowidx = i + 1;
				colidx = j;
				demds.getLandPoint(rowidx, colidx, demdata(i + 1, j + 1), p0);
				demds.getLandPoint(rowidx - 1, colidx, demdata(i - 1 + 1, j), p1);
				demds.getLandPoint(rowidx, colidx - 1, demdata(i + 1, j - 1), p2);
				demds.getLandPoint(rowidx + 1, colidx, demdata(i + 1 + 1, j), p3);
				demds.getLandPoint(rowidx, colidx + 1, demdata(i + 1, 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;

			}
		}
		omp_set_lock(&lock);
		demsloperxyz.saveImage(demsloper_x, start_ids , 0, 1);
		demsloperxyz.saveImage(demsloper_y, start_ids , 0, 2);
		demsloperxyz.saveImage(demsloper_z, start_ids , 0, 3);
		demsloperxyz.saveImage(demsloper_angle, start_ids , 0, 4);
		linecount = linecount + temp_line_invert;
		qDebug() << "dem -> Sloper [" << linecount * 100.0 / demds.height << "] %";
		omp_unset_lock(&lock);	 // 锟酵放伙拷斤拷
	
	}


	omp_destroy_lock(&lock); // 劫伙拷斤拷
	return ErrorCode::SUCCESS;
}



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);
}








/**​
* @brief 将VRT文件转换为ENVI格式（生成.dat和.hdr文件）
* @param inputVrtPath 输入VRT文件路径
* @param outputEnviPath 输出ENVI数据文件路径（不含扩展名，自动生成.dat和.hdr）
* @param papszOptions GDAL创建选项（如压缩参数、存储顺序等，默认NULL）
* @return 成功返回true，失败返回false
*/
bool ConvertVrtToEnvi(QString vrtPath, QString outPath) {

	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES"); // 注绞斤拷斤拷锟?1锟?7
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");

	const char* inputVrtPath = vrtPath.toUtf8().constData();
	const char* outputEnviPath = outPath.toUtf8().constData();

	// 注册所有GDAL驱动[7](@ref)
	GDALAllRegister();

	// 打开输入VRT文件[1,3](@ref)
	GDALDataset* poVRTDataset = (GDALDataset*)GDALOpen(inputVrtPath, GA_ReadOnly);
	if (!poVRTDataset) {
		qDebug() << "错误：无法打开VRT文件 " << inputVrtPath ;
		return false;
	}

	// 获取ENVI驱动[4,7](@ref)
	GDALDriver* poENVIDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
	if (!poENVIDriver) {
		qDebug() << "错误：ENVI驱动未加载" ;
		GDALClose(poVRTDataset);
		return false;
	}

	// 构造输出文件名（强制添加.dat后缀）[4](@ref)
	std::string outputPath = std::string(outputEnviPath) + ".dat";

	// 执行格式转换（自动生成.hdr头文件）[4,7](@ref)
	GDALDataset* poENVIDataset = poENVIDriver->CreateCopy(
		outputPath.c_str(),    // 输出路径
		poVRTDataset,           // 输入数据集
		FALSE,                  // 非严格模式（允许格式调整）
		nullptr,           // 用户自定义选项（如压缩参数）
		nullptr, nullptr        // 进度条和参数（暂不启用）
	);

	// 检查转换结果
	bool success = (poENVIDataset != nullptr);
	if (!success) {
		qDebug() << "错误：ENVI文件创建失败" ;
	}
	else {
		GDALClose(poENVIDataset);
	}

	// 关闭输入数据集
	GDALClose(poVRTDataset);
	return success;
}




// 遥感类
ErrorCode ResampleDEM(QString indemPath, QString outdemPath, double gridx, double gridy)
{
	double gridlat = gridy / 110000.0;
	double gridlon = gridx / 100000.0;

	long espgcode = GetEPSGFromRasterFile(indemPath.toUtf8().constData());
	if (espgcode == 4326) {
		resampleRaster(indemPath.toUtf8().constData(), outdemPath.toUtf8().constData(), gridlon, gridlat);
		return ErrorCode::SUCCESS;
	}
	else {
		QMessageBox::warning(nullptr, u8"警告", u8"请输入WGS84坐标的DEM影像");
		return ErrorCode::FAIL;
	}
}

void BASECONSTVARIABLEAPI CloseAllGDALRaster()
{
	GDALDestroyDriverManager();
	return ;
}



ErrorCode Complex2AmpRaster(QString inComplexPath, QString outRasterPath)
{
	gdalImageComplex inimg(inComplexPath);
	gdalImage ampimg = CreategdalImage(outRasterPath, inimg.height, inimg.width, inimg.band_num, inimg.gt, inimg.projection, true, true);

	long blocklines = Memory1GB / 8 / inimg.width;
	blocklines = blocklines < 100 ? 100 : blocklines;
	omp_lock_t lock;
	omp_init_lock(&lock);
#pragma omp parallel for
 	for (int64_t startrow = 0; startrow < inimg.height; startrow = startrow + blocklines) {

		Eigen::MatrixXd imgArrb1 = ampimg.getData(startrow, 0, blocklines, inimg.width, 1);
		Eigen::MatrixXcd imgArr = inimg.getDataComplex(startrow, 0, blocklines, inimg.width, 1);
		imgArrb1 = imgArr.array().abs();
		omp_set_lock(&lock);
		ampimg.saveImage(imgArrb1, startrow, 0, 1);
		omp_unset_lock(&lock);	 //  
	}

	omp_destroy_lock(&lock); //  
	qDebug() << u8"影像写入到：" << outRasterPath;
	return ErrorCode::SUCCESS;
}

ErrorCode Complex2PhaseRaster(QString inComplexPath, QString outRasterPath)
{
	gdalImageComplex inimg(inComplexPath);
	gdalImage ampimg = CreategdalImage(outRasterPath, inimg.height, inimg.width, inimg.band_num, inimg.gt, inimg.projection, true, true);


	long blocklines = Memory1MB  / 8 / inimg.width*200;
	blocklines = blocklines < 100 ? 100 : blocklines;
 ;
	omp_lock_t lock;
	omp_init_lock(&lock);
#pragma omp parallel for
 	for (int64_t startrow = 0; startrow < inimg.height; startrow = startrow + blocklines) {

		Eigen::MatrixXd imgArrb1 = ampimg.getData(startrow, 0, blocklines, inimg.width, 1);
		Eigen::MatrixXcd  imgArr = inimg.getData(startrow, 0, blocklines, inimg.width, 1);
		imgArrb1 = imgArr.array().arg();
		omp_set_lock(&lock);
		ampimg.saveImage(imgArrb1, startrow, 0, 1);
		omp_unset_lock(&lock);	 // 
	}
	omp_destroy_lock(&lock); //  
	qDebug() << "影像写入到：" << outRasterPath;
	return ErrorCode::SUCCESS;
}

ErrorCode Complex2dBRaster(QString inComplexPath, QString outRasterPath)
{
	gdalImageComplex inimg(inComplexPath);
	gdalImage ampimg = CreategdalImage(outRasterPath, inimg.height, inimg.width, inimg.band_num, inimg.gt, inimg.projection, true, true);


	long blocklines = Memory1MB / 8 / inimg.width * 200;
	blocklines = blocklines < 100 ? 100 : blocklines;
 
	omp_lock_t lock;
	omp_init_lock(&lock);
#pragma omp parallel for
	for (int64_t startrow = 0; startrow < inimg.height; startrow = startrow + blocklines) {
		Eigen::MatrixXd imgArrb1 = ampimg.getData(startrow, 0, blocklines, inimg.width, 1);
		Eigen::MatrixXcd imgArr = inimg.getData(startrow, 0, blocklines, inimg.width, 1);
		imgArrb1 = imgArr.array().abs().log10() * 20.0;
		omp_set_lock(&lock);
		ampimg.saveImage(imgArrb1, startrow, 0, 1);
		omp_unset_lock(&lock);	 // 
	}
	omp_destroy_lock(&lock); // 
	qDebug() << "影像写入到：" << outRasterPath;
	return ErrorCode::SUCCESS;
}



ErrorCode BASECONSTVARIABLEAPI amp2dBRaster(QString inPath, QString outRasterPath)
{
	gdalImage inimg(inPath);
	gdalImage dBimg = CreategdalImage(outRasterPath, inimg.height, inimg.width, inimg.band_num, inimg.gt, inimg.projection, true, true);

	long blocklines = Memory1MB / 8 / inimg.width * 200;
	blocklines = blocklines < 100 ? 100 : blocklines;

	omp_lock_t lock;
	omp_init_lock(&lock);
#pragma omp parallel for
	for (int64_t startrow = 0; startrow < inimg.height; startrow = startrow + blocklines) {

		Eigen::MatrixXd imgArrdB = dBimg.getData(startrow, 0, blocklines, inimg.width, 1);
		Eigen::MatrixXd imgArr = inimg.getData(startrow, 0, blocklines, inimg.width, 1);
		imgArrdB = imgArr.array().log10() * 20.0;
		omp_set_lock(&lock);
		dBimg.saveImage(imgArrdB, startrow, 0, 1);
		omp_unset_lock(&lock);	 // 
	}
	omp_destroy_lock(&lock); // 
	qDebug() << "影像写入到：" << outRasterPath;
	return ErrorCode::SUCCESS;
}



void MultiLookRaster(QString inRasterPath, QString outRasterPath, long looklineNumrow, long looklineNumCol)
{
	gdalImage inRaster(inRasterPath);

	int outRows = inRaster.height / looklineNumrow;
	int outCols = inRaster.width / looklineNumCol;
	int bandnum = 1;
	QString project = inRaster.projection;
	Eigen::MatrixXd inRasterGt = inRaster.gt;
	GDALDataType datetype = inRaster.getDataType();

	gdalImage outRaster = CreategdalImage(outRasterPath, outRows, outCols, 1, inRasterGt, project, true, true, true, datetype);

	Eigen::MatrixXd inRasterArr = inRaster.getData(0, 0, inRaster.height, inRaster.width, 1);
	Eigen::MatrixXd outRasterArr = outRaster.getData(0, 0, outRows, outCols, 1);


	// 多视处理
#pragma omp parallel for collapse(2)
	for (int row = 0; row < outRows; ++row) {
		for (int col = 0; col < outCols; ++col) {
			// 计算输入矩阵的窗口位置[2,3](@ref)
			const int inputRow = row * looklineNumrow;
			const int inputCol = col * looklineNumCol;

			// 动态计算实际窗口大小（处理边界情况）[3](@ref)
			const int actualRows = (row == outRows - 1) ?
				(inRaster.height - inputRow) : looklineNumrow;
			const int actualCols = (col == outCols - 1) ?
				(inRaster.width - inputCol) : looklineNumCol;

			// 提取当前窗口数据块[7](@ref)
			Eigen::MatrixXd window = inRasterArr.block(inputRow, inputCol, actualRows, actualCols);

			// 计算多视平均值（自动处理数据类型转换）[2,7](@ref)
			double sum = 0.0;
			if constexpr (std::is_integral_v<typename Eigen::MatrixXd::Scalar>) {
				sum = window.cast<double>().sum(); // 整型数据转换为双精度计算
			}
			else {
				sum = window.sum();
			}
			outRasterArr(row, col) = sum / (actualRows * actualCols);
		}
	}

	// 可选：处理残留边界（当输入尺寸不是整数倍时）
	if ((inRaster.height % looklineNumrow != 0) || (inRaster.width % looklineNumCol != 0)) {
		// 底部边界处理[3](@ref)
		const int residualRowStart = outRows * looklineNumrow;
		const int residualRowSize = inRaster.height - residualRowStart;
		if (residualRowSize > 0) {
#pragma omp parallel for
			for (int col = 0; col < outCols; ++col) {
				const int inputCol = col * looklineNumCol;
				const int actualCols = (col == outCols - 1) ?
					(inRaster.width - inputCol) : looklineNumCol;

				Eigen::MatrixXd window = inRasterArr.block(residualRowStart, inputCol,
					residualRowSize, actualCols);
				outRasterArr(outRows - 1, col) = window.mean();
			}
		}

		// 右侧边界处理[3](@ref)
		const int residualColStart = outCols * looklineNumCol;
		const int residualColSize = inRaster.width - residualColStart;
		if (residualColSize > 0) {
#pragma omp parallel for
			for (int row = 0; row < outRows; ++row) {
				const int inputRow = row * looklineNumrow;
				const int actualRows = (row == outRows - 1) ?
					(inRaster.height - inputRow) : looklineNumrow;

				Eigen::MatrixXd window = inRasterArr.block(inputRow, residualColStart,
					actualRows, residualColSize);
				outRasterArr(row, outCols - 1) = window.mean();
			}
		}
	}

	// 保存结果
	outRaster.saveImage(outRasterArr, 0, 0, 1);

}













/* 启动下面的函数，会导致编译器错误，可能与Eigen 库 本身的bug相关，不建议排查，太费时间，而且不一定能排查出来
*
***

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;
}


****/