#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 <assert.h>

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

//#define PROBABILISTIC_COST	//for EM_ALGO, i use either deviation from isometry for matching i to j or abs-log'ed probability cost of matching i to j

#define SUBDIV_TO_5_BASES	//subdivides each patch in a given level (at all times) to ~5 tmpBases; since # tmpBases is low, i can match them w/ combinations/permutations, i.e. w/o blossom stuff

#define SAMPLING_V2	//samples all bases of current level over the whole surface via COES and then subset of close bases are used as the newcomers of the parent/generator patch; if not defined, use computeTmpBaseVerts2() which uses closeToPatchBorder() (which might leave uncovered regions and therefore bad)

#define OPTIMAL_R_BY_FORMULA //use a formula that gives optimal r given the biggest-area over all patches; this r is expected to distribute ~5 bases on each patch, at least should not exceed ~5 to keep comb/permutations fast; if not defined, make a binary search w/ trial and error on changing r values

#define MAX_M_FOR_COMBINATIONS		10 //C(m, 5) combinations tried for the traversalList creation (for isometryCosts[][]) where m cannot exceed MAX_M_FOR_COMBINATIONS for efficiency

//#define ENFORCE_SMOOTH_MAPPING	//defined: tmpSafeCorrespTry tried by getErrOfCorrespListCP()'ll be discouraged if its worst match is too big (even if it has a small original isometryCostTry), hence enforcing smooth mappings

//#define DBG_DETAIL	//define to print lots of detail that are useful during (part of) my agonizing debugging process

#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							90000000000000.0f
#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

//theInput constants
#define DANCING_JUMPING				1 //hybrid
#define DOG_WOLF					2 //hybrid
#define BALLERINA					3 //hybrid [this now corresponds to huge Watertight (SHREC'07) dataset which involves ballerina pairs in human class]
#define JUMPING						6 //not hybrid, i.e. from the same animation sequence
#define DANCING						7 //not hybrid, i.e. from the same animation sequence
#define HORSE						8 //not hybrid, i.e. from the same animation sequence
#define WOLF_DANCING				9 //hybrid
#define BALLERINA_JUMPING			10 //hybrid
#define WOLF_HORSE					11 //hybrid
#define DOG_HORSE					12 //hybrid
#define CAT							13 //hybrid
#define DOG							14 //hybrid
#define LION						15 //hybrid
#define HORSE2						16 //hybrid
#define CAT_DOG						17 //hybrid
#define CAT_WOLF					18 //hybrid
#define WOLF_HORSE2					19 //hybrid
#define DAVID						20 //hybrid
#define VICTORIA					21 //hybrid
#define MICHEAL						22 //hybrid
#define DAVID_VICTORIA				23 //hybrid
#define MICHAEL_VICTORIA			24 //hybrid
#define DAVID_MICHAEL				25 //hybrid
#define SCAPE						26
#define SHREC11						27
#define TOSCA						28 //high-reso grd-truth versions of nonrigid world data, e.g. CAT, DOG, .. stuff
//isometry deviation, best to worst: jumping (0.030), ballerina (0.034), dancing (0.036), horse (0.039), dancing-jumping (0.046), dog-wolf (0.060), where max dist over whole mesh is normalized to 1
//within animals 13..16 (0.044), within humans 20..22 (0.014), across humans 23..25 (0.026)

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 2+ times for patch area computation
	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;

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

	//k-D spectral coords (1st 3 coords of drawn for visualization)
	float* spectralK;

	//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
	vector< int > triNeighbors;

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

	//d of Dijkstra in Algo textbook, i.e. the shortest path estimate (from source baseVerts[i] to this v, during the computation of baseVerts[i])
	float d, d2; //also used for BFS, shortest dist from source to this vert, assuming dist = # of edges, i.e. unweighted shortest path; when a 2nd fnctn updates d, the caller fnctn's d suck; so, use d2 in such 2nd fnctncs such as tmpBaseCandidateFalse()

	//geodesic-based function: (averaged) sum of geodesic distances from this vert to all other verts
	float Mu;

	//similar to d above, but that d is valid for verts within a neighborhood of baseVerts[i]; this dSpanned spans whole mesh
	float dSpanned;
	//dSpanning spans whole mesh so that [9].dSpanning[4] gives g(verts[9], verts[4]) where verts[9] is base, i.e. geodesic distance b/w v9 (which is base) and v4 (which may be nonbase); note that, dSpanning is initialized only for baseVerts and will be normalized later
	float* dSpanning, * dSpanningUnn; //Unn keeps the unnormalized values at all times (recall, dSpanning is normalized)
//	vector< float >* geodesics, * geodesicsUnn; //[2*i]=7, [2*i + 1]=0.4 means normalized geodesic dist from this vert to verts[7] is 0.4; so, every two consec elements is related, idx and geodesic dist to vert at that idx; *Unn is the unnormalized version
//vector.clear() not dealloc memo (it just calls destructors of its elements, which makes no good 'cos my elements here are int/float anyway); so, use vector ptr w/ new and delete it when done for memo efficiency

	float avgGd2Roots; //when roots has 1 element, this reduces to geodesic distance b/w this vert and root vertex (a function value like Mu, but more stable for symmetric parts)
					   //when roots has 2+ elements, instead of 1 root, use 2+ roots and therefore average geodesic distance to roots [descriptor]

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

	//some verts will be base for each level of computation; our curvature-oriented evenly-spaced (COES) sampling decides bases
	bool baseCandidate, base, //is this a baseCandidate or a baseVert?
		 tmpBaseCandidate, //tmp bases that are uniformly distributed within a patch, for pairUpPatchesByGraphMatching2()
		 inPatch, //tmp bases better not be selected from ~borders of to-be-subdivided patch; picking close to borders bad 'cos the next bv coming computeTmpBaseVerts2() may share such a newcomer c, where bv-c dist is ~0, i.e. uniformity sucks
		 closeToBorder, //v is close to patch border if its new patch (of radius r) stays within current patch p as a whole, where p is the big one getting subdivided in the ctf algo
		 residue, //in matchSmallTmpBases(), e.g. for the 3 vs. 2 patch case, 2 vs. 2 match is set via permutations, and 1 residue in patch1 gets its match w/o permutations
		 useInSafeList; //some matches of prev level are marked as n/a to be used in safe list (for matching of the current level) due to their high iso cost compared to general D_iso (all costs computed over trimmed/optimal/close-to-1-to-1 corresp list)

	float patchArea; //the approximated area of the patch represented by this base vert
	int patchIdx; //call each patchArea a patch, and this is just to color each patch in different color; not maintained for patches of tmpBases

	//becomes true only when this v is added to safeCorresp list
	bool safe;

	vector< int > baseNghbs, //inited only for baseVerts; list of surrounding base neighbs of this base vert
				 * baseNghbhood, //inited only for baseVerts; involves verts in the neighborhood of this baseVert (trying to prevent duplicates took big time in baseVert creation; so, use vector w/ duplicates); ptr to delete it for memo during levels
				 * tmpBaseNghbhood; //for pairUpPatchesByGraphMatching2() that uses tmp patches within a base patch; ptr to delete it for memo during levels

	int matchIdx, matchIdxCpy, //matchIdx = 7 means this vert matches w/ m2.vert7
		nMatches, //defined for mesh2 verts only; maintaining matchIdxs[] is hard; nMatches is easier (always 1 for mesh1.baseVerts) and 1+ for mesh2.baseVerts; it is the # of matches of this vert and updated only in greedyForBaseUpdate(), only place for 1-to-many matchings
		matchIdxCoarse; //defined for mesh1.bases only for dense matching; remember the coarse results to roll-back after force1to1 effect
	vector< float > matchIdxs; //saves matches for v from each touching patch, as 2 consecutive elements matchIdx & matchCost, and may decide the best later if i want to have 1-to-1 mapping from many-to-many
	vector< float > matchIndexes; //even idx are same as matchIdxs but the followig odd idxs hold the isometry cost w.r.t. whole mesh, i.e. safeList spans the whole mesh, unlike matchIdxs which is w.r.t. patches, i.e. safeList is around patch being matched (in matchSmallTmpBases())
	float matchCostI, matchCostICpy, //isometry, i.e. geodesic consistency, cost
		  matchCostG, //ground-truth geodesic distance cost
		  matchCostDense, //force1to1 cases in make1to1() uses this costs for correct operation
		  matchProbI, //probability of matching this vert to its match; same as matchCostI but this is prob = 1 - cost; to be linearly proportional to cost, i use abs-log (log for EM framework and abs for linear proportion); see setCorrespondenceErrs() for an example
		  matchCostIResidue; //matchCost of the residue to set the min cost during the bottom of matchSmallTmpBases()

	bool processed, drawMatchLine, //2nd is to draw the matches of a given patch only
		 hasBadIsometryCost, //bottom 9 of coarse result file is regarded as bad
		 top10, //top10 match w.r.t. isometryCost
		 topCurv10, //top10 bases w.r.t. curvature (nonbases not made topCurv10 even if they have higher curvature than bases that are made topCurv10=true)
		 topMu10; //top10 bases w.r.t. mu (nonbases not made topMu10 even if they have higher mu than bases that are made topMu10=true)

	int idxInBaseVerts, //-1 for non-baseVerts; used in painting to fetch correct diffuseColor
		bam, //greedy's matchIdx results of this vert before entering alternating-optimization, i.e. just after greedy, just to count # of correspondence updates made by alternating-optimization (isometryBasedGraphMatching())
		bam2, //initial spectral matching's matchIdx results of this vert before entering alternating-optimization, i.e. same as bam but this holds initial spectral match's results
		worstIin2, worstGin2; //since i don't have matchIdx values available for many-to-many case, i.e. c2f algo, save some special matches, e.g. worstMatchCostI in mesh2 side, explicitly for drawing purposes
		
	
	float avgGd2Bases, //another descriptor that looks like Mun, but this one is the avg geodesic to all neighboring bases unlike Mun which is avg geodesic to all verts
		  gc; //Gaussian curvature as another [descriptor]

	vector< float > isometryCosts, bestIsometryCosts, //remember in this format for use in matchTmpBasesByGraph/matchSmallTmpBases() as well as 1-to-1 conversion: for 4471->4388 w/ cost 0.019232, verts[4471]->isometryCosts[x] = 4388, [x+1] = 0.019232
					geo2coarses; //for interpolation of coarse mapping into full dense mapping, i need geo dists from all verts to matched coarse corresp verts; so, geo2coarses is defined for all verts v
	float isometryCostChosenByBlossom; //blossom picks 4471->4388 while its isometryCost in vertexBasedIso is 0.019232; so, remember verts[4471].isometryCostChosenByBlossom = 0.019232
	int //tmpMatchIdx, //for the tmpBases that make into top5 of matchTmpBasesByGraph(), need this to reach the mesh2 counterpart of this vert; hence, defined for mesh1 verts that act as tmpBase only; [replaced by easy-to-maintain 2*ss+1 access to otherside]
		tmpBaseLevel; //at which level this vert became a tmpBase? e.g. 0 for the initial ~10 baseVerts matched in the spectral domain; also = -1 means this is not a tmpBase

	Vertex(int i, float* c) : idx(i), coords(c), baseCandidate(true), base(false), Mu(0.0f), safe(false), drawMatchLine(false), top10(false), topCurv10(false), topMu10(false),
							   matchIdx(-1), matchCostI(INF), matchCostG(INF), idxInBaseVerts(-1), patchIdx(-1), tmpBaseLevel(-1), dSpanning(NULL), baseNghbhood(NULL), tmpBaseNghbhood(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

		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(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
		}
	}

	/*inline float getGeodesic(int v)
	{
		//returns normalized ([0,1]) geodesic dist from this vert to v

		int gSize = (int) geodesics->size() / 2;
		for (int i = 0; i < gSize; i++)
			if ((*geodesics)[2*i] == v)
				return (*geodesics)[2*i + 1];
		//if normalized geodesic n/a, then return special -1 to let me launch a new original/notshortcut dijkstra that touches all verts (see farThrshld declaration comments)
		return -1.0f;
	}
	inline float getGeodesicUnn(int v)
	{
		//returns geodesic dist from this vert to v

		int gSize = (int) geodesicsUnn->size() / 2;
		for (int i = 0; i < gSize; i++)
			if ((*geodesicsUnn)[2*i] == v)
				return (*geodesicsUnn)[2*i + 1];
		return INF; //it may validly happen for unnormalized case and must return a big number as explained in computeTmpBaseVerts2().info comment
	}*/
};

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; }
};

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
	vector< int > baseVerts, //base verts that are distributed almost equally on the surface
				  sepBases, //bases separated sufficiently in spectral domain, i.e. none in this list is too close to another, hence no spectral collision
				  tmpBaseVerts, //tmp bases that are uniformly distributed within a patch
				  tmpBaseVertsAll; //tmp bases that are uniformly distributed over the whole mesh

	vector< int > safeCorresp, //safe correspondences, i.e. desired output, will be added to this list in the format of i <- v1 and i+1 <- corresp(v1)
				  safeCorrespAll, //keep all matches from coarseStage to process them in pairUpPatches(); use top10 for isometryCosts of costs[][] however				  
//				  tmpSafeCorresp, //safe corresps b/w tmpBaseVerts, i.e. the verts within given 2 patches, one from source one from target mesh
//				  safeCorrespsAdaptive, //adaptive/dynamic safe list that changes for each vertex, consisting of e.g, geodesically close bases only
				  fallens; //defined for bases only; in spectral domain, other bases falling on top of this base go here
//				  initialCorresp, //result of initial matching (has ~5 corresps)
//				  adaptList; //correspondence list from patch of previous level to be used for the initial correspondence of this level
	vector< int >* tmpSafeCorresp, * safeCorrespsAdaptive, * initialCorresp, * adaptList; //make frequently used ones delete'able for memo efficiency since clear() not dealloc memo

	int nCorresps; //# of correspondences computed

	//edge information for display only
	float minEdgeLen, maxEdgeLen, edgeLenTotal, minSedgeLen, maxSedgeLen, sedgeLenTotal;

	//principalAxes[0] is the eigvector/principalAxis corresponding to largest eigvalue, [1] & [2] follows that in order
	float* principalAxes;

	//holds the 100% accurate shortest distances b/w all pairs of base verts, not all vertices; e.g. [2][5] is dist b/w baseVert2 & 5
	float** baseGeoAffinity;

	//this is mesh of which jumping man frame?
	int frameID;

	//is this the target/last mesh of whole frame sequence, e.g. for 2-frame seq, true for 2nd mesh
	bool targetMesh;

	int worstMatchI, bestMatchI, //display the worst & best match w.r.t. isometry error, i.e. geodesic consistency
		worstMatchG, bestMatchG, //display the worst & best match w.r.t. ground-truth geodesic distance error (always -1 for theInput)
		theInput, //models from which data set?
		fileForSymProblem, //explained at top of coarseToFineMatching()
		levelToEvaluate, //correspondence stored in level&.dat file is getting evaluated
		trackTillLevel, //top10 maps will be tracked until level trackTillLevel, which is typically 4
		nInitialSamples, //non -1 value forces 1-to-1 initial mapping of size nInitialSamples (not many-to-one unlike upcoming levels because 8 initial extreme samples are very widely separated and demanding many-to-1 may cause unrecoverable errors)
		levelToInterpolate, //any existing level on which the interpolation is made
		maxLevelSuggested; //proceed until maxLevel, use -1 for automatic maxLevel leading to the full dense map, other values (< 14) for intermediate results
	
	//overall isometry and ground-truth cost of (whole) correspondence
	float isometryCost, groundTruthCost, maxCostIg, maxCostGg, isometryCostPass1,
		  minGeoDistG, maxGeoDistG, //min/max geo dist around the whole mesh ('cos computed once w/ the initial ~10 bases that span the whole surface)
		  MIN_TRIANGLE_AREA_ALLOWED; //no triangle w/ area below this value exists on this mesh

	int* gcOrder, //for nonarbitrary next base selection, i.e. high curvatures first
	   * muOrder, //for nonarbitrary next base selection, i.e. high mu values first
	   * nOverlaps; //nOverlaps[v] is the # of patches that involves v'th vertex, i.e. 1 for base verts and 1+ for others

	bool rootBasedSymSoln, //use root based gd2r1/2 to align the pairs that exhibit symmetry problem otherwise
		 lessBases, //true makes ~30, false makes ~80 bases for jumping, and similar for other
		 batchRun; //is this a run via batch file?

	//some interesting matches to remember in order to draw w/ a sphere during visualization
	vector< int > interestingMatches;

	//[i] = radius used at level i; to reach parent (or grand-grand-.. parent), keep all radii used in completed levels
	vector< float > radii,
					safeCorrespCosts; //half the size of safeCorresp 'cos for each match (pair of verts) in safeCorresp, i save its matchCostI in safeCorrespCosts; used only by trimCorresps() pass # 3
	
	float r1, r2, //radius of patches of source and target
		  normalizedR1, //= r1/maxGeoDist1; isometryCost is in [0,1]; so, div it by this to get an isometryCost in terms of r accurately
		  totalPatchSize; //to print the avg patch size

	bool denseMatchByGraphMatching, //do the dense matching by min-weight perfect matching, e.g. blossom, or by interpolation
		 useKernel1, //use either kernel1 or kernel2 to embed pnts into spectral domain via geodesic affinity formed in computeBaseGeodesicAffinityMtrx()
		 arrayInsteadOfHeap, //true uses simple arrays instead of fast fibonacci heap in closeToPatchBorder(); heap is fast only in the first level and then becomes inefficient due to excessive redundant inserts; so, switch back to arrays after level 1
		 C2Fmode, //true means EM_ALGO undefined and therefore i use C2F algo of SGP instead of the EM/alternating algo of CVPR/PAMI
		 finalLevel, //true if this is the very last level of c2f process
		 executionTimeOnly, //false: don't do informative setCorrespondenceErrsTmp() and levelInfoToFile() to save time
		 shortcutEdges, //use shortcut edges to improve geodesic distances; increases runtime very slightly and crashes at weird places for x64 runs; so, i can disable them here using false
		 collapsed; //input is collapsed model either for computeCorrespondenceErrorGMDS() or for coarseToFineMatching()

	Mesh(bool trg) : nCorresps(0), targetMesh(trg), nInitialSamples(8), lessBases(false), denseMatchByGraphMatching(false), baseGeoAffinity(NULL), executionTimeOnly(false), collapsed(false),
					 rootBasedSymSoln(false), totalPatchSize(0.0f), useKernel1(false), MIN_TRIANGLE_AREA_ALLOWED(0.001f), fileForSymProblem(2), levelToEvaluate(4),
					 tmpSafeCorresp(NULL), safeCorrespsAdaptive(NULL), initialCorresp(NULL), adaptList(NULL), perms2(NULL), finalLevel(false), muOrder(NULL) {};

	int addVertex(float* v);
	int addVertex(float x, float y, float z);
	void addVertexND(float* c);
	void addTriangle(int v1i, int v2i, int v3i, bool temporary0Area);
	void updateNormals();
	void moreBasesForDrawing();

	int initialSpectralMatching(Mesh* mesh2);
	int initialEuclideanMatching(Mesh* mesh2);
	int greedyOpt(Mesh* mesh2, bool alternatingOpt, int k);
	void finalizeMatching(Mesh* mesh2, int k, bool detailsToFile);

	void correspFromFile(Mesh* mesh2);
	void computeCorrespondenceError(Mesh* mesh2);
	void computeCorrespondenceErrorC2F(Mesh* mesh2, bool compareWithZhang, bool justDisplay = true);
	void computeCorrespondenceErrorBIM(Mesh* mesh2, bool justDisplay = true);
	void computeCorrespondenceErrorGMDS(Mesh* mesh2, bool interpolateCoarseMap);
	void computeInputQuality(Mesh* mesh2);

	void interpolateCoarseCorresp(Mesh* mesh2, bool loadInterpolatedFile);

	void denseMatching(Mesh* mesh2);

	void coarseToFineMatching(Mesh* mesh2); //-1 means no file loading or whatsoever
private:
	//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;

	//shortest distance b/w all pairs of verts held in D; before allPairsShortestPaths() call, D maintains adjacency matrix
	float** geoAffinity; //again holds the shortest distance b/w all pairs of verts, but approximation (hence faster than D)

	vector< int > roots; //use 1+ roots for avgGd2Roots value

	float** eigen_vectors, //nxn eigvecs of geodesic affinity mrtx in descending order of eigvals; eigen_vectors.column1 = eigvector1 (corresponding to largest eigval), col2 = eigvec2, and so on..
		 * eigen_values, //corresponding eigvalues (for geodesic affinity mrtx)
		 ** spectralCoords; //remember this mtrx computed in PCAonAffinity() to be used in distanceBasedTriangulation() as lkPseudoinvTranspose

	FILE* errfPtr; //errs*.dat; used for greedy algo, not for c2f algo

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

	char outFolder[100]; //symHeal mode outputs fwritten into this folder

	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 moveCOMtoOrigin();
	void prepareForGeodesicComputation();	

	void hybridModelsToSameCoordFrame(Mesh* mesh2);
	
	///////// EM_ALGO.coarseCorresp ////////
	float computeBaseVerts(bool dijkstraFromBases, bool COES, bool computeMu = false);
	void computeBaseVerts(float r, bool COES);
	float setRadius(float wholeSurfaceArea, bool initialBasesByMu);
	int getBaseVertIdx(bool COES);
	inline float computePatchArea(int bv, float r);

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

	void computeGaussianCurvatures(); //avg of patch for bases not done anymore; not a big deal anyway (delLaterrrrrr this comment)
	void computeBaseNeighbs(float r);
	void sortCurvatures();
	void sortMus();

	void initRoots();
	void computeGeodesicDistsToRoot(bool print = true);

	void computeBaseGeodesicAffinityMtrx(int k, bool toFile = false);
	void PCAonAffinity(int k, bool pseudoinverseTranspose = false);
	void translateAndScaleSpectralBasesK(int k, bool translateOnly = false);

	void findBestAffinityEigvectors(Mesh* mesh2, int k, bool pseudoinverseTranspose = false);
	void PCAalignmentK(float** &pntSet1, float** &pntSet2, int n1, int n2, int k);
	inline float hybridCost(float* v1Aligned, float* v2Aligned, int k, Vertex* v1, Vertex* v2);
	void distsAmongKDSets(int* idxSet1, int* idxSet2, int n1, int n2, int k, Mesh* mesh2);
	float computeMSEKD(float** pntSet1, float** pntSet2, int n1, int n2, int k, bool div = true);

	float setCorrespondenceErrs(Mesh* mesh2, float r, int when = 2);
	void refillSafeCorresp(Mesh* mesh2, int when = 0);
	void sortSafeCorresp(Mesh* mesh2, char* fName);
	void decideTop10Iso();
	void decideTop10Curv();
	void decideTop10Mu(bool evenlySpaced);
	void setExtremeCorrespondences(float r, bool lightPrint, bool manyToMany = false);
	int adjustBaseVerts();

	void greedyForBaseUpdate(Mesh* mesh2, float r, int iter);
	void isometryBasedGraphMatching(Mesh* mesh2, float r, int* vIdx2, int iter, int k);
	void duplicateFreePatchesToFile(Mesh* mesh2, char* fName);
	void setOutFolder();
	///////// EM_ALGO.coarseCorresp ends ////////

	///////// EM_ALGO.denseCorresp ////////
	void fillBaseGeoAffinity();
	void distanceBasedTriangulation(int k);

	void pairUpPatches(Mesh* mesh2);
	void make1to1(Mesh* mesh2, FILE* fPtrU, int force1to1, bool updateMatchCostI = false, bool baseFlagValid = true);
	void allowManyTo1(Mesh* mesh2, int currLevel, bool updateMatchCostI = true);
	void allow1ToMany(Mesh* mesh2, int currLevel);

	void pairUpPatchesByInterpolation(Mesh* mesh2, FILE* fPtrU, int bv1, int bv2);

	bool adaptiveListForIsometryCost(int bv1, int bv2, int level = -1, Mesh* mesh2 = NULL);
	void pairUpPatchesByGraphMatching(Mesh* mesh2, FILE* fPtrU, int bv1, int bv2);

	float initialEuclideanMatchingFromFile(Mesh* mesh2, int nthMapping, bool firstCall, bool toFile);
	void computeTmpBaseVerts(float r, int bv, bool discardBV);
	bool closeToPatchBorder(int w, float r);
	void computeTmpBaseVerts2(float r, int bv, float rOfPrevLevel, int currentLevel);
	void computeTmpBaseVerts3(int bv);
	void computeTmpBaseVerts4(float r, int bv, float rOfPrevLevel, int currentLevel);
	void computeTmpBaseVerts4Eff(float r, int bv, float rOfPrevLevel, int currentLevel);
	void basesOverWholeMesh(float r, int level);
	int countTmpBaseVerts(float r, int bv);
	void tmpBaseCandidateFalse(int sourceVert, float far);
	void matchTmpBasesByGraph(Mesh* mesh2, int bv1, int bv2, bool useTopnForTraversalList, int level = 0);
	void matchTmpBasesInitial(Mesh* mesh2, int bv1, int bv2, int level);
	//void matchTmpBasesTopn(Mesh* mesh2, int bv1, int bv2);
	void matchTmpBasesC10_5(Mesh* mesh2, int bv1, int bv2);
	float computeTmpSafeCorrespErr(Mesh* mesh2);
	void pairUpPatchesByGraphMatching2(Mesh* mesh2, int bv1, int bv2);
	///////// EM_ALGO.denseCorresp ends ////////

	void refillSafeCorrespTmp(Mesh* mesh2, int currentLevel, bool isoCostPass1 = false);
	void setCorrespondenceErrsTmp(Mesh* mesh2, bool reducedCorresp);
	void setCorrespondenceErrsTmpEff(Mesh* mesh2, bool reducedCorresp, int currentLevel);
	void setCorrespondenceErrsTmpManyTo1(Mesh* mesh2);
	void trimCorresps(Mesh* mesh2, int currentLevel, bool pass1Only, bool isoCostPass1);
	void clearVects2();
	void matchSmallPatches(Mesh* mesh2, int level);
	void matchSmallTmpBases(Mesh* mesh2, bool forPostprocessing = true, int bv1 = -1, int bv2 = -1, int level = 0);
	
	inline float getErrOfCorrespList(Mesh* mesh2, vector< int > tmpSafeCorrespTry, bool saveIsoCosts, bool smoothnessTest = false);
	inline float getErrOfCorrespListCP(Mesh* mesh2, vector< int > tmpSafeCorrespTry, vector< int > parentCorresp, bool saveIsoCosts = false);
	void initialCorrespC10_5(Mesh* mesh2);
	void safeCorrespC10_5(Mesh* mesh2);
	
	void levelInfoToFile(int currentLevel, Mesh* mesh2);
	void resumeFromLevel(int levelToLoad, Mesh* mesh2);
	void fillPermutationMtrxs();
	int refineInitialBasesByGraph(Mesh* mesh2, int k);

//	void readAffinityFromFile(char* fileAddr);
	void geodesicsBetweenBasesToFile(float r, int finalLevel);
	void printUnmatched();
};

#endif
