#include "Mesh.h"

#include <string>
#include <vector>


#include <iostream>
#include <fstream>
#include <sstream>


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


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

inline float distanceBetween(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]));
}

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

	//update the list of edges of which verts[v1i] is a member
	//for geodesic computations in dijkstraShortestPaths() i will discard interior edges (v.interiorEdgeList) as geodesic is defined as surface distance and interior edges are just noises distracting the true geodesics; interiorEdges however are necessary for volumetric regulation term in spring system and deformations
	if (interiorEdge)
	{
		verts[v1i]->interiorEdgeList.push_back(newcomer->idx);
		verts[v2i]->interiorEdgeList.push_back(newcomer->idx); //do same thing for [v2i]
	}
	else
	{
		verts[v1i]->edgeList.push_back(newcomer->idx);
		verts[v2i]->edgeList.push_back(newcomer->idx); //do same thing for [v2i]
	}

	edges.push_back(newcomer);

	//information
	if (newcomer->length < minEdgeLen)
		minEdgeLen = newcomer->length;
	if (newcomer->length > maxEdgeLen)
	{
		maxEdgeLen = newcomer->length;
		//cout << v1i << " " << v2i << "\t" << newcomer->length << endl;
	}
	edgeLenTotal += newcomer->length;
}

int* Mesh::getTrianglesSharedBy(Edge* bigEdge, int v1i, int v2i)
{
	//returns 2 tris adjacent to bigEdge (where endpnts of bigEdge is v1i and v2i)

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

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

	//assert(shTris[0] != shTris[1]);
	return shTris;
}


void Mesh::loadObj_face(char* meshFile){ //with tricky face data

	cout << "Mesh initializing (to " << meshFile << ")..\n";

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

	char type;	//type of line: vertex line (v) or face line (f)
	float a, b, c;	//for face lines casting these to int will suffice
	float temp1, temp2, temp3;//for reading faces
	minEdgeLen = INF;
	maxEdgeLen = -INF;
	edgeLenTotal = 0.0f;

	//skip comment

	while (1){

		char lineHeader[512];
		// read the first word of the line
		int res = fscanf(fPtr, "%s", lineHeader);
		if (res == EOF)
			break; // EOF = End Of File. Quit the loop.

		if (strcmp(lineHeader, "v") == 0){
			{
				fscanf(fPtr, "%c %f %f %f\n", &type, &a, &b, &c);
				//assume none of verts are duplicated initially
				//a *= SCALE_FACTOR;
				//b *= SCALE_FACTOR;
				//c *= SCALE_FACTOR;
				float* coords = new float[3];

				coords[0] = a;
				coords[1] = b;
				coords[2] = c;
				addVertexND(coords); //ND: no duplicate check
			}
		}
		else if (strcmp(lineHeader, "f") == 0){
			//f 8317//8317 8424//8424 8423//8423 format

			fscanf(fPtr, "%c %f//%f %f//%f %f//%f\n", &type, &a, &temp1, &b, &temp2, &c, &temp3);

			addTriangle((int)a - 1, (int)b - 1, (int)c - 1); //-1: obj indices start at 1
			//verts[(int)a - 1]->interior = verts[(int)b - 1]->interior = verts[(int)c - 1]->interior = false;

		}


	}
	avgEdgeLen = edgeLenTotal / ((int)edges.size());
	//computeBoundingBox();
	fclose(fPtr);
	//nBoundaryVerts = (int)verts.size();
	//nBoundaryTris = (int)tris.size();
	//nBoundaryEdges = (int)edges.size();
	//cout << "Mesh has " << (int)tris.size() << " tris, " << (int)verts.size() << " verts, " << (int)edges.size() << " edges\nInitialization done\n";
}

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

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

	char type;	//type of line: vertex line (v) or face line (f)
	float a, b, c;	//for face lines casting these to int will suffice
	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')
		{
			//assume none of verts are duplicated initially
			//	a *= SCALE_FACTOR;
			//	b *= SCALE_FACTOR;
			//	c *= SCALE_FACTOR;
			float* coords = new float[3];
			float* regulator = new float[3];

			coords[0] = a;
			coords[1] = b;
			coords[2] = c;
			addVertexND(coords); //ND: no duplicate check
		}
		else //if (type == 'f')
		{
			addTriangle((int)a - 1, (int)b - 1, (int)c - 1); //-1: obj indices start at 1

		}
	////	if (alpha > 0.0f)
	//for (int v = 0; v < (int)verts.size(); v++)
	//	if (verts[v]->triNeighbors.empty())
	//	{
	//	cout << "2D Dirichlet/ARAP regularization behaves weirdly w/ isolated vertices: get rid of v" << v << endl;
	//	exit(0);
	//	}
	avgEdgeLen = edgeLenTotal / ((int)edges.size());


	fclose(fPtr);

	//cout << "Mesh has " << (int)tris.size() << " tris, " << (int)verts.size() << " verts, " << (int)edges.size() << " edges\nInitialization done\n";

}


void Mesh::loadOff(char* meshFile)
{
	cout << "Mesh initializing (to " << meshFile << ")..\n";

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

	//cout << "reading  ..." << 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 verts coord section
	{
		fscanf(fPtr, "%f %f %f\n", &a, &b, &c);

		float* coords = new float[3];
		coords[0] = a;
		coords[1] = b;
		coords[2] = c;

		addVertexND(coords);

	}

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

	//color = new int[3];
	//	color[0] = 1;
	//color[1] = 1;
	//color[2] = 0;

	float totalLength = 0;
	for (unsigned int i = 0; i < edges.size(); i++){
		totalLength += edges[i]->length;
	}
	//avgEdgeLength = totalLength / edges.size();
	//draggedCoords = new float[3];

	//	avgEdgeLen = edgeLenTotal / ((int) edges.size());
	//	computeBoundingBox();
	fclose(fPtr);

	//cout << "Mesh has " << (int)tris.size() << " tris, " << (int)verts.size() << " verts, " << (int)edges.size() << " edges\nInitialization done\n";
}



void Mesh::readID(int id, string foldername, string format)
{
	string meshname = foldername;
	if (format == "hand" ){
		string meshname = foldername;
		meshname += to_string(id);
		meshname += ".off";
		//cout << "meshname: " << meshname << endl;

		char *cstr = &meshname[0u];
		this->loadOff(cstr);
	}
	else if (format == "cat" || format == "horse"){
		string meshname = foldername;
		
		meshname += to_string(id);
		meshname += ".obj";
		//cout << "meshname: " << meshname << endl;

		char *cstr = &meshname[0u];
		this->loadObj_face(cstr);
	}
	else if (format == "dancer"){
		string meshname = foldername;
		meshname += "mesh_0";
		if (id < 10){
			meshname += "00";
		}
		else if (id < 100){
			meshname += "0";
		}

		meshname += to_string(id);
		meshname += ".obj";
		//cout << "meshname: " << meshname << endl;

		char *cstr = &meshname[0u];
		this->loadObj(cstr);
	}
	else if (format == "scape"){
		meshname += "mesh0";
		if (id < 10){
			meshname += "0";
		}
		if (id > 50){
			meshname += to_string(id + 1);
		}
		else
			meshname += to_string(id);

		meshname += ".off";
		char *cstr = &meshname[0u];
		this->loadOff(cstr);
	}
	
}

void Mesh::addVertexND(float* c)
{
	//fastly add verts w/ coords c to mesh; ND: no duplication check

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


void  Mesh::triNormal(Triangle* tri, float(&normal)[3]){

	//Nf = (B - A) ×(C - A)
	float v21[3] = { 0, 0, 0 };
	float v31[3] = { 0, 0, 0 };
	v21[0] = verts[tri->v2i]->coords[0] - verts[tri->v1i]->coords[0];
	v21[1] = verts[tri->v2i]->coords[1] - verts[tri->v1i]->coords[1];
	v21[2] = verts[tri->v2i]->coords[2] - verts[tri->v1i]->coords[2];

	v31[0] = verts[tri->v3i]->coords[0] - verts[tri->v1i]->coords[0];
	v31[1] = verts[tri->v3i]->coords[1] - verts[tri->v1i]->coords[1];
	v31[2] = verts[tri->v3i]->coords[2] - verts[tri->v1i]->coords[2];

	//float normal[3];
	mycross(v21, v31, normal);

	float length = sqrt(normal[0] * normal[0] + normal[1] * normal[1] + normal[2] * normal[2]);

	normal[0] = normal[0] / length;
	normal[1] = normal[1] / length;
	normal[2] = normal[2] / length;
}



float dot(float * v1, float* v2){
	return (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]);
}

void mycross(float* v1, float* v2, float(&returnArray)[3]){

	returnArray[0] = v1[1] * v2[2] - v1[2] * v2[1];
	returnArray[1] = v1[2] * v2[0] - v1[0] * v2[2];
	returnArray[2] = v1[0] * v2[1] - v1[1] * v2[0];

}

float length(Vertex* v1){
	return sqrt(v1->coords[0] * v1->coords[0] + v1->coords[1] * v1->coords[1] + v1->coords[2] * v1->coords[2]);
}

void Mesh::LaplacianSmoothing(){
	//for each vertex
	for (unsigned int i = 0; i < verts.size(); i++){
		//check all the neighbors

		double noofneigbors = verts[i]->vertNeighborsSet.size();
		double temp[3] = { 0, 0, 0 };
		for (std::set<int>::iterator it = verts[i]->vertNeighborsSet.begin(); it != verts[i]->vertNeighborsSet.end(); ++it) {
			int vertexid = *it;

			temp[0] += verts[vertexid]->coords[0];
			temp[1] += verts[vertexid]->coords[1];
			temp[2] += verts[vertexid]->coords[2];
		}
		verts[i]->coords[0] = (1 / noofneigbors) * temp[0];
		verts[i]->coords[1] = (1 / noofneigbors) * temp[1];
		verts[i]->coords[2] = (1 / noofneigbors) * temp[2];

	}
}

void Mesh::inVertexNormal(Vertex* myvertex, float(&vertexnormal)[3]){

	for (unsigned int i = 0; i < myvertex->triNeighbors.size(); i++){
		float normaldir[3];
		//get tri normals
		triNormal(tris[myvertex->triNeighbors[i]], normaldir);
		vertexnormal[0] += normaldir[0];
		vertexnormal[1] += normaldir[1];
		vertexnormal[2] += normaldir[2];
	}
	//normalize
	float length = sqrt(vertexnormal[0] * vertexnormal[0] + vertexnormal[1] * vertexnormal[1] + vertexnormal[2] * vertexnormal[2]);

	//ray direction is opposite of vertex normal
	vertexnormal[0] = -vertexnormal[0] / length;
	vertexnormal[1] = -vertexnormal[1] / length;
	vertexnormal[2] = -vertexnormal[2] / length;

}




void normalize(float* f){
	float len = sqrt(f[0] * f[0] + f[1] * f[1] + f[2] * f[2]);
	f[0] /= len;
	f[1] /= len;
	f[2] /= len;
}


unsigned long long int factorial(int k){
	unsigned long long int result = 1;
	for (int i = 1; i <= k; i++){
		result *= i;
	}
	return result;

}

unsigned long long int combination(int k, int i){
	unsigned long long int comb = factorial(k) / ((factorial(k - i)* factorial(i)));
	//if (bezierprob)	cout << "comb:" << comb << " k:" <<k<<" i: "<<i <<endl;
	return comb;
}

void bezierNmesh(vector<Mesh*> meshlist, float amount, Mesh* meshinter, bool overwrite){
	int noofnodes = meshlist.size() - 1;
	////cout << "noofnodes: " << noofnodes << endl;
	//for (int n = 0; n <= noofnodes; n++){
	//	float factor = combination(noofnodes, n)*pow((1 - amount), (noofnodes - n))*pow(amount, n);
	//	//cout << "factor for: " << n << ": " << factor << endl;
	//	//cout << "combination: " << combination(noofnodes, n) << endl;
	//	//cout << "nth mesh: n: " << n << endl; 
	//	//cout << "pow((1 - amount), (noofnodes - n)) : (1-t)^:" << (noofnodes - n) << " t^" << n << endl;
	//}


	for (unsigned int i = 0; i < meshinter->verts.size(); i++){

		float newcoords[3] = { 0, 0, 0 };

		for (int n = 0; n <= noofnodes; n++){
			float factor = combination(noofnodes, n)     *    pow((1 - amount), (noofnodes - n))     *      pow(amount, n);
	
			newcoords[0] += factor*meshlist[n]->verts[i]->coords[0];
			newcoords[1] += factor*meshlist[n]->verts[i]->coords[1];
			newcoords[2] += factor*meshlist[n]->verts[i]->coords[2];

			/*if (newcoords[0] <= 0.0000001){
			cout << "sth is really wrong" << endl;
			}*/

		}
		if (overwrite == false){
			meshinter->verts[i]->coords[0] = newcoords[0];
			meshinter->verts[i]->coords[1] = newcoords[1];
			meshinter->verts[i]->coords[2] = newcoords[2];
		}
		else{
			meshinter->verts[i]->coords[0] = 0.5*(meshinter->verts[i]->coords[0] + newcoords[0]);
			meshinter->verts[i]->coords[1] = 0.5*(meshinter->verts[i]->coords[1] + newcoords[1]);
			meshinter->verts[i]->coords[2] = 0.5*(meshinter->verts[i]->coords[2] + newcoords[2]);
		}

	}

}
void Mesh::alignHelper(const vector<vector<int>>& partvertlist, int  stillpart, int movingpart){
	//cout << "now in still part:" << stillpart << " and moving:" << movingpart << endl;

	vector<int> edgelist1 = this->FindEdgeVertices(partvertlist[movingpart], partvertlist[stillpart]);
	vector<int> edgelist2 = this->FindEdgeVertices(partvertlist[stillpart], partvertlist[movingpart]);
	
	this->alignTwoParts(edgelist1, edgelist2, partvertlist[movingpart]);
	//this->alignRotation(edgelist1, edgelist2, partvertlist[movingpart]);


	vector <int> edges_concat;
	edges_concat.insert(edges_concat.end(), edgelist1.begin(), edgelist1.end());
	edges_concat.insert(edges_concat.end(), edgelist2.begin(), edgelist2.end());

	//cout << "edge fairing" << endl;
	for (int n = 0; n < 10; n++) {
		this->fairing_vert_set(edges_concat);
	}

}

vector<int> Mesh::FindEdgeVertices(const vector<int>& vertlist1, const vector<int>& vertlist2){
	//cout << "in fn findedgevertices" << endl;
	vector<int> edgevertices;
	bool foundone = false;
	for (unsigned int i = 0; i < vertlist1.size(); i++){
		for (std::set<int>::iterator it = verts[vertlist1[i]]->vertNeighborsSet.begin(); it != verts[vertlist1[i]]->vertNeighborsSet.end(); ++it) {
			int neighid = *it;
			//check each neighbor, if it is in the other list, vertex is an edge vertex
			for (unsigned int j = 0; j < vertlist2.size(); j++){
				if (neighid == vertlist2[j]){
					edgevertices.push_back(vertlist1[i]);
					foundone = true;
					break;
				}
			}
			if (foundone){//if there is one, dont bother checking others
				foundone = false;
				break;
			}
		}
	}
//	cout << "size of edge vert list in fn: " << edgevertices.size() << endl;
	if (edgevertices.size() == 0) cout << "!!!!!!!!!!!!!not found any edges!!!!!!!!!!!!!!" << endl;
//	cout << "got edges" << endl;
	return edgevertices;

}


void Mesh::AlignParts(const vector<vector<int>>& partvertlist, int noofparts, string format){
	if (format == "cat"){
		this->alignHelper(partvertlist, 3, 0);
		this->alignHelper(partvertlist, 3, 1);
		this->alignHelper(partvertlist, 3, 2);
		this->alignHelper(partvertlist, 3, 4);
	}
	else if (format == "horse"){
		this->alignHelper(partvertlist, 5, 0);
		this->alignHelper(partvertlist, 5, 1);
		this->alignHelper(partvertlist, 5, 2);
		this->alignHelper(partvertlist, 5, 3);
		this->alignHelper(partvertlist, 5, 4);
	}
	else if (format == "dancer"){
		this->alignHelper(partvertlist, 2, 0);
		this->alignHelper(partvertlist, 2, 1);
		this->alignHelper(partvertlist, 2, 3);
		this->alignHelper(partvertlist, 2, 4);
		this->alignHelper(partvertlist, 4, 5);
		this->alignHelper(partvertlist, 4, 6);
	}
	else if (noofparts == 5){
		/*		0- 2 - 1
				  / \
				 3   4
								*/
		//align 0 -2
		this->alignHelper(partvertlist, 2, 0);
		//align 1 -2
		this->alignHelper(partvertlist, 2, 1);
		//align 3 -2
		this->alignHelper(partvertlist, 2, 3);
		//align 4 -2
		this->alignHelper(partvertlist, 2, 4);
		//align 3 -4
		this->alignHelper(partvertlist, 3, 4);

	}
	else if (noofparts == 9){
		//keep torso (part 4), align others
		/*
				  0-7- 4 -6-5
					  /\
				   1 -   2
				   /	  \
			      3	      8

				   */

		////align 7 -4
		//this->alignHelper(partvertlist, 4, 7);
		////align 0 -7
		//this->alignHelper(partvertlist, 7, 0);

		////align 6 -5
		//this->alignHelper(partvertlist, 4, 6);
		////align 6 -4
		//this->alignHelper(partvertlist, 6, 5);

		////align 1-4
		//this->alignHelper(partvertlist, 4, 1);
		////align 1-3
		//this->alignHelper(partvertlist, 1, 3);

		////align 2-4
		//this->alignHelper(partvertlist, 4, 2);
		////align 2-8
		//this->alignHelper(partvertlist, 2, 8);

		////align 1-2
		//this->alignHelper(partvertlist, 1, 2);

		//align 7 -4
		this->alignHelper(partvertlist, 4, 7);
		//align 0 -7
		this->alignHelper(partvertlist, 7, 0);

		//align 6 -5
		this->alignHelper(partvertlist, 4, 6);
		//align 6 -4
		this->alignHelper(partvertlist, 6, 5);

		this->alignHelper(partvertlist, 4, 1);

		//align 2-4
		this->alignHelper(partvertlist, 4, 2);

		this->alignHelper(partvertlist, 1, 3);

		//align 2-8
		this->alignHelper(partvertlist, 2, 8);




	}
	else if (noofparts == 2){
		this->alignHelper(partvertlist, 0, 1);
	}
	else if (noofparts == 6){ //handdb
 		//  4 3 2 1
		//  \ | |/
		//     5  - 0
		//
		this->alignHelper(partvertlist, 5, 0);
		this->alignHelper(partvertlist, 5, 1);
		this->alignHelper(partvertlist, 5, 2);
		this->alignHelper(partvertlist, 5, 3);
		this->alignHelper(partvertlist, 5, 4);
	}
	else if (noofparts == 11){ //handdb
		// 8  6 4   2
		// \  | |  /
		// 10 7 9 3
		//  \ | |/
		//     5 - 1- 0
		//

		this->alignHelper(partvertlist, 5, 1);
		this->alignHelper(partvertlist, 1, 0);

		this->alignHelper(partvertlist, 5, 3);
		this->alignHelper(partvertlist, 3, 2);

		this->alignHelper(partvertlist, 5, 9);
		this->alignHelper(partvertlist, 9, 4);

		this->alignHelper(partvertlist, 5, 7);
		this->alignHelper(partvertlist, 7, 6);

		this->alignHelper(partvertlist, 5, 10);
		this->alignHelper(partvertlist, 10, 8);
	}

}
void Mesh::alignTwoParts(const vector<int>& edgetobealigned,const  vector<int>& edgestill,const vector<int>& movingpart){
//	cout << "aligning" << endl;

	int N = edgestill.size();
	//compute center of mass of not moving parts edges

	float q_[3] = { 0, 0, 0 };


	for (int i = 0; i < N; i++){
		//W.insert(i, i) = 1;//diagonal matrix n*n
		q_[0] += verts[edgestill[i]]->coords[0];
		q_[1] += verts[edgestill[i]]->coords[1];
		q_[2] += verts[edgestill[i]]->coords[2];
	}

	q_[0] /= N;
	q_[1] /= N;
	q_[2] /= N;

	//compute second center of mass of the moving part

	float p_[3] = { 0, 0, 0 };
	int N2 = edgetobealigned.size();

	//cout << "N2:" << N2 << endl;

	for (int i = 0; i < N2; i++){
		p_[0] += verts[edgetobealigned[i]]->coords[0];
		p_[1] += verts[edgetobealigned[i]]->coords[1];
		p_[2] += verts[edgetobealigned[i]]->coords[2];
	}
	p_[0] /= N2;
	p_[1] /= N2;
	p_[2] /= N2;

//	cout << "got center of masses" << endl;
//	cout << "q_[0]: " << q_ << "p_[0]: " << p_ << endl;


	for (int i = 0; i < movingpart.size(); i++){
		verts[movingpart[i]]->coords[0] += q_[0] - p_[0];
		verts[movingpart[i]]->coords[1] += q_[1] - p_[1];
		verts[movingpart[i]]->coords[2] += q_[2] - p_[2];
	}

}

void Mesh::fairing_N(int n){
	for (int i = 0; i < n; i++){
		this->fairing();
	}
}

void Mesh::fairing(){
	cout << "fairing" << endl;
	float lambda = 0.6307f, mu = -0.6732f;	//mu < -lambda and k = 1/lambda + 1/mu > 0, where k = .1, must hold


	for (unsigned int v = 0; v < verts.size(); v++)	//not v = 0; v < numOfVertsInUse; v++!!!
	{

	//	cout << "v: " << v << endl;
		float displacement[3] = { 0, 0, 0 };

		//float temp[3] = { 0, 0, 0 };
		//verts[v]->regulator = temp;	//totalDisplacement in can's code
		for (int i = 0; i < 3; i++) {
			verts[v]->regulator[i] = 0.0f;
		}

		float howManyNeighb = 0.0f, distance, totalWeight = 0.0f, weight;
		set< int, std::less< int > >::const_iterator neighbVertIdx;

		for (neighbVertIdx = verts[v]->vertNeighborsSet.begin(); neighbVertIdx != verts[v]->vertNeighborsSet.end(); neighbVertIdx++)
		{
			//if (! verts[(*neighbVertIdx)]->removed)
			if (verts[(*neighbVertIdx)] != 0)	//delete'd; so, prevent null ptr exception
			{

				//compute the displacement between these two vertices (v and neighb)
				for (int i = 0; i < 3;i++)
					displacement[i] = verts[(*neighbVertIdx)]->coords[i] - verts[v]->coords[i];
				//compute the length of this displacement
				distance = sqrt(dot(displacement, displacement));//length 

				//if distance is 0 then displacement is nullified
				if (distance == 0.0f)
				{
					//float temp2[3] = { 0, 0, 0 };
					//verts[v]->regulator = temp2;
					for (int i = 0; i < 3; i++)
						verts[v]->regulator[i] = 0.0f;
					break;
				}
				//compute the individual weight for this neighboring vertex
				weight = 1 / distance;
				//add this weight to the total weight
				totalWeight += weight;
				//multiply the weight with the displacement
				for (int i = 0; i < 3; i++)
					displacement[i] *= weight;
				//...and add the weighted displacement to the total displacement
				for (int i = 0; i < 3; i++)
					verts[v]->regulator[i] += displacement[i];	//totalDisplacement of can

				howManyNeighb++;
			}
		}
		//each non-removed vert must have 1+ neighb
		//assert(howManyNeighb != 0.0f); deceptive 'cos loop may break in 1st vert although v has 1+ vneighb
		//divide the final displacement value by the total weight.
		for (int i = 0; i < 3; i++) {
		//	cout << "v: " << v << " regulator: " << verts[v]->regulator[0] << " " << verts[v]->regulator[1] << " " << verts[v]->regulator[2] << endl;
			if (totalWeight < pow(10,-4)) cout << "0 total weight" << endl;
			verts[v]->regulator[i] /= totalWeight;
		//	cout << "totalweight: " << totalWeight << endl;
			

		}
		//cout << "v: " << v << " regulator: " << verts[v]->regulator[0] << " " << verts[v]->regulator[1] << " " << verts[v]->regulator[2] << endl;
	} //end of v
	/*for (unsigned int v = 0; v < verts.size(); v++) {
		cout << "v: " << v << " regulator: " << verts[v]->regulator[0] << " " << verts[v]->regulator[1] << " " << verts[v]->regulator[2] << endl;
	}*/

	//ApplyAllDisplacements() of can that uses verts.regulators set above
	for (unsigned int v = 0; v < verts.size(); v++)	//not v = 0; v < numOfVertsInUse; v++!!!
	{
		//for index consistency, verts are kept even if they're removed (ignore removed ones)
		//if (verts[v]->removed)	continue;
		if (verts[v] == 0)	continue;	//delete'd; so, prevent null ptr exception
		
	//	cout << "v: " << v << " regulator: " << verts[v]->regulator[0] << " " << verts[v]->regulator[1] << " " << verts[v]->regulator[2] << endl;


		verts[v]->coords[0] += lambda*verts[v]->regulator[0];	//fairing part1: vNew = v + lambda*deltaV
		verts[v]->coords[1] += lambda*verts[v]->regulator[1];	//fairing part1: vNew = v + lambda*deltaV
		verts[v]->coords[2] += lambda*verts[v]->regulator[2];	//fairing part1: vNew = v + lambda*deltaV
	}

	//nonshrinking effect by the same job above but with mu instead of lambda (part2)
	for (unsigned int v = 0; v < verts.size(); v++)	//not v = 0; v < numOfVertsInUse; v++!!!
	{

		if (verts[v] == 0)	continue;	//delete'd; so, prevent null ptr exception

		float displacement[3] = { 0,0,0 };

		//float temp3[3] = { 0, 0, 0 };
		//verts[v]->regulator = temp3;  //totalDisplacement in can's code

		for (int i = 0; i < 3; i++) {

			verts[v]->regulator[i] = 0.0f;
		}

		float howManyNeighb = 0.0f, distance, totalWeight = 0.0f, weight;
		//create a constant iterator
		set< int, std::less< int > >::const_iterator neighbVertIdx;
		for (neighbVertIdx = verts[v]->vertNeighborsSet.begin(); neighbVertIdx != verts[v]->vertNeighborsSet.end(); neighbVertIdx++)
		{
			//if (! verts[(*neighbVertIdx)]->removed)
			if (verts[(*neighbVertIdx)] != 0)	//delete'd; so, prevent null ptr exception
			{
				//assert((*neighbVertIdx) != v);	//i do not add 7 to the vertNeighborsSet of v7 in contract

				//compute the displacement between these two vertices (v and neighb)
				for (int i = 0; i < 3; i++)
					displacement[i] = verts[(*neighbVertIdx)]->coords[i] - verts[v]->coords[i];

				//compute the length of this displacement
				distance = sqrt(dot(displacement, displacement));//length 

				//if distance is 0 then displacement is nullified
				if (distance == 0.0f)
				{
				//	float temp4[3] = { 0, 0, 0 };
				//	verts[v]->regulator = temp4;

					for (int i = 0; i < 3; i++)
						verts[v]->regulator[i] = 0.0f;


					break;
				}
				//compute the individual weight for this neighboring vertex
				weight = 1 / distance;
				//add this weight to the total weight
				totalWeight += weight;
				//multiply the weight with the displacement
				for (int i = 0; i < 3; i++)
					displacement[i] *= weight;
				//...and add the weighted displacement to the total displacement
				for (int i = 0; i < 3; i++)
					verts[v]->regulator[i] += displacement[i];	//totalDisplacement of can


				howManyNeighb++;
			}
		}
		//each non-removed vert must have 1+ neighb
		//assert(howManyNeighb != 0.0f); deceptive 'cos loop may break in 1st vert although v has 1+ vneighb
		//divide the final displacement value by the total weight.


		for (int i = 0; i < 3; i++)
			verts[v]->regulator[i] /= totalWeight;

	} //end of v


	//ApplyAllDisplacements() of can that uses verts.regulators set above

	for (unsigned v = 0; v < verts.size(); v++)	//not v = 0; v < numOfVertsInUse; v++!!!
	{
		//for index consistency, verts are kept even if they're removed (ignore removed ones)
		//if (verts[v]->removed)	continue;



		if (verts[v] == 0)	continue;	//delete'd; so, prevent null ptr exception

		verts[v]->coords[0] += mu*verts[v]->regulator[0];	//fairing part2: vNew = v + mu*deltaV
		verts[v]->coords[1] += mu*verts[v]->regulator[1];	//fairing part2: vNew = v + mu*deltaV
		verts[v]->coords[2] += mu*verts[v]->regulator[2];	//fairing part2: vNew = v + mu*deltaV


	}
}





void Mesh::fairing_vert_set(const vector<int>& vert_set){

	//cout << "fairing" << endl;
	float lambda = 0.6307f, mu = -0.6732f;	//mu < -lambda and k = 1/lambda + 1/mu > 0, where k = .1, must hold


	for (unsigned int t = 0; t < vert_set.size(); t++)	//not v = 0; v < numOfVertsInUse; v++!!!
	{
		int v = vert_set[t];

		float displacement[3] = { 0, 0, 0 };

		for (int i = 0; i < 3; i++) {
			verts[v]->regulator[i] = 0.0f;
		}
			//totalDisplacement in can's code
		float howManyNeighb = 0.0f, distance, totalWeight = 0.0f, weight;
		set< int, std::less< int > >::const_iterator neighbVertIdx;

		for (neighbVertIdx = verts[v]->vertNeighborsSet.begin(); neighbVertIdx != verts[v]->vertNeighborsSet.end(); neighbVertIdx++)
		{
			//if (! verts[(*neighbVertIdx)]->removed)
			if (verts[(*neighbVertIdx)] != 0)	//delete'd; so, prevent null ptr exception
			{

				//compute the displacement between these two vertices (v and neighb)
				for (int i = 0; i < 3; i++)
					displacement[i] = verts[(*neighbVertIdx)]->coords[i] - verts[v]->coords[i];
				//compute the length of this displacement
				distance = sqrt(dot(displacement, displacement));//length 

				//if distance is 0 then displacement is nullified
				if (distance == 0.0f)
				{
					for (int i = 0; i < 3; i++) {
						verts[v]->regulator[i] = 0.0f;
					}
					break;
				}
				//compute the individual weight for this neighboring vertex
				weight = 1 / distance;
				//add this weight to the total weight
				totalWeight += weight;
				//multiply the weight with the displacement
				for (int i = 0; i < 3; i++)
					displacement[i] *= weight;
				//...and add the weighted displacement to the total displacement
				for (int i = 0; i < 3; i++)
					verts[v]->regulator[i] += displacement[i];	//totalDisplacement of can

				howManyNeighb++;
			}
		}
		//each non-removed vert must have 1+ neighb
		//assert(howManyNeighb != 0.0f); deceptive 'cos loop may break in 1st vert although v has 1+ vneighb
		//divide the final displacement value by the total weight.
		for (int i = 0; i < 3; i++)
			verts[v]->regulator[i] /= totalWeight;
	} //end of v

	//ApplyAllDisplacements() of can that uses verts.regulators set above
	for (unsigned int t = 0; t < vert_set.size(); t++)	//not v = 0; v < numOfVertsInUse; v++!!!
	{
		int v = vert_set[t];

		//for index consistency, verts are kept even if they're removed (ignore removed ones)
		//if (verts[v]->removed)	continue;
		if (verts[v] == 0)	continue;	//delete'd; so, prevent null ptr exception

		verts[v]->coords[0] += lambda * verts[v]->regulator[0];	//fairing part1: vNew = v + lambda*deltaV
		verts[v]->coords[1] += lambda * verts[v]->regulator[1];	//fairing part1: vNew = v + lambda*deltaV
		verts[v]->coords[2] += lambda * verts[v]->regulator[2];	//fairing part1: vNew = v + lambda*deltaV
	}

	//nonshrinking effect by the same job above but with mu instead of lambda (part2)
	for (unsigned int t = 0; t < vert_set.size(); t++)	//not v = 0; v < numOfVertsInUse; v++!!!
	{
		int v = vert_set[t];

		if (verts[v] == 0)	continue;	//delete'd; so, prevent null ptr exception

		float displacement[3] = { 0,0,0 };

		for (int i = 0; i < 3; i++) {
			verts[v]->regulator[i] = 0.0f;
		} //totalDisplacement in can's code
		float howManyNeighb = 0.0f, distance, totalWeight = 0.0f, weight;
		//create a constant iterator
		set< int, std::less< int > >::const_iterator neighbVertIdx;
		for (neighbVertIdx = verts[v]->vertNeighborsSet.begin(); neighbVertIdx != verts[v]->vertNeighborsSet.end(); neighbVertIdx++)
		{
			//if (! verts[(*neighbVertIdx)]->removed)
			if (verts[(*neighbVertIdx)] != 0)	//delete'd; so, prevent null ptr exception
			{
				//assert((*neighbVertIdx) != v);	//i do not add 7 to the vertNeighborsSet of v7 in contract

				//compute the displacement between these two vertices (v and neighb)
				for (int i = 0; i < 3; i++)
					displacement[i] = verts[(*neighbVertIdx)]->coords[i] - verts[v]->coords[i];

				//compute the length of this displacement
				distance = sqrt(dot(displacement, displacement));//length 

				//if distance is 0 then displacement is nullified
				if (distance == 0.0f)
				{
					for (int i = 0; i < 3; i++) {
						verts[v]->regulator[i] = 0.0f;
					}
					break;
				}
				//compute the individual weight for this neighboring vertex
				weight = 1 / distance;
				//add this weight to the total weight
				totalWeight += weight;
				//multiply the weight with the displacement
				for (int i = 0; i < 3; i++)
					displacement[i] *= weight;
				//...and add the weighted displacement to the total displacement
				for (int i = 0; i < 3; i++)
					verts[v]->regulator[i] += displacement[i];	//totalDisplacement of can


				howManyNeighb++;
			}
		}
		//each non-removed vert must have 1+ neighb
		//assert(howManyNeighb != 0.0f); deceptive 'cos loop may break in 1st vert although v has 1+ vneighb
		//divide the final displacement value by the total weight.
		for (int i = 0; i < 3; i++)
			verts[v]->regulator[i] /= totalWeight;

	} //end of v

	//ApplyAllDisplacements() of can that uses verts.regulators set above

	for (unsigned int t = 0; t < vert_set.size(); t++)	//not v = 0; v < numOfVertsInUse; v++!!!
	{
		int v = vert_set[t];

		//for index consistency, verts are kept even if they're removed (ignore removed ones)
		//if (verts[v]->removed)	continue;
		if (verts[v] == 0)	continue;	//delete'd; so, prevent null ptr exception

		verts[v]->coords[0] += mu * verts[v]->regulator[0];	//fairing part2: vNew = v + mu*deltaV
		verts[v]->coords[1] += mu * verts[v]->regulator[1];	//fairing part2: vNew = v + mu*deltaV
		verts[v]->coords[2] += mu * verts[v]->regulator[2];	//fairing part2: vNew = v + mu*deltaV


	}
}

void Mesh::linearinterpolation(Mesh* mesh1, Mesh* mesh2, float amount){

	for (unsigned int i = 0; i < verts.size(); i++){
		verts[i]->coords[0] = amount*mesh1->verts[i]->coords[0] + (1 - amount)*mesh2->verts[i]->coords[0];
		verts[i]->coords[1] = amount*mesh1->verts[i]->coords[1] + (1 - amount)*mesh2->verts[i]->coords[1];
		verts[i]->coords[2] = amount*mesh1->verts[i]->coords[2] + (1 - amount)*mesh2->verts[i]->coords[2];
	}

}


void linearint(Mesh* meshlin1, Mesh* mesh1, Mesh*  mesh2, float  amount){
	for (unsigned int i = 0; i < mesh1->verts.size(); i++){
		meshlin1->verts[i]->coords[0] = (1 - amount)*mesh1->verts[i]->coords[0] + amount*mesh2->verts[i]->coords[0];
		meshlin1->verts[i]->coords[1] = (1 - amount)*mesh1->verts[i]->coords[1] + amount*mesh2->verts[i]->coords[1];
		meshlin1->verts[i]->coords[2] = (1 - amount)*mesh1->verts[i]->coords[2] + amount*mesh2->verts[i]->coords[2];
	}
}
void Mesh::WriteOff(string filename) {
	std::ofstream out;

	out.open(filename, std::ios::app);
	if (out.is_open())
	{
		out << "OFF\n";
		out << verts.size() << " " << tris.size() << " " << edges.size() << "\n";
		for (unsigned int i = 0; i < verts.size(); i++) {
			out << verts[i]->coords[0] << " " << verts[i]->coords[1] << " " << verts[i]->coords[2] << "\n";
		}
		for (unsigned int i = 0; i < tris.size(); i++) {
			out << "3 " << tris[i]->v1i << " " << tris[i]->v2i << " " << tris[i]->v3i << "\n";
		}
	}
	out.close();
}