// 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>


using namespace std;

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, 3>, 3> QMatrix;

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

	CMesh(VectArr3i& triangles)
	{
		int i = 0;
		int j = 0;
		map< Array2i, int> mapOfEdges;
		map< Array2i, int>::iterator imaoe;
		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);
				}

				imaoe = mapOfEdges.find({ (*it)[k], (*it)[(k + 1) % 3] });
				if (imaoe == mapOfEdges.end())
				{
					imaoe = mapOfEdges.find({ (*it)[(k + 1) % 3],(*it)[k] });
					if (imaoe == 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.push_back(f);
						t.Edges[k] = j;
						j++;
					}
					else
					{
						EdgesWithInfo[(*imaoe).second].IndexOfFace[1] = i;
						t.Edges[k] = (*imaoe).second;
					}
				}
				else
				{
					EdgesWithInfo[(*imaoe).second].IndexOfFace[1] = i;
					t.Edges[k] = (*imaoe).second;
				}
			}
			CTriangles.push_back(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, MapVertex& 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, VectArr3d& 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;
			i++;
			uniqueVertexes.push_back(v);
		}
		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, crossProduct;;
	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<0.0001)
	{
		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, VectorVertex& 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 };
		//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];
			}
			bool err = false;
			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;
	}
}
void decimate(string path, string name, double reduction)
{
	VectorVertex vertices;
	VectArr3i triangles;
	VectArr3d uniqueVertexes;
	VectEdges edges;

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

	if (isAscii(sourceFile))
	{
		asciiFile(sourceFile, vertices);
	}
	else binaryFile(sourceFile, vertices);
	getUniqVertAndTrianglesList(vertices, triangles, uniqueVertexes);
	CMesh mesh(triangles);
	getVertexMatrices(mesh, vertices);

}


int main()
{
	decimate("C:\\Users\\Coreli\\Desktop\\Hanka\\práce\\B&R\\H\\Script\\TestMesh\\huge\\EiffelTower.stl", "object.stl", 70.0 );
    return 0;
}


