LAMPCAE/src/LAMPTool/BaseToollib/ImageOperatorBase.cpp

#include "ImageOperatorBase.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.h>
#include <memory>
#include <iostream>
#include "FileOperator.h"

using namespace std;
using namespace Eigen;

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






std::shared_ptr<GDALDataset> OpenDataset(const QString& in_path,GDALAccess rwmode)
{
	GDALAllRegister();
	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();
	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();
	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<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<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<complex<float>> 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<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();

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


gdalImage::gdalImage()
{
	this->height = 0;
	this->width = 0;
	this->data_band_ids = 1;
	this->start_row = 0;
	this->start_col = 0;
}

/// <summary>
/// 斤拷图饺∮帮拷锟?1锟?7
/// </summary>
/// <param name="dem_path"></param>
gdalImage::gdalImage(const QString& raster_path)
{
	omp_lock_t lock;
	omp_init_lock(&lock); // 锟斤拷始斤拷斤拷
	omp_set_lock(&lock); //锟斤拷没斤拷锟?1锟?7
	this->img_path = raster_path;

	GDALAllRegister();// 注绞斤拷斤拷锟?1锟?7
	// 锟斤拷DEM影锟斤拷
	GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(raster_path.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷
	this->width = rasterDataset->GetRasterXSize();
	this->height = rasterDataset->GetRasterYSize();
	this->band_num = rasterDataset->GetRasterCount();

	double* gt = new double[6];
	// 锟斤拷梅斤拷锟斤拷
	rasterDataset->GetGeoTransform(gt);
	this->gt = Eigen::MatrixXd(2, 3);
	this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5];

	this->projection = rasterDataset->GetProjectionRef();
	// 斤拷斤拷
	//double* inv_gt = new double[6];;
	//GDALInvGeoTransform(gt, inv_gt); // 斤拷斤拷
	// 斤拷投影
	GDALFlushCache((GDALDatasetH)rasterDataset);
	GDALClose((GDALDatasetH)rasterDataset);
	rasterDataset = NULL;// 指矫匡拷
	this->InitInv_gt();
	delete[] gt;
	////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	omp_unset_lock(&lock); //锟酵放伙拷斤拷
	omp_destroy_lock(&lock); //劫伙拷斤拷
}

gdalImage::~gdalImage()
{
}

void gdalImage::setHeight(int height)
{
	this->height = height;
}

void gdalImage::setWidth(int width)
{
	this->width = width;
}

void gdalImage::setTranslationMatrix(Eigen::MatrixXd gt)
{
	this->gt = gt;
}

void gdalImage::setData(Eigen::MatrixXd, int data_band_ids)
{
	this->data = data;
	this->data_band_ids = data_band_ids;
}

 

Eigen::MatrixXd gdalImage::getData(int start_row, int start_col, int rows_count, int cols_count, int band_ids = 1)
{
	omp_lock_t lock;
	omp_init_lock(&lock);
	omp_set_lock(&lock);
	GDALAllRegister();

	GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷

	GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
	GDALRasterBand* demBand = rasterDataset->GetRasterBand(band_ids);

	rows_count = start_row + rows_count <= this->height ? rows_count : this->height - start_row;
	cols_count = start_col + cols_count <= this->width ? cols_count : this->width - start_col;

	MatrixXd datamatrix(rows_count, cols_count);


	if (gdal_datatype == GDT_Byte) {
		char* temp = new char[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);
		for (int i = 0; i < rows_count; i++) {
			for (int j = 0; j < cols_count; j++) {
				datamatrix(i, j) = temp[i * cols_count + j];
			}
		}
		delete[] temp;
	}
	else if (gdal_datatype == 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);
		for (int i = 0; i < rows_count; i++) {
			for (int j = 0; j < cols_count; j++) {
				datamatrix(i, j) = temp[i * cols_count + j];
			}
		}
		delete[] temp;
	}
	else if (gdal_datatype == 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);
		for (int i = 0; i < rows_count; i++) {
			for (int j = 0; j < cols_count; j++) {
				datamatrix(i, j) = temp[i * cols_count + j];
			}
		}
		delete[] temp;
	}
	else if (gdal_datatype == GDT_UInt32) {
		unsigned int* temp = new unsigned int[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);
		for (int i = 0; i < rows_count; i++) {
			for (int j = 0; j < cols_count; j++) {
				datamatrix(i, j) = temp[i * cols_count + j];
			}
		}
		delete[] temp;
	}
	else if (gdal_datatype == GDT_Int32) {
		int* temp = new int[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);
		for (int i = 0; i < rows_count; i++) {
			for (int j = 0; j < cols_count; j++) {
				datamatrix(i, j) = temp[i * cols_count + j];
			}
		}
		delete[] temp;
	}
	//else if (gdal_datatype == GDT_UInt64) {
	//	unsigned long* temp = new unsigned long[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);
	//	for (int i = 0; i < rows_count; i++) {
	//		for (int j = 0; j < cols_count; j++) {
	//			datamatrix(i, j) = temp[i * cols_count + j];
	//		}
	//	}
	//	delete[] temp;
	//}
	//else if (gdal_datatype == GDT_Int64) {
	//	long* temp = new long[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);
	//	for (int i = 0; i < rows_count; i++) {
	//		for (int j = 0; j < cols_count; j++) {
	//			datamatrix(i, j) = temp[i * cols_count + j];
	//		}
	//	}
	//	delete[] temp;
	//}
	else if (gdal_datatype == 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);

		for (int i = 0; i < rows_count; i++) {
			for (int j = 0; j < cols_count; j++) {
				datamatrix(i, j) = temp[i * cols_count + j];
			}
		}
		delete[] temp;
	}
	else if (gdal_datatype == GDT_Float64) {
		double* temp = new double[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);

		for (int i = 0; i < rows_count; i++) {
			for (int j = 0; j < cols_count; j++) {
				datamatrix(i, j) = temp[i * cols_count + j];
			}
		}
		delete[] temp;
	}
	else {}
	GDALClose((GDALDatasetH)rasterDataset);
	omp_unset_lock(&lock); //锟酵放伙拷斤拷
	omp_destroy_lock(&lock); //劫伙拷斤拷
	//GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	return datamatrix;

}

Eigen::MatrixXd gdalImage::getGeoTranslation()
{
	return this->gt;
}

GDALDataType gdalImage::getDataType()
{
	GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly));
	GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
	return gdal_datatype;
}

/// <summary>
/// 写锟斤拷遥锟斤拷影锟斤拷
/// </summary>
/// <param name="data"></param>
/// <param name="start_row"></param>
/// <param name="start_col"></param>
/// <param name="band_ids"></param>
void gdalImage::saveImage(Eigen::MatrixXd data, int start_row = 0, int start_col = 0, int band_ids = 1)
{
	omp_lock_t lock;
	omp_init_lock(&lock);
	omp_set_lock(&lock);
	if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
		QString tip = u8"file path: " + this->img_path;
		qDebug() << tip;
		throw exception(tip.toUtf8().constData());
	}
	GDALAllRegister();
	GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset* poDstDS = nullptr;
	if (exists_test(this->img_path)) {
		poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
	}
	else {
		poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height, this->band_num, GDT_Float32, NULL); // 斤拷锟斤拷
		poDstDS->SetProjection(this->projection.toUtf8().constData());


		double gt_ptr[6];
		for (int i = 0; i < this->gt.rows(); i++) {
			for (int j = 0; j < this->gt.cols(); j++) {
				gt_ptr[i * 3 + j] = this->gt(i, j);
			}
		}
		poDstDS->SetGeoTransform(gt_ptr);
		delete[] gt_ptr;
	}

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

	float* databuffer = new float[datarows * datacols];// (float*)malloc(datarows * datacols * sizeof(float));

	for (int i = 0; i < datarows; i++) {
		for (int j = 0; j < datacols; j++) {
			float temp = float(data(i, j));
			databuffer[i * datacols + j] = temp;
		}
	}
	//poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols, datarows, GDT_Float32,band_ids, num,0,0,0);
	poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows, databuffer, datacols, datarows, GDT_Float32, 0, 0);

	GDALFlushCache(poDstDS);
	GDALClose((GDALDatasetH)poDstDS);
	//GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	delete[] databuffer;
	omp_unset_lock(&lock); //锟酵放伙拷斤拷
	omp_destroy_lock(&lock); //劫伙拷斤拷
}

void gdalImage::saveImage()
{
	this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids);
}

void gdalImage::setNoDataValue(double nodatavalue = -9999, int band_ids = 1)
{
	GDALAllRegister();// 注绞斤拷斤拷锟?1锟?7
	//GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset* poDstDS = (GDALDataset*)(GDALOpen(img_path.toUtf8().constData(), GA_Update));
	poDstDS->GetRasterBand(band_ids)->SetNoDataValue(nodatavalue);
	GDALFlushCache((GDALDatasetH)poDstDS);
	GDALClose((GDALDatasetH)poDstDS);
}

int gdalImage::InitInv_gt()
{
	//1 lon lat = x
	//1 lon lat = y
	Eigen::MatrixXd temp=Eigen::MatrixXd::Zero(2, 3);
	this->inv_gt = temp;
	double a = this->gt(0, 0);
	double b = this->gt(0, 1);
	double c = this->gt(0, 2);
	double d = this->gt(1, 0);
	double e = this->gt(1, 1);
	double f = this->gt(1, 2);
	double g = 1;
	double det_gt = b * f - c * e;
	if (det_gt == 0) {
		return 0;
	}
	this->inv_gt(0, 0) = (c * d - a * f) / det_gt; //2
	this->inv_gt(0, 1) = f / det_gt; // lon
	this->inv_gt(0, 2) = -c / det_gt; // lat
	this->inv_gt(1, 0) = (a * e - b * d) / det_gt; //1
	this->inv_gt(1, 1) = -e / det_gt; // lon
	this->inv_gt(1, 2) = b / det_gt; // lat
	return 1;
}

Landpoint gdalImage::getRow_Col(double lon, double lat)
{
	Landpoint p{ 0,0,0 };
	p.lon = this->inv_gt(0, 0) + lon * this->inv_gt(0, 1) + lat * this->inv_gt(0, 2); //x
	p.lat = this->inv_gt(1, 0) + lon * this->inv_gt(1, 1) + lat * this->inv_gt(1, 2); //y
	return p;
}

Landpoint gdalImage::getLandPoint(double row, double col, double ati = 0)
{
	Landpoint p{ 0,0,0 };
	p.lon = this->gt(0, 0) + col * this->gt(0, 1) + row * this->gt(0, 2); //x
	p.lat = this->gt(1, 0) + col * this->gt(1, 1) + row * this->gt(1, 2); //y
	p.ati = ati;
	return p;
}

double gdalImage::mean(int bandids)
{
	double mean_value = 0;
	double count = this->height * this->width;
	int line_invert = 100;
	int start_ids = 0;
	do {
		Eigen::MatrixXd	sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
		mean_value = mean_value + sar_a.sum() / count;
		start_ids = start_ids + line_invert;
	} while (start_ids < this->height);
	return mean_value;
}

double gdalImage::max(int bandids)
{
	double max_value = 0;
	bool state_max = true;
	int line_invert = 100;
	int start_ids = 0;
	double temp_max = 0;
	do {
		Eigen::MatrixXd	sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
		if (state_max) {
			state_max = false;
			max_value = sar_a.maxCoeff();
		}
		else {
			temp_max = sar_a.maxCoeff();
			if (max_value < temp_max) {
				max_value = temp_max;
			}
		}
		start_ids = start_ids + line_invert;
	} while (start_ids < this->height);
	return max_value;
}

double gdalImage::min(int bandids)
{
	double min_value = 0;
	bool state_min = true;
	int line_invert = 100;
	int start_ids = 0;
	double temp_min = 0;
	do {
		Eigen::MatrixXd	sar_a = this->getData(start_ids, 0, line_invert, this->width, bandids);
		if (state_min) {
			state_min = false;
			min_value = sar_a.minCoeff();
		}
		else {
			temp_min = sar_a.minCoeff();
			if (min_value < temp_min) {
				min_value = temp_min;
			}
		}
		start_ids = start_ids + line_invert;
	} while (start_ids < this->height);
	return min_value;
}

GDALRPCInfo gdalImage::getRPC()
{
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "NO");
	CPLSetConfigOption("GDAL_DATA", "./data");
	GDALAllRegister();//注斤拷锟斤拷
	//斤拷锟斤拷
	GDALDataset* pDS = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
	//锟斤拷元斤拷锟叫伙拷取RPC锟斤拷息
	char** papszRPC = pDS->GetMetadata("RPC");

	//斤拷取锟斤拷RPC锟斤拷息斤拷山峁癸拷锟?1锟?7
	GDALRPCInfo oInfo;
	GDALExtractRPCInfo(papszRPC, &oInfo);

	GDALClose((GDALDatasetH)pDS);

	return oInfo;
}

Eigen::MatrixXd gdalImage::getLandPoint(Eigen::MatrixXd points)
{
	if (points.cols() != 3) {
		throw new exception("the size of points is equit 3!!!");
	}

	Eigen::MatrixXd result(points.rows(), 3);
	result.col(2) = points.col(2);// 锟竭筹拷
	points.col(2) = points.col(2).array() * 0 + 1;
	points.col(0).swap(points.col(2));// 斤拷
	Eigen::MatrixXd gts(3, 2);
	gts.col(0) = this->gt.row(0);
	gts.col(1) = this->gt.row(1);

	result.block(0, 0, points.rows(), 2) = points * gts;
	return result;
}

Eigen::MatrixXd gdalImage::getHist(int bandids)
{
	GDALAllRegister();// 注绞斤拷斤拷锟?1锟?7
	// 锟斤拷DEM影锟斤拷
	GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷

	GDALDataType gdal_datatype = rasterDataset->GetRasterBand(1)->GetRasterDataType();
	GDALRasterBand* xBand = rasterDataset->GetRasterBand(bandids);


	double dfMin = this->min(bandids);
	double dfMax = this->max(bandids);
	int count = int((dfMax - dfMin) / 0.01);
	count = count > 255 ? count : 255;
	GUIntBig* panHistogram = new GUIntBig[count];
	xBand->GetHistogram(dfMin, dfMax, count, panHistogram, TRUE, FALSE, NULL, NULL);
	Eigen::MatrixXd result(count, 2);
	double delta = (dfMax - dfMin) / count;
	for (int i = 0; i < count; i++) {
		result(i, 0) = dfMin + i * delta;
		result(i, 1) = double(panHistogram[i]);
	}
	delete[] panHistogram;
	GDALClose((GDALDatasetH)rasterDataset);
	return result;
}


gdalImage CreategdalImage(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt,bool overwrite) {

	if (exists_test(img_path.toUtf8().constData())) {
		if (overwrite) {
			gdalImage result_img(img_path);
			return result_img;
		}
		else {
		throw "file has exist!!!";
		exit(1);
		}
	}
	GDALAllRegister();// 注锟斤拷锟绞斤拷锟斤拷锟斤拷锟?1锟?7
	GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_Float32, NULL); // 锟斤拷锟斤拷锟斤拷
	if (need_gt) {
		poDstDS->SetProjection(projection.toUtf8().constData());

		// 锟斤拷锟斤拷转锟斤拷锟斤拷锟斤拷
		double gt_ptr[6] = { 0 };
		for (int i = 0; i < gt.rows(); i++) {
			for (int j = 0; j < gt.cols(); j++) {
				gt_ptr[i * 3 + j] = gt(i, j);
			}
		}
		poDstDS->SetGeoTransform(gt_ptr);
	}
	for (int i = 1; i <= band_num; i++) {
		poDstDS->GetRasterBand(i)->SetNoDataValue(-9999);
	}
	GDALFlushCache((GDALDatasetH)poDstDS);
	GDALClose((GDALDatasetH)poDstDS);
	////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	gdalImage result_img(img_path);
	return result_img;
}

gdalImageComplex CreategdalImageComplex(const QString& img_path, int height, int width, int band_num, Eigen::MatrixXd gt, QString projection, bool need_gt, bool overwrite)
{
	if (exists_test(img_path.toUtf8().constData())) {
		if (overwrite) {
			gdalImageComplex result_img(img_path);
			return result_img;
		}
		else {
			throw "file has exist!!!";
			exit(1);
		}
	}
	GDALAllRegister();
	GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
	GDALDataset* poDstDS = poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_CFloat64, NULL);
	if (need_gt) {
		poDstDS->SetProjection(projection.toUtf8().constData());

		// 锟斤拷锟斤拷转锟斤拷锟斤拷锟斤拷
		double gt_ptr[6] = { 0 };
		for (int i = 0; i < gt.rows(); i++) {
			for (int j = 0; j < gt.cols(); j++) {
				gt_ptr[i * 3 + j] = gt(i, j);
			}
		}
		poDstDS->SetGeoTransform(gt_ptr);
	}
	for (int i = 1; i <= band_num; i++) {
		poDstDS->GetRasterBand(i)->SetNoDataValue(-9999);
	}
	GDALFlushCache((GDALDatasetH)poDstDS);
	GDALClose((GDALDatasetH)poDstDS);
	////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	gdalImageComplex result_img(img_path);
	return result_img;
}


int ResampleGDAL(const char* pszSrcFile, const char* pszOutFile, double* gt, int new_width, int new_height, GDALResampleAlg eResample)
{
	GDALAllRegister();
	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 " << endl;
	void* hTransformArg;
	hTransformArg = GDALCreateGenImgProjTransformer((GDALDatasetH)pDSrc, pszSrcWKT, NULL, pszSrcWKT, FALSE, 0.0, 1);
	qDebug() << "no proj " << endl;
	//(没锟斤拷投影锟斤拷影锟斤拷锟斤拷锟斤拷锟斤拷卟锟酵?拷锟?1锟?7)  
	if (hTransformArg == NULL)
	{
		GDALClose((GDALDatasetH)(GDALDatasetH)pDSrc);
		return -3;
	}
	qDebug() << "has proj " << 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);

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

	//psWo->papszWarpOptions = CSLDuplicate(NULL);    
	psWo->eWorkingDataType = dataType;
	psWo->eResampleAlg = eResample;

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

	qDebug() << "GDALCreateGenImgProjTransformer has created" << 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 << endl;
	oWo.ChunkAndWarpMulti(0, 0, new_width, new_height);
	GDALFlushCache(pDDst);
	qDebug() << "ChunkAndWarpImage over" << 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", "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 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;
}

 
gdalImageComplex::gdalImageComplex(const QString& raster_path)
{
	omp_lock_t lock;
	omp_init_lock(&lock);
	omp_set_lock(&lock); 
	this->img_path = raster_path;

	GDALAllRegister(); 
	GDALDataset* rasterDataset = (GDALDataset*)(GDALOpen(raster_path.toUtf8().constData(), GA_ReadOnly));//锟斤拷只斤拷式锟斤拷取斤拷影锟斤拷
	this->width = rasterDataset->GetRasterXSize();
	this->height = rasterDataset->GetRasterYSize();
	this->band_num = rasterDataset->GetRasterCount();

	double* gt = new double[6];
	rasterDataset->GetGeoTransform(gt);
	this->gt = Eigen::MatrixXd(2, 3);
	this->gt << gt[0], gt[1], gt[2], gt[3], gt[4], gt[5];


	double a = this->gt(0, 0);
	double b = this->gt(0, 1);
	double c = this->gt(0, 2);
	double d = this->gt(1, 0);
	double e = this->gt(1, 1);
	double f = this->gt(1, 2);

	this->projection = rasterDataset->GetProjectionRef();
 
	// 斤拷投影
	GDALFlushCache((GDALDatasetH)rasterDataset);
	GDALClose((GDALDatasetH)rasterDataset);
	rasterDataset = NULL;// 指矫匡拷
	this->InitInv_gt();
	delete[] gt;
	////GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	omp_unset_lock(&lock); //锟酵放伙拷斤拷
	omp_destroy_lock(&lock); //劫伙拷斤拷


}

gdalImageComplex::~gdalImageComplex()
{
}

void gdalImageComplex::setData(Eigen::MatrixXcd data)
{
	this->data = data;
}

void gdalImageComplex::saveImage(Eigen::MatrixXcd data, int start_row, int start_col, int band_ids)
{
	omp_lock_t lock;
	omp_init_lock(&lock);
	omp_set_lock(&lock);
	if (start_row + data.rows() > this->height || start_col + data.cols() > this->width) {
		QString tip = u8"file path: " + this->img_path;
		qDebug() << tip;
		throw exception(tip.toUtf8().constData());
	}
	GDALAllRegister();
	GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
	GDALDataset* poDstDS = nullptr;
	if (exists_test(this->img_path)) {
		poDstDS = (GDALDataset*)(GDALOpen(this->img_path.toUtf8().constData(), GA_Update));
	}
	else {
		poDstDS = poDriver->Create(this->img_path.toUtf8().constData(), this->width, this->height, this->band_num, GDT_CFloat64, NULL); // 斤拷锟斤拷
		poDstDS->SetProjection(this->projection.toUtf8().constData());


		double gt_ptr[6];
		for (int i = 0; i < this->gt.rows(); i++) {
			for (int j = 0; j < this->gt.cols(); j++) {
				gt_ptr[i * 3 + j] = this->gt(i, j);
			}
		}
		poDstDS->SetGeoTransform(gt_ptr);
		delete[] gt_ptr;
	}

	int datarows = data.rows();
	int datacols = data.cols();
 
	double* databuffer = new double[data.size() * 2];
	for (int i = 0; i < data.rows(); i++) {
		for (int j = 0; j < data.cols(); j++) {
			databuffer[i * data.cols() * 2 + j * 2] = data(i, j).real();
			databuffer[i * data.cols() * 2 + j * 2 + 1] = data(i, j).imag();
		}
	}
 
	//poDstDS->RasterIO(GF_Write,start_col, start_row, datacols, datarows, databuffer, datacols, datarows, GDT_Float32,band_ids, num,0,0,0);
	poDstDS->GetRasterBand(band_ids)->RasterIO(GF_Write, start_col, start_row, datacols, datarows, databuffer, datacols, datarows, GDT_CFloat64, 0, 0);
	GDALFlushCache(poDstDS);
	delete databuffer;
	GDALClose((GDALDatasetH)poDstDS);
	//GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
 	omp_unset_lock(&lock); //锟酵放伙拷斤拷
	omp_destroy_lock(&lock); //劫伙拷斤拷

}

Eigen::MatrixXcd gdalImageComplex::getDataComplex(int start_row, int start_col, int rows_count, int cols_count, int band_ids)
{


	GDALDataset* poDataset;
	GDALAllRegister();

	// 打开TIFF文件
	poDataset = (GDALDataset*)GDALOpen(this->img_path.toUtf8().constData(), GA_ReadOnly);
	if (poDataset == nullptr) {
		QMessageBox::warning(nullptr, u8"错误", u8"无法打开文件：" + this->img_path);
		qDebug() << u8"无法打开文件：" + this->img_path;
	}

	// 获取数据集的第一个波段
	GDALRasterBand* poBand;
	poBand = poDataset->GetRasterBand(1);

	// 读取波段信息，假设是复数类型
	int nXSize = poBand->GetXSize();
	int nYSize = poBand->GetYSize();
	

	double* databuffer = new double[nXSize * nYSize * 2];
	poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count, rows_count, GDT_CFloat64, 0, 0);
	GDALClose((GDALDatasetH)poDataset);
	Eigen::MatrixXcd rasterData(nYSize, nXSize); // 使用Eigen的MatrixXcd
	for (size_t i = 0; i < nYSize; i++)
	{
		for (size_t j = 0; j < nXSize; j++) {
			rasterData(i, j) = std::complex<double>(databuffer[i * nXSize * 2 + j * 2], databuffer[i * nXSize * 2 + j * 2 + 1]);
		}
	}

	delete databuffer;
	return rasterData;
 
}

void gdalImageComplex::saveImage()
{
	this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids);
}