#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 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 BIGi						2147483647 //biggest int
#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 take as 0s undesiredly; this prevents possibility of such fake 0s)
#define POTENTIAL_THRESHOLD			0.15f //two pnt sets make a potential triplet/quadruple/.. if geodesic consistency b/w the two is below this (make it INF to use all combinations, i.e. no geodesic pruning)
#define	JUMP_BASED_CONVERGENCE //if defined, worst-removal iterations in conquering step (bottom) of matchBasesViaVoting2() stops w/o a threshold parameter using auto-detected jump value

//theInput constants
#define JUMPING_DANCING				1
#define WOLF_CAT					2
#define BALLERINA					3
#define JUMPING						6
#define DANCING						7
#define WOLF_DANCING				9
#define BALLERINA_JUMPING			10
#define BALLERINA_DANCING			11
#define CENTAUR_HORSE				12
#define BOY_GIRL					13
#define CENTAUR_GIRL				14
#define SHREC11						15
#define VICTORIA_SHREC11			16 //some TOSCA.victoria vs. (manually scaled) partial class of SHREC11
#define HORSE						17 //low-reso horses from non-rigid world
#define WATERTIGHT					18 //all 20 humans and their partial models (total of 40)
#define CENTAUR_DAVID				19
#define SEAHORSE_HORSE				20
#define SEAHORSE_SHARK				21
#define GORILLA_GIRL				22
#define SCAPE						23
#define HORSE_TOSCA					24 //high-reso horses from tosca dataset
#define DOG_TOSCA					25 //high-reso dogs from tosca dataset

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,
		 * coords2, //original coords before laplacian smoothing (in case i want to roll back to originals later)
		 * regulator; //displacement for laplacian smoothing

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

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

	//color of this matched vert's sphere, for consistency even in non 1-1 case
	float sphereColorR, sphereColorG, sphereColorB;

	//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; gives clues for extremities (higher mu = more extreme, e.g. at finger tip)
	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;
	//dSpanningUnn spans whole mesh so that [9].dSpanningUnn[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 never normalized
	float* dSpanningUnn, //min/maxGeoDistG here in partial matching case not compatible b/w mesh1 and mesh2 (e.g. long-tail dog has a higher maxGeoDistG than a short-tail one), so no normalized dSpanning[]
		 * dSpanningUnnSS; //properly scaled version of sMesh.dSpanningUnn to be used in matchBasesViaVoting2(); will be inited for those who has dSpanningUnn's defined, i.e. only for the extreme ~10 bases
//	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

	//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, rbase; //is this a baseCandidate, a baseVert, a regional base vert?

	//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
				  baseNghbsK, //inited only for baseVerts; list of k closest base neighbs of this base vert; sorted list s.t. [0] is the closest, [1] next closest, ..
				 * 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

	int matchIdx, //matchIdx = 7 means this mesh1 vert matches w/ mesh2.vert7
		matchIdx2; //adaptiveTraversalList needs matchIdx values around a given base vert; matchIdx gets -1 when it's excluded, but the most recent non -1 value is kept in matchIdx2 to be used
	
	vector< float > matchIdxs, //saves matches for v, as 2 consecutive elements matchIdx & matchCost
					matchIdxs2, //same as matchIdxs but this is for bMesh.regionBaseVerts to be blended into 1-to-1 later
					matchIdxs3, //saves matches for bMesh.v, as 3 consecutive elements matchIdx (in sMesh) & distortion & # of votes, and decides the best later (for matchBasesViaVoting(), 2 consecutive elements per match: matchIdx & # of votes)
					matchIdxs3FirstJump; //good match candidates for this bMesh.v that falls to the left of the first significant jump, i.e. the most confident ones; same format of matchIdxs3

	float matchCostI, //isometry, i.e. geodesic consistency, cost
		  matchCostG, //ground-truth geodesic distance cost
		  matchCostDense, //force1to1 cases in make1to1() uses this costs for correct operation
		  confidence; //confidence of this vert v to be added in entry v-match(v) of vote/confidence mtrx; defined for bMesh only
	//float avgBlossomCost, minBlossomCost, maxBlossomCost; //info about blossom costs[][] entries used in extendToDense2()

	bool processed,
		 hasBadIsometryCost, //bottom 9 of coarse result file is regarded as bad
		 topN, //topN matches w.r.t. isometryCost (= matchCostI)
		 topCurv10, //top10 bases w.r.t. curvature (nonbases not made topCurv10 even if they have higher curvature than bases that are made topCurv10=true)
		 additionalBase, //true if this base is added via computeAdditionalBases() to improve matchBasesViaNghbs()
		 maxMaximaMu, minMinimaMu, //is this a maximaMu vert w/ the max mu value; and similarly for minimaMu vert w/ the min mu value
		 disconnected; //seahorse input has disconnected verts; the ones in the small disconnected are marked to be ignored thorugh all algorithm

	int	 lMaxMinMu, //1 if this base is a local maxima of its 1-ring mu values, 2 if it is local minima, anything else if this is neither maxima nor minima, i.e. non-interesting vert that does not make into baseVerts to be pruned
		 lMaxMinCurv; //1 if this base is a local maxima of its 1-ring curvature values, 2 if it is local minima, anything else if this is neither maxima nor minima, i.e. non-interesting vert that does not make into baseVerts to be pruned

	int idxInBaseVerts, //-1 for non-baseVerts; used in painting to fetch correct diffuseColor
		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 patchArea, //the approximated area of the patch represented by this base vert
		  gc; //Gaussian curvature as another [descriptor]
	int patchIdx; //call each patchArea a patch, and this is just to color each patch in different color

	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; //geo dists from this regional base to all regional bases in the ordering suggested by blossom mapping matchList
	float distortion7, distortion2; //distortion7 tends to be in [0,1], and distortion2 be in [0,0.36], neither can exceed [0,1] interval; i'll not normalize them further 'cos normalizing 3 distortions of 3 generating extremes in [0,1] not good 'cos even if all 3 are good, i'd force one of them to be 1, which is not desired
	
	Vertex(int i, float* c) : idx(i), coords(c), coords2(NULL), baseCandidate(true), base(false), rbase(false), Mu(0.0f), safe(false), worstIin2(-1), worstGin2(-1), additionalBase(false), lMaxMinMu(-1), lMaxMinCurv(-1), disconnected(false),
							   matchIdx(-1), matchCostI(INF), matchCostG(INF), idxInBaseVerts(-1), dSpanningUnn(NULL), dSpanningUnnSS(NULL), baseNghbhood(NULL), topN(false), topCurv10(false), maxMaximaMu(false), minMinimaMu(false) {};

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

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, //sampled base verts on the surface
				  baseVerts2; //copy of the original fps base verts
	vector< int > safeCorresp, //safe correspondences, i.e. desired output, will be added to this list in the format of [i] ~ mesh1.v and [i+1] ~ corresp(mesh1.v)
				  worsts, //worst of each blossom in partialMatching3() is accumulated in this list
				  potentialTriplets, potentialQuadruples, potentialQuintuples, //all potential triplets, quadruples, and quintuples of matches (potentially/quite probably) from the same part
				  denseMap; //coarse map b/w shape extremes can be extended to a dense one; same format in safeCorresp, i.e. 2 consecutive items make 1 matching pair, left from mesh1, right from mesh2
	vector< vector< int > > candidateMaps; //candidate maps to be tracked, i.e. to be dense matched; typically 5 maps stored here, each of size 5

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

	//this is mesh of which jumping man frame?
	int frameID;
	char name[100];

	//is this the target/last mesh of whole frame sequence, e.g. for 2-frame seq, true for 2nd mesh
	bool targetMesh,
		 adaptiveTraversalList, //if true, each entry of blossom mtrx will be formed due to a specific traversalList of nearby bases [never true 'cos u never keep it local while filling all columns in a row, e.g. source.foot to target.hand entry is local only for source.foot part]
		 laplSmoothing, //smooth out the mesh before computing local maxima/minima mu values		 
		 dense40, //use 40 regional bases for dense matching used in confidence computations; if false, 120 bases used (observed no execution time difference b/w 40 and 120 bases on GORILLA_GIRL runs which has low-reso and high-reso meshes)		 
		 batchRun, //output file name convention differs when batch runs are being made
		 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
		 scaleInvDemo; //true when SCALE_INV_DISTORTION_DEMO, i.e. the easier isometric part matching solution

	bool fpsOnMaximaMinima, //true makes fps algo select from the pool of baseVerts which have local maxima/minima mu verts only; false makes fps algo select from the pool of all verts, hence even better coverage but likely to miss salient maxima/minima points at tips/centers
		 distortionByNormalizedGeodesics, //for computation of bestCorresp, true normalizes geodesics by local max geodesic distance and use the second-order distortion (cvpr'10, pami, sgp distortion); false applies scale-inveariant distortion directly on unnormalized geodesics
		 trackCombinatorialMaps, //true to perform dense map on 1+ combinatorial (5x5) maps and decide the best 5x5 by this tracking; false to simple use the min-distortion 5x5 as the combinatorial map
		 noScaleAdjustment; //meshes are already at ~same scale in input; so no scale normalization if true

	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?
		pruningDistDenom; //user provides this in [8,20]
	
	//overall isometry and ground-truth cost of (whole) correspondence
	float isometryCost, groundTruthCost, maxCostIg, maxCostGg,
		  minGeoDistG, maxGeoDistG, //min/max geo dist around the whole mesh
		  MIN_TRIANGLE_AREA_ALLOWED, //no triangle w/ area below this value exists on this mesh
		  radius, //radius yielding ~50 bases on this mesh
		  avgMu, //average of all mu, integral geodesic distance function, values
		  weight7, //weight distortion7-based confidence w/ this value in [0,1] and distortion2-based w/ 1 - weight7
		  xZeta; //multiplier of jump factor in outlier removal

	int* gcOrder, //for nonarbitrary next base selection, i.e. high curvatures first
	   * muOrder; //for nonarbitrary next base selection, i.e. high mu values first
	
	Mesh(bool trg) : targetMesh(trg), MIN_TRIANGLE_AREA_ALLOWED(0.001f), perms2(NULL), worstMatchI(-1), worstMatchG(-1), muOrder(NULL), gcOrder(NULL), laplSmoothing(false), dense40(false), weight7(0.0f), batchRun(false), shortcutEdges(true),
				     scaleInvDemo(false), fpsOnMaximaMinima(false), distortionByNormalizedGeodesics(false), noScaleAdjustment(false), trackCombinatorialMaps(true) {};

	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();

	void computeCorrespondenceErrorMobius(Mesh* mesh2);

	void partialMatching(Mesh* mesh2, float arbitraryScale, int approach);
	void partialMatching2(Mesh* mesh2);
	void partialMatching3(Mesh* mesh2, bool extendCoarseMapToDense = false);
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;

	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;

	//normalizers to be applied to distortions in conf = w7*(1-distortion7) + (1-w7)*(1-distortion2) formulation for compatible weighting, defined on bMesh verts
//	float maxDistortion7, maxDistortion2, minDistortion7, minDistortion2; normalizing 3 distortions of 3 generating extremes in [0,1] not good 'cos even if all 3 are good, i'd force one of them to be 1, which is not desired

	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 scaleBy(float s);
	void prepareForGeodesicComputation();

	int computeBaseVerts(bool COESforDenseSampling = false, bool localMaxMinMu = false, bool localMaxMinCurv = false);
	void computeBaseVerts(float r, bool additionalBases = false, bool print = true);
	vector< int > computeRegionBaseVerts(float r, vector< int > regionVerts, int denseDistoComputation);
	vector< int > computeRegionBaseVerts2(float r, vector< int > regionVerts);
	int computeAdditionalBases();
	void pruneBases(int nRing);
	void pruneBases2(float dist, bool populateTill10 = true);
	int getBaseVertIdx(bool mu, bool COES);
	inline float computePatchArea(int bv, float r);
	inline float computeRegionArea(vector< int > regionVertsB);

	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 sortCurvatures();
	void sortMus();

	void sortSafeCorresp1to1(Mesh* mesh2, char* fName, int n, float avgErr);

	inline float getCorrespListDistortion(Mesh* mesh2, vector< int > tmpSafeCorrespTry, float maxDist1 = 1.0f, float maxDist2 = 1.0f);
	inline float getCorrespListDistortionSS(Mesh* bMesh, vector< int > tmpSafeCorrespTry, float maxDistS, float maxDistB, bool avg, bool sameConfToAll, float w25 = -1.0f);
	inline float getIndividualDistortion(Mesh* mesh2, int si, int bj, vector< int > tmpSafeCorrespTry, float maxDist1 = 1.0f, float maxDist2 = 1.0f);
	inline float getCorrespListScaleInvDistortion(Mesh* mesh2, vector< int > tmpSafeCorrespTry);
	inline float getIndividualScaleInvDistortion(Mesh* bMesh, int si, int bj, vector< int > traversalListTry);
	
	void fillPermutationMtrxs();
	
	void geodesicsBetweenBasesToFile(float r, int finalLevel);

	float isometricPartMatching(Mesh* mesh2);
	void computeBaseNeighbs();
	void computeBaseNeighbsK(int n);
	int trackCandidateMaps(Mesh* mesh2);

	void initial3Matches(Mesh* mesh2);
	void computePotentialTriplets(Mesh* mesh2, bool bestTriplet = false);
	void initial4Matches(Mesh* mesh2);
	void computePotentialQuadruples(Mesh* mesh2, bool bestQuadruple = false);
	void computePotentialQuintuples(Mesh* mesh2, bool bestQuintuple = false);
	void computeAllPotentialsIncremental(Mesh* mesh2, int bestTuple = -1);
	void initial2Matches(Mesh* mesh2); //in use
	void initial3MatchesV2(Mesh* mesh2);
	void initial50Matches(Mesh* mesh2); //in use
	float blossomOnBases(Mesh* mesh2, vector< int > traversalListTry);
	float blossomOnBases2(Mesh* bMesh, Mesh* sMesh);
	int excludeFromBases(int v);
	inline int leastConfidentExtreme(int e1, int e2, int e3);

	void matchBasesViaAll(Mesh* bMesh, Mesh* sMesh);
	void matchBasesViaAllEff(Mesh* bMesh, Mesh* sMesh);	
	void matchBasesViaNghbs(Mesh* bMesh, Mesh* sMesh);
	bool matchBasesViaVoting(Mesh* bMesh, Mesh* sMesh, char* fName);
	bool matchBasesViaVoting2(Mesh* bMesh, Mesh* sMesh, char* fName, bool oneToOneGreedy);
	void extendToDense(Mesh* sMesh, bool bMeshIsThis);
	void extendToDense2(Mesh* sMesh, bool bMeshIsThis);

	void laplacianSmoothing();
};

#endif
