// MeshDecEKavg.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <array>
#include <string>
#include <algorithm>
//#include <chrono>
//#include <windows.h>

using namespace std;
//using namespace std::chrono;

typedef array<int, 3> Array3i;
typedef array<int, 2> Array2i;
typedef array<double, 3> Array3D;
typedef vector<array<int, 3>> VectArr3i;
typedef vector<array<double, 3>> VectArr3d;
typedef map<int, vector<int>> Map1iVecti;
typedef array<array<double, 4>, 4> QMatrix;
typedef vector< pair <double, int> > VectPair;


struct CVertex
{
	double X;
	double Y;
	double Z;
	QMatrix Q;
};
typedef vector<CVertex> VectorVertex;
typedef map<int, CVertex > MapVertex;
struct CEdge
{
	Array2i IndexOfFace;
	Array2i Verts;
	CVertex Vc;
	double Error;
};
typedef map<int, CEdge> MapEdges;
class CTriangle
{
public:
	Array3i Edges;
};
typedef map<int, CTriangle> MapTriangle;
class CMesh
{
public:
	MapEdges EdgesWithInfo;
	MapTriangle CTriangles;
	Map1iVecti triAdjToPoint, edgAdjToPoint;

	CMesh(VectArr3i& triangles)
	{
		int i = 0;
		int j = 0;
		map< Array2i, int> mapOfEdges;
		map< Array2i, int>::iterator imaoe1, imaoe2;
		Map1iVecti::iterator ieop;

		for (VectArr3i::iterator it = triangles.begin(); it != triangles.end(); ++it)
		{
			CTriangle t;
			for (int k = 0; k < 3; k++)
			{
				ieop = triAdjToPoint.find((*it)[k]);
				if (ieop == triAdjToPoint.end())
				{
					triAdjToPoint[(*it)[k]] = { i };
				}
				else
				{
					(*ieop).second.push_back(i);
				}

				imaoe1 = mapOfEdges.find({ (*it)[k], (*it)[(k + 1) % 3] });
				if (imaoe1 == mapOfEdges.end())
				{
					imaoe2 = mapOfEdges.find({ (*it)[(k + 1) % 3],(*it)[k] });
					if (imaoe2 == mapOfEdges.end())
					{
						mapOfEdges[{ (*it)[k], (*it)[(k + 1) % 3] }] = j;
						CEdge f;
						f.Verts = { (*it)[k], (*it)[(k + 1) % 3] };
						f.IndexOfFace[0] = i;
						EdgesWithInfo[j] = f;
						edgAdjToPoint[(*it)[k]].push_back(j);
						edgAdjToPoint[(*it)[(k + 1) % 3]].push_back(j);
						t.Edges[k] = j;
						j++;
					}
					else
					{
						EdgesWithInfo[(*imaoe2).second].IndexOfFace[1] = i;
						t.Edges[k] = (*imaoe2).second;
					}
				}
				else
				{
					EdgesWithInfo[(*imaoe1).second].IndexOfFace[1] = i;
					t.Edges[k] = (*imaoe1).second;
				}
			}
			CTriangles[i] = t;
			i++;
		}
	};
};
bool isAscii(ifstream& sourceFile)
{
	char a[4];
	sourceFile.seekg(80);
	sourceFile.read(a, 4);

	unsigned int numberOfTriangles = *reinterpret_cast<unsigned int*>(a);

	int fileSizeCalculated = 84 + (numberOfTriangles * 50);

	sourceFile.seekg(0);
	streampos beginOfFile, endOfFile;
	beginOfFile = sourceFile.tellg();
	sourceFile.seekg(0, ios::end);
	endOfFile = sourceFile.tellg();
	int fileSize = endOfFile - beginOfFile;

	if (fileSizeCalculated == 0)
	{
		return false;
	}
	if (fileSize == fileSizeCalculated)
	{
		return false;
	}
	sourceFile.seekg(0);
	return true;
}
CVertex getVertexValues(string lineOfFile, int positionInLine)
{
	string numberFromFile;
	double value;
	CVertex vertex;

	lineOfFile = lineOfFile + ' ';
	positionInLine = positionInLine + 7;
	for (int i = 0; i < 3; i++)
	{
		while (lineOfFile.at(positionInLine) != ' ')
		{
			numberFromFile = numberFromFile + lineOfFile.at(positionInLine);
			positionInLine = positionInLine + 1;
		}
		value = stod(numberFromFile);

		if (value< FLT_EPSILON && value > -FLT_EPSILON)
		{
			value = 0.0;
		}

		switch (i)
		{
		case 0: vertex.X = value;
			break;
		case 1: vertex.Y = value;
			break;
		case 2: vertex.Z = value;
			break;
		}

		numberFromFile.clear();
		positionInLine = positionInLine + 1;
	}
	return vertex;
}
void asciiFile(ifstream& sourceFile, VectorVertex& vertexes)
{
	int positionOfStrVertex;
	string line;
	while (getline(sourceFile, line))
	{
		positionOfStrVertex = line.find("vertex");
		if (positionOfStrVertex != -1)
		{
			vertexes.push_back(getVertexValues(line, positionOfStrVertex));
		}
	};
	sourceFile.seekg(0);
}
void binaryFile(ifstream& sourceFile, VectorVertex& vertexes)
{
	char a[4];
	CVertex vertex;
	sourceFile.seekg(80);
	sourceFile.read(a, 4);
	unsigned int numberOfTriangles = *reinterpret_cast<unsigned int*>(a);
	for (int i = 0; i < numberOfTriangles; i++)
	{
		for (int j = 0; j < 3; j++)
		{
			sourceFile.read(a, 4);
		}
		for (int j = 0; j < 3; j++)
		{
			for (int k = 0; k < 3; k++)
			{
				sourceFile.read(a, 4);
				double value = *reinterpret_cast<float*>(a);

				switch (k)
				{
				case 0:vertex.X = value;
					break;
				case 1:vertex.Y = value;
					break;
				case 2:vertex.Z = value;
					break;
				}
			}
			vertexes.push_back(vertex);
		}
		sourceFile.read(a, 2);
	}
	sourceFile.seekg(0);
}
void getUniqVertAndTrianglesList(VectorVertex vertexes, VectArr3i& listOfTriangles, MapVertex& uniqueVertexes)
{
	Array3D v;
	Array3i pos;
	map<Array3D, int> a;
	map<Array3D, int>::iterator iti;
	int i = 0;
	int j = 0;

	for (VectorVertex::iterator it = vertexes.begin(); it != vertexes.end(); ++it, j++)
	{
		v = { (*it).X ,(*it).Y ,(*it).Z };
		iti = a.find(v);

		if (iti == a.end())
		{
			a[v] = i;
			pos[j % 3] = i;
			uniqueVertexes[i] = (*it);
			i++;
		}
		else
		{
			pos[j % 3] = (*iti).second;
		}

		if (j % 3 == 2)
		{
			listOfTriangles.push_back(pos);
		}
	}
}
Array3D normal(CVertex A, CVertex B, CVertex C)
{
	Array3D normal, v1, v2;
	double lenght;

	v1 = { A.X - B.X, A.Y - B.Y, A.Z - B.Z };
	v2 = { C.X - B.X, C.Y - B.Y, C.Z - B.Z };

	normal[0] = v1[1] * v2[2] - v1[2] * v2[1];
	normal[1] = v1[2] * v2[0] - v1[0] * v2[2];
	normal[2] = v1[0] * v2[1] - v1[1] * v2[0];

	lenght = sqrt(normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2]);

	if (lenght >= DBL_EPSILON)
	{
		normal[0] = normal[0] / lenght;
		normal[1] = normal[1] / lenght;
		normal[2] = normal[2] / lenght;
		return normal;
	}

	return{ 0.0,0.0,0.0 };
}
void getVertexMatrices(CMesh& mesh, MapVertex& vertices)
{
	//iterace pres body
	for (Map1iVecti::iterator it = mesh.triAdjToPoint.begin(); it != mesh.triAdjToPoint.end(); ++it)
	{
		QMatrix Q;
		Q[0] = { 0.0,0.0,0.0 };
		Q[1] = { 0.0,0.0,0.0 };
		Q[2] = { 0.0,0.0,0.0 };
		Q[3] = { 0.0,0.0,0.0 };
		//iterace pres prilehle trojuhelniky
		for (vector<int>::iterator i = (*it).second.begin(); i != (*it).second.end(); i++)
		{
			int a = mesh.EdgesWithInfo[mesh.CTriangles[(*i)].Edges[0]].Verts[0];
			int b = mesh.EdgesWithInfo[mesh.CTriangles[(*i)].Edges[0]].Verts[1];
			int c = mesh.EdgesWithInfo[mesh.CTriangles[(*i)].Edges[1]].Verts[0];
			CVertex A = vertices[a];
			CVertex B = vertices[b];
			CVertex C;
			Array3D n;

			if (a == c || b == c)
			{
				C = vertices[mesh.EdgesWithInfo[mesh.CTriangles[(*i)].Edges[1]].Verts[1]];
			}
			else
			{
				C = vertices[c];
			}
			n = normal(A, B, C);
			double d = -n[0] * A.X - n[1] * A.Y - n[2] * A.Z;
			array<double, 4 > plEq = { n[0] ,n[1], n[2],d };
			for (int k = 0; k < 4; k++)
			{
				for (int j = 0; j < 4; j++)
				{
					Q[k][j] = Q[k][j] + plEq[k] * plEq[j];
				}
			}
		}
		vertices[(*it).first].Q = Q;
	}
}
bool flipCheck(int edg, MapEdges& mEdg, Map1iVecti& edgAdjToPoint)
{
	//vsechny hrany
	bool flip = false;
	vector<int>e1 = edgAdjToPoint[mEdg[edg].Verts[0]];
	vector<int>e2 = edgAdjToPoint[mEdg[edg].Verts[1]];
	vector<int>p1, p2;
	for (vector<int>::iterator it = e1.begin(); it != e1.end(); it++)
	{
		if (mEdg[(*it)].Verts[0] != mEdg[edg].Verts[0] && mEdg[(*it)].Verts[0] != mEdg[edg].Verts[1])
		{
			p1.push_back(mEdg[(*it)].Verts[0]);
		}
		else if (mEdg[(*it)].Verts[1] != mEdg[edg].Verts[0] && mEdg[(*it)].Verts[1] != mEdg[edg].Verts[1])
		{
			p1.push_back(mEdg[(*it)].Verts[1]);
		}
	}
	for (vector<int>::iterator it = e2.begin(); it != e2.end(); it++)
	{
		if (mEdg[(*it)].Verts[0] != mEdg[edg].Verts[0] && mEdg[(*it)].Verts[0] != mEdg[edg].Verts[1])
		{
			p2.push_back(mEdg[(*it)].Verts[0]);
		}
		else if (mEdg[(*it)].Verts[1] != mEdg[edg].Verts[0] && mEdg[(*it)].Verts[1] != mEdg[edg].Verts[1])
		{
			p2.push_back(mEdg[(*it)].Verts[1]);
		}
	}
	int k = 0;
	for (vector<int>::iterator it = p1.begin(); it != p1.end(); it++)
	{
		for (vector<int>::iterator it2 = p2.begin(); it2 != p2.end(); it2++)
		{
			if ((*it) == (*it2))
			{
				k++;
			}
		}
	}
	if (k > 2)
	{
		flip = true;
	}
	return flip;
}
double det(QMatrix& m)
{
	double value;
	value =
		m[0][3] * m[1][2] * m[2][1] * m[3][0] - m[0][2] * m[1][3] * m[2][1] * m[3][0] - m[0][3] * m[1][1] * m[2][2] * m[3][0] + m[0][1] * m[1][3] * m[2][2] * m[3][0] +
		m[0][2] * m[1][1] * m[2][3] * m[3][0] - m[0][1] * m[1][2] * m[2][3] * m[3][0] - m[0][3] * m[1][2] * m[2][0] * m[3][1] + m[0][2] * m[1][3] * m[2][0] * m[3][1] +
		m[0][3] * m[1][0] * m[2][2] * m[3][1] - m[0][0] * m[1][3] * m[2][2] * m[3][1] - m[0][2] * m[1][0] * m[2][3] * m[3][1] + m[0][0] * m[1][2] * m[2][3] * m[3][1] +
		m[0][3] * m[1][1] * m[2][0] * m[3][2] - m[0][1] * m[1][3] * m[2][0] * m[3][2] - m[0][3] * m[1][0] * m[2][1] * m[3][2] + m[0][0] * m[1][3] * m[2][1] * m[3][2] +
		m[0][1] * m[1][0] * m[2][3] * m[3][2] - m[0][0] * m[1][1] * m[2][3] * m[3][2] - m[0][2] * m[1][1] * m[2][0] * m[3][3] + m[0][1] * m[1][2] * m[2][0] * m[3][3] +
		m[0][2] * m[1][0] * m[2][1] * m[3][3] - m[0][0] * m[1][2] * m[2][1] * m[3][3] - m[0][1] * m[1][0] * m[2][2] * m[3][3] + m[0][0] * m[1][1] * m[2][2] * m[3][3];
	return value;
}
bool inverse(QMatrix& m, QMatrix& Q)
{
	double deter = det(m);
	if (deter <= FLT_EPSILON)
	{
		return false;
	}
	else
	{
		Q[0][0] = (m[1][2] * m[2][3] * m[3][1] - m[1][3] * m[2][2] * m[3][1] + m[1][3] * m[2][1] * m[3][2] - m[1][1] * m[2][3] * m[3][2] - m[1][2] * m[2][1] * m[3][3] + m[1][1] * m[2][2] * m[3][3]) / deter;
		Q[0][1] = (m[0][3] * m[2][2] * m[3][1] - m[0][2] * m[2][3] * m[3][1] - m[0][3] * m[2][1] * m[3][2] + m[0][1] * m[2][3] * m[3][2] + m[0][2] * m[2][1] * m[3][3] - m[0][1] * m[2][2] * m[3][3]) / deter;
		Q[0][2] = (m[0][2] * m[1][3] * m[3][1] - m[0][3] * m[1][2] * m[3][1] + m[0][3] * m[1][1] * m[3][2] - m[0][1] * m[1][3] * m[3][2] - m[0][2] * m[1][1] * m[3][3] + m[0][1] * m[1][2] * m[3][3]) / deter;
		Q[0][3] = (m[0][3] * m[1][2] * m[2][1] - m[0][2] * m[1][3] * m[2][1] - m[0][3] * m[1][1] * m[2][2] + m[0][1] * m[1][3] * m[2][2] + m[0][2] * m[1][1] * m[2][3] - m[0][1] * m[1][2] * m[2][3]) / deter;
		Q[1][0] = (m[1][3] * m[2][2] * m[3][0] - m[1][2] * m[2][3] * m[3][0] - m[1][3] * m[2][0] * m[3][2] + m[1][0] * m[2][3] * m[3][2] + m[1][2] * m[2][0] * m[3][3] - m[1][0] * m[2][2] * m[3][3]) / deter;
		Q[1][1] = (m[0][2] * m[2][3] * m[3][0] - m[0][3] * m[2][2] * m[3][0] + m[0][3] * m[2][0] * m[3][2] - m[0][0] * m[2][3] * m[3][2] - m[0][2] * m[2][0] * m[3][3] + m[0][0] * m[2][2] * m[3][3]) / deter;
		Q[1][2] = (m[0][3] * m[1][2] * m[3][0] - m[0][2] * m[1][3] * m[3][0] - m[0][3] * m[1][0] * m[3][2] + m[0][0] * m[1][3] * m[3][2] + m[0][2] * m[1][0] * m[3][3] - m[0][0] * m[1][2] * m[3][3]) / deter;
		Q[1][3] = (m[0][2] * m[1][3] * m[2][0] - m[0][3] * m[1][2] * m[2][0] + m[0][3] * m[1][0] * m[2][2] - m[0][0] * m[1][3] * m[2][2] - m[0][2] * m[1][0] * m[2][3] + m[0][0] * m[1][2] * m[2][3]) / deter;
		Q[2][0] = (m[1][1] * m[2][3] * m[3][0] - m[1][3] * m[2][1] * m[3][0] + m[1][3] * m[2][0] * m[3][1] - m[1][0] * m[2][3] * m[3][1] - m[1][1] * m[2][0] * m[3][3] + m[1][0] * m[2][1] * m[3][3]) / deter;
		Q[2][1] = (m[0][3] * m[2][1] * m[3][0] - m[0][1] * m[2][3] * m[3][0] - m[0][3] * m[2][0] * m[3][1] + m[0][0] * m[2][3] * m[3][1] + m[0][1] * m[2][0] * m[3][3] - m[0][0] * m[2][1] * m[3][3]) / deter;
		Q[2][2] = (m[0][1] * m[1][3] * m[3][0] - m[0][3] * m[1][1] * m[3][0] + m[0][3] * m[1][0] * m[3][1] - m[0][0] * m[1][3] * m[3][1] - m[0][1] * m[1][0] * m[3][3] + m[0][0] * m[1][1] * m[3][3]) / deter;
		Q[2][3] = (m[0][3] * m[1][1] * m[2][0] - m[0][1] * m[1][3] * m[2][0] - m[0][3] * m[1][0] * m[2][1] + m[0][0] * m[1][3] * m[2][1] + m[0][1] * m[1][0] * m[2][3] - m[0][0] * m[1][1] * m[2][3]) / deter;
		Q[3][0] = (m[1][2] * m[2][1] * m[3][0] - m[1][1] * m[2][2] * m[3][0] - m[1][2] * m[2][0] * m[3][1] + m[1][0] * m[2][2] * m[3][1] + m[1][1] * m[2][0] * m[3][2] - m[1][0] * m[2][1] * m[3][2]) / deter;
		Q[3][1] = (m[0][1] * m[2][2] * m[3][0] - m[0][2] * m[2][1] * m[3][0] + m[0][2] * m[2][0] * m[3][1] - m[0][0] * m[2][2] * m[3][1] - m[0][1] * m[2][0] * m[3][2] + m[0][0] * m[2][1] * m[3][2]) / deter;
		Q[3][2] = (m[0][2] * m[1][1] * m[3][0] - m[0][1] * m[1][2] * m[3][0] - m[0][2] * m[1][0] * m[3][1] + m[0][0] * m[1][2] * m[3][1] + m[0][1] * m[1][0] * m[3][2] - m[0][0] * m[1][1] * m[3][2]) / deter;
		Q[3][3] = (m[0][1] * m[1][2] * m[2][0] - m[0][2] * m[1][1] * m[2][0] + m[0][2] * m[1][0] * m[2][1] - m[0][0] * m[1][2] * m[2][1] - m[0][1] * m[1][0] * m[2][2] + m[0][0] * m[1][1] * m[2][2]) / deter;
		return true;
	}
}
double errFcn(CVertex Vc, QMatrix Q)
{
	double error;
	error = (Vc.X*Q[0][0] + Vc.Y*Q[1][0] + Vc.Z*Q[2][0] + Q[3][0])*Vc.X + (Vc.X*Q[0][1] + Vc.Y*Q[1][1] + Vc.Z*Q[2][1] + Q[3][1])*Vc.Y + (Vc.X*Q[0][2] + Vc.Y*Q[1][2] + Vc.Z*Q[2][2] + Q[3][2])*Vc.Z + (Vc.X*Q[0][3] + Vc.Y*Q[1][3] + Vc.Z*Q[2][3] + Q[3][3]);
	return error;
}
void getErrors(CMesh& mesh, MapVertex& vertices, VectPair& VP)
{
	for (MapEdges::iterator it = mesh.EdgesWithInfo.begin(); it != mesh.EdgesWithInfo.end(); it++)
	{
		if (!flipCheck((*it).first, mesh.EdgesWithInfo, mesh.edgAdjToPoint))
		{
			QMatrix Qc, Qcd, Qcdi;
			for (int i = 0; i < 4; i++)
			{
				for (int j = 0; j < 4; j++)
				{
					Qc[i][j] = vertices[(*it).second.Verts[0]].Q[i][j] + vertices[(*it).second.Verts[1]].Q[i][j];
					if (i < 3)
					{
						Qcd[i][j] = Qc[i][j];
					}
				}
			}
			Qcd[3] = { 0.0, 0.0, 0.0, 1.0 };
			if (inverse(Qcd, Qcdi))
			{
				(*it).second.Vc.X = Qcdi[0][3];
				(*it).second.Vc.Y = Qcdi[1][3];
				(*it).second.Vc.Z = Qcdi[2][3];
				(*it).second.Vc.Q = Qc;
				(*it).second.Error = errFcn((*it).second.Vc, Qc);

				
				VP.push_back(make_pair((*it).second.Error, (*it).first));
			}
			else
			{
				CVertex avg, A, B;
				A = vertices[(*it).second.Verts[0]];
				B = vertices[(*it).second.Verts[1]];

				avg.X = (A.X + B.X) / 2;	avg.Y = (A.Y + B.Y) / 2;	avg.Z = (A.Z + B.Z) / 2;
				double a = errFcn(A, Qc);
				double b = errFcn(B, Qc);
				double c = errFcn(avg, Qc);
				double m = min(min(a, b), c);

				if (m == a)
				{
					(*it).second.Vc = A;
					(*it).second.Vc.Q = Qc;
					(*it).second.Error = a;
					VP.push_back(make_pair(a, (*it).first));
				}
				else if (m == b)
				{
					(*it).second.Vc = B;
					(*it).second.Vc.Q = Qc;
					(*it).second.Error = b;
					VP.push_back(make_pair(b, (*it).first));
				}
				else
				{
					(*it).second.Vc = avg;
					(*it).second.Vc.Q = Qc;
					(*it).second.Error = c;
					VP.push_back(make_pair(c, (*it).first));
				}

			}
		}
	}
	sort(VP.begin(), VP.end());
}
void createEdge(CMesh& mesh, VectPair VP, int deletTri, Array2i deletVert, CEdge& E1, int indexH, Array2i& edgI)
{

	vector<int> faces, edgInd;
	Array3i::iterator it;
	// pokud hrana trojuhelniku neni mazana
	if (mesh.CTriangles[deletTri].Edges[0] != VP[0].second)
	{

		//pokud bod teto hrany neni mazany
		if (mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[0]].Verts[0] != deletVert[0] && mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[0]].Verts[0] != deletVert[1])
		{
			E1.Verts[1] = mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[0]].Verts[0];
		}
		else
		{
			E1.Verts[1] = mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[0]].Verts[1];
		}

		// pokud prilehly trojuhelnik neni mazan !! pro kazdou hranu

		if (mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[0]].IndexOfFace[0] != deletTri)
		{
			faces.push_back(mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[0]].IndexOfFace[0]);
		}
		else
		{
			faces.push_back(mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[0]].IndexOfFace[1]);
		}
		edgInd.push_back(mesh.CTriangles[deletTri].Edges[0]);

		if (mesh.CTriangles[deletTri].Edges[1] != VP[0].second)
		{

			if (mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[1]].IndexOfFace[0] != deletTri)
			{
				faces.push_back(mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[1]].IndexOfFace[0]);
			}
			else
			{
				faces.push_back(mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[1]].IndexOfFace[1]);
			}
			edgInd.push_back(mesh.CTriangles[deletTri].Edges[1]);
		}
		else
		{

			if (mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[2]].IndexOfFace[0] != deletTri)
			{
				faces.push_back(mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[2]].IndexOfFace[0]);

			}
			else
			{
				faces.push_back(mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[2]].IndexOfFace[1]);

			}
			edgInd.push_back(mesh.CTriangles[deletTri].Edges[2]);
		}
	}
	else
	{
		//pokud bod teto hrany neni mazany
		if (mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[1]].Verts[0] != deletVert[0] && mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[1]].Verts[0] != deletVert[1])
		{
			E1.Verts[1] = mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[1]].Verts[0];
		}
		else
		{
			E1.Verts[1] = mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[1]].Verts[1];
		}

		if (mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[1]].IndexOfFace[0] != deletTri)
		{
			faces.push_back(mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[1]].IndexOfFace[0]);
		}
		else
		{
			faces.push_back(mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[1]].IndexOfFace[1]);
		}

		edgInd.push_back(mesh.CTriangles[deletTri].Edges[1]);
		if (mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[2]].IndexOfFace[0] != deletTri)
		{
			faces.push_back(mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[2]].IndexOfFace[0]);
		}
		else
		{
			faces.push_back(mesh.EdgesWithInfo[mesh.CTriangles[deletTri].Edges[2]].IndexOfFace[1]);
		}
		edgInd.push_back(mesh.CTriangles[deletTri].Edges[2]);
	}


	edgI[0] = edgInd[0]; edgI[1] = edgInd[1];
	E1.IndexOfFace[0] = faces[0]; E1.IndexOfFace[1] = faces[1];
	//nahrazeni starych hran v trojuhelnicich hranami novymi
	for (int i = 0; i < 3; i++)
	{
		if (mesh.CTriangles[E1.IndexOfFace[0]].Edges[i] == edgI[0] || mesh.CTriangles[E1.IndexOfFace[0]].Edges[i] == edgI[1])
		{
			mesh.CTriangles[E1.IndexOfFace[0]].Edges[i] = indexH;
		}
		if (mesh.CTriangles[E1.IndexOfFace[1]].Edges[i] == edgI[0] || mesh.CTriangles[E1.IndexOfFace[1]].Edges[i] == edgI[1])
		{
			mesh.CTriangles[E1.IndexOfFace[1]].Edges[i] = indexH;
		}
	}
}
void updateEdges(CMesh& mesh, vector<int> edges, MapVertex& vertices, VectPair& VP)
{
	for (vector<int>::iterator it = edges.begin(); it != edges.end(); it++)
	{
		for (VectPair::iterator ip = VP.begin(); ip != VP.end();)
		{
			if ((*ip).second == (*it))
			{
				ip = VP.erase(ip);
			}
			else
			{
				ip++;
			}
		}
		if (!flipCheck((*it), mesh.EdgesWithInfo, mesh.edgAdjToPoint))
		{
			QMatrix Qc, Qcd, Qcdi;
			for (int i = 0; i < 4; i++)
			{
				for (int j = 0; j < 4; j++)
				{

					Qc[i][j] = vertices[mesh.EdgesWithInfo[(*it)].Verts[0]].Q[i][j] + vertices[mesh.EdgesWithInfo[(*it)].Verts[1]].Q[i][j];
					if (i < 3)
					{
						Qcd[i][j] = Qc[i][j];
					}
				}
			}
			Qcd[3] = { 0.0, 0.0, 0.0, 1.0 };
			if (inverse(Qcd, Qcdi))
			{
				mesh.EdgesWithInfo[(*it)].Vc.X = Qcdi[0][3];
				mesh.EdgesWithInfo[(*it)].Vc.Y = Qcdi[1][3];
				mesh.EdgesWithInfo[(*it)].Vc.Z = Qcdi[2][3];
				mesh.EdgesWithInfo[(*it)].Vc.Q = Qc;
				mesh.EdgesWithInfo[(*it)].Error = errFcn(mesh.EdgesWithInfo[(*it)].Vc, Qc);
				VP.push_back(make_pair(mesh.EdgesWithInfo[(*it)].Error, (*it)));
			}
			else
			{
				CVertex avg, A, B;
				A = vertices[mesh.EdgesWithInfo[(*it)].Verts[0]];
				B = vertices[mesh.EdgesWithInfo[(*it)].Verts[1]];

				avg.X = (A.X + B.X) / 2;	avg.Y = (A.Y + B.Y) / 2;	avg.Z = (A.Z + B.Z) / 2;
				double a = errFcn(A, Qc);
				double b = errFcn(B, Qc);
				double c = errFcn(avg, Qc);
				double m = min(min(a, b), c);

				if (m == a)
				{
					mesh.EdgesWithInfo[(*it)].Vc = A;
					mesh.EdgesWithInfo[(*it)].Vc.Q = Qc;
					mesh.EdgesWithInfo[(*it)].Error = a;
					VP.push_back(make_pair(a, (*it)));
				}
				else if (m == b)
				{
					mesh.EdgesWithInfo[(*it)].Vc = B;
					mesh.EdgesWithInfo[(*it)].Vc.Q = Qc;
					mesh.EdgesWithInfo[(*it)].Error = b;
					VP.push_back(make_pair(b, (*it)));
				}
				else
				{
					mesh.EdgesWithInfo[(*it)].Vc = avg;
					mesh.EdgesWithInfo[(*it)].Vc.Q = Qc;
					mesh.EdgesWithInfo[(*it)].Error = c;
					VP.push_back(make_pair(c, (*it)));
				}

			}
		}
	}
	sort(VP.begin(), VP.end());
}
void dec(CMesh& mesh, VectPair VP, double level, MapVertex& vertices)
{
	int iter = round(((mesh.CTriangles.size() / 100.0)*(100.0 - level)) / 2);
	cout << "puvodni pocet trojuhelniku:  " << mesh.CTriangles.size() << endl;
	cout << "iter " << iter << endl;
	for (int i = 0; i < iter; i++)
	{
		if (VP.size() != 0)
		{
			for (VectPair::iterator it = VP.begin(); it != VP.end(); it++)
			{
				//cout << "chyba: " << (*it).first << "  "<< (*it).second << endl;
			}

			//1) novy bod
			MapVertex::iterator mv = vertices.end();
			MapEdges::iterator me = mesh.EdgesWithInfo.end();
			mv--; me--;
			int indexv = (*mv).first + 1; int indexe1 = (*me).first + 1; int indexe2 = (*me).first + 2;
			vertices[indexv] = mesh.EdgesWithInfo[VP[0].second].Vc;

			//indexy mazane hrany, bodu, trojuhelniku
			Array2i deletVert = mesh.EdgesWithInfo[VP[0].second].Verts;
			//cout << "body ke smazani:  "<< deletVert[0]<<"  "<< deletVert[1] << endl;

			Array2i deletTri = mesh.EdgesWithInfo[VP[0].second].IndexOfFace;
			Array2i dePom;
			array<int, 5> deletEdg;
			deletEdg[0] = VP[0].second;

			//2) nove hrany
			CEdge E1, E2;
			E1.Verts[0] = indexv;
			E2.Verts[0] = indexv;

			createEdge(mesh, VP, deletTri[0], deletVert, E1, indexe1, dePom);
			deletEdg[1] = dePom[0]; deletEdg[2] = dePom[1];

			createEdge(mesh, VP, deletTri[1], deletVert, E2, indexe2, dePom);
			deletEdg[3] = dePom[0]; deletEdg[4] = dePom[1];
			mesh.EdgesWithInfo[indexe1] = E1;
			mesh.EdgesWithInfo[indexe2] = E2;

			//smazani hran z edgAdjToPoint + vlozeni novych
			for (array<int, 5>::iterator it = deletEdg.begin(); it != deletEdg.end(); it++)
			{
				vector<int>::iterator ie;
				ie = find(mesh.edgAdjToPoint[mesh.EdgesWithInfo[(*it)].Verts[0]].begin(), mesh.edgAdjToPoint[mesh.EdgesWithInfo[(*it)].Verts[0]].end(), (*it));
				if (ie != mesh.edgAdjToPoint[mesh.EdgesWithInfo[(*it)].Verts[0]].end())
				{
					mesh.edgAdjToPoint[mesh.EdgesWithInfo[(*it)].Verts[0]].erase(ie);
				}

				ie = find(mesh.edgAdjToPoint[mesh.EdgesWithInfo[(*it)].Verts[1]].begin(), mesh.edgAdjToPoint[mesh.EdgesWithInfo[(*it)].Verts[1]].end(), (*it));
				if (ie != mesh.edgAdjToPoint[mesh.EdgesWithInfo[(*it)].Verts[1]].end())
				{
					mesh.edgAdjToPoint[mesh.EdgesWithInfo[(*it)].Verts[1]].erase(ie);
				}
				for (VectPair::iterator ip = VP.begin(); ip != VP.end();)
				{
					if ((*ip).second == (*it))
					{
						ip = VP.erase(ip);
					}
					else
					{
						ip++;
					}
				}
			}
			mesh.edgAdjToPoint[E1.Verts[0]].push_back(indexe1);
			mesh.edgAdjToPoint[E1.Verts[1]].push_back(indexe1);
			mesh.edgAdjToPoint[E2.Verts[0]].push_back(indexe2);
			mesh.edgAdjToPoint[E2.Verts[1]].push_back(indexe2);

			//nahrazeni puvodnich bodu v prilehlych hranach bodem novym
			for (vector<int>::iterator it = mesh.edgAdjToPoint[deletVert[0]].begin(); it != mesh.edgAdjToPoint[deletVert[0]].end(); it++)
			{
				if (mesh.EdgesWithInfo[(*it)].Verts[0] == deletVert[0])
				{
					mesh.EdgesWithInfo[(*it)].Verts[0] = indexv;
				}
				else if (mesh.EdgesWithInfo[(*it)].Verts[1] == deletVert[0])
				{
					mesh.EdgesWithInfo[(*it)].Verts[1] = indexv;
				}
			}
			for (vector<int>::iterator it = mesh.edgAdjToPoint[deletVert[1]].begin(); it != mesh.edgAdjToPoint[deletVert[1]].end(); it++)
			{
				if (mesh.EdgesWithInfo[(*it)].Verts[0] == deletVert[1])
				{
					mesh.EdgesWithInfo[(*it)].Verts[0] = indexv;
				}
				else if (mesh.EdgesWithInfo[(*it)].Verts[1] == deletVert[1])
				{
					mesh.EdgesWithInfo[(*it)].Verts[1] = indexv;
				}
			}

			// prilehle hrany novemu bodu
			vector<int> v1 = mesh.edgAdjToPoint[deletVert[0]];
			for (vector<int>::iterator it = mesh.edgAdjToPoint[deletVert[1]].begin(); it != mesh.edgAdjToPoint[deletVert[1]].end(); it++)
			{
				vector<int>::iterator ite = find(v1.begin(), v1.end(), (*it));
				if (ite == v1.end() && (*it) != deletEdg[0])	///////////
				{
					v1.push_back((*it));
				}
			}
			v1.push_back((indexe1)); v1.push_back((indexe2));
			mesh.edgAdjToPoint[indexv] = v1;
			//mazani hrany, body, trojuhelniky
			mesh.edgAdjToPoint.erase(deletVert[0]); mesh.edgAdjToPoint.erase(deletVert[1]);

			mesh.EdgesWithInfo.erase(deletEdg[0]); mesh.EdgesWithInfo.erase(deletEdg[1]); mesh.EdgesWithInfo.erase(deletEdg[2]); mesh.EdgesWithInfo.erase(deletEdg[3]); mesh.EdgesWithInfo.erase(deletEdg[4]);

			vertices.erase(deletVert[0]);
			vertices.erase(deletVert[1]);

			mesh.CTriangles.erase(deletTri[0]);
			mesh.CTriangles.erase(deletTri[1]);
			updateEdges(mesh, mesh.edgAdjToPoint[indexv], vertices, VP);
		}
		else
		{
			break;
		}
	}
	cout << "finalni pocet trojuhelniku:  " << mesh.CTriangles.size() << endl;
}
void writeToSTL(CMesh mesh, MapVertex& vertices, ofstream& stlFile)
{

	stlFile << "solid" << '\n';
	cout << "pocet trojuhelniku " << mesh.CTriangles.size() << endl;
	for (MapTriangle::iterator it = mesh.CTriangles.begin(); it != mesh.CTriangles.end(); it++)
	{
		CVertex A, B, C;
		Array3D n;
		int a = mesh.EdgesWithInfo[(*it).second.Edges[0]].Verts[0];
		int b = mesh.EdgesWithInfo[(*it).second.Edges[0]].Verts[1];
		A = vertices[a];
		B = vertices[b];

		if (mesh.EdgesWithInfo[(*it).second.Edges[1]].Verts[0] == a || mesh.EdgesWithInfo[(*it).second.Edges[1]].Verts[0] == b)
		{
			C = vertices[mesh.EdgesWithInfo[(*it).second.Edges[1]].Verts[1]];
		}
		else
		{
			C = vertices[mesh.EdgesWithInfo[(*it).second.Edges[1]].Verts[0]];
		}
		n = normal(A, B, C);

		stlFile << "facet normal " << n[0] << " " << n[1] << " " << n[2] << '\n';
		stlFile << "outer loop" << '\n';
		stlFile << "vertex " << A.X << " " << A.Y << " " << A.Z << '\n';
		stlFile << "vertex " << B.X << " " << B.Y << " " << B.Z << '\n';
		stlFile << "vertex " << C.X << " " << C.Y << " " << C.Z << '\n';
		stlFile << "endloop" << '\n';
		stlFile << "endfacet" << '\n';
	}
	stlFile << "endsolid";
	stlFile.close();
}
void decimate(string path, string name, double level)
{
	VectorVertex vertices;
	VectArr3i triangles;
	MapVertex uniqueVertices;
	MapEdges edges;
	VectPair VP;

	ifstream sourceFile(path, ios::binary);
	ofstream outputFile(name);

	if (isAscii(sourceFile))
	{
		asciiFile(sourceFile, vertices);
	}
	else binaryFile(sourceFile, vertices);
	getUniqVertAndTrianglesList(vertices, triangles, uniqueVertices);
	cout << "pocet bodu orig: " << uniqueVertices.size() << endl;
	CMesh mesh(triangles);

	getVertexMatrices(mesh, uniqueVertices);
	getErrors(mesh, uniqueVertices, VP);
	dec(mesh, VP, level, uniqueVertices);

	writeToSTL(mesh, uniqueVertices, outputFile);
	sourceFile.close();

}


int main()
{


	decimate("..\\elephant.stl", "object1.stl", 10.0);



	cin.get();
	return 0;
}


