Merge pull request #690 from Gjacquenot/curvature

Update principal_curvature.cpp

Former-commit-id: ef98e0f5269ca3c61d6affac5f5e2bd5235c6c25

include/igl/principal_curvature.cpp

@@ -105,13 +105,12 @@ public:
     }
 
 
-    IGL_INLINE static Quadric fit(std::vector<Eigen::Vector3d> &VV, bool zeroDetCheck, bool svd)
+    IGL_INLINE static Quadric fit(const std::vector<Eigen::Vector3d> &VV)
     {
-      using namespace std;
       assert(VV.size() >= 5);
       if (VV.size() < 5)
       {
-        cerr << "ASSERT FAILED!" << endl;
+        std::cerr << "ASSERT FAILED! fit function requires at least 5 points: Only " << VV.size() << " were given." << std::endl;
         exit(0);
       }
 
@@ -160,8 +159,6 @@ public:
   bool localMode; /* Use local mode */
   bool projectionPlaneCheck; /* Check collected vertices on tangent plane */
   bool montecarlo;
-  bool svd; /* Use svd calculation instead of pseudoinverse */
-  bool zeroDetCheck; /* Check if the determinant is close to zero */
   unsigned int montecarloN;
 
   searchType st; /* Use either a sphere search or a k-ring search */
@@ -179,23 +176,23 @@ public:
   IGL_INLINE CurvatureCalculator();
   IGL_INLINE void init(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F);
 
-  IGL_INLINE void finalEigenStuff (int, std::vector<Eigen::Vector3d>, Quadric );
-  IGL_INLINE void fitQuadric (Eigen::Vector3d, std::vector<Eigen::Vector3d> ref, const  std::vector<int>& , Quadric *);
-  IGL_INLINE void applyProjOnPlane(Eigen::Vector3d, std::vector<int>, std::vector<int>&);
+  IGL_INLINE void finalEigenStuff(int, const std::vector<Eigen::Vector3d>&, Quadric&);
+  IGL_INLINE void fitQuadric(const Eigen::Vector3d&, const std::vector<Eigen::Vector3d>& ref, const std::vector<int>& , Quadric *);
+  IGL_INLINE void applyProjOnPlane(const Eigen::Vector3d&, const std::vector<int>&, std::vector<int>&);
   IGL_INLINE void getSphere(const int, const double, std::vector<int>&, int min);
   IGL_INLINE void getKRing(const int, const double,std::vector<int>&);
-  IGL_INLINE Eigen::Vector3d project(Eigen::Vector3d, Eigen::Vector3d, Eigen::Vector3d);
-  IGL_INLINE void computeReferenceFrame(int, Eigen::Vector3d, std::vector<Eigen::Vector3d>&);
-  IGL_INLINE void getAverageNormal(int, std::vector<int>, Eigen::Vector3d&);
-  IGL_INLINE void getProjPlane(int, std::vector<int>, Eigen::Vector3d&);
-  IGL_INLINE void applyMontecarlo(std::vector<int>&,std::vector<int>*);
+  IGL_INLINE Eigen::Vector3d project(const Eigen::Vector3d&, const Eigen::Vector3d&, const Eigen::Vector3d&);
+  IGL_INLINE void computeReferenceFrame(int, const Eigen::Vector3d&, std::vector<Eigen::Vector3d>&);
+  IGL_INLINE void getAverageNormal(int, const std::vector<int>&, Eigen::Vector3d&);
+  IGL_INLINE void getProjPlane(int, const std::vector<int>&, Eigen::Vector3d&);
+  IGL_INLINE void applyMontecarlo(const std::vector<int>&,std::vector<int>*);
   IGL_INLINE void computeCurvature();
-  IGL_INLINE void printCurvature(std::string outpath);
+  IGL_INLINE void printCurvature(const std::string& outpath);
   IGL_INLINE double getAverageEdge();
 
-  IGL_INLINE static int rotateForward (float *v0, float *v1, float *v2)
+  IGL_INLINE static int rotateForward (double *v0, double *v1, double *v2)
   {
-    float t;
+    double t;
 
     if (std::abs(*v2) >= std::abs(*v1) && std::abs(*v2) >= std::abs(*v0))
       return 0;
@@ -208,9 +205,9 @@ public:
     return 1 + rotateForward (v0, v1, v2);
   }
 
-  IGL_INLINE static void rotateBackward (int nr, float *v0, float *v1, float *v2)
+  IGL_INLINE static void rotateBackward (int nr, double *v0, double *v1, double *v2)
   {
-    float t;
+    double t;
 
     if (nr == 0)
       return;
@@ -223,18 +220,18 @@ public:
     rotateBackward (nr - 1, v0, v1, v2);
   }
 
-  IGL_INLINE static Eigen::Vector3d chooseMax (Eigen::Vector3d n, Eigen::Vector3d abc, float ab)
+  IGL_INLINE static Eigen::Vector3d chooseMax (Eigen::Vector3d n, Eigen::Vector3d abc, double ab)
   {
-    int i, max_i;
-    float max_sp;
+    int max_i;
+    double max_sp;
     Eigen::Vector3d nt[8];
 
     n.normalize ();
     abc.normalize ();
 
-    max_sp = - std::numeric_limits<float>::max();
+    max_sp = - std::numeric_limits<double>::max();
 
-    for (i = 0; i < 4; i++)
+    for (int i = 0; i < 4; ++i)
     {
       nt[i] = n;
       if (ab > 0)
@@ -290,7 +287,6 @@ public:
         max_i = i;
       }
     }
-
     return nt[max_i];
   }
 
@@ -315,8 +311,6 @@ IGL_INLINE CurvatureCalculator::CurvatureCalculator()
   this->montecarlo=false;
   this->montecarloN=0;
   this->kRing=3;
-  this->svd=true;
-  this->zeroDetCheck=true;
   this->curvatureComputed=false;
   this->expStep=true;
 }
@@ -337,7 +331,7 @@ IGL_INLINE void CurvatureCalculator::init(const Eigen::MatrixXd& V, const Eigen:
   igl::per_vertex_normals(V, F, face_normals, vertex_normals);
 }
 
-IGL_INLINE void CurvatureCalculator::fitQuadric (Eigen::Vector3d v, std::vector<Eigen::Vector3d> ref, const std::vector<int>& vv, Quadric *q)
+IGL_INLINE void CurvatureCalculator::fitQuadric(const Eigen::Vector3d& v, const std::vector<Eigen::Vector3d>& ref, const std::vector<int>& vv, Quadric *q)
 {
   std::vector<Eigen::Vector3d> points;
   points.reserve (vv.size());
@@ -354,17 +348,25 @@ IGL_INLINE void CurvatureCalculator::fitQuadric (Eigen::Vector3d v, std::vector<
     double z = vTang.dot(ref[2]);
     points.push_back(Eigen::Vector3d (x,y,z));
   }
-  *q = Quadric::fit (points, zeroDetCheck, svd);
+  if (points.size() < 5)
+  {
+    std::cerr << "ASSERT FAILED! fit function requires at least 5 points: Only " << points.size() << " were given." << std::endl;
+    *q = Quadric(0,0,0,0,0);
+  }
+  else
+  {
+    *q = Quadric::fit (points);
+  }
 }
 
-IGL_INLINE void CurvatureCalculator::finalEigenStuff (int i, std::vector<Eigen::Vector3d> ref, Quadric q)
+IGL_INLINE void CurvatureCalculator::finalEigenStuff(int i, const std::vector<Eigen::Vector3d>& ref, Quadric& q)
 {
 
-  double a = q.a();
-  double b = q.b();
-  double c = q.c();
-  double d = q.d();
-  double e = q.e();
+  const double a = q.a();
+  const double b = q.b();
+  const double c = q.c();
+  const double d = q.d();
+  const double e = q.e();
 
 //  if (fabs(a) < 10e-8 || fabs(b) < 10e-8)
 //  {
@@ -457,7 +459,7 @@ IGL_INLINE void CurvatureCalculator::getKRing(const int start, const double r, s
     vv.push_back(toVisit);
     if (distance<(int)r)
     {
-      for (unsigned int i=0; i<vertex_to_vertices[toVisit].size(); i++)
+      for (unsigned int i=0; i<vertex_to_vertices[toVisit].size(); ++i)
       {
         int neighbor=vertex_to_vertices[toVisit][i];
         if (!visited[neighbor])
@@ -488,13 +490,13 @@ IGL_INLINE void CurvatureCalculator::getSphere(const int start, const double r,
     int toVisit=queue->front();
     queue->pop_front();
     vv.push_back(toVisit);
-    for (unsigned int i=0; i<vertex_to_vertices[toVisit].size(); i++)
+    for (unsigned int i=0; i<vertex_to_vertices[toVisit].size(); ++i)
     {
       int neighbor=vertex_to_vertices[toVisit][i];
       if (!visited[neighbor])
       {
         Eigen::Vector3d neigh=vertices.row(neighbor);
-        float distance=(me-neigh).norm();
+        double distance=(me-neigh).norm();
         if (distance<r)
           queue->push_back(neighbor);
         else if ((int)vv.size()<min)
@@ -508,13 +510,13 @@ IGL_INLINE void CurvatureCalculator::getSphere(const int start, const double r,
     std::pair<int, double> cand=extra_candidates->top();
     extra_candidates->pop();
     vv.push_back(cand.first);
-    for (unsigned int i=0; i<vertex_to_vertices[cand.first].size(); i++)
+    for (unsigned int i=0; i<vertex_to_vertices[cand.first].size(); ++i)
     {
       int neighbor=vertex_to_vertices[cand.first][i];
       if (!visited[neighbor])
       {
         Eigen::Vector3d neigh=vertices.row(neighbor);
-        float distance=(me-neigh).norm();
+        double distance=(me-neigh).norm();
         extra_candidates->push(std::pair<int,double>(neighbor,distance));
         visited[neighbor]=true;
       }
@@ -525,16 +527,15 @@ IGL_INLINE void CurvatureCalculator::getSphere(const int start, const double r,
   free(visited);
 }
 
-IGL_INLINE Eigen::Vector3d CurvatureCalculator::project(Eigen::Vector3d v, Eigen::Vector3d  vp, Eigen::Vector3d ppn)
+IGL_INLINE Eigen::Vector3d CurvatureCalculator::project(const Eigen::Vector3d& v, const Eigen::Vector3d& vp, const Eigen::Vector3d& ppn)
 {
   return (vp - (ppn * ((vp - v).dot(ppn))));
 }
 
-IGL_INLINE void CurvatureCalculator::computeReferenceFrame(int i, Eigen::Vector3d normal, std::vector<Eigen::Vector3d>& ref )
+IGL_INLINE void CurvatureCalculator::computeReferenceFrame(int i, const Eigen::Vector3d& normal, std::vector<Eigen::Vector3d>& ref )
 {
 
-  Eigen::Vector3d longest_v=Eigen::Vector3d::Zero();
-  longest_v=Eigen::Vector3d(vertices.row(vertex_to_vertices[i][0]));
+  Eigen::Vector3d longest_v=Eigen::Vector3d(vertices.row(vertex_to_vertices[i][0]));
 
   longest_v=(project(vertices.row(i),longest_v,normal)-Eigen::Vector3d(vertices.row(i))).normalized();
 
@@ -546,25 +547,25 @@ IGL_INLINE void CurvatureCalculator::computeReferenceFrame(int i, Eigen::Vector3
   ref[2]=normal;
 }
 
-IGL_INLINE void CurvatureCalculator::getAverageNormal(int j, std::vector<int> vv, Eigen::Vector3d& normal)
+IGL_INLINE void CurvatureCalculator::getAverageNormal(int j, const std::vector<int>& vv, Eigen::Vector3d& normal)
 {
   normal=(vertex_normals.row(j)).normalized();
   if (localMode)
     return;
 
-  for (unsigned int i=0; i<vv.size(); i++)
+  for (unsigned int i=0; i<vv.size(); ++i)
   {
     normal+=vertex_normals.row(vv[i]).normalized();
   }
   normal.normalize();
 }
 
-IGL_INLINE void CurvatureCalculator::getProjPlane(int j, std::vector<int> vv, Eigen::Vector3d& ppn)
+IGL_INLINE void CurvatureCalculator::getProjPlane(int j, const std::vector<int>& vv, Eigen::Vector3d& ppn)
 {
   int nr;
-  float a, b, c;
-  float nx, ny, nz;
-  float abcq;
+  double a, b, c;
+  double nx, ny, nz;
+  double abcq;
 
   a = b = c = 0;
 
@@ -605,9 +606,9 @@ IGL_INLINE double CurvatureCalculator::getAverageEdge()
   double sum = 0;
   int count = 0;
 
-  for (int i = 0; i<faces.rows(); i++)
+  for (int i = 0; i<faces.rows(); ++i)
   {
-    for (short unsigned j=0; j<3; j++)
+    for (short unsigned j=0; j<3; ++j)
     {
       Eigen::Vector3d p1=vertices.row(faces.row(i)[j]);
       Eigen::Vector3d p2=vertices.row(faces.row(i)[(j+1)%3]);
@@ -623,14 +624,14 @@ IGL_INLINE double CurvatureCalculator::getAverageEdge()
 }
 
 
-IGL_INLINE void CurvatureCalculator::applyProjOnPlane(Eigen::Vector3d ppn, std::vector<int> vin, std::vector<int> &vout)
+IGL_INLINE void CurvatureCalculator::applyProjOnPlane(const Eigen::Vector3d& ppn, const std::vector<int>& vin, std::vector<int> &vout)
 {
-  for (std::vector<int>::iterator vpi = vin.begin(); vpi != vin.end(); ++vpi)
-    if (vertex_normals.row(*vpi) * ppn > 0.0f)
-      vout.push_back (*vpi);
+  for (std::vector<int>::const_iterator vpi = vin.begin(); vpi != vin.end(); ++vpi)
+    if (vertex_normals.row(*vpi) * ppn > 0.0)
+      vout.push_back(*vpi);
 }
 
-IGL_INLINE void CurvatureCalculator::applyMontecarlo(std::vector<int>& vin, std::vector<int> *vout)
+IGL_INLINE void CurvatureCalculator::applyMontecarlo(const std::vector<int>& vin, std::vector<int> *vout)
 {
   if (montecarloN >= vin.size ())
   {
@@ -639,24 +640,22 @@ IGL_INLINE void CurvatureCalculator::applyMontecarlo(std::vector<int>& vin, std:
   }
 
   float p = ((float) montecarloN) / (float) vin.size();
-  for (std::vector<int>::iterator vpi = vin.begin(); vpi != vin.end(); ++vpi)
+  for (std::vector<int>::const_iterator vpi = vin.begin(); vpi != vin.end(); ++vpi)
   {
     float r;
     if ((r = ((float)rand () / RAND_MAX)) < p)
     {
-      vout->push_back (*vpi);
+      vout->push_back(*vpi);
     }
   }
 }
 
 IGL_INLINE void CurvatureCalculator::computeCurvature()
 {
-  using namespace std;
-
   //CHECK che esista la mesh
-  size_t vertices_count=vertices.rows() ;
+  const size_t vertices_count=vertices.rows();
 
-  if (vertices_count <=0)
+  if (vertices_count ==0)
     return;
 
   curvDir=std::vector< std::vector<Eigen::Vector3d> >(vertices_count);
@@ -691,10 +690,9 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
         return;
     }
 
-    std::vector<Eigen::Vector3d> ref(3);
     if (vv.size()<6)
     {
-      std::cerr << "Could not compute curvature of radius " << scaledRadius << endl;
+      std::cerr << "Could not compute curvature of radius " << scaledRadius << std::endl;
       return;
     }
 
@@ -704,8 +702,7 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
       vvtmp.reserve (vv.size ());
       applyProjOnPlane (vertex_normals.row(i), vv, vvtmp);
       if (vvtmp.size() >= 6 && vvtmp.size()<vv.size())
-	  vv = vvtmp;
-
+        vv = vvtmp;
     }
 
 
@@ -723,7 +720,7 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
     }
     if (vv.size()<6)
     {
-      std::cerr << "Could not compute curvature of radius " << scaledRadius << endl;
+      std::cerr << "Could not compute curvature of radius " << scaledRadius << std::endl;
       return;
     }
     if (montecarlo)
@@ -737,6 +734,7 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
 
     if (vv.size()<6)
       return;
+    std::vector<Eigen::Vector3d> ref(3);
     computeReferenceFrame(i,normal,ref);
 
     Quadric q;
@@ -748,7 +746,7 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
   curvatureComputed=true;
 }
 
-IGL_INLINE void CurvatureCalculator::printCurvature(std::string outpath)
+IGL_INLINE void CurvatureCalculator::printCurvature(const std::string& outpath)
 {
   using namespace std;
   if (!curvatureComputed)
@@ -765,7 +763,7 @@ IGL_INLINE void CurvatureCalculator::printCurvature(std::string outpath)
 
   int vertices_count=vertices.rows();
   of << vertices_count << endl;
-  for (int i=0; i<vertices_count; i++)
+  for (int i=0; i<vertices_count; ++i)
   {
     of << curv[i][0] << " " << curv[i][1] << " " << curvDir[i][0][0] << " " << curvDir[i][0][1] << " " << curvDir[i][0][2] << " " <<
     curvDir[i][1][0] << " " << curvDir[i][1][1] << " " << curvDir[i][1][2] << endl;
@@ -792,12 +790,10 @@ IGL_INLINE void igl::principal_curvature(
   unsigned radius,
   bool useKring)
 {
-  using namespace std;
-
   if (radius < 2)
   {
     radius = 2;
-    cout << "WARNING: igl::principal_curvature needs a radius >= 2, fixing it to 2." << endl;
+    std::cout << "WARNING: igl::principal_curvature needs a radius >= 2, fixing it to 2." << std::endl;
   }
 
   // Preallocate memory
@@ -823,10 +819,8 @@ IGL_INLINE void igl::principal_curvature(
   cc.computeCurvature();
 
   // Copy it back
-  for (unsigned i=0; i<V.rows(); i++)
+  for (unsigned i=0; i<V.rows(); ++i)
   {
-    Eigen::Vector3d d1;
-    Eigen::Vector3d d2;
     PD1.row(i) << cc.curvDir[i][0][0], cc.curvDir[i][0][1], cc.curvDir[i][0][2];
     PD2.row(i) << cc.curvDir[i][1][0], cc.curvDir[i][1][1], cc.curvDir[i][1][2];
     PD1.row(i).normalize();
@@ -843,7 +837,7 @@ IGL_INLINE void igl::principal_curvature(
 
     if (PD1.row(i) * PD2.row(i).transpose() > 10e-6)
     {
-      cerr << "PRINCIPAL_CURVATURE: Something is wrong with vertex: i" << endl;
+      std::cerr << "PRINCIPAL_CURVATURE: Something is wrong with vertex: " << i << std::endl;
       PD1.row(i) *= 0;
       PD2.row(i) *= 0;
     }