Merge branch 'vskhitkov-principal_curvature_fix' into dev

Former-commit-id: ad42bf1d2cc1c99c30cfca46fca043a56aef2604

include/igl/principal_curvature.cpp

@@ -692,8 +692,8 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
 
     if (vv.size()<6)
     {
-      std::cerr << "Could not compute curvature of radius " << scaledRadius << std::endl;
-      return;
+      //std::cerr << "Could not compute curvature of radius " << scaledRadius << std::endl;
+      continue;
     }
 
 
@@ -720,8 +720,8 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
     }
     if (vv.size()<6)
     {
-      std::cerr << "Could not compute curvature of radius " << scaledRadius << std::endl;
-      return;
+      //std::cerr << "Could not compute curvature of radius " << scaledRadius << std::endl;
+      continue;
     }
     if (montecarlo)
     {
@@ -773,6 +773,92 @@ IGL_INLINE void CurvatureCalculator::printCurvature(const std::string& outpath)
 
 }
 
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedPD1,
+  typename DerivedPD2,
+  typename DerivedPV1,
+  typename DerivedPV2,
+  typename Index>
+IGL_INLINE void igl::principal_curvature(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedPD1>& PD1,
+  Eigen::PlainObjectBase<DerivedPD2>& PD2,
+  Eigen::PlainObjectBase<DerivedPV1>& PV1,
+  Eigen::PlainObjectBase<DerivedPV2>& PV2,
+  std::vector<Index>& bad_vertices,
+  unsigned radius,
+  bool useKring)
+{
+
+  if (radius < 2)
+  {
+    radius = 2;
+    std::cout << "WARNING: igl::principal_curvature needs a radius >= 2, fixing it to 2." << std::endl;
+  }
+
+  // Preallocate memory
+  PD1.resize(V.rows(),3);
+  PD2.resize(V.rows(),3);
+
+  // Preallocate memory
+  PV1.resize(V.rows(),1);
+  PV2.resize(V.rows(),1);
+
+  // Precomputation
+  CurvatureCalculator cc;
+  cc.init(V.template cast<double>(),F.template cast<int>());
+  cc.sphereRadius = radius;
+
+  if (useKring)
+  {
+    cc.kRing = radius;
+    cc.st = K_RING_SEARCH;
+  }
+
+  // Compute
+  cc.computeCurvature();
+
+  // Copy it back
+  for (unsigned i=0; i<V.rows(); ++i)
+  {
+    if (!cc.curv[i].empty())
+    {
+      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();
+      PD2.row(i).normalize();
+
+      if (std::isnan(PD1(i,0)) || std::isnan(PD1(i,1)) || std::isnan(PD1(i,2)) || std::isnan(PD2(i,0)) || std::isnan(PD2(i,1)) || std::isnan(PD2(i,2)))
+      {
+        PD1.row(i) << 0,0,0;
+        PD2.row(i) << 0,0,0;
+      }
+
+      PV1(i) = cc.curv[i][0];
+      PV2(i) = cc.curv[i][1];
+
+      if (PD1.row(i) * PD2.row(i).transpose() > 10e-6)
+      {
+        bad_vertices.push_back((Index)i);
+
+        PD1.row(i) *= 0;
+        PD2.row(i) *= 0;
+      }
+    } else {
+      bad_vertices.push_back((Index)i);
+
+      PV1(i) = 0;
+      PV2(i) = 0;
+      PD1.row(i) << 0,0,0;
+      PD2.row(i) << 0,0,0;
+    }
+  }
+
+}
+
 template <
   typename DerivedV,
   typename DerivedF,
@@ -851,4 +937,5 @@ IGL_INLINE void igl::principal_curvature(
 template void igl::principal_curvature<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, unsigned int, bool);
 template void igl::principal_curvature<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, unsigned int, bool);
 template void igl::principal_curvature<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, unsigned int, bool);
-#endif
+template void igl::principal_curvature<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, int>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, std::vector<int, std::allocator<int> >&, unsigned int, bool);
+#endif

include/igl/principal_curvature.h

@@ -12,6 +12,8 @@
 #include <Eigen/Geometry>
 #include <Eigen/Dense>
 
+#include <vector>
+
 #include "igl_inline.h"
 //#include <igl/cotmatrix.h>
 //#include <igl/writeOFF.h>
@@ -35,6 +37,9 @@ namespace igl
   //   PV1 #V by 1 maximal curvature value for each vertex.
   //   PV2 #V by 1 minimal curvature value for each vertex.
   //
+  // Return value:
+  //   Function returns vector of indices of bad vertices if any.
+  //
   // See also: average_onto_faces, average_onto_vertices
   //
   // This function has been developed by: Nikolas De Giorgis, Luigi Rocca and Enrico Puppo.
@@ -43,11 +48,11 @@ namespace igl
   // Daniele Panozzo, Enrico Puppo, Luigi Rocca
   // GraVisMa, 2010
 template <
-  typename DerivedV, 
+  typename DerivedV,
   typename DerivedF,
-  typename DerivedPD1, 
-  typename DerivedPD2, 
-  typename DerivedPV1, 
+  typename DerivedPD1,
+  typename DerivedPD2,
+  typename DerivedPV1,
   typename DerivedPV2>
 IGL_INLINE void principal_curvature(
   const Eigen::PlainObjectBase<DerivedV>& V,
@@ -58,6 +63,26 @@ IGL_INLINE void principal_curvature(
   Eigen::PlainObjectBase<DerivedPV2>& PV2,
   unsigned radius = 5,
   bool useKring = true);
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedPD1,
+  typename DerivedPD2,
+  typename DerivedPV1,
+  typename DerivedPV2,
+  typename Index>
+IGL_INLINE void principal_curvature(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedPD1>& PD1,
+  Eigen::PlainObjectBase<DerivedPD2>& PD2,
+  Eigen::PlainObjectBase<DerivedPV1>& PV1,
+  Eigen::PlainObjectBase<DerivedPV2>& PV2,
+  std::vector<Index>& bad_vertices,
+  unsigned radius = 5,
+  bool useKring = true);
+
 }
 
 

tutorial/203_CurvatureDirections/main.cpp

@@ -39,7 +39,8 @@ int main(int argc, char *argv[])
   // Compute curvature directions via quadric fitting
   MatrixXd PD1,PD2;
   VectorXd PV1,PV2;
-  igl::principal_curvature(V,F,PD1,PD2,PV1,PV2);
+  std::vector<int> Z;
+  igl::principal_curvature(V,F,PD1,PD2,PV1,PV2,Z);
   // mean curvature
   H = 0.5*(PV1+PV2);
 

tutorial/CMakeLists.txt

@@ -32,13 +32,13 @@ else()
 endif()
 
 ### Choose which chapters to compile
-option(TUTORIALS_CHAPTER1 "Compile chapter 1" ON)
+option(TUTORIALS_CHAPTER1 "Compile chapter 1" OFF)
 option(TUTORIALS_CHAPTER2 "Compile chapter 2" ON)
-option(TUTORIALS_CHAPTER3 "Compile chapter 3" ON)
-option(TUTORIALS_CHAPTER4 "Compile chapter 4" ON)
-option(TUTORIALS_CHAPTER5 "Compile chapter 5" ON)
-option(TUTORIALS_CHAPTER6 "Compile chapter 6" ON)
-option(TUTORIALS_CHAPTER7 "Compile chapter 7" ON)
+option(TUTORIALS_CHAPTER3 "Compile chapter 3" OFF)
+option(TUTORIALS_CHAPTER4 "Compile chapter 4" OFF)
+option(TUTORIALS_CHAPTER5 "Compile chapter 5" OFF)
+option(TUTORIALS_CHAPTER6 "Compile chapter 6" OFF)
+option(TUTORIALS_CHAPTER7 "Compile chapter 7" OFF)
 
 # Store location of tutorial/shared directory
 set(TUTORIAL_SHARED_PATH ${CMAKE_CURRENT_SOURCE_DIR}/shared CACHE PATH "location of shared tutorial resources")