Added first version of quadric surface fitting ala Taubin (incomplete).

cmath3d/TriangleMesh.cpp

 #include <iostream>
 #include <fstream>
 #include <map>
+#include <lapacke.h>
 
 #include "TriangleMesh.h"
+#include "../cmath3d_v/TriangleMesh_v.h"
 
 //'multiple' should be ideally 10^desired_decimal_accuracy
 int inline RoundTo(const float val, const float multiple=1000.f)
@@ -391,3 +393,104 @@ void ActiveMesh::ScaleMesh(const Vector3F& scale)
 		Pos[i].z*=scale.z;
 	}
 }
+
+
+// ============== surface fitting ============== 
+void 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)
+	for (unsigned int n=0; n < neigs.size(); ++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 < 10; ++j)  //for column
+		 for (int i=0; i < 10; ++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*10 +i; //Fortran
+
+			M[off]+= V[j]* V[i];
+			N[off]+=V1[j]*V1[i];
+			N[off]+=V2[j]*V2[i];
+			N[off]+=V3[j]*V3[i];
+		 }
+	}
+
+	//now, solve the generalized eigenvector of the matrix pair (M,N):
+	//MC = nNC
+	//
+	//M,N are 10x10 matrices constructed above,
+	//C is vector (infact, the coeff), n is scalar Lagrange multiplier
+	//
+	//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=10;
+	float w[10];
+	float work[1024];
+	lapack_int lwork=1024;
+	lapack_int info;
+	LAPACK_ssygv(&itype,&jobz,&uplo,&n,M,&n,N,&n,w,work,&lwork,&info);
+
+	std::cout << "info=" << info << "\n";
+	std::cout << "work(1)=" << work[0] << "\n";
+	
+
+
+	
+
+	
+
+
+
+
+}
+	                              
+
+
+
+
+
+
+
+
+

cmath3d/TriangleMesh.h

@@ -25,6 +25,25 @@ class ActiveMesh
                                       //the array length should be 3x the number of triangles in the mesh
 	std::vector<Vector3F> norm;        //list of vectors that are normal/orthogonal to the planes given by the individual triangles
 
+
+	/**
+	 * Determine coeffs that represent the local surface around vertex \e vertexID.
+	 *
+	 * Any point [x,y,z] on the curve satisfies the equation:
+	 * Scalar product of \e coeffs and [1,x,y,z,xy,xz,yz,x^2,y^2,z^2] = 0.
+	 * This also explains the order (and meaning) of the \e coeffs array.
+	 * Order was adopted from [1], eq (2), page 3.
+	 * 
+	 * 
+	 * [1]: Dong-Ming Yan, Wenping Wang, Yang Liu, Zhouwang Yang.
+	 * Variational mesh segmentation via quadric surface fitting.
+	 * Elsevier: Computer-Aided Design 44 (2012), pp. 1072-1082.
+	 */
+	void CalcQuadricSurface_Taubin(const int vertexID,
+	                               float (&coeffs)[10]);
+
+
+
 	///input is filename
 	///returns 0 on success, otherwise non-zero
 	int ImportSTL(const char* filename);

src/main.cpp

@@ -16,8 +16,8 @@ int main(void)
 	ParamsSetup();
 
 	//load mesh
-	char filename[]="../sample_input/samplecell.stl";
-	//char filename[]="../sample_input/torus_100_30.stl";
+	//char filename[]="../sample_input/samplecell.stl";
+	char filename[]="../sample_input/torus_100_30.stl";
 
 	int retval=mesh.ImportSTL(filename);
 	if (retval) 
@@ -34,6 +34,11 @@ int main(void)
 	mesh.CenterMesh(params.sceneCentre);
 
 	//work with mesh
+	float surf_coeff[10];
+	mesh.CalcQuadricSurface_Taubin(10,surf_coeff);
+
+/*
+	//render mesh
 	i3d::Image3d<i3d::GRAY16> mask;
 	mask.SetOffset(i3d::Offset(0.0,0.0,0.0));
 	//mask.SetResolution(i3d::Resolution(20.0,20.0,20.0)); //for torus
@@ -43,10 +48,13 @@ int main(void)
 	mask.GetVoxelData()=0;
 	mesh.RenderMask(mask);
 	mask.SaveImage("mesh.ics");
+*/
 
+/*
 	initializeGL();
 	loopGL();
 	closeGL();
+*/
 
 	return(0);
 }
@@ -56,8 +64,8 @@ void ParamsSetup(void)
 {
 	//set up the environment
    params.sceneOffset=Vector3d<float>(0.f);
-   //params.sceneSize=Vector3d<float>(15.f,15.f,15.f); //for torus
-   params.sceneSize=Vector3d<float>(150.f,150.f,150.f); //for sample cell
+   params.sceneSize=Vector3d<float>(15.f,15.f,15.f); //for torus
+   //params.sceneSize=Vector3d<float>(150.f,150.f,150.f); //for sample cell
 
 	params.sceneCentre=params.sceneSize;
 	params.sceneCentre/=2.0f;