#include "Mesh.h"

Mesh::Mesh(int n, int type, int inputDB, int modelID)
{
	//2nd constructor to be called by database models
	sketch = false;
	z = 0.0f;
	N = n;
	int dbType = 2;//1;//2;//2 is very cool for mcgill.snake and spectacle
	if (inputDB == WATERTIGHT)
	{
		if (type == 0) //human class
			frameID = (dbType == 1 ? 9 : 5);//4/*10*/);
		else if (type == 1) //cup class
			frameID = (dbType == 1 ? 21 : 25);
		else if (type == 2) //glasses class
			frameID = (dbType == 1 ? 44 : 52/*54*/);
		else if (type == 3) //airplane class
			frameID = (dbType == 1 ? 63 : 74);
		else if (type == 4) //ant class
			frameID = (dbType == 1 ? 86 : 88);
		else if (type == 5) //chair class
			frameID = (dbType == 1 ? 116 : 104);
		else if (type == 6) //octopus class
			frameID = (dbType == 1 ? 122 : 130);
		else if (type == 7) //table class
			frameID = (dbType == 1 ? 156 : 160);
		else if (type == 8) //teddy class
			frameID = (dbType == 1 ? 167 : 169);
		else if (type == 9) //hand class
			frameID = (dbType == 1 ? 194 : 189);
		else if (type == 10) //plier class
			frameID = (dbType == 1 ? 214 : 212);
		else if (type == 11) //fish class
			frameID = (dbType == 1 ? 221 : 228);
		else if (type == 12) //bird class
			frameID = (dbType == 1 ? 259 : 260);
		else if (type == 13) //armadillo class
			frameID = (dbType == 1 ? 283/*282*/ : 287);
		else if (type == 14) //bust class
			frameID = (dbType == 1 ? 318 : 319);
		else if (type == 15) //mech class
			frameID = (dbType == 1 ? 325 : 328);
		else if (type == 16) //bearing class
			frameID = (dbType == 1 ? 341 : 345);
		else if (type == 17) //vase class
			frameID = (dbType == 1 ? 370 : 373);
		else if (type == 18) //fourleg class
			frameID = (dbType == 1 ? 386/*383*/ : 398/*394*/);//392
		else if (type == 19) //rope class
			frameID = (dbType == 1 ? 280 : 261);
	}
	else if (inputDB == MCGILL)
	{
		if (type == 0) //ant
			frameID = (dbType == 1 ? 7 : 27);
		else if (type == 1) //crab
			frameID = (dbType == 1 ? 32 : 59);
		else if (type == 2) //hand
			frameID = (dbType == 1 ? 67 : 74);
		else if (type == 3) //human
			frameID = (dbType == 1 ? 105 : 108);
		else if (type == 4) //octopus
			frameID = (dbType == 1 ? 113 : 130);
		else if (type == 5) //plier
			frameID = (dbType == 1 ? 136 : 147);
		else if (type == 6) //snake
			frameID = (dbType == 1 ? 158 : 178);
		else if (type == 7) //spectacle
			frameID = (dbType == 1 ? 185 : 202);
		else if (type == 8) //spider
			frameID = (dbType == 1 ? 207 : 231);
		else if (type == 9) //teddy
			frameID = (dbType == 1 ? 238 : 254);
	}
	//TOSCA
	else if (type < 0) //TOSCA cat class (type == -1), TOSCA centaur class (type == -2), ..
		frameID = modelID;
	else
	{
		cout << "undefined type\n";
		exit(0);
	}
}
int Mesh::addVertex(float x, float y, float z)
{
	float v[3] = {x, y, z};
	return addVertex(v);
}
int Mesh::addVertex(float* v)
{
	//add v to verts, if not duplicated

	bool duplicateExists = false;
	int i = (int) verts.size();
	
	//duplicate prevention
	float ROUNDING_ERR = 0.01f;
	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(SbVec3f 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));
}

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

	return triInd;
}

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

	for (unsigned int vie = 0; vie < verts[vi]->edgeList.size(); vie++)
		for (unsigned int vje = 0; vje < 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::computeNormalsAndAreas()
{
	//compute normals of tris and verts in mesh as well as the areas of the tris

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

		if (tris[t]->area < 0) //should never happen unlike the tetrahedron case
		{
			cout << "flipped/inverted tri w/ area = " << tris[t]->area << "\t" << t << "\t" << v1i << " " << v2i << " " << v3i << "\n";
			continue;
		}
		else if (tris[t]->area == 0) //mInv'll not make any sense, i.e. all entries 1.#INF due to div by 0 so discard it in energy computations later
		{
			cout << "degenerate tri w/ area = " << tris[t]->area << "\t" << t << "\t" << v1i << " " << v2i << " " << v3i << "\n";
			continue;
		}
	}

	//vertices
	for (unsigned int v = 0; v < (int) verts.size(); v++)
	{
		SbVec3f normalTotal(0.0f, 0.0f, 0.0f);
		int howManyTriNeighb = 0;
		for (unsigned int tneighb = 0; tneighb < 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
			{
				normalTotal += tris[ verts[v]->triNeighbors[tneighb] ]->normal; //maby a better way is to weight tri-normals w/ their areas for v.normal computation
								//tris[ verts[v]->triNeighbors[tneighb] ]->normal * tris[ verts[v]->triNeighbors[tneighb] ]->area;
				howManyTriNeighb++;
			}

		if (howManyTriNeighb != 0) //if 0, go on w/ the ex-normal
			verts[v]->setNormal(normalTotal / (float) howManyTriNeighb);
		//else
		//	nDegen++;
	}
}

inline float distanceBetween(SbVec3f v1, SbVec3f v2)
{
	//Euclidean distance b/w v1 & v2
	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 distance2Between(SbVec3f v1, SbVec3f v2)
{
	//squared Euclidean distance b/w v1 & v2
	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 distBetween(float* v1, float* v2)
{
	//Euclidean distance b/w v1 & v2
	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 dist2Between(float* v1, float* v2)
{
	//squared Euclidean distance b/w v1 & v2
	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 distanBetween(SbVec3f v1, float* v2)
{
	//Euclidean distance b/w v1 & v2
	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 distanceBetween2D(SbVec2f v1, SbVec2f v2)
{
	//Euclidean distance b/w v1 & v2
	return (float) sqrt( (v2[0]-v1[0])*(v2[0]-v1[0]) + (v2[1]-v1[1])*(v2[1]-v1[1]) );
}
inline float distance2Between2D(SbVec2f v1, SbVec2f v2)
{
	//Euclidean distance b/w v1 & v2
	return (float) (v2[0]-v1[0])*(v2[0]-v1[0]) + (v2[1]-v1[1])*(v2[1]-v1[1]);
}

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

	SbVec3f p_r = v2 - v1;
	SbVec3f q_r = v3 - v1;
	SbVec3f cross = p_r.cross(q_r);

	return 0.5f * cross.length();
}

void Mesh::addEdge(int v1i, int v2i, float dist)
{
	//create an edge b/w v1i and v2i; dist=-1 means this edge is b/w sketch pnts[], otherwise b/w mesh verts[]

	int initial = (int) edges.size();
	Edge* newcomer = new Edge(initial, v1i, v2i, (dist == -1.0f ? distanceBetween(verts[v1i]->coords, verts[v2i]->coords) : dist));

	if (dist == -1.0f)
	{
		verts[v1i]->edgeList.push_back(newcomer->idx);
		verts[v2i]->edgeList.push_back(newcomer->idx); //do same thing for [v2i]
	}
	else
	{
		pnts[v1i]->edgeList.push_back(newcomer->idx);
		pnts[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::getBoundaryTrianglesSharedBy(Edge* bigEdge, int* shTris, int v1i, int v2i)
{
	//writes 2 noninterior/boundary tris adjacent to bigEdge (or to an arbitrary edge defined by v1i and v2i) to shTris[]

	//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;	
	for (unsigned int t1 = 0; t1 < verts[v1i]->triNeighbors.size(); t1++)
	{
		tr1 = verts[v1i]->triNeighbors[t1];
		for (unsigned int t2 = 0; t2 < verts[v2i]->triNeighbors.size(); t2++)
		{
			tr2 = verts[v2i]->triNeighbors[t2];
			if  (tr1 == tr2)
				shTris[a++] = tr1; //or tr2
		}
	}
	//a is 1 if bigEdge is a border edge
	if (a == 1)
	{
		shTris[1] = -1;
		cout << v1i << "," << v2i << " is a border edge found in my closed-mesh :(\n\n"; //exit(0);}
	}
	if (shTris[0] == shTris[1])
	{
		cout << "wrong\n";
		exit(0);
	}
}

#define LINESIZE 1024
void Mesh::loadSketch(int sid)
{
	//load the sketch from xml file w/ existing stroke info; i may segment/subdivide these preknown strokes and may also add virtual strokes

	sketchID = sid;
	char fName[250], line[LINESIZE];
	sprintf(classStr, "human");
/*sprintf(classStr, "hand");
sprintf(classStr, "bird");
sprintf(classStr, "dog");
sprintf(classStr, "glasses");
sprintf(classStr, "fish");
sprintf(classStr, "chair");
sprintf(classStr, "table");
sprintf(classStr, "centaur");//*/
	sprintf(classStr, "ant");
	sprintf(classStr, "crab");
	sprintf(classStr, "octopus");
	sprintf(classStr, "plier");
	sprintf(classStr, "snake");
	sprintf(classStr, "spider");
	sprintf(classStr, "teddy");

	sprintf(fName, "tools\\IUILabUtilities\\DataCollection\\Csharp\\Exe\\subjects\\%s\\%d_%s.xml", classStr, sketchID, classStr);
	cout << "Sketch initializing (to " << fName << " file)...\n";
	FILE* fPtr;
	if (! (fPtr = fopen(fName, "r")))
	{
		cout << fName << " not found\n";
		exit(0);
	}
	int lineNum = 1, nActivePnts = 0;
    while (fgets(line, LINESIZE, fPtr))
	{
        if (lineNum++ < 3)
			continue; //first 2 lines have no valuable info
		//cout << line;exit(0);

		//process the current line
		char* strID, * token = strtok(line, "\""); //delimiter is "
		float x, y, tokNum = 1;
		bool strokeStarting = false;
		while (token != NULL)
		{
			if (strcmp(token, "  <stroke id=") == 0)
			{
				strokeStarting = true;
				break;
			}
			if (tokNum == 2)
			{
				//strID = token; sucks and i do my manual string copy
				strID = new char[(int) strlen(token)+1]; //+1 for null terminator \0 [mandatory]
				for (int c = 0; c < (int) strlen(token); c++)
					strID[c] = token[c];
				strID[(int) strlen(token)] = '\0'; //mandatory
			}
			else if (tokNum == 4)
				x = (float) atof(token);
			else if (tokNum == 6)
				y = (float) atof(token);
			tokNum++;
			token = strtok(NULL, "\"");
		}

		if (strokeStarting)
		{
			//create a new stroke consisting of points fread before
			fgets(line, LINESIZE, fPtr);
			Edge* stroke = new Edge((int) strokes.size(), false);
			int idx = -1;
			while (strcmp(line, "  </stroke>\n") != 0)
			{
				token = strtok(line, "<>/"); //delimiters are < and > and /
				tokNum = 1;
				while (token != NULL)
				{
					if (tokNum++ == 3)
					{
						idx = getIdxByName(token);
						stroke->pnts.push_back(idx);
						pnts[idx]->strokeID = (int) strokes.size();
						if ((int) stroke->pnts.size() > 1)
							stroke->length += distanceBetween2D(pnts[ stroke->pnts[(int) stroke->pnts.size() - 1] ]->coords, pnts[ stroke->pnts[(int) stroke->pnts.size() - 2] ]->coords);
						else if ((int) stroke->pnts.size() == 1)
							stroke->v1i = idx; //first point (idx to pnts[])
					}
					token = strtok(NULL, "<>/");
				}
				fgets(line, LINESIZE, fPtr);
			}
			stroke->v2i = idx; //last point (idx to pnts[])
			pnts[ stroke->v1i ]->strokeList.push_back(stroke->idx);
			pnts[ stroke->v2i ]->strokeList.push_back(stroke->idx);
			strokes.push_back(stroke);
			nActivePnts += (int) stroke->pnts.size();
cout << "strokes" << (int) strokes.size()-1 << " w/ " << (int) stroke->pnts.size() << " points and length = " << stroke->length << endl;
		}
		else
			pnts.push_back(new SketchPnt(strID, SbVec2f(x, y), (int) pnts.size()));
	}
    fclose(fPtr);
	//center of mass computation for better visualization
	com[0] = com[1] = com[2] = 0.0f;
	for (int i = 0; i < (int) pnts.size(); i++)
	{
		//additional reflection around x-axis for sketch canvas and open inventor screen compatibility
		pnts[i]->coords[1] = -pnts[i]->coords[1]; //x-coord remains the same, y changes sign
		for (int c = 0; c < 3; c++)
			com[c] += pnts[i]->coords[c];
	}
	com /= ((float) pnts.size());
	nOriginalStrokes = (int) strokes.size();
	//some pnts[] are not referenced by any strokes[] so use nActivePnts instead of pnts.size()
	cout << nOriginalStrokes << " strokes loaded\ntotal # of sketch pnts = " << (int) pnts.size() << "; (# referred pnts = " << nActivePnts << ")\n\n";
	
	//set initial endpnts[] which may later be extended during segmentStrokes()
	for (int s = 0; s < (int) strokes.size(); s++)
	{
		endpnts.push_back(pnts[ strokes[s]->v1i ]);
		strokes[s]->ev1i = (int) endpnts.size() - 1; //v1i of this strokes[s] as a pointer to endpnts[]
		endpnts.push_back(pnts[ strokes[s]->v2i ]);
		strokes[s]->ev2i = (int) endpnts.size() - 1; //v2i of this strokes[s] as a pointer to endpnts[]
		pnts[ strokes[s]->v1i ]->endpnt = pnts[ strokes[s]->v2i ]->endpnt = true; //info: endpnts[(int)endpnts.size()-2 and -1]->endpnt becomes true after update here
	}




//increase # of sketch pnts by uniformly filling the gap b/w consecutive sketch points w/ new sketch pnts
densifySketch(); //this creates a file w/ dense pnt info and nothing else, i.e. no stroke info




	//for geodesic distance computation over the entire sketch, segment nonrigid strokes into smaller rigid ones separated by joint points endpnts
	segmentStrokes();

	//for geodesic distance computation over the entire sketch, decide which stroke is adjacent to which stroke by introducing virtual strokes
	setStrokeConnectivity();
	
	int nStrokesNonremoved = 0, nAdds = 0;
	for (int s = 0; s < (int) strokes.size(); s++)
		nStrokesNonremoved += (strokes[s]->removed ? 0 : 1);
	cout << nStrokesNonremoved << " final # of strokes\n";

	//determine avg stroke length
	edgeLenTotal = 0.0f;
	minEdgeLen = INF;
	maxEdgeLen = -INF;
	for (int s = 0; s < (int) strokes.size(); s++)
	{
		if (strokes[s]->removed || strokes[s]->virtua)
			continue; //ignore removed strokes (they're kept physically for index consistency); also ignore virtual strokes for avg length computations
		if (strokes[s]->length < minEdgeLen)
			minEdgeLen = strokes[s]->length;
		if (strokes[s]->length > maxEdgeLen)
			maxEdgeLen = strokes[s]->length;
		edgeLenTotal += strokes[s]->length;
		nAdds++;
	}
	avgEdgeLen = edgeLenTotal / nAdds;

	cout << (int) endpnts.size() << " endpnts exist over all strokes\n";
}
void Mesh::densifySketch()
{
	//make sketch denser by adding new sketch pnts b/w distant existing sketch pnts

	float threshold = 1.0f; //existing pnts[] pair p w/ distance(p.v1, p.v2) > threshold are considered distant and i add new sketch pnts on the line connecting p.v1 and p.v2; new pnts are separated by distance threshold, hence uniform sampling
	char fName[250], pName[50];
	sprintf(fName, "tools\\IUILabUtilities\\DataCollection\\Csharp\\Exe\\subjects\\%s\\dense\\%d_pnts[].dat", classStr, sketchID);
	FILE* fPtr = fopen(fName, "w");
	int nNewPnts = 0; //# of new pnts added
	for (int s = 0; s < (int) strokes.size(); s++) //for each stroke
	{
		for (int p = 0; p < (int) strokes[s]->pnts.size()-1; p++) //for each point of this stroke; -1 'cos i use p+1 below
		{
			float dist = distanceBetween2D(pnts[ strokes[s]->pnts[p] ]->coords, pnts[ strokes[s]->pnts[p+1] ]->coords); //distance b/w [p] and [p+1]
			fprintf(fPtr, "%f\t%f\t%d\n", pnts[ strokes[s]->pnts[p] ]->coords[0], pnts[ strokes[s]->pnts[p] ]->coords[1], 1); //last 1 means this is an existing sketch pnt from xml file
			if (dist > threshold)
			{
				SbVec2f dir = pnts[ strokes[s]->pnts[p+1] ]->coords - pnts[ strokes[s]->pnts[p] ]->coords; //from p to p+1
				dir.normalize(); //unit direction				
				float threshold2 = threshold / 2.0f; //make new samples very close to each other by decreasing existing threshold
				int n = (int) (dist / threshold2);
//cout << n << "\n";
				for (float i = 1.0f; i < n; i++)
				{
//cout << (i/n)*dir[0] << " " << (i/n)*dir[1] << " is added to p= " << pnts[ strokes[s]->pnts[p] ]->coords[0] << " " << pnts[ strokes[s]->pnts[p] ]->coords[1] << endl; 
					sprintf(pName, "newPnt%d", nNewPnts++);
					SbVec2f newPnt = pnts[ strokes[s]->pnts[p] ]->coords + (i/n)*dist*dir;
					pnts.push_back(new SketchPnt(pName, newPnt, (int) pnts.size()));
					fprintf(fPtr, "%f\t%f\t%d\n", newPnt[0], newPnt[1], -1); //last -1 means this is a new sketch pnt sampled b/w p and p+1
				}
			}
		}
		fprintf(fPtr, "%f\t%f\t%d\n", pnts[ strokes[s]->pnts[ (int) strokes[s]->pnts.size()-1 ] ]->coords[0], pnts[ strokes[s]->pnts[ (int) strokes[s]->pnts.size()-1 ] ]->coords[1], 1); //last 1 means this is an existing sketch pnt from xml file, namely this is last pnt of this stroke
//exit(0);
	}
	fclose(fPtr);

	cout << "sketch is more dense now (existing + new pnts = " << (int) pnts.size() - nNewPnts << " + " << nNewPnts << " = " << (int) pnts.size() << "\nuse these sketch pnts fprinted at: " << fName << endl;
	exit(0);
}
void Mesh::loadSketchDense_notInUse(int sid)
{
	//load the sketch w/ increased/dense # of pnts from dat file w/o^^ any stroke info; i'll first decide strokes by corner detection and then may segment/subdivide these computed strokes and may also add virtual strokes
	//^^ i'm implicitly using stroke info in presegment=true mode by using consecutiveness; that's why this function no longer in use

	sketchID = sid;
	char fName[250], pName[50];
	sprintf(fName, "tools\\IUILabUtilities\\DataCollection\\Csharp\\Exe\\subjects\\human\\dense\\%d_pnts[].dat", sketchID);
	cout << "Sketch initializing (to " << fName << " file)...\n";
	FILE* fPtr;
	if (! (fPtr = fopen(fName, "r")))
	{
		cout << fName << " not found\n";
		exit(0);
	}
	
	int nPnts = 0, c; //# of all pnts in dat file
	float a, b;
	while (! feof(fPtr))
	{
		fscanf(fPtr, "%f\t%f\t%d\n", &a, &b, &c);
		sprintf(pName, "pnt%d", nPnts++);
		pnts.push_back(new SketchPnt(pName, SbVec2f(a, b), (int) pnts.size()));
	}
	cout << (int) pnts.size() << " sketch pnts loaded\n";
	fclose(fPtr);
	
	//center of mass computation for better visualization
	com[0] = com[1] = com[2] = 0.0f;
	for (int i = 0; i < (int) pnts.size(); i++)
		//pnts[i]->coords[1] = -pnts[i]->coords[1]; //additional reflection around x-axis for sketch canvas and open inventor screen compatibility [no! 'cos %d_pnts[].dat file created by densifySketch() already has the reflected coords]
		for (int c = 0; c < 3; c++)
			com[c] += pnts[i]->coords[c];
	com /= ((float) pnts.size());

	//my auto segmentation on the loaded sketch pnts requires 1+ strokes to segment; so put all pnts[] into strokes[] by initiating a new stroke whenever 2 consecutive pnts are far enough
	bool presegment = true;
	float tolerance = 50.0f;
	int p = 0;
	while (p < (int) pnts.size())
	{
		Edge* stroke = new Edge((int) strokes.size(), false);
		for ( ; p < (int) pnts.size(); p++)
		{
			stroke->pnts.push_back(p);
			pnts[p]->strokeID = (int) strokes.size();
			if ((int) stroke->pnts.size() > 1)
			{
//cout << "dist= " << distanceBetween2D(pnts[ stroke->pnts[(int) stroke->pnts.size() - 1] ]->coords, pnts[ stroke->pnts[(int) stroke->pnts.size() - 2] ]->coords) << endl << endl;
				float consecDist = distanceBetween2D(pnts[ stroke->pnts[(int) stroke->pnts.size() - 1] ]->coords, pnts[ stroke->pnts[(int) stroke->pnts.size() - 2] ]->coords);
				if (consecDist <= tolerance || ! presegment) //if presegmentation not wanted then i increase length, i.e. keep p a part of this segment, no matter what consecDist is
					stroke->length += consecDist;
				else
				{
					//finish this stroke and start a new one by breaking from this inner loop
					//this p is the v1i (beginning) of the next stroke; so remove it from this stroke's pnts[]
					stroke->pnts.erase(stroke->pnts.begin() + (int) stroke->pnts.size()-1); //removes entry [last]
					pnts[p]->strokeID = -1; //redundant 'cos it will be set to the correct value in the next outer loop iteration
					break; //p++ not executed so i'll treat this p again in the next outer loop iteration
				}
			}
			else if ((int) stroke->pnts.size() == 1)
				stroke->v1i = p; //first point (idx to pnts[])
		}
		//finish and add this stroke
		stroke->v2i = p - 1; //last point (idx to pnts[])
		pnts[ stroke->v1i ]->strokeList.push_back(stroke->idx);
		pnts[ stroke->v2i ]->strokeList.push_back(stroke->idx);
		strokes.push_back(stroke);
		cout << "strokes" << (int) strokes.size()-1 << " w/ " << (int) stroke->pnts.size() << " points and length = " << stroke->length << endl;
	}
	cout << (int) strokes.size() << " strokes after pre-segmentation\n\n";

	//set initial endpnts[] which may later be extended during segmentStrokes()
	for (int s = 0; s < (int) strokes.size(); s++)
	{
		endpnts.push_back(pnts[ strokes[s]->v1i ]);
		strokes[s]->ev1i = (int) endpnts.size() - 1; //v1i of this strokes[s] as a pointer to endpnts[]
		endpnts.push_back(pnts[ strokes[s]->v2i ]);
		strokes[s]->ev2i = (int) endpnts.size() - 1; //v2i of this strokes[s] as a pointer to endpnts[]
		pnts[ strokes[s]->v1i ]->endpnt = pnts[ strokes[s]->v2i ]->endpnt = true; //info: endpnts[(int)endpnts.size()-2 and -1]->endpnt becomes true after update here
	}

	//my auto segmentation on the loaded sketch pnts 
	segmentStrokes(); //for geodesic distance computation over the entire sketch, segment nonrigid strokes into smaller rigid ones separated by joint points endpnts

	//for geodesic distance computation over the entire sketch, decide which stroke is adjacent to which stroke by introducing virtual strokes
	setStrokeConnectivity();
	
	int nStrokesNonremoved = 0, nAdds = 0;
	for (int s = 0; s < (int) strokes.size(); s++)
		nStrokesNonremoved += (strokes[s]->removed ? 0 : 1);
	cout << nStrokesNonremoved << " final # of strokes\n";

	//determine avg stroke length
	edgeLenTotal = 0.0f;
	minEdgeLen = INF;
	maxEdgeLen = -INF;
	for (int s = 0; s < (int) strokes.size(); s++)
	{
		if (strokes[s]->removed || strokes[s]->virtua)
			continue; //ignore removed strokes (they're kept physically for index consistency); also ignore virtual strokes for avg length computations
		if (strokes[s]->length < minEdgeLen)
			minEdgeLen = strokes[s]->length;
		if (strokes[s]->length > maxEdgeLen)
			maxEdgeLen = strokes[s]->length;
		edgeLenTotal += strokes[s]->length;
		nAdds++;
	}
	avgEdgeLen = edgeLenTotal / nAdds;

	cout << (int) endpnts.size() << " endpnts exist over all strokes\n";
}
void Mesh::loadSketchDense()
{
	//load the sketch w/ increased/dense # of pnts from dat file w/o any stroke info; no stroke decision is necessary; each sketch pnt is connected to close pnts by edges; these edges are used in geodesic computations

	bool continuationBasedEdgeCreation = false, //true: create edges to make good continuations [not in use]; false: create edges b/w close-enough pnts
		 setPntConnectivity = false; //true: sets neighborhood of each pnt in this function; false: just loads pnts w/o any connectivity computation and later setSampleConnetivity() handles topology
	char fName[250], pName[50];
	sprintf(classStr, "human");
//sprintf(classStr, "hand");
//sprintf(classStr, "bird");
//sprintf(classStr, "dog");
//sprintf(classStr, "glasses");
//sprintf(classStr, "fish");
//sprintf(classStr, "chair");
//sprintf(classStr, "table");
//sprintf(classStr, "centaur");
	sprintf(classStr, "ant");
	sprintf(classStr, "crab");
	sprintf(classStr, "hand");
	sprintf(classStr, "human");
	sprintf(classStr, "octopus");
	sprintf(classStr, "plier");
	sprintf(classStr, "snake");
//	sprintf(classStr, "glasses");
//	sprintf(classStr, "spider");
//	sprintf(classStr, "teddy");

	if (strcmp(classStr, "human") == 0)
		sketchID = 14;//12;//20;//12;//13;//14;//4;//0;//8, 2
	else if (strcmp(classStr, "hand") == 0)
		sketchID = 5;//7;//6;
	else if (strcmp(classStr, "bird") == 0)
		sketchID = 11;//10;//1;//2;0;
	else if (strcmp(classStr, "dog") == 0)
		sketchID = 7;//3;//1;//2;//0;4;5;
	else if (strcmp(classStr, "glasses") == 0)
		sketchID = 0;//2;//3;//2;//1;//0;
	else if (strcmp(classStr, "fish") == 0)
		sketchID = 5;//6;//0;//4;//3;//2;//1;//0;
	else if (strcmp(classStr, "chair") == 0)
		sketchID = 6;//3;//1;
	else if (strcmp(classStr, "table") == 0)
		sketchID = 0;//1;
	else if (strcmp(classStr, "centaur") == 0)
		sketchID = 0;//1;//2;
	else if (strcmp(classStr, "ant") == 0)
		sketchID = 1;//0;//2;
	else if (strcmp(classStr, "crab") == 0)
		sketchID = 1;//0;//2;
	else if (strcmp(classStr, "octopus") == 0)
		sketchID = 2;//1;//0;
	else if (strcmp(classStr, "plier") == 0)
		sketchID = 1;//2;//0;
	else if (strcmp(classStr, "snake") == 0)
		sketchID = 2;//0;//1;
	else if (strcmp(classStr, "spider") == 0)
		sketchID = 0;//1;//2;
	else if (strcmp(classStr, "teddy") == 0)
		sketchID = 8;//5;//7;//5;//6;//3;//4;//1;//2;//0;
	else
	{
		cout << "undefined class\n";
		exit(0);
	}

	sprintf(fName, "tools\\IUILabUtilities\\DataCollection\\Csharp\\Exe\\subjects\\%s\\dense\\%d_pnts[].dat", classStr, sketchID); //use %d_pnts[] coarse.dat to see the original sparse set of pnts
	cout << "Sketch initializing (to " << fName << " file)...\n";
	FILE* fPtr;
	if (! (fPtr = fopen(fName, "r")))
	{
		cout << fName << " not found\n";
		exit(0);
	}
	
	int nPnts = 0, c; //# of all pnts in dat file
	float a, b;
	while (! feof(fPtr))
	{
		fscanf(fPtr, "%f\t%f\t%d\n", &a, &b, &c);
		sprintf(pName, "pnt%d", nPnts++);
		pnts.push_back(new SketchPnt(pName, SbVec2f(a, b), (int) pnts.size()));
		pnts[(int) pnts.size() - 1]->processed = false; //true means this is incident to an edge
	}
	cout << (int) pnts.size() << " sketch pnts loaded\n";
	fclose(fPtr);
	
	//center of mass computation for better visualization
	com[0] = com[1] = com[2] = 0.0f;
	for (int i = 0; i < (int) pnts.size(); i++)
		//pnts[i]->coords[1] = -pnts[i]->coords[1]; //additional reflection around x-axis for sketch canvas and open inventor screen compatibility [no! 'cos %d_pnts[].dat file created by densifySketch() already has the reflected coords]
		for (int c = 0; c < 3; c++)
			com[c] += pnts[i]->coords[c];
	com /= ((float) pnts.size());
	//bounding box computation for spacing distance computation to be used in sampling
	float minX, minY, maxX, maxY, w, h;
	for (int p = 0; p < (int) pnts.size(); p++)
		if (p == 0)
		{
			minX = maxX = pnts[p]->coords[0];	minY = maxY = pnts[p]->coords[1];
		}
		else
		{
			if (pnts[p]->coords[0] < minX)	minX = pnts[p]->coords[0]; else if (pnts[p]->coords[0] > maxX)	maxX = pnts[p]->coords[0];
			if (pnts[p]->coords[1] < minY)	minY = pnts[p]->coords[1]; else if (pnts[p]->coords[1] > maxY)	maxY = pnts[p]->coords[1];
//if (pnts[p]->coords[0] > 315 && pnts[p]->coords[0] < 350) cout << p << "p\n";
		}
	w = maxX-minX;
	h = maxY-minY;
	spacing = sqrt(w*w + h*h) / 40.0f; //diagonal of bbox divided by a constant factor, say 40
	cout << "sketch bbox.(w,h): (" << w << ", " << h << ")\tspacing = " << spacing << "\n";
//cout << "sketch bbox.(w,h): (" << minX << ", " << maxX << ")\tspacing = " << spacing << "\n";
	if (! setPntConnectivity)
	{
		cout << "\n";
		return;
	}

	//connectivity of the loaded sketch pnt
	float maxEdgeLenAllowed = (sketchID > 3 ? 8.0f : 2.0f), //2 pnts closer than this value will (continuationBasedEdgeCreation=false) or may (continuationBasedEdgeCreation=true) have an edge in between
		  maxRatioAllowed = 2.0f; //extension is a good continuation if generates a linenessRatio below this threshold
	if (continuationBasedEdgeCreation)
	{
		//select the first neighbor of p such that i can define a line (continuation start point) from p to that neighbor
		int p = 0, clStartPnt = p; //properties of continuation line (cl), i.e. the starting pnt of this cl
		float clPath; //properties of continuation line (cl), i.e. the path/geodesic distance accumulated on it
		////////////////////////// initialization ////////////////////////////////////
		for (; p < (int) pnts.size(); p++)
		{
			int closestNghb = -1;
			float minDist = INF;
			for (int r = 0; r < (int) pnts.size(); r++)
				if (p != r && ! pnts[r]->processed) //no edge (of length 0) b/w the same pnts p=r
				{
					float dist = distanceBetween2D(pnts[p]->coords, pnts[r]->coords);
					if (dist < minDist)
					{
						minDist = dist;
						closestNghb = r;
					}
				}
			//create a new edge b/w p and closestNghb
			if (pnts[p]->addPntNeighbor(closestNghb)) //addVertNeighbor returns true if verts[v1i] is not neighbor of [v2i]; so they can form a new edge
				addEdge(p, closestNghb, minDist);
			pnts[closestNghb]->addPntNeighbor(p);	//p closestNghb edge is (probably) added just-above; do not add closestNghb p edge!!!
			pnts[closestNghb]->processed = true; //closestNghb is incident to an edge now
			if (minDist == 0.0f) //p & closestNghb overlap; so no line defined; hence select the next closestNghb for this p
				p--;
			else
			{
				//continuation line (cl) is now defined w/ the following properties
				clPath = minDist;
				break;
			}
		}
		////////////////////////// extension that considers candidates within maxEdgeLenAllowed distance ////////////////////////////////////
		for (; p < (int) pnts.size(); p++)
		{
			//verts close to p are stored along with their linenessRatio acccording to the current continuation line (cl); 3 consecutive items are [i] = close pnt, [i+1] = distance, [i+2] = ratio
			vector< float > candids;
			float minRatio, theDist;
			int minRatioIdx;
			for (int r = 0; r < (int) pnts.size(); r++)
				if (p != r && ! pnts[r]->processed) //no edge (of length 0) b/w the same pnts p=r
				{
					float dist = distanceBetween2D(pnts[p]->coords, pnts[r]->coords);
					if (dist <= maxEdgeLenAllowed)
					{
						float pathLen = clPath + dist, //geodesic length consisting of 2 distances: accumulated clPath distance from clStartPnt to here + current dist
							  eucLen = distanceBetween2D(pnts[clStartPnt]->coords, pnts[r]->coords), //euclidean distance between clStartPnt and here
							  linenessRatio = (eucLen != 0.0f ? pathLen / eucLen : 0.0f); //pathLen / eucLen ratio is >= 1 for sure
						if (linenessRatio < maxRatioAllowed)
						{
							candids.push_back((float) r);
							candids.push_back(dist);
							candids.push_back(linenessRatio);
							if (linenessRatio < minRatio)
							{
								minRatio = linenessRatio;
								minRatioIdx = (int) candids.size() - 3;
								theDist = dist;
							}
						}
					}
				}
			//check all candidates to see whether there is a closer pnt w/ a similar linenessRatio; use that pnt as the neighbor of p, if exists
			float threshold = 0.1f;
			for (int c = 0; c < (int) candids.size(); c += 3)
				if (candids[c+2] >= minRatio-threshold && candids[c+2] <= minRatio+threshold && //current ratio is similar to minRatio
					candids[c+1] <= theDist) //distance from p to this candidate is smaller than theDist
				{
					//update the neighbor of p to this current pnt
					minRatio = candids[c+2];
					minRatioIdx = c;
					theDist = candids[c+1];
				}
			//finally create a new edge b/w p and best candidate
			if (pnts[p]->addPntNeighbor((int) candids[minRatioIdx])) //addVertNeighbor returns true if verts[v1i] is not neighbor of [v2i]; so they can form a new edge
				addEdge(p, (int) candids[minRatioIdx], candids[minRatioIdx+1]);
			pnts[ (int) candids[minRatioIdx] ]->addPntNeighbor(p);	//p bestCandid edge is (probably) added just-above; do not add bestCandid p edge!!!
		}
	}
	else //each pnt generates new edges to its closest pnts
		for (int p = 0; p < (int) pnts.size(); p++)
			for (int r = 0; r < (int) pnts.size(); r++)
				if (p != r) //no edge (of length 0) b/w the same pnts p=r
				{
					float dist = distanceBetween2D(pnts[p]->coords, pnts[r]->coords);
//if (p == 159 && r == 160) cout << dist << " euc\n"; //see euclidean distance b/w these 2 pnts
					if (dist <= maxEdgeLenAllowed)
					{
						//create a new edge b/w p and r				
						if (pnts[p]->addPntNeighbor(r)) //addVertNeighbor returns true if verts[v1i] is not neighbor of [v2i]; so they can form a new edge
							addEdge(p, r, dist);
						pnts[r]->addPntNeighbor(p);	//p r edge is (probably) added just-above; do not add r p edge!!!
					}
				}
	avgEdgeLen = edgeLenTotal / (int) edges.size();
#ifdef PRINT_DETAILS
	cout << (int) edges.size() << " edges on this sketch w/ an average length of " << avgEdgeLen << "\n\n";
#endif
}

void Mesh::loadObj(int fid)
{
	//loads the mesh in obj format

	frameID = fid;
	char meshFile[250];
	sprintf(meshFile, "input\\%d.obj", frameID);
	cout << "Mesh initializing (to " << meshFile << " file)...\n";

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

	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
	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')
			//none of verts are removed/duplicated initially 'cos i read .obj
			addVertexND(SbVec3f(a, b, c)); //ND: no duplicate check
		else //if (type == 'f')
			addTriangle((int) a - 1, (int) b - 1, (int) c - 1); //-1: obj indices start at 1
	avgEdgeLen = edgeLenTotal/((int) edges.size());
	color = SbVec3f(0.75, 0.75f, 0.75f);
}
void Mesh::loadOff(int inputDB, int type)
{
	//loads the mesh in off format

	//frameID is already set in the constructor
	char meshFile[250];
	if (inputDB == TOSCA)
	{
		char typeStr[50];
		sprintf(typeStr, "cat");
		if (type == 1)		sprintf(typeStr, "centaur"); //TOSCA centaur
		else if (type == 2)	sprintf(typeStr, "david"); //TOSCA david
		else if (type == 3)	sprintf(typeStr, "dog"); //TOSCA dog
		else if (type == 4)	sprintf(typeStr, "gorilla"); //TOSCA gorilla
		else if (type == 5)	sprintf(typeStr, "horse"); //TOSCA horse
		else if (type == 6)	sprintf(typeStr, "victoria"); //TOSCA victoria
		else if (type == 7)	sprintf(typeStr, "wolf"); //TOSCA wolf
		sprintf(meshFile, "..\\..\\..\\PhD\\code\\input\\toscahires-asci\\off\\%s%d.off", typeStr, frameID);
	}
	else if (inputDB == MCGILL)
	{
		char typeStr[50];
		sprintf(typeStr, "ant");
		int startID = 0;
		if (type == 1) {		sprintf(typeStr, "crab"); startID = 30; }
		else if (type == 2)	{	sprintf(typeStr, "hand"); startID = 60; }
		else if (type == 3)	{	sprintf(typeStr, "human"); startID = 80; }
		else if (type == 4)	{	sprintf(typeStr, "octopus"); startID = 109; }
		else if (type == 5)	{	sprintf(typeStr, "plier"); startID = 134; }
		else if (type == 6)	{	sprintf(typeStr, "snake"); startID = 154; }
		else if (type == 7)	{	sprintf(typeStr, "spectacle"); startID = 179; }
		else if (type == 8)	{	sprintf(typeStr, "spider"); startID = 204; }
		else if (type == 9)	{	sprintf(typeStr, "teddy"); startID = 235; }
		sprintf(meshFile, "..\\..\\..\\PhD\\code\\input\\McGill\\new_McGill_benchmark\\McGillDBAll\\%s\\b%d.off", typeStr, frameID);
	}
	else //watertight
	{
	sprintf(meshFile, "..\\..\\..\\PhD\\code\\input\\MeshsegBenchmark-1.0\\MeshsegBenchmark-1.0 (everything)\\data\\off\\%d.off", frameID);
/*sprintf(meshFile, "..\\..\\..\\PhD\\code\\input\\psb_v1\\benchmark\\db\\5\\m586\\m586.off"); //glasses from psb benchmark
sprintf(meshFile, "..\\..\\..\\PhD\\code\\input\\psb_v1\\benchmark\\db\\1\\m156\\m156.off"); //human from psb benchmark
sprintf(meshFile, "..\\..\\..\\PhD\\code\\input\\psb_v1\\benchmark\\db\\0\\m67\\m67.off"); //fish from psb benchmark
sprintf(meshFile, "..\\..\\..\\PhD\\code\\input\\psb_v1\\benchmark\\db\\17\\m1762\\m1762.off"); //phone from psb benchmark*/
	}

	//cout << "Mesh initializing (to " << frameID << ".off (short name) file)..\n";
	cout << "Mesh initializing (to " << meshFile << ")..\n";

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

	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;

	minEdgeLen = INF;
	maxEdgeLen = -INF;
	edgeLenTotal = 0.0f;
	while (v++ < nVerts) //go till the end of coords section
	{
		fscanf(fPtr, "%f %f %f\n", &a, &b, &c);
		//none of verts are removed/duplicated initially 'cos i read .obj
		addVertexND(SbVec3f(a, b, c)); //ND: no duplicate check
	}

	//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
	avgEdgeLen = edgeLenTotal/((int) edges.size());
	color = SbVec3f(0.75, 0.75f, 0.75f);
	additionalRotations();
//	cout << "nVerts = " << (int) verts.size() << endl;
}
void Mesh::loadPly(int fid)
{
	//loads the mesh in ply format

	frameID = fid;
	char meshFile[250];
	sprintf(meshFile, "input\\SCAPE\\%d.ply", frameID);
	cout << "Mesh initializing (to " << 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", &ply);
	if (strcmp(ply, "comment") == 0) //supporting 0 or 1 line comments only
	{
		fscanf(fPtr, "%s %s\n", &ply, &nVertsStr); //skip "comment VCGLIB generated" at next line
		fscanf(fPtr, "%s %s %s\n", &ply, &ply, &nVertsStr); //skip "element vertex 12500" at next line
	}
	else
		fscanf(fPtr, "%s %s\n", &ply, &nVertsStr); //skip "element vertex 12500" at next line; element already fread
	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
	minEdgeLen = INF;
	maxEdgeLen = -INF;
	edgeLenTotal = 0.0f;
	while (v++ < nVerts) //go till the end of coords section
	{
		fscanf(fPtr, "%f %f %f\n", &a, &b, &c);
		//none of verts are removed/duplicated initially 'cos i read .obj
		addVertexND(SbVec3f(a, b, c)); //ND: no duplicate check
	}

	//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);
	avgEdgeLen = edgeLenTotal/((int) edges.size());
	color = SbVec3f(0.75, 0.75f, 0.75f);
//	additionalRotations();
}
void Mesh::computeCenterOfMass()
{
	//center of mass of this mesh

	com[0] = com[1] = com[2] = 0.0f;
	for (int v = 0; v < (int) verts.size(); v++)
	{						
		com[0] += verts[v]->coords[0];
		com[1] += verts[v]->coords[1];
		com[2] += verts[v]->coords[2];
	}
	com /= ((float) verts.size());
	//cout << "center-of-mass: " << com[0] << ", " << com[1] << ", " << com[2] << " for mesh" << frameID << "\n";
}

int Mesh::getIdxByName(char* name)
{
	//returns the idx of the sketch point (pnts[idx]) identified by name

	for (int idx = 0; idx < (int) pnts.size(); idx++)
		if (strcmp(pnts[idx]->strID, name) == 0)
			return idx;
cout<<"shouldn't reach here\n";
	return -1;
}

float minOf(float a, float b, float c, float d)
{
	//returns the min of the parameters

	if (a <= b && a <= c && a <= d)
		return a;
	if (b <= a && b <= c && b <= d)
		return b;
	if (c <= a && c <= b && c <= d)
		return c;
	return d;
}

float nCalls = 0.0f;
void Mesh::setStrokeConnectivity()
{
	//sets the neighborhood info of the strokes of this sketch; introduce virtual strokes of ~0 length to connect close/neighbor strokes

	float close = (1.0f + nCalls) * spacing; //strokes' endpoint-to-endpnt distance tolerance; 2 strokes to be connected should be close enough to prevent weird connections b/w strokes that would mess up pairwise distance computation
	nCalls += ((strcmp(classStr, "human") == 0 && sketchID == 14) || (strcmp(classStr, "dog") == 0 && sketchID == 5) || (strcmp(classStr, "hand") == 0 && sketchID == 7) || strcmp(classStr, "fish") == 0 || strcmp(classStr, "chair") == 0 || strcmp(classStr, "table") || strcmp(classStr, "teddy") == 0 ? 2.5f : 0.5f); //for the next call


#ifdef PRINT_DETAILS
	cout << "closeness threshold for virtual stroke creation: " << close <<"\n";
#endif
	//remove virtual strokes of (possible) previous calls
	for (int s = 0; s < (int) strokes.size(); s++)
		if (strokes[s]->virtua)
			strokes[s]->removed = true;
	int initialStrokeSize = (int) strokes.size(); //use this one in for-loop headers as the strokes[] will be appended by new virtual strokes
	for (int s = 0; s < initialStrokeSize; s++) //for each stroke
	{
		if (strokes[s]->removed)
			continue; //ignore removed strokes (they're kept physically for index consistency)
		
		SbVec2f sv1 = pnts[ strokes[s]->v1i ]->coords, sv2 = pnts[ strokes[s]->v2i ]->coords;
		for (int o = s+1; o < initialStrokeSize; o++) //test all other strokes as the neighbor of s; +1 not to test stroke s vs. s and o=s+1 not to test 2 vs. 1 since 1 vs. 2 already tested before for sure
		{
			if (strokes[o]->removed)
				continue; //ignore removed strokes (they're kept physically for index consistency)

			//amongst 4 distances b/w sv1/2 and ov1/2 pick the smallest one 
			SbVec2f ov1 = pnts[ strokes[o]->v1i ]->coords, ov2 = pnts[ strokes[o]->v2i ]->coords;

			float d1 = distanceBetween2D(sv1, ov1), d2 = distanceBetween2D(sv1, ov2), d3 = distanceBetween2D(sv2, ov1), d4 = distanceBetween2D(sv2, ov2);
//cout << strokes[s]->idx << "\t" << strokes[o]->idx << "\t" << d1 << " " << d2 << " " << d3 << " " << d4 << endl;
			//virtual stroke consisting of 2 points and a short (~0) distance b/w them			
			if (d1 <= close && minOf(d1, d2, d3, d4) == d1) //sv1 & ov1 are closest and they are close enough to make strokes s and o neighbors
			{
				Edge* virtualStroke = new Edge((int) strokes.size(), true);
				virtualStroke->pnts.push_back(strokes[s]->v1i);
				virtualStroke->pnts.push_back(strokes[o]->v1i);
				virtualStroke->length = d1;
				virtualStroke->v1i = strokes[s]->v1i; //first point
				virtualStroke->v2i = strokes[o]->v1i; //second/last point
				virtualStroke->ev1i = getIdxToEndpnts(virtualStroke->v1i); //v1i of this virtualStroke as a pointer to endpnts[]
				virtualStroke->ev2i = getIdxToEndpnts(virtualStroke->v2i); //v2i of this virtualStroke as a pointer to endpnts[]
				pnts[ virtualStroke->v1i ]->strokeList.push_back(virtualStroke->idx);
				pnts[ virtualStroke->v2i ]->strokeList.push_back(virtualStroke->idx);
				strokes.push_back(virtualStroke);
#ifdef PRINT_DETAILS
				cout << "strokes" << (int) strokes.size()-1 << " is virtual b/w pnts" << virtualStroke->v1i << " & " << virtualStroke->v2i << " (strokes" << s << " & " << o << ") w/ length = " << virtualStroke->length << "\n";
#endif
			}
			else if (d2 <= close && minOf(d1, d2, d3, d4) == d2) //sv1 & ov2 are closest and they are close enough to make strokes s and o neighbors
			{
				Edge* virtualStroke = new Edge((int) strokes.size(), true);
				virtualStroke->pnts.push_back(strokes[s]->v1i);
				virtualStroke->pnts.push_back(strokes[o]->v2i);
				virtualStroke->length = d2;
				virtualStroke->v1i = strokes[s]->v1i; //first point
				virtualStroke->v2i = strokes[o]->v2i; //second/last point
				virtualStroke->ev1i = getIdxToEndpnts(virtualStroke->v1i); //v1i of this virtualStroke as a pointer to endpnts[]
				virtualStroke->ev2i = getIdxToEndpnts(virtualStroke->v2i); //v2i of this virtualStroke as a pointer to endpnts[]
				pnts[ virtualStroke->v1i ]->strokeList.push_back(virtualStroke->idx);
				pnts[ virtualStroke->v2i ]->strokeList.push_back(virtualStroke->idx);
				strokes.push_back(virtualStroke);
#ifdef PRINT_DETAILS
				cout << "strokes" << (int) strokes.size()-1 << " is virtual b/w pnts" << virtualStroke->v1i << " & " << virtualStroke->v2i << " (strokes" << s << " & " << o << ") w/ length = " << virtualStroke->length << "\n";
#endif
			}
			else if (d3 <= close && minOf(d1, d2, d3, d4) == d3) //sv2 & ov1 are closest and they are close enough to make strokes s and o neighbors
			{
				Edge* virtualStroke = new Edge((int) strokes.size(), true);
				virtualStroke->pnts.push_back(strokes[s]->v2i);
				virtualStroke->pnts.push_back(strokes[o]->v1i);
				virtualStroke->length = d3;
				virtualStroke->v1i = strokes[s]->v2i; //first point
				virtualStroke->v2i = strokes[o]->v1i; //second/last point
				virtualStroke->ev1i = getIdxToEndpnts(virtualStroke->v1i); //v1i of this virtualStroke as a pointer to endpnts[]
				virtualStroke->ev2i = getIdxToEndpnts(virtualStroke->v2i); //v2i of this virtualStroke as a pointer to endpnts[]
				pnts[ virtualStroke->v1i ]->strokeList.push_back(virtualStroke->idx);
				pnts[ virtualStroke->v2i ]->strokeList.push_back(virtualStroke->idx);
				strokes.push_back(virtualStroke);
#ifdef PRINT_DETAILS
				cout << "strokes" << (int) strokes.size()-1 << " is virtual b/w pnts" << virtualStroke->v1i << " & " << virtualStroke->v2i << " (strokes" << s << " & " << o << ") w/ length = " << virtualStroke->length << "\n";
#endif
			}
			else if (d4 <= close && minOf(d1, d2, d3, d4) == d4) //sv2 & ov2 are closest and they are close enough to make strokes s and o neighbors
			{
				Edge* virtualStroke = new Edge((int) strokes.size(), true);
				virtualStroke->pnts.push_back(strokes[s]->v2i);
				virtualStroke->pnts.push_back(strokes[o]->v2i);
				virtualStroke->length = d4;
				virtualStroke->v1i = strokes[s]->v2i; //first point
				virtualStroke->v2i = strokes[o]->v2i; //second/last point
				virtualStroke->ev1i = getIdxToEndpnts(virtualStroke->v1i); //v1i of this virtualStroke as a pointer to endpnts[]
				virtualStroke->ev2i = getIdxToEndpnts(virtualStroke->v2i); //v2i of this virtualStroke as a pointer to endpnts[]
				pnts[ virtualStroke->v1i ]->strokeList.push_back(virtualStroke->idx);
				pnts[ virtualStroke->v2i ]->strokeList.push_back(virtualStroke->idx);
				strokes.push_back(virtualStroke);
#ifdef PRINT_DETAILS
				cout << "strokes" << (int) strokes.size()-1 << " is virtual b/w pnts" << virtualStroke->v1i << " & " << virtualStroke->v2i << " (strokes" << s << " & " << o << ") w/ length = " << virtualStroke->length << "\n";
#endif
			}
		}
	}
	int nStrokes = 0, nViStrokes = 0, nEnds = 0;
	for (int s = 0; s < (int) strokes.size(); s++)
	{
		nStrokes += (strokes[s]->removed ? 0 : 1);
		nViStrokes += (! strokes[s]->removed && strokes[s]->virtua ? 1 : 0);
	}
	for (int p = 0; p < (int) pnts.size(); p++)
		nEnds += (pnts[p]->endpnt ? 1 : 0);
#ifdef PRINT_DETAILS
	cout << (nStrokes-nViStrokes) << " main + " << nViStrokes << " virtual = " << nStrokes << " strokes w/ " << nEnds << " endpnts total\n\n";
#endif
//	if (nViStrokes > 3*(nStrokes-nViStrokes))
//		cout << "\n\n!!!!!!!!WARNING!!!!!!!! too many virtual strokes; try reducing `close' variable\n\n";
}

void Mesh::segmentStrokes()
{
	//segments nonrigid strokes into smaller rigid ones separated by joint points endpnts

/*	//i started implementing my own segmentation idea but then switched to well-known douglas-peucker algo; so canceling this block
	float farr = 10.0f, dist, maxDist = -INF;
	int initialStrokeSize = (int) strokes.size(), tmp; //use this one in for-loop headers as the strokes[] will be appended by new strokes
	for (int s = 0; s < initialStrokeSize; s++) //for each stroke
		for (int p = 1; p < (int) strokes[s]->pnts.size()-1; p++) //-1 'cos i access pnts[p-1, p, and p+1] so break early when p = size-1
		{
			SketchPnt* pLeft = pnts[ strokes[s]->pnts[p-1] ], * pMid = pnts[ strokes[s]->pnts[p] ], * pRight = pnts[ strokes[s]->pnts[p+1] ];
			//pMid is a joint point if its distance to the line defined by pLeft and pRight is big enough
			SbVec2f dir = pRight->coords - pLeft->coords; //direction of the line from pLeft to pRight
			dir.normalize(); //need unit vector
			//line equation is x = pLeft + t*dir; test point is pMid; distance b/w this line and pMid is (from wikipedia)
			dist = ((pLeft->coords - pMid->coords) - (pLeft->coords - pMid->coords).dot(dir) * dir).length();
cout << dist << "**" << strokes[s]->pnts[p] << "\t\t";
			if (dist > farr)
			{
				cout << "dist from pnts" << strokes[s]->pnts[p] << " to the line = " << dist << " at stroke" << s << endl;
			}
			if (dist > maxDist)
			{
				maxDist = dist;
				tmp = strokes[s]->pnts[p];
			}
//break;
		}
	cout << "\nmaxDist = " << maxDist << " for pnts" << tmp << endl;*/

	int initialStrokeSize = (int) strokes.size(); //use this one in for-loop headers as the strokes[] will be appended by new strokes
	for (int s = 0; s < initialStrokeSize; s++) //for each stroke
	{
		vector< double > polyline; //original polyline, assume not empty
		for (int p = 0; p < (int) strokes[s]->pnts.size(); p++)
		{
			//prepare polyline for douglas-peucker algo
			polyline.push_back(pnts[ strokes[s]->pnts[p] ]->coords[0]);
			polyline.push_back(pnts[ strokes[s]->pnts[p] ]->coords[1]);
		}
		////////////////////// douglas-peucker line simplification algo, which acts as a corner-detector on my input stroke(s) //////////////////////
		double tolerance = 10.0;   //point-to-edge distance tolerance
		double* result = new double [ (int) polyline.size() ]; //make sure the result array is large enough
		int nCorners = 0;

		//simplify the 2d polyline [see tools\doug-peucker simplification\ w/ a better code and explanation]
		psimpl::simplify_douglas_peucker<2>(polyline.begin(), polyline.end(), tolerance, result, nCorners); //result[] has the desired corner coords via call-by-ref; nCorners keeps # corners detected

/*		cout<<"\ninputttttttt\n";
		for (int v = 0; v < (int) polyline.size(); v += 2)
			cout << polyline[v] << ", " << polyline[v+1] << "\n";
		cout<<"outtttttttt\tnCorners = " << nCorners << "\n";
		for (int v = 0; v < nCorners; v++)
			cout << result[2*v] << ", " << result[2*v+1] << "\n";//*/

		//process the output, which is stored in result[] as x1, y1, x2, y2, x3, y3, .. where 2 consecutive 2d coords make 1 corner vertex in 2D
		for (int c = 0; c < nCorners; c++)
		{
//cout << c << " &&&&&&&&& " << nCorners << endl; 
			if (c > 0 && c < nCorners-1) //first and last corner are the endpnts of this stroke for sure; so do nothing for them
			{
				/////////////////////////////// creation of 2 substrokes /////////////////////////////
				//result[v], result[v+1] is the next new corner (nc), which is an existing SketchPoint; find this SketchPoint of interest
				SketchPnt* nc = NULL;
				for (int p = 0; p < (int) strokes[s]->pnts.size(); p++)
				{
					if (pnts[ strokes[s]->pnts[p] ]->coords[0] == result[2*c] && pnts[ strokes[s]->pnts[p] ]->coords[1] == result[2*c+1])
					{
						cout << "new corner detected: strokes[s]->pnts[" << p << "] = pnts[" << strokes[s]->pnts[p] << "]\n";
						nc = pnts[ strokes[s]->pnts[p] ];
						nc->endpnt = true; //info: pnts[ strokes[s]->pnts[p] ]->endpnt becomes true after nc update here
						endpnts.push_back(nc);
						break;
					}
				}
				//let nc split strokes[s] into 2 sub-strokes whose 1 endpnt is nc and the other is found by traversing along strokes[s]
				/////////////// time-backward substroke ////////////////////////
				Edge* substroke1 = new Edge((int) strokes.size(), false);
				int idx = nc->idx;
//if (c==11) cout << "ss\t" << idx << endl; 
				while (true) //populate pnts[] of this new substroke 'till an existing endpnts[] is reached; go backwards in time (idx--)
				{
					substroke1->pnts.push_back(idx);
					pnts[idx]->strokeID = (int) strokes.size();
					if ((int) substroke1->pnts.size() > 1)
						substroke1->length += distanceBetween2D(pnts[ substroke1->pnts[(int) substroke1->pnts.size() - 1] ]->coords, pnts[ substroke1->pnts[(int) substroke1->pnts.size() - 2] ]->coords);
					else if ((int) substroke1->pnts.size() == 1)
					{
						substroke1->v1i = idx; //first point (idx to pnts[])
						substroke1->ev1i = (int) endpnts.size() - 1; //v1i of this substroke1 as a pointer to endpnts[]
					}
					idx--;
					if (pnts[idx]->endpnt) //endpnt reached
					{
						//repeat stuff above for this last point and break
						substroke1->pnts.push_back(idx);
						pnts[idx]->strokeID = (int) strokes.size();
						substroke1->length += distanceBetween2D(pnts[ substroke1->pnts[(int) substroke1->pnts.size() - 1] ]->coords, pnts[ substroke1->pnts[(int) substroke1->pnts.size() - 2] ]->coords);
						break;
					}
				}
				substroke1->v2i = idx; //last point (idx to pnts[])
				substroke1->ev2i = getIdxToEndpnts(substroke1->v2i); //v2i of this substroke1 as a pointer to endpnts[]
				pnts[ substroke1->v1i ]->strokeList.push_back(substroke1->idx);
				pnts[ substroke1->v2i ]->strokeList.push_back(substroke1->idx);
				strokes.push_back(substroke1);
				/////////////// time-backward substroke ends ////////////////////////
				/////////////// time-forward substroke ////////////////////////
				Edge* substroke2 = new Edge((int) strokes.size(), false);
				idx = nc->idx;
				while (true) //populate pnts[] of this new substroke 'till an existing endpnts[] is reached; go forwards in time (idx++)
				{
					substroke2->pnts.push_back(idx);
					pnts[idx]->strokeID = (int) strokes.size();
					if ((int) substroke2->pnts.size() > 1)
						substroke2->length += distanceBetween2D(pnts[ substroke2->pnts[(int) substroke2->pnts.size() - 1] ]->coords, pnts[ substroke2->pnts[(int) substroke2->pnts.size() - 2] ]->coords);
					else if ((int) substroke2->pnts.size() == 1)
					{
						substroke2->v1i = idx; //first point (idx to pnts[])
						substroke2->ev1i = (int) endpnts.size() - 1; //v1i of this substroke2 as a pointer to endpnts[]
					}
					idx++;
					if (pnts[idx]->endpnt) //endpnt reached
					{
						//repeat stuff above for this last point and break
						substroke2->pnts.push_back(idx);
						pnts[idx]->strokeID = (int) strokes.size();
						substroke2->length += distanceBetween2D(pnts[ substroke2->pnts[(int) substroke2->pnts.size() - 1] ]->coords, pnts[ substroke2->pnts[(int) substroke2->pnts.size() - 2] ]->coords);
						break;
					}
				}
				substroke2->v2i = idx; //last point (idx to pnts[])
				substroke2->ev2i = getIdxToEndpnts(substroke2->v2i); //v2i of this substroke1 as a pointer to endpnts[]
				pnts[ substroke2->v1i ]->strokeList.push_back(substroke2->idx);
				pnts[ substroke2->v2i ]->strokeList.push_back(substroke2->idx);
				strokes.push_back(substroke2);
				/////////////// time-forward substroke ends ////////////////////////
				/////////////////////////////// creation of 2 substrokes ends /////////////////////////////
			}
		}
//cout << "hi there\n\n\n\n\n";
		delete [] result;
//exit(0);
		////////////////////// douglas-peucker line simplification algo, which acts as a corner-detector on my input stroke(s) ends //////////////////////
	}

	/////////////////////////////// deletion of splitted big stroke(s) /////////////////////////////
	//addition of new substrokes above makes splitted strokes extra and bad; so remove a stroke if it has 3+ endpnts on it
	for (int s = 0; s < (int) strokes.size(); s++)
	{
		int nEndpnts = 0;
		for (int p = 0; p < (int) strokes[s]->pnts.size(); p++)
			nEndpnts += (pnts[ strokes[s]->pnts[p] ]->endpnt ? 1 : 0);

/*if (strokes[s]->idx == 5)
{
	cout << "eeeeeeeeeeeeeeee " << nEndpnts << "\t" << strokes[s]->ev1i << "--" << strokes[s]->ev2i << "\n" << strokes[s]->v1i << "--" << strokes[s]->v2i << endl;
	for (int p = 0; p < (int) strokes[s]->pnts.size(); p++) cout << pnts[ strokes[s]->pnts[p] ]->idx << "x" << pnts[ strokes[s]->pnts[p] ]->endpnt<<"\t";cout<<endl;
	cout<<endl;
	for (int p = 0; p < (int) strokes[0]->pnts.size(); p++) cout << pnts[ strokes[0]->pnts[p] ]->idx << "y" << pnts[ strokes[0]->pnts[p] ]->endpnt<<"\t";cout<<endl;
}*/
		if (nEndpnts > 2)
		{
			strokes[s]->removed = true; //explicit deletion w/ erase() causes lots of index shifts; so keep strokes[s] physically but ignore throughout the program
cout << "splitted and to-be-ignored: strokes" << s << " w/ " << (int) strokes[s]->pnts.size() << " points and length = " << strokes[s]->length << "\n";
		}
	}
	/////////////////////////////// deletion of splitted big stroke(s) ends /////////////////////////////
	cout << "strokes after segmentation:\n";
	for (int s = 0; s < (int) strokes.size(); s++)
		if (! strokes[s]->removed)
			cout << "strokes" << strokes[s]->idx << ": endpnts" << strokes[s]->ev1i << "--" << strokes[s]->ev2i << endl;
}

void Mesh::dijkstraShortestPaths(int sourceVert, vector< int >* savePathToTargets)
{
	//computes shortest paths from sourceVert to all other verts

	if (verts[sourceVert]->dSpanningUnn)
		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);
		verts[v]->predecessorV = verts[v]->predecessorE = -1;
		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 surface/boundary edge (not the interior edges) 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
				verts[va]->predecessorV = minDidx;
				verts[va]->predecessorE = e;
				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 (savePathToTargets) //if (savePathToTargetV != -1) this was from the simpler 1-element savePathToTargetV, which is extended to multiple elements via savePathToTargets
		//in addition to desired sourceVert.dSpanningUnn[] values above, save the generator shortest/geodesic path from sourveVert to savePathToTargetV
		for (int i = 0; i < (int) (*savePathToTargets).size(); i++)
		{
			int savePathToTargetV = (*savePathToTargets)[i];
			bool edgeSet = true; //save shortest/geodesic path as a set of edges (true) or verts (false)
			int predecessor = (edgeSet ? verts[savePathToTargetV]->predecessorE : savePathToTargetV);
			vector< int > geoPath;
			while (predecessor != -1)
			{
				geoPath.push_back(predecessor);
				savePathToTargetV = verts[savePathToTargetV]->predecessorV;
//				if (edgeSet)
//					edges[predecessor]->constrainedByGeoPath = true;
				predecessor = (edgeSet ? verts[savePathToTargetV]->predecessorE : savePathToTargetV);				
			}
			verts[sourceVert]->geoPaths.push_back(geoPath); //i'th geoPath is set, i.e. verts[sourceVert]->geoPaths[i]
		}
		
	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	
}
void Mesh::dijkstraShortestPathsSketchE(int sourcePnt)
{
	//same as dijkstraShortestPaths() except sketch calls here and therefore use SketchPnts and strokeList instead of the usual verts[] and the edgeList for meshes (dijkstraShortestPaths() case)

	if (endpnts[sourcePnt]->dSpanningUnn)
		return;
	endpnts[sourcePnt]->dSpanning = new float[ (int) endpnts.size() ]; //in [0,1] scale (normalized)
	endpnts[sourcePnt]->dSpanningUnn = new float[ (int) endpnts.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) endpnts.size(); v++)
	{
		endpnts[v]->dSpanned = (v == sourcePnt ? 0.0f : INF);
		FibHeapNode* hn = new FibHeapNode(endpnts[v]->dSpanned, v); //v is a handle from headNode to endpnts,
		endpnts[v]->heapNode = hn;				  				    //and hn is an handle from endpnts 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 endpnts w/ unknown shortest dist; so, we're done
		int minDidx = minNode->element;
		float minD = minNode->key;
		endpnts[sourcePnt]->dSpanningUnn[minDidx] = minD; //remember original value before normalization in the end of this function
//if (sourcePnt==8) {cout << sourcePnt << "\t\t" << minDidx << "\t\t" << minD << endl;}
		delete minNode; //new hn above is deleted here to save memo [info: delete minNode; automatically calls minNode->~FibHeapNode(), i.e. its destructor]

		//relax each surface/boundary stroke incident to minDidx
		for (int ve = 0; ve < (int) endpnts[minDidx]->strokeList.size(); ve++)
		{
			int e = endpnts[minDidx]->strokeList[ve];
			if (strokes[e]->removed)
				continue; //ignore removed strokes (they're kept physically for index consistency)
			int va = (strokes[e]->ev1i == minDidx ? strokes[e]->ev2i : strokes[e]->ev1i);
			if (endpnts[minDidx]->dSpanned + strokes[e]->length < endpnts[va]->dSpanned)
			{
				endpnts[va]->dSpanned = endpnts[minDidx]->dSpanned + strokes[e]->length; //relaxation
				h->DecreaseKey(endpnts[va]->heapNode, endpnts[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) endpnts.size(); i++)
		endpnts[sourcePnt]->dSpanning[ i ] = (endpnts[sourcePnt]->dSpanningUnn[i] - minGeoDist) / denom; //don't forget to change minGeoDist to minGeoDistG if you change denom	
}
void Mesh::dijkstraShortestPathsSketch(int sourcePnt)
{
	//same as dijkstraShortestPathsSketchE() except this one uses pnts[] and edgeList instead of endpnts[]

	if (pnts[sourcePnt]->dSpanningUnn)
		return;
	pnts[sourcePnt]->dSpanning = new float[ (int) pnts.size() ]; //in [0,1] scale (normalized)
	pnts[sourcePnt]->dSpanningUnn = new float[ (int) pnts.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) pnts.size(); v++)
	{
		pnts[v]->dSpanned = (v == sourcePnt ? 0.0f : INF);
		FibHeapNode* hn = new FibHeapNode(pnts[v]->dSpanned, v); //v is a handle from headNode to pnts,
		pnts[v]->heapNode = hn;				  				     //and hn is an handle from pnts 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 pnts w/ unknown shortest dist; so, we're done
		int minDidx = minNode->element;
		float minD = minNode->key;
		pnts[sourcePnt]->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]

		//relax each surface/boundary stroke incident to minDidx
		for (int ve = 0; ve < (int) pnts[minDidx]->edgeList.size(); ve++)
		{
			int e = pnts[minDidx]->edgeList[ve];
			//if (edges[e]->removed) continue; //ignore removed edges (they're kept physically for index consistency) [edges never removed]
			int va = (edges[e]->v1i == minDidx ? edges[e]->v2i : edges[e]->v1i);
			if (pnts[minDidx]->dSpanned + edges[e]->length < pnts[va]->dSpanned)
			{
				pnts[va]->dSpanned = pnts[minDidx]->dSpanned + edges[e]->length; //relaxation
				h->DecreaseKey(pnts[va]->heapNode, pnts[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) pnts.size(); i++)
		pnts[sourcePnt]->dSpanning[ i ] = (pnts[sourcePnt]->dSpanningUnn[i] - minGeoDist) / denom; //don't forget to change minGeoDist to minGeoDistG if you change denom	
}
#ifdef USE_GOOD_CONTINUATION
void Mesh::dijkstraShortestPathsSketchGC(int sourcePnt)
{
	//same as dijkstraShortestPathsSketch() except this one minimizes deviation from good continuation instead of the distances (standard dijkstraShortestPathsSketch() case)

	if (pnts[sourcePnt]->dSpanningUnn)
		return;
	pnts[sourcePnt]->dSpanning = new float[ (int) pnts.size() ]; //in [0,1] scale (normalized)
	pnts[sourcePnt]->dSpanningUnn = new float[ (int) pnts.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) pnts.size(); v++)
	{
		pnts[v]->dSpanned = (v == sourcePnt ? 0.0f : INF);
		pnts[v]->predecessorP = (v == sourcePnt ? v : -1); //for sourcePnt, predecessorP is -1 naturally; to prevent crashes below, consider this sourcePnt as the predecessorP of itself
//		pnts[v]->processed = false; //true means pnts[sourcePnt]->dSpanningUnn[v] is finalized and will never change (redundant for standard dijkstra but here mandatory 'cos i'm doing a different relaxation) [not in use]
		FibHeapNode* hn = new FibHeapNode(pnts[v]->dSpanned, v); //v is a handle from headNode to pnts,
		pnts[v]->heapNode = hn;				  				     //and hn is an handle from pnts to heapNode
		h->Insert(hn);
		//no delete hn; here 'cos it'll be extracted, used, and deleted below
	}

//float maxL = -INF, minL = INF;int p1min, p1max, p2min, p2max;
	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 pnts w/ unknown shortest dist; so, we're done
		int minDidx = minNode->element;
		float minD = minNode->key;
		pnts[sourcePnt]->dSpanningUnn[minDidx] = minD; //remember original value before normalization in the end of this function
//		pnts[minDidx]->processed = true; //[sourcePnt]->dSpanningUnn[minDidx] finalized just-above
		delete minNode; //new hn above is deleted here to save memo [info: delete minNode; automatically calls minNode->~FibHeapNode(), i.e. its destructor]

		//decide the edge to be used in relaxation below [no decision; directly relax each surface/boundary edge incident to minDidx]
		float minLinenessRatio = INF;
		int minE = -1; //edge giving the minLinenessRatio
//cout<<"queue" << minDidx << " (" << minD << ") preceded by " << pnts[minDidx]->predecessorP << "\t" << (int) pnts[minDidx]->edgeList.size() << endl;
		for (int ve = 0; ve < (int) pnts[minDidx]->edgeList.size(); ve++)
		{
//if (minDidx == 379) cout << ve << "ve\n";
			int e = pnts[minDidx]->edgeList[ve];
//if (minDidx == 379) cout << e << "eeeeeeeeeeeeeeeee\n";
			//if (edges[e]->removed) continue; //ignore removed edges (they're kept physically for index consistency) [edges never removed]
			int va = (edges[e]->v1i == minDidx ? edges[e]->v2i : edges[e]->v1i);
//if (minDidx == 379) cout << va << " " << pnts[va]->processed << " " << pnts[minDidx]->predecessorP << " " << pnts[minDidx]->dSpanned << " vaaaaaaa\n";
//			if (pnts[va]->processed) //va is a bad candidate 'cos its dSpanningUnn to sourcePnt is already decided
//				continue;
			float pathLen = distanceBetween2D(pnts[ pnts[minDidx]->predecessorP ]->coords, pnts[minDidx]->coords) + edges[e]->length, //geodesic length consisting of 2 distances: from predecessor to minDidx + from minDidx to va
				  eucLen = distanceBetween2D(pnts[ pnts[minDidx]->predecessorP ]->coords, pnts[va]->coords), //euclidean distance between predecessor and va
				  linenessRatio = (eucLen != 0.0f ? pathLen / eucLen : 0.0f); //pathLen / eucLen ratio is >= 1 for sure
//cout<<linenessRatio<<"\t";
/*if (edges[e]->v1i == 350 && edges[e]->v2i == 362) cout << "?????? " << pathLen << " / " << eucLen << "::::::: " << pnts[minDidx]->predecessorP << " " << minDidx << " " << va << endl;
if (linenessRatio > maxL) {maxL = linenessRatio;p1max=edges[e]->v1i;p2max=edges[e]->v2i;}
if (linenessRatio < minL) {minL = linenessRatio;p1min=edges[e]->v1i;p2min=edges[e]->v2i;}*/
			if (pnts[minDidx]->dSpanned + linenessRatio < pnts[va]->dSpanned)
			{
				pnts[va]->dSpanned = pnts[minDidx]->dSpanned + linenessRatio; //relaxation
				pnts[va]->predecessorP = minDidx;
				h->DecreaseKey(pnts[va]->heapNode, pnts[va]->dSpanned); //worst O(lgV), amortized O(1)
			}
//cout << linenessRatio << " --------- " << minLinenessRatio << "\n";
/*			if (linenessRatio < minLinenessRatio)
			{
				minLinenessRatio = linenessRatio;
				minE = e;
			}*/
		}
//cout<<minE << "a\n";
/*		if (minE != -1) //prevent array out of bounds error; true means minDidx not have any unfinalized neighbors so no relaxation at all
		{
			//relaxation: relax minE edge incident to minDidx
			int va = (edges[minE]->v1i == minDidx ? edges[minE]->v2i : edges[minE]->v1i);
//cout << minDidx << "uuuuu " << pnts[minDidx]->dSpanned + edges[minE]->length << " vs. " << pnts[va]->dSpanned << "\n\n";
			if (pnts[minDidx]->dSpanned + edges[minE]->length < pnts[va]->dSpanned)
			{
				pnts[va]->dSpanned = pnts[minDidx]->dSpanned + edges[minE]->length; //relaxation
				pnts[va]->predecessorP = minDidx;
				h->DecreaseKey(pnts[va]->heapNode, pnts[va]->dSpanned); //worst O(lgV), amortized O(1)
			}
		}*/
	} //end of while (true)
	/////// original Dijkstra to compute shortest paths from sourceVert to all others ///////
//exit(0);
//cout << sourcePnt << "s\t" << minL << " b/w " << p1min << " & " << p2min << endl;cout << sourcePnt << "s\t" << maxL << " b/w " << p1max << " & " << p2max << endl << endl;

	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) pnts.size(); i++)
		pnts[sourcePnt]->dSpanning[ i ] = (pnts[sourcePnt]->dSpanningUnn[i] - minGeoDist) / denom; //don't forget to change minGeoDist to minGeoDistG if you change denom

	if (sourcePnt == 362)
	{
		int b = 553;
		cout << "geodesic path from pnts" << sourcePnt << " to pnts" << b << ":\n";
		int predecessor = b;
		while (predecessor != -1) //never false in this project
		{
			cout << predecessor << "\t";
			b = pnts[b]->predecessorP;
			predecessor = b;
			//below is necessary 'cos i don't use predecessor=-1 for source pnts for technical reasons
			if (predecessor == sourcePnt)
				break;
		}
		cout << predecessor << "\n\n";
	}
}
#endif

void Mesh::computeEvenlySpacedSamples()
{
	//sampling of N evenly-spaced FPS samples (that naturally/automatically land on extremes)

	if (! samples.empty())
		return; //already computed so do nothing

	cout << "Evenly-spaced sampling of " << N << " extremities [mesh]......\n";

	if (N > (int) verts.size() || N <= 0)
	{
		cout << "# samples requested > # verts in mesh or nonpositive request\n";
		exit(0);
	}

	minGeoDist = INF;
	maxGeoDist = -INF;
	//fixed number (N) of fps (farthest point sampling) samples goes to samples[] first of which is an arbitrary vertex
	int extVert = getAnExtremeVert();
	for (int v = 0; v < (int) verts.size(); v++)
		verts[v]->sample = verts[v]->processed = false; //initialize this to prevent duplications in samples; true means added to samples	
	samples.push_back(extVert);
	dijkstraShortestPaths(extVert); //dSpanningUnn[] to be used below is now ready, dSpanning is bad due to unhealthy min/maxGeoDist values (to be corrected at the bottom)
	verts[extVert]->sample = verts[extVert]->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)
	//as a result, 26 secs for N=500 on a V=12.5K scape model, just 1- for N=9
	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 select s as a sample 'cos it's already added or is an interior vert

			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]->sample = 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

//if (frameID < 6 && (samples[7] > 2141 && samples[7] < 2250)) {samples[7] = 6575;dijkstraShortestPaths(samples[7]);} do this for centaur only to land a sample on head instead of on back

	//set min/maxGeoDist properly which in turn gives me a chance to set the dSpanning[] values properly
	minGeoDist = INF;
	maxGeoDist = -INF;
	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);}
		}
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 << " .. ";
	}
#ifdef PRINT_DETAILS
	cout << "\ngeodesics:\n";
	for (int i = 0; i < (int) samples.size(); i++)
		for (int j = i+1; j < (int) samples.size(); j++)
			cout << verts[ samples[i] ]->idx << "-" << verts[ samples[j] ]->idx << "\t" << verts[ samples[i] ]->dSpanning[ samples[j] ] << "\t" << verts[ samples[i] ]->dSpanningUnn[ samples[j] ] << endl;
	cout << "\nmaxGeoDist on mesh = " << maxGeoDist << "\n\n";
#endif
}
bool Mesh::computeEvenlySpacedSamplesSketch()
{
	//same as computeEvenlySpacedSamples() except this one works on pnts[] of this sketch (or on endpnts[] of sketch strokes)

	//samples[] already filled by computeEvenlySpacedSamplesSketchEuc() w/ dense uniform pnts; start from scratch
	for (int v = 0; v < (int) samples.size(); v++)
		pnts[ samples[v] ]->sample = false;
	samples.clear();
	
	usePnts = false; //true: uses pnts[]; false: uses endpnts[] throughout the extremity sampling and matching
	if (usePnts)
	{
		//after setStrokes(), each stroke is represented by its endpnts[]; so use those endpnts[] from now on
		cout << "must use endpnts[], not pnts[]\n";
		exit(0);
	}

	cout << "Evenly-spaced sampling of " << N << " extremities [sketch]......\n";
	if (N <= 0 || (usePnts && N > (int) pnts.size()) || (! usePnts && N > (int) endpnts.size()))
	{
		cout << "# samples requested > # stroke endpnts[] (" << (int) endpnts.size() << ") or sketch pnts[] in sketch or nonpositive request\n";
		exit(0);
	}
	
	minGeoDist = INF;
	maxGeoDist = -INF;
	//fixed number (N) of fps (farthest point sampling) samples goes to samples[] first of which is an arbitrary vertex
	int extVert = getAnExtremeVertSketch();
	samples.push_back(extVert);
	if (usePnts) //pnts[]
	{
		for (int v = 0; v < (int) pnts.size(); v++)
			pnts[v]->sample = pnts[v]->processed = false; //initialize this to prevent duplications in samples; true means added to samples
#ifdef USE_GOOD_CONTINUATION
		dijkstraShortestPathsSketchGC(extVert); //dSpanningUnn[] to be used below is now ready, dSpanning is bad due to unhealthy min/maxGeoDist values (to be corrected at the bottom);
#else
		dijkstraShortestPathsSketch(extVert); //dSpanningUnn[] to be used below is now ready, dSpanning is bad due to unhealthy min/maxGeoDist values (to be corrected at the bottom)
#endif
		pnts[extVert]->sample = pnts[extVert]->processed = true; //to prevent duplications in samples; true means added to samples
	}
	else //endpnts[]
	{
		for (int v = 0; v < (int) endpnts.size(); v++)
			endpnts[v]->sample = endpnts[v]->processed = false; //initialize this to prevent duplications in samples; true means added to samples
		dijkstraShortestPathsSketchE(extVert); //dSpanningUnn[] to be used below is now ready, dSpanning is bad due to unhealthy min/maxGeoDist values (to be corrected at the bottom)
		endpnts[extVert]->sample = endpnts[extVert]->processed = true; //to prevent duplications in samples; true means added to samples	
	}
	
//cout << endpnts[extVert]->idx << "\n\n\n\n\n\n";return false;

	//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)
	//as a result, 26 secs for N=500 on a V=12.5K scape model, just 1- for N=9
	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;
		if (usePnts) //pnts[]
		{
			for (int s = 0; s < (int) pnts.size(); s++) //fps selects from amongst all vertices
			{
				int e2 = -1;
				if (pnts[s]->processed)
					continue; //no chance to select s as a sample 'cos it's already added or is an interior vert

				float minDist = pnts[ 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 (pnts[ samples[e] ]->dSpanningUnn[s] <= minDist)
					{
						minDist = pnts[ samples[e] ]->dSpanningUnn[s];
						e2 = samples[e];
					}
				if (pnts[e2]->dSpanningUnn[s] > maxDist)
				{
					maxDist = pnts[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);
			pnts[nextSample]->sample = pnts[nextSample]->processed = true; //added, so not add again later
#ifdef USE_GOOD_CONTINUATION
			dijkstraShortestPathsSketchGC(nextSample); //upcoming iteration needs samples.dSpanningUnn which is not currently ready when selecting from the pool of all verts;
#else
			dijkstraShortestPathsSketch(nextSample); //upcoming iteration needs samples.dSpanningUnn which is not currently ready when selecting from the pool of all verts
#endif
		}
		else //endpnts[]
		{
			for (int s = 0; s < (int) endpnts.size(); s++) //fps selects from amongst all vertices
			{
				int e2 = -1;
				if (endpnts[s]->processed)
					continue; //no chance to select s as a sample 'cos it's already added or is an interior vert

				float minDist = endpnts[ 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 (endpnts[ samples[e] ]->dSpanningUnn[s] <= minDist)
					{
						minDist = endpnts[ samples[e] ]->dSpanningUnn[s];
						e2 = samples[e];
					}
				if (endpnts[e2]->dSpanningUnn[s] > maxDist)
				{
					maxDist = endpnts[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);
			endpnts[nextSample]->sample = endpnts[nextSample]->processed = true; //added, so not add again later		
			dijkstraShortestPathsSketchE(nextSample); //upcoming iteration needs samples.dSpanningUnn which is not currently ready when selecting from the pool of all verts
		}
	} //end of while N
//if (target) {samples[4] = 493;dijkstraShortestPaths(samples[4]);}

	//set min/maxGeoDist properly which in turn gives me a chance to set the dSpanning[] values properly
	minGeoDist = INF;
	maxGeoDist = -INF;
	bool disconnection = false; //becomes true if there are disconnected extremities, which triggers me to re-do the setStrokeConnectivity() w/ a larger `close' threshold
	if (usePnts) //pnts[]
	{
		for (int i = 0; i < (int) samples.size(); i++)
			for (int v = 0; v < (int) pnts.size(); v++)
			{
				if (pnts[ samples[i] ]->dSpanningUnn[ v ] < minGeoDist)
					minGeoDist = pnts[ samples[i] ]->dSpanningUnn[ v ];
				if (pnts[ samples[i] ]->dSpanningUnn[ v ] > maxGeoDist)
					maxGeoDist = pnts[ samples[i] ]->dSpanningUnn[ v ];
				if (pnts[ samples[i] ]->dSpanningUnn[ v ] == INF)
				{
					cout<<"pnts" << samples[i] << " " << v << " not connected by a path; adding more virtual strokes by auto-increasing `close' variable\n\n";
					//cout << "euc = " << distanceBetween2D(pnts[ samples[i] ]->coords, pnts[v]->coords) << endl;
					disconnection = true;//exit(0);
				}
			}
		if (disconnection) //disconnected sample(s) should be canceled
		{
			//clean stuff as this function'll be called again
			for (int i = 0; i < (int) samples.size(); i++)
			{
				if (pnts[ samples[i] ]->dSpanning) //may already be deleted in for-loop above
					delete pnts[ samples[i] ]->dSpanning;
				if (pnts[ samples[i] ]->dSpanningUnn)
					delete pnts[ samples[i] ]->dSpanningUnn; //arrays allocated with new[] must be deallocated with delete[] so i should've added [] here but not a big deal
				pnts[ samples[i] ]->dSpanning = pnts[ samples[i] ]->dSpanningUnn = NULL;
			}
			return true; //disconnection detected
		}

		//clean old samples.dSpanning values
		for (int i = 0; i < (int) samples.size(); i++)
		{
			if (pnts[ samples[i] ]->dSpanning) //may already be deleted in for-loop above
				delete pnts[ samples[i] ]->dSpanning;
			if (pnts[ samples[i] ]->dSpanningUnn)
				delete pnts[ samples[i] ]->dSpanningUnn; //arrays allocated with new[] must be deallocated with delete[] so i should've added [] here but not a big deal
			pnts[ samples[i] ]->dSpanning = pnts[ samples[i] ]->dSpanningUnn = NULL;

			//dSpanning[] good now 'cos min/maxGeoDist ready; dSpanning[]'ll be used during distortion computations of maps b/w shape samples
#ifdef USE_GOOD_CONTINUATION
			dijkstraShortestPathsSketchGC(samples[i]);
#else
			dijkstraShortestPathsSketch(samples[i]);
#endif	
			if (i < 5 || i >= N-5) //print the first and last 5 sample ids
				cout << samples[i] << " ";
			else if (i == 5)
				cout << " .. ";
		}
#ifdef PRINT_DETAILS
		cout << "\ngeodesics:\n";
		for (int i = 0; i < (int) samples.size(); i++)
			for (int j = i+1; j < (int) samples.size(); j++)
				cout << pnts[ samples[i] ]->idx << "-" << pnts[ samples[j] ]->idx << "\t" << pnts[ samples[i] ]->dSpanning[ samples[j] ] << "\t" << pnts[ samples[i] ]->dSpanningUnn[ samples[j] ] << endl;
#endif
	}
	else //endpnts[]
	{
		for (int i = 0; i < (int) samples.size(); i++)
			for (int v = 0; v < (int) endpnts.size(); v++)
			{
				if (endpnts[ samples[i] ]->dSpanningUnn[ v ] < minGeoDist)
					minGeoDist = endpnts[ samples[i] ]->dSpanningUnn[ v ];
				if (endpnts[ samples[i] ]->dSpanningUnn[ v ] > maxGeoDist)
					maxGeoDist = endpnts[ samples[i] ]->dSpanningUnn[ v ];
				if (endpnts[ samples[i] ]->dSpanningUnn[ v ] == INF)
				{
					cout<<"endpnts" << samples[i]<< " " << v << " (= pnts " << endpnts[ samples[i] ]->idx << " " << endpnts[v]->idx << ") not connected by a path; adding more virtual strokes by auto-increasing `close' variable\n\n";
					//cout << "euc = " << distanceBetween2D(endpnts[ samples[i] ]->coords, pnts[v]->coords) << endl;
					disconnection = true;//exit(0);
					i = (int) samples.size();
					break;
				}
			}
		if (disconnection) //disconnected sample(s) should be canceled
		{
			//clean stuff as this function'll be called again
			for (int i = 0; i < (int) samples.size(); i++)
			{
				if (endpnts[ samples[i] ]->dSpanning) //may already be deleted in for-loop above
					delete endpnts[ samples[i] ]->dSpanning;
				if (endpnts[ samples[i] ]->dSpanningUnn)
					delete endpnts[ samples[i] ]->dSpanningUnn; //arrays allocated with new[] must be deallocated with delete[] so i should've added [] here but not a big deal
				endpnts[ samples[i] ]->dSpanning = endpnts[ samples[i] ]->dSpanningUnn = NULL;
			}
			return true; //disconnection detected
		}

		//clean old samples.dSpanning values
		for (int i = 0; i < (int) samples.size(); i++)
		{
			if (endpnts[ samples[i] ]->dSpanning) //may already be deleted in for-loop above
				delete endpnts[ samples[i] ]->dSpanning;
			if (endpnts[ samples[i] ]->dSpanningUnn)
				delete endpnts[ samples[i] ]->dSpanningUnn; //arrays allocated with new[] must be deallocated with delete[] so i should've added [] here but not a big deal
			endpnts[ samples[i] ]->dSpanning = endpnts[ samples[i] ]->dSpanningUnn = NULL;

			//dSpanning[] good now 'cos min/maxGeoDist ready; dSpanning[]'ll be used during distortion computations of maps b/w shape samples
			dijkstraShortestPathsSketchE(samples[i]);
		
			if (i < 5 || i >= N-5) //print the first and last 5 sample ids
				cout << endpnts[ samples[i] ]->idx << " ";
			else if (i == 5)
				cout << " .. ";
		}
#ifdef PRINT_DETAILS
		cout << "\ngeodesics:\n";
		for (int i = 0; i < (int) samples.size(); i++)
			for (int j = i+1; j < (int) samples.size(); j++)
				cout << endpnts[ samples[i] ]->idx << "-" << endpnts[ samples[j] ]->idx << "\t" << endpnts[ samples[i] ]->dSpanning[ samples[j] ] << "\t" << endpnts[ samples[i] ]->dSpanningUnn[ samples[j] ] << endl;
#endif
	}
#ifdef PRINT_DETAILS
	cout << "\nmaxGeoDist on sketch = " << maxGeoDist << "\n\n";
#else
	cout << "\n\n";
#endif
	return false;
}

void Mesh::preprocessSketch()
{
	//uniform resampling of dense sketch pnts, connectivity of those samples, strokes as clustered samples, connectivity of those strokes

	computeEvenlySpacedSamplesSketchEuc();
	setSampleConnetivity();
	setStrokes();
	setStrokeConnectivity();
//cout << distanceBetween2D(pnts[134]->coords, pnts[87]->coords) << " euc\n\n\n\n";
}
void Mesh::computeEvenlySpacedSamplesSketchEuc()
{
	//same as computeEvenlySpacedSamplesSketch() except this one is based on euclidean distances, and not the dijkstra geodesic distances

	if (! samples.empty())
		return; //already computed so do nothing

	cout << "\nEvenly-spaced (" << spacing << ") dense sampling based on Euclidean distances......\n";

	for (int v = 0; v < (int) pnts.size(); v++)
		pnts[v]->sample = pnts[v]->processed = false; //sample=true means added to samples (prevent duplications in samples) and processed=true means marked as unavailable for being the next sample
	int samplePnt = 0;
	do
	{
		samples.push_back(samplePnt);
		pnts[samplePnt]->sample = true;
		for (int v = 0; v < (int) pnts.size(); v++)
			if (distanceBetween2D(pnts[v]->coords, pnts[samplePnt]->coords) <= spacing)
				pnts[v]->processed = true; //euclidean-neighborhood of new sample samplePnt is marked as unavailable
		//next sample pnt
		samplePnt = -1;
		for (int v = 0; v < (int) pnts.size(); v++)
			if (! pnts[v]->processed)
				samplePnt = v;
	} while (samplePnt != -1);
	cout << (int) samples.size() << " samples generated\n";
}
void Mesh::setSampleConnetivity()
{
	//sets the neighborhood info of the samples of this sketch; it makes few neighbors for each sample (1 if it is on the tip) which are the closest pnts; essentially 
	//create an edge b/w sample s and its closest neighbor which is not connected to s yet (first closest might have already connected in a prev iter, which causes 2+ neighbors for s)

	float dspacing = 2*spacing;
	for (int s = 0; s < (int) samples.size(); s++)	   //set neighbor(s) of s
	{
		float minDist = INF;
		int minT = -1;
		for (int t = 0; t < (int) samples.size(); t++) //by checking all other samples t
		{
			if (s == t) //skip t=s case as a sample cannot be the neighbor of itself
				continue;
			float dist = distanceBetween2D(pnts[ samples[s] ]->coords, pnts[ samples[t] ]->coords);
			if (dist < minDist) //closestPntNghb=true which means i am selecting the closest && not-yet-neighbor sample as the neighbor
			{
				bool alreadyNghb = false; //not-yet-neighbor test
				for (int e = 0; e < (int) pnts[ samples[s] ]->edgeList.size(); e++) //check for the existence of t in the incident edges of s (same as checking for the existence of s in the incident edges of t, due to duality)
				{
					int p = (edges[ pnts[ samples[s] ]->edgeList[e] ]->v1i == samples[s] ? edges[ pnts[ samples[s] ]->edgeList[e] ]->v2i : edges[ pnts[ samples[s] ]->edgeList[e] ]->v1i);
					if (p == samples[t])
					{
						alreadyNghb = true;
						break;
					}
				}
				if (! alreadyNghb)
				{
					minDist = dist;
					minT = samples[t];
				}
			}
		}
		if (minDist <= dspacing) //s and minT should be close enough to prevent weird pairings from one tip to another
		{
			//create a new edge b/w s and minT
			if (pnts[ samples[s] ]->addPntNeighbor(minT)) //addVertNeighbor returns true if s is not neighbor of minT; so they can form a new edge
				addEdge(samples[s], minT, minDist);
			pnts[minT]->addPntNeighbor(samples[s]);	//s minT edge is (probably) added just-above; do not add minT s edge!!!
		}
	}
#ifdef PRINT_DETAILS
//	for (int s = 0; s < (int) samples.size(); s++) 
//		cout << samples[s] << " has " << (int) pnts[ samples[s] ]->edgeList.size() << " neighboring samples\n";
//if (samples[s] == 0) cout << edges[ pnts[ samples[s] ]->edgeList[0] ]->v2i << ", " << edges[ pnts[ samples[s] ]->edgeList[1] ]->v1i << ", " << edges[ pnts[ samples[s] ]->edgeList[2] ]->v1i << "\n\n";
#endif
	avgEdgeLen = edgeLenTotal / (int) edges.size();
#ifdef PRINT_DETAILS
	cout << (int) edges.size() << " edges on this sketch w/ an average length of " << avgEdgeLen << "\n\n";
#endif
}
void Mesh::setStrokes()
{
	//sets the strokes of this sketch; a stroke is a set of ordered points (i.e. polyline) w/ similar slopes, where slope of pnt i is computed using the line from i-1 to i+1

	bool first = true, slopeBasedContinuation = false, //true makes slope-based continuation in strokes, false makes geo/euc-ratio-based continuation [better]
		 doDougPeuckerPostprocessing = true; //true launches doug-peucker corner detection on strokes set in do-while loop below in order to segment non-line strokes further using a more global/robust method; false skips doug-peucker stuff
	for (int s = 0; s < (int) samples.size(); s++)
		pnts[ samples[s] ]->processed = false; //true means this is a member of a stroke
	int startingDegree = 1; //vertices of degree startingDegree will be used to start a new stroke
	do //decide next startingDegree value
	{
		while (true)//(int) strokes.size() < 7) //strokes starting from startingDegree-degree vertices
		{
			int p = -1;
			//start from a new/unprocessed tip (or part) of the object, which is detected as the pnt w/ only 1 neighbor (w/ startingDegree neighbors)
			for (int s = 0; s < (int) samples.size(); s++)
				if (! pnts[ samples[s] ]->processed && (int) pnts[ samples[s] ]->edgeList.size() == startingDegree) //e.g. handle all strokes starting w/ a 1-degree vertex when startingDegree=1
				{
					p = samples[s];
					break;
				}
			if (p == -1)
				break; //no more starting point, i.e. tip, is left so quit the while-loop

			//generate the stroke starting at p
			//stroke is a polyline consisting of idxs to pnts[]; stroke.pnts[0]=v1i and stroke.pnts[last]=v2i can be used to find its length; or i can accumulate small steps from v1i to v2i for a more accurate length [in use]
			Edge* stroke = new Edge((int) strokes.size(), false);
			float tiny = 0.000001f, denom, maxSlopeDiffAllowed = 1.0f, maxRatioAllowed = 1.05f;//(sketchID == 7 ? 1.021f : 1.05f); //extension is a good continuation if it generates a slope-difference or ratio below these thresholds, respectively
																						 //1.05f is based on shortstraw paper; enable cout<<ratio below to come up w/ the best maxRatioAllowed value; or even better use doug-peucker corner detection to segment strokes further
			while (p != -1)
			{
#ifdef PRINT_DETAILS
				cout << "pnts[" << p << "] as the next pnt of strokes[" << (int) strokes.size() << "]\n";
#endif
				stroke->pnts.push_back(p);
				int n = (int) stroke->pnts.size(); //pnts[p] and pnts[ stroke.pnts[n-1] ] are the same, i.e. p = stroke.pnts[n-1]
				pnts[p]->strokeID = (int) strokes.size();
				pnts[p]->processed = true; //added to this stroke
				if (n > 1)
				{
					stroke->length += distanceBetween2D(pnts[ stroke->pnts[n-1] ]->coords, pnts[ stroke->pnts[n-2] ]->coords);
					if (n == 2)
					{
						if (slopeBasedContinuation)
						{
							denom = (pnts[ stroke->pnts[n-1] ]->coords[0] - pnts[ stroke->pnts[n-2] ]->coords[0] != 0.0f ? pnts[ stroke->pnts[n-1] ]->coords[0] - pnts[ stroke->pnts[n-2] ]->coords[0] : tiny);
if (denom==tiny){cout<<"0000 case; may be good as is\n";exit(0);}
							//newly added point of the next iteration (nextP below) must respect this slope value, which is approximated from the current local part of this stroke
							stroke->slope = (pnts[ stroke->pnts[n-1] ]->coords[1] - pnts[ stroke->pnts[n-2] ]->coords[1]) / denom; //slope = y2-y1 / x2-x1, that is rise/run, change in ys / change in xs
						}
					}
				}
				else if ((int) stroke->pnts.size() == 1)
					stroke->v1i = p; //first point (idx to pnts[])

				//select the next p
				int nextP = -1;
				for (int e = 0; e < (int) pnts[p]->edgeList.size(); e++) //pnts[p] = pnts[ stroke->pnts[n-1] ], the lastly added point
				{
					int pa = (edges[ pnts[p]->edgeList[e] ]->v1i == p ? edges[ pnts[p]->edgeList[e] ]->v2i : edges[ pnts[p]->edgeList[e] ]->v1i);
					if (pnts[pa]->processed)
						continue; //pa is already a member of a stroke so skip it
					if ((int) stroke->pnts.size() > 1)
					{
						if (slopeBasedContinuation)
						{
							//slope test using 3 pnts: pnts[p].prev (n-2), pnts[p] (n-1), pnts[pa]; actually i use the slope of the line from [p].prev to [pa]
							float denom = (pnts[pa]->coords[0] - pnts[ stroke->pnts[n-2] ]->coords[0] != 0.0f ? pnts[pa]->coords[0] - pnts[ stroke->pnts[n-2] ]->coords[0] : tiny), slope = (pnts[pa]->coords[1] - pnts[ stroke->pnts[n-2] ]->coords[1]) / denom;
if (denom==tiny){cout<<"0 case; may be good as is\n";exit(0);}
//cout << slope << " vs " << stroke->slope << ":: " << fabs(slope - stroke->slope) << "   " << "\t\t" << stroke->pnts[n-2] << " " << p << " " << pa << "\n\n";
							if (fabs(slope - stroke->slope) < maxSlopeDiffAllowed)
							{
								nextP = pa;
								stroke->slope = slope; //newly added point of the next iteration must respect this slope value, which is approximated from the current local part of this stroke
								break;
							}
						}
						else //geo/euc ratio based continuation
						{
							//ratio test using 3 pnts: pnts[p].prev (n-2), pnts[p] (n-1), pnts[pa]
							float geoLen = distanceBetween2D(pnts[ stroke->pnts[n-2] ]->coords, pnts[ stroke->pnts[n-1] ]->coords) + distanceBetween2D(pnts[ stroke->pnts[n-1] ]->coords, pnts[pa]->coords), //geodesic length consisting of 2 distances: from p.prev to p + from p to pa
								  eucLen = distanceBetween2D(pnts[ stroke->pnts[n-2] ]->coords, pnts[pa]->coords), //euclidean distance between p.prev and pa
								  ratio = (eucLen != 0.0f ? geoLen / eucLen : 0.0f); //geoLen / eucLen ratio is >= 1 for sure
//cout << ratio << "   " << "\t\t" << stroke->pnts[n-2] << " " << p << " " << pa << "\n\n";
							if (ratio < maxRatioAllowed)
							{
								nextP = pa;
								break;
							}
						}
					}
					else
					{
						//no slope test can be done 'cos i need at least 2 pnts in this stroke to perform a slope test; so add nextP directly
						nextP = pa;
						break;
					}
				}
				p = nextP;
			}
			if ((int) stroke->pnts.size() == 1) //rare case where this stroke has only 1 pnt; this crashes shortest path algo so add a new/2nd pnt w/ the same coords
			{
				pnts.push_back(new SketchPnt("virtualPnt", pnts[ stroke->v1i ]->coords, (int) pnts.size()));
				int newp = (int) pnts.size() - 1;
#ifdef PRINT_DETAILS
				cout << "pnts[" << newp << "] as the next (virtual) pnt of strokes[" << (int) strokes.size() << "]\n";
#endif
				stroke->pnts.push_back(newp);
				pnts[newp]->strokeID = (int) strokes.size();
				pnts[newp]->processed = true; //added to this stroke
			}
			//stroke finishes w/ this last pnt
			stroke->v2i = stroke->pnts[ (int) stroke->pnts.size()-1 ]; //last point (idx to pnts[])
			pnts[ stroke->v1i ]->strokeList.push_back(stroke->idx);
			pnts[ stroke->v2i ]->strokeList.push_back(stroke->idx);
			strokes.push_back(stroke);
#ifdef PRINT_DETAILS
			cout << "strokes[" << (int) strokes.size()-1 << "] w/ " << (int) stroke->pnts.size() << " points and length = " << stroke->length << "\n\n";
#endif
		}
		//next startingDegree to start the remaining strokes with
		if (startingDegree == 1 && first)
		{
			startingDegree = 5+1; //start from this new degree and decrease 1-by-1; +1 is for the immediate -- just-below
			first = false;
		}
		if (startingDegree == 1 && ! first) //i hit 1 again but this time due to 1-by-1 decrements, which means time to stop
			startingDegree = -8; //flag to quit the outer loop
		startingDegree--; //when i first come here, it means try new starting points p above (in next iter) that has startingDegree neighbors
	} while (startingDegree != -9);

	//postprocessing to segment strokes[] further via doug-peucker, which has a more global point of view; i could omit this step by adjusting maxRatioAllowed value but this is a more principle/scientific solution
	if (doDougPeuckerPostprocessing)
	{
		int initialStrokeSize = (int) strokes.size(); //use this one in for-loop headers as the strokes[] will be appended by new strokes
		for (int s = 0; s < initialStrokeSize; s++) //for each stroke
		{
			vector< double > polyline; //original polyline, assume not empty
			for (int p = 0; p < (int) strokes[s]->pnts.size(); p++)
			{
				//prepare polyline for douglas-peucker algo
				polyline.push_back(pnts[ strokes[s]->pnts[p] ]->coords[0]);
				polyline.push_back(pnts[ strokes[s]->pnts[p] ]->coords[1]);
			}
			////////////////////// douglas-peucker line simplification algo, which acts as a corner-detector on my input stroke(s) //////////////////////
			double tolerance = 1.2f*spacing;   //point-to-edge distance tolerance
			double* result = new double [ (int) polyline.size() ]; //make sure the result array is large enough
			int nCorners = 0;

			//simplify the 2d polyline [see tools\doug-peucker simplification\ w/ a better code and explanation]
			psimpl::simplify_douglas_peucker<2>(polyline.begin(), polyline.end(), tolerance, result, nCorners); //result[] has the desired corner coords via call-by-ref; nCorners keeps # corners detected

			/*cout<<"\ninputttttttt\n";
			for (int v = 0; v < (int) polyline.size(); v += 2)
				cout << polyline[v] << ", " << polyline[v+1] << "\t(pnts[" << strokes[s]->pnts[v/2] << "])\n";
			cout<<"outtttttttt\tnCorners = " << nCorners << "\n";
			for (int v = 0; v < nCorners; v++)
				cout << result[2*v] << ", " << result[2*v+1] << "\n";//*/

			//process the output, which is stored in result[] as x1, y1, x2, y2, x3, y3, .. where 2 consecutive 2d coords make 1 corner vertex in 2D
			for (int c = 0; c < nCorners; c++)
			{
				if (c > 0 && c < nCorners-1) //first and last corner are the endpnts of this stroke for sure; so do nothing for them
				{
					/////////////////////////////// creation of 2 substrokes /////////////////////////////
					pnts[ strokes[s]->pnts[0] ]->endpnt = pnts[ strokes[s]->pnts[ (int) strokes[s]->pnts.size()-1 ] ]->endpnt = true;
					//result[v], result[v+1] is the next new corner (nc), which is an existing SketchPoint; find this SketchPoint of interest
					SketchPnt* nc = NULL;
					for (int p = 0; p < (int) strokes[s]->pnts.size(); p++)
					{
						if (pnts[ strokes[s]->pnts[p] ]->coords[0] == result[2*c] && pnts[ strokes[s]->pnts[p] ]->coords[1] == result[2*c+1])
						{
#ifdef PRINT_DETAILS
							cout << "new corner detected: strokes[s]->pnts[" << p << "] = pnts[" << strokes[s]->pnts[p] << "]\n";
#endif
							nc = pnts[ strokes[s]->pnts[p] ];
							nc->spidx = p;
							nc->endpnt = true; //info: pnts[ strokes[s]->pnts[p] ]->endpnt becomes true after nc update here
							break;
						}
					}
					if (nc->spidx == 0 || nc->spidx == (int) strokes[s]->pnts.size())
					{
						cout << "new corner idx must be intermediate in (0, last-1), exclusive\n";
						exit(0);
					}
					//let nc split strokes[s] into 2 sub-strokes whose 1 endpnt is nc and the other is found by traversing along strokes[s] 'till an endpnt is hit
					/////////////// time-backward substroke ////////////////////////
					Edge* substroke1 = new Edge((int) strokes.size(), false);
					int sp = nc->spidx; //idx to strokes[s]->pnts[]
//cout << sp << "sp\n";
					while (true) //populate pnts[] of this new substroke 'till an existing endpnt=true set above is reached; go backwards in time (idx--)
					{
						substroke1->pnts.push_back( strokes[s]->pnts[sp] );
//cout << strokes[s]->pnts[sp] << " added\n";
						pnts[ strokes[s]->pnts[sp] ]->strokeID = (int) strokes.size();
						if ((int) substroke1->pnts.size() > 1)
							substroke1->length += distanceBetween2D(pnts[ substroke1->pnts[(int) substroke1->pnts.size() - 1] ]->coords, pnts[ substroke1->pnts[(int) substroke1->pnts.size() - 2] ]->coords);
						else if ((int) substroke1->pnts.size() == 1)
							substroke1->v1i = strokes[s]->pnts[sp]; //first point (idx to pnts[])
						sp--;
//cout << sp << "iiiiiiiiiiiiiiiisp\n\n";
						if (pnts[ strokes[s]->pnts[sp] ]->endpnt) //endpnt reached
						{
							//repeat stuff above for this last point and break
							substroke1->pnts.push_back(strokes[s]->pnts[sp]);
							pnts[ strokes[s]->pnts[sp] ]->strokeID = (int) strokes.size();
							substroke1->length += distanceBetween2D(pnts[ substroke1->pnts[(int) substroke1->pnts.size() - 1] ]->coords, pnts[ substroke1->pnts[(int) substroke1->pnts.size() - 2] ]->coords);
							break;
						}
					}
					substroke1->v2i = strokes[s]->pnts[sp]; //last point (idx to pnts[])
					pnts[ substroke1->v1i ]->strokeList.push_back(substroke1->idx);
					pnts[ substroke1->v2i ]->strokeList.push_back(substroke1->idx);
					pnts[ substroke1->v1i ]->endpnt = pnts[ substroke1->v2i ]->endpnt = true;
					strokes.push_back(substroke1);
					/////////////// time-backward substroke ends ////////////////////////
					/////////////// time-forward substroke ////////////////////////
					Edge* substroke2 = new Edge((int) strokes.size(), false);
					sp = nc->spidx + 1; //+1 'cos sp = nc->spidx already added to substroke1
					while (true) //populate pnts[] of this new substroke 'till an existing endpnt=true set above is reached; go forwards in time (idx++)
					{
						substroke2->pnts.push_back(strokes[s]->pnts[sp]);
						pnts[ strokes[s]->pnts[sp] ]->strokeID = (int) strokes.size();
						if ((int) substroke2->pnts.size() > 1)
							substroke2->length += distanceBetween2D(pnts[ substroke2->pnts[(int) substroke2->pnts.size() - 1] ]->coords, pnts[ substroke2->pnts[(int) substroke2->pnts.size() - 2] ]->coords);
						else if ((int) substroke2->pnts.size() == 1)
							substroke2->v1i = strokes[s]->pnts[sp]; //first point (idx to pnts[])
						sp++;
						if (pnts[ strokes[s]->pnts[sp] ]->endpnt) //endpnt reached
						{
							//repeat stuff above for this last point and break
							substroke2->pnts.push_back(strokes[s]->pnts[sp]);
							pnts[ strokes[s]->pnts[sp] ]->strokeID = (int) strokes.size();
							substroke2->length += distanceBetween2D(pnts[ substroke2->pnts[(int) substroke2->pnts.size() - 1] ]->coords, pnts[ substroke2->pnts[(int) substroke2->pnts.size() - 2] ]->coords);
							break;
						}
					}
					substroke2->v2i = strokes[s]->pnts[sp]; //last point (idx to pnts[])
					pnts[ substroke2->v1i ]->strokeList.push_back(substroke2->idx);
					pnts[ substroke2->v2i ]->strokeList.push_back(substroke2->idx);
					pnts[ substroke2->v1i ]->endpnt = pnts[ substroke2->v2i ]->endpnt = true;
					strokes.push_back(substroke2);
					/////////////// time-forward substroke ends ////////////////////////
					/////////////////////////////// creation of 2 substrokes ends /////////////////////////////
				}
			}
			delete [] result;
			////////////////////// douglas-peucker line simplification algo, which acts as a corner-detector on my input stroke(s) ends //////////////////////
		}
		/////////////////////////////// deletion of splitted big stroke(s) /////////////////////////////
		//addition of new substrokes above makes splitted strokes extra and bad; so remove a stroke if it has 3+ endpnts on it
		for (int s = 0; s < (int) strokes.size(); s++)
		{
			int nEndpnts = 0;
			for (int p = 0; p < (int) strokes[s]->pnts.size(); p++)
				nEndpnts += (pnts[ strokes[s]->pnts[p] ]->endpnt ? 1 : 0);
			if (nEndpnts > 2)
			{
				strokes[s]->removed = true; //explicit deletion w/ erase() causes lots of index shifts; so keep strokes[s] physically but ignore throughout the program
#ifdef PRINT_DETAILS
				cout << "splitted and to-be-ignored: strokes" << s << " w/ " << (int) strokes[s]->pnts.size() << " points and length = " << strokes[s]->length << "\n";
#endif
			}
		}
		/////////////////////////////// deletion of splitted big stroke(s) ends /////////////////////////////
	}

	//set initial endpnts[]; clear endpnt=true setting made during doug-peucker
	for (int s = 0; s < (int) pnts.size(); s++)
		pnts[s]->endpnt = false;
#ifdef PRINT_DETAILS
	cout << endl;
#endif
	for (int s = 0; s < (int) strokes.size(); s++)
	{
		if (strokes[s]->removed)
			continue; //big strokes splitted and made removed=true above
		endpnts.push_back(pnts[ strokes[s]->v1i ]);
		strokes[s]->ev1i = (int) endpnts.size() - 1; //v1i of this strokes[s] as a pointer to endpnts[]
		endpnts.push_back(pnts[ strokes[s]->v2i ]);
		strokes[s]->ev2i = (int) endpnts.size() - 1; //v2i of this strokes[s] as a pointer to endpnts[]
		pnts[ strokes[s]->v1i ]->endpnt = pnts[ strokes[s]->v2i ]->endpnt = true; //info: endpnts[(int)endpnts.size()-2 and -1]->endpnt becomes true after update here
#ifdef PRINT_DETAILS
		cout << "strokes" << strokes[s]->idx << ": endpnts" << strokes[s]->ev1i << "--" << strokes[s]->ev2i << " = pnts" << strokes[s]->v1i << "--" << strokes[s]->v2i << endl;
#endif
	}
}
int Mesh::getAnExtremeVert()
{
	//returns idx of an extreme vertex, e.g. on arm tip or toe of the mesh

	//vertex that is farthest away from arbitraryVert is definitely an extreme, where arbitraryVert can be anywhere on mesh (center, tips, head, ..)
	int arbitraryVert = 0, result = -1;
	dijkstraShortestPaths(arbitraryVert); //arbitraryVert.dSpanningUnn[] is now ready, dSpanning is bad due to unhealthy min/maxGeoDist values
	float maxDist = -INF;
	for (int s = 0; s < (int) verts.size(); s++) //fps selects from amongst all vertices
		if (verts[arbitraryVert]->dSpanningUnn[s] > maxDist)
		{
			maxDist = verts[arbitraryVert]->dSpanningUnn[s];
			result = s;
		}
	return result;
}
int Mesh::getAnExtremeVertSketch()
{
	//returns idx of an extreme vertex, e.g. on arm tip or toe of the sketch

	//vertex that is farthest away from arbitraryVert is definitely an extreme, where arbitraryVert can be anywhere on mesh (center, tips, head, ..)
	int arbitraryVert = 0, result = -1;
	float maxDist = -INF;
	if (usePnts) //pnts[]
	{
#ifdef USE_GOOD_CONTINUATION
		dijkstraShortestPathsSketchGC(arbitraryVert); //arbitraryVert.dSpanningUnn[] is now ready, dSpanning is bad due to unhealthy min/maxGeoDist values
#else
		dijkstraShortestPathsSketch(arbitraryVert); //arbitraryVert.dSpanningUnn[] is now ready, dSpanning is bad due to unhealthy min/maxGeoDist values
#endif
		for (int s = 0; s < (int) pnts.size(); s++) //fps selects from amongst all vertices
			if (pnts[arbitraryVert]->dSpanningUnn[s] > maxDist)
			{
				maxDist = pnts[arbitraryVert]->dSpanningUnn[s];
				result = s;
			}
	}
	else //endpnts[]
	{
		dijkstraShortestPathsSketchE(arbitraryVert); //arbitraryVert.dSpanningUnn[] is now ready, dSpanning is bad due to unhealthy min/maxGeoDist values
		for (int s = 0; s < (int) endpnts.size(); s++) //fps selects from amongst all vertices
			if (endpnts[arbitraryVert]->dSpanningUnn[s] > maxDist)
			{
				maxDist = endpnts[arbitraryVert]->dSpanningUnn[s];
				result = s;
			}
	}

	return result;
}
int Mesh::getIdxToEndpnts(int p)
{
	//returns index to endpnts[] that locates p

	for (int index = 0; index < (int) endpnts.size(); index++)
		if (endpnts[index]->idx == p) //->idx is an index to pnts and i'm looking for pnts[*] == p case
			return index;
	cout << "shouldn't reach here\n";
	exit(0);
	//return -1;
}

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;
}
float Mesh::matchWith(Mesh* sketch)
{
	//finds a good mapping b/w this mesh and sketch

	//scale invariance by making the maxGeoDist of this mesh equal to that of sketch's maxGeoDist
	float s = sketch->maxGeoDist / maxGeoDist;
#ifdef USE_GOOD_CONTINUATION
	s *= 8; //this manual update is necessary for visual scaling 'cos sketch->maxGeoDist is based on continuation costs whereas this->maxGeoDist is based on distance costs; so apple and orange scenario (not a problem during matching 'cos i'm using normalized dSpanning[] values there)
#endif
	scaleBy(s);
//dSpanning[] will not change w/ these new scaled coords/edgelengths but dSpanningUnn[]'ll change; since i don't use dSpanningUnn[] i'm skipping dijkstraShortestPaths() calls here
	maxGeoDist *= s;
#ifdef PRINT_DETAILS
	cout << "maxGeoDist on mesh after scale-normalization = " << maxGeoDist << "\navg edge lengths on sketch & mesh = " << sketch->avgEdgeLen << " & " << avgEdgeLen << "\n\n";
#endif
	//for better visualization translate com of the mesh to com of the sketch
	computeCenterOfMass(); //compute center of mass of this mesh
	translateBy(sketch->com - com); //geodesics not affected by this call
	computeCenterOfMass(); //redundant; update com info
	
	///////////// permutations ///////////////////////
	//no combinations below 'cos i use all N of N extremes for sure
	fillPermutationMtrxs(N);
	if ((int) samples.size() != N || (int) sketch->samples.size() != N)
	{
		cout << N << " vs. " << ((int) samples.size() != N ? (int) samples.size() : (int) sketch->samples.size()) << " incompatible extreme sizes; apply choose m from n trick (combinations)\n";
		exit(0);
	}

	int factR = factorial(N), twoR = 2*N;
	vector< int > minDistoMap;
	float minIsometryCost = INF;
	for (int r = 0; r < factR; r++) //r! reordering of mesh2 samples to be tested w/ the fixed mesh1 (this) quadruple
	{
		//refill tmpSafeCorrespTry according to the current reordering (one of k! permutations) of mesh2 samples
		vector< int > tmpSafeCorrespTry;
		int idx1 = 0, idx2 = 0;
		for (int y = 0; y < twoR; y++)
			if (y % 2 == 0)
				tmpSafeCorrespTry.push_back(samples[idx1++]); //even idxs from stable mesh1 side
			else //odd idxs for permutations from mesh2 side
			{
				if (N == 4)
					tmpSafeCorrespTry.push_back(sketch->samples[ perms4[r][idx2++] ]); 
				else if (N == 5)
					tmpSafeCorrespTry.push_back(sketch->samples[ perms5[r][idx2++] ]);
				else if (N == 6)
					tmpSafeCorrespTry.push_back(sketch->samples[ perms6[r][idx2++] ]);
				else if (N == 7)
					tmpSafeCorrespTry.push_back(sketch->samples[ perms7[r][idx2++] ]);
				else if (N == 8)
					tmpSafeCorrespTry.push_back(sketch->samples[ perms8[r][idx2++] ]);
				else //if (N == 9)
					tmpSafeCorrespTry.push_back(sketch->samples[ 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 = getIsometricDistortion(sketch, tmpSafeCorrespTry);
		if (isometryCostTry < minIsometryCost) //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
		{
#ifdef DENSE_MATCHING
			//new isometryCostTry based on dense matching will be used in enumaration of candidateMapsUnsorted; i still enumarate coarse maps b/w extremes, i.e. tmpSafeCorrespTry, but based on this new dense-friendly isometryCostTry
			isometryCostTry = getIsometricDistortionDense(mesh2, tmpSafeCorrespTry);
#endif
			minIsometryCost = isometryCostTry;
			minDistoMap = tmpSafeCorrespTry;
		}
	} //end of r
	cout << "\nfinal map isometric distortion = " << minIsometryCost << "\n\n";

	//transfer map in minDistoMap to matchIdxs which lets me visualize/evaluate the map
	for (int i = 0; i < (int) minDistoMap.size(); i += 2)
	{
		verts[ minDistoMap[i] ]->matchIdx = minDistoMap[i+1];
//cout << minDistoMap[i] << "<->" << verts[ minDistoMap[i] ]->matchIdx << endl;
	}
/*int a = 2, b = 6, tmp = verts[ minDistoMap[a] ]->matchIdx; //use 2x indices (from getMatchingSpheresSep) as minDistoMap[] holds 2 elements per match
verts[ minDistoMap[a] ]->matchIdx = verts[ minDistoMap[b] ]->matchIdx;
verts[ minDistoMap[b] ]->matchIdx = tmp;//*/

	return minIsometryCost;
}
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);
	}
}
void Mesh::scaleBy(float s)
{
	//scales this mesh by s

	for (int i = 0; i < (int) verts.size(); i++)
		verts[i]->coords *= s;
	//edge lengths also need to be updated since all verts have been scaled above
	updateEdgeLengths();
	//updateNormals();
}
void Mesh::updateEdgeLengths()
{
	//called after coords scaling which changes edge lengths

	edgeLenTotal = 0.0f;
	minEdgeLen = INF;
	maxEdgeLen = -INF;
	for (unsigned int e = 0; e < edges.size(); e++)
	{
		edges[e]->length = distanceBetween(verts[edges[e]->v1i]->coords, verts[edges[e]->v2i]->coords);
		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());
}
void Mesh::translateBy(SbVec3f t)
{
	//scales this mesh by s

	for (int i = 0; i < (int) verts.size(); i++)
		verts[i]->coords += t;
	//updateNormals();
}

void Mesh::additionalRotations()
{
	//apply additional x- y- and/or z-axis rotations to mesh

//	cout << "additional rotations for visual convenience..\n";
	float x = 0, y = 0, z = 0, theta; //radian rotation angle
	
	//y-axis rotation by theta degree
	theta = -90.0f * (PI / 180.0f);
	for (int i = 0; i < (int) verts.size(); i++)
	{
		x = verts[i]->coords[2]*sin(theta) + verts[i]->coords[0]*cos(theta);
		z = verts[i]->coords[2]*cos(theta) - verts[i]->coords[0]*sin(theta);
		verts[i]->coords[0] = x;
		verts[i]->coords[2] = z;
	}//*/

/*	//x-axis rotation by theta degree
	theta = -90.0f * (PI / 180.0f);
	for (int i = 0; i < (int) verts.size(); i++)
	{
		y = verts[i]->coords[1]*cos(theta) - verts[i]->coords[2]*sin(theta);
		z = verts[i]->coords[1]*sin(theta) + verts[i]->coords[2]*cos(theta);
		verts[i]->coords[1] = y;
		verts[i]->coords[2] = z;
	}//*/
	
/*	//z-axis rotation by theta degree
	theta = 90.0f * (PI / 180.0f);
	for (int i = 0; i < (int) verts.size(); i++)
	{
		float x = verts[i]->coords[0]*cos(theta) - verts[i]->coords[1]*sin(theta);
		float y = verts[i]->coords[0]*sin(theta) + verts[i]->coords[1]*cos(theta);
		verts[i]->coords[0] = x;
		verts[i]->coords[1] = y;
	}//*/

	//edge lengths are preserved after rotations but bounding box info might change
	//computeBoundingBox();
}

inline float Mesh::getIsometricDistortion(Mesh* sketch, 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 = (sketch->usePnts ? sketch->pnts[bvb]->dSpanning[s2] : sketch->endpnts[bvb]->dSpanning[s2]); //getIsometricDistortion() 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;
}
