#include "Mesh.h"

Mesh::Mesh(int fid, char* fileFormat, int mIdx)
{
	theInput = SCAPE; //not in use in this public release
	ordering = 0; //0 for the original ordering of the mesh files, not in use in this public release

	frameID = fid;
	shortcutEdges = false;//due to incredibly weird occasional crash in addShortcutEdge() for x64 release runs, e.g. for scape 2.ply; so never use shortcuts
	affinity = NULL;
	char dName[250], fName[250], orderingName[250] = "";
//note that dName and orderingName are unused in this public release
	if (theInput == SCAPE)
		sprintf(dName, "SCAPE\\");
	else if (theInput == TOSCA_CAT)
		sprintf(dName, "TOSCA\\cat\\");
	else if (theInput == TOSCA_CENTAUR)
		sprintf(dName, "TOSCA\\centaur\\");
	else if (theInput == TOSCA_DAVID)
		sprintf(dName, "TOSCA\\david\\");
	else if (theInput == TOSCA_DOG)
		sprintf(dName, "TOSCA\\dog\\");
	else if (theInput == TOSCA_GORILLA)
		sprintf(dName, "TOSCA\\gorilla\\");
	else if (theInput == TOSCA_HORSE)
		sprintf(dName, "TOSCA\\horse\\");
	else if (theInput == TOSCA_VICTORIA)
		sprintf(dName, "TOSCA\\victoria\\");
	else if (theInput == TOSCA_WOLF)
		sprintf(dName, "TOSCA\\wolf\\");
	else if (theInput == TOSCA_DAVIDVICTORIA)
		sprintf(dName, "TOSCA\\david-victoria\\");
	if (ordering > 0)
		sprintf(orderingName, "another order%d\\", ordering);

	//sprintf(fName, "input\\%s%s%d.%s", dName, orderingName, frameID, fileFormat);
	sprintf(fName, "input\\%d.%s", frameID, fileFormat); //public release format

	if (strcmp(fileFormat, "obj") == 0)
		loadObj(fName);
	else if (strcmp(fileFormat, "off") == 0)
		loadOff(fName);
	else if (strcmp(fileFormat, "ply") == 0)
		loadPly(fName);	

	//ground-truth matches are generally in i-i format; but for TOSCA_DAVIDVICTORIA i set them manually and fread from file; for consistency in code i use grdTruthOnMesh1 for i-i as well for inputs != TOSCA_DAVIDVICTORIA
	if (theInput != TOSCA_DAVIDVICTORIA)
		for (int v = 0; v < (int) verts.size(); v++)
			verts[v]->grdTruthOnMesh1 = v; //for i-i case
	else //if (theInput == TOSCA_DAVIDVICTORIA)
	{
		FILE* dvPtr, * vdPtr;
		if (! (dvPtr = fopen("input\\TOSCA\\david-victoria ground-truth.dat", "r")) || ! (vdPtr = fopen("input\\TOSCA\\victoria-david ground-truth.dat", "r")))
		{
			cout << "manual ground-truth correspondence file(s) not found (required for evaluation only)\n";
			exit(0);
		}
		//even mesh idxs (mIdx) for david, odd for victoria; so, if this is a david mesh (mIdx is even) then find ground-truth matches on victoria mesh for each david extreme via vdPtr; similarly for victoria mesh use dvPtr
		fscanf(dvPtr, "david's grd-truth match (left) for victoria0-1-2-3-5-6-7 extreme (right):\n");
		fscanf(vdPtr, "victoria's grd-truth match (left) for david0-1-2-3-4-5-6 extreme (right):\n");
		int vicIdx, davIdx;
		if (mIdx % 2 == 0) //current mesh is david
			while (! feof(vdPtr))
			{
				fscanf(vdPtr, "%d\t%d\n", &vicIdx, &davIdx);
				verts[davIdx]->grdTruthOnMesh1 = vicIdx; //here mesh1 represents victoria
			}
		else //current mesh is victoria
			while (! feof(dvPtr))
			{				
				fscanf(dvPtr, "%d\t%d\n", &davIdx, &vicIdx);
				verts[vicIdx]->grdTruthOnMesh1 = davIdx; //here mesh1 represents david
			}
		fclose(vdPtr);
		fclose(dvPtr);
	}	

	updateNormals();
	//more accurate shortest paths
	prepareForGeodesicComputation();
	cout << "Mesh w/ ID = " << frameID << " in collection has " << (int) tris.size() << " tris, " << (int) verts.size() << " verts, " << (int) edges.size() << " edges, " << (int) sedges.size() << " shortcut edges\n\n";
}

void Mesh::loadObj(char* meshFile)
{
	cout << "Mesh initializing (to model in " << meshFile << " file)..\n";

	FILE* fPtr = fopen(meshFile, "r");
	if (fPtr == 0)
	{
		cout << "cannot read " << meshFile << endl;
		exit(0);
	}

	char type;	//type of line: vertex line (v) or face line (f)
	float a, b, c;	//for face lines casting these to int will suffice
	bool first = true;
	float minX, minY, minZ, maxX, maxY, maxZ; //bbox to be used in center-of-bbox calculation
	minEdgeLen = INF;	maxEdgeLen = -INF; edgeLenTotal = 0.0f;
	while (fscanf(fPtr, "%c %f %f %f\n", &type, &a, &b, &c) != EOF) //go till the end of file
		if (type == 'v')
		{
			//in original .obj file, tri.areas ~= .1f and, e.g., split of them deceives >= .01f assertion in addTria wrongly
			a *= SCALE_FACTOR;	b *= SCALE_FACTOR;	c *= SCALE_FACTOR;

			//none of verts are removed/duplicated initially 'cos i read .obj
			float* coords = new float[3];
			coords[0] = a;
			coords[1] = b;
			coords[2] = c;
			addVertexND(coords); //ND: no duplicate check

			//bbox computation
			if ((int) verts.size() == 1)
			{
				minX = maxX = a;	minY = maxY = b;	minZ = maxZ = c;
			}
			else
			{
				if (a < minX)	minX = a; else if (a > maxX)	maxX = a;
				if (b < minY)	minY = b; else if (b > maxY)	maxY = b;
				if (c < minZ)	minZ = c; else if (c > maxZ)	maxZ = c;
			}
		}
		else //if (type == 'f')
			addTriangle((int) a - 1, (int) b - 1, (int) c - 1); //-1: obj indices start at 1

	float cobX = (minX+maxX) / 2.0f, cobY = (minY+maxY) / 2.0f, cobZ = (minZ+maxZ) / 2.0f;
	cout << "bbox.center: (" << cobX << ", " << cobY << ", " << cobZ << "); "
		 << "bbox.(w,h,d): (" << maxX-minX << ", " << maxY-minY << ", " << maxZ-minZ << ")\n"
		 << "min <= avg <= maxEdgeLen: " << minEdgeLen << " <= " << edgeLenTotal/((int) edges.size()) << " <= " << maxEdgeLen << "\n";
	avgEdgeLen = edgeLenTotal/((int) edges.size());
	fclose(fPtr);
}
void Mesh::loadOff(char* meshFile)
{
	cout << "Mesh initializing (to model in " << meshFile << " file)..\n";

	FILE* fPtr = fopen(meshFile, "r");
	if (fPtr == 0)
	{
		cout << "cannot read " << meshFile << endl;
		exit(0);
	}

	char off[25];
	fscanf(fPtr, "%s\n", &off); //cout << off << " type file\n";

	float a, b, c, d;	//for face lines and the 2nd line (that gives # of verts, faces, and edges) casting these to int will suffice
	fscanf(fPtr, "%f %f %f\n", &a, &b, &c);
	int nVerts = (int) a, v = 0;

	float minX, minY, minZ, maxX, maxY, maxZ; //bbox to be used in center-of-bbox calculation
	minEdgeLen = INF;
	maxEdgeLen = -INF;
	edgeLenTotal = 0.0f;
	while (v++ < nVerts) //go till the end of verts coord section
	{
		fscanf(fPtr, "%f %f %f\n", &a, &b, &c);

		//in original .off file, tri.areas ~= .1f and, e.g., split of them deceives >= .01f assertion in addTria wrongly
		a *= SCALE_FACTOR;	b *= SCALE_FACTOR;	c *= SCALE_FACTOR;

		//none of verts are removed/duplicated initially 'cos i read .off
		float* coords = new float[3];
		coords[0] = a, coords[1] = b, coords[2] = c;
		addVertexND(coords); //ND: no duplicate check

		//bbox computation
		if ((int) verts.size() == 1)
		{
			minX = maxX = a;	minY = maxY = b;	minZ = maxZ = c;
		}
		else
		{
			if (a < minX)	minX = a; else if (a > maxX)	maxX = a;
			if (b < minY)	minY = b; else if (b > maxY)	maxY = b;
			if (c < minZ)	minZ = c; else if (c > maxZ)	maxZ = c;
		}
	}

	//verts ready, time to fill triangles
	while (fscanf(fPtr, "%f %f %f %f\n", &d, &a, &b, &c) != EOF) //go till the end of file
		addTriangle((int) a, (int) b, (int) c); //no -1 'cos idxs start from 0 for off files

	float cobX = (minX+maxX) / 2.0f, cobY = (minY+maxY) / 2.0f, cobZ = (minZ+maxZ) / 2.0f;
	cout << "bbox.center: (" << cobX << ", " << cobY << ", " << cobZ << "); "
		 << "bbox.(w,h,d): (" << maxX-minX << ", " << maxY-minY << ", " << maxZ-minZ << ")\n"
		 << "min <= avg <= maxEdgeLen: " << minEdgeLen << " <= " << edgeLenTotal/((int) edges.size()) << " <= " << maxEdgeLen << "\n";
	avgEdgeLen = edgeLenTotal/((int) edges.size());
	fclose(fPtr);
}
void Mesh::loadPly(char* meshFile)
{
	cout << "Mesh initializing (to model in " << meshFile << " file)...\n";

	FILE* fPtr = fopen(meshFile, "r");
	if (fPtr == 0)
	{
		cout << "cannot read " << meshFile << endl;
		exit(0);
	}

	char ply[25], nVertsStr[25];
	fscanf(fPtr, "%s\n", &ply); //skip "ply" at line 1
	fscanf(fPtr, "%s %s %s\n", &ply, &ply, &ply); //skip "format ascii 1.0" at next line
	fscanf(fPtr, "%s %s %s\n", &ply, &ply, &nVertsStr); //skip "element vertex 12500" at next line
	int nVerts = atoi(nVertsStr), v = 0;
	fscanf(fPtr, "%s %s %s\n", &ply, &ply, &ply); //skip "property float x" at next line
	fscanf(fPtr, "%s %s %s\n", &ply, &ply, &ply); //skip "property float y" at next line
	fscanf(fPtr, "%s %s %s\n", &ply, &ply, &ply); //skip "property float z" at next line
	fscanf(fPtr, "%s %s %s\n", &ply, &ply, &ply); //skip "element face 25000" at next line
	fscanf(fPtr, "%s %s %s %s %s\n", &ply, &ply, &ply, &ply, &ply); //skip "property list uchar int vertex_indices" at next line
	fscanf(fPtr, "%s\n", &ply); //skip "end_header" at next line

	float a, b, c, d;	//for face lines and the 2nd line (that gives # of verts, faces, and edges) casting these to int will suffice
	float minX, minY, minZ, maxX, maxY, maxZ; //bbox to be used in center-of-bbox calculation
	minEdgeLen = INF;
	maxEdgeLen = -INF;
	edgeLenTotal = 0.0f;
	while (v++ < nVerts) //go till the end of file
	{
		fscanf(fPtr, "%f %f %f\n", &a, &b, &c);

		//in original .off file, tri.areas ~= .1f and, e.g., split of them deceives >= .01f assertion in addTria wrongly
		a *= SCALE_FACTOR;	b *= SCALE_FACTOR;	c *= SCALE_FACTOR;

		//none of verts are removed/duplicated initially 'cos i read .ply
		float* coords = new float[3];
		coords[0] = a, coords[1] = b, coords[2] = c;
		addVertexND(coords); //ND: no duplicate check

		//bbox computation
		if ((int) verts.size() == 1)
		{
			minX = maxX = a;	minY = maxY = b;	minZ = maxZ = c;
		}
		else
		{
			if (a < minX)	minX = a; else if (a > maxX)	maxX = a;
			if (b < minY)	minY = b; else if (b > maxY)	maxY = b;
			if (c < minZ)	minZ = c; else if (c > maxZ)	maxZ = c;
		}
	}

	//verts ready, time to fill triangles
	while (fscanf(fPtr, "%f %f %f %f\n", &d, &a, &b, &c) != EOF) //go till the end of file
		addTriangle((int) a, (int) b, (int) c);

	float cobX = (minX+maxX) / 2.0f, cobY = (minY+maxY) / 2.0f, cobZ = (minZ+maxZ) / 2.0f;
	cout << "bbox.center: (" << cobX << ", " << cobY << ", " << cobZ << "); "
		 << "bbox.(w,h,d): (" << maxX-minX << ", " << maxY-minY << ", " << maxZ-minZ << ")\n"
		 << "min <= avg <= maxEdgeLen: " << minEdgeLen << " <= " << edgeLenTotal/((int) edges.size()) << " <= " << maxEdgeLen << "\n";
	fclose(fPtr);
}

/*int Mesh::addVertex(float* v)
{
	//add v to verts, if not duplicated

	bool duplicateExists = false;
	int i = (int) verts.size();

	//duplicate prevention
	for (i = (int) verts.size() - 1; i >= 0; i--)
	{
		if (fabs(verts[i]->coords[0] - v[0]) < ROUNDING_ERR && //since i scale stuff initially w/ SCALE_FACTOR = 1000 just in case, ROUNDING_ERR = 0.01f should be good here
			fabs(verts[i]->coords[1] - v[1]) < ROUNDING_ERR &&
			fabs(verts[i]->coords[2] - v[2]) < ROUNDING_ERR)
		{
			duplicateExists = true;
			break;	//more than efficiency: return the correct i by this break
		}
	}
	if (! duplicateExists)
	{
		int idx = (int) verts.size();
		verts.push_back(new Vertex(idx, v));
		return idx; //pointer to this new item (last index)
	}

	//this duplicate vert v is ith element of verts
	return i;
}*/
void Mesh::addVertexND(float* c)
{
	//fastly add verts w/ coords c to mesh; ND: no duplication check

	int vSize = (int) verts.size(); //size before push_back just-below
	verts.push_back(new Vertex(vSize, c));
}

void Mesh::addTriangle(int v1i, int v2i, int v3i)
{
	//add tri v1i-v2i-v3i to the triInd'th element of tris

	int triInd = (int) tris.size();

	tris.push_back( new Triangle(triInd, v1i, v2i, v3i) );

	//update each vertex with its new tri and vert neighbors
	verts[v1i]->triNeighbors.push_back(triInd);
	//addVertNeighbor returns true if verts[v1i] is not neighbor of [v2i]; so they can form a new edge
	if (verts[v1i]->addVertNeighbor(v2i))	addEdge(v1i, v2i);
	if (verts[v1i]->addVertNeighbor(v3i))	addEdge(v1i, v3i);
	//in other words, a 'false' returned by addVertNeighbor prevents adding the same edge if it already exists

	verts[v2i]->triNeighbors.push_back(triInd);
	verts[v2i]->addVertNeighbor(v1i);	//v1i v2i edge is (probably) added above; do not add v2i v1i edge!!!
	if (verts[v2i]->addVertNeighbor(v3i))	addEdge(v2i, v3i);	//v3i v2i is not added above, so add this if ok

	verts[v3i]->triNeighbors.push_back(triInd);
	verts[v3i]->addVertNeighbor(v1i);	//no addEdge risk, so no if control
	verts[v3i]->addVertNeighbor(v2i);
}

int Mesh::getEdgeSharedBy(int vi, int vj)
{
	//return idx of edge[vi, vj]

	for (int vie = 0; vie < (int) verts[vi]->edgeList.size(); vie++)
		for (int vje = 0; vje < (int) verts[vj]->edgeList.size(); vje++)
			if (verts[vj]->edgeList[vje] == verts[vi]->edgeList[vie])
				return verts[vj]->edgeList[vje];

	//edge not exist
	return -99;
}

void Mesh::updateNormals()
{
	//update normals of tris and verts in mesh

	//triangles
	for (int t = 0; t < (int) tris.size(); t++)
	{
		int v1i = tris[t]->v1i, v2i = tris[t]->v2i, v3i = tris[t]->v3i;

		tris[t]->setNormal(verts[v1i]->coords, verts[v2i]->coords, verts[v3i]->coords);
		tris[t]->setArea(verts[v1i]->coords, verts[v2i]->coords, verts[v3i]->coords);
	}

	//vertices
	for (int v = 0; v < (int) verts.size(); v++)
	{
		float* normalTotal = new float[3];
		normalTotal[0] = normalTotal[1] = normalTotal[2] = 0.0f;
		int howManyTriNeighb = 0;
		for (int tneighb = 0; tneighb < (int) verts[v]->triNeighbors.size(); tneighb++)		
			if (tris[ verts[v]->triNeighbors[tneighb] ]->area >= 0.001f/*MIN_TRIANGLE_AREA_ALLOWED*/)	//ignore temporary 0-area tris which occur by i.e. split
			{
				for (int i = 0; i < 3; i++)
					normalTotal[i] += tris[ verts[v]->triNeighbors[tneighb] ]->normal[i]; //maby a better way is to weight tri-normals w/ their areas for v.normal computation
				howManyTriNeighb++;
			}
		if (howManyTriNeighb != 0) //if 0, go on w/ the ex-normal
		{
			for (int i = 0; i < 3; i++)
				normalTotal[i] /= (float) howManyTriNeighb;
			verts[v]->setNormal(normalTotal);
		}
	}
}

inline float distanceBetween(float* v1, float* v2)
{
	//In fact, the keyword inline is not necessary. If the function is defined with its body directly and the function
	//has a smaller block of code, it will be automatically treated as inline by the compiler
	return (float) sqrt( (v2[0]-v1[0])*(v2[0]-v1[0]) + (v2[1]-v1[1])*(v2[1]-v1[1]) + (v2[2]-v1[2])*(v2[2]-v1[2]));
}
inline float distanceBetween2(float* v1, float* v2)
{
	//squared distance which skips sqrt part of distanceBetween() for efficiency; e.g. use this for comparison-based distances
	return (float) (v2[0]-v1[0])*(v2[0]-v1[0]) +  (v2[1]-v1[1])*(v2[1]-v1[1]) + (v2[2]-v1[2])*(v2[2]-v1[2]);
}

inline float Mesh::calcTriArea(int t)
{
	float* v1 = verts[ tris[t]->v1i ]->coords,
		 * v2 = verts[ tris[t]->v2i ]->coords,
		 * v3 = verts[ tris[t]->v3i ]->coords;

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

	//When it comes time to de-allocate them, I had them all in a single delete[] statement: delete[] a, b, c, d, e; this appears to only be de-allocating a and completely ignoring b, c, d, and e
	delete p_r;
	delete q_r;
	delete cross; //arrays allocated with new[] must be deallocated with delete[] so i should've added [] here but not a big deal

	return 0.5f * length;
}

void Mesh::addEdge(int v1i, int v2i)
{
	int initial = (int) edges.size();
	Edge* newcomer = new Edge(initial, v1i, v2i, distanceBetween(verts[v1i]->coords, verts[v2i]->coords));
	edgeMultiSet.insert(MultiSetContent(newcomer->idx, newcomer->length));

	//update the list of edges of which verts[v1i] is a member
	verts[v1i]->edgeList.push_back(newcomer->idx);
	verts[v2i]->edgeList.push_back(newcomer->idx); //do same thing for [v2i]

	edges.push_back(newcomer);

	//information
	if (newcomer->length < minEdgeLen)
		minEdgeLen = newcomer->length;
	if (newcomer->length > maxEdgeLen)
		maxEdgeLen = newcomer->length;
	edgeLenTotal += newcomer->length;
}

void Mesh::prepareForGeodesicComputation()
{
	//prepares the mesh to upcoming geodesic computation process

	//1) fine distribution of verts needed; so, use split to resample verts until all edge lengths < threshold
//	edgeSplit();
//	updateNormals(); //newcomer tris/verts should learn their normals

	//2) edges should be isotropic, i.e. uniform in all directions, for good geodesic approximation
	genShortcutEdges();
}

bool Mesh::vertsCollinear(float* v1, float* v2, float* v3)
{
	//3 points may be on the same line which makes area of tri 0; for convexity check, 4th vert
	//projected onto the plane (no plane actually) of such 3 collinear verts always yield a convex shape

	// 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; //arrays allocated with new[] must be deallocated with delete[] so i should've added [] here but not a big deal

	return (area == 0.0f); //area ~ 0 means that verts are collinear

/*	if (theInput != JUMPING && theInput != DANCING && area == 0.0f)
		return true;
	else if (theInput != JUMPING && theInput != DANCING && area != 0.0f)
		return false;
	else //not hybrid run, so hybridModelsToSameCoorFrame() not scaled coords and triangle areas not violated
		return (area < MIN_TRIANGLE_AREA_ALLOWED ? true : false);	//.01 margin for variable precision*/
}

float* Mesh::projectVertOntoPlane(int toBeProjected, float* planeUnormal, float distOfPlaneFromOrigin)
{
	//using the projectiveMatrix from www.euclideanspace.com/maths/geometry/elements/plane/

	float* prjctdVert = new float[3], * v = verts[toBeProjected]->coords;

	//using the projectiveMatrix from www.euclideanspace.com/maths/geometry/elements/plane/
	prjctdVert[0] = (planeUnormal[1]*planeUnormal[1] + planeUnormal[2]*planeUnormal[2])*v[0] - planeUnormal[0]*planeUnormal[1]*v[1] - planeUnormal[0]*planeUnormal[2]*v[2];
	prjctdVert[1] = -planeUnormal[1]*planeUnormal[0]*v[0] + (planeUnormal[0]*planeUnormal[0] + planeUnormal[2]*planeUnormal[2])*v[1] - planeUnormal[1]*planeUnormal[2]*v[2];
	prjctdVert[2] = -planeUnormal[2]*planeUnormal[0]*v[0] - planeUnormal[2]*planeUnormal[1]*v[1] + (planeUnormal[0]*planeUnormal[0] + planeUnormal[1]*planeUnormal[1])*v[2];

	//this crucial part is added after my visualization; thanx to my test friendly style!!!
	//result just above is projected onto plane passing through origin; move this point to the
	//actual plane by moving it along the plane planeUnormal direction with an amount of D
	for (int i = 0; i < 3; i++)
		prjctdVert[i] += distOfPlaneFromOrigin*planeUnormal[i];

	return prjctdVert;
}

bool lineSegmentsIntersect(float* s1p1, float* s1p2, float* s2p1, float* s2p2)
{
	//returns true if segment1 (s1) defined by 2 pnts (p1, p2) intersects segment2; z-coords are ignored, having the
	//projection effect on x-y plane (safe 'cos s1 & s2 are guaranteed to be on same plane, not necessarily x-y tough)

	//algo from local.wasp.uwa.edu.au/~pbourke/geometry/lineline2d/
	float denom = (s2p2[1]-s2p1[1])*(s1p2[0]-s1p1[0]) - (s2p2[0]-s2p1[0])*(s1p2[1]-s1p1[1]);
	if (denom == 0.0f)
		return false; //parallel lines
	float numer1 = (s2p2[0]-s2p1[0])*(s1p1[1]-s2p1[1]) - (s2p2[1]-s2p1[1])*(s1p1[0]-s2p1[0]),
		  numer2 = (s1p2[0]-s1p1[0])*(s1p1[1]-s2p1[1]) - (s1p2[1]-s1p1[1])*(s1p1[0]-s2p1[0]),
		  u1 = numer1 / denom, u2 = numer2 / denom;
	return (u1 >= 0.0f && u1 <= 1.0f		&&		u2 >= 0.0f && u2 <= 1.0f);
}

void Mesh::genShortcutEdges()
{
	//introduces additional edges so that resulting mesh has edges that are uniform in all directions

	//shortcut sometimes crash below on x64 release mode only; better not use them at all
	if (! shortcutEdges || verts.size() > 40000) //40K+ are dense meshes anyway; so, need for path smoothings
		return; //no shortcut 'cos it slows down the base computation as well as dijkstra (but shortcut'ed more accurate; so, do it)

//	cout << "shortcut edges being generated....\n";

	minSedgeLen = INF;
	maxSedgeLen = -INF;
	sedgeLenTotal = 0.0f;
	for (int t = 0; t < (int) tris.size(); t++)
	{
		//get the 3 triangles adjacent to t (edge-adjacency)
		int v1i = tris[t]->v1i, v2i = tris[t]->v2i, v3i = tris[t]->v3i;
		float* v1ic = verts[v1i]->coords, * v2ic = verts[v2i]->coords, * v3ic = verts[v3i]->coords;
		if ( vertsCollinear(v1ic, v2ic, v3ic) )
			continue; //no plane defined for this degenarate tri

		int e1 = getEdgeSharedBy(v1i, v3i), e2 = getEdgeSharedBy(v1i, v2i), e3 = getEdgeSharedBy(v2i, v3i);
		if (e1 < 0 || e2 < 0 || e3 < 0)
			continue; //not expected to happen but anyway
		int* shrdTris = getTrianglesSharedBy(edges[e1]);
		int t1 = (shrdTris[0] == t ? shrdTris[1] : shrdTris[0]);
		shrdTris = getTrianglesSharedBy(edges[e2]);
		int t2 = (shrdTris[0] == t ? shrdTris[1] : shrdTris[0]);
		shrdTris = getTrianglesSharedBy(edges[e3]);
		int t3 = (shrdTris[0] == t ? shrdTris[1] : shrdTris[0]);

		if (t1 == -1 || t2 == -1 || t3 == -1) //border edge on (non-closed) mesh makes one triangle neighbor to an edge; so, the 2nd tri may become -1
			continue;
		int u1, u2, u3; //u1 is the uncommon vert of the 1st adj tri (t1), and so on
		if (! (tris[t1]->v1i == v1i || tris[t1]->v1i == v3i)) //v2i ignored 'cos t1 is created via v1i-v3i edge
			u1 = tris[t1]->v1i;
		else if (! (tris[t1]->v2i == v1i || tris[t1]->v2i == v3i))
			u1 = tris[t1]->v2i;
		else
			u1 = tris[t1]->v3i;

		if (! (tris[t3]->v1i == v3i || tris[t3]->v1i == v2i))
			u2 = tris[t3]->v1i;
		else if (! (tris[t3]->v2i == v3i || tris[t3]->v2i == v2i))
			u2 = tris[t3]->v2i;
		else
			u2 = tris[t3]->v3i;
		
		if (! (tris[t2]->v1i == v2i || tris[t2]->v1i == v1i))
			u3 = tris[t2]->v1i;
		else if (! (tris[t2]->v2i == v2i || tris[t2]->v2i == v1i))
			u3 = tris[t2]->v2i;
		else
			u3 = tris[t2]->v3i;

		//now i've v1i, u1, v2i, u2, v3i, u3; need to project u's to the plane of v1i-v2i-v3i first
		float distOfPlaneFromOrigin = 0; //tris[t]->unormal.dot(v1ic);
		for (int i = 0; i < 3; i++)
			distOfPlaneFromOrigin += tris[t]->unormal[i]*v1ic[i];
		//u1 projected to the plane of t, and projected coord is kept in u1c
		float*  u1c = projectVertOntoPlane(u1, tris[t]->unormal, distOfPlaneFromOrigin),
				//project u2 as u2c
			 *	u2c = projectVertOntoPlane(u2, tris[t]->unormal, distOfPlaneFromOrigin),
				//project u3 as u3c		
			 *	u3c = projectVertOntoPlane(u3, tris[t]->unormal, distOfPlaneFromOrigin);


		//now i've v1i, u1c, v2i, u2c, v3i, u3c; create desired shortcut edges b/w corresponding pairs
		//test 1 corresponding pair: v1i ~ u2c (v1i ~ u2c line segment is related to edge e3 segment)
		if ( lineSegmentsIntersect(v1ic, u2c,	v2ic, v3ic) )
			addShortcutEdge(v1i, u2, v1ic, u2c);
		//test other corresponding pairs in each if-statement
		if ( lineSegmentsIntersect(v2ic, u1c,	v1ic, v3ic) )
			addShortcutEdge(v2i, u1, v2ic, u1c);
		if ( lineSegmentsIntersect(v3ic, u3c,	v1ic, v2ic) )
			addShortcutEdge(v3i, u3, v3ic, u3c);
		//shortcut edges from u1
		if ( lineSegmentsIntersect(u1c,	u3c,	v1ic, v3ic) )
			addShortcutEdge(u1, u3, u1c, u3c);
		if ( lineSegmentsIntersect(u1c,	u2c,	v1ic, v3ic) )
			addShortcutEdge(u1, u2, u1c, u2c);
		//shortcut edges from u2 (note that, u1-u2 already checked above, so no u2-u1 here
		if ( lineSegmentsIntersect(u2c,	u3c,	v2ic, v3ic) )
			addShortcutEdge(u2, u3, u2c, u3c);
	}

//	if (((int) sedges.size()) != 0)
//		cout << "shortcutEdges: min <= avg <= maxSedgeLen: " << minSedgeLen << " <= " << sedgeLenTotal/((int) sedges.size()) << " <= " << maxSedgeLen << "\n";
//	cout << "Done! # of generated shortcut edges: " << (int) sedges.size() << ", for " << (int) verts.size() << " verts\n\n";
}
inline void Mesh::addShortcutEdge(int v1i, int v2i, float* v1, float* v2Prjctd)
{
	//shortcut edge b/w v1i & v2i added if it has not been previously generated

	//check if virtual connection already exists
	for (int se = 0; se < (int) verts[v1i]->sedgeList.size(); se++)
	{
		int sidx = verts[v1i]->sedgeList[se];
		if (sidx < 0 || sidx >= (int) sedges.size())
			return; //paranoia
		Edge* sedge = sedges[sidx];
		int otherVert = (sedge->v1i == v1i ? sedge->v2i : sedge->v1i);
		if (otherVert == v2i) //virtual connection already exists
			return;
	}

	//add the shortcut edge for the first time b/w these 2 verts, where edge length is from projected domain
	float plen = distanceBetween(v1, v2Prjctd); //still crash below w/o fnctn calls (float) sqrt( (v2Prjctd[0]-v1[0])*(v2Prjctd[0]-v1[0]) + (v2Prjctd[1]-v1[1])*(v2Prjctd[1]-v1[1]) + (v2Prjctd[2]-v1[2])*(v2Prjctd[2]-v1[2]));

	int sSize = (int) sedges.size(); //size before push_back just-below
	sedges.push_back(new Edge(sSize, v1i, v2i, plen));
	verts[v1i]->sedgeList.push_back(sSize);

//if (v1i == 511 && v2i == 435) cout << "c3\tvs = " << (int) verts.size() << "\n\n" << v2i << "\n\n" << verts[v2i]->sedgeList.size() << "\n\n" << sSize << "\n\n"; weird: this prints logical values all in bounds

	verts[v2i]->sedgeList.push_back(sSize); //do same thing for [v2i]

//if (v1i == 511 && v2i == 435) cout << "c4" << endl; yet can't reach here in x64 release mode in scape model 2.ply; why is this weird crash???????????

	if (plen < minSedgeLen)
		minSedgeLen = plen;
	if (plen > maxSedgeLen)
		maxSedgeLen = plen;
	sedgeLenTotal += plen;
}


int* Mesh::getTrianglesSharedBy(Edge* bigEdge, int v1i, int v2i)
{
	//tri t is shared by bigEdge iff it exists in the triNeigbors of verts[v1i] & [v2i]
	if (v1i == -1) //learn v1i and v2i from bigEdge
	{
		v1i = bigEdge->v1i;
		v2i = bigEdge->v2i;
	}
	int a = 0, tr1, tr2;
	int* shTris = new int[2];	//portable version (memo-friendly as well 'cos i delete shrdTris in split/flip())
	//seeing this took my 6 hours!!! prevent removed triangles from being selected as shared!
	for (int t1 = 0; t1 < (int) verts[v1i]->triNeighbors.size(); t1++)
	{
		tr1 = verts[v1i]->triNeighbors[t1];
		for (int t2 = 0; t2 < (int) verts[v2i]->triNeighbors.size(); t2++)
		{
			tr2 = verts[v2i]->triNeighbors[t2];
			if  (tris[tr2] != 0 && (tr1 == tr2))
				shTris[a++] = tr1; //or tr2
		}
	}
	//if called from split/flip, a is certainly incremented by 1; but it remains at 1 (no more ++) 
	//if bigEdge is a border edge (not in my case 'cos my mesh is closed); take precautions anyway;
	if (a == 1) shTris[1] = -1;	//number -1 must not change 'cos used by somewhere else
if (shTris[1] == -1) cout << v1i << "," << v2i << " is a border edge found in my closed-mesh :(\n\n"; //exit(0);}

//if called from contract, a may remain 0, which tells me to remove that edgesAdded[cea] neighbor to 2 removed tris
if (a == 0)	shTris[0] = -1;	//since that part of contract() is inactive, don't care here

	return shTris;
}

inline float maxOf(float a, float b) { return (a > b ? a : b); }
inline float minOf(float a, float b) { return (a < b ? a : b); }

void Mesh::updateEdgeLengths()
{
	//called after coords scaling which changes edge lengths change

	edgeMultiSet.clear(); //re-fill this w/ new edge lengths (although order not change after scale, i want to maintain latest length values in edgeMultiSet (for no reason)
	edgeLenTotal = 0.0f;
	minEdgeLen = INF;
	maxEdgeLen = -INF;
	for (int e = 0; e < (int) edges.size(); e++)
	{
		edges[e]->length = distanceBetween(verts[edges[e]->v1i]->coords, verts[edges[e]->v2i]->coords);
		edgeMultiSet.insert(MultiSetContent(edges[e]->idx, edges[e]->length));
		if (edges[e]->length < minEdgeLen)
			minEdgeLen = edges[e]->length;
		if (edges[e]->length > maxEdgeLen)
			maxEdgeLen = edges[e]->length;
		edgeLenTotal += edges[e]->length;
	}
	avgEdgeLen = edgeLenTotal/((int) edges.size());
//	cout << "after updateEdgeLengths(): min <= avg <= maxEdgeLen: " << minEdgeLen << " <= " << avgEdgeLen << " <= " << maxEdgeLen << "\n";
}
void Mesh::updateSedgeLengths()
{
	//similar to updateEdgeLengths() but this updates shortcut edges

	sedgeLenTotal = 0.0f;
	for (int e = 0; e < (int) sedges.size(); e++)
	{
		sedges[e]->length = distanceBetween(verts[sedges[e]->v1i]->coords, verts[sedges[e]->v2i]->coords);
		sedgeLenTotal += sedges[e]->length;
	}
//	cout << "after updateSedgeLengths(): avg shortcut edge len = " << sedgeLenTotal/((int) sedges.size()) << "\n";
}

void Mesh::dijkstraShortestPaths(int sourceVert)//, float far)
{
	//computes shortest paths from sourceVert to all other verts that are at distance far or less

	if (verts[sourceVert]->dSpanningUnn) //don't recompute if already computed
		return;
	verts[sourceVert]->dSpanning = new float[ (int) verts.size() ]; //in [0,1]
	verts[sourceVert]->dSpanningUnn = new float[ (int) verts.size() ]; //in original scale

	//use Fibonacci heap to get O(VlgV + e); minHeap takes O(VlgV + elgV), simple array takes O(V^2)
	FibHeap* h = new FibHeap();

	/////// original Dijkstra to compute shortest paths from sourceVert to all others ///////
	//initialization
	for (int v = 0; v < (int) verts.size(); v++)
	{
		verts[v]->dSpanned = (v == sourceVert ? 0.0f : INF);
		FibHeapNode* hn = new FibHeapNode(verts[v]->dSpanned, v); //v is a handle from headNode to verts,
		verts[v]->heapNode = hn;								  //and hn is an handle from verts to heapNode
		h->Insert(hn);
		//no delete hn; here 'cos it'll be extracted, used, and deleted below
	}

	while (true)
	{
		//extract vert w/ min dSpanned value
//not doing the brute force O(V) algo to extract min (useful for dense graph, i.e. simple array usage, which is not my case)
		FibHeapNode* minNode = h->ExtractMin(); //fast O(lgV) priority queue algo
		if (minNode == NULL)
			break; //no more verts w/ unknown shortest dist; so, we're done
		int minDidx = minNode->element;
		float minD = minNode->key;
		verts[sourceVert]->dSpanningUnn[minDidx] = minD; //remember original value before normalization in the end of this function
		delete minNode; //new hn above is deleted here to save memo [info: delete minNode; automatically calls minNode->~FibHeapNode(), i.e. its destructor]

		//minD is an ascending value; so, when i reach far enough verts, i can safely skip the rest since all verts visited so far have finalized their shortest paths,
		//which means any vertex that i visit after this pnt will be have dist >= minD to sourceVert
		//if (minD > far)  //far enough vert reached
		//	break;

		//relax each edge, including shortcuts, incident to minDidx
		for (int ve = 0; ve < (int) verts[minDidx]->edgeList.size(); ve++)
		{
			int e = verts[minDidx]->edgeList[ve];
			int va = (edges[e]->v1i == minDidx ? edges[e]->v2i : edges[e]->v1i);
			if (verts[minDidx]->dSpanned + edges[e]->length < verts[va]->dSpanned)
			{
				verts[va]->dSpanned = verts[minDidx]->dSpanned + edges[e]->length; //relaxation
				h->DecreaseKey(verts[va]->heapNode, verts[va]->dSpanned); //worst O(lgV), amortized O(1)
			}
		}
		for (int ve = 0; ve < (int) verts[minDidx]->sedgeList.size(); ve++)
		{
			int e = verts[minDidx]->sedgeList[ve];
			int va = (sedges[e]->v1i == minDidx ? sedges[e]->v2i : sedges[e]->v1i);
			if (verts[minDidx]->dSpanned + sedges[e]->length < verts[va]->dSpanned)
			{
				verts[va]->dSpanned = verts[minDidx]->dSpanned + sedges[e]->length; //relaxation
				h->DecreaseKey(verts[va]->heapNode, verts[va]->dSpanned); //worst: O(lgV), amortized O(1)
			}
		}//*/
	} //end of while (true)
	/////// original Dijkstra to compute shortest paths from sourceVert to all others ///////

	delete h;

	if (maxGeoDist == -INF) //don't compute bad dSpanning[] w/ unhealthy min/maxGeoDist values in vain
		return;
	//normalized dSpanning[] from dSpanningUnn[]
	float denom = maxGeoDist - minGeoDist; //use the global min/maxGeoDist that spans the whole surface
	for (int i = 0; i < (int) verts.size(); i++)
		verts[sourceVert]->dSpanning[ i ] = (verts[sourceVert]->dSpanningUnn[i] - minGeoDist) / denom; //don't forget to change minGeoDist to minGeoDistG if you change denom	
}

inline float angleBetween(float* v1, float* v2)
{
	//returns smallest radian angle b/w normalized v1 & v2 via v1 DOT v2 = v1Len * v2Len * cosTheta (len = 1 for normalized, hence acos(dot) is the theta)

	float cosTheta = 0; //v1.dot(v2);
	for (int i = 0; i < 3; i++)
		cosTheta += v1[i]*v2[i];
	if (cosTheta < -1.0f)	cosTheta = -1.0f;
	if (cosTheta > 1.0f)	cosTheta = 1.0f;
	return acos(cosTheta); //* (180/PI); no radian to float to save time (besides computeGaussianCurvature() expects radian output)
}

inline int factorial(int k)
{
	if (k > 10)
		return -1; //don't bother big numbers
	int res = 1;
	for (int i = k; i > 0; i--)
		res *= i;
	return res;
}
void Mesh::fillPermutationMtrxs(int N)
{
	//fills permutation mtrxs perms2/3/4/5/6/7, e.g. perms3[1] = {2,0,1} and it has 3! = 6 such rows/arrays of size 3

	//permutations to be read from file to perms*[][] mtrxs
	FILE* fPtr;
	char fName[150];
	//2! permutations in perms2[][]
	int r = 2, factR = factorial(r);
	sprintf(fName, "input\\permutations\\0-start\\%d!.txt", r);
	if (! (fPtr = fopen(fName, "r")))
	{
		cout << fName << " not found\n";
		exit(0);
	}
	perms2 = new int*[factR];
	for (int i = 0; i < factR; i++) //r! reorderings/permutations
	{
		perms2[i] = new int[r];
		for (int j = 0; j < r; j++) //each of size r
			fscanf(fPtr, "%d ", &perms2[i][j]);
	}
	fclose(fPtr);

	//3! permutations in perms3[][]
	r = 3;
	factR = factorial(r);
	sprintf(fName, "input\\permutations\\0-start\\%d!.txt", r);
	if (! (fPtr = fopen(fName, "r")))
	{
		cout << fName << " not found\n";
		exit(0);
	}
	perms3 = new int*[factR];
	for (int i = 0; i < factR; i++) //r! reorderings/permutations
	{
		perms3[i] = new int[r];
		for (int j = 0; j < r; j++) //each of size r
			fscanf(fPtr, "%d ", &perms3[i][j]);
	}
	fclose(fPtr);

	//4! permutations in perms4[][]
	r = 4;
	factR = factorial(r);
	sprintf(fName, "input\\permutations\\0-start\\%d!.txt", r);
	if (! (fPtr = fopen(fName, "r")))
	{
		cout << fName << " not found\n";
		exit(0);
	}
	perms4 = new int*[factR];
	for (int i = 0; i < factR; i++) //r! reorderings/permutations
	{
		perms4[i] = new int[r];
		for (int j = 0; j < r; j++) //each of size r
			fscanf(fPtr, "%d ", &perms4[i][j]);
	}
	fclose(fPtr);

	//5! permutations in perms5[][]
	r = 5;
	factR = factorial(r);
	sprintf(fName, "input\\permutations\\0-start\\%d!.txt", r);
	if (! (fPtr = fopen(fName, "r")))
	{
		cout << fName << " not found\n";
		exit(0);
	}
	perms5 = new int*[factR];
	for (int i = 0; i < factR; i++) //r! reorderings/permutations
	{
		perms5[i] = new int[r];
		for (int j = 0; j < r; j++) //each of size r
			fscanf(fPtr, "%d ", &perms5[i][j]);
	}
	fclose(fPtr);

	//6! permutations in perms6[][]
	r = 6;
	factR = factorial(r);
	sprintf(fName, "input\\permutations\\0-start\\%d!.txt", r);
	if (! (fPtr = fopen(fName, "r")))
	{
		cout << fName << " not found\n";
		exit(0);
	}
	perms6 = new int*[factR];
	for (int i = 0; i < factR; i++) //r! reorderings/permutations
	{
		perms6[i] = new int[r];
		for (int j = 0; j < r; j++) //each of size r
			fscanf(fPtr, "%d ", &perms6[i][j]);
	}
	fclose(fPtr);

	//7! permutations in perms7[][]
	r = 7;
	factR = factorial(r);
	sprintf(fName, "input\\permutations\\0-start\\%d!.txt", r);
	if (! (fPtr = fopen(fName, "r")))
	{
		cout << fName << " not found\n";
		exit(0);
	}
	perms7 = new int*[factR];
	for (int i = 0; i < factR; i++) //r! reorderings/permutations
	{
		perms7[i] = new int[r];
		for (int j = 0; j < r; j++) //each of size r
			fscanf(fPtr, "%d ", &perms7[i][j]);
	}
	fclose(fPtr);

	if (N == 8) //don't allocate 8!=40320 bytes in vain for memo efficiency
	{
		//8! permutations in perms8[][]
		r = 8;
		factR = factorial(r);
		sprintf(fName, "input\\permutations\\0-start\\%d!.txt", r);
		if (! (fPtr = fopen(fName, "r")))
		{
			cout << fName << " not found\n";
			exit(0);
		}
		perms8 = new int*[factR];
		for (int i = 0; i < factR; i++) //r! reorderings/permutations
		{
			perms8[i] = new int[r];
			for (int j = 0; j < r; j++) //each of size r
				fscanf(fPtr, "%d ", &perms8[i][j]);
		}
		fclose(fPtr);
	}

	if (N == 9) //don't allocate 9!=362880 bytes in vain for memo efficiency
	{
		//9! permutations in perms9[][]
		r = 9;
		factR = factorial(r);
		sprintf(fName, "input\\permutations\\0-start\\%d!.txt", r);
		if (! (fPtr = fopen(fName, "r")))
		{
			cout << fName << " not found\n";
			exit(0);
		}
		perms9 = new int*[factR];
		for (int i = 0; i < factR; i++) //r! reorderings/permutations
		{
			perms9[i] = new int[r];
			for (int j = 0; j < r; j++) //each of size r
				fscanf(fPtr, "%d ", &perms9[i][j]);
		}
		fclose(fPtr);
	}
}
inline float Mesh::getCorrespListDistortion(Mesh* mesh2, vector< int > tmpSafeCorrespTry)
{
	//returns D_iso = distortion, i.e. subtraction of normalized geodesics, of the correspondence lying in the tmpSafeCorrespTry list using tmpSafeCorrespTry itself as the traversalList

	float isometryCostTry = 0.0f;
	int n = (int) tmpSafeCorrespTry.size() / 2;
	for (int s = 0; s < n; s++)
	{
		int bvs = tmpSafeCorrespTry[2*s], bvb = tmpSafeCorrespTry[2*s + 1]; //tmpSafeCorrespTry[even idxs] ~ this and [odd] ~ mesh2
		float vertexBasedIso = 0.0f;
		for (int ss = 0; ss < n; ss++) //traversal
		{
			int s1 = tmpSafeCorrespTry[2*ss], s2 = tmpSafeCorrespTry[2*ss + 1];
			if (s1 == bvs) //when bvs-bvb is a bad match, s1-s2 will wrongly award it by += 0 effect when s1=bvs (and hence s2=bvb); to prevent this, skip s1=bvs case; if bvs-bvb is a good match then this += 0 has no harm
				continue;

			//if bvs vs. bvb is a good correspondence, then geoDist(bvs, s1) = geoDist(bvb, s2) for all s1 & s2 in tmpSafeCorrespTry
			float t1s1 = verts[bvs]->dSpanning[s1], t2s2 = mesh2->verts[bvb]->dSpanning[s2]; //getCorrespListDistortion() called over this mesh; so, verts[bvs]->dSpanning is indeed this.verts[bvs]->dSpanning
			vertexBasedIso += fabs(t1s1 - t2s2); //normalized difference here is small if tbv1 vs. tbv2 is a good correspondence
		}
		if (n != 1) //no div by zero
			vertexBasedIso /= (n-1); //in computation of vertexBasedIso of bvs-bvb, 1 match, bvs-bvb itself, is skipped during traversal, and therefore n-1 += performed
		else
			vertexBasedIso = 0.0f; //due to symmetry problem, all costs will be the same and hence best selected arbitrarily
		isometryCostTry += vertexBasedIso;
		verts[bvs]->diso = vertexBasedIso; //individual distortion of this vert for interpCoarseToDense mode
	} //end of s
	//overall isometry cost of this tmpSafeCorrespTry combination
	isometryCostTry /= n;

	return isometryCostTry;
}
OutputPkt Mesh::matchWith(Mesh* mesh2, int N, int topK)
{
	//all possible 1-to-1 maps b/w extremes of this mesh and mesh2
	
	cout << "\nAll pairwise maps b/w mesh" << frameID << " & mesh" << mesh2->frameID << "......\n";

	///////////// permutations ///////////////////////
	//no combinations below 'cos i use all N of N extremes for sure
	fillPermutationMtrxs(N);
	//initially unosrted maps w/ their distortions; sorted distortions and corresponding maps'll be overwritten on these
	vector< vector< int > > candidateMapsUnsorted;
	vector< float > allDistortions;
	if ((int) extremes.size() != N || (int) mesh2->extremes.size() != N)
	{
		cout << N << " vs. " << ((int) extremes.size() != N ? (int) extremes.size() : (int) mesh2->extremes.size()) << " incompatible extreme sizes\n";
		exit(0);
	}

	int factR = factorial(N), twoR = 2*N;
	float tooBigDistortion = 0.14f;//, avgIsoCost = 0.0f; //don't bother sorting very bad mappings w/ tooBigDistortions at all; avg of all isometric distortion values computed in loop below gives an idea for tooBigDistortion; for scape0-scape1 pair avg = 0.265298
	vector< int > fillTheRest;
	for (int r = 0; r < factR; r++) //r! reordering of mesh2 extremes to be tested w/ the fixed mesh1 (this) quadruple
	{
		//refill tmpSafeCorrespTry according to the current reordering (one of k! permutations) of mesh2 extremes
		vector< int > tmpSafeCorrespTry;
		int idx1 = 0, idx2 = 0;
		for (int y = 0; y < twoR; y++)
			if (y % 2 == 0)
				tmpSafeCorrespTry.push_back(extremes[idx1++]); //even idxs from stable mesh1 side
			else //odd idxs for permutations from mesh2 side
			{
				if (N == 4)
					tmpSafeCorrespTry.push_back(mesh2->extremes[ perms4[r][idx2++] ]); 
				else if (N == 5)
					tmpSafeCorrespTry.push_back(mesh2->extremes[ perms5[r][idx2++] ]);
				else if (N == 6)
					tmpSafeCorrespTry.push_back(mesh2->extremes[ perms6[r][idx2++] ]);
				else if (N == 7)
					tmpSafeCorrespTry.push_back(mesh2->extremes[ perms7[r][idx2++] ]);
				else if (N == 8)
					tmpSafeCorrespTry.push_back(mesh2->extremes[ perms8[r][idx2++] ]);
				else //if (N == 9)
					tmpSafeCorrespTry.push_back(mesh2->extremes[ perms9[r][idx2++] ]);//*/
			}

		//compute isometry cost of this N x N mapping and store it along w/ the map itself, if it's not a very bad map
		float isometryCostTry = getCorrespListDistortion(mesh2, tmpSafeCorrespTry);
		if (isometryCostTry < tooBigDistortion) //very big distortions very unlikely to be ground-truth maps; so, don't even consider them in sorted list of candidate maps to be tracked
		{ 
			candidateMapsUnsorted.push_back(tmpSafeCorrespTry); //a candidate map for selection during shortest paths on graph
			allDistortions.push_back(isometryCostTry); //distortion of this candidate map
		}
		else //in the very rare case where allDistortions.size does not reach topK due to tooBigDistortion test above, i fill the remaining slots w/ this dummy unused map to prevent crashes in the loop on topK at the bottom
			fillTheRest = tmpSafeCorrespTry;
		//avgIsoCost += isometryCostTry;
	} //end of r
	//cout << "avg of all isometric distortions produced: " << avgIsoCost/factR << endl;
	while ((int) allDistortions.size() < topK) //very rare
	{
		candidateMapsUnsorted.push_back(fillTheRest); //this dummy map won't be selected by Graph.solve() due to its high distortion and i put it here just to guarantee that there's topK verts in each vertical layer of my graph
		allDistortions.push_back(0.99f); //dummy distortion
		cout << "dummy map to preserve the graph structure; no harm at all\n";
	}
	///////////// permutations ends ///////////////////////
	
	//sort allDistortions and the corresponding candidate maps; if there are k distortions (passing tooBigDistortion test), then there are k candidateMaps each of size N; in the end, both allDistortions & candidateMapsUnsorted will be sorted
	for (int i = 0; i < (int) allDistortions.size(); i++)
	{
		//get the min from the right side of [i] 'cos left of [i] already sorted and fixed
		float minDisto = INF;
		int minDistoIdx = -1;
		for (int j = i; j < (int) allDistortions.size(); j++)
			if (allDistortions[j] < minDisto)
			{
				minDisto = allDistortions[j];
				minDistoIdx = j;
			}
		//replace [i] & [minDistoIdx]
		float tmp = allDistortions[i];
		allDistortions[i] = minDisto;
		allDistortions[minDistoIdx] = tmp;
		//replace corresponding maps
		vector< int > tmp2 = candidateMapsUnsorted[i]; //remember i'th map
		candidateMapsUnsorted[i] = candidateMapsUnsorted[minDistoIdx]; //update i'th map to be this next min-distortion map
		candidateMapsUnsorted[minDistoIdx] = tmp2; //remembered i'th map goes to the slot of the next min-distortion used just-above
	}

	//put the first topK entry of sorted candidateMapsUnsorted along w/ theire distortions into result struct to be returned
	OutputPkt result;
	for (int i = 0; i < topK; i++)
	{
		result.candidateMapsSorted.push_back(candidateMapsUnsorted[i]);
		result.distortionsSorted.push_back(allDistortions[i]);
	}
	cout << "Done!\n\n";
	return result;
}

void Mesh::computeExtremes(int N)
{
	//sampling of N extreme samples; start w/ the most extreme (max mu) and follow by N-1 farthest-point extremes

	//this mesh already knows its extremes, related geodesics, .. due to processing in the previous iter, e.g., in the current iter of mesh4-mesh5, mesh4 already knows stuff due to previous mesh3-mesh4 processing
	if (! extremes.empty())
		return;

	cout << "\nExtremity sampling......\n";
	//////////////// ~40 bases for mu computation for mostExtreme decision ////////////////
	//some # of baseVerts will be created to find average geodesic distance, mu, values that gives us the most extreme sample later
	float wholeSurfaceArea = 0.0f, r;
	for (int t = 0; t < (int) tris.size(); t++)
		wholeSurfaceArea += tris[t]->area;

	//~40 evenly-spaced baseVerts for more accurate mu values
	r = 0.158f * sqrt(wholeSurfaceArea / PI) * 1.8f; //A = 40 pi r^2 ==> r = 0.1 sqrt(A/pi) gives 40 on a disk/patch; *~2 to make 40 on a closed surface (as 2 same patches connected to make the closed surface)
	//bases that are at least r apart from each other
	computeBaseVerts(r);

	minGeoDist = INF;
	maxGeoDist = -INF;
//	cout << "original Dijkstra from each baseVert being computed....\n";
	for (int i = 0; i < (int) baseVerts.size(); i++)
		dijkstraShortestPaths(baseVerts[i]); //dSpanning[] bad now 'cos min/maxGeoDist not ready yet (will be ready based on extremes[])
//	cout << "done!\n\n";

	//compute mu for extremities using healthy dSpanningUnn[], not dSpanning[]
	//for (int i = 0; i < (int) baseVerts.size(); i++) verts[ baseVerts[i] ]->patchArea = computePatchArea(baseVerts[i], r); //to be used in mu computation below
	//use baseVerts to determine the v.Mu values via v.Mu = SUM_i [g(v, bi) * bi.area] for each vert v
//	cout << "Mu (sum of geodesics from v to all verts) being computed....\n";
	int mostExtreme = -1;
	float maxMu = -INF;
	for (int v = 0; v < (int) verts.size(); v++)
	{
		for (int i = 0; i < (int) baseVerts.size(); i++)
			verts[v]->Mu += verts[ baseVerts[i] ]->dSpanningUnn[v];// * verts[ baseVerts[i] ]->patchArea; //patchArea weight
		//not normalizing v.Mu w/ (int) baseVerts.size() to save time 'cos Mu values across meshes not to be compared; i'll use these values just to compute the most extreme on this shape for which normalization is redundant
		if (verts[v]->Mu > maxMu)
		{
			maxMu = verts[v]->Mu;
			mostExtreme = v;
		}
	}
//	cout << "done!\n\nextremes being created..\n";
	//////////////// ~40 bases for mu computation for mostExtreme decision ends ////////////////

	//////////////// fps of N extremes starting from the mostExtreme ////////////////
	//fixed number (N) of fps (farthest point sampling) samples goes to extremes[] first of which is the mostExtreme
	for (int v = 0; v < (int) verts.size(); v++)
		verts[v]->processed = false; //initialize this to prevent duplications in extremes; true means added to extremes
	extremes.push_back(mostExtreme);
	dijkstraShortestPaths(mostExtreme); //dSpanningUnn[] to be used below is now ready, dSpanning is bad due to unhealthy min/maxGeoDist values
	verts[mostExtreme]->processed = true; //to prevent duplications in extremes; true means added to extremes
	//fill the remaining N-1 items via fps for uniform distribution that lets nice coverage on models
	//complexity of FPS below is O(NV + NVlogV), where V is for farthest point search, and VlogV is for dijkstra in the end that fills dSpanningUnn[] to be used in the next iter; one can reduce the first V to logV by efficient searching via kd-tree or sth, hence O(NlogV + NVlogV)
	while ((int) extremes.size() < N)
	{
		//in fps, the next sample candidate x_i \in all vertices finds its closest existing sample x_clo and remembers d_i = d(x_i, x_clo); if d_i is the max among all other candidates closest remembered dists, then x_i is selected as the next sample;
		//here my candidates are restricted to nothing, i.e. i use all verts as candidates; so, find the closest existing extreme e2 \in extremes to s \in verts and add s to extremes as nextSample if that closest dist is very big, i.e. the biggest among all closest dists for different s's
		float maxDist = -INF;
		int nextSample = -1;				
		for (int s = 0; s < (int) verts.size(); s++) //fps selects from amongst all vertices
		{
			int e2 = -1;
			if (verts[s]->processed)
				continue; //no chance to add s to iBaseVerts 'cos it's already added
			float minDist = verts[ extremes[0] ]->dSpanningUnn[s]; //=INF is tricky 'cos objects may be scaled to those high coordinates			
			for (int e = 0; e < (int) extremes.size(); e++) //find the closest existing base e \in extremes to s \in verts
				if (verts[ extremes[e] ]->dSpanningUnn[s] <= minDist)
				{
					minDist = verts[ extremes[e] ]->dSpanningUnn[s];
					e2 = extremes[e];
				}
			if (verts[e2]->dSpanningUnn[s] > maxDist)
			{
				maxDist = verts[e2]->dSpanningUnn[s];
				nextSample = s;
			}
		} //end of for s
if (nextSample == -1) {cout << "can't happen in fps\n"; exit(0);}
		extremes.push_back(nextSample);
		verts[nextSample]->processed = true; //added, so not add again later
		dijkstraShortestPaths(nextSample); //upcoming iteration needs extremes.dSpanningUnn which is not currently ready when selecting from the pool of all verts
	} //end of while N

	//dSpanning[] values of baseVerts and extremes are bad now due to undecided min/maxGeoDist above; make them null to let upcoming dijkstraShortestPaths() find healthy dSpanning[] for them, for extremes actually 'cos baseVerts won't be used anymore, thanks to the healthy mi/maxGeoDist values set here
	for (int i = 0; i < (int) extremes.size(); i++) //use extremes not baseVerts 'cos it is expected to represent max geo better
		for (int v = 0; v < (int) verts.size(); v++)
		{
			if (verts[ extremes[i] ]->dSpanningUnn[ v ] < minGeoDist)
				minGeoDist = verts[ extremes[i] ]->dSpanningUnn[ v ];
			if (verts[ extremes[i] ]->dSpanningUnn[ v ] > maxGeoDist)
				maxGeoDist = verts[ extremes[i] ]->dSpanningUnn[ v ];
if (verts[ extremes[i] ]->dSpanningUnn[ v ] == INF){cout<<extremes[i]<< " " << v << " not connected by a path\n"; exit(0);}
		}
//	cout << "min/max geodesic on the surface (for normalization): " << minGeoDist << " / " << maxGeoDist << "\nextremes to be matched: ";
	cout << frameID << "'s extremes to be matched: ";
if (maxGeoDist == INF) {cout << "disconnected vert(s) should be canceled\n";exit(0);} //guarantees that mesh is not disconnected (especially manually cut meshes can be bad); i may not crush if bases i'll use not hit these disconnected ones but be safe and stop here
	//clean old bad baseVerts and extremes
	for (int i = 0; i < (int) baseVerts.size(); i++)
	{
		delete verts[ baseVerts[i] ]->dSpanning;
		delete verts[ baseVerts[i] ]->dSpanningUnn; //arrays allocated with new[] must be deallocated with delete[] so i should've added [] here but not a big deal
		verts[ baseVerts[i] ]->dSpanning = verts[ baseVerts[i] ]->dSpanningUnn = NULL;
	}
	for (int i = 0; i < (int) extremes.size(); i++)
	{
		if (verts[ extremes[i] ]->dSpanning) //may already be deleted in for-loop above
			delete verts[ extremes[i] ]->dSpanning;
		if (verts[ extremes[i] ]->dSpanningUnn)
			delete verts[ extremes[i] ]->dSpanningUnn; //arrays allocated with new[] must be deallocated with delete[] so i should've added [] here but not a big deal
		verts[ extremes[i] ]->dSpanning = verts[ extremes[i] ]->dSpanningUnn = NULL;

		//dSpanning[] good now 'cos min/maxGeoDist ready; dSpanning[]'ll be used during distortion computations of maps b/w shape extremes
		dijkstraShortestPaths(extremes[i]);

		cout << extremes[i] << " ";
//for (int v = 0; v < (int) verts.size(); v++) if (verts[ extremes[i] ]->dSpanning[v] < 0 || verts[ extremes[i] ]->dSpanning[v] > 1) {cout << verts[ extremes[i] ]->dSpanning[v] << " wwwwwrong"; exit(0); }//delLaterrrrrrrrrr
	}
	cout << "\n";
	//////////////// fps of N extremes starting from the mostExtreme ends ////////////////
}
void Mesh::computeBaseVerts(float r)
{
	//evenly-spaced baseVerts samples that are at least r apart from each other

//	cout << frameID << "'s base vertices for upcoming mu computation being prepared (r = " << r << ")....\n";

	//get the next available base vertex
	int baseVertIdx = -1;
	for (int v = 0; v < (int) verts.size(); v++)
		if (verts[v]->baseCandidate)
			baseVertIdx = v;

	vector< int >* vlist = new vector<int>; //first a baseVert enters this list, and then verts in the area/patch of that baseVert

	//output dependent: O(k * V), where k = # of baseVerts generated, and V is # of verts
	while (baseVertIdx != -1) //there is still a baseVert whose area to be expanded below
	{
		/////// shortcut Dijkstra w/ slow arrays to create 1 baseVert ///////
		//step 1: initialization
		for (int v = 0; v < (int) verts.size(); v++)
			verts[v]->d = INF;

		//step 2: select a base vert
		verts[baseVertIdx]->d = 0.0f; //so that this baseVert is initially selected from vlist in the loop below
		vlist->push_back(baseVertIdx);		//used below
		baseVerts.push_back(baseVertIdx);	//used later; whole aim is to fill this vector indeed
		verts[baseVertIdx]->baseCandidate = false;
		while (! vlist->empty())
		{
			//step 3: extract vert w/ min d value
			float minD = INF;
			int minDidx, v;
			for (int i = 0; i < (int) vlist->size(); i++)
				if (verts[ (*vlist)[i] ]->d < minD) //dereferencing
				{
					v = (*vlist)[i];
					minD = verts[v]->d;
					minDidx = i;
				}
			//erase the element w/ min d value
			vlist->erase(vlist->begin() + minDidx); //ve.erase(begin + 4) erases ve[4]

			//step 4: relaxation loop that updates d values and expands baseVert.area as long as d is small enough
			for (int ve = 0; ve < (int) verts[v]->edgeList.size(); ve++)
			{
				int e = verts[v]->edgeList[ve];
				int va = (edges[e]->v1i == v ? edges[e]->v2i : edges[e]->v1i);
				if (verts[v]->d + edges[e]->length < verts[va]->d)
				{
					verts[va]->d = verts[v]->d + edges[e]->length;
					if (verts[va]->d < r)
					{
						vlist->push_back(va);
						verts[va]->baseCandidate = false;
					}
				}
			}
			for (int ve = 0; ve < (int) verts[v]->sedgeList.size(); ve++)
			{
				int e = verts[v]->sedgeList[ve];
				int va = (sedges[e]->v1i == v ? sedges[e]->v2i : sedges[e]->v1i);
				if (verts[v]->d + sedges[e]->length < verts[va]->d)
				{
					verts[va]->d = verts[v]->d + sedges[e]->length;
					if (verts[va]->d < r)
					{
						vlist->push_back(va);
						verts[va]->baseCandidate = false;
					}
				}
			}
		} //end of vlist.empty()
		/////// shortcut Dijkstra w/ slow arrays to create 1 baseVert ends ///////
#ifdef DBG_DETAIL
//		cout << "base" << baseVerts.size()-1 << ": " << baseVertIdx << endl;
#endif
		//get the next available base vertex
		baseVertIdx = -1;
		for (int v = 0; v < (int) verts.size(); v++)
			if (verts[v]->baseCandidate)
				baseVertIdx = v;
	} //end of baseVertIdx

	delete vlist;

//	cout << "done! # of base vertices: " << (int) baseVerts.size() << ", for " << (int) verts.size() << " verts\n\n";
}

/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//////////////////////////// everything below is for dense matching of dense samples[] instead of coarse extremes[] /////////////////////////////////
void Mesh::computeEvenlySpacedSamples(int N, int D)
{
	//sampling of N evenly-spaced FPS samples that are not necessarily on extremes; since N is big here nonextreme is not an issue for our matching performance

	//this mesh already knows its samples, related geodesics, .. due to processing in the previous iter, e.g., in the current iter of mesh4-mesh5, mesh4 already knows stuff due to previous mesh3-mesh4 processing
	if (! samples.empty())
		return;

	cout << "\nEvenly-spaced sampling......\n";

	//fixed number (N) of fps (farthest point sampling) samples goes to samples[] first of which is an arbitrary vertex
	for (int v = 0; v < (int) verts.size(); v++)
		verts[v]->processed = false; //initialize this to prevent duplications in samples; true means added to samples
	int arbitVert = 0;
	samples.push_back(arbitVert);
	dijkstraShortestPaths(arbitVert); //dSpanningUnn[] to be used below is now ready, dSpanning is bad due to unhealthy min/maxGeoDist values
	verts[arbitVert]->processed = true; //to prevent duplications in samples; true means added to samples
	//fill the remaining N-1 items via fps for uniform distribution that lets nice coverage on models
	//complexity of FPS below is O(NV + NVlogV), where V is for farthest point search, and VlogV is for dijkstra in the end that fills dSpanningUnn[] to be used in the next iter; one can reduce the first V to logV by efficient searching via kd-tree or sth, hence O(NlogV + NVlogV)
	while ((int) samples.size() < N)
	{
		//in fps, the next sample candidate x_i \in all vertices finds its closest existing sample x_clo and remembers d_i = d(x_i, x_clo); if d_i is the max among all other candidates closest remembered dists, then x_i is selected as the next sample;
		//here my candidates are restricted to nothing, i.e. i use all verts as candidates; so, find the closest existing sample e2 \in samples to s \in verts and add s to samples as nextSample if that closest dist is very big, i.e. the biggest among all closest dists for different s's
		float maxDist = -INF;
		int nextSample = -1;				
		for (int s = 0; s < (int) verts.size(); s++) //fps selects from amongst all vertices
		{
			int e2 = -1;
			if (verts[s]->processed)
				continue; //no chance to add s to iBaseVerts 'cos it's already added
			float minDist = verts[ samples[0] ]->dSpanningUnn[s]; //=INF is tricky 'cos objects may be scaled to those high coordinates			
			for (int e = 0; e < (int) samples.size(); e++) //find the closest existing base e \in samples to s \in verts
				if (verts[ samples[e] ]->dSpanningUnn[s] <= minDist)
				{
					minDist = verts[ samples[e] ]->dSpanningUnn[s];
					e2 = samples[e];
				}
			if (verts[e2]->dSpanningUnn[s] > maxDist)
			{
				maxDist = verts[e2]->dSpanningUnn[s];
				nextSample = s;
			}
		} //end of for s
if (nextSample == -1) {cout << "can't happen in fps\n"; exit(0);}
		samples.push_back(nextSample);
		verts[nextSample]->processed = true; //added, so not add again later
		dijkstraShortestPaths(nextSample); //upcoming iteration needs samples.dSpanningUnn which is not currently ready when selecting from the pool of all verts
	} //end of while N

	//dSpanning[] values of samples are bad now due to undecided min/maxGeoDist above; make them null to let upcoming dijkstraShortestPaths() find healthy dSpanning[] for them
	for (int i = 0; i < (int) samples.size(); i++)
		for (int v = 0; v < (int) verts.size(); v++)
		{
			if (verts[ samples[i] ]->dSpanningUnn[ v ] < minGeoDist)
				minGeoDist = verts[ samples[i] ]->dSpanningUnn[ v ];
			if (verts[ samples[i] ]->dSpanningUnn[ v ] > maxGeoDist)
				maxGeoDist = verts[ samples[i] ]->dSpanningUnn[ v ];
if (verts[ samples[i] ]->dSpanningUnn[ v ] == INF){cout<<samples[i]<< " " << v << " not connected by a path\n"; exit(0);}
		}
//	cout << "min/max geodesic on the surface (for normalization): " << minGeoDist << " / " << maxGeoDist << "\nsamples to be matched: ";
	cout << frameID << "'s samples to be matched: ";
if (maxGeoDist == INF) {cout << "disconnected vert(s) should be canceled\n";exit(0);} //guarantees that mesh is not disconnected (especially manually cut meshes can be bad); i may not crush if bases i'll use not hit these disconnected ones but be safe and stop here
	//clean old samples.dSpanning values
	for (int i = 0; i < (int) samples.size(); i++)
	{
		if (verts[ samples[i] ]->dSpanning) //may already be deleted in for-loop above
			delete verts[ samples[i] ]->dSpanning;
		if (verts[ samples[i] ]->dSpanningUnn)
			delete verts[ samples[i] ]->dSpanningUnn; //arrays allocated with new[] must be deallocated with delete[] so i should've added [] here but not a big deal
		verts[ samples[i] ]->dSpanning = verts[ samples[i] ]->dSpanningUnn = NULL;

		//dSpanning[] good now 'cos min/maxGeoDist ready; dSpanning[]'ll be used during distortion computations of maps b/w shape samples
		dijkstraShortestPaths(samples[i]);

		if (i < 5 || i >= N-5) //print the first and last 5 sample ids
			cout << samples[i] << " ";
		else if (i == 5)
			cout << " .. ";
	}
	cout << "\n";

	//geodesic affinity matrix for spectral embedding
	if (D != -1) //call from interpCoarseToDense not need specrtal stuff
		computeAffinityMtrx(N, D);
}

void negateK(float* a, int D) { for (int i = 0; i < D; i++) a[i] = -a[i]; }
void mtrxMult(float** result, int m1Row, int m1Col, int m2Row, int m2Col, float** m1, float** m2)
{
	//result = m1 * m2
	
	//init result to 0
	for (int r = 0; r < m1Row; r++)
		for (int c = 0; c < m2Col; c++)
			result[r][c] = 0.0f;

	/*pseudocode:
	for each row i in matrix1 do:
            for each column j in matrix 2 do:
                for each column k in matrix1 do:
                    increase  matrix3[i,j] by matrix1[i,k]*matrix2[k,j] */
	for (int m1r = 0; m1r < m1Row; m1r++)
		for (int m2c = 0; m2c < m2Col; m2c++)
			for (int m1c = 0; m1c < m1Col; m1c++)
				result[m1r][m2c] += m1[m1r][m1c] * m2[m1c][m2c]; //call-by-ref to result
}
inline float distanceBetweenKD(float* v1, float* v2, int k)
{
	//L2 distance b/w two k-D pnts
	float ss = 0;
	for (int i = 0; i < k; i++)
		ss += ((v2[i]-v1[i])*(v2[i]-v1[i]));
	return (float) sqrt(ss);
}
void translateAndScaleK(float** &pnts, int N, int D)
{
	//translates D-D pnts so that center of mass is at the origin (0,0,0) and scale them so that avg dist to origin is 1

	float* m1 = new float[D], * m2 = new float[D]; //mass ~ mean names used interchangably
	float* sum = new float[D];
	for (int i = 0; i < D; i++)
		sum[i] = 0.0f;
	for (int i = 0; i < N; i++)
		for (int c = 0; c < D; c++)
			sum[c] += pnts[i][c];
	for (int c = 0; c < D; c++)
		m1[c] = sum[c] / (float) N;
	//initial translation (t) that moves shape centers to (0,0,0,0,0) when D = 5
	for (int i = 0; i < N; i++)
		for (int c = 0; c < D; c++)
			pnts[i][c] -= m1[c];
	//now 2 models are optimally translation-aligned
	for (int c = 0; c < D; c++)
		m1[c] = m2[c] = 0.0f; //m1 and m2 moved to center after translation above
	//initial scaling (s) that makes mean variance from origin 1, i.e. scale each model to an appropriate size via avgDistToOrigin
	float avgDistToOrigin = 0.0f;
	for (int i = 0; i < N; i++)
		//taking pow2 below makes sum of means from origin to be 1; w/o pow2 it makes sum of distances from origin 1, which i want
		avgDistToOrigin += distanceBetweenKD(pnts[i], m1, D); //m1 = (0,0,0,0) when D = 4
	for (int i = 0; i < N; i++)
		for (int c = 0; c < D; c++)
			pnts[i][c] /= avgDistToOrigin;
	//now 2 models are optimally scale-aligned by Horn; m1 & m2 cannot change after scaling above; so, no redundant update here
	delete [] m1; delete [] m2; delete [] sum;
}
OutputPkt Mesh::spectralMatchWith(Mesh* mesh2, int N, int D, int topK)
{
	//2^D possible 1-to-1 min-weight perfect matchings b/w samples of this mesh and mesh2 for 2^D configurations suggested by eigenvector sign flips to create D-dimensional spectral embedding
	
	cout << "\nSome pairwise maps b/w mesh" << frameID << " & mesh" << mesh2->frameID << "......\n";

	//initially unsorted maps w/ their distortions; sorted distortions and corresponding maps'll be overwritten on these
	vector< vector< int > > candidateMapsUnsorted;
	vector< float > allDistortions;
	if ((int) samples.size() != N || (int) mesh2->samples.size() != N)
	{
		cout << N << " vs. " << ((int) samples.size() != N ? (int) samples.size() : (int) mesh2->samples.size()) << " incompatible extreme sizes\n";
		exit(0);
	}

	//////////////////////////////// map creation ////////////////////////////////
	//v.spectralD and eigen_vectors/values[] made ready by a prev computeAffinityMtrx() calls
	//eigvecs of geodesic affinity matrices may suffer from reordering and sign flips; each such eigen config gives me 1 map to be inserted in candidateMapsUnsorted along w/ its distortion	
	int n1 = (int) samples.size(), n2 = (int) mesh2->samples.size();

	int* vIdx1 = new int[n1], * vIdx2 = new int[n2]; //if r,c = 3,7 is a match in costs[r][c], i.e. [3][7]=0, then samples[3].matchIdx = vIdx2[7], hence vIdx2[7] is an idx to verts, unlike c which is an idx to samples[]
	//spectral coords and hence initial pntSet(D)1/2 are already centered at origin and scaled so that avg dist to origin is 1
	float** pntSetK1 = new float*[n1]; //for mesh1's (this) spectralCoords (of samples)
	for (int i = 0; i < n1; i++)
	{
		pntSetK1[i] = new float[D];
		vIdx1[i] = samples[i];
		for (int c = 0; c < D; c++)
			pntSetK1[i][c] = verts[ samples[i] ]->spectralD[c]; //fixed; so fill it now and update never
	}
	float** pntSetK2 = new float*[n2];
	for (int i = 0; i < n2; i++)
	{
		pntSetK2[i] = new float[D]; //changes according to the current eigen config	
		vIdx2[i] = mesh2->samples[i];
	}

//eigen_values in descending order as we go from [0] to [N]; eigen_vectors.column1 = eigvector1 (corresponding to largest eigval), col2 = eigvec2, and so on..
//so, eigen_vectors[1..7][i] is 1st 7 components, e.g. x,y,z,.., of ith eigenvector
	float** eigvecs1 = new float*[D]; //D leading eigvecs (each of which has N components) of mesh1.affinityMtrx; each row is an eigvec
	for (int i = 1; i <= D; i++) //i=0 is 1st eigvec w/ huge eigval and therefore ignored
	{
		eigvecs1[i-1] = new float[n1]; //kxn1 mtrx
		for (int j = 0; j < n1; j++)
			eigvecs1[i-1][j] = eigen_vectors[j][i]; //eigen_vectors[1..n1][0] is ignored 'cos it is constant corresponding to the largest eigvalue
	}
	float* eigvals1 = new float[D]; //corresponding eigvalues
	//fixed eigvecs and eigvals for mesh1 filled once and never updated as opposed to mesh2 stuff below
	for (int i = 1; i <= D; i++)
		eigvals1[i-1] = eigen_values[i]; //eigen_values[0] will never be accessed (huge eigval)
	//eigvecs/vals for mesh2 to be updated inside the loop
	float** eigvecs2 = new float*[D]; //D leading eigvecs of mesh2.affinityMtrx (each row is an eigvec)
	for (int i = 1; i <= D; i++)
	{
		eigvecs2[i-1] = new float[n2]; //kxn2 mtrx
		for (int j = 0; j < n2; j++) //init here in case i skip reordering below
			eigvecs2[i-1][j] = mesh2->eigen_vectors[j][i];
	}
	float* eigvals2 = new float[D]; //corresponding eigvalues to be filled and refilled inside the loop below
	for (int i = 1; i <= D; i++) //init here in case i skip reordering below
		eigvals2[i-1] = mesh2->eigen_values[i];
	//auxiliary arrays to get the new mesh2.spectrals
	float** spectralCoords2 = new float*[D]; //desired spectrals according to the current eigen config
	for (int i = 0; i < D; i++) //alloc but not fill yet
		spectralCoords2[i] = new float[n2]; //kxn2 mtrx
	float** diagK2 = new float*[D]; //kxk
	for (int i = 1; i <= D; i++) //use 2nd, 3rd, 4th eigvecs as suggested by paper (for 3-D)
		diagK2[i-1] = new float[D];	

//	for (int r = 0; r < D!; r++) //D! reordering; skipping reordering stuff 'cos my dynamic programming graph proecssing working over the entire shape collection alleviates the slight problems due to this choice
	//arbitrary sign flip handling
	int maxJ = (int) pow(2.0f, D); //= 16 for D = 4, 32 for D = 5
	for (int j = 0; j < maxJ; j++) //try each orientation of mesh2.eigen_vectors
	{
		///////// +- eigvector handling ////////////
		//one of 2^D permutations
		if (j == 1 || j == 9 || j == 17 || j == 25 || j == 33 || j == 41 || j == 49 || j == 57) //do nothing for j = 0
			negateK(eigvecs2[2], n2); //negate effect, i.e. eigvecs2[2] = -eigvecs2[2]; ++-++- for j = 33, ++-+-- for j = 49, ++---- for j = 57
		else if (j == 2 || j == 10 || j == 18 || j == 26 || j == 34 || j == 42 || j == 50 || j == 58) //+-++++ for j = 2 || +-+-++ for j = 10 || +-++-+ for j = 18 || +-+--+ for j = 26 || +-+++- for j = 34 || +-+-+- for j = 42 || +-++-- for j = 50 || +-+--- for j = 58
		{
			negateK(eigvecs2[1], n2);
			negateK(eigvecs2[2], n2);
		}
		else if (j == 3 || j == 11 || j == 19 || j == 27 || j == 35 || j == 43 || j == 51 || j == 59) //+--+++ || +---++ || +--+-+ || +----+ || +--++- || +---+- || +--+-- || +-----
			negateK(eigvecs2[2], n2);
		else if (j == 4 || j == 12 || j == 20 || j == 28 || j == 36 || j == 44 || j == 52 || j == 60) //-+++++ || -++-++ || -+++-+ || -++--+ || -++++- || -++-+- || -+++-- || -++---
		{
			negateK(eigvecs2[0], n2);
			negateK(eigvecs2[1], n2);
			negateK(eigvecs2[2], n2);
		}
		else if (j == 5 || j == 13 || j == 21 || j == 29 || j == 37 || j == 45 || j == 53 || j == 61) //-+-+++ || -+--++ || -+-+-+ || -+---+ || -+-++- || -+--+- || -+-+-- || -+----
			negateK(eigvecs2[2], n2);
		else if (j == 6 || j == 14 || j == 22 || j == 30 || j == 38 || j == 46 || j == 54 || j == 62) //--++++ || --+-++ || --++-+ || --+--+ || --+++- || --+-+- || --++-- || --+---
		{
			negateK(eigvecs2[1], n2);
			negateK(eigvecs2[2], n2);
		}
		else if (j == 7 || j == 15 || j == 23 || j == 31 || j == 39 || j == 47 || j == 55 || j == 63) //---+++ || ----++ || ---+-+ || -----+ || ---++- || ----+- || ---+-- || ------
			negateK(eigvecs2[2], n2);
		else if (j == 8 || j == 16 || j == 24 || j == 32 || j == 40 || j == 48 || j == 56) //+++-++ || ++++-+ || +++--+ || +++++- || +++-+- || ++++-- || +++---
		{
			negateK(eigvecs2[0], n2);
			negateK(eigvecs2[1], n2);
			negateK(eigvecs2[2], n2); //reset to original +++ for first 3 axes
			negateK(eigvecs2[3], n2); //and flip 4th axes for the next 7 iters
			if (j == 16 || j == 32 || j == 48) //if we enter here, it means D = 5 (or 6), therefore we have 5 (or 6) principal axes
				negateK(eigvecs2[4], n2); //and flip the 5th axes for the next 16 iters
			if (j == 32) //if we enter here, it means D = 6, therefore we have 6 principal axes
				negateK(eigvecs2[5], n2); //and flip the 6th axes for the next 16 iters
		}
		///////// +- eigvector handling ends ////////////

		//obtain new mesh2.spectralCoords w/ these orders and signs of eigvectors
		//compute 2 matrices that need to be multiplied to reach spectralCoords
		for (int i = 0; i < D; i++) //use 2nd, 3rd, 4th eigvecs as suggested by paper (for 3-D)
			for (int j = 0; j < D; j++)
				diagK2[i][j] = (i == j ? sqrt(eigvals2[i]) : 0.0f); //use sqrt(eigval) as suggested by literature, i.e. for better eigval scaling

		//desired result (of size kxn2) filled by call-by-ref
		mtrxMult(spectralCoords2, D, D, D, n2, diagK2, eigvecs2); //eigvecs2 is the mtrx whose rows are eigvecs; so, same as the one referred as ekt in PCAonAffinity()

		for (int i = 0; i < n2; i++)
			for (int c = 0; c < D; c++)
				pntSetK2[i][c] = spectralCoords2[c][i];

		//since pntSet1 is translated to origin and scaled so that avg dist to origin is 1 during PCAonAffinity(), do same for pntSet2 for a good comparison
		translateAndScaleK(pntSetK2, n2, D);

		///////// blossom on current spectral coords ////////////
		//before performing O(N^2 lgN) blossom, i may skip some very bad configs by O(N^2) closest-pnt-differences-sum test (tooBigDistortion = computeMSEKD(pntSetK1, pntSetK2, n1, n2, D, false);) but not doing that
		//compute isometry cost of this upcoming N x N blossom mapping and store it along w/ the map itself
		int nRows = (n1 > n2 ? n1 : n2), nCols = nRows; //square mtrx
		float discourage = INF, ** costs = new float*[nRows]; //[3][8]: costs of assigning samples[3] to mesh2.samples[8]
		for (int r = 0; r < nRows; r++)
		{
			costs[r] = new float[nCols];
			for (int c = 0; c < nCols; c++)
				costs[r][c] = ((r < n1 && c < n2) ? distanceBetweenKD(pntSetK1[r], pntSetK2[c], D) : discourage); //since n1=n2 for this project discourage will never be used but anyway
				//remember this entry in matchIdxs to access it in isometry cost computation of the resulting blossom match //no not needed anymore: verts[ samples[r] ]->matchIdxDists.push_back((float) mesh2->samples[c]); verts[ samples[r] ]->matchIdxDists.push_back(costs[r][c]);
		}
		//////////// blossom input ///////////////
		int nVertNum = nRows*2, nEdgeNum = nRows*nRows, e = 0, * edges = new int[2*nEdgeNum];
		double* weightsD = new double[nEdgeNum];
		for (int i = 0; i < nRows; i++)
			for (int j = 0; j < nCols; j++)
			{
				edges[2*e] = i;
				edges[2*e+1] = nRows+j;
				weightsD[e++] = (double) costs[i][j];
			}
		//////////// blossom input ends ///////////////
		Blossom* blossom = new Blossom();
		int* matches = new int[nVertNum]; //see this as lines of matches where 2 consecutive items make a matching line, e.g. 369 694
		blossom->match(nVertNum, nEdgeNum, edges, weightsD, matches); //matches update by call-by-ref
		delete blossom;
		delete [] edges;
		delete [] weightsD;
		vector< int > tmpSafeCorrespTry;		
		for (int m = 0; m < nRows; m++)
		{
			int m1 = matches[2*m]; //left side of the (m+1)'th matching line
			if (m1 < n1)
			{
				int m2 = matches[2*m + 1] - nRows; //mathches[m1 + 1], right side of the line, is in [nRows..2*nRows-1], hence nRows entry at mesh2
				int mi = (m2 < n2 ? vIdx2[m2] : -1); //if virtual column(s) exists in initial costs, some rows will not have valid matches, hence -1
				if (mi != -1)
				{
					//match
					tmpSafeCorrespTry.push_back(vIdx1[m1]); //base1 from this mesh1; be sure even idxs ~ mesh1, odds ~ mesh2
					tmpSafeCorrespTry.push_back(mi); //base2 from mesh2
					/*//and its cost by learning the individual distortion used in blossom mtrx
					float vertexBasedIso = -1.0f;
					for (int j = 0; j < (int) verts[ vIdx1[m1] ]->matchIdxDists.size(); j += 2)
						if (verts[ vIdx1[m1] ]->matchIdxDists[j] == mi)
						{
							vertexBasedIso = verts[ vIdx1[m1] ]->matchIdxDists[j+1];
							break;
						}
if (vertexBasedIso == -1.0f){cout<<"\n\n\n\n\nsth wrong\n";exit(0);}
					//spectral coords are scaled so that avg dist to origin is 1; due to this scaling vertexBasedIso is expected to be less than 1 but nevertheless += path->overallDistortion in Graph.cpp will not be that compatible w/ this L2-based initial distortion values; so, instead of using this costs[][]-based distortion here, feed tmpSafeCorrespTry to getCorrespListDistortion() which makes initial distortions (b/w meshes[0] & [1]) compatible w/ Graph.cpp
					isometryCostTry += vertexBasedIso;*/
				}
			}
		} //end of m
		delete [] matches;
		float isometryCostTry = getCorrespListDistortion(mesh2, tmpSafeCorrespTry); //as explained above, += vertexBasedIso tactic is wrong 'cos that normalized L2-based dist is not compatible w/ normalized geodesic-based distortion; so do this way
		//store isometry cost of this N x N mapping along w/ the map itself
		candidateMapsUnsorted.push_back(tmpSafeCorrespTry); //a candidate map for selection during shortest paths on graph
		allDistortions.push_back(isometryCostTry); //distortion of this candidate map
		///////// blossom on current spectral coords ////////////
	} //end of j

	//re-capture memo
	//delete [] eigen_values;
	//delete [] mesh2->eigen_values;
	//for (int i = 0; i < n1; i++) delete [] eigen_vectors[i]; delete [] eigen_vectors;
	//for (int i = 0; i < n2; i++) delete [] mesh2->eigen_vectors[i]; delete [] mesh2->eigen_vectors; don't clear these 'cos consecutive mesh in shape collection will access these
	for (int i = 0; i < D; i++) delete [] eigvecs1[i]; delete [] eigvecs1;
	for (int i = 0; i < D; i++) delete [] eigvecs2[i]; delete [] eigvecs2;
	delete [] eigvals1;	delete [] eigvals2;
	for (int i = 0; i < D; i++) delete [] spectralCoords2[i]; delete [] spectralCoords2;
	for (int i = 0; i < D; i++) delete [] diagK2[i]; delete [] diagK2;
	for (int i = 0; i < n1; i++) delete [] pntSetK1[i];	delete [] pntSetK1;
	for (int i = 0; i < n2; i++) delete [] pntSetK2[i];	delete [] pntSetK2;
	//////////////////////////////// map creation ends ////////////////////////////////
	
	//sort allDistortions and the corresponding candidate maps; if there are k distortions (passing tooBigDistortion test), then there are k candidateMaps each of size N; in the end, both allDistortions & candidateMapsUnsorted will be sorted
	for (int i = 0; i < (int) allDistortions.size(); i++)
	{
		//get the min from the right side of [i] 'cos left of [i] already sorted and fixed
		float minDisto = INF;
		int minDistoIdx = -1;
		for (int j = i; j < (int) allDistortions.size(); j++)
			if (allDistortions[j] < minDisto)
			{
				minDisto = allDistortions[j];
				minDistoIdx = j;
			}
		//replace [i] & [minDistoIdx]
		float tmp = allDistortions[i];
		allDistortions[i] = minDisto;
		allDistortions[minDistoIdx] = tmp;
		//replace corresponding maps
		vector< int > tmp2 = candidateMapsUnsorted[i]; //remember i'th map
		candidateMapsUnsorted[i] = candidateMapsUnsorted[minDistoIdx]; //update i'th map to be this next min-distortion map
		candidateMapsUnsorted[minDistoIdx] = tmp2; //remembered i'th map goes to the slot of the next min-distortion used just-above
	}

	//put the first topK entry of sorted candidateMapsUnsorted along w/ theire distortions into result struct to be returned
	OutputPkt result;
	for (int i = 0; i < topK; i++)
	{
		result.candidateMapsSorted.push_back(candidateMapsUnsorted[i]);
		result.distortionsSorted.push_back(allDistortions[i]);
	}
	cout << "Done!\n\n";
	return result;
}

////////////// eigenvector/value finder ////////////////
//adapted from Tina4 (tina-vision.net/tina4/doxygen/html/mat__eigen_8c-source.html)
#define SIGN(a,b) ((b)<0 ? -fabs(a) : fabs(a))
static void tred2(int n, float **a, float *d, float *e)
{
	int l, k, j, i;
	float scale, hh, h, g, f;

	//start arrays at 0 not 1
	for (i = 0; i < n; ++i)
		(a[i])--;
	a--;
	d--;
	e--;
	
	for (i = n; i >= 2; i--)
	{
		l = i - 1;
		h = scale = 0.0;
		if (l > 1)
		{
			for (k = 1; k <= l; k++)
				scale += fabs(a[i][k]);
			if (scale == 0.0)
				e[i] = a[i][l];
			else
			{
				for (k = 1; k <= l; k++)
				{
					a[i][k] /= scale;
					h += a[i][k] * a[i][k];
				}
				f = a[i][l];
				g = f > 0 ? -sqrt(h) : sqrt(h);
				e[i] = scale * g;
				h -= f * g;
				a[i][l] = f - g;
				f = 0.0;
				for (j = 1; j <= l; j++)
				{
					a[j][i] = a[i][j] / h;
					g = 0.0;
					for (k = 1; k <= j; k++)
						g += a[j][k] * a[i][k];
					for (k = j + 1; k <= l; k++)
						g += a[k][j] * a[i][k];
					e[j] = g / h;
					f += e[j] * a[i][j];
				}
				hh = f / (h + h);
				for (j = 1; j <= l; j++)
				{
					f = a[i][j];
					e[j] = g = e[j] - hh * f;
					for (k = 1; k <= j; k++)
						a[j][k] -= (f * e[k] + g * a[i][k]);
				}
			}
		}
		else
			e[i] = a[i][l];
		d[i] = h;
	}
	d[1] = 0.0;
	e[1] = 0.0;
	for (i = 1; i <= n; i++)
	{
		l = i - 1;
		if (d[i])
		{
			for (j = 1; j <= l; j++)
			{
					g = 0.0;
				for (k = 1; k <= l; k++)
					g += a[i][k] * a[k][j];
				for (k = 1; k <= l; k++)
					a[k][j] -= g * a[k][i];
			}
		}
		d[i] = a[i][i];
		a[i][i] = 1.0;
		for (j = 1; j <= l; j++)
			a[j][i] = a[i][j] = 0.0;
	}

	a++;
	d++;
	e++;
	for (i = 0; i < n; ++i)
		(a[i])++;
}
static int tqli(int n, float **z, float *d, float *e)
{
	int m, l, iter, i, k;
	float s, r, p, g, f, dd, c, b;

	//start arrays at 0 not 1
	for (i = 0; i < n; ++i)
		(z[i])--;
	z--;
	d--;
	e--;

	for (i = 2; i <= n; i++)
		e[i - 1] = e[i];
	e[n] = 0.0;
	for (l = 1; l <= n; l++)
	{
		iter = 0;
		do
		{
			for (m = l; m <= n - 1; m++)
			{
				dd = fabs(d[m]) + fabs(d[m + 1]);
				if (fabs(e[m]) + dd == dd)
					break;
			}
			if (m != l)
			{
				if (iter++ == 30)
					return (0);
				g = (d[l + 1] - d[l]) / (2.0f * e[l]);
				r = sqrt((g * g) + 1.0f);
				g = d[m] - d[l] + e[l] / (g + SIGN(r, g));
				s = c = 1.0;
				p = 0.0;
				for (i = m - 1; i >= l; i--)
				{
					f = s * e[i];
					b = c * e[i];
					if (fabs(f) >= fabs(g))
					{
						c = g / f;
						r = sqrt((c * c) + 1.0f);
						e[i + 1] = f * r;
						c *= (s = 1.0f / r);
					}
					else
					{
						s = f / g;
						r = sqrt((s * s) + 1.0f);
						e[i + 1] = g * r;
						s *= (c = 1.0f / r);
					}
					g = d[i + 1] - p;
					r = (d[i] - g) * s + 2.0f * c * b;
					p = s * r;
					d[i + 1] = g + p;
					g = c * r - b;
					for (k = 1; k <= n; k++)
					{
						f = z[k][i + 1];
						z[k][i + 1] = s * z[k][i] + c * f;
						z[k][i] = c * z[k][i] - s * f;
					}
				}
				d[l] = d[l] - p;
				e[l] = g;
				e[m] = 0.0;
			}
		}
		while (m != l);
	}

	z++;
	d++;
	e++;
	for (i = 0; i < n; ++i)
		(z[i])++;

	return (1);
}
static int eigen(int n, float **a, float **evec, float *eval)
{
	//Find eigenvalues and eigenvectors of symmetric n by n matrix a

	int i, j, k, success;
	float *d = new float[n]; //tvector_alloc(0, n, float);
	float *e = new float[n]; //tvector_alloc(0, n, float);

	for(i = 0; i < n; i++)
		for(j = 0; j < n; j++)
			evec[i][j] = a[i][j];
	
	tred2(n, evec, eval, e);
	success = tqli(n, evec, eval, e);

	if(!success)
	{
		delete [] d; //tvector_free(d, 0, float);
		delete [] e; //tvector_free(e, 0, float);
		return(success);
	}

	for (j = 0; j < n - 1; j++)
	{
		float emax = eval[j];
		int jmax = j;
		for (k = j + 1; k < n; k++)
		if (emax < eval[k])
		{
	          emax = eval[k];
			jmax = k;
		}
        if (jmax != j)
		{
			float tmp = eval[j]; eval[j] = eval[jmax]; eval[jmax] = tmp; //SWAP(float, eval[j], eval[jmax]);
			for(i = 0; i < n; i++)
			{
              float tmp = evec[i][j]; evec[i][j] = evec[i][jmax]; evec[i][jmax] = tmp; //SWAP(float, evec[i][j], evec[i][jmax]);
			}
		}
	}
	delete [] d; //tvector_free(d, 0, float);
	delete [] e; //tvector_free(e, 0, float);
	return (success);
}
////////////// eigenvector/value finder ends ////////////////
void Mesh::computeAffinityMtrx(int N, int D)
{
	//computes NxN geodesic affinity mtrx b/w samples and sends it to k-D spectral domain via PCAonAffinity()

//	cout << "geodesic affinity mtrx (b/w pairs of samples only) being computed..\n";

	//use Gaussian kernel'ed version of geodesic distances to normalize mtrx elements in (0, 1), i.e. scale invariant,
	//and to decrease the importance of far-away vertices
	float twoSigmaSquare = 2.0f; //kernel(geoDist(p, q)) = e^{- geoDist(p, q) / twoSigmaSquare}
//Jain's paper says sigma = maxGeoDist; since maxGeoDist is always 1 after normalization below, set twoSigmaSquare = 2 directly//twoSigmaSquare = 0.3f; //so that far away pnts, i.e. big geodesics, map to 0.04; if twoSigmaSquare = 2.0f, they map to as low as 0.61; closest pnts, i.e. geoDist = 0, always map to 

	//ready to use normalized dSpanning/geoDist values to compute affinity mtrx
	affinity = new float*[N];
	float sigma = maxGeoDist / 2.0f, a = 2.0f / (sqrt(2*PI) * sigma); //true 0-mean gaussian kernel parameters that are explained below
	for (int i = 0; i < N; i++)
	{
		affinity[i] = new float[N];
		for (int j = 0; j < N; j++)
		{
			if (j < i) //save time
				affinity[i][j] = affinity[j][i]; //symmetric mtrx for undirected graph (copy precomputed result)
			else
			{
				float geoDistUnn = verts[ samples[i] ]->dSpanningUnn[ samples[j] ];
				//unnormalized geo dists range from 0 to say 666; i want to map [0,666] -> [0,1] where values close to 0 map close to 1 and close to 666 map to 0: 0 -> 1, 666 -> ~0
				//this effect is obtained via gaussian kernel; mu is the value i want to favor the most: 0; sigma^2 (variance) is the width of the kernel; values after 3sigma is regarded as 0; so,
				//3sigma = 666 --> sigma = 666/3=222, making variance = 222^2 = 49284; recall 0-mean gausian kernel: (1 / (sqrt(2pi) x sigma)) x e^(-x^2 / (2 x sigma^2)); since my values are +ve,
				//i need the right half of this gaussian and to make the integral 1 in this interval, use 2 x above, i.e. to make the summation of weights 1, use 2 x above; so, kernel works like this:
				affinity[i][j] = a * pow(EXP, -(geoDistUnn*geoDistUnn) / (2 * sigma * sigma)); //[kernel2] [PCAonAffinity() can use the 1st eigvec as well; no, still can't use 'cos i still use a big gaussian width (sigma) and therefore 1st eigvec is still a bad consant vector]
			}
		} //end of j
	} //end of i
//	cout << "\n";
	
	PCAonAffinity(N, D);
}
void Mesh::PCAonAffinity(int N, int D)
{
	//principal component analysis on the affinity matrix (to embed it into D-dimensional spectral domain; to do this, it computes the eigenvals/vecs of affinity matrix, and then selects the D best eigvecs (corresponding to D-largest eigvals); hence D-dimensional spectral domain will be used to represent N-dimensional affinity mtrx we have here

//	cout << "PCA on affinity matrix to get the spectral coords..\n";

//info: eigenvectors corresponding to symmetric positive definite matrices are orthogonal; so, they can be used as basis-vectors of our new D-dimensional space

	// solve for eigenvalues and eigenvectors.
	// eigen values are sorted in descending order, 
	// eigen vectors are sorted in accordance.

	eigen_values = new float[N];
	eigen_vectors = new float*[N]; //n eigenvectors w/ n components each, i.e. n-dimensions, each (in PCA, i'll choose the best/largest k of them)
	for (int i = 0; i < N; i++)
		eigen_vectors[i] = new float[N];
	eigen(N, affinity, eigen_vectors, eigen_values); //call-by-reference effect to fill eigen_* (eigvectors come normalized, i.e. length = 1)

//eigen_values in descending order as we go from [0] to [n]; eigen_vectors.column1 = eigvector1 (corresponding to largest eigval), col2 = eigvec2, and so on..

	//compute 2 matrices that need to be multiplied to reach spectralCoords
	float** diagK = new float*[D]; //DxD
	for (int i = 1; i <= D; i++) //use D eigenvectors starting from the second 'cos 1st is constant
	{
		diagK[i-1] = new float[D];
		for (int j = 1; j <= D; j++)
			diagK[i-1][j-1] = (i == j ? sqrt(eigen_values[i]) : 0.0f); //use sqrt(eigval) as suggested by literature, i.e. for better eigval scaling
	}
	for (int i = 0; i <= D; i++)
		if (eigen_values[i] < 0.0f)
			eigen_values[i] = 0.0f;
//remember: literature uses 2nd, 3rd, and 4th eigvecs (for 3-D)!! ppt: First eigenvector ignored as it is almost a constant 'cos row-sums of affinity mtrx are almost constant

	float** ekt = new float*[D]; //DxN
	for (int i = 1; i <= D; i++) //use D eigenvectors starting from the second 'cos 1st is constant
	{
		ekt[i-1] = new float[N];
		for (int j = 0; j < N; j++)
			ekt[i-1][j] = eigen_vectors[j][i];
	}
	
	//let's get the desired spectral coords
	float** spectralCoords = new float*[D];
	for (int i = 0; i < D; i++) //alloc but not fill yet
		spectralCoords[i] = new float[N]; //kxn mtrx
	mtrxMult(spectralCoords, D, D, D, N, diagK, ekt); //desired result (of size kxn) filled by call-by-ref

	//set the desired spectral coords
	for (int i = 0; i < N; i++)
	{
		if (D <= 6)
		{
			verts[ samples[i] ]->spectralD = new float[D];
			for (int j = 0; j < D; j++)
				verts[ samples[i] ]->spectralD[j] = spectralCoords[j][i];
//note that, for k > 3, quaternion-based icp algo not work; use svd-based algo in that case (literature\others\icp\coding (svd instead of quaternions).m)
		}
		else
		{
			cout << "need some work to enumarate principal axes for 7+ dimension; so, currently can perform spectral embedding up to 6-D\n";
			exit(0);
		}
	}

	//translate spectral coords so that center of mass is at the origin (0,0,0) and scale them so that avg dist to origin is 1
	translateAndScaleSpectralBasesK(N, D);

	//recapture memo
	//delete [] eigen_values; //not here
	//for (int i = 0; i < n; i++) delete [] eigen_vectors[i]; delete [] eigen_vectors; //not here 'cos used in findB..ors()
	for (int i = 0; i < D; i++)
	{
		delete [] diagK[i];
		delete [] ekt[i];
		delete [] spectralCoords[i];
	}
	delete [] diagK;
	delete [] ekt;
	delete [] spectralCoords;

	//cout << "\n";
}
void Mesh::translateAndScaleSpectralBasesK(int N, int D)
{
	//translated all D-D spectral samples s.t. their center at origin (0,0,0) and scale'em s.t. samples' avg dist to origin is 1

	float* m1 = new float[D], * m2 = new float[D]; //mass ~ mean names used interchangably
	float* sum = new float[D];
	for (int i = 0; i < D; i++)
		sum[i] = 0.0f;
	for (int i = 0; i < N; i++)
		for (int c = 0; c < D; c++)
			sum[c] += verts[ samples[i] ]->spectralD[c];
	for (int c = 0; c < D; c++)
		m1[c] = sum[c] / (float) N;
	//initial translation (t) that moves shape centers to (0,0,0,0,0) when D = 5
	for (int i = 0; i < N; i++)
		for (int c = 0; c < D; c++)
			verts[ samples[i] ]->spectralD[c] -= m1[c];
	//now 2 models are optimally translation-aligned
	for (int c = 0; c < D; c++)
		m1[c] = m2[c] = 0.0f; //m1 and m2 moved to center after translation above

	//initial scaling (s) that makes mean variance from origin 1, i.e. scale each model to an appropriate size via avgDistToOrigin
	float avgDistToOrigin = 0.0f;
	for (int i = 0; i < N; i++)
		//taking pow2 below makes sum of means from origin to be 1; w/o pow2 it makes sum of distances from origin 1, which i want
		avgDistToOrigin += distanceBetweenKD(verts[ samples[i] ]->spectralD, m1, D); //m1 = (0,0,0,0) when D = 4
	for (int i = 0; i < N; i++)
		for (int c = 0; c < D; c++)
			verts[ samples[i] ]->spectralD[c] /= avgDistToOrigin;
	
	//now 2 models are optimally scale-aligned by Horn; m1 & m2 cannot change after scaling above; so, no redundant update here
	delete [] m1; delete [] m2; delete [] sum;
}

vector< int > Mesh::denseMatchWith(Mesh* mesh2, int N, vector< int > goodMatches)
{
	//min-weight perfect matching w/ this mesh and mesh2 using coarse goodMatches

	int nRows = N, nCols = nRows;
	int* vIdx1 = new int[N], * vIdx2 = new int[N];
	float** costs = new float*[nRows];
	int nVertNum = nRows*2, nEdgeNum = nRows*nRows, e, * edges = new int[2*nEdgeNum], * matches = new int[nVertNum];
	double* weightsD = new double[nEdgeNum];
	int row = 0;
	for (int i = 0; i < N; i++)
	{
		costs[row] = new float[nCols];
		vIdx1[row] = samples[i];
		int col = 0;
		for (int j = 0; j < N; j++)
			vIdx2[col++] = mesh2->samples[j];
		row++;
	}
	row = 0;
	for (int i = 0; i < N; i++)
	{
		int base1 = samples[i];
		int col = 0;
		for (int j = 0; j < N; j++)
		{
			int base2 = mesh2->samples[j];
			//bool combinatorialMatch = false; for (int y = 0; y < (int) safeCorresp.size(); y += 2) if (safeCorresp[y] == base1 && safeCorresp[y+1] == base2) combinatorialMatch = true;
			float vertexBasedIso = 0.0f;
			int nAdds = 0;
			for (int ss = 0; ss < (int) goodMatches.size(); ss += 2)
			{
				int s1 = goodMatches[ss], s2 = goodMatches[ss + 1];
				//if (s1 == base1) continue; not done here 'cos unlike getCorrespListDistortion() i don't just compute distortion of a given mapping here; rather, i create all possible pairwise distortions using a safe list in safeCorresp which requires no continue here
				float t1s1 = verts[s1]->dSpanning[base1], t2s2 = mesh2->verts[s2]->dSpanning[base2]; //t1s1 = verts[s1]->dSpanningUnn[base1] / maxExtremeDist1, t2s2 = mesh2->verts[s2]->dSpanningUnnSS[base2] / maxExtremeDist2;
				vertexBasedIso += fabs(t1s1 - t2s2);
				nAdds++;				
			}
			vertexBasedIso /= nAdds;
			//combinatorial matches will have vertexBasedIso=0 for sure'cos both t1s1 and t2s2'll be 0 above; so, they're encouraged already
			costs[row][col++] = vertexBasedIso; //(! combinatorialMatch ? vertexBasedIso : -INF); //encourage combinatorial matches (initial ~5 samples in baseVerts2) to be selected by blossom
		}
		row++;
	}
	e = 0;
	for (int i = 0; i < nRows; i++)
		for (int j = 0; j < nCols; j++)
		{
			edges[2*e] = i;
			edges[2*e+1] = nRows+j;
			weightsD[e++] = (double) costs[i][j];
		}
	Blossom* blossom = new Blossom();
	blossom->match(nVertNum, nEdgeNum, edges, weightsD, matches);
	delete blossom;
	vector< int > denseMap; //resulting dense map to be returned; same format in safeCorresp, i.e. 2 consecutive items make 1 matching pair, left from this mesh, right from mesh2
	for (int m = 0; m < nRows; m++)
	{
		int m1 = matches[2*m]; //left side of the (m+1)'th matching line
		int m2 = matches[2*m + 1] - nRows; //mathches[m1 + 1], right side of the line, is in [nRows..2*nRows-1], hence nRows entry at mesh2
		int mi = vIdx2[m2];
		denseMap.push_back(vIdx1[m1]); //base1 from this mesh1; be sure even idxs ~ mesh1, odds ~ mesh2
		denseMap.push_back(mi); //base2 from mesh2
	} //end of m
	//memo re-capture
	delete [] matches;
	delete [] edges;
	delete [] weightsD;
	for (int i = 0; i < nRows; i++) delete [] costs[i];		delete [] costs;
	delete [] vIdx1;
	delete [] vIdx2;

	return denseMap;
}
