/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#include "trackline.h"

double node2dDistance(node2d n1,node2d n2)
{
	return sqrt(pow(n2.x-n1.x,2)+pow(n2.y-n1.y,2));
}

int trackLine::getInputNode(char* para)
{
	char order[1024] = "0,1,2";
	char filename[1024] = "NULL";
	int orders[3];
	double temp_d;
	std::vector<double> tempRow;
	string temp_str;
	stringstream temp_ss;

	//首先尝试分离文件名和+d命令
	if (2 != sscanf(para,"%[^+]+d%s",filename,order))
	{
		//分离失败 直接将para解释为文件名
		strcpy(filename,para);
	}

	ifstream datain;
	gctl::open_infile(datain,filename);

	if (3 != sscanf(order,"%d,%d,%d",&orders[0],&orders[1],&orders[2]))
	{
		cout << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "wrong order parameters for the file: " << filename << endl;
		return -1;
	}

	node2d tempNode;
	while(getline(datain,temp_str))
	{
		if(*(temp_str.begin()) == '#') continue;
		temp_ss.clear();
		temp_ss << temp_str;
		//解析数据行
		if(!tempRow.empty()) tempRow.clear();
		while (temp_ss >> temp_d)
		{
			tempRow.push_back(temp_d);
		}

		//读入指定位置的数据
		tempNode.x = tempRow[orders[0]];
		tempNode.y = tempRow[orders[1]];
		tempNode.val = tempRow[orders[2]];
		tempNode.id = inputNode.size();
		inputNode.push_back(tempNode);
	}
	datain.close();
	return 0;
}

int trackLine::initBox(char* para)
{
	//定位二维点位的中心位置
	xmin = ymin = GCTL_BDL_MAX;
	xmax = ymax = GCTL_BDL_MIN;
	for (int i = 0; i < inputNode.size(); i++)
	{
		if (inputNode[i].x < xmin) xmin = inputNode[i].x;
		if (inputNode[i].x > xmax) xmax = inputNode[i].x;
		if (inputNode[i].y < ymin) ymin = inputNode[i].y;
		if (inputNode[i].y > ymax) ymax = inputNode[i].y;
	}

	double cx,cy;
	cx = 0.5*(xmin+xmax); cy = 0.5*(ymin+ymax);
	//确定盒子的边长
	double pointArea;
	/*这里我们发现使用任意两点间的最短距离作为自动网格化大小并不能取得很好的效果
	经常可能会造成插值点无值的状况 因此我们这里将使用平均点位面积估算网格大小 效果出众 而且这样做仅需要很少的运算*/
	if (2 != sscanf(para,"%lf/%lf",&dx,&dy))
	{
		pointArea = (xmax-xmin)*(ymax-ymin)/inputNode.size();
		dx = dy = sqrt(pointArea);
	}
	//确定盒子在x和y方向的半盒子数 向上取整
	int halfNum_x = ceil((xmax-cx)/dx);
	int halfNum_y = ceil((ymax-cy)/dy);
	//重新确定盒子的范围 向外推出一层
	xmin = cx - halfNum_x*dx - dx;
	xmax = cx + halfNum_x*dx + dx;
	ymin = cy - halfNum_y*dy - dy;
	ymax = cy + halfNum_y*dy + dy;
	//确定盒子在x和y方向的盒子数
	xnum = int ((xmax-xmin)/dx);
	ynum = int ((ymax-ymin)/dy);
	//初始化盒子上的节点
	node2d tempNode;
	for (int j = 0; j < ynum+1; j++)
	{
		for (int k = 0; k < xnum+1; k++)
		{
			tempNode.x = xmin + k*dx;
			tempNode.y = ymin + j*dy;
			tempNode.val = GCTL_BDL_MAX;
			tempNode.id = boxNode.size();
			boxNode.push_back(tempNode);
		}
	}
	//将输入点位放置在盒子中
	int temp_m,temp_n;
	boxIndex.resize(xnum*ynum);
	for (int i = 0; i < inputNode.size(); i++)
	{
		temp_m = floor((inputNode[i].x - xmin)/dx);
		temp_n = floor((inputNode[i].y - ymin)/dy);
		boxIndex[temp_n*xnum+temp_m].push_back(inputNode[i].id);
	}
	//查看每个盒子节点周围四个盒子中的输入点 插值确定盒子节点上的值
	double tempRes,oneDist,totalDist;
	int rectIndex;
	int nodeIndex;
	for (int j = 0; j < ynum-1; j++)
	{
		for (int k = 0; k < xnum-1; k++)
		{
			tempRes = 0;
			totalDist = 0;
			nodeIndex = (j+1)*(xnum+1)+k+1;
			//查看四个格子内的顶点 并利用距离平方反比插值
			for (int n = 0; n < 2; n++)
			{
				for (int p = 0; p < 2; p++)
				{
					rectIndex = (j+n)*xnum+k+p;
					if (!boxIndex[rectIndex].empty())
					{
						for (int b = 0; b < boxIndex[rectIndex].size(); b++)
						{
							totalDist += 
							1.0/(1e-30
								+pow(inputNode[boxIndex[rectIndex][b]].x - boxNode[nodeIndex].x,2)
								+pow(inputNode[boxIndex[rectIndex][b]].y - boxNode[nodeIndex].y,2));
						}
					}
				}
			}
			if (totalDist == 0)
			{
				boxNode[nodeIndex].val = GCTL_BDL_MAX;
				continue;
			}
			//计算插值
			for (int n = 0; n < 2; n++)
			{
				for (int p = 0; p < 2; p++)
				{
					rectIndex = (j+n)*xnum+k+p;
					if (!boxIndex[rectIndex].empty())
					{
						for (int b = 0; b < boxIndex[rectIndex].size(); b++)
						{
							oneDist = 
							1.0/(1e-30
								+pow(inputNode[boxIndex[rectIndex][b]].x - boxNode[nodeIndex].x,2)
								+pow(inputNode[boxIndex[rectIndex][b]].y - boxNode[nodeIndex].y,2));
							tempRes += inputNode[boxIndex[rectIndex][b]].val*oneDist/totalDist;
						}
					}
				}
			}
			boxNode[nodeIndex].val = tempRes;
		}
	}
	return 0;
}

int trackLine::initLine(char* para)
{
	char order[1024] = "0,1";
	char filename[1024] = "NULL";
	int orders[2];
	double temp_d;
	std::vector<double> tempRow;
	string temp_str;
	stringstream temp_ss;

	node2d tempNode;
	double length;
	double xinterval;
	node2d n1,n2,snode,enode;
	//首先尝试从参数初始化待插值点
	if (5 == sscanf(para,"%lf/%lf/%lf/%lf/%lf",&n1.x,&n1.y,&n2.x,&n2.y,&xinterval))
	{
		snode = n1; enode = n2;
		//计算水平方向剖面长度
		double plen = sqrt(pow(enode.x-snode.x,2)+pow(enode.y-snode.y,2));

		length = 0;
		while(length <= plen)
		{
			tempNode.x = snode.x+(enode.x-snode.x)*length/plen;
			tempNode.y = snode.y+(enode.y-snode.y)*length/plen;
			tempNode.px = length;
			tempNode.id = lineNode.size();
			lineNode.push_back(tempNode);
			length += xinterval;
		}
	}
	//首先尝试分离文件名和+d命令
	else if (2 == sscanf(para,"%[^+]+d%s",filename,order))
	{
		ifstream datain;
		gctl::open_infile(datain,filename);

		if (2 != sscanf(order,"%d,%d",&orders[0],&orders[1]))
		{
			cout << GCTL_BOLDRED << "error ==> " << GCTL_RESET << "wrong order parameters for the file: " << filename << endl;
			return -1;
		}

		while(getline(datain,temp_str))
		{
			if(*(temp_str.begin()) == '#') continue;
			temp_ss.clear();
			temp_ss << temp_str;
			//解析数据行
			if(!tempRow.empty()) tempRow.clear();
			while (temp_ss >> temp_d)
			{
				tempRow.push_back(temp_d);
			}

			//读入指定位置的数据
			tempNode.x = tempRow[orders[0]];
			tempNode.y = tempRow[orders[1]];
			tempNode.id = lineNode.size();
			lineNode.push_back(tempNode);
		}
		datain.close();

		//计算待插值点在剖面上的累积长度
		length = 0;
		lineNode[0].px = 0.0;
		for (int i = 1; i < lineNode.size(); i++)
		{
			length += node2dDistance(lineNode[i-1],lineNode[i]);
			lineNode[i].px = length;
		}
	}
	//分离失败 直接将para解释为文件名
	else
	{
		strcpy(filename,para);

		ifstream datain;
		gctl::open_infile(datain,filename);

		while(getline(datain,temp_str))
		{
			if(*(temp_str.begin()) == '#') continue;
			temp_ss.clear();
			temp_ss << temp_str;
			//解析数据行
			if(!tempRow.empty()) tempRow.clear();
			while (temp_ss >> temp_d)
			{
				tempRow.push_back(temp_d);
			}

			//读入指定位置的数据
			tempNode.x = tempRow[orders[0]];
			tempNode.y = tempRow[orders[1]];
			tempNode.id = lineNode.size();
			lineNode.push_back(tempNode);
		}
		datain.close();

		//计算待插值点在剖面上的累积长度
		length = 0;
		lineNode[0].px = 0.0;
		for (int i = 1; i < lineNode.size(); i++)
		{
			length += node2dDistance(lineNode[i-1],lineNode[i]);
			lineNode[i].px = length;
		}
	}
	return 0;
}

int trackLine::interLine()
{
	int temp_m,temp_n;
	for (int i = 0; i < lineNode.size(); i++)
	{
		//检查待插值点的位置是否在盒子内
		if (lineNode[i].x < xmin + dx || 
			lineNode[i].x > xmax - dx || 
			lineNode[i].y < ymin + dy || 
			lineNode[i].y > ymax - dy)
		{
			lineNode[i].val = GCTL_BDL_MAX;
			continue;
		}
		temp_m = floor((lineNode[i].x - xmin)/dx);
		temp_n = floor((lineNode[i].y - ymin)/dy);
		lineNode[i].val = gctl::rect_interpolate(xmin+temp_m*dx,ymin+temp_n*dy,dx,dy,
											lineNode[i].x,lineNode[i].y,
											boxNode[temp_n*(xnum+1)+temp_m].val,
											boxNode[temp_n*(xnum+1)+temp_m+1].val,
											boxNode[(temp_n+1)*(xnum+1)+temp_m+1].val,
											boxNode[(temp_n+1)*(xnum+1)+temp_m].val);
	}
	return 0;
}

int trackLine::outLine(char* filename)
{
	time_t now = time(0);
	char* dt = ctime(&now);
	ofstream lineOut;
	gctl::open_outfile(lineOut,filename);
	lineOut << "# This file is generated by trackLine on " << dt;
	lineOut << "# For more information please contact the author via: zhangyi.cugwuhan@gmail.com" << endl;
	lineOut << "# XonLine XonPlane YonPlane Val" << endl;
	for (int i = 0; i < lineNode.size(); i++)
	{
		if (lineNode[i].val == GCTL_BDL_MAX)
		{
			lineOut << lineNode[i].px << " " << lineNode[i].x << " " << lineNode[i].y << " NaN" << endl;
		}
		else lineOut << lineNode[i].px << " " << lineNode[i].x << " " << lineNode[i].y << " " << lineNode[i].val << endl;
	}
	lineOut.close();
	return 0;
}