#ifndef _MESH_H_
#define _MESH_H_

#define _CRT_SECURE_NO_DEPRECATE //i won't use fopen_s; so, don't warn me

#include <iostream>
#include <vector>
#include <set>
#include <queue>
using namespace std;

#include "FibHeap.h"
#include "Blossom/Blossom.h" //efficient min-weight perfect matching, 15secs for 693 vs. 693 verts

#define PI							3.141592653f
#define	ToRadian(degree)			(degree * PI/180.0f) //sin(30 degree) = sin(0.52 radian) = 0.5
#define ROUNDING_ERR				0.01f //e.g. duplicate coordinate detection
#define INF							9000000000000000000000000000000.0f //largest float is 3.4x10^38
#define EXP							2.718281828f //exponential e to be used in Gaussian kernel
#define SCALE_FACTOR				1000 //scale original .obj files to deal w/ big area tris (due to precision weirdnesses, non-zeros may be taken as 0s undesiredly; this prevents possibility of such fake 0-area triangles)
//theInput constants
#define SCAPE						1
#define TOSCA_CAT					2
#define TOSCA_CENTAUR				3
#define TOSCA_DAVID					4
#define TOSCA_DOG					5
#define TOSCA_GORILLA				6
#define TOSCA_HORSE					7
#define TOSCA_VICTORIA				8
#define TOSCA_WOLF					9
#define TOSCA_DAVIDVICTORIA			10

struct Triangle
{
	//tris[idx] is this tri
	int idx;

	//idx to verts forming this tri
	int v1i, v2i, v3i;

	//normal and unit normal
	float* normal, * unormal;

	//unsigned area
	float area;

	//not to use the same tri multiple times for the area computation of the same patch
	bool usedInPatchAreaComputation;

	Triangle(int i, int i1, int i2, int i3) : idx(i), v1i(i1), v2i(i2), v3i(i3) {};

	void setNormal(float* v1, float* v2, float* v3)
	{
		normal = new float[3];
		unormal = new float[3];

		//normal is the A, B, C triple of Ax + By + Cz + D = 0; and unormal is the unitized version (direction)
		normal[0] = v1[1]*(v2[2]-v3[2]) + v2[1]*(v3[2]-v1[2]) + v3[1]*(v1[2]-v2[2]);
		normal[1] = v1[2]*(v2[0]-v3[0]) + v2[2]*(v3[0]-v1[0]) + v3[2]*(v1[0]-v2[0]);
		normal[2] = v1[0]*(v2[1]-v3[1]) + v2[0]*(v3[1]-v1[1]) + v3[0]*(v1[1]-v2[1]);
		//D = -v1[0]*(v2[1]*v3[2] - v3[1]*v2[2]) - v2[0]*(v3[1]*v1[2] - v1[1]*v3[2]) - v3[0]*(v1[1]*v2[2] - v2[1]*v1[2]);
		//D is needed only for distance from a point to tri computation; not needed in snake 3d

		//length of this vector
		float length = sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);

		//note that, contractToV2 or middle may produce temporary 0-area (due to collinearity) tris
		if (length >= ROUNDING_ERR) //0.00001f thrshld is bad 'cos 0.0000854606 and, undesiredly, 0.0000854606 >= 0.00001f holds
		{
			unormal[0] = normal[0] / length;
			unormal[1] = normal[1] / length;
			unormal[2] = normal[2] / length; //frequent case
		}
		else
			unormal[0] = unormal[1] = unormal[2] = 0.0f; //rare case
	}

	void setArea(float* v1, float* v2, float* v3)
	{
		// If a triangle is defined by 3 points, say p, q and r, then
		// its area is 0.5 * length of ((p - r) cross (q - r)) [See Real-Time Rendering book, Appendix A]
/*		SbVec3f p_r = v2 - v1;
		SbVec3f q_r = v3 - v1;
		SbVec3f cross = p_r.cross(q_r);*/
		float* p_r = new float[3];
		for (int i = 0; i < 3; i++)
			p_r[i] = v2[i] - v1[i];

		float* q_r = new float[3];
		for (int i = 0; i < 3; i++)
			q_r[i] = v3[i] - v1[i];

		float* cross = new float[3];
		cross[0] = p_r[1]*q_r[2] - p_r[2]*q_r[1];
		cross[1] = p_r[2]*q_r[0] - p_r[0]*q_r[2];
		cross[2] = p_r[0]*q_r[1] - p_r[1]*q_r[0];

		//length of this vector
		float length = sqrt(cross[0]*cross[0] + cross[1]*cross[1] + cross[2]*cross[2]);

		area = 0.5f * length;

		//delete a, b, c; ignores b and c
		delete p_r;
		delete q_r;
		delete cross;
	}
};

struct Vertex
{
	//verts[idx] is this vert
	int idx, grdTruthOnMesh1; //generally ground-truth matches are i-i so this is redundant but for the pairs where i manually set the ground-truths i need grdTruthOnMesh1 trick

	//spatial coordinates (on traditional Cartesian frame of xyz-axes)
	float* coords;

	//normal and unit normal
	float* normal, * unormal;

	//idx to vert neighbors of this vert (use set for efficient duplicate entry prevention)
	typedef set< int, std::less< int > > vertNeighborsIS; vertNeighborsIS vertNeighborsSet;
	
	//idx to sedges of this vert; in addition to vertNeighborsSet, this vert is virtually connected to other verts, decided by sedges
	vector< int > sedgeList; //shortcut edges to smooth out geodesic computation; not keeping virtualVertNeighbors sort of thing 'cos i need sedge length in projected domain, which cannot be computed after genShortcutEdges()

	//idx to tri neighbs and to edges incident to this vertex
	vector< int > triNeighbors, edgeList;

	//dSpanningUnn spans whole mesh so that [9].dSpanningUnn[4] gives g(verts[9], verts[4]) in the original scale; normalized into [0,1] version is dSpanning[] s.t. max geodesic over the surface is 1
	float* dSpanningUnn, * dSpanning, dSpanned, d;

	//corresponding heapNode for this vert to be used in DecreaseKey(heapNode, newKey)
	FibHeapNode* heapNode;

	//true if this vert is available to be an evenly-spaced base to be used in mu computation
	bool baseCandidate, processed;

	//average geodesic distance from this vert to all other verts
	float Mu, * spectralD, diso; //latter for spectral embedding and last is the individual isometric distortion of the match ending at this vert

	Vertex(int i, float* c) : idx(i), coords(c), baseCandidate(true), Mu(0.0f), dSpanningUnn(NULL) {};

	bool addVertNeighbor(int vertInd)
	{
		//add vertInd as a vert neighb of this vert if it won't cause a duplication

		pair< vertNeighborsIS::const_iterator, bool > pa;
		//O(log N) search for duplicates by the use of red-black tree (balanced binary search tree)
		pa = vertNeighborsSet.insert(vertInd);
		//pa.second is false if item already exists, and true if insertion succeeds
		return pa.second;
	}
	
	void removeVertNeighbor(int toBeRemoved)
	{
		//size_type erase (const key_type& x): Deletes all the elements matching x. Returns the 
		//number of elements erased. Since a set supports unique keys, erase will always return 1 or 0

		vertNeighborsSet.erase(toBeRemoved);
	}

	void removeTriNeighbor(int t)
	{
//this fnctn could've been used by contract/split/flip as well but i wrote it later; so, only createSphere() uses it

		int removeIdx;
		for (removeIdx = 0; removeIdx < (int) triNeighbors.size(); removeIdx++)
			if (triNeighbors[removeIdx] == t)
				break;

		if (removeIdx == triNeighbors.size()) { cout << "tri" << t << " not in vert" << idx << ".triNeighbs!\n"; exit(0); }

		triNeighbors.erase(triNeighbors.begin() + removeIdx);
	}

	void setNormal(float* nnormal)
	{
		normal = new float[3];
		unormal = new float[3];
		for (int i = 0; i < 3; i++)
			normal[i] = nnormal[i];

		//length of this vector
		float length = sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
		//if (length >= ROUNDING_ERR) //0.00001f thrshld is bad 'cos 0.0000854606 and, undesiredly, 0.0000854606 >= 0.00001f holds
		if (length > 0)
		{
			unormal[0] = normal[0] / length;
			unormal[1] = normal[1] / length;
			unormal[2] = normal[2] / length; //frequent case
		}
		else
		{
			unormal[0] = unormal[1] = unormal[2] = 0.0f; //rare case
			cout << idx << " has 0000 vertex normal!\t" << length << endl; //happens if this vert has only 2 tri neighbs which are temporary 0-area
		}
	}
};

struct Edge
{
	//edges[idx] is this edge
	int idx;

	//idx to endpnt verts of this edge
	int v1i, v2i;

	//distance between v1i & v2i
	float length;

	Edge(int i, int i1, int i2, float l) : idx(i), v1i(i1), v2i(i2), length(l) {};
};

struct MultiSetContent
{
	//identify each edge by their unique idx field
	int idx;

	//the key for sorting
	float length;

	MultiSetContent(int i, float l) : idx(i), length(l) {};

	//operator< is used to order content by their lengths
	bool operator<(const MultiSetContent c) const {	return length < c.length; }
};

struct OutputPkt
{
	//matchWith()'ll return so complicated objects that i need a struct to package them
	vector< vector< int > > candidateMapsSorted;
	vector< float > distortionsSorted;
};

class Mesh
{
public:
	//use vector to allow random access of elements in O(1) time (to maintain index consistency, never physically erase elements)
	vector< Triangle* > tris;
	vector< Vertex* > verts;
	vector< Edge* > edges;
	vector< Edge* > sedges; //shortcut edges used for better geodesic distance computation

	//id of this mesh and which dataset this mesh is from and in which order
	int frameID, theInput, ordering;
	//min/max geodesics over the surface for geodesic normalization into [0,1]
	float minGeoDist, maxGeoDist;
	//geodesic affinity matrix for spectraql embedding
	float** affinity;

	Mesh(int fid, char* format, int mIdx);// : shortcutEdges(false) {}; can't do here or it complains Mesh already has a body

	void loadObj(char* meshFile);
	void loadObj2(char* meshFile);
	void loadOff(char* meshFile);
	void loadPly(char* meshFile);

	void computeExtremes(int N);
	void computeEvenlySpacedSamples(int N, int D);
	OutputPkt matchWith(Mesh* mesh2, int N, int topK);
	OutputPkt spectralMatchWith(Mesh* mesh2, int N, int D, int topK);
	inline float getCorrespListDistortion(Mesh* mesh2, vector< int > tmpSafeCorrespTry);
	vector< int > denseMatchWith(Mesh* mesh2, int N, vector< int > goodMatches);
private:
	//edge information for display only
	float minEdgeLen, maxEdgeLen, edgeLenTotal, minSedgeLen, maxSedgeLen, sedgeLenTotal, avgEdgeLen;

	//use multiset to place all edge indices in a heap keyed on length (cost) with the min length at the top
	typedef multiset< MultiSetContent, std::less< MultiSetContent > > theEdgeMS; theEdgeMS edgeMultiSet;

	//permutation mtrxs, e.g. perms3[1] = {2,0,1} and it has 3! = 6 such rows/arrays of size 3
	int** perms2, ** perms3, ** perms4, ** perms5, ** perms6, ** perms7, ** perms8, ** perms9;

	//use extra shortcut edges to smooth out dijkstra distance computations over the surface
	bool shortcutEdges;

	//evenly-spaced bases used in mu computation and few extreme points to be matched; samples is for evenly-spaced samples to be matched for another scenario
	vector< int > baseVerts,  extremes, samples;

	//spectral embedding stuff
	float* eigen_values, ** eigen_vectors;

	void addVertexND(float* c);
	void addTriangle(int v1i, int v2i, int v3i);
	void addEdge(int v1i, int v2i);
	inline float calcTriArea(int t);
	int* getTrianglesSharedBy(Edge* bigEdge, int v1i = -1, int v2i = -1);
	bool vertsCollinear(float* v1, float* v2, float* v3);
	int getEdgeSharedBy(int vi, int vj);
	void genShortcutEdges();
	float* projectVertOntoPlane(int toBeProjected, float* planeUnormal, float distOfPlaneFromOrigin);
	inline void addShortcutEdge(int v1i, int v2i, float* v1, float* v2);
	void updateEdgeLengths();
	void updateSedgeLengths();
	void prepareForGeodesicComputation();
	void updateNormals();

	void computeBaseVerts(float r);
	inline float computePatchArea(int bv, float r);
	
	void computeAffinityMtrx(int N, int D);
	void PCAonAffinity(int N, int D);
	void translateAndScaleSpectralBasesK(int N, int D);

	void dijkstraShortestPaths(int sourceVert);//, float far = INF);

	void fillPermutationMtrxs(int N);	
};

#endif
