TriangleMesh.cpp

int ActiveMesh::CalcQuadricSurface_Taubin(const int vertexID,
                                          float (&coeffs)[10])
{
	//determine some reasonable number of nearest neighbors
	std::vector< std::vector<size_t> > neigsLeveled;
	ulm::getVertexNeighbours(*this,vertexID,2,neigsLeveled);

	//make it flat...
	std::vector<size_t> neigs;
	for (unsigned int l=0; l < neigsLeveled.size(); ++l)
		for (unsigned int i=0; i < neigsLeveled[l].size(); ++i)
			neigs.push_back(neigsLeveled[l][i]);
	neigsLeveled.clear();

	//V be the vector of [x,y,z] combinations, a counterpart to coeffs
	float V[10],V1[10],V2[10],V3[10];

	//M be the matrix holding initially sum of (square matrices) V*V'
	//N is similar, just a sum of three different Vs is used
	float M[100],N[100];
	for (int i=0; i < 100; ++i) M[i]=N[i]=0.f;

	//over all neighbors (including the centre vertex itself),
	//
	//this order garuantees that array V will be relevant for input
	//vertexID after the cycles are over -- will become handy later
	for (signed int n=(int)neigs.size()-1; n >= 0; --n)
	{
		//shortcut to the current point
		const Vector3FC& v=Pos[neigs[n]];

		//get Vs for the given point 'v'
		V[0]=1.f;       V1[0]=0.f;      V2[0]=0.f;      V3[0]=0.f;
		V[1]=v.x;       V1[1]=1.f;      V2[1]=0.f;      V3[1]=0.f;
		V[2]=v.y;       V1[2]=0.f;      V2[2]=1.f;      V3[2]=0.f;
		V[3]=v.z;       V1[3]=0.f;      V2[3]=0.f;      V3[3]=1.f;
		V[4]=v.x*v.y;   V1[4]=v.y;      V2[4]=v.x;      V3[4]=0.f;
		V[5]=v.x*v.z;   V1[5]=v.z;      V2[5]=0.f;      V3[5]=v.x;
		V[6]=v.y*v.z;   V1[6]=0.f;      V2[6]=v.z;      V3[6]=v.y;
		V[7]=v.x*v.x;   V1[7]=2.f*v.x;  V2[7]=0.f;      V3[7]=0.f;
		V[8]=v.y*v.y;   V1[8]=0.f;      V2[8]=2.f*v.y;  V3[8]=0.f;
		V[9]=v.z*v.z;   V1[9]=0.f;      V2[9]=0.f;      V3[9]=2.f*v.z;

		//construct V*V' and add it to M and N
		for (int j=0; j < 9; ++j)  //for column
		 for (int i=0; i < 9; ++i) //for row;  note the order optimal for Fortran
		 {
			//C order                 (row-major):  M_i,j -> M[i][j] -> &M +i*STRIDE +j
			//Lapack/Fortran order (column-major):  M_i,j -> M[i][j] -> &M +j*STRIDE +i
			//const int off=i*10 +j; //C
			const int off=j*9 +i; //Fortran

			M[off]+= V[j+1]* V[i+1];
			N[off]+=V1[j+1]*V1[i+1];
			N[off]+=V2[j+1]*V2[i+1];
			N[off]+=V3[j+1]*V3[i+1];
		 }
	}

	//now, solve the generalized eigenvector of the matrix pair (M,N):
	//MC = nNC
	//
	//M,N are (reduced) 9x9 matrices constructed above,
	//C is vector (infact, the coeff), n is scalar Lagrange multiplier
	//
	//if M,N were (full-size) 10x10 matrices, the N would be singular
	//as the 1st row would contain only zeros,
	//it is therefore reduced to 9x9 sacrifing the first row
	//
	//according to netlib (Lapack) docs, http://www.netlib.org/lapack/lug/node34.html
	//type 1, Az=lBz -- A=M, B=N, z=C
	//function: SSYGV

/*

	lapack_int itype=1;
	char jobz='V';
	char uplo='U';
	lapack_int n=9;
	float w[10];
	float work[512];
	lapack_int lwork=512;
	lapack_int info;
	LAPACK_ssygv(&itype,&jobz,&uplo,&n,M,&n,N,&n,w,work,&lwork,&info);

	std::cout << "vertices considered: " << neigs.size() << "\n";
	std::cout << "info=" << info << " (0 is OK)\n";
	std::cout << "work(1)=" << work[0] << " (should be below 512)\n";
	
	//if some error, report it to the caller
	if (info != 0) return info;

	//M is now matrix of eigenvectors
	//it should hold (according to Lapack docs):
	//Z^T N Z = I where Z is one eigenvector, I is identity matrix
	//
	//w holds eigenvalues in ascending order

	//our result c[1]...c[9] is the eigenvector
	//corresponding to the smallest non-negative eigenvalue, so the j-th eigenvector
	int j=0;
	while (j < 9 && w[j] < 0.f) ++j;
	//have we found some non-negative eigenvalue?
	if (j == 9) return(-9999);

	//also:
	//the last missing coefficient c[0] we will determine by submitting
	//the given input vertex to the algebraic expresion of the surface
	//(given with coeffs) and equating it to zero:
	coeffs[0]=0.f;
	for (int i=0; i < 9; ++i)
	{
		coeffs[i+1]=M[j*9 +i];          //copy eigenvector
		coeffs[0]-=coeffs[i+1]*V[i+1];  //determine c[0]
	}

	std::cout << "w(j)=" << w[j] << ", j=" << j << "\n";

*/

	return 0;
}


bool ActiveMesh::GetPointOnQuadricSurface(const float x,const float y,
                                          float &z1, float &z2,
                                          const float (&coeffs)[10])
{
	const float a=coeffs[9];
	const float b=coeffs[3] +coeffs[5]*x +coeffs[6]*y;
	const float c=coeffs[0] +coeffs[1]*x +coeffs[2]*y
	             +coeffs[4]*x*y +coeffs[7]*x*x +coeffs[8]*y*y;

	const float sqArg=b*b - 4*a*c;
	if (sqArg < 0.f) return false;
	if (a == 0.f) return false;

	z1=(-b + sqrtf(sqArg)) / (2.f*a);
	z2=(-b - sqrtf(sqArg)) / (2.f*a);

	return true;
}


float ActiveMesh::GetClosestPointOnQuadricSurface(Vector3F& point,
                                                  const float (&coeffs)[10])
{
	//backup original input coordinate
	const float x=point.x;
	const float y=point.y;
	const float z=point.z;
	float tmp1,tmp2;

	//list of possible coordinates
	std::vector<Vector3F> pointAdepts;

	//took a pair of coordinates, calculate the third one
	//and make it an adept...
	if (GetPointOnQuadricSurface(x,y,tmp1,tmp2,coeffs))
	{
		pointAdepts.push_back(Vector3F(x,y,tmp1));
		pointAdepts.push_back(Vector3F(x,y,tmp2));
	}
	if (GetPointOnQuadricSurface(x,z,tmp1,tmp2,coeffs))
	{
		pointAdepts.push_back(Vector3F(x,tmp1,z));
		pointAdepts.push_back(Vector3F(x,tmp2,z));
	}
	if (GetPointOnQuadricSurface(y,z,tmp1,tmp2,coeffs))
	{
		pointAdepts.push_back(Vector3F(tmp1,y,z));
		pointAdepts.push_back(Vector3F(tmp2,y,z));
	}

	//are we doomed?
	if (pointAdepts.size() == 0)
		return (-999999.f);

	//find the closest
	int closestIndex=ChooseClosestPoint(pointAdepts,point);
	
	//calc distance to it
	point-=pointAdepts[closestIndex];
	tmp1=point.Len();

	//adjust the input/output point
	point=pointAdepts[closestIndex];

	return (tmp1);
}


int ActiveMesh::ChooseClosestPoint(const std::vector<Vector3F>& points,
                                   const Vector3F& point)
{
	int minIndex=-1;
	float minSqDist=9999999999999.f;

	Vector3F p;
	for (unsigned int i=0; i < points.size(); ++i)
	{
		p=point;
		p-=points[i];
		if (p.LenQ() < minSqDist)
		{
			minIndex=i;
			minSqDist=p.LenQ();
		}
	}

	return minIndex;
}