#include "stdafx.h"
#include "FltDensityAnalyzer.h"
#include <sstream>
#include <fstream>
#include <numeric>
CFltDensityAnalyzer::CFltDensityAnalyzer()
	:m_fDx(1.0)
    ,m_fDy(1.0)
	,m_aveDensity(0.0)
	, m_minFltLength(0.0)
	, m_bCalcWithRealArea(TRUE)
	, m_fileLastOpen(_T(""))
{
}


CFltDensityAnalyzer::~CFltDensityAnalyzer()
{
	ClearAll();
}
void CFltDensityAnalyzer::ClearAll(void)
{
	ClearFlts();
	ClearPolygons();
}
void CFltDensityAnalyzer::ClearPolygons(void)
{
	for (int i = 0; i < m_cellPolygons.size(); i++)
		delete m_cellPolygons[i];
	m_cellPolygons.clear();
}
void CFltDensityAnalyzer::ClearFlts(void)
{
	FLTITER i = m_flts.begin();

	//for (; i != m_flts.end(); i ++)
	//	delete (*i);
	m_flts.clear();
	m_fltRects.clear();
}
void CFltDensityAnalyzer::SetOutputPaths(/*CString strDenAve,*/
	CString strDenDist, CString strCellPolygons)
{
	//m_strFltInput = strFlt;
	//m_strBorder = strBd;
	//m_strOutputDensityAverage = strDenAve;
	m_strOuputDensityDistribution = strDenDist;
	m_strOutputCellPolygons = strCellPolygons;


}

bool CFltDensityAnalyzer::ReadFltFile(CString strFlt)
{
	m_strFltInput = strFlt;
	m_fileLastOpen = strFlt;
	m_fltMap.ClearAll();
	if (0 == m_fltMap.ReadFile(m_strFltInput.GetBuffer()))
		return false;

	return true;
}
bool CFltDensityAnalyzer::ReadBorderFile(CString strBd)
{
	m_strBorder = strBd;
	return ReadBorder();
}

void CFltDensityAnalyzer::SetGrid(float dx, float dy)
{
	if (dx < 1e-6)
		dx = 1.0;
	if (dy < 1e-6)
		dy = 1.0;
	m_fDx = dx;
	m_fDy = dy;
}

//写出平均密度
bool CFltDensityAnalyzer::WriteAverageDensity(double dden, CString strAveDen)
{
	FILE* fw = fopen((LPSTR)(LPCTSTR)strAveDen, "w");
	if (0 == fw)
		return false;
	fprintf(fw, "断层平均密度: %.8f 条/平方公里\n", dden);
	fclose(fw);
	return true;

}
double CFltDensityAnalyzer::CalcDensity()
{
	RemoveInvalidFlts();
	CreateGridRects(); //生成网格单元矩形
	CalcAverageDensity();
	return m_aveDensity;
}

bool CFltDensityAnalyzer::Execute(void)
{
	RemoveInvalidFlts();
	CreateGridRects(); //生成网格单元矩形
	//平均密度
	/*double aveDensity = CalcAverageDensity();*/
	CalcAverageDensity();
	WriteAverageDensity(m_aveDensity, m_strOutputDensityAverage);
	
	WriteGridFltsInfo();

	WriteCellPolygons();

	return true;
}

void   CFltDensityAnalyzer::WriteGridFltsInfo(void)
{
	std::ifstream file((LPSTR)(LPCTSTR)m_strOuputDensityDistribution);
	if (file.good())
	{
		// 文件存在，关闭文件并删除文件
		file.close();
		std::remove(m_strOuputDensityDistribution);
	}
	FILE* fw = fopen((LPSTR)(LPCTSTR)m_strOuputDensityDistribution, "w");
	if (0 == fw)
		return ;
	//fprintf(fw, "断层平均密度: %.6f 条/平方公里\n", dden);
	GRect8 grt;
	vector<GPline*> gflts;
	int n = 0;
	GPline testline;
	GPoint3D cpt;
	fprintf(fw, "中心点X,中心点Y,条数,面积(Km^2),密度(条/Km^2)\n");
	double area = m_fDy*m_fDx / 1e6;
	for (int i = 0; i < m_cellPolygons.size(); i++)
	{
		grt = m_grdRects[i];
		//gflts.clear();
		//n = GetRectFlts(grt, gflts);		
		////grt.GetCenter(cpt.x0,cpt.y0);
		//cpt = m_grdCenters[i];
		//fprintf(fw, "%.12g,%.12g,%d\n", cpt.x0, cpt.y0, n);
		if (m_bCalcWithRealArea)
		{
			cpt = m_grdCenters[i];
			n = m_grdFltNums[i];
			 area = m_cellPolygons[i]->Area() / 1e6;
			fprintf(fw, "%.12g,%.12g,%d,%.6g,%.6g\n", cpt.x0, cpt.y0, n, area,
				double(n) / area);
		}
		else
		{
			m_grdRects[i].GetCenter(cpt.x0,cpt.y0);
			n = m_grdFltNums[i];
			
			fprintf(fw, "%.12g,%.12g,%d,%.6g,%.6g\n", cpt.x0, cpt.y0, n, area,
				double(n) / area);
		}
	}	
	fclose(fw);
	
}

void CFltDensityAnalyzer::WriteCellPolygons(void)
{
	std::ifstream file((LPSTR)(LPCTSTR)m_strOutputCellPolygons);
	if (file.good())
	{
		// 文件存在，关闭文件并删除文件
		file.close();
		std::remove(m_strOutputCellPolygons);
	}
	FILE* fw = fopen((LPSTR)(LPCTSTR)m_strOutputCellPolygons, "w");
	if (0 == fw)
		return;
	
	CString str; 
	int n = 0;
	GPline* pc = 0;
	GPoint3D cpt;
	double dFontHeight = m_fDx / 5.0;
	double dFontWidth = dFontHeight*0.4;
	fprintf(fw, "HowToViewPoint\n");
	fprintf(fw, "Font 0 13 0 0 0 0 0 1 0 0 0 0 Times New Roman\n");
	//fprintf(fw, "Text 200.0000 80.0000 0.0000 0.0000 0.0000 82\n");
	fprintf(fw, "Text %lf %lf 0.0000 0.0000 0.0000 82\n", dFontHeight, dFontWidth);
	fprintf(fw, "Color 0 0\n");
	fprintf(fw, "FontZ 0 13 0 0 0 0 0 1 0 0 0 0 Times New Roman\n");
	fprintf(fw, "TextZ 0.0000 0.0000 0.0000 0.0000 0.0000 100\n");
	fprintf(fw, "ColorZ 0 0\n");
	fprintf(fw, "Mark 0.0000 0.0000 0.0000\n");
	fprintf(fw, "Mode 2048\n");
	fprintf(fw, "\n");
	fprintf(fw, "Layer M 0\n");
	fprintf(fw, "Layer HowToViewPoint\n");
	fprintf(fw, "HowToViewCurve\n");
	fprintf(fw, "\n");
	fprintf(fw, "HowToViewPoint\n");
	fprintf(fw, "\n");
	fprintf(fw, "State 10\n");
	fprintf(fw, "Color 0\n");
	fprintf(fw, "Layer M 0\n");

	if (m_bCalcWithRealArea)
	{
		for (int i = 0; i < m_cellPolygons.size(); i++)
		{
			n = m_grdFltNums[i];
			str.Format("%d", n);
			pc = m_cellPolygons[i];
			pc->SetName(str.GetBuffer());
			cpt = m_grdCenters[i];
			cpt.SetName(str.GetBuffer());
			pc->Write(fw);
			cpt.Write(fw);


			fprintf(fw, "\n");

		}
	}
	else
	{
		GRect8 grt;
		GPline pc;
		for (int i = 0; i < m_grdRects.size(); i++)
		{
			grt = m_grdRects[i];
			n = m_grdFltNums[i];
			str.Format("%d", n);
			pc.FromRect(grt);


			pc.SetName(str.GetBuffer());
			grt.GetCenter(cpt.x0, cpt.y0);
			cpt.SetName(str.GetBuffer());
			pc.Write(fw);
			cpt.Write(fw);


			fprintf(fw, "\n");

		}
	}

	



	fclose(fw);
}

float CFltDensityAnalyzer::GetFltLength(GPline * line)
{
	double len = line->GetLength();
	if (line->IsClosed(1e-6))
		len /= 2.0;


	return len;


}

int CFltDensityAnalyzer::GetLayers(vector<CString>& strlayers)
{
	
	vector<string> veclayers;
	m_fltMap.GetLayers(veclayers);
	for (int i = 0; i < veclayers.size(); i++)
	{
		strlayers.push_back(veclayers[i].c_str());
	}

	return veclayers.size();


}

bool CFltDensityAnalyzer::CalcAverageDensity(void) //返回平均密度
{
	int N = m_grdFltNums.size();
	if (m_grdFltNums.empty())
		return false;
	double sum = 0.0;
	for (int i = 0; i < m_cellPolygons.size(); i++)
	{
		sum += m_grdFltNums[i];
		//sum += m_cellPolygons[i]->Area() * m_grdFltNums[i]/1e6;
	}
	double s = 0;
	if (m_bCalcWithRealArea)
		s = m_border.Area() / 1e6;
	else
		s = m_grdRects.size() * m_fDx*m_fDy / 1e6;
	double d = sum / s;

	
	m_aveDensity = d;
	return true;
	/*double s = m_border.Area();
	if (s < 1e-8)
		return 0;
	double d = m_flts.size()*1e6/ s;

	return d;*/
}

int CFltDensityAnalyzer::GetRectFlts(const GRect8& rect, vector<GPline*>& dstflts)
{
	
	GPline rcur; rcur.FromRect(rect);
	GRect8 fltrect;
	FLTITER iter = m_flts.begin();
	GPline* p1, *p2;
	GPline* pFlt = 0;
	int count = 0;
	for (; iter != m_flts.end(); iter++)
	{
		p1 = &rcur;
		p2 = *iter;
		pFlt = p2;
		fltrect = m_fltRects[p2];
		if (!rect.IsIntersected(fltrect))
			continue;
		if (p2->IsClosed(1e-6) && p1->Area() < p2->Area())
			std::swap(p1, p2);
		
		if (0 == p1->IsInside(*p2))
			continue;

		dstflts.push_back(pFlt);
		count++;
	}
	return count;
}

int CFltDensityAnalyzer::CreateGridRects(void)
{
	m_grdRects.clear();
	ClearPolygons();
	m_grdCenters.clear();
	m_grdFltNums.clear();
	GRect8 bdrect = m_border.GetRect();
	GRect8 rt;
	int xnum = ceil(bdrect.Width() / m_fDx);
	int ynum = ceil(bdrect.Height() / m_fDy);

	double xx0 = bdrect.left;
	double yy0 = bdrect.bottom;
	double x0, y0;
	int flag = 0;
	double cx, cy;
	GPline* pgnInside = 0;
	GPline* pgnRect = 0;
	m_grdRects.reserve(xnum*ynum);
	vector<GPline*> relFlts;
	for (int j = 0; j < ynum; j++)
	{
		y0 = yy0 + j*m_fDy;
		for (int i = 0; i < xnum; i++)
		{
			x0 = xx0 + i * m_fDx;
			rt.Create(x0, y0, m_fDx, m_fDy);
			flag = m_border.IsInside(rt);
			if(0 == flag) //不相交
				continue;
			


			if (1 == flag) //在边界内部
			{

				rt.GetCenter(cx, cy);
				pgnRect = new GPline();
				pgnRect->FromRect(rt);
				
				//pgnRect.FromRect(rt);
				m_grdRects.push_back(rt);
				m_grdCenters.push_back(GPoint3D(cx, cy, 0));
				m_cellPolygons.push_back(pgnRect);
			}
			else
			{
				pgnRect = new GPline();
				pgnInside = new GPline();
				pgnRect->FromRect(rt);
				if (!m_border.GetIntersectedPolygon(*pgnRect, *pgnInside))
				{
					delete pgnRect;
					delete pgnInside;
					continue;
				}
				pgnInside->GetGravityCenter(cx, cy);


				////begin test
				//if(0 == pgnInside.IsInside(cx,cy))
				//{
				//	FILE* fw = fopen("f:/rectandpgn.dfd", "w");
				//	pgnRect.Write(fw);
				//	pgnInside.Write(fw);
				//	GPoint3D ppt(cx, cy,0);
				//	ppt.Write(fw);
				//	fclose(fw);
				//	printf("");
				//}

				////end test

				m_grdRects.push_back(rt);
				m_grdCenters.push_back(GPoint3D(cx, cy, 0));

				m_cellPolygons.push_back(pgnInside);
				delete pgnRect;

			}
			relFlts.clear();
			int nflt = GetRectFlts(rt, relFlts);
			m_grdFltNums.push_back(nflt);

			//if (2 == flag) //相交 ,判断中心点是否在边界内
			//{
			//	rt.GetCenter(cx, cy);
			//	if (0 == m_border.IsInside(cx, cy))
			//		continue; //中心点不在边界内，跳过
			//}
			//m_grdRects.push_back(rt);
		}
	}


	return m_grdRects.size();
}

//bool CFltDensityAnalyzer::ReadFltFile(void)
//{
//
//	
//
//	////obsolete code
//
//	//ClearFlts();
//
//	//string line;
//	//stringstream sstr;
//
//	//ifstream ifs(m_strFltInput, fstream::in);
//	//if (!ifs)
//	//	return false;
//	//int idx = 0;
//	//GRect8 rect;
//	//int i = 0;
//	//string str;
//	//while (getline(ifs, line))
//	//{
//	//	if (line.size() < 3)
//	//		continue;
//	//	idx = line.find("Pline");
//	//	if (idx < 0)
//	//		continue;
//	//	GPline* pc = new GPline();
//	//	i = line.find_first_of('.', 5);
//	//	
//	//	pc->SetName(line.substr(i + 1).c_str());
//
//	//	if (pc->ReadFile(ifs) < 2 || GetFltLength(pc) < m_minFltLength)
//	//	{
//
//	//		delete pc;
//	//	}
//	//	else
//	//	{
//
//	//		m_flts.push_back(pc);
//	//		rect = pc->GetRect();
//	//		m_fltRects[pc] = rect;
//	//	}
//	//}
//	//return true;
//
//}

bool CFltDensityAnalyzer::ReadFlts(CString strL, bool bIncludingSubLayers /*= true*/)
{
	ClearFlts();
	std::list<GBaseObj*> sels;
	string strLayer(strL.GetBuffer());
	if (strL.IsEmpty())
		sels = m_fltMap.GetElements();
	else
		m_fltMap.GetElements(sels, strLayer, bIncludingSubLayers);

	if (sels.empty())
		return false;
	GPline* pc = 0;
	std::list<GBaseObj*>::iterator iter = sels.begin();
	for (; iter != sels.end(); iter++)
	{
		if (eLine != (*iter)->GetType())
			continue;
		pc = (GPline*)(*iter);


		if (GetFltLength(pc) < m_minFltLength)
			continue;


		m_flts.push_back(pc);
		GRect8 rect = pc->GetRect();
		m_fltRects[pc] = rect;

	}

	return true;
}

//读取指定层位边界
bool CFltDensityAnalyzer::ReadBorder(CString strL, bool bIncludingSubLayer)
{

	m_border.Clear();

	if (strL.IsEmpty())
		return false;
	std::list<GBaseObj*> sels;
	string strLayer(strL.GetBuffer());



	m_fltMap.GetElements(sels, strLayer, bIncludingSubLayer);

	GPline* pc = 0;
	std::list<GBaseObj*>::iterator iter = sels.begin();
	for (; iter != sels.end(); iter++)
	{
		if (eLine != (*iter)->GetType())
			continue;
		m_border = *(GPline*)(*iter);


		GRect8 trt = m_border.GetRect();
		m_fltRects[&m_border] = trt;	

		break;

	}

	return true;



}


bool CFltDensityAnalyzer::ReadBorder(void)
{
	string line;
	std::ifstream ifs(m_strBorder, fstream::in);
	if (!ifs)
		return false;
	int idx = 0;
	GRect8 trt;
	while (getline(ifs, line))
	{
		if (line.size() < 3)
			continue;
		idx = line.find("Pline");
		if (idx < 0)
			continue;
		
		if (m_border.ReadFile(ifs) < 3)
			continue;
		 trt = m_border.GetRect();
		m_fltRects[&m_border] = trt;


		return true;
	}


	return true;
}

int CFltDensityAnalyzer::IsInside(GPline & pc1, GPline & pc2)
{
	GRect8 rt1, rt2;
	rt1 = m_fltRects[&pc1];
	rt2 = m_fltRects[&pc2];
	//首先快速排斥
	if (!rt1.IsIntersected(rt2))
		return false;


	return pc1.IsInside(pc2);
}

void CFltDensityAnalyzer::RemoveInvalidFlts(void)
{
	std::list<GPline*>::iterator iter = m_flts.begin();
	vector<GPline*> invalidplines;
	for (iter = m_flts.begin(); iter != m_flts.end(); iter++)
	{
		if (0 == m_border.IsInside(*(*iter)))
			invalidplines.push_back(*iter);
	}
	
	for (int i = 0; i < invalidplines.size(); i++)
	{
		m_flts.remove(invalidplines[i]);
		//delete invalidplines[i];
	}


}