geometric cotangent, volume, face areas, edge lengths for tet meshes

Former-commit-id: 5ff193ae71315779e5ec0ebfe26696aae43edef7

documentation/implemented-papers.tex

@@ -38,6 +38,17 @@ choosing skinning transformations \cite{Jacobson:FAST:2012}.
 ``skeletal subspace deformation'', or ``enveloping''. This technique is often
 attributed to \cite{Magnenat-Thalmann:1988:JLD}.
 
+  % appears in the appendix of: ``Interactive Topology-aware Surface
+  % Reconstruction,'' by Sharf, A. et al
+  %
+  % versus
+  %
+  % ND derivation given in "A MONOTONE FINITE ELEMENT SCHEME FOR
+  % CONVECTION-DIFFUSION EQUATIONS" [Xu & ZIKATANOV 1999]
+  %
+  % 3D derivation given in "Aspects of unstructured grids and finite-volume
+  % solvers for the Euler and Navier-Stokes equations" [Barth 1992]
+
 \bibliographystyle{acmsiggraph}
 \bibliography{references} 
 

include/igl/cotangent.cpp

@@ -6,13 +6,20 @@
 // 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 "cotangent.h"
+#include "doublearea.h"
 #include "edge_lengths.h"
+#include "face_areas.h"
+#include "volume.h"
+#include "dihedral_angles.h"
 
 #include "verbose.h"
 
 
-template <class MatV, class MatF, class MatC>
-IGL_INLINE void igl::cotangent(const MatV & V, const MatF & F, MatC & C)
+template <typename DerivedV, typename DerivedF, typename DerivedC>
+IGL_INLINE void igl::cotangent(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedC>& C)
 {
   using namespace igl;
   using namespace std;
@@ -23,127 +30,145 @@ IGL_INLINE void igl::cotangent(const MatV & V, const MatF & F, MatC & C)
   int m = F.rows();
 
   // Law of cosines + law of sines
-  if(simplex_size == 3)
+  switch(simplex_size)
   {
-    // Triangles
-    //Matrix<typename MatC::Scalar,Dynamic,3> l;
-    //edge_lengths(V,F,l);
-    // edge lengths numbered same as opposite vertices
-    Matrix<typename MatC::Scalar,Dynamic,3> l(m,3);
-    // loop over faces
-    for(int i = 0;i<m;i++)
+    case 3:
     {
-      l(i,0) = sqrt((V.row(F(i,1))-V.row(F(i,2))).array().pow(2).sum());
-      l(i,1) = sqrt((V.row(F(i,2))-V.row(F(i,0))).array().pow(2).sum());
-      l(i,2) = sqrt((V.row(F(i,0))-V.row(F(i,1))).array().pow(2).sum());
+      // Triangles
+      //Matrix<typename DerivedC::Scalar,Dynamic,3> l;
+      //edge_lengths(V,F,l);
+      // edge lengths numbered same as opposite vertices
+      Matrix<typename DerivedC::Scalar,Dynamic,3> l;
+      igl::edge_lengths(V,F,l);
+      // double area
+      Matrix<typename DerivedC::Scalar,Dynamic,1> dblA;
+      doublearea(l,dblA);
+      // cotangents and diagonal entries for element matrices
+      // correctly divided by 4 (alec 2010)
+      C.resize(m,3);
+      for(int i = 0;i<m;i++)
+      {
+        C(i,0) = (l(i,1)*l(i,1) + l(i,2)*l(i,2) - l(i,0)*l(i,0))/dblA(i)/4.0;
+        C(i,1) = (l(i,2)*l(i,2) + l(i,0)*l(i,0) - l(i,1)*l(i,1))/dblA(i)/4.0;
+        C(i,2) = (l(i,0)*l(i,0) + l(i,1)*l(i,1) - l(i,2)*l(i,2))/dblA(i)/4.0;
+      }
+      break;
     }
-    // semiperimeters
-    Matrix<typename MatC::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
-    assert(s.rows() == m);
-    // Heron's formula for area
-    Matrix<typename MatC::Scalar,Dynamic,1> dblA(m);
-    for(int i = 0;i<m;i++)
+    case 4:
     {
-      dblA(i) = 2.0*sqrt(s(i)*(s(i)-l(i,0))*(s(i)-l(i,1))*(s(i)-l(i,2)));
-    }
-    // cotangents and diagonal entries for element matrices
-    // correctly divided by 4 (alec 2010)
-    C.resize(m,3);
-    for(int i = 0;i<m;i++)
-    {
-      C(i,0) = (l(i,1)*l(i,1) + l(i,2)*l(i,2) - l(i,0)*l(i,0))/dblA(i)/4.0;
-      C(i,1) = (l(i,2)*l(i,2) + l(i,0)*l(i,0) - l(i,1)*l(i,1))/dblA(i)/4.0;
-      C(i,2) = (l(i,0)*l(i,0) + l(i,1)*l(i,1) - l(i,2)*l(i,2))/dblA(i)/4.0;
-    }
-  }else if(simplex_size == 4)
-  {
-    // Tetrahedra
-    typedef Matrix<typename MatV::Scalar,3,1> Vec3;
-    typedef Matrix<typename MatV::Scalar,3,3> Mat3;
-    typedef Matrix<typename MatV::Scalar,3,4> Mat3x4;
-    typedef Matrix<typename MatV::Scalar,4,4> Mat4x4;
 
-    // preassemble right hand side
-    // COLUMN-MAJOR ORDER FOR LAPACK
-    Mat3x4 rhs;
-    rhs <<
-      1,0,0,-1,
-      0,1,0,-1,
-      0,0,1,-1;
+      // edge lengths numbered same as opposite vertices
+      Matrix<typename DerivedC::Scalar,Dynamic,6> l;
+      edge_lengths(V,F,l);
+      Matrix<typename DerivedC::Scalar,Dynamic,4> s;
+      face_areas(l,s);
+      Matrix<typename DerivedC::Scalar,Dynamic,6> cos_theta,theta;
+      dihedral_angles_intrinsic(l,s,theta,cos_theta);
 
-    bool diag_all_pos = true;
-    C.resize(m,6);
+      // volume
+      Matrix<typename DerivedC::Scalar,Dynamic,1> vol;
+      volume(l,vol);
 
-    // loop over tetrahedra
-    for(int j = 0;j<F.rows();j++)
-    {
-      // points a,b,c,d make up the tetrahedra
-      size_t a = F(j,0);
-      size_t b = F(j,1);
-      size_t c = F(j,2);
-      size_t d = F(j,3);
-      //const std::vector<double> & pa = vertices[a];
-      //const std::vector<double> & pb = vertices[b];
-      //const std::vector<double> & pc = vertices[c];
-      //const std::vector<double> & pd = vertices[d];
-      Vec3 pa = V.row(a);
-      Vec3 pb = V.row(b);
-      Vec3 pc = V.row(c);
-      Vec3 pd = V.row(d);
-
-      // Following definition that appears in the appendix of: ``Interactive
-      // Topology-aware Surface Reconstruction,'' by Sharf, A. et al
-      // http://www.cs.bgu.ac.il/~asharf/Projects/InSuRe/Insure_siggraph_final.pdf
 
-      // compute transpose of jacobian Jj
-      Mat3 JTj;
-      JTj.row(0) = pa-pd;
-      JTj.row(1) = pb-pd;
-      JTj.row(2) = pc-pd;
+      // Law of sines
+      // http://mathworld.wolfram.com/Tetrahedron.html
+      Matrix<typename DerivedC::Scalar,Dynamic,6> sin_theta(m,6);
+      sin_theta.col(0) = vol.array() / ((2./(3.*l.col(0).array())).array() * s.col(1).array() * s.col(2).array());
+      sin_theta.col(1) = vol.array() / ((2./(3.*l.col(1).array())).array() * s.col(2).array() * s.col(0).array());
+      sin_theta.col(2) = vol.array() / ((2./(3.*l.col(2).array())).array() * s.col(0).array() * s.col(1).array());
+      sin_theta.col(3) = vol.array() / ((2./(3.*l.col(3).array())).array() * s.col(3).array() * s.col(0).array());
+      sin_theta.col(4) = vol.array() / ((2./(3.*l.col(4).array())).array() * s.col(3).array() * s.col(1).array());
+      sin_theta.col(5) = vol.array() / ((2./(3.*l.col(5).array())).array() * s.col(3).array() * s.col(2).array());
 
-      // compute abs(determinant of JTj)/6 (volume of tet)
-      // determinant of transpose of A equals determinant of A
-      double volume = fabs(JTj.determinant())/6.0;
-      //printf("volume[%d] = %g\n",j+1,volume);
 
-      // solve Jj' * Ej = [-I -1], for Ej
-      // in other words solve JTj * Ej = [-I -1], for Ej
-      Mat3x4 Ej = JTj.inverse() * rhs;
-      // compute Ej'*Ej
-      Mat4x4 EjTEj = Ej.transpose() * Ej;
+      // http://arxiv.org/pdf/1208.0354.pdf Page 18
+      C = (1./6.) * l.array() * cos_theta.array() / sin_theta.array();
 
-      // Kj =  det(JTj)/6 * Ej'Ej 
-      Mat4x4 Kj = EjTEj*volume;
-      diag_all_pos &= ((Kj(0,0)>0) & (Kj(1,1)>0)) & ((Kj(2,2)>0) & (Kj(3,3)>0));
-      C(j,0) = Kj(1,2);
-      C(j,1) = Kj(2,0);
-      C(j,2) = Kj(0,1);
-      C(j,3) = Kj(3,0);
-      C(j,4) = Kj(3,1);
-      C(j,5) = Kj(3,2);
+// LEGACY
+//      // Tetrahedra
+//      typedef Matrix<typename MatV::Scalar,3,1> Vec3;
+//      typedef Matrix<typename MatV::Scalar,3,3> Mat3;
+//      typedef Matrix<typename MatV::Scalar,3,4> Mat3x4;
+//      typedef Matrix<typename MatV::Scalar,4,4> Mat4x4;
+//
+//      // preassemble right hand side
+//      // COLUMN-MAJOR ORDER FOR LAPACK
+//      Mat3x4 rhs;
+//      rhs <<
+//        1,0,0,-1,
+//        0,1,0,-1,
+//        0,0,1,-1;
+//
+//      bool diag_all_pos = true;
+//      C.resize(m,6);
+//
+//      // loop over tetrahedra
+//      for(int j = 0;j<F.rows();j++)
+//      {
+//        // points a,b,c,d make up the tetrahedra
+//        size_t a = F(j,0);
+//        size_t b = F(j,1);
+//        size_t c = F(j,2);
+//        size_t d = F(j,3);
+//        //const std::vector<double> & pa = vertices[a];
+//        //const std::vector<double> & pb = vertices[b];
+//        //const std::vector<double> & pc = vertices[c];
+//        //const std::vector<double> & pd = vertices[d];
+//        Vec3 pa = V.row(a);
+//        Vec3 pb = V.row(b);
+//        Vec3 pc = V.row(c);
+//        Vec3 pd = V.row(d);
+//
+//        // Following definition that appears in the appendix of: ``Interactive
+//        // Topology-aware Surface Reconstruction,'' by Sharf, A. et al
+//        // http://www.cs.bgu.ac.il/~asharf/Projects/InSuRe/Insure_siggraph_final.pdf
+//
+//        // compute transpose of jacobian Jj
+//        Mat3 JTj;
+//        JTj.row(0) = pa-pd;
+//        JTj.row(1) = pb-pd;
+//        JTj.row(2) = pc-pd;
+//
+//        // compute abs(determinant of JTj)/6 (volume of tet)
+//        // determinant of transpose of A equals determinant of A
+//        double volume = fabs(JTj.determinant())/6.0;
+//        //printf("volume[%d] = %g\n",j+1,volume);
+//
+//        // solve Jj' * Ej = [-I -1], for Ej
+//        // in other words solve JTj * Ej = [-I -1], for Ej
+//        Mat3x4 Ej = JTj.inverse() * rhs;
+//        // compute Ej'*Ej
+//        Mat4x4 EjTEj = Ej.transpose() * Ej;
+//
+//        // Kj =  det(JTj)/6 * Ej'Ej 
+//        Mat4x4 Kj = EjTEj*volume;
+//        diag_all_pos &= ((Kj(0,0)>0) & (Kj(1,1)>0)) & ((Kj(2,2)>0) & (Kj(3,3)>0));
+//        C(j,0) = Kj(1,2);
+//        C(j,1) = Kj(2,0);
+//        C(j,2) = Kj(0,1);
+//        C(j,3) = Kj(3,0);
+//        C(j,4) = Kj(3,1);
+//        C(j,5) = Kj(3,2);
+//      }
+//      if(diag_all_pos)
+//      {
+//#ifdef VERBOSE 
+//        verbose("cotangent.h: Flipping sign of cotangent, so that cots are positive\n");
+//#endif
+//        C *= -1.0;
+//      }
+      break;
     }
-    if(diag_all_pos)
+    default:
     {
-#ifdef VERBOSE 
-      verbose("cotangent.h: Flipping sign of cotangent, so that cots are positive\n");
-#endif
-      C *= -1.0;
+      fprintf(stderr,
+          "cotangent.h: Error: Simplex size (%d) not supported\n", simplex_size);
+      assert(false);
     }
-  }else
-  {
-    fprintf(stderr,
-      "cotangent.h: Error: Simplex size (%d) not supported\n", simplex_size);
-    assert(false);
   }
 }
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
-// generated by autoexplicit.sh
-template void igl::cotangent<Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-// generated by autoexplicit.sh
-template void igl::cotangent<Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-// generated by autoexplicit.sh
-template void igl::cotangent<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> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-template void igl::cotangent<Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::PlainObjectBase<Eigen::Matrix<int, -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::Matrix<double, -1, -1, 0, -1, -1>&);
+template void igl::cotangent<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::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> >&);
 #endif

include/igl/cotangent.h

@@ -8,6 +8,7 @@
 #ifndef IGL_COTANGENT_H
 #define IGL_COTANGENT_H
 #include "igl_inline.h"
+#include <Eigen/Core>
 namespace igl
 {
   // COTANGENT compute the cotangents of each angle in mesh (V,F)
@@ -22,12 +23,17 @@ namespace igl
   // Outputs:
   //   C  #F by {3|6} list of cotangents corresponding angles
   //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
-  //     for tets, columns correspond to edges [1,2],[2,0],[0,1],[3,0],[3,1],[3,2]
+  //     for tets, columns correspond to edges
+  //     [1,2],[2,0],[0,1],[3,0],[3,1],[3,2] **times corresponding edge
+  //     lengths**
   //
   // Known bug:
   //   This computes 0.5*cotangent
-  template <class MatV, class MatF, class MatC>
-  IGL_INLINE void cotangent(const MatV & V, const MatF & F, MatC & C);
+  template <typename DerivedV, typename DerivedF, typename DerivedC>
+  IGL_INLINE void cotangent(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedC>& C);
 }
 
 #ifdef IGL_HEADER_ONLY

include/igl/cotmatrix.cpp

@@ -17,8 +17,8 @@
 
 template <typename DerivedV, typename DerivedF, typename Scalar>
 IGL_INLINE void igl::cotmatrix(
-  const Eigen::MatrixBase<DerivedV> & V, 
-  const Eigen::MatrixBase<DerivedF> & F, 
+  const Eigen::PlainObjectBase<DerivedV> & V, 
+  const Eigen::PlainObjectBase<DerivedF> & F, 
   Eigen::SparseMatrix<Scalar>& L)
 {
   using namespace igl;
@@ -81,5 +81,5 @@ IGL_INLINE void igl::cotmatrix(
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
-template void igl::cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/cotmatrix.h

@@ -43,8 +43,8 @@ namespace igl
   // arithmetic order in cotangent.h: C(i,e) = (arithmetic)/dblA/4
   template <typename DerivedV, typename DerivedF, typename Scalar>
   IGL_INLINE void cotmatrix(
-    const Eigen::MatrixBase<DerivedV> & V, 
-    const Eigen::MatrixBase<DerivedF> & F, 
+    const Eigen::PlainObjectBase<DerivedV> & V, 
+    const Eigen::PlainObjectBase<DerivedF> & F, 
     Eigen::SparseMatrix<Scalar>& L);
 }
 

include/igl/dihedral_angles.cpp

@@ -0,0 +1,94 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "dihedral_angles.h"
+#include <cassert>
+
+template <
+  typename DerivedV, 
+  typename DerivedT, 
+  typename Derivedtheta,
+  typename Derivedcos_theta>
+IGL_INLINE void igl::dihedral_angles(
+  Eigen::PlainObjectBase<DerivedV>& V,
+  Eigen::PlainObjectBase<DerivedT>& T,
+  Eigen::PlainObjectBase<Derivedtheta>& theta,
+  Eigen::PlainObjectBase<Derivedcos_theta>& cos_theta)
+{
+  using namespace Eigen;
+  assert(T.cols() == 4);
+  Matrix<typename Derivedtheta::Scalar,Dynamic,6> l;
+  edge_lengths(V,T,l);
+  Matrix<typename Derivedtheta::Scalar,Dynamic,4> s;
+  face_areas(l,s);
+  return dihedral_angles_intrinsic(l,s,theta,cos_theta);
+}
+
+template <
+  typename DerivedL, 
+  typename DerivedA, 
+  typename Derivedtheta,
+  typename Derivedcos_theta>
+IGL_INLINE void igl::dihedral_angles_intrinsic(
+  Eigen::PlainObjectBase<DerivedL>& L,
+  Eigen::PlainObjectBase<DerivedA>& A,
+  Eigen::PlainObjectBase<Derivedtheta>& theta,
+  Eigen::PlainObjectBase<Derivedcos_theta>& cos_theta)
+{
+  using namespace Eigen;
+  const int m = L.rows();
+  assert(m == A.rows());
+  // Law of cosines
+  // http://math.stackexchange.com/a/49340/35376
+  Matrix<typename Derivedtheta::Scalar,Dynamic,6> H_sqr(m,6);
+  H_sqr.col(0) = (1./16.) * (4. * L.col(3).array().square() * L.col(0).array().square() - 
+                                ((L.col(1).array().square() + L.col(4).array().square()) -
+                                 (L.col(2).array().square() + L.col(5).array().square())).square());
+  H_sqr.col(1) = (1./16.) * (4. * L.col(4).array().square() * L.col(1).array().square() - 
+                                ((L.col(2).array().square() + L.col(5).array().square()) -
+                                 (L.col(3).array().square() + L.col(0).array().square())).square());
+  H_sqr.col(2) = (1./16.) * (4. * L.col(5).array().square() * L.col(2).array().square() - 
+                                ((L.col(3).array().square() + L.col(0).array().square()) -
+                                 (L.col(4).array().square() + L.col(1).array().square())).square());
+  H_sqr.col(3) = (1./16.) * (4. * L.col(0).array().square() * L.col(3).array().square() - 
+                                ((L.col(4).array().square() + L.col(1).array().square()) -
+                                 (L.col(5).array().square() + L.col(2).array().square())).square());
+  H_sqr.col(4) = (1./16.) * (4. * L.col(1).array().square() * L.col(4).array().square() - 
+                                ((L.col(5).array().square() + L.col(2).array().square()) -
+                                 (L.col(0).array().square() + L.col(3).array().square())).square());
+  H_sqr.col(5) = (1./16.) * (4. * L.col(2).array().square() * L.col(5).array().square() - 
+                                ((L.col(0).array().square() + L.col(3).array().square()) -
+                                 (L.col(1).array().square() + L.col(4).array().square())).square());
+  cos_theta.resize(m,6);
+  cos_theta.col(0) = (H_sqr.col(0).array() - 
+      A.col(1).array().square() - A.col(2).array().square()).array() / 
+      (-2.*A.col(1).array() * A.col(2).array());
+  cos_theta.col(1) = (H_sqr.col(1).array() - 
+      A.col(2).array().square() - A.col(0).array().square()).array() / 
+      (-2.*A.col(2).array() * A.col(0).array());
+  cos_theta.col(2) = (H_sqr.col(2).array() - 
+      A.col(0).array().square() - A.col(1).array().square()).array() / 
+      (-2.*A.col(0).array() * A.col(1).array());
+  cos_theta.col(3) = (H_sqr.col(3).array() - 
+      A.col(3).array().square() - A.col(0).array().square()).array() / 
+      (-2.*A.col(3).array() * A.col(0).array());
+  cos_theta.col(4) = (H_sqr.col(4).array() - 
+      A.col(3).array().square() - A.col(1).array().square()).array() / 
+      (-2.*A.col(3).array() * A.col(1).array());
+  cos_theta.col(5) = (H_sqr.col(5).array() - 
+      A.col(3).array().square() - A.col(2).array().square()).array() / 
+      (-2.*A.col(3).array() * A.col(2).array());
+
+  theta = cos_theta.array().acos();
+
+  cos_theta.resize(m,6);
+
+}
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::dihedral_angles_intrinsic<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<double, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
+#endif

include/igl/dihedral_angles.h

@@ -0,0 +1,55 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_DIHEDRAL_ANGLES_H
+#define IGL_DIHEDRAL_ANGLES_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // DIHEDRAL_ANGLES Compute dihedral angles for all tets of a given tet mesh
+  // (V,T)
+  //
+  // theta = dihedral_angles(V,T)
+  // theta = dihedral_angles(V,T,'ParameterName',parameter_value,...)
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   T  #V by 4 list of tet indices
+  // Outputs:
+  //   theta  #T by 6 list of dihedral angles (in radians)
+  //   cos_theta  #T by 6 list of cosine of dihedral angles (in radians)
+  //
+  template <
+    typename DerivedV, 
+    typename DerivedT, 
+    typename Derivedtheta,
+    typename Derivedcos_theta>
+  IGL_INLINE void dihedral_angles(
+    Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedT>& T,
+    Eigen::PlainObjectBase<Derivedtheta>& theta,
+    Eigen::PlainObjectBase<Derivedcos_theta>& cos_theta);
+  template <
+    typename DerivedL, 
+    typename DerivedA, 
+    typename Derivedtheta,
+    typename Derivedcos_theta>
+  IGL_INLINE void dihedral_angles_intrinsic(
+    Eigen::PlainObjectBase<DerivedL>& L,
+    Eigen::PlainObjectBase<DerivedA>& A,
+    Eigen::PlainObjectBase<Derivedtheta>& theta,
+    Eigen::PlainObjectBase<Derivedcos_theta>& cos_theta);
+
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "dihedral_angles.cpp"
+#endif
+
+#endif
+

include/igl/doublearea.cpp

@@ -101,6 +101,8 @@ IGL_INLINE void igl::doublearea(
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<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::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> >&);
 template void igl::doublearea<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<double, -1, 1, 0, -1, 1> >&);

include/igl/edge_lengths.cpp

@@ -6,6 +6,7 @@
 // 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 "edge_lengths.h"
+#include <iostream>
 
 template <typename DerivedV, typename DerivedF, typename DerivedL>
 IGL_INLINE void igl::edge_lengths(
@@ -13,20 +14,61 @@ IGL_INLINE void igl::edge_lengths(
   const Eigen::PlainObjectBase<DerivedF>& F,
   Eigen::PlainObjectBase<DerivedL>& L)
 {
-  assert(F.cols() == 3);
-  L.resize(F.rows(),3);
-  // loop over faces
-  for(int i = 0;i<F.rows();i++)
+  using namespace std;
+  switch(F.cols())
   {
-    L(i,0) = sqrt((V.row(F(i,1))-V.row(F(i,2))).array().pow(2).sum());
-    L(i,1) = sqrt((V.row(F(i,2))-V.row(F(i,0))).array().pow(2).sum());
-    L(i,2) = sqrt((V.row(F(i,0))-V.row(F(i,1))).array().pow(2).sum());
+    case 3:
+    {
+      L.resize(F.rows(),3);
+      // loop over faces
+      for(int i = 0;i<F.rows();i++)
+      {
+        L(i,0) = sqrt((V.row(F(i,1))-V.row(F(i,2))).array().pow(2).sum());
+        L(i,1) = sqrt((V.row(F(i,2))-V.row(F(i,0))).array().pow(2).sum());
+        L(i,2) = sqrt((V.row(F(i,0))-V.row(F(i,1))).array().pow(2).sum());
+      }
+      break;
+    }
+    case 4:
+    {
+      const int m = F.rows();
+      L.resize(m,6);
+      // loop over faces
+      for(int i = 0;i<m;i++)
+      {
+        L(i,0) = sqrt((V.row(F(i,3))-V.row(F(i,0))).array().pow(2).sum());
+        L(i,1) = sqrt((V.row(F(i,3))-V.row(F(i,1))).array().pow(2).sum());
+        L(i,2) = sqrt((V.row(F(i,3))-V.row(F(i,2))).array().pow(2).sum());
+        L(i,3) = sqrt((V.row(F(i,1))-V.row(F(i,2))).array().pow(2).sum());
+        L(i,4) = sqrt((V.row(F(i,2))-V.row(F(i,0))).array().pow(2).sum());
+        L(i,5) = sqrt((V.row(F(i,0))-V.row(F(i,1))).array().pow(2).sum());
+      }
+      break;
+    }
+    default:
+    {
+      cerr<< "edge_lengths.h: Error: Simplex size ("<<F.cols()<<
+        ") not supported"<<endl;
+      assert(false);
+    }
   }
 }
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<float, -1, 2, 0, -1, 2>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> >&);
+// generated by autoexplicit.sh
+template void igl::edge_lengths<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::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> >&);
+// generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
+// generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
+// generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(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, 6, 0, -1, 6> >&);
+// generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(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, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
 template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
 template void igl::edge_lengths<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::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> >&);
 #endif

include/igl/edge_lengths.h

@@ -14,7 +14,8 @@
 namespace igl
 {
   // Constructs a list of lengths of edges opposite each index in a face
-  // (triangle) list
+  // (triangle/tet) list
+  //
   // Templates:
   //   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
   //   DerivedF derived from face indices matrix type: i.e. MatrixXi
@@ -22,10 +23,14 @@ namespace igl
   // Inputs:
   //   V  eigen matrix #V by 3
   //   F  #F by 3 list of mesh faces (must be triangles)
+  //    or
+  //   T  #T by 4 list of mesh elements (must be tets)
   // Outputs:
-  //   E #E by 2 list of edges in no particular order
+  //   L  #F by {3|6} list of edge lengths 
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  //     for tets, columns correspond to edges
+  //     [1,2],[2,0],[0,1],[3,0],[3,1],[3,2]
   //
-  // See also: adjacency_matrix
   template <typename DerivedV, typename DerivedF, typename DerivedL>
   IGL_INLINE void edge_lengths(
     const Eigen::PlainObjectBase<DerivedV>& V,

include/igl/face_areas.cpp

@@ -0,0 +1,48 @@
+#include "face_areas.h"
+#include "edge_lengths.h"
+#include "doublearea.h"
+
+template <typename DerivedV, typename DerivedT, typename DerivedA>
+IGL_INLINE void igl::face_areas(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedT>& T,
+  Eigen::PlainObjectBase<DerivedA>& A)
+{
+  assert(T.cols() == 4);
+  typename Eigen::PlainObjectBase<DerivedA> L;
+  edge_lengths(V,T,L);
+  return face_areas(L,A);
+}
+
+template <typename DerivedL, typename DerivedA>
+IGL_INLINE void igl::face_areas(
+  const Eigen::PlainObjectBase<DerivedL>& L,
+  Eigen::PlainObjectBase<DerivedA>& A)
+{
+  using namespace Eigen;
+  assert(L.cols() == 6);
+  const int m = L.rows();
+  // (unsigned) face Areas (opposite vertices: 1 2 3 4)
+  Matrix<typename DerivedA::Scalar,Dynamic,1> 
+    A0(m,1), A1(m,1), A2(m,1), A3(m,1);
+  Matrix<typename DerivedA::Scalar,Dynamic,3> 
+    L0(m,3), L1(m,3), L2(m,3), L3(m,3);
+  L0<<L.col(1),L.col(2),L.col(3);
+  L1<<L.col(0),L.col(2),L.col(4);
+  L2<<L.col(0),L.col(1),L.col(5);
+  L3<<L.col(3),L.col(4),L.col(5);
+  doublearea(L0,A0);
+  doublearea(L1,A1);
+  doublearea(L2,A2);
+  doublearea(L3,A3);
+  A.resize(m,4);
+  A.col(0) = 0.5*A0;
+  A.col(1) = 0.5*A1;
+  A.col(2) = 0.5*A2;
+  A.col(3) = 0.5*A3;
+}
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::face_areas<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<double, -1, 4, 0, -1, 4> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&);
+#endif

include/igl/face_areas.h

@@ -0,0 +1,46 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_FACE_AREAS_H
+#define IGL_FACE_AREAS_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+
+namespace igl
+{
+  // Constructs a list of face areas of faces opposite each index in a tet list
+  //
+  // Templates:
+  //   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
+  //   DerivedT derived from face indices matrix type: i.e. MatrixXi
+  //   DerivedL derived from edge lengths matrix type: i.e. MatrixXd
+  // Inputs:
+  //   V  eigen matrix #V by 3
+  //   T  #T by 3 list of mesh faces (must be triangles)
+  // Outputs:
+  //   E #E by 2 list of edges in no particular order
+  //
+  // See also: adjacency_matrix
+  template <typename DerivedV, typename DerivedT, typename DerivedA>
+  IGL_INLINE void face_areas(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedT>& T,
+    Eigen::PlainObjectBase<DerivedA>& A);
+  template <typename DerivedL, typename DerivedA>
+  IGL_INLINE void face_areas(
+    const Eigen::PlainObjectBase<DerivedL>& L,
+    Eigen::PlainObjectBase<DerivedA>& A);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "face_areas.cpp"
+#endif
+
+#endif
+
+

include/igl/sort.cpp

@@ -213,6 +213,12 @@ IGL_INLINE void igl::sort(
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<float, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > > const&, bool, std::vector<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
@@ -225,4 +231,5 @@ template void igl::sort<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<in
 template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
+template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/volume.cpp

@@ -0,0 +1,72 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "volume.h"
+template <
+  typename DerivedV, 
+  typename DerivedT, 
+  typename Derivedvol>
+IGL_INLINE void igl::volume(
+  Eigen::PlainObjectBase<DerivedV>& V,
+  Eigen::PlainObjectBase<DerivedT>& T,
+  Eigen::PlainObjectBase<Derivedvol>& vol)
+{
+  using namespace Eigen;
+  const int m = T.rows();
+  vol.resize(m,1);
+  for(int t = 0;t<m;t++)
+  {
+    const auto & a = V.row(T(t,0));
+    const auto & b = V.row(T(t,1));
+    const auto & c = V.row(T(t,2));
+    const auto & d = V.row(T(t,3));
+    vol(t) = -(a-d).dot((b-d).cross(c-d))/6.;
+  }
+}
+
+template <
+  typename DerivedL, 
+  typename Derivedvol>
+IGL_INLINE void igl::volume(
+  Eigen::PlainObjectBase<DerivedL>& L,
+  Eigen::PlainObjectBase<Derivedvol>& vol)
+{
+  using namespace Eigen;
+  const int m = L.rows();
+  typedef typename Derivedvol::Scalar ScalarS;
+  vol.resize(m,1);
+  for(int t = 0;t<m;t++)
+  {
+    const ScalarS u = L(t,0);
+    const ScalarS v = L(t,1);
+    const ScalarS w = L(t,2);
+    const ScalarS U = L(t,3);
+    const ScalarS V = L(t,4);
+    const ScalarS W = L(t,5);
+    const ScalarS X = (w - U + v)*(U + v + w);
+    const ScalarS x = (U - v + w)*(v - w + U);
+    const ScalarS Y = (u - V + w)*(V + w + u);
+    const ScalarS y = (V - w + u)*(w - u + V);
+    const ScalarS Z = (v - W + u)*(W + u + v);
+    const ScalarS z = (W - u + v)*(u - v + W);
+    const ScalarS a = sqrt(x*Y*Z); 
+    const ScalarS b = sqrt(y*Z*X); 
+    const ScalarS c = sqrt(z*X*Y); 
+    const ScalarS d = sqrt(x*y*z); 
+    vol(t) = sqrt(
+       (-a + b + c + d)*
+       ( a - b + c + d)*
+       ( a + b - c + d)*
+       ( a + b + c - d))/
+       (192.*u*v*w);
+  }
+}
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::volume<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/volume.h

@@ -0,0 +1,51 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_VOLUME_H
+#define IGL_VOLUME_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // VOLUME Compute volume for all tets of a given tet mesh
+  // (V,T)
+  //
+  // vol = volume(V,T)
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   T  #V by 4 list of tet indices
+  // Outputs:
+  //   vol  #T list of dihedral angles (in radians)
+  //
+  template <
+    typename DerivedV, 
+    typename DerivedT, 
+    typename Derivedvol>
+  IGL_INLINE void volume(
+    Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedT>& T,
+    Eigen::PlainObjectBase<Derivedvol>& vol);
+  // Intrinsic version:
+  //
+  // Inputs:
+  //   L  #V by 6 list of edge lengths (see edge_lengths)
+  template <
+    typename DerivedL, 
+    typename Derivedvol>
+  IGL_INLINE void volume(
+    Eigen::PlainObjectBase<DerivedL>& L,
+    Eigen::PlainObjectBase<Derivedvol>& vol);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "volume.cpp"
+#endif
+
+#endif
+
+