#include <opencv/highgui.h>
#include <math.h>
#include "GeoTiffMerger.hxx"
#include "GeoTiffReader.hxx"
#include "GeoTiffImage.hxx"
#include "GeoTiffWriter.hxx"
#include "FileUtils.hxx"

GeoTiffMerger::GeoTiffMerger()
{
	_sourceDir[0] = '\0';
	_destFileName[0] = '\0';

	_minModelCoordX = 10000.0f;
	_maxModelCoordX = -10000.0f;
	_minModelCoordY = 10000.0f;
	_maxModelCoordY = -10000.0f;

	_avgXScale = 1.0;
	_avgYScale = 1.0;

	_imageScaleFactor = 0.1;

	_depth = 32;
	_channels = 4;

	_maxVal = -10000.0;
	_minVal = 10000.0f;

	_targetImage = nullptr;

}

void GeoTiffMerger::estimateParameters()
{
	_targetImage = new GeoTiffImage();

	// iterate tiff files and get dimensions;
	//QDir dir(_sourceDir);
	//dir.setFilter(QDir::Files | QDir::Hidden | QDir::NoSymLinks);
	//dir.setSorting(QDir::Name | QDir::Reversed);

	//QStringList filters;
	//filters << "*.tif";
	//dir.setNameFilters(filters);

	//QFileInfoList list = dir.entryInfoList();
	//QString str;
        std::vector<std::string> list;
        ListFiles(_sourceDir, "tif", list);

	double sumScaleX = 0.0f;
	double sumScaleY = 0.0f;

	for (int i = 0; i < list.size(); i++)
	{
		//QFileInfo fileInfo = list.at(i);
		//str = fileInfo.fileName();
                std::string str = list.at(i);

		char fname[256];
		//sprintf(fname, "%s/%s", _sourceDir, str.toStdString().c_str());
                sprintf(fname,"%s", str.c_str());

		printf("[ %04d/%04d ]:%s\n", i + 1, list.size(), fname);

		GeoTiffImage geoImage;
		GeoTiffReader reader;


		reader.setGeoTiffFileName(fname);
		reader.setGeoTiffImage(&geoImage);
		reader.setJustReadHeader();

		reader.readImage();

		double sx, sy, sz;
		double mx, my, mz;

		geoImage.getGeoScale(sx, sy, sz);
		printf("GeoScale:%2.15f,%2.15f,%2.15f\n", sx, sy, sz);
		geoImage.getGeoModelCoord(mx, my, mz);
		printf("GeoModePoint:%f,%f,%f\n", mx, my, mz);

		if (i == 0)
		{
			_depth = geoImage.getDepth();
			_channels = geoImage.getChannels();

			_targetImage->setDepth(_depth);
			_targetImage->setChannels(_channels);
			// _targetImage->setGeoKeyDirectory(geoImage.getGeoKeyDirectory());
			// _targetImage->setGeoDoubleParams(geoImage.getGeoDoubleParams());
			// _targetImage->setGeoAsciiParams(geoImage.getGeoAsciiParams());
			// _targetImage->setGeoNoData(geoImage.getGeoNoData());
		}
		else
		{
			assert(_depth == geoImage.getDepth());
			assert(_channels = geoImage.getChannels());
		}

		sumScaleX += sx;
		sumScaleY += sy;

		double tMinX = mx;
		double tMaxX = mx + sx * geoImage.getWidth();

		double tMinY = my;
		double tMaxY = my + sy * geoImage.getHeight();

		if (tMinX < _minModelCoordX) _minModelCoordX = tMinX;
		if (tMaxX > _maxModelCoordX) _maxModelCoordX = tMaxX;

		if (tMinY < _minModelCoordY) _minModelCoordY = tMinY;
		if (tMaxY > _maxModelCoordY) _maxModelCoordY = tMaxY;

		printf("Model Range X:%03.15f,%03.15f\n", _minModelCoordX, _maxModelCoordX);
		printf("Model Range Y:%03.15f,%03.15f\n", _minModelCoordY, _maxModelCoordY);
	}

	_avgXScale = sumScaleX / list.size();
	_avgYScale = sumScaleY / list.size();

	printf("Avg Scale: %03.15f\t%03.15f\n", _avgXScale, _avgYScale);

	int mergedWidth = (_maxModelCoordX - _minModelCoordX) / _avgXScale;
	int mergedHeight = (_maxModelCoordY - _minModelCoordY) / _avgYScale;
	printf("Merged Image Widtd: %f\n", mergedWidth );
	printf("Merged Image Height: %f\n",mergedHeight );

	double controlPoint[6];
	controlPoint[3] = _minModelCoordX;
	controlPoint[4] = _maxModelCoordY;
	controlPoint[0] = controlPoint[1] = controlPoint[2] = controlPoint[5] = 0.0f;

	_targetImage->setWidth(mergedWidth);
	_targetImage->setHeight(mergedHeight);
	_targetImage->setGeoScale(_avgXScale, _avgYScale, 0.0);
	_targetImage->setGeoControlPoint(controlPoint);

}

void GeoTiffMerger::mergeFiles()
{
	// Initialize target image
	int targetWidth = (_maxModelCoordX - _minModelCoordX) / _avgXScale * _imageScaleFactor + 0.5f;
	int targetHeight = (_maxModelCoordY - _minModelCoordY) / _avgYScale * _imageScaleFactor + 0.5f;

	double targetScaleX = _avgXScale / _imageScaleFactor;
	double targetScaleY = _avgYScale / _imageScaleFactor;

	printf("Target Image Size: %d, %d\n", targetWidth, targetHeight);

	_targetImage->setWidth(targetWidth);
	_targetImage->setHeight(targetHeight);
	_targetImage->setGeoScale(targetScaleX, targetScaleY, 0.0f);


	IplImage* targetImage = cvCreateImage(cvSize(targetWidth, targetHeight), _depth, _channels);
	//cvZero(targetImage);
	cvSet(targetImage, cvScalar(-9999, -9999, -9999, -9999));

	// QDir dir(_sourceDir);
	// dir.setFilter(QDir::Files | QDir::Hidden | QDir::NoSymLinks);
	// dir.setSorting(QDir::Name | QDir::Reversed);

	// QStringList filters;
	// filters << "*.tif";
	// dir.setNameFilters(filters);

	// QFileInfoList list = dir.entryInfoList();
	// QString str;
	std::vector<std::string> fileList;
	ListFiles(_sourceDir,"tif",fileList);

	for (int i = 0; i < fileList.size(); i++)
	{
		//QFileInfo fileInfo = list.at(i);
		//str = fileInfo.fileName();

		std::string filePath = fileList.at(i);
		//char fname[256];
		//sprintf(fname, "%s/%s", _sourceDir, str.toStdString().c_str());

		printf("[ %04d/%04d ]:%s\n", i + 1, fileList.size(), filePath.c_str());

		GeoTiffImage geoImage;
		GeoTiffReader reader;


		reader.setGeoTiffFileName(filePath.c_str());
		reader.setGeoTiffImage(&geoImage);
		reader.setJustReadHeader();

		reader.readImage();

		double sx, sy, sz;
		double mx, my, mz;


		// mx, my are left bottom coordinate of geo coordinate.
		geoImage.getGeoScale(sx, sy, sz);
		printf("GeoScale:%2.15f,%2.15f,%2.15f\n", sx, sy, sz);
		geoImage.getGeoModelCoord(mx, my, mz);
		printf("GeoModePoint:%f,%f,%f\n", mx, my, mz);

		// transform my to top coordinate of geo coordinate.
		my = my + geoImage.getHeight() * sy;

		int imgPosX = floor((mx - _minModelCoordX) / targetScaleX);
		//int imgPosY = floor((my - _minModelCoordY) / targetScaleY);
		int imgPosY = ceil((_maxModelCoordY - my) / targetScaleY);

		if (imgPosY < 0) imgPosY = 0;

		//int tmpImageWidth = ceil(geoImage.getImageWidth() * sx / targetScaleX);
		//int tmpImageHeight = ceil(geoImage.getImageHeight() * sy / targetScaleY);
		int tmpImageWidth = geoImage.getWidth() * sx / targetScaleX + 0.5;
		int tmpImageHeight = geoImage.getHeight() * sy / targetScaleY+ 0.5;

		if (imgPosX + tmpImageWidth >= targetImage->width)
			tmpImageWidth = targetImage->width - imgPosX - 1;

		if (imgPosY + tmpImageHeight >= targetImage->height)
			tmpImageHeight = targetImage->height - imgPosY - 1;

		IplImage* tmpImage = cvCreateImage(cvSize(tmpImageWidth, tmpImageHeight), _depth , _channels);
		IplImage* tifImage = cvLoadImage(filePath.c_str(), CV_LOAD_IMAGE_UNCHANGED);
		
		cvResize(tifImage, tmpImage);	

		assert(((imgPosX + tmpImageWidth) < targetImage->width) && ((imgPosY + tmpImageHeight) < targetImage->height) && "Out of target range!!!");

		cvSetImageROI(targetImage, cvRect(imgPosX, imgPosY, tmpImageWidth, tmpImageHeight));

		blendImageBlock(tmpImage, targetImage);

		cvResetImageROI(targetImage);

		cvReleaseImage(&tmpImage);
		cvReleaseImage(&tifImage);
	}

	//cvSaveImage("F:/TargetImage.png", targetImage);

	//IplImage* tmpImage = cvCreateImage(cvGetSize(targetImage), IPL_DEPTH_8U, 1);
	//cvConvertScale(targetImage, tmpImage, 255.0f / _maxVal);

	//cvSaveImage(_destFileName, tmpImage);
	//exportTxtDEM(targetImage);
	if (_channels == 1)
	{
		GeoTiffWriter writer;
		writer.setGeoTiffImage(_targetImage);
		writer.setFileName(_destFileName);
		writer.write(targetImage, _targetImage);
	}
	else
	{
		//IplImage* tmpImage = cvCreateImage(cvGetSize(targetImage), IPL_DEPTH_8U, 1);
		//cvConvertScale(targetImage, tmpImage, 255.0f / _maxVal);

		//cvSaveImage(_destFileName, tmpImage);
		cvSaveImage(_destFileName, targetImage);
	}

	cvReleaseImage(&targetImage);
}

void GeoTiffMerger::blendImageBlock(IplImage* source, IplImage* target)
{
	if (_channels == 4 && _depth == 8)
		blendRGBAImage(source, target);
	else if (_channels == 1 && _depth == 32)
		blendGrayImage(source, target);
}

void GeoTiffMerger::blendRGBAImage(IplImage* source, IplImage* target)
{
	int w = source->width;
	int h = source->height;

	IplImage* sR = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 1);
	IplImage* sG = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 1);
	IplImage* sB = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 1);
	IplImage* sA = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 1);

	IplImage* sF = cvCreateImage(cvSize(w, h), IPL_DEPTH_32F, 1);

	IplImage* tR = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 1);
	IplImage* tG = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 1);
	IplImage* tB = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 1);
	IplImage* tA = cvCreateImage(cvSize(w, h), IPL_DEPTH_8U, 1);

	IplImage* tF = cvCreateImage(cvSize(w, h), IPL_DEPTH_32F, 1);

	IplImage* asum = cvCreateImage(cvSize(w, h), IPL_DEPTH_32F, 1);

	IplImage* tmpChannel = cvCreateImage(cvSize(w, h), IPL_DEPTH_32F, 1);


	cvSplit(source, sB, sG, sR, sA);
	cvSplit(target, tB, tG, tR, tA);

	cvConvertScale(sA, sF, 1.0f / 255.0f);
	cvConvertScale(tA, tF, 1.0f / 255.0f);

	cvAdd(sF, tF, asum);
	cvDiv(sF, asum, sF);
	cvDiv(tF, asum, tF);

	cvConvertScale(sB, tmpChannel);
	cvMul(tmpChannel, sF, tmpChannel);
	cvConvertScale(tmpChannel, sB);

	cvConvertScale(sG, tmpChannel);
	cvMul(tmpChannel, sF, tmpChannel);
	cvConvertScale(tmpChannel, sG);

	cvConvertScale(sR, tmpChannel);
	cvMul(tmpChannel, sF, tmpChannel);
	cvConvertScale(tmpChannel, sR);

	cvMerge(sB, sG, sR, sA, source);

	cvConvertScale(tB, tmpChannel);
	cvMul(tmpChannel, tF, tmpChannel);
	cvConvertScale(tmpChannel, tB);

	cvConvertScale(tG, tmpChannel);
	cvMul(tmpChannel, tF, tmpChannel);
	cvConvertScale(tmpChannel, tG);

	cvConvertScale(tR, tmpChannel);
	cvMul(tmpChannel, tF, tmpChannel);
	cvConvertScale(tmpChannel, tR);

	cvMerge(tB, tG, tR, tA, target);

	cvAdd(source, target, target);

	cvReleaseImage(&sR);
	cvReleaseImage(&sB);
	cvReleaseImage(&sG);
	cvReleaseImage(&sA);

	cvReleaseImage(&tR);
	cvReleaseImage(&tB);
	cvReleaseImage(&tG);
	cvReleaseImage(&tA);

	cvReleaseImage(&sF);
	cvReleaseImage(&tF);

	cvReleaseImage(&asum);
	cvReleaseImage(&tmpChannel);
}

void GeoTiffMerger::blendGrayImage(IplImage* source, IplImage* target)
{
	for (int y = 0; y < source->height; y++)
	{
		for (int x = 0; x < source->width; x++)
		{
			float* sval = (float*)(source->imageData + y * source->widthStep + x * (_depth / 8));

			IplROI* roi = target->roi;

			float* tval = (float*)(target->imageData +  ( roi->yOffset + y ) * target->widthStep + ( roi->xOffset + x ) * (_depth / 8));

			if (*sval < 0.0f)
			{
				if (*tval <= 0.0f) *tval = -9999;

				continue;
			}

			if (*tval <= 0.0f)
				*tval = *sval;
			else
				*tval = (*tval + *sval) / 2.0f;

			if (*tval < _minVal) _minVal = *tval;
			if (*tval > _maxVal) _maxVal = *tval;			

			if (x == source->width - 1)
			{
				if (x + roi->xOffset + 1 < target->width)
					*(tval + 1) = *(tval);
			}

			if (y == source->height - 1)
			{
				if (y + roi->yOffset + 1 < target->height)
				{
					float* vval = (float*)(target->imageData + (roi->yOffset + y + 1) * target->widthStep + (roi->xOffset + x) * (_depth / 8));
					*vval = *tval;
				}

			}

			if (*tval > 85.0f)
				*tval = *tval - 85.0f;
		}
	}
}

void GeoTiffMerger::mergeTiffFiles()
{
	estimateParameters();
	mergeFiles();
}

void GeoTiffMerger::exportTxtDEM(IplImage* image)
{
	FILE* file = fopen("F:/HBH_DEM_TXT.txt", "w");

	fprintf(file, "ncols         %d\n", image->width);
	fprintf(file, "nrows         %d\n", image->height);

	double sx, sy, sz;
	double mx, my, mz;

	_targetImage->getGeoGridCoord(mx, my, mz);
	_targetImage->getGeoScale(sx, sy, sz);

	fprintf(file, "xllcorner     %3.12f\r\n", mx);
	fprintf(file, "yllcorner     %3.12f\r\n", my);

	fprintf(file, "cellsize       %3.12f\r\n", sx);
	fprintf(file, "NODATA_value  -9999\n");

	for (int y = 0; y < image->height; y++)
	{
		for (int x = 0; x < image->width; x++)
		{
			float* fval = (float*)(image->imageData + image->widthStep * y + x * sizeof(float));

			int val = -9999;
			if (*fval > 0.0f)
				val = *fval + 0.5;

			if (x < image->width - 1)
				fprintf(file, "%d ", val);
			else
				fprintf(file, "%d\n", val);
		}
	}

	fclose(file);
}