#ifndef _MESH_H_
#define _MESH_H_

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

//#define HAVE_SINGLEPRECISION_MATH //Eigen package does #include <complex> which conflicts with Inventor/C/basic.h definitions; so prevent them with this in files using include Inventor

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

#include <Inventor/SbLinear.h>
#include <Inventor/nodes/SoSeparator.h>
#include <Inventor/nodes/SoCoordinate3.h>
#include <assert.h>
#include <windows.h> //for main.get_time()

#include "FibHeap.h"
#include "psimpl.h"

#define PI							3.141592653f
#define INF							9000000000000000000000000000000.0f //largest float is 3.4x10^38
//#define USE_GOOD_CONTINUATION //define this to use good continuation rule in geodesic distance computations [bad]; undefine to use standard euclidean edge-length distances
#define RETRIEVAL_APP //define this to launch shape retrieval application, i.e. match sketch w/ lots of 3d models; undefine this to match sketch w/ one 3d model only
//#define PRINT_DETAILS //define this to print program flow which is useful in debugging/development; undefine to reduce prints


#define TOSCA	1 //use TOSCA dataset instead of the defaut Watertight dataset
#define MCGILL	2 //use McGill articulated dataset instead of the defaut Watertight dataset
#define WATERTIGHT 3 //use Watertight dataset

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

	//idx to verts forming this tri
	int v1i, v2i, v3i;
	
	//normal and unit normal
	SbVec3f normal, unormal;

	//unsigned area
	float area;

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

	void setNormal(SbVec3f v1, SbVec3f v2, SbVec3f v3)
	{
		//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

		//note that, contractToV2 or middle may produce temporary 0-area (due to collinearity) tris
		//if (normal.length() >= ROUNDING_ERR) //0.00001f thrshld is bad 'cos 0.0000854606 and, undesiredly, 0.0000854606 >= 0.00001f holds
		if (normal.length() > 0)
			unormal = normal / normal.length(); //frequent case
		else
		{
			unormal = SbVec3f(0.0f, 0.0f, 0.0f); //rare case
			//cout << idx << " has 0000 trrrrrrrrrrrrrrrrrrrrrrrriangle normal!\t" << normal.length() << endl;
		}
	}

	void setArea(SbVec3f v1, SbVec3f v2, SbVec3f 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);

		area = 0.5f * cross.length();
	}
};

struct Vertex
{
	//verts[idx] is this vert
	int idx,
		predecessorV, predecessorE, //predecessor vert/edge of this vert in the current shortest/geodesic path
		matchIdx; //matching sketch endpnt for this mesh vert

	//3D coordinate
	SbVec3f coords;

	//normal and unit normal
	SbVec3f 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 triangle (tris) and tetrahedron (tets) neighbs
	vector< int > triNeighbors, tetNeighbors;

	//idxs of the verts (or edges) for the shortest/geodesic path from this vert to some other verts, e.g. verts in savePathToTargets; idxs stored in reverse order, i.e. geoPaths[2][0]=savePathToTargetV, .., geoPaths[2][last]=this for 3rd ([2]) geoPath
	vector< vector< int > > geoPaths;

	//idx to edges incident to this vertex
	vector< int > edgeList;

	//dSpanningUnn spans whole mesh so that [9].dSpanningUnn[4] gives geodesicDistanceBetween(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;
	//corresponding heapNode for this vert to be used in DecreaseKey(heapNode, newKey)
	FibHeapNode* heapNode;

	bool processed, sample; //is this a sampled vertex, i.e. in samples[]

	Vertex(int i, SbVec3f c) : idx(i), coords(c), dSpanningUnn(NULL), dSpanning(NULL), matchIdx(-1) {};

	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

		unsigned int removeIdx;
		for (removeIdx = 0; removeIdx < 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(SbVec3f nnormal)
	{
		normal = nnormal;
		float vecLen = normal.length();	//sqrt(x*x + y*y + z*z)
		//if (vecLen >= ROUNDING_ERR) //0.00001f thrshld is bad 'cos 0.0000854606 and, undesiredly, 0.0000854606 >= 0.00001f holds
		if (vecLen > 0)
			unormal = normal / vecLen;
		else
		{
				unormal = SbVec3f(0.0f, 0.0f, 0.0f);
//				cout << idx << " has 0000 vertex normal!\t" << vecLen << endl; //happens if this vert has only 2 tri neighbs which are temporary 0-area
		}
	}
};

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

	//idx to endpnt verts[] of this edge; or idx to endpnt pnts[] (SketchPnts) of this stroke
	int v1i, v2i,
		ev1i, ev2i; //or idx to endpnt endpnts[] (SketchPnts) of this stroke, i.e. for strokes, these values are pointers to endpnts[] but v1i/v2i are ptrs to pnts[]

	//distance between v1i & v2i
	float length, slope; //slope of this stroke approximated from a local part of the whole stroke

	//set of pointers to SketchPnts that make this stroke
	vector< int > pnts;
	
	//is this a virtual stroke or a real stroke?
	bool virtua, removed; //is this stroke removed and therefore to be ignroed throughout the program?

	Edge(int i, int i1, int i2, float l) : idx(i), v1i(i1), v2i(i2), length(l) {}; //for mesh
	Edge(int i, bool v) : idx(i), virtua(v), length(0.0f), removed(false) {}; //for sketch
};

struct SketchPnt
{
	//unique string ID of this sketch point
	char* strID;

	//pnts[idx] is this pnt; stroke.pnts[spidx] is this pnt for a given stroke
	int idx, spidx, strokeID, //unique^^ stroke ID of this sketch point, i.e. which non-virtual stroke has this pnt? ^^endpnts of splitted strokes may appear on 1+ non-virtual strokes [not in use]
		predecessorP; //predecessor point of this point in the current shortest/geodesic path

	//is this pnt an endpnt?
	bool endpnt;

	//2d coords of this sketch point
	SbVec2f coords;

	//idx to strokes incident to this sketch point (defined only for pnts that are at the ends of strokes)
	vector< int > strokeList, edgeList; //idx to edges incident to this pnt

	//dSpanningUnn spans whole mesh so that [9].dSpanningUnn[4] gives geodesicDistanceBetween(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;
	//corresponding heapNode for this vert to be used in DecreaseKey(heapNode, newKey)
	FibHeapNode* heapNode;

	bool processed, sample; //is this a sampled vertex, i.e. in samples[]

	//idx to sketch point neighbors of this sketch pnt (use set for efficient duplicate entry prevention)
	typedef set< int, std::less< int > > pntNeighborsIS; pntNeighborsIS pntNeighborsSet;
	
	SketchPnt(char* strid, SbVec2f c, int i) : strID(strid), coords(c), idx(i), dSpanningUnn(NULL), dSpanning(NULL), strokeID(-1), endpnt(false) {};

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

		pair< pntNeighborsIS::const_iterator, bool > pa;
		//O(log N) search for duplicates by the use of red-black tree (balanced binary search tree)
		pa = pntNeighborsSet.insert(pntInd);
		//pa.second is false if item already exists, and true if insertion succeeds
		return pa.second;
	}
};

class Mesh
{
public:
	//triangular mesh stuff
	vector< Triangle* > tris;
	vector< Vertex* > verts;
	vector< Edge* > edges;

	//sketch stuff
	vector< SketchPnt* > pnts, endpnts; //sketch points (ends + intermediates) and sketch endpnts (ends only)
	vector< Edge* > strokes; //sketch strokes

	vector< int > samples; //evenly-spaced samples from verts[] (if this is a surface mesh) or from endpnts[] (if this is a sketch)

	//color and center of mass
	SbVec3f color, com;
	//each scr obj needs a separator to be drawn to screen
	SoSeparator* sep;

	bool sketch, //is this 2.5D sketch or 3D model
		 usePnts; //true: uses pnts[]; false: uses endpnts[] throughout the extremity sampling and matching
	
	//min/max geodesics over the surface for geodesic normalization into [0,1]
	float minGeoDist, maxGeoDist, avgEdgeLen, spacing, z; //put sketch's xy plane on this z coord

	//# of original/initial strokes
	int nOriginalStrokes;
	char classStr[250];

	Mesh(int n, bool s = false) : N(n), sketch(s), z(0.0f), usePnts(true) {};
	Mesh(int n, int type, int inputDB, int modelID = -1);
	void loadSketch(int sketchID);
	void loadSketchDense();
	void loadSketchDense_notInUse(int sid);
	void loadPly(int fid);
	void loadObj(int fid);
	void loadOff(int inputDB, int type = -1);
	void computeEvenlySpacedSamples();
	bool computeEvenlySpacedSamplesSketch();
	void preprocessSketch();
	void setStrokeConnectivity();
	float matchWith(Mesh* sketch);
private:
	int sketchID, frameID, //ids of sketch and mesh
		** perms2, ** perms3, ** perms4, ** perms5, ** perms6, ** perms7, ** perms8, ** perms9, //permutation mtrxs
		N; //# of samples to match (FPS samples)
	float minEdgeLen, maxEdgeLen, edgeLenTotal;
	int addVertex(float* v);
	int addVertex(float x, float y, float z);
	void addVertexND(SbVec3f c);
	int addTriangle(int v1i, int v2i, int v3i);
	void addEdge(int v1i, int v2i, float dist = -1.0f);
	inline float calcTriArea(int t);
	void getBoundaryTrianglesSharedBy(Edge* bigEdge, int* shTris, int v1i = -1, int v2i = -1);
	int getEdgeSharedBy(int vi, int vj);
	void computeNormalsAndAreas();
	int getIdxByName(char* name);
	int getIdxToEndpnts(int p);
	void segmentStrokes();
	void dijkstraShortestPaths(int sourceVert, vector< int >* savePathToTargets = NULL);
	void dijkstraShortestPathsSketchE(int sourcePnt);
	void dijkstraShortestPathsSketch(int sourcePnt);
	void dijkstraShortestPathsSketchGC(int sourcePnt);
	int getAnExtremeVertSketch();
	int getAnExtremeVert();
	void scaleBy(float s);
	void translateBy(SbVec3f t);
	void updateEdgeLengths();
	void computeCenterOfMass();
	void additionalRotations();
	void fillPermutationMtrxs(int N);
	inline float getIsometricDistortion(Mesh* mesh2, vector< int > tmpSafeCorrespTry);
	void densifySketch();
	void computeEvenlySpacedSamplesSketchEuc();
	void setSampleConnetivity();
	void setStrokes();
};

#endif
