Merge branch 'master' of github.com:libigl/libigl

Former-commit-id: 337e7d078a5b8ead5e0b6995242655a28ef9e689

include/igl/boundary_vertices_sorted.cpp

@@ -0,0 +1,101 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Stefan Brugger <stefanbrugger@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "boundary_vertices_sorted.h"
+
+#include "igl/tt.h"
+#include "igl/vf.h"
+
+IGL_INLINE void igl::boundary_vertices_sorted(
+    const Eigen::MatrixXd& V,
+    const Eigen::MatrixXi& F,
+    Eigen::VectorXi& b)
+{
+  std::vector<int> bnd;
+  bnd.clear();
+  std::vector<bool> isVisited(V.rows(),false);
+
+  Eigen::MatrixXi TT,TTi;
+  std::vector<std::vector<int> > VF, VFi;
+  igl::tt(V,F,TT,TTi);
+  igl::vf(V,F,VF,VFi);
+
+  // Extract one boundary edge
+  bool done = false;
+  for (int i = 0; i < TT.rows() && !done; i++)
+  {
+    for (int j = 0; j < TT.cols(); j++)
+    {
+      if (TT(i,j) < 0)
+      {
+        int idx1, idx2;
+        idx1 = F(i,j);
+        idx2 = F(i,(j+1) % F.cols());
+        bnd.push_back(idx1);
+        bnd.push_back(idx2);
+        isVisited[idx1] = true;
+        isVisited[idx2] = true;
+        done = true;
+        break;
+      }
+    }
+  }
+
+  // Traverse boundary
+  while(1)
+  {
+    bool changed = false;
+    int lastV;
+    lastV = bnd[bnd.size()-1];
+
+    for (int i = 0; i < VF[lastV].size(); i++)
+    {
+      int curr_neighbor = VF[lastV][i];
+
+      if (TT.row(curr_neighbor).minCoeff() < 0.) // Face contains boundary edge
+      {
+        int idx_lastV_in_face;
+        if (F(curr_neighbor,0) == lastV) idx_lastV_in_face = 0;
+        if (F(curr_neighbor,1) == lastV) idx_lastV_in_face = 1;
+        if (F(curr_neighbor,2) == lastV) idx_lastV_in_face = 2;
+
+        int idx_prev = (idx_lastV_in_face + F.cols()-1) % F.cols();
+        int idx_next = (idx_lastV_in_face + 1) % F.cols();
+        bool isPrevVisited = isVisited[F(curr_neighbor,idx_prev)];
+        bool isNextVisited = isVisited[F(curr_neighbor,idx_next)];
+
+        bool gotBndEdge = false;
+        int next_bnd_vertex;
+        if (!isNextVisited && TT(curr_neighbor,idx_lastV_in_face) < 0)
+        {
+          next_bnd_vertex = idx_next;
+          gotBndEdge = true;
+        }
+        else if (!isPrevVisited && TT(curr_neighbor,(idx_lastV_in_face+2) % F.cols()) < 0)
+        {
+          next_bnd_vertex = idx_prev;
+          gotBndEdge = true;
+        }
+
+        if (gotBndEdge)
+        {
+          changed = true;
+          bnd.push_back(F(curr_neighbor,next_bnd_vertex));
+          isVisited[F(curr_neighbor,next_bnd_vertex)] = true;
+          break;
+        }
+      }
+    }
+
+    if (!changed)
+      break;
+  }
+
+  b.resize(bnd.size());
+  for(unsigned i=0;i<bnd.size();++i)
+    b(i) = bnd[i];
+}

include/igl/boundary_vertices_sorted.h

@@ -0,0 +1,35 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Stefan Brugger <stefanbrugger@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_BOUNDARYVERTICESSORTED_H
+#define IGL_BOUNDARYVERTICESSORTED_H
+#include <igl/igl_inline.h>
+
+#include <Eigen/Dense>
+#include <vector>
+
+namespace igl
+{
+
+  // Compute sorted list of boundary vertices for mesh with single boundary.
+  //
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #V by dim list of mesh faces
+  // Outputs:
+  //   bnd   sorted list of boundary vertex indices
+  IGL_INLINE void boundary_vertices_sorted(
+  	const Eigen::MatrixXd& V, 
+  	const Eigen::MatrixXi& F, 
+    Eigen::VectorXi& bnd);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "boundary_vertices_sorted.cpp"
+#endif
+
+#endif

include/igl/map_vertices_to_circle.cpp

@@ -0,0 +1,59 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Stefan Brugger <stefanbrugger@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "map_vertices_to_circle.h"
+
+#include <Eigen/Sparse>
+
+#include "igl/cotmatrix.h"
+#include "igl/boundary_vertices_sorted.h"
+
+IGL_INLINE void igl::map_vertices_to_circle(
+  const Eigen::MatrixXd& V,
+  const Eigen::MatrixXi& F,
+  const Eigen::VectorXi& bnd,
+  Eigen::MatrixXd& UV)
+{
+  // Get sorted list of boundary vertices
+  std::vector<int> interior,map_ij;
+  map_ij.resize(V.rows());
+
+  std::vector<bool> isOnBnd(V.rows(),false);
+  for (int i = 0; i < bnd.size(); i++)
+  {
+    isOnBnd[bnd[i]] = true;
+    map_ij[bnd[i]] = i;
+  }
+
+  for (int i = 0; i < isOnBnd.size(); i++)
+  {
+    if (!isOnBnd[i])
+    {
+      map_ij[i] = interior.size();
+      interior.push_back(i);
+    }
+  }
+
+  // Map boundary to unit circle
+  std::vector<double> len(bnd.size());
+  len[0] = 0.;
+
+  for (int i = 1; i < bnd.size(); i++)
+  {
+    len[i] = len[i-1] + (V.row(bnd[i-1]) - V.row(bnd[i])).norm();
+  }
+  double total_len = len[len.size()-1] + (V.row(bnd[0]) - V.row(bnd[bnd.size()-1])).norm();
+
+  UV.resize(bnd.size(),2);
+  for (int i = 0; i < bnd.size(); i++)
+  {
+    double frac = len[i] * 2. * M_PI / total_len;
+    UV.row(map_ij[bnd[i]]) << cos(frac), sin(frac);
+  }
+
+}

include/igl/map_vertices_to_circle.h

@@ -0,0 +1,37 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Stefan Brugger <stefanbrugger@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_MAPVERTICESTOCIRCLE_H
+#define IGL_MAPVERTICESTOCIRCLE_H
+#include <igl/igl_inline.h>
+
+#include <Eigen/Dense>
+#include <vector>
+
+namespace igl
+{
+
+  // Map the vertices whose indices are in b on the unit circle.
+  //
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #V by dim list of mesh faces
+  //   b  #W list of vertex ids
+  // Outputs:
+  //   UV   #W by 2 list of 2D position on the unit circle for the vertices in b
+  IGL_INLINE void map_vertices_to_circle(
+  	const Eigen::MatrixXd& V,
+  	const Eigen::MatrixXi& F,
+    const Eigen::VectorXi& b,
+  	Eigen::MatrixXd& UV);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "map_vertices_to_circle.cpp"
+#endif
+
+#endif

include/igl/min_quad_with_fixed.h

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
+//
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_MIN_QUAD_WITH_FIXED_H
 #define IGL_MIN_QUAD_WITH_FIXED_H
@@ -55,9 +55,9 @@ namespace igl
     );
 
   // Solves a system previously factored using min_quad_with_fixed_precompute
-  // 
+  //
   // Template:
-  //   T  type of sparse matrix (e.g. double) 
+  //   T  type of sparse matrix (e.g. double)
   //   DerivedY  type of Y (e.g. derived from VectorXd or MatrixXd)
   //   DerivedZ  type of Z (e.g. derived from VectorXd or MatrixXd)
   // Inputs:
@@ -71,9 +71,9 @@ namespace igl
   // Returns true on success, false on error
   template <
     typename T,
-    typename DerivedB, 
+    typename DerivedB,
     typename DerivedY,
-    typename DerivedBeq, 
+    typename DerivedBeq,
     typename DerivedZ,
     typename Derivedsol>
   IGL_INLINE bool min_quad_with_fixed_solve(
@@ -86,9 +86,9 @@ namespace igl
   // Wrapper without sol
   template <
     typename T,
-    typename DerivedB, 
+    typename DerivedB,
     typename DerivedY,
-    typename DerivedBeq, 
+    typename DerivedBeq,
     typename DerivedZ>
   IGL_INLINE bool min_quad_with_fixed_solve(
     const min_quad_with_fixed_data<T> & data,

include/igl/principal_curvature.cpp

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// Copyright (C) 2013 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "principal_curvature.h"
 #include <iostream>
@@ -51,12 +51,12 @@ public:
   class Quadric
   {
   public:
-    
+
     IGL_INLINE Quadric ()
     {
       a() = b() = c() = d() = e() = 1.0;
     }
-    
+
     IGL_INLINE Quadric(double av, double bv, double cv, double dv, double ev)
     {
       a() = av;
@@ -65,46 +65,46 @@ public:
       d() = dv;
       e() = ev;
     }
-    
+
     IGL_INLINE double& a() { return data[0];}
     IGL_INLINE double& b() { return data[1];}
     IGL_INLINE double& c() { return data[2];}
     IGL_INLINE double& d() { return data[3];}
     IGL_INLINE double& e() { return data[4];}
-    
+
     double data[5];
-    
+
     IGL_INLINE double evaluate(double u, double v)
     {
       return a()*u*u + b()*u*v + c()*v*v + d()*u + e()*v;
     }
-    
+
     IGL_INLINE double du(double u, double v)
     {
       return 2.0*a()*u + b()*v + d();
     }
-    
+
     IGL_INLINE double dv(double u, double v)
     {
       return 2.0*c()*v + b()*u + e();
     }
-    
+
     IGL_INLINE double duv(double u, double v)
     {
       return b();
     }
-    
+
     IGL_INLINE double duu(double u, double v)
     {
       return 2.0*a();
     }
-    
+
     IGL_INLINE double dvv(double u, double v)
     {
       return 2.0*c();
     }
-    
-    
+
+
     IGL_INLINE static Quadric fit(std::vector<Eigen::Vector3d> &VV, bool zeroDetCheck, bool svd)
     {
       using namespace std;
@@ -118,67 +118,67 @@ public:
       Eigen::MatrixXd A(VV.size(),5);
       Eigen::MatrixXd b(VV.size(),1);
       Eigen::MatrixXd sol(5,1);
-      
+
       for(unsigned int c=0; c < VV.size(); ++c)
       {
         double u = VV[c][0];
         double v = VV[c][1];
         double n = VV[c][2];
-        
+
         A(c,0) = u*u;
         A(c,1) = u*v;
         A(c,2) = v*v;
         A(c,3) = u;
         A(c,4) = v;
-        
+
         b(c) = n;
       }
-      
+
       sol=A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
 
       return Quadric(sol(0),sol(1),sol(2),sol(3),sol(4));
     }
   };
-  
+
 public:
-  
+
   Eigen::MatrixXd vertices;
   // Face list of current mesh    (#F x 3) or (#F x 4)
   // The i-th row contains the indices of the vertices that forms the i-th face in ccw order
   Eigen::MatrixXi faces;
-  
+
   std::vector<std::vector<int> > vertex_to_vertices;
   std::vector<std::vector<int> > vertex_to_faces;
   std::vector<std::vector<int> > vertex_to_faces_index;
   Eigen::MatrixXd face_normals;
   Eigen::MatrixXd vertex_normals;
-  
+
   /* Size of the neighborhood */
   double sphereRadius;
   int kRing;
-  
+
   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 */
   normalType nt;
-  
+
   double lastRadius;
   double scaledRadius;
   std::string lastMeshName;
-  
+
   /* Benchmark related variables */
   bool expStep; /* True if we want the radius to increase exponentially */
   int step;  /* If expStep==false, by how much rhe radius increases on every step */
   int maxSize; /* The maximum limit of the radius in the benchmark */
-  
+
   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>&);
@@ -192,48 +192,48 @@ public:
   IGL_INLINE void computeCurvature();
   IGL_INLINE void printCurvature(std::string outpath);
   IGL_INLINE double getAverageEdge();
-  
+
   IGL_INLINE static int rotateForward (float *v0, float *v1, float *v2)
   {
     float t;
-    
+
     if (abs(*v2) >= abs(*v1) && abs(*v2) >= abs(*v0))
       return 0;
-    
+
     t = *v0;
     *v0 = *v2;
     *v2 = *v1;
     *v1 = t;
-    
+
     return 1 + rotateForward (v0, v1, v2);
   }
-  
+
   IGL_INLINE static void rotateBackward (int nr, float *v0, float *v1, float *v2)
   {
     float t;
-    
+
     if (nr == 0)
       return;
-    
+
     t = *v2;
     *v2 = *v0;
     *v0 = *v1;
     *v1 = t;
-    
+
     rotateBackward (nr - 1, v0, v1, v2);
   }
-  
+
   IGL_INLINE static Eigen::Vector3d chooseMax (Eigen::Vector3d n, Eigen::Vector3d abc, float ab)
   {
     int i, max_i;
     float max_sp;
     Eigen::Vector3d nt[8];
-    
+
     n.normalize ();
     abc.normalize ();
-    
+
     max_sp = - std::numeric_limits<float>::max();
-    
+
     for (i = 0; i < 4; i++)
     {
       nt[i] = n;
@@ -243,16 +243,16 @@ public:
         {
           case 0:
             break;
-            
+
           case 1:
             nt[i][2] = -n[2];
             break;
-            
+
           case 2:
             nt[i][0] = -n[0];
             nt[i][1] = -n[1];
             break;
-            
+
           case 3:
             nt[i][0] = -n[0];
             nt[i][1] = -n[1];
@@ -267,33 +267,33 @@ public:
           case 0:
             nt[i][0] = -n[0];
             break;
-            
+
           case 1:
             nt[i][1] = -n[1];
             break;
-            
+
           case 2:
             nt[i][0] = -n[0];
             nt[i][2] = -n[2];
             break;
-            
+
           case 3:
             nt[i][1] = -n[1];
             nt[i][2] = -n[2];
             break;
         }
       }
-      
+
       if (nt[i].dot(abc) > max_sp)
       {
         max_sp = nt[i].dot(abc);
         max_i = i;
       }
     }
-    
+
     return nt[max_i];
   }
-  
+
 };
 
 class comparer
@@ -325,11 +325,11 @@ IGL_INLINE void CurvatureCalculator::init(const Eigen::MatrixXd& V, const Eigen:
 {
   // Normalize vertices
   vertices = V;
-  
+
 //  vertices = vertices.array() - vertices.minCoeff();
 //  vertices = vertices.array() / vertices.maxCoeff();
 //  vertices = vertices.array() * (1.0/igl::avg_edge_length(V,F));
-  
+
   faces = F;
   igl::adjacency_list(F, vertex_to_vertices);
   igl::vf(V, F, vertex_to_faces, vertex_to_faces_index);
@@ -341,14 +341,14 @@ IGL_INLINE void CurvatureCalculator::fitQuadric (Eigen::Vector3d v, std::vector<
 {
   std::vector<Eigen::Vector3d> points;
   points.reserve (vv.size());
-  
+
   for (unsigned int i = 0; i < vv.size(); ++i) {
-    
+
     Eigen::Vector3d  cp = vertices.row(vv[i]);
-    
+
     // vtang non e` il v tangente!!!
     Eigen::Vector3d  vTang = cp - v;
-    
+
     double x = vTang.dot(ref[0]);
     double y = vTang.dot(ref[1]);
     double z = vTang.dot(ref[2]);
@@ -359,13 +359,13 @@ IGL_INLINE void CurvatureCalculator::fitQuadric (Eigen::Vector3d v, std::vector<
 
 IGL_INLINE void CurvatureCalculator::finalEigenStuff (int i, 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();
-  
+
 //  if (fabs(a) < 10e-8 || fabs(b) < 10e-8)
 //  {
 //    std::cout << "Degenerate quadric: " << i << std::endl;
@@ -374,52 +374,52 @@ IGL_INLINE void CurvatureCalculator::finalEigenStuff (int i, std::vector<Eigen::
   double E = 1.0 + d*d;
   double F = d*e;
   double G = 1.0 + e*e;
-  
+
   Eigen::Vector3d n = Eigen::Vector3d(-d,-e,1.0).normalized();
-  
+
   double L = 2.0 * a * n[2];
   double M = b * n[2];
   double N = 2 * c * n[2];
-  
-  
+
+
   // ----------------- Eigen stuff
   Eigen::Matrix2d m;
   m << L*G - M*F, M*E-L*F, M*E-L*F, N*E-M*F;
   m = m / (E*G-F*F);
   Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> eig(m);
-  
+
   Eigen::Vector2d c_val = eig.eigenvalues();
   Eigen::Matrix2d c_vec = eig.eigenvectors();
 
   // std::cerr << "c_val:" << c_val << std::endl;
   // std::cerr << "c_vec:" << c_vec << std::endl;
-  
+
   // std::cerr << "c_vec:" << c_vec(0) << " "  << c_vec(1) << std::endl;
 
   c_val = -c_val;
-  
+
   Eigen::Vector3d v1, v2;
   v1[0] = c_vec(0);
   v1[1] = c_vec(1);
   v1[2] = 0; //d * v1[0] + e * v1[1];
-  
+
   v2[0] = c_vec(2);
   v2[1] = c_vec(3);
   v2[2] = 0; //d * v2[0] + e * v2[1];
-  
-  
+
+
   // v1 = v1.normalized();
   // v2 = v2.normalized();
-  
+
   Eigen::Vector3d v1global = ref[0] * v1[0] + ref[1] * v1[1] + ref[2] * v1[2];
   Eigen::Vector3d v2global = ref[0] * v2[0] + ref[1] * v2[1] + ref[2] * v2[2];
-  
+
   v1global.normalize();
   v2global.normalize();
-  
+
   v1global *= c_val(0);
   v2global *= c_val(1);
-  
+
   if (c_val[0] > c_val[1])
   {
     curv[i]=std::vector<double>(2);
@@ -532,12 +532,12 @@ IGL_INLINE Eigen::Vector3d CurvatureCalculator::project(Eigen::Vector3d v, Eigen
 
 IGL_INLINE void CurvatureCalculator::computeReferenceFrame(int i, 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]));
-  
+
   longest_v=(project(vertices.row(i),longest_v,normal)-Eigen::Vector3d(vertices.row(i))).normalized();
-  
+
   /* L'ultimo asse si ottiene come prodotto vettoriale tra i due
    * calcolati */
   Eigen::Vector3d y_axis=(normal.cross(longest_v)).normalized();
@@ -551,7 +551,7 @@ IGL_INLINE void CurvatureCalculator::getAverageNormal(int j, std::vector<int> vv
   normal=(vertex_normals.row(j)).normalized();
   if (localMode)
     return;
-  
+
   for (unsigned int i=0; i<vv.size(); i++)
   {
     normal+=vertex_normals.row(vv[i]).normalized();
@@ -565,9 +565,9 @@ IGL_INLINE void CurvatureCalculator::getProjPlane(int j, std::vector<int> vv, Ei
   float a, b, c;
   float nx, ny, nz;
   float abcq;
-  
+
   a = b = c = 0;
-  
+
   if (localMode)
   {
     for (unsigned int i=0; i<vertex_to_faces.at(j).size(); ++i)
@@ -594,7 +594,7 @@ IGL_INLINE void CurvatureCalculator::getProjPlane(int j, std::vector<int> vv, Ei
   nz = sqrt (1 - nx*nx - ny*ny);
   rotateBackward (nr, &a, &b, &c);
   rotateBackward (nr, &nx, &ny, &nz);
-  
+
   ppn = chooseMax (Eigen::Vector3d(nx, ny, nz), Eigen::Vector3d (a, b, c), a * b);
   ppn.normalize();
 }
@@ -604,21 +604,21 @@ IGL_INLINE double CurvatureCalculator::getAverageEdge()
 {
   double sum = 0;
   int count = 0;
-  
+
   for (int i = 0; i<faces.rows(); i++)
   {
     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]);
-      
+
       double l = (p1-p2).norm();
-      
+
       sum+=l;
       ++count;
     }
   }
-  
+
   return (sum/(double)count);
 }
 
@@ -637,7 +637,7 @@ IGL_INLINE void CurvatureCalculator::applyMontecarlo(std::vector<int>& vin, std:
     *vout = vin;
     return;
   }
-  
+
   float p = ((float) montecarloN) / (float) vin.size();
   for (std::vector<int>::iterator vpi = vin.begin(); vpi != vin.end(); ++vpi)
   {
@@ -652,27 +652,27 @@ IGL_INLINE void CurvatureCalculator::applyMontecarlo(std::vector<int>& vin, std:
 IGL_INLINE void CurvatureCalculator::computeCurvature()
 {
   using namespace std;
-  
+
   //CHECK che esista la mesh
   size_t vertices_count=vertices.rows() ;
-  
+
   if (vertices_count <=0)
     return;
-  
+
   curvDir=std::vector< std::vector<Eigen::Vector3d> >(vertices_count);
   curv=std::vector<std::vector<double> >(vertices_count);
-  
-  
-  
+
+
+
   scaledRadius=getAverageEdge()*sphereRadius;
-  
+
   std::vector<int> vv;
   std::vector<int> vvtmp;
   Eigen::Vector3d normal;
-  
+
   //double time_spent;
   //double searchtime=0, ref_time=0, fit_time=0, final_time=0;
-  
+
   for (size_t i=0; i<vertices_count; ++i)
   {
     vv.clear();
@@ -690,7 +690,7 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
         fprintf(stderr,"Error: search type not recognized");
         return;
     }
-    
+
     std::vector<Eigen::Vector3d> ref(3);
     if (vv.size()<6)
     {
@@ -709,7 +709,7 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
         fprintf(stderr,"Error: normal type not recognized");
         return;
     }
-    
+
     if (projectionPlaneCheck)
     {
       vvtmp.reserve (vv.size ());
@@ -717,7 +717,7 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
       if (vvtmp.size() >= 6)
         vv = vvtmp;
     }
-    
+
     if (vv.size()<6)
     {
       std::cerr << "Could not compute curvature of radius " << scaledRadius << endl;
@@ -731,16 +731,16 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
       applyMontecarlo(vv,&vvtmp);
       vv=vvtmp;
     }
-    
+
     if (vv.size()<6)
       return;
     computeReferenceFrame(i,normal,ref);
-    
+
     Quadric q;
     fitQuadric (me, ref, vv, &q);
     finalEigenStuff(i,ref,q);
   }
-  
+
   lastRadius=sphereRadius;
   curvatureComputed=true;
 }
@@ -750,16 +750,16 @@ IGL_INLINE void CurvatureCalculator::printCurvature(std::string outpath)
   using namespace std;
   if (!curvatureComputed)
     return;
-  
+
   std::ofstream of;
   of.open(outpath.c_str());
-  
+
   if (!of)
   {
     fprintf(stderr, "Error: could not open output file %s\n", outpath.c_str());
     return;
   }
-  
+
   int vertices_count=vertices.rows();
   of << vertices_count << endl;
   for (int i=0; i<vertices_count; i++)
@@ -767,9 +767,9 @@ IGL_INLINE void CurvatureCalculator::printCurvature(std::string outpath)
     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;
   }
-  
+
   of.close();
-  
+
 }
 
 template <typename DerivedV, typename DerivedF>
@@ -778,51 +778,61 @@ IGL_INLINE void igl::principal_curvature(
                                          const Eigen::PlainObjectBase<DerivedF>& F,
                                          Eigen::PlainObjectBase<DerivedV>& PD1,
                                          Eigen::PlainObjectBase<DerivedV>& PD2,
+                                         Eigen::PlainObjectBase<DerivedV>& PV1,
+                                         Eigen::PlainObjectBase<DerivedV>& PV2,
                                          unsigned radius,
                                          bool useKring
                                          )
 {
   using namespace std;
-  
+
   // 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++)
   {
     Eigen::Vector3d d1;
     Eigen::Vector3d d2;
-    d1 << cc.curvDir[i][0][0], cc.curvDir[i][0][1], cc.curvDir[i][0][2];
-    d2 << cc.curvDir[i][1][0], cc.curvDir[i][1][1], cc.curvDir[i][1][2];
-    d1.normalize();
-    d2.normalize();
-    d1 *= cc.curv[i][0];
-    d2 *= cc.curv[i][1];
-    PD1.row(i) = d1;
-    PD2.row(i) = 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();
+    PD2.row(i).normalize();
     
+    if (isnan(PD1(i,0)) || isnan(PD1(i,1)) || isnan(PD1(i,2)) || isnan(PD2(i,0)) || isnan(PD2(i,1)) || 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)
     {
-      cerr << "Something is wrong, vertex: i" << endl;
+      cerr << "PRINCIPAL_CURVATURE: Something is wrong with vertex: i" << endl;
       PD1.row(i) *= 0;
       PD2.row(i) *= 0;
     }
   }
-  
-}
 
+}

include/igl/principal_curvature.h

@@ -1,9 +1,9 @@
 // This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+//
+// Copyright (C) 2013 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
 // obtain one at http://mozilla.org/MPL/2.0/.
 #ifndef IGL_PRINCIPAL_CURVATURE_H
 #define IGL_PRINCIPAL_CURVATURE_H
@@ -34,8 +34,9 @@ namespace igl
   //
   // Outputs:
   //   PD1 #V by 3 maximal curvature direction for each vertex.
-  //   PD2 #V by 3 minimal curvature direction for each vertex.  
-  //  Note that to maximal/minimal curvature value is the rowise norm of PD1/PD2
+  //   PD2 #V by 3 minimal curvature direction for each vertex.
+  //   PV1 #V by 1 maximal curvature value for each vertex.
+  //   PV2 #V by 1 minimal curvature value for each vertex.
   //
   // See also: moveVF, moveFV
   //
@@ -51,6 +52,8 @@ IGL_INLINE void principal_curvature(
   const Eigen::PlainObjectBase<DerivedF>& F,
   Eigen::PlainObjectBase<DerivedV>& PD1,
   Eigen::PlainObjectBase<DerivedV>& PD2,
+  Eigen::PlainObjectBase<DerivedV>& PV1,
+  Eigen::PlainObjectBase<DerivedV>& PV2,
   unsigned radius = 5,
   bool useKring = true);
 }

include/igl/viewer/Viewer.cpp

@@ -1616,7 +1616,7 @@ namespace igl
 
   void Viewer::OpenGL_state::draw_overlay_lines()
   {
-    glDrawElements(GL_LINES, 2*lines_F_vbo.cols(), GL_UNSIGNED_INT, 0);
+    glDrawElements(GL_LINES, lines_F_vbo.cols(), GL_UNSIGNED_INT, 0);
   }
 
   void Viewer::OpenGL_state::draw_overlay_points()
@@ -2151,10 +2151,22 @@ namespace igl
   void Viewer::set_mesh(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F)
   {
     using namespace std;
+
+    Eigen::MatrixXd V_temp;
+    
+    // If V only has two columns, pad with a column of zeros
+    if (V.cols() == 2)
+    {
+      V_temp = Eigen::MatrixXd::Zero(V.rows(),3);
+      V_temp.block(0,0,V.rows(),2) = V;
+    }
+    else
+      V_temp = V;
+
     if (data.V.rows() == 0 && data.F.rows() == 0)
     {
       clear_mesh();
-      data.V = V;
+      data.V = V_temp;
       data.F = F;
 
       compute_normals();
@@ -2169,7 +2181,7 @@ namespace igl
     {
       if (data.V.rows() == V.rows() && data.F.rows() == F.rows())
       {
-        data.V = V;
+        data.V = V_temp;
         data.F = F;
         alignCameraCenter();
       }
@@ -2249,7 +2261,7 @@ namespace igl
 
   }
 
-  void Viewer::set_UV(const Eigen::MatrixXd& UV)
+  void Viewer::set_uv(const Eigen::MatrixXd& UV)
   {
     using namespace std;
     if (UV.rows() == data.V.rows())
@@ -2262,7 +2274,7 @@ namespace igl
     data.dirty |= DIRTY_UV;
   }
 
-  void Viewer::set_UV(const Eigen::MatrixXd& UV_V, const Eigen::MatrixXi& UV_F)
+  void Viewer::set_uv(const Eigen::MatrixXd& UV_V, const Eigen::MatrixXi& UV_F)
   {
     set_face_based(true);
     data.V_uv = UV_V;
@@ -2282,23 +2294,27 @@ namespace igl
     data.dirty |= DIRTY_TEXTURE;
   }
 
-  void Viewer::draw_points(const Eigen::MatrixXd& P,  const Eigen::MatrixXd& C)
+  void Viewer::add_points(const Eigen::MatrixXd& P,  const Eigen::MatrixXd& C)
   {
     int lastid = data.points.rows();
     data.points.conservativeResize(data.points.rows() + P.rows(),6);
     for (unsigned i=0; i<P.rows(); ++i)
       data.points.row(lastid+i) << P.row(i), i<C.rows() ? C.row(i) : C.row(C.rows()-1);
+
+    data.dirty |= DIRTY_OVERLAY_POINTS;
   }
 
-  void Viewer::draw_edges (const Eigen::MatrixXd& P1, const Eigen::MatrixXd& P2, const Eigen::MatrixXd& C)
+  void Viewer::add_edges(const Eigen::MatrixXd& P1, const Eigen::MatrixXd& P2, const Eigen::MatrixXd& C)
   {
     int lastid = data.lines.rows();
     data.lines.conservativeResize(data.lines.rows() + P1.rows(),9);
     for (unsigned i=0; i<P1.rows(); ++i)
       data.lines.row(lastid+i) << P1.row(i), P2.row(i), i<C.rows() ? C.row(i) : C.row(C.rows()-1);
+
+    data.dirty |= DIRTY_OVERLAY_LINES;
   }
 
-  void Viewer::draw_label (const Eigen::VectorXd& P,  const std::string& str)
+  void Viewer::add_label(const Eigen::VectorXd& P,  const std::string& str)
   {
     int lastid = data.labels_positions.rows();
     data.labels_positions.conservativeResize(lastid+1, 3);

include/igl/viewer/Viewer.h

@@ -357,16 +357,16 @@ namespace igl
     // Inputs:
     //   C  #V|#F|1 by 3 list of colors
     void set_colors(const Eigen::MatrixXd &C);
-    void set_UV(const Eigen::MatrixXd& UV);
-    void set_UV(const Eigen::MatrixXd& UV_V, const Eigen::MatrixXi& UV_F);
+    void set_uv(const Eigen::MatrixXd& UV);
+    void set_uv(const Eigen::MatrixXd& UV_V, const Eigen::MatrixXi& UV_F);
     void set_texture(
                       const Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic>& R,
                       const Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic>& G,
                       const Eigen::Matrix<char,Eigen::Dynamic,Eigen::Dynamic>& B);
 
-    void draw_points(const Eigen::MatrixXd& P,  const Eigen::MatrixXd& C);
-    void draw_edges (const Eigen::MatrixXd& P1, const Eigen::MatrixXd& P2, const Eigen::MatrixXd& C);
-    void draw_label (const Eigen::VectorXd& P,  const std::string& str);
+    void add_points(const Eigen::MatrixXd& P,  const Eigen::MatrixXd& C);
+    void add_edges (const Eigen::MatrixXd& P1, const Eigen::MatrixXd& P2, const Eigen::MatrixXd& C);
+    void add_label (const Eigen::VectorXd& P,  const std::string& str);
 
     // Save the OpenGL transformation matrices used for the previous rendering pass
     Eigen::Matrix4f view;

tutorial/104_Events/main.cpp

@@ -15,12 +15,12 @@ bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
     // Draw_mesh creates or updates the vertices and faces of the displayed mesh.
     // If a mesh is already displayed, draw_mesh returns an error if the given V and
     // F have size different than the current ones
-    viewer.draw_mesh(V1, F1);
+    viewer.set_mesh(V1, F1);
   }
   else if (key == '2')
   {
     viewer.clear_mesh();
-    viewer.draw_mesh(V2, F2);
+    viewer.set_mesh(V2, F2);
   }
 
   return false;
@@ -37,6 +37,6 @@ int main(int argc, char *argv[])
   // Register a keyboard callback that allows to switch between
   // the two loaded meshes
   viewer.callback_key_down = &key_down;
-  viewer.draw_mesh(V1, F1);
+  viewer.set_mesh(V1, F1);
   viewer.launch();
 }

tutorial/106_Overlays/main.cpp

@@ -45,14 +45,14 @@ int main(int argc, char *argv[])
 
   // Plot the mesh
   igl::Viewer viewer;
-  viewer.draw_mesh(V, F);
+  viewer.set_mesh(V, F);
 
   // Plot the corners of the bounding box as points
-  viewer.draw_points(V_box,Eigen::RowVector3d(1,0,0));
+  viewer.add_points(V_box,Eigen::RowVector3d(1,0,0));
 
   // Plot the edges of the bounding box
   for (unsigned i=0;i<E_box.rows(); ++i)
-    viewer.draw_edges
+    viewer.add_edges
     (
       V_box.row(E_box(i,0)),
       V_box.row(E_box(i,1)),
@@ -65,10 +65,10 @@ int main(int argc, char *argv[])
   // Plot labels with the coordinates of bounding box vertices
   std::stringstream l1;
   l1 << m(0) << ", " << m(1) << ", " << m(2);
-  viewer.draw_label(m,l1.str());
+  viewer.add_label(m,l1.str());
   std::stringstream l2;
   l2 << M(0) << ", " << M(1) << ", " << M(2);
-  viewer.draw_label(M,l2.str());
+  viewer.add_label(M,l2.str());
 
   // Launch the viewer
   viewer.launch();

tutorial/203_CurvatureDirections/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(203_CurvatureDirections)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${CMAKE_PROJECT_NAME} ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${CMAKE_PROJECT_NAME} ${SHARED_LIBRARIES})

tutorial/203_CurvatureDirections/main.cpp

@@ -0,0 +1,37 @@
+#include <igl/readOFF.h>
+#include <igl/viewer/Viewer.h>
+#include <igl/per_vertex_normals.h>
+#include <igl/per_face_normals.h>
+#include <igl/per_corner_normals.h>
+#include <igl/avg_edge_length.h>
+#include <igl/barycenter.h>
+#include <igl/principal_curvature.h>
+
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+
+int main(int argc, char *argv[])
+{
+  // Load a mesh in OFF format
+  igl::readOFF("../shared/fertility.off", V, F);
+
+  // Compute curvature directions via quadric fitting
+  Eigen::MatrixXd PD1,PD2,PV1,PV2;
+  igl::principal_curvature(V,F,PD1,PD2,PV1,PV2);
+
+  // Plot the mesh
+  igl::Viewer viewer;
+  viewer.set_mesh(V, F);
+
+  // Find the average edge length
+  double avg = igl::avg_edge_length(V,F);
+  
+  // Draw a red segment on each vertex parallel to the minimal curvature direction
+  viewer.add_edges(V + PD1*avg, V - PD1*avg, Eigen::RowVector3d(0,0,1));
+
+  // Draw a blue segment on each vertex parallel to the maximal curvature direction
+  viewer.add_edges(V + PD2*avg, V - PD2*avg, Eigen::RowVector3d(1,0,0));
+
+  // Launch the viewer
+  viewer.launch();
+}

tutorial/501_HarmonicParam/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(501_HarmonicParam)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${CMAKE_PROJECT_NAME} ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${CMAKE_PROJECT_NAME} ${SHARED_LIBRARIES})

tutorial/501_HarmonicParam/main.cpp

@@ -0,0 +1,58 @@
+#include <igl/readOFF.h>
+#include <igl/viewer/Viewer.h>
+#include <igl/boundary_vertices_sorted.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/harmonic.h>
+
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+Eigen::MatrixXd V_uv;
+
+bool key_down(igl::Viewer& viewer, unsigned char key, int modifier)
+{
+  if (key == '1')
+    // Plot the 3D mesh
+    viewer.set_mesh(V,F);
+  else if (key == '2')
+    // Plot the mesh in 2D using the UV coordinates as vertex coordinates
+    viewer.set_mesh(V_uv,F);
+
+  viewer.compute_normals();
+
+  return false;
+}
+
+int main(int argc, char *argv[])
+{
+  // Load a mesh in OFF format
+  igl::readOFF("../shared/camelhead.off", V, F);
+
+  // Find the open boundary
+  Eigen::VectorXi bnd;
+  igl::boundary_vertices_sorted(V,F,bnd);
+
+  // Map the boundary to a circle, preserving edge proportions
+  Eigen::MatrixXd bnd_uv;
+  igl::map_vertices_to_circle(V,F,bnd,bnd_uv);
+
+  // Harmonic parametrization for the internal vertices
+  igl::harmonic(V,F,bnd,bnd_uv,1,V_uv);
+
+  // Scale UV to make the texture more clear
+  V_uv *= 5;
+
+  // Plot the mesh
+  igl::Viewer viewer;
+  viewer.set_mesh(V, F);
+  viewer.set_uv(V_uv);
+  viewer.callback_key_down = &key_down;
+
+  // Disable wireframe
+  viewer.options.show_lines = false;
+
+  // Draw checkerboard texture
+  viewer.options.show_texture = true;
+
+  // Launch the viewer
+  viewer.launch();
+}

tutorial/CMakeLists.shared

@@ -9,12 +9,15 @@ set(CMAKE_COLOR_MAKEFILE ON)
 
 find_package(ANTTWEAKBAR REQUIRED)
 find_package(OpenGL REQUIRED)
-find_package(GLEW REQUIRED)
 find_package(EIGEN REQUIRED)
 find_package(LIBIGL REQUIRED)
 find_package(TINYXML2 REQUIRED)
 find_package(GLFW REQUIRED)
 
+if(WIN32)
+	find_package(GLEW REQUIRED)
+endif(WIN32)
+
 if(APPLE)
 	set(CMAKE_CXX_LINK_FLAGS "-framework OpenGL -framework Cocoa -framework IOKit -framework CoreVideo")
 endif (APPLE) #APPLE
@@ -23,8 +26,10 @@ include_directories( ${ANT_TWEAK_BAR_INCLUDE_DIR} )
 
 # message(FATAL_ERROR ${ANT_TWEAK_BAR_INCLUDE_DIR})
 
+if(WIN32)
+	include_directories( ${GLEW_INCLUDE_DIR} )
+endif(WIN32)
 
-include_directories( ${GLEW_INCLUDE_DIR} )
 include_directories( ${EIGEN_INCLUDE_DIR})
 include_directories( ${LIBIGL_INCLUDE_DIR})
 include_directories( ${TINYXML2_INCLUDE_DIR})
@@ -50,7 +55,11 @@ if ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang" OR "${CMAKE_CXX_COMPILER_ID}" ST
 	set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-deprecated-register")
 endif()
 
-set(SHARED_SOURCES ${LIBIGL_SOURCES} ${GLEW_SOURCES})
+set(SHARED_SOURCES ${LIBIGL_SOURCES})
+
+if(WIN32)
+	set(SHARED_SOURCES ${SHARED_SOURCES} ${GLEW_SOURCES})
+endif(WIN32)
 
 # message(FATAL_ERROR ${LIBIGL_SOURCES})
 

tutorial/shared/camelhead.off.REMOVED.git-id

@@ -0,0 +1 @@
+9bfd194b1bc4dd4b10e394f41b05e91643a64779

tutorial/shared/lion.off.REMOVED.git-id

@@ -0,0 +1 @@
+6aa7cd79a038749b98f9f2602e4829037295bcc7