/*

	CageMesh.cpp


*/

#include "CageMesh.h"
#include <iostream>


inline double CageMesh::calcGLCint(Point3DD p, Point3DD t1, Point3DD t2, Point3DD eta)
{
	
	Point3DD t2mt1 = t2-t1;
	Point3DD pmt1 = p-t1;
	Point3DD pmt2 = p-t2;

	double Spmt1 = pmt1.Size();
	double  Spmt2 = pmt2.Size();
	
	double  tempval =   ( Point3DD::ScalarProduct(t2mt1,pmt1))/ ((t2mt1.Size()*pmt1.Size()) );
	if (abs(tempval)>(1.0-0.00000000001))
		return 0.0;
	double  alpha = acos(tempval );

	if ( abs(alpha - 3.14159265358979) < 0.00000000001 || abs(alpha) < 0.00000000001 )
		return 0.0;


	tempval =   (  Point3DD::ScalarProduct(pmt1,pmt2))/(  (Spmt1*Spmt2));
	if (abs(tempval)>(1.0 - 0.00000000001))
		return 0.0;

	double  beta    = acos( tempval );

	double  Lam = pmt1.SquareOfSize() *sin(alpha)*sin(alpha);

	double  delta = 3.14159265358979-alpha;
	Point3DD pmeta = p-eta;
	double  c = pmeta.SquareOfSize();
	double  sqrtc = sqrt(c);
	// direct calculation
	double 	a=Lam;
	double  sqrta = sqrt(a);
	double  sx = sin(delta);
	double  cx = cos(delta);



	double  I=0.0;
	I=-0.5*Sign(sx)* ( 2*sqrtc*atan((sqrtc*cx) / (sqrt(a+c*sx*sx) ) )+sqrta*log(((sqrta*(1-2*c*cx/(c*(1+cx)+a+sqrta*sqrt(a+c*sx*sx)))))*(2*sx*sx/pow((1-cx),2))));

	sx = sin((delta-beta));
	cx = cos((delta-beta));
	
	double  II=0.0;
	II=-0.5*Sign(sx)* ( 2*sqrtc*atan((sqrtc*cx) / (sqrt(a+c*sx*sx) ) )+sqrta*(log((sqrta*(1-2*c*cx/(c*(1+cx)+a+sqrta*sqrt(a+c*sx*sx))))*(2*sx*sx/pow((1-cx),2)))));

	double  myInt = (-1.0/(4*3.14159265358979) )*abs( I - II -sqrtc*beta);
	return myInt;

}

bool CageMesh::calcGLCclosed(TexturedMesh	*objectMesh)
{
	if ( objectMesh->m_vertices.size() > 0 )
	{
		GLCv.resize(objectMesh->m_vertices.size());
		GLCf.resize(objectMesh->m_vertices.size());
		GLCfactors.resize(objectMesh->m_vertices.size());

		t_circum.resize(objectMesh->m_faces.size());

		for (int i = 0; i < (int)objectMesh->m_vertices.size(); ++i ) {
			GLCv[i].resize(m_vertices.size());
			GLCf[i].resize(m_faces.size());
		}
	}

			
	std::vector<Tripplet>::iterator it;
	std::cout << "\n";
	UINT cnt=-1;
	for ( it = m_faces.begin(); 
		it != m_faces.end(); 
		++it )
	{
		std::cout << cnt+2 << " of " << m_faces.size() << '\r';
		Point3D ot1t = m_vertices[it->GetFirstIndex()];
		Point3D ot2t = m_vertices[it->GetSecondIndex()];
		Point3D ot3t = m_vertices[it->GetThirdIndex()];
		Point3DD ot1(ot1t.GetX(),ot1t.GetY(),ot1t.GetZ());
		Point3DD ot2(ot2t.GetX(),ot2t.GetY(),ot2t.GetZ());
		Point3DD ot3(ot3t.GetX(),ot3t.GetY(),ot3t.GetZ());

		cnt++;
		Point3D tnt = m_face_normals[cnt];
		Point3DD tn(tnt.GetX(),tnt.GetY(),tnt.GetZ());
		real tA = m_faces_area[cnt];
		
		
		for (int i = 0; i < (int)objectMesh->m_vertices.size(); ++i ) 	{

			Point3D oetat = objectMesh->m_vertices[i];
			Point3DD oeta(oetat.GetX(),oetat.GetY(),oetat.GetZ());

		
			Point3DD t1 = ot1 - oeta;
			Point3DD t2 = ot2 - oeta;
			Point3DD t3 = ot3 - oeta;
			Point3DD eta(0,0,0);

			Point3DD p = Point3DD::ScalarProduct(t1-eta,tn)*tn;
			
		
			
			double  myInt1 = 0;
			double  mysign1 = Point3DD::ScalarProduct(Point3DD::VectorProduct(t1-p,t2-p),tn);
			mysign1 = Sign(mysign1);
			myInt1 = calcGLCint(p,t1,t2,eta);
			
			
			double  myInt2 = 0;
			double  mysign2 = Point3DD::ScalarProduct(Point3DD::VectorProduct(t2-p,t3-p),tn);
			mysign2=Sign(mysign2);
			myInt2 = calcGLCint(p,t2,t3,eta);
			
		
			double  myInt3 = 0;
			double  mysign3 = Point3DD::ScalarProduct(Point3DD::VectorProduct(t3-p,t1-p),tn);// det( [p t3 t1]  );
			mysign3 = Sign( mysign3);
			myInt3 = calcGLCint(p,t3,t1,eta);
			
			

			double I_t = -abs( mysign1*myInt1 + mysign2*myInt2 + mysign3*myInt3);
			double II = -I_t;//had minus

		
			//for triangle t' in N(t)
			double I_tp1 = (calcGLCint(eta,t2,t1,eta));
			double I_tp2 = (calcGLCint(eta,t3,t2,eta));
			double I_tp3 = (calcGLCint(eta,t1,t3,eta));


			Point3DD n1 = Point3DD::VectorProduct(t2-eta,t1-eta);  
			double sn1 = n1.Size();
			if (sn1 < 0.0000001){
				I_tp1=0.0;
			}
			else
				n1 = n1*(1.0/sn1);//  sqrt(mydot(n1,n1));
			
			Point3DD n2 = Point3DD::VectorProduct(t3-eta,t2-eta);
			double sn2 = n2.Size();
			if (sn2 < 0.0000001){
				I_tp2=0.0;
			}
			else
				n2 = n2*(1.0/sn2); // sqrt(mydot(n2,n2));


			Point3DD n3 = Point3DD::VectorProduct(t1-eta,t3-eta);  
			double sn3 = n3.Size();
			if (sn3 < 0.0000001){
				I_tp3=0.0;
			}
			else			
				n3 = n3*(1.0/sn3); //sqrt(mydot(n3,n3));

			
	
			Point3DD w = eta +  tn*I_t +(n1*I_tp1 + n2*I_tp2 + n3*I_tp3);

			// check if t1,t2,t3 are on the same plane - if so dont contribute to psi
			double ws = w.Size();
			if ( ws < 0.00001 ){ 
					if (mysign1 > 0 && mysign2> 0 && mysign3 >0){
						//the point is inside the face, remember factor 2
						GLCfactors[i]=1.0+cnt; //add the face index for later to recognize the face
					}
				}else{
		
	
				double I3 = Point3DD::ScalarProduct(n1,w)/Point3DD::ScalarProduct(n1,t3);
				double I1 = Point3DD::ScalarProduct(n2,w)/Point3DD::ScalarProduct(n2,t1);
				double I2 = Point3DD::ScalarProduct(n3,w)/Point3DD::ScalarProduct(n3,t2);

		
				GLCv[i][it->GetFirstIndex()]	=	GLCv[i][it->GetFirstIndex()]  + I1;
				GLCv[i][it->GetSecondIndex()]	=	GLCv[i][it->GetSecondIndex()] + I2;
				GLCv[i][it->GetThirdIndex()]	=	GLCv[i][it->GetThirdIndex()]  + I3;
			}
	

				//not using any normalization here only at deformation time
				GLCf[i][cnt] =  II;
			

		}
		
		
	}



	std::cout << "\n";

	std::cout << "Finished!" << "\n";

	return true;
}