Added testing of surface fitting on a sphere. Seems to work well.

cmath3d/TriangleMesh.cpp

@@ -526,6 +526,7 @@ bool ActiveMesh::GetPointOnQuadricSurface(const float x,const float 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);

cmath3d/TriangleMesh.h

@@ -46,6 +46,23 @@ class ActiveMesh
 	 */
 	int CalcQuadricSurface_Taubin(const int vertexID,
 	                              float (&coeffs)[10]);
+	int CalcQuadricSurface_sphere(const float radius,
+	                              const Vector3F& centre,
+	                              float (&coeffs)[10])
+	{
+		const float R=radius*radius;
+		coeffs[0]=-1.f + (centre.x*centre.x +centre.y*centre.y +centre.z*centre.z)/R;
+		coeffs[1]=-2.f*centre.x /R;
+		coeffs[2]=-2.f*centre.y /R;
+		coeffs[3]=-2.f*centre.z /R;
+		coeffs[4]=0.f;
+		coeffs[5]=0.f;
+		coeffs[6]=0.f;
+		coeffs[7]=1.f/R;
+		coeffs[8]=1.f/R;
+		coeffs[9]=1.f/R;
+		return(0);
+	}
 
 	/**
 	 * Calculate the 3rd coordinate (e.g., z) to complete a 3D point

src/graphics.cpp

@@ -149,7 +149,18 @@ void ActiveMesh::displayQuadricSurface(void)
 {
 	//get surface params
 	float surf_coeff[10];
-	CalcQuadricSurface_Taubin(VertexID,surf_coeff);
+	//CalcQuadricSurface_Taubin(VertexID,surf_coeff);
+	CalcQuadricSurface_sphere(5.0f,params.sceneCentre,surf_coeff);
+	std::cout << "\nsurface: " << surf_coeff[0] << " + "
+	          << surf_coeff[1] << "*x + "
+	          << surf_coeff[2] << "*y + "
+	          << surf_coeff[3] << "*z + "
+	          << surf_coeff[4] << "*xy + "
+	          << surf_coeff[5] << "*xz + "
+	          << surf_coeff[6] << "*yz + "
+	          << surf_coeff[7] << "*x*x + "
+	          << surf_coeff[8] << "*y*y + "
+	          << surf_coeff[9] << "*z*z = 0\n";
 
 	glPointSize(2.0f);
 	glBegin(GL_POINTS);
@@ -165,12 +176,16 @@ void ActiveMesh::displayQuadricSurface(void)
 
 	//FOR CYCLE BEGIN
 	//over all triangles
-	for (unsigned int t=0; t < neigsT.size(); ++t)
+	//for (unsigned int t=0; t < neigsT.size(); ++t)
+	for (unsigned int t=0; t < 2; ++t)
 	{
 		//determine triangle vertices
 		const Vector3FC& v1=Pos[ID[3*neigsT[t] +0]];
 		const Vector3FC& v2=Pos[ID[3*neigsT[t] +1]];
 		const Vector3FC& v3=Pos[ID[3*neigsT[t] +2]];
+		std::cout << "v1: (" << v1.x << "," << v1.y << "," << v1.z << ")\n";
+		std::cout << "v2: (" << v2.x << "," << v2.y << "," << v2.z << ")\n";
+		std::cout << "v3: (" << v3.x << "," << v3.y << "," << v3.z << ")\n";
 
 		//sweep (coarsely!) across current triangle
 		for (float c=0.1f; c <= 0.9f; c += 0.1f)
@@ -190,10 +205,14 @@ void ActiveMesh::displayQuadricSurface(void)
 			GetClosestPointOnQuadricSurface(point,surf_coeff);
 
 			//display the point
+			glColor3f(0.0f,1.0f,0.0f);
 			glVertex3f(point.x,point.y,point.z);
 			std::cout << "  (" << origPoint.x << "," << origPoint.y << "," << origPoint.z
 			          << ") -> ("
 			          << point.x << "," << point.y << "," << point.z << ")\n";
+
+			glColor3f(0.0f,1.0f,1.0f);
+			glVertex3f(origPoint.x,origPoint.y,origPoint.z);
 		}
 	}
 	//FOR CYCLE END
@@ -203,7 +222,7 @@ void ActiveMesh::displayQuadricSurface(void)
 
 void ActiveMesh::displayMesh(void)
 {
-	glColor3f(0.8f,0.8f,0.8f);
+	glColor3f(0.6f,0.6f,0.6f);
 
 	long unsigned int* id=ID.data();
 	for (unsigned int i=0; i < ID.size(); i+=3)
@@ -503,15 +522,15 @@ void reshape(int w, int h)
 void initializeGL(void)
 {
 	 int Argc=1;
-	 char msg[]="TRAgen: A Tool for Generation of Synthetic Time-Lapse Image Sequences of Living Cells";
+	 char msg[]="meshSurface GUI";
 	 char *Argv[10];
 	 Argv[0]=msg;
 
     glutInit(&Argc, Argv);
 
     glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA);
-    glutInitWindowSize(600, 600); 
-    glutInitWindowPosition(100, 100);
+    glutInitWindowSize(800, 800); 
+    glutInitWindowPosition(600, 100);
     glutCreateWindow(Argv[0]);
 
     glutKeyboardFunc(keyboard);

src/main.cpp

@@ -17,7 +17,8 @@ int main(void)
 
 	//load mesh
 	//char filename[]="../sample_input/samplecell.stl";
-	char filename[]="../sample_input/torus_100_30.stl";
+	//char filename[]="../sample_input/torus_100_30.stl";
+	char filename[]="../sample_input/sphere.stl";
 
 	int retval=mesh.ImportSTL(filename);
 	if (retval)