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









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

	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	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 std::exception(tip.toUtf8().constData());
	}
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	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();

	if (this->getDataType() == GDT_CFloat64)
	{

		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;

	}
	else if (this->getDataType() == GDT_CFloat32) {

		float* databuffer = new float[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] = float(data(i, j).real());
				databuffer[i * data.cols() * 2 + j * 2 + 1] =float( 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_CFloat32, 0, 0);
		GDALFlushCache(poDstDS);
		delete databuffer;
	}
	else {
		throw std::exception("gdalImageComplex::saveImage: data type error");
	}



	GDALClose((GDALDatasetH)poDstDS);
	GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	omp_unset_lock(&lock);	 //
	omp_destroy_lock(&lock); //
}

void gdalImageComplex::saveImage(std::shared_ptr<std::complex<double>> data, long start_row, long start_col, long rowCount, long colCount, int band_ids)
{
	omp_lock_t lock;
	omp_init_lock(&lock);
	omp_set_lock(&lock);
	if (start_row + rowCount > this->height || start_col + colCount > this->width) {
		QString tip = u8"file path: " + this->img_path;
		qDebug() << tip;
		throw std::exception(tip.toUtf8().constData());
		return;
	}
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	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;
	}

	double* databuffer = new double[rowCount * colCount * 2];
	for (long i = 0; i < rowCount; i++) {
		for (long j = 0; j < colCount; j++) {
			databuffer[i * colCount * 2 + j * 2] = data.get()[i * colCount + j].real();
			databuffer[i * colCount * 2 + j * 2 + 1] = data.get()[i * colCount + 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, colCount, rowCount,
		databuffer, colCount, rowCount, 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); //

}

void gdalImageComplex::saveImage(std::complex<double>* data, long start_row, long start_col, long rowCount, long colCount, int band_ids)
{
	omp_lock_t lock;
	omp_init_lock(&lock);
	omp_set_lock(&lock);
	if (start_row + rowCount > this->height || start_col + colCount > this->width) {
		QString tip = u8"file path: " + this->img_path;
		qDebug() << tip;
		throw std::exception(tip.toUtf8().constData());
		return;
	}
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	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;
	}

	double* databuffer = new double[rowCount * colCount * 2];
	for (long i = 0; i < rowCount; i++) {
		for (long j = 0; j < colCount; j++) {
			databuffer[i * colCount * 2 + j * 2] = data[i * colCount + j].real();
			databuffer[i * colCount * 2 + j * 2 + 1] = data[i * colCount + 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, colCount, rowCount,
		databuffer, colCount, rowCount, 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();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");

	// 打开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 = cols_count; poBand->GetXSize();
	int		nYSize = rows_count; poBand->GetYSize();
	Eigen::MatrixXcd rasterData(nYSize, nXSize); // 使用Eigen的MatrixXcd
	if (this->getDataType() == GDT_CFloat64)
	{
		long long pixelCount =long long( nXSize) *long long( nYSize);
		double* databuffer = new double[pixelCount * 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);

		for (long long i = 0; i < nYSize; i++) {
			for (long long 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;
	}
	else if(this->getDataType()==GDT_CFloat32)
	{
		long long pixelCount = long long(nXSize) * long long(nYSize);
		float* databuffer = new float[pixelCount * 2];
		poBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, databuffer, cols_count, rows_count, GDT_CFloat32, 0, 0);
		GDALClose((GDALDatasetH)poDataset);

		for (long long i = 0; i < nYSize; i++) {
			for (long long 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;
}

std::shared_ptr<std::complex<double>> gdalImageComplex::getDataComplexSharePtr(int start_row, int start_col, int rows_count, int cols_count, int band_ids)
{
	GDALDataset* poDataset;
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");

	// 打开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);

	// 读取波段信息，假设是复数类型


	double* databuffer = new double[rows_count * cols_count * 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);

	std::shared_ptr<std::complex<double>> rasterData(new std::complex<double>[rows_count * cols_count], delArrPtr);

	for (size_t i = 0; i < rows_count; i++) {
		for (size_t j = 0; j < cols_count; j++) {
			rasterData.get()[i * cols_count + j] = std::complex<double>(databuffer[i * cols_count * 2 + j * 2],
				databuffer[i * cols_count * 2 + j * 2 + 1]);
		}
	}

	delete[] databuffer;
	return rasterData;
}

void gdalImageComplex::saveComplexImage()
{
	this->saveImage(this->data, this->start_row, this->start_col, this->data_band_ids);
}
void gdalImageComplex::savePreViewImage()
{
	qDebug() << "void gdalImageComplex::savePreViewImage()";
	Eigen::MatrixXd data_abs = Eigen::MatrixXd::Zero(this->height, this->width);
	data_abs = (this->data.array().real().pow(2) + this->data.array().imag().pow(2))
		.array()
		.log10() * 10.0; // 计算振幅

	double				   min_abs = data_abs.minCoeff();		 // 最大值
	double				   max_abs = data_abs.maxCoeff();		 // 最小值
	double				   delta = (max_abs - min_abs) / 1000; // 1000分位档
	Eigen::MatrixX<size_t> data_idx =
		((data_abs.array() - min_abs).array() / delta).array().floor().cast<size_t>();

	// 初始化
	double hist[1001];
	for (size_t i = 0; i < 1001; i++) {
		hist[i] = 0; // 初始化
	}
	for (size_t i = 0; i < this->height; i++) {
		for (size_t j = 0; j < this->width; j++) {
			hist[data_idx(i, j)]++;
		}
	}

	// 统计
	size_t count = this->height * this->width;
	double precent = 0;
	size_t curCount = 0;
	double pre2 = 0;
	bool   findprec_2 = true;
	double pre98 = 0;
	bool   findprec_98 = true;
	for (size_t i = 0; i < 1001; i++) {
		precent = precent + hist[i];
		if (findprec_2 & precent / count > 0.02) {
			pre2 = i * delta + min_abs;
			findprec_2 = false;
		}
		if (findprec_98 & precent / count > 0.98) {
			pre98 = (i - 1) * delta + min_abs;
			findprec_98 = false;
		}
	}
	// 拉伸
	delta = (pre98 - pre2) / 200;
	data_idx =
		((data_abs.array() - pre2).array() / delta).array().floor().cast<size_t>();

	for (size_t i = 0; i < this->height; i++) {
		for (size_t j = 0; j < this->width; j++) {
			if (data_idx(i, j) < 0) {
				data_idx(i, j) = 0;
			}
			else if (data_idx(i, j) > 255) {
				data_idx(i, j) = 255;
			}
			else {

			}
		}
	}

	// 赋值
	QString	  filePath = this->img_path;
	QFile	  file(filePath);
	QFileInfo fileInfo(file);

	QString	  directory = fileInfo.absolutePath();
	qDebug() << "文件所在目录：" << directory;
	QString baseName = fileInfo.completeBaseName();
	qDebug() << "无后缀文件名：" << baseName;

	// 创建文件路径
	QString			previewImagePath = JoinPath(directory, baseName + "_preview.png");
	cv::Mat			img(this->height, this->width, CV_8U, cv::Scalar(0));

	for (size_t i = 0; i < this->height; i++) {
		for (size_t j = 0; j < this->width; j++) {
			img.at<uchar>(i, j) = (uchar)(data_idx(i, j));
		}
	}
	//std::string outimgpath=previewImagePath.toUtf8().data();
	cv::imwrite(previewImagePath.toUtf8().data(), 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();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
	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(0); // 回波部分
	//}
	GDALFlushCache((GDALDatasetH)poDstDS);
	GDALClose((GDALDatasetH)poDstDS);
	GDALDestroy(); // or, DllMain at DLL_PROCESS_DETACH
	gdalImageComplex result_img(img_path);
	return result_img;
}

gdalImageComplex  CreategdalImageComplexNoProj(const QString& img_path, int height, int width, int band_num, bool overwrite)
{

	// 创建新文件
	omp_lock_t lock;
	omp_init_lock(&lock);
	omp_set_lock(&lock);
	GDALAllRegister();
	CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");

	GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
	GDALDataset* poDstDS = (poDriver->Create(img_path.toUtf8().constData(), width, height, band_num, GDT_CFloat64, NULL));
	GDALFlushCache((GDALDatasetH)poDstDS);
	GDALClose((GDALDatasetH)poDstDS);
	//poDstDS.reset();
	omp_unset_lock(&lock);	 //  
	omp_destroy_lock(&lock); //  

	gdalImageComplex result_img(img_path);
	return result_img;
}

gdalImageComplex CreateEchoComplex(const QString& img_path, int height, int width, int band_num)
{
	// 创建图像
	Eigen::MatrixXd gt = Eigen::MatrixXd::Zero(2, 3);
	//Xgeo = GeoTransform[0] + Xpixel * GeoTransform[1] + Ypixel * GeoTransform[2]
	//Ygeo = GeoTransform[3] + Xpixel * GeoTransform[4] + Ypixel * GeoTransform[5]
	// X
	gt(0, 0) = 0; gt(0, 2) = 1; gt(0, 2) = 0;
	gt(1, 0) = 0; gt(1, 1) = 0; gt(1, 2) = 1;
	// Y
	QString projection = "";
	gdalImageComplex echodata = CreategdalImageComplex(img_path, height, width, 1, gt, projection, false, true);
	return echodata;

}