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

ErrorCode MergeRasterProcess(QVector<QString> filepaths, QString outfileptah, QString mainString, MERGEMODE  mergecode, bool isENVI, ShowProessAbstract* dia)
{
	// 参数检查
	if (!isExists(mainString)) {
		qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::FILENOFOUND)) << "\t" << mainString;
		return ErrorCode::FILENOFOUND;
	}
	else {}
	gdalImage mainimg(mainString);
	QVector<gdalImage> imgdslist(filepaths.count());
	for (long i = 0; i < filepaths.count(); i++) {
		if (!isExists(filepaths[i])) {
			qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::FILENOFOUND)) << "\t" << filepaths[i];
			return ErrorCode::FILENOFOUND;
		}
		else {
			imgdslist[i] = gdalImage(filepaths[i]);
			if (imgdslist[i].band_num != mainimg.band_num) {
				qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::RASTERBAND_NOTEQUAL)) << "\t" << imgdslist[i].band_num << " != " << mainimg.band_num;
				return ErrorCode::RASTERBAND_NOTEQUAL;
			}
		}
	}

	// 检查坐标系是否统一
	long EPSGCode = GetEPSGFromRasterFile(mainString);
	long tempCode = 0;
	for (long i = 0; i < filepaths.count(); i++) {
		tempCode = GetEPSGFromRasterFile(filepaths[i]);
		if (EPSGCode != tempCode) {
			qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::EPSGCODE_NOTSAME)) << "\t" << EPSGCode << "!=" << tempCode;
			return ErrorCode::EPSGCODE_NOTSAME;
		}
	}

	// 检查影像类型是否统一
	GDALDataType mainType = mainimg.getDataType();
	for (long i = 0; i < imgdslist.count(); i++) {
		if (mainType != imgdslist[i].getDataType()) {
			qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::RASTER_DATETYPE_NOTSAME)) << "\t" << mainType << "!=" << imgdslist[i].getDataType();
			return ErrorCode::RASTER_DATETYPE_NOTSAME;
		}
	}

	Eigen::MatrixXd maingt = mainimg.getGeoTranslation();
	Eigen::MatrixXd rgt = Eigen::MatrixXd::Zero(2, 3);
	RasterExtend mainExtend = mainimg.getExtend();
	rgt(0, 1) = (mainExtend.max_x - mainExtend.min_x) / (mainimg.width - 1); //dx
	rgt(1, 2) = -1 * std::abs(((mainExtend.max_y - mainExtend.min_y) / (mainimg.height - 1)));//dy
	QVector<RasterExtend> extendlist(imgdslist.count());
	for (long i = 0; i < imgdslist.count(); i++) {
		extendlist[i] = imgdslist[i].getExtend();
		mainExtend.min_x = mainExtend.min_x < extendlist[i].min_x ? mainExtend.min_x : extendlist[i].min_x;
		mainExtend.max_x = mainExtend.max_x > extendlist[i].max_x ? mainExtend.max_x : extendlist[i].max_x;
		mainExtend.min_y = mainExtend.min_y < extendlist[i].min_y ? mainExtend.min_y : extendlist[i].min_y;
		mainExtend.max_y = mainExtend.max_y > extendlist[i].max_y ? mainExtend.max_y : extendlist[i].max_y;
	}

	rgt(0, 0) = mainExtend.min_x;
	rgt(1, 0) = mainExtend.max_y;

	// 计算数量

	long width = std::ceil((mainExtend.max_x - mainExtend.min_x) / rgt(0, 1) + 1);
	long height = std::ceil(std::abs((mainExtend.min_y - mainExtend.max_y) / rgt(1, 2)) + 1);
	OGRSpatialReference oSRS;
	if (oSRS.importFromEPSG(EPSGCode) != OGRERR_NONE) {
		qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::EPSGCODE_NOTSUPPORT)) << "\t" << EPSGCode;
		return ErrorCode::EPSGCODE_NOTSUPPORT;
	}

	gdalImage resultImage = CreategdalImage(outfileptah, height, width, mainimg.band_num, rgt, EPSGCode, mainType, true, true, isENVI);

	QString resultMaskString = addMaskToFileName(outfileptah, QString("_MASK"));
	gdalImage maskImage = CreategdalImage(resultMaskString, height, width, 1, rgt, EPSGCode, GDT_Int32, true, true, isENVI);

	// 初始化
	long resultline = Memory1MB * 500 / 8 / resultImage.width;
	resultline = resultline < 10000 ? resultline : 10000; // 最多100行
	resultline = resultline > 0 ? resultline : 2;
	long bandnum = resultImage.band_num + 1;
	long starti = 0;
	long rasterCount = imgdslist.count();



	QProgressDialog progressDialog(u8"初始化影像", u8"终止", 0, resultImage.height);
	progressDialog.setWindowTitle(u8"初始化影像");
	progressDialog.setWindowModality(Qt::WindowModal);
	progressDialog.setAutoClose(true);
	progressDialog.setValue(0);
	progressDialog.setMaximum(resultImage.height);
	progressDialog.setMinimum(0);
	progressDialog.show();


	for (starti = 0; starti < resultImage.height; starti = starti + resultline) {
		long blocklines = resultline;
		blocklines = starti + blocklines < resultImage.height ? blocklines : resultImage.height - starti;
		for (long b = 1; b < bandnum; b++) {
			Eigen::MatrixXd data = resultImage.getData(starti, 0, blocklines, resultImage.width, b);
			Eigen::MatrixXi maskdata = maskImage.getDatai(starti, 0, blocklines, resultImage.width, b);
			data = data.array() * 0;
			maskdata = maskdata.array() * 0;
			resultImage.saveImage(data, starti, 0, b);
			maskImage.saveImage(maskdata, starti, 0, b);
		}
		if (nullptr != dia) {
			dia->showProcess(starti * 1.0 / resultImage.height, u8"初始化影像数据");
		}
		progressDialog.setValue(starti + blocklines);
	}
	progressDialog.close();



	switch (mergecode)
	{
	case MERGE_GEOCODING:
		return MergeRasterInGeoCoding(imgdslist, resultImage, maskImage, dia);
	default:
		break;
	}


	return ErrorCode::SUCCESS;
}

ErrorCode MergeRasterInGeoCoding(QVector<gdalImage> imgdslist, gdalImage resultimg, gdalImage maskimg, ShowProessAbstract* dia)
{
	omp_set_num_threads(Paral_num_thread);
	// 逐点合并计算
	QVector<RasterExtend> extendlist(imgdslist.count());
	for (long i = 0; i < imgdslist.count(); i++) {
		extendlist[i] = imgdslist[i].getExtend();
		imgdslist[i].InitInv_gt();
	}

	// 分块计算
	long resultline = Memory1MB * 1000 / 8 / resultimg.width;
	resultline = resultline < 300 ? resultline : 300; // 最多100行

	long bandnum = resultimg.band_num + 1;
	long starti = 0;
	long rasterCount = imgdslist.count();

	long processNumber = 0;
	QProgressDialog progressDialog(u8"合并影像", u8"终止", 0, resultimg.height);
	progressDialog.setWindowTitle(u8"合并影像");
	progressDialog.setWindowModality(Qt::WindowModal);
	progressDialog.setAutoClose(true);
	progressDialog.setValue(0);
	progressDialog.setMaximum(resultimg.height);
	progressDialog.setMinimum(0);
	progressDialog.show();
	omp_lock_t lock;
	omp_init_lock(&lock);

#pragma omp parallel for
	for (starti = 0; starti < resultimg.height; starti = starti + resultline) {
		long blocklines = resultline;
		blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti;

		long rid = starti;
		long cid = 0;

		Landpoint pp = { 0,0,0 };
		Landpoint lp = { 0,0,0 };

		for (long ir = 0; ir < rasterCount; ir++) {// 影像
			long minRid = imgdslist[ir].height;
			long maxRid = 0;

			Eigen::MatrixXd ridlist = resultimg.getData(starti, 0, blocklines, resultimg.width, 1);
			ridlist = ridlist.array() * 0;
			Eigen::MatrixXd cidlist = ridlist.array() * 0;

			for (long i = 0; i < blocklines; i++) {// 行号
				rid = starti + i;
				for (long j = 0; j < resultimg.width; j++) {// 列号
					cid = j;
					resultimg.getLandPoint(rid, cid, 0, pp);
					lp = imgdslist[ir].getRow_Col(pp.lon, pp.lat); // 获取点坐标
					ridlist(i, j) = lp.lat;
					cidlist(i, j) = lp.lon;
				}
			}

			//ImageShowDialogClass* dialog = new ImageShowDialogClass;
			//dialog->show();
			//dialog->load_double_MatrixX_data(cidlist, u8"");

			//dialog->exec();


			if (ridlist.maxCoeff() < 0 || ridlist.minCoeff() >= imgdslist[ir].height) {
				continue;
			}

			if (cidlist.maxCoeff() < 0 || cidlist.minCoeff() >= imgdslist[ir].width) {
				continue;
			}

			minRid = std::floor(ridlist.minCoeff());
			maxRid = std::ceil(ridlist.maxCoeff());
			minRid = minRid < 0 ? 0 : minRid;
			maxRid = maxRid < imgdslist[ir].height ? maxRid : imgdslist[ir].height - 1;

			long rowlen = maxRid - minRid + 1;
			if (rowlen <= 0) {
				continue;
			}
			// 获取分配代码
			Landpoint p0{ 0,0,0 }, p11{ 0,0,0 }, p21{ 0,0,0 }, p12{ 0,0,0 }, p22{ 0,0,0 }, p{ 0,0,0 };


			long rowcount = 0;
			long colcount = 0;
			double ridtemp = 0, cidtemp = 0;

			long lastr = 0, nextr = 0;
			long lastc = 0, nextc = 0;

			double r0 = 0, c0 = 0;

			for (long b = 1; b < bandnum; b++) {
				Eigen::MatrixXd resultdata = resultimg.getData(starti, 0, blocklines, resultimg.width, b);
				Eigen::MatrixXi resultmask = maskimg.getDatai(starti, 0, blocklines, resultimg.width, b);
				Eigen::MatrixXd data = imgdslist[ir].getData(minRid, 0, rowlen, imgdslist[ir].width, b);

				double nodata = imgdslist[ir].getNoDataValue(b);
				for (long ii = 0; ii < data.rows(); ii++) {
					for (long jj = 0; jj < data.cols(); jj++) {
						if (std::abs(data(ii, jj) - nodata) < 1e-6) {
							data(ii, jj) = 0;
						}
					}
				}
				rowcount = ridlist.rows();
				colcount = ridlist.cols();
				double Bileanervalue = 0;
				for (long i = 0; i < rowcount; i++) {
					for (long j = 0; j < colcount; j++) {
						ridtemp = ridlist(i, j);
						cidtemp = cidlist(i, j);

						lastr = std::floor(ridtemp);
						nextr = std::ceil(ridtemp);
						lastc = std::floor(cidtemp);
						nextc = std::ceil(cidtemp);

						if (lastr < 0 || lastr >= imgdslist[ir].height
							|| nextr < 0 || nextr >= imgdslist[ir].height
							|| lastc < 0 || lastc >= imgdslist[ir].width
							|| nextc < 0 || nextc >= imgdslist[ir].width) {
							continue;
						}
						else {}

						r0 = ridtemp - std::floor(ridtemp);
						c0 = cidtemp - std::floor(cidtemp);

						lastr = lastr - minRid;
						nextr = nextr - minRid;

						p0 = Landpoint{ c0,r0,0 };
						p11 = Landpoint{ 0,0,data(lastr,lastc) };
						p21 = Landpoint{ 0,1,data(nextr,lastc) };
						p12 = Landpoint{ 1,0,data(lastr,nextc) };
						p22 = Landpoint{ 1,1,data(nextr,nextc) };
						Bileanervalue = Bilinear_interpolation(p0, p11, p21, p12, p22);
						if (std::abs(Bileanervalue) < 1e-6 || resultmask(i, j) > 0) {
							continue;
						}
						resultdata(i, j) = resultdata(i, j) + Bileanervalue;
						resultmask(i, j) = resultmask(i, j) + 1;
					}
				}
				resultimg.saveImage(resultdata, starti, 0, b);
				maskimg.saveImage(resultmask, starti, 0, b);
			}
		}

		omp_set_lock(&lock);
		processNumber = processNumber + blocklines;
		qDebug() << "\rprocess bar:\t" << processNumber * 100.0 / resultimg.height << " % " << "\t\t\t";
		if (nullptr != dia) {
			dia->showProcess(processNumber * 1.0 / resultimg.height, u8"合并图像");
		}
		if (progressDialog.maximum() <= processNumber) {
			processNumber = progressDialog.maximum() - 1;
		}
		progressDialog.setValue(processNumber);
		omp_unset_lock(&lock);
	}
	omp_destroy_lock(&lock);

	progressDialog.setWindowTitle(u8"影像掩膜");
	progressDialog.setLabelText(u8"影像掩膜");
	for (starti = 0; starti < resultimg.height; starti = starti + resultline) {
		long blocklines = resultline;
		blocklines = starti + blocklines < resultimg.height ? blocklines : resultimg.height - starti;
		for (long b = 1; b < bandnum; b++) {
			Eigen::MatrixXd data = resultimg.getData(starti, 0, blocklines, resultimg.width, b);
			Eigen::MatrixXd maskdata = maskimg.getData(starti, 0, blocklines, maskimg.width, b);

			for (long i = 0; i < data.rows(); i++) {
				for (long j = 0; j < data.cols(); j++) {
					if (maskdata(i, j) == 0) {
						data(i, j) = -9999;
						continue;
					}
					data(i, j) = data(i, j) / maskdata(i, j);
				}
			}

			resultimg.saveImage(data, starti, 0, b);
			maskimg.saveImage(maskdata, starti, 0, b);
		}
		if (nullptr != dia) {
			dia->showProcess((starti + blocklines) * 1.0 / resultimg.height, u8"影像掩膜");
		}
		progressDialog.setValue(starti + blocklines);
	}


	progressDialog.close();
	return ErrorCode::SUCCESS;
}


void MergeTiffs(QList<QString> inputFiles, QString outputFile) {
	GDALAllRegister();

	if (inputFiles.isEmpty()) {
		fprintf(stderr, "No input files provided.\n");
		return;
	}

	// Open the first file to determine the data type and coordinate system
	GDALDataset* poFirstDS = (GDALDataset*)GDALOpen(inputFiles.first().toUtf8().constData(), GA_ReadOnly);
	if (poFirstDS == nullptr) {
		fprintf(stderr, "Failed to open the first file %s\n", inputFiles.first().toUtf8().constData());
		return;
	}

	double adfGeoTransform[6];
	CPLErr eErr = poFirstDS->GetGeoTransform(adfGeoTransform);
	if (eErr != CE_None) {
		fprintf(stderr, "Failed to get GeoTransform for the first file %s\n", inputFiles.first().toUtf8().constData());
		GDALClose(poFirstDS);
		return;
	}

	int nXSize = 0;
	int nYSize = 0;
	double minX = std::numeric_limits<double>::max();
	double minY = std::numeric_limits<double>::max();
	double maxX = std::numeric_limits<double>::lowest();
	double maxY = std::numeric_limits<double>::lowest();
	double pixelWidth = adfGeoTransform[1];
	double pixelHeight = adfGeoTransform[5];

	// Determine the bounding box and size of the output raster
	for (const QString& inputFile : inputFiles) {
		GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(inputFile.toUtf8().constData(), GA_ReadOnly);
		if (poSrcDS == nullptr) {
			fprintf(stderr, "Failed to open %s\n", inputFile.toUtf8().constData());
			continue;
		}

		double adfThisTransform[6];
		eErr = poSrcDS->GetGeoTransform(adfThisTransform);
		if (eErr != CE_None) {
			fprintf(stderr, "Failed to get GeoTransform for %s\n", inputFile.toUtf8().constData());
			GDALClose(poSrcDS);
			continue;
		}

		minX = std::min(minX, adfThisTransform[0]);
		minY = std::min(minY, adfThisTransform[3] + adfThisTransform[5] * poSrcDS->GetRasterYSize());
		maxX = std::max(maxX, adfThisTransform[0] + adfThisTransform[1] * poSrcDS->GetRasterXSize());
		maxY = std::max(maxY, adfThisTransform[3]);

		GDALClose(poSrcDS);
	}

	nXSize = static_cast<int>(std::ceil((maxX - minX) / pixelWidth));
	nYSize = static_cast<int>(std::ceil((maxY - minY) / (-pixelHeight)));

	GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (poDriver == nullptr) {
		fprintf(stderr, "GTiff driver not available.\n");
		GDALClose(poFirstDS);
		return;
	}

	char** papszOptions = nullptr;
	GDALDataset* poDstDS = poDriver->Create(outputFile.toUtf8().constData(), nXSize, nYSize, 1, poFirstDS->GetRasterBand(1)->GetRasterDataType(), papszOptions);
	if (poDstDS == nullptr) {
		fprintf(stderr, "Creation of output file failed.\n");
		GDALClose(poFirstDS);
		return;
	}

	poDstDS->SetGeoTransform(adfGeoTransform);

	const OGRSpatialReference* oSRS = poFirstDS->GetSpatialRef();
	poDstDS->SetSpatialRef(oSRS);

	float fillValue = std::numeric_limits<float>::quiet_NaN();
	void* pafScanline = CPLMalloc(GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * nXSize);
	memset(pafScanline, 0, GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * nXSize);

	// Initialize all pixels to NaN
	for (int iY = 0; iY < nYSize; ++iY) {
		GDALRasterBand* poBand = poDstDS->GetRasterBand(1);
		poBand->RasterIO(GF_Write, 0, iY, nXSize, 1, pafScanline, nXSize, 1, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);
	}

	CPLFree(pafScanline);

	// Read each source image and write into the destination image
	for (const QString& inputFile : inputFiles) {
		GDALDataset* poSrcDS = (GDALDataset*)GDALOpen(inputFile.toUtf8().constData(), GA_ReadOnly);
		if (poSrcDS == nullptr) {
			fprintf(stderr, "Failed to open %s\n", inputFile.toUtf8().constData());
			continue;
		}

		double adfThisTransform[6];
		poSrcDS->GetGeoTransform(adfThisTransform);

		int srcXSize = poSrcDS->GetRasterXSize();
		int srcYSize = poSrcDS->GetRasterYSize();

		int dstXOffset = static_cast<int>(std::round((adfThisTransform[0] - minX) / pixelWidth));
		int dstYOffset = static_cast<int>(std::round((maxY - adfThisTransform[3]) / (-pixelHeight)));

		GDALRasterBand* poSrcBand = poSrcDS->GetRasterBand(1);
		GDALRasterBand* poDstBand = poDstDS->GetRasterBand(1);

		void* pafBuffer = CPLMalloc(GDALGetDataTypeSizeBytes(poFirstDS->GetRasterBand(1)->GetRasterDataType()) * srcXSize * srcYSize);
		poSrcBand->RasterIO(GF_Read, 0, 0, srcXSize, srcYSize, pafBuffer, srcXSize, srcYSize, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);

		poDstBand->RasterIO(GF_Write, dstXOffset, dstYOffset, srcXSize, srcYSize, pafBuffer, srcXSize, srcYSize, poFirstDS->GetRasterBand(1)->GetRasterDataType(), 0, 0);

		CPLFree(pafBuffer);
		GDALClose(poSrcDS);
	}

	GDALClose(poDstDS);
	GDALClose(poFirstDS);

}