diff --git a/include/igl/principal_curvature.cpp b/include/igl/principal_curvature.cpp
index 5712130a0..aede8191e 100644
--- a/include/igl/principal_curvature.cpp
+++ b/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)
     {
@@ -780,7 +780,7 @@ template <
   typename DerivedPD2,
   typename DerivedPV1,
   typename DerivedPV2>
-IGL_INLINE void igl::principal_curvature(
+IGL_INLINE std::vector<size_t> igl::principal_curvature(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
   Eigen::PlainObjectBase<DerivedPD1>& PD1,
@@ -790,6 +790,8 @@ IGL_INLINE void igl::principal_curvature(
   unsigned radius,
   bool useKring)
 {
+  std::vector<size_t> bad_vertices;
+
   if (radius < 2)
   {
     radius = 2;
@@ -821,34 +823,47 @@ IGL_INLINE void igl::principal_curvature(
   // Copy it back
   for (unsigned i=0; i<V.rows(); ++i)
   {
-    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)))
+    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((size_t)i);
+
+        PD1.row(i) *= 0;
+        PD2.row(i) *= 0;
+      }
+    } else {
+      bad_vertices.push_back((size_t)i);
+
+      PV1(i) = 0;
+      PV2(i) = 0;
       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)
-    {
-      std::cerr << "PRINCIPAL_CURVATURE: Something is wrong with vertex: " << i << std::endl;
-      PD1.row(i) *= 0;
-      PD2.row(i) *= 0;
-    }
   }
 
+  return bad_vertices;
+
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
-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);
+template std::vector<size_t> 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 std::vector<size_t> 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 std::vector<size_t> 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
diff --git a/include/igl/principal_curvature.h b/include/igl/principal_curvature.h
index 958e58418..cb14e4fb4 100644
--- a/include/igl/principal_curvature.h
+++ b/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,13 +48,13 @@ 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(
+IGL_INLINE std::vector<size_t> principal_curvature(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
   Eigen::PlainObjectBase<DerivedPD1>& PD1,