cleanup new code mess enough to compile

Former-commit-id: d35a417d6c65061032c723eca386035b510f65cb

include/igl/adjacency_matrix.cpp

@@ -12,24 +12,30 @@
 // Bug in unsupported/Eigen/SparseExtra needs iostream first
 #include <iostream>
 #include <unsupported/Eigen/SparseExtra>
+#warning "old includes"
+#include <vector>
 
 template <typename T>
 IGL_INLINE void igl::adjacency_matrix(
   const Eigen::MatrixXi & F, 
   Eigen::SparseMatrix<T>& A)
 {
+  using namespace std;
   Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> 
     dyn_A(F.maxCoeff()+1, F.maxCoeff()+1);
   switch(F.cols())
   {
     case 3:
-      dyn_A.reserve(6*(F.maxCoeff()+1));
+      dyn_A.reserve(8*(F.maxCoeff()+1));
       break;
     case 4:
       dyn_A.reserve(26*(F.maxCoeff()+1));
       break;
   }
 
+  //typedef Triplet<T> IJV;
+  //vector<IJV > ijv;
+  //ijv.reserve(F.size()*2);
   // Loop over faces
   for(int i = 0;i<F.rows();i++)
   {

include/igl/cross_field_missmatch.cpp

@@ -9,143 +9,143 @@
 #include "tt.h"
 
 
-namespace igl {
-  template <typename DerivedV, typename DerivedF, typename DerivedO>
-  class MissMatchCalculator
-  {
-  public:
-    
-    const Eigen::PlainObjectBase<DerivedV> &V;
-    const Eigen::PlainObjectBase<DerivedF> &F;
-    const Eigen::PlainObjectBase<DerivedV> &PD1;
-    const Eigen::PlainObjectBase<DerivedV> &PD2;
-    Eigen::PlainObjectBase<DerivedV> N;
-    
-  private:
-    // internal
-    std::vector<bool> V_border; // bool
-    std::vector<std::vector<int> > VF;
-    std::vector<std::vector<int> > VFi;
-    Eigen::PlainObjectBase<DerivedF> TT;
-    Eigen::PlainObjectBase<DerivedF> TTi;
-    
-    
-  private:
-    
-    ///return true if a vertex is singluar by looking at initialized missmatches
-    // possible bugs, verify deleted flag vs IsD()
-    // not sorted vf, but should not make a difference
-    // olga: TODO: this returns the index modulo 4.
-    int oneRingMissMatch(const int vid)
-    {
-      ///check that is on border..
-      if (V_border[vid])
-        return 0;
-      
-      int missmatch=0;
-      for (unsigned int i=0;i<VF[vid].size();i++)
-      {
-        // look for the vertex
-        int j=-1;
-        for (unsigned z=0; z<3; ++z)
-          if (F(VF[vid][i],z) == vid)
-            j=z;
-        assert(j!=-1);
-        
-        missmatch+=Handle_MMatch(VF[vid][i],j);
-      }
-      
-      missmatch=missmatch%4;
-      return missmatch;
-    }
-    
-    
-    ///compute the mismatch between 2 faces
-    int MissMatchByCross(const int f0,
-                         const int f1)
-    {
-      Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0 = PD1.row(f0);
-      Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1 = PD1.row(f1);
-      Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n0 = N.row(f0);
-      Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n1 = N.row(f1);
-      
-      Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1Rot = Comb<DerivedV, DerivedF>::Rotate(n1,n0,dir1);
-      dir1Rot.normalize();
-      
-      // TODO: this should be equivalent to the other code below, to check!
-      // Compute the angle between the two vectors
-      //    double a0 = atan2(dir0.dot(B2.row(f0)),dir0.dot(B1.row(f0)));
-      //    double a1 = atan2(dir1Rot.dot(B2.row(f0)),dir1Rot.dot(B1.row(f0)));
-      //
-      //    double angle_diff = a1-a0;   //VectToAngle(f0,dir1Rot);
-      
-      double angle_diff = atan2(dir1Rot.dot(PD2.row(f0)),dir1Rot.dot(PD1.row(f0)));
-      
-      //    std::cerr << "Dani: " << dir0(0) << " " << dir0(1) << " " << dir0(2) << " " << dir1Rot(0) << " " << dir1Rot(1) << " " << dir1Rot(2) << " " << angle_diff << std::endl;
-      
-      double step=M_PI/2.0;
-      int i=(int)floor((angle_diff/step)+0.5);
-      int k=0;
-      if (i>=0)
-        k=i%4;
-      else
-        k=(-(3*i))%4;
-      return k;
-    }
-    
-    
-public:
-  MissMatchCalculator(const Eigen::PlainObjectBase<DerivedV> &_V,
-                      const Eigen::PlainObjectBase<DerivedF> &_F,
-                      const Eigen::PlainObjectBase<DerivedV> &_PD1,
-                      const Eigen::PlainObjectBase<DerivedV> &_PD2
-                      ):
-  V(_V),
-  F(_F),
-  PD1(_PD1),
-  PD2(_PD2)
-  {
-    igl::per_face_normals(V,F,N);
-    V_border = igl::is_border_vertex(V,F);
-    igl::vf(V,F,VF,VFi);
-    igl::tt(V,F,TT,TTi);
-  }
-  
-  void calculateMissmatch(Eigen::PlainObjectBase<DerivedO> &Handle_MMatch)
-  {
-    Handle_MMatch.setConstant(F.rows(),3,-1);
-    for (unsigned int i=0;i<F.rows();i++)
-    {
-      for (int j=0;j<3;j++)
-      {
-        if (i==TT(i,j) || TT(i,j) == -1)
-          Handle_MMatch(i,j)=0;
-        else
-          Handle_MMatch(i,j) = MissMatchByCross(i,TT(i,j));
-      }
-    }
-  }
-  
-};
-}
-template <typename DerivedV, typename DerivedF, typename DerivedO>
-IGL_INLINE void igl::cross_field_missmatch(const Eigen::PlainObjectBase<DerivedV> &V,
-                                           const Eigen::PlainObjectBase<DerivedF> &F,
-                                           const Eigen::PlainObjectBase<DerivedV> &PD1,
-                                           const Eigen::PlainObjectBase<DerivedV> &PD2,
-                                           const bool isCombed,
-                                           Eigen::PlainObjectBase<DerivedO> &missmatch)
-{
-  Eigen::PlainObjectBase<DerivedV> PD1_combed;
-  Eigen::PlainObjectBase<DerivedV> PD2_combed;
-  
-  if (!isCombed)
-    igl::comb_cross_field(V,F,PD1,PD2,PD1_combed,PD2_combed);
-  else
-  {
-    PD1_combed = PD1;
-    PD2_combed = PD2;
-  }
-  igl::MissMatchCalculator<DerivedV, DerivedF, DerivedO> sf(V, F, PD1_combed, PD2_combed);
-  sf.calculateMissmatch(missmatch);
-}
+//namespace igl {
+//  template <typename DerivedV, typename DerivedF, typename DerivedO>
+//  class MissMatchCalculator
+//  {
+//  public:
+//    
+//    const Eigen::PlainObjectBase<DerivedV> &V;
+//    const Eigen::PlainObjectBase<DerivedF> &F;
+//    const Eigen::PlainObjectBase<DerivedV> &PD1;
+//    const Eigen::PlainObjectBase<DerivedV> &PD2;
+//    Eigen::PlainObjectBase<DerivedV> N;
+//    
+//  private:
+//    // internal
+//    std::vector<bool> V_border; // bool
+//    std::vector<std::vector<int> > VF;
+//    std::vector<std::vector<int> > VFi;
+//    Eigen::PlainObjectBase<DerivedF> TT;
+//    Eigen::PlainObjectBase<DerivedF> TTi;
+//    
+//    
+//  private:
+//    
+//    ///return true if a vertex is singluar by looking at initialized missmatches
+//    // possible bugs, verify deleted flag vs IsD()
+//    // not sorted vf, but should not make a difference
+//    // olga: TODO: this returns the index modulo 4.
+//    int oneRingMissMatch(const int vid)
+//    {
+//      ///check that is on border..
+//      if (V_border[vid])
+//        return 0;
+//      
+//      int missmatch=0;
+//      for (unsigned int i=0;i<VF[vid].size();i++)
+//      {
+//        // look for the vertex
+//        int j=-1;
+//        for (unsigned z=0; z<3; ++z)
+//          if (F(VF[vid][i],z) == vid)
+//            j=z;
+//        assert(j!=-1);
+//        
+//        missmatch+=Handle_MMatch(VF[vid][i],j);
+//      }
+//      
+//      missmatch=missmatch%4;
+//      return missmatch;
+//    }
+//    
+//    
+//    ///compute the mismatch between 2 faces
+//    int MissMatchByCross(const int f0,
+//                         const int f1)
+//    {
+//      Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir0 = PD1.row(f0);
+//      Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1 = PD1.row(f1);
+//      Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n0 = N.row(f0);
+//      Eigen::Matrix<typename DerivedV::Scalar, 3, 1> n1 = N.row(f1);
+//      
+//      Eigen::Matrix<typename DerivedV::Scalar, 3, 1> dir1Rot = Comb<DerivedV, DerivedF>::Rotate(n1,n0,dir1);
+//      dir1Rot.normalize();
+//      
+//      // TODO: this should be equivalent to the other code below, to check!
+//      // Compute the angle between the two vectors
+//      //    double a0 = atan2(dir0.dot(B2.row(f0)),dir0.dot(B1.row(f0)));
+//      //    double a1 = atan2(dir1Rot.dot(B2.row(f0)),dir1Rot.dot(B1.row(f0)));
+//      //
+//      //    double angle_diff = a1-a0;   //VectToAngle(f0,dir1Rot);
+//      
+//      double angle_diff = atan2(dir1Rot.dot(PD2.row(f0)),dir1Rot.dot(PD1.row(f0)));
+//      
+//      //    std::cerr << "Dani: " << dir0(0) << " " << dir0(1) << " " << dir0(2) << " " << dir1Rot(0) << " " << dir1Rot(1) << " " << dir1Rot(2) << " " << angle_diff << std::endl;
+//      
+//      double step=M_PI/2.0;
+//      int i=(int)floor((angle_diff/step)+0.5);
+//      int k=0;
+//      if (i>=0)
+//        k=i%4;
+//      else
+//        k=(-(3*i))%4;
+//      return k;
+//    }
+//    
+//    
+//public:
+//  MissMatchCalculator(const Eigen::PlainObjectBase<DerivedV> &_V,
+//                      const Eigen::PlainObjectBase<DerivedF> &_F,
+//                      const Eigen::PlainObjectBase<DerivedV> &_PD1,
+//                      const Eigen::PlainObjectBase<DerivedV> &_PD2
+//                      ):
+//  V(_V),
+//  F(_F),
+//  PD1(_PD1),
+//  PD2(_PD2)
+//  {
+//    igl::per_face_normals(V,F,N);
+//    V_border = igl::is_border_vertex(V,F);
+//    igl::vf(V,F,VF,VFi);
+//    igl::tt(V,F,TT,TTi);
+//  }
+//  
+//  void calculateMissmatch(Eigen::PlainObjectBase<DerivedO> &Handle_MMatch)
+//  {
+//    Handle_MMatch.setConstant(F.rows(),3,-1);
+//    for (unsigned int i=0;i<F.rows();i++)
+//    {
+//      for (int j=0;j<3;j++)
+//      {
+//        if (i==TT(i,j) || TT(i,j) == -1)
+//          Handle_MMatch(i,j)=0;
+//        else
+//          Handle_MMatch(i,j) = MissMatchByCross(i,TT(i,j));
+//      }
+//    }
+//  }
+//  
+//};
+//}
+//template <typename DerivedV, typename DerivedF, typename DerivedO>
+//IGL_INLINE void igl::cross_field_missmatch(const Eigen::PlainObjectBase<DerivedV> &V,
+//                                           const Eigen::PlainObjectBase<DerivedF> &F,
+//                                           const Eigen::PlainObjectBase<DerivedV> &PD1,
+//                                           const Eigen::PlainObjectBase<DerivedV> &PD2,
+//                                           const bool isCombed,
+//                                           Eigen::PlainObjectBase<DerivedO> &missmatch)
+//{
+//  Eigen::PlainObjectBase<DerivedV> PD1_combed;
+//  Eigen::PlainObjectBase<DerivedV> PD2_combed;
+//  
+//  if (!isCombed)
+//    igl::comb_cross_field(V,F,PD1,PD2,PD1_combed,PD2_combed);
+//  else
+//  {
+//    PD1_combed = PD1;
+//    PD2_combed = PD2;
+//  }
+//  igl::MissMatchCalculator<DerivedV, DerivedF, DerivedO> sf(V, F, PD1_combed, PD2_combed);
+//  sf.calculateMissmatch(missmatch);
+//}

include/igl/list_to_matrix.cpp

@@ -121,5 +121,4 @@ template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<int,
 template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > >(std::vector<std::vector<double, std::allocator<double> >, std::allocator<std::vector<double, std::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template bool igl::list_to_matrix<int, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 template bool igl::list_to_matrix<double, Eigen::Matrix<double, 1, -1, 1, 1, -1> >(std::vector<double, std::allocator<double> > const&, Eigen::Matrix<double, 1, -1, 1, 1, -1>&);
-template bool igl::list_to_matrix<double, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > >(std::__1::vector<std::__1::vector<double, std::__1::allocator<double> >, std::__1::allocator<std::__1::vector<double, std::__1::allocator<double> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
 #endif

include/igl/planarize_quad_mesh.cpp

@@ -1,6 +1,7 @@
 #include "planarize_quad_mesh.h"
 #include "quad_planarity.h"
 #include <Eigen/Sparse>
+#include <iostream>
 
 namespace igl
 {
@@ -170,16 +171,17 @@ inline void igl::PlanarizerShapeUp<DerivedV, DerivedF>::assembleP()
       CC.col(i) = Vi.segment(3*i, 3);
     Eigen::Matrix<typename DerivedV::Scalar, 3, 3> C = CC*CC.transpose();
     
-    Eigen::EigenSolver<Eigen::Matrix<typename DerivedV::Scalar, 3, 3>> es(C);
-    // the real() is for compilation purposes
-    Eigen::Matrix<typename DerivedV::Scalar, 3, 1> lambda = es.eigenvalues().real();
-    Eigen::Matrix<typename DerivedV::Scalar, 3, 3> U = es.eigenvectors().real();
-    int min_i;
-    lambda.cwiseAbs().minCoeff(&min_i);
-    U.col(min_i).setZero();
-    Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> PP = U*U.transpose()*CC;
-    for (int i = 0; i <ni; ++i)
-      P.segment(3*ni*fi+3*i, 3) =  weightsSqrt[fi]*PP.col(i);
+    // Alec: Doesn't compile
+    //Eigen::EigenSolver<Eigen::Matrix<typename DerivedV::Scalar, 3, 3>> es(C);
+    //// the real() is for compilation purposes
+    //Eigen::Matrix<typename DerivedV::Scalar, 3, 1> lambda = es.eigenvalues().real();
+    //Eigen::Matrix<typename DerivedV::Scalar, 3, 3> U = es.eigenvectors().real();
+    //int min_i;
+    //lambda.cwiseAbs().minCoeff(&min_i);
+    //U.col(min_i).setZero();
+    //Eigen::Matrix<typename DerivedV::Scalar, Eigen::Dynamic, Eigen::Dynamic> PP = U*U.transpose()*CC;
+    //for (int i = 0; i <ni; ++i)
+    //  P.segment(3*ni*fi+3*i, 3) =  weightsSqrt[fi]*PP.col(i);
     
   }
 }
@@ -227,4 +229,4 @@ IGL_INLINE void igl::planarize_quad_mesh(const Eigen::PlainObjectBase<DerivedV>
 {
   PlanarizerShapeUp<DerivedV, DerivedF> planarizer(Vin, Fin, maxIter, threshold);
   planarizer.planarize(Vout);
-}
+}

include/igl/readOBJ.cpp

@@ -388,5 +388,4 @@ template bool igl::readOBJ<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matr
 template bool igl::readOBJ<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 2, 1, -1, 2> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 1, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >&);
 template bool igl::readOBJ<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template bool igl::readOBJ<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::string, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
-template bool igl::readOBJ<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 #endif

include/igl/readSTL.cpp

@@ -171,7 +171,7 @@ IGL_INLINE bool igl::readSTL(
     V.resize(num_tri*3,vector<TypeV >(3,0));
     N.resize(num_tri,vector<TypeN >(3,0));
     F.resize(num_tri,vector<TypeF >(3,0));
-    for(int t = 0;t<num_tri;t++)
+    for(int t = 0;t<(int)num_tri;t++)
     {
       // Read normal
       float n[3];

tutorial/203_CurvatureDirections/main.cpp

@@ -7,6 +7,9 @@
 #include <igl/barycenter.h>
 #include <igl/principal_curvature.h>
 #include <igl/jet.h>
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/invert_diag.h>
 
 Eigen::MatrixXd V;
 Eigen::MatrixXi F;
@@ -17,17 +20,29 @@ int main(int argc, char *argv[])
   // Load a mesh in OFF format
   igl::readOFF("../shared/fertility.off", V, F);
 
+  // Alternative discrete mean curvature
+  MatrixXd HN;
+  SparseMatrix<double> L,M,Minv;
+  igl::cotmatrix(V,F,L);
+  igl::massmatrix(V,F,igl::MASSMATRIX_VORONOI,M);
+  igl::invert_diag(M,Minv);
+  HN = -Minv*(L*V);
+  H = (HN.rowwise().squaredNorm()).array().sqrt();
+
   // Compute curvature directions via quadric fitting
   MatrixXd PD1,PD2;
-  VectorXd PV1,PV2;
+  VectorXd PV1,PV2,H;
   igl::principal_curvature(V,F,PD1,PD2,PV1,PV2);
+  // mean curvature
+  H = 0.5*(PV1+PV2);
 
   igl::Viewer viewer;
   viewer.set_mesh(V, F);
 
+
   // Compute pseudocolor
   MatrixXd C;
-  igl::jet(PV1,true,C);
+  igl::jet(H,true,C);
   viewer.set_colors(C);
 
   // Average edge length for sizing

tutorial/readme.md

@@ -68,7 +68,7 @@ on a triangle mesh is via a vertex's _angular deficit_:
  $k_G(v_i) = 2π - \sum\limits_{j\in N(i)}θ_{ij},$
 
 where $N(i)$ are the triangles incident on vertex $i$ and $θ_{ij}$ is the angle
-at vertex $i$ in triangle $j$.
+at vertex $i$ in triangle $j$. (**Alec: cite Meyer or something**)
 
 Just like the continuous analog, our discrete Gaussian curvature reveals
 elliptic, hyperbolic and parabolic vertices on the domain.
@@ -76,8 +76,56 @@ elliptic, hyperbolic and parabolic vertices on the domain.
 ![The `GaussianCurvature` example computes discrete Gaussian curvature and visualizes it in
 pseudocolor.](images/bumpy-gaussian-curvature.jpg)
 
+## <a id=principal></a> Curvature Directions
+The two principal curvatures $(k_1,k_2)$ at a point on a surface measure how much the
+surface bends in different directions. The directions of maximum and minimum
+(signed) bending are call principal directions and are always
+orthogonal.
+
+Mean curvature is defined simply as the average of principal curvatures:
+
+ $H = \frac{1}{2}(k_1 + k_2).$ 
+
+One way to extract mean curvature is by examining the Laplace-Beltrami operator
+applied to the surface positions. The result is a so-called mean-curvature
+normal:
+
+  $-\Delta \mathbf{x} = H \mathbf{n}.$ 
+
+It is easy to compute this on a discrete triangle mesh in libigl using the cotangent
+Laplace-Beltrami operator (**Alec: cite Meyer**):
+
+```cpp
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/invert_diag.h>
+...
+MatrixXd HN;
+SparseMatrix<double> L,M,Minv;
+igl::cotmatrix(V,F,L);
+igl::massmatrix(V,F,igl::MASSMATRIX_VORONOI,M);
+igl::invert_diag(M,Minv);
+HN = -Minv*(L*V);
+H = (HN.rowwise().squaredNorm()).array().sqrt();
+```
+
+Combined with the angle defect definition of discrete Gaussian curvature, one
+can define principal curvatures and use least squares fitting to find
+directions (**Alec: cite meyer**).
+
+Alternatively, a robust method for determining principal curvatures is via
+quadric fitting (**Alec: cite whatever we're using**). In the neighborhood
+around every vertex, a best-fit quadric is found and principal curvature values
+and directions are sampled from this quadric. With these in tow, one can
+compute mean curvature and Gaussian curvature as sums and products
+respectively.
+
+![The `CurvatureDirections` example computes principal curvatures via quadric
+fitting and visualizes mean curvature in pseudocolor and principal directions
+with a cross field.](images/fertility-principal-curvature.jpg)
 
 This is an example of syntax highlighted code:
+
 ```cpp
 #include <foo.html>
 int main(int argc, char * argv[])

tutorial/style.css

@@ -17,6 +17,36 @@
    on the URL you're *sending* from, not the URL where the recipient views it.
 */
 /* .markdown-here-wrapper[data-md-url*="mail.yahoo."] ul { color: red; } */
+html
+{
+  margin: 0px;
+  padding: 0px;
+  /*background-color: #bbb;*/
+  background-image: url(images/background.gif);
+  font-weight: 300;
+  height: 100%;
+}
+body
+{
+  font-size: 11pt;
+  background-color: #fff;
+  padding: 10px;
+  margin:0px auto;
+  border-left: 1px solid #888;
+  border-right: 1px solid #888;
+  width: 765px;
+  text-align: left;
+}
+
+img
+{
+  max-width: 650px;
+}
+
+figcaption
+{
+  font-style: italic;
+}
 
 pre, code {
   font-size: 0.85em;