squared_edge_lengths (clean)

Former-commit-id: b16bf37e330d31a98960d5fdc4159245834c0cf5

include/igl/cotmatrix.h

@@ -43,8 +43,6 @@ namespace igl
   // therefore in general negative and the matrix is **negative** semi-definite
   // (immediately, -L is **positive** semi-definite)
   //
-  // Known bugs: off by 1e-16 on regular grid. I think its a problem of
-  // arithmetic order in cotmatrix_entries.h: C(i,e) = (arithmetic)/dblA/4
   template <typename DerivedV, typename DerivedF, typename Scalar>
   IGL_INLINE void cotmatrix(
     const Eigen::PlainObjectBase<DerivedV> & V, 

include/igl/cotmatrix_entries.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "cotmatrix_entries.h"
 #include "doublearea.h"
+#include "squared_edge_lengths.h"
 #include "edge_lengths.h"
 #include "face_areas.h"
 #include "volume.h"
@@ -34,11 +35,13 @@ IGL_INLINE void igl::cotmatrix_entries(
     case 3:
     {
       // Triangles
-      //Matrix<typename DerivedC::Scalar,Dynamic,3> l;
-      //edge_lengths(V,F,l);
-      // edge lengths numbered same as opposite vertices
+      //Compute Squared Edge lenghts 
+      Matrix<typename DerivedC::Scalar,Dynamic,3> l2;
+      igl::squared_edge_lengths(V,F,l2);
+      //Compute Edge lenghts 
       Matrix<typename DerivedC::Scalar,Dynamic,3> l;
-      igl::edge_lengths(V,F,l);
+      l = l2.array().sqrt();
+      
       // double area
       Matrix<typename DerivedC::Scalar,Dynamic,1> dblA;
       doublearea(l,dblA);
@@ -47,9 +50,9 @@ IGL_INLINE void igl::cotmatrix_entries(
       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;
+        C(i,0) = (l2(i,1) + l2(i,2) - l2(i,0))/dblA(i)/4.0;
+        C(i,1) = (l2(i,2) + l2(i,0) - l2(i,1))/dblA(i)/4.0;
+        C(i,2) = (l2(i,0) + l2(i,1) - l2(i,2))/dblA(i)/4.0;
       }
       break;
     }

include/igl/doublearea.cpp

@@ -146,10 +146,10 @@ IGL_INLINE void igl::doublearea(
   //
   // "Miscalculating Area and Angles of a Needle-like Triangle"
   // https://people.eecs.berkeley.edu/~wkahan/Triangle.pdf
-  sort(ul,2,false,l,_);
+  igl::sort(ul,2,false,l,_);
   // semiperimeters
-  Matrix<typename Derivedl::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
-  assert((Index)s.rows() == m);
+  //Matrix<typename Derivedl::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
+  //assert((Index)s.rows() == m);
   // resize output
   dblA.resize(l.rows(),1);
   parallel_for(

include/igl/edge_lengths.cpp

@@ -6,69 +6,18 @@
 // 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 "parallel_for.h"
-#include <iostream>
+#include "squared_edge_lengths.h"
 
 template <typename DerivedV, typename DerivedF, typename DerivedL>
 IGL_INLINE void igl::edge_lengths(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
   Eigen::PlainObjectBase<DerivedL>& L)
-{
-  using namespace std;
-  const int m = F.rows();
-  switch(F.cols())
   {
-    case 2:
-    {
-      L.resize(F.rows(),1);
-      for(int i = 0;i<F.rows();i++)
-      {
-        L(i,0) = (V.row(F(i,1))-V.row(F(i,0))).norm();
-      }
-      break;
-    }
-    case 3:
-    {
-      L.resize(m,3);
-      // loop over faces
-      parallel_for(
-        m,
-        [&V,&F,&L](const int i)
-        {
-          L(i,0) = (V.row(F(i,1))-V.row(F(i,2))).norm();
-          L(i,1) = (V.row(F(i,2))-V.row(F(i,0))).norm();
-          L(i,2) = (V.row(F(i,0))-V.row(F(i,1))).norm();
-        },
-        1000);
-      break;
-    }
-    case 4:
-    {
-      L.resize(m,6);
-      // loop over faces
-      parallel_for(
-        m,
-        [&V,&F,&L](const int i)
-        {
-          L(i,0) = (V.row(F(i,3))-V.row(F(i,0))).norm();
-          L(i,1) = (V.row(F(i,3))-V.row(F(i,1))).norm();
-          L(i,2) = (V.row(F(i,3))-V.row(F(i,2))).norm();
-          L(i,3) = (V.row(F(i,1))-V.row(F(i,2))).norm();
-          L(i,4) = (V.row(F(i,2))-V.row(F(i,0))).norm();
-          L(i,5) = (V.row(F(i,0))-V.row(F(i,1))).norm();
-        },
-        1000);
-      break;
-    }
-    default:
-    {
-      cerr<< "edge_lengths.h: Error: Simplex size ("<<F.cols()<<
-        ") not supported"<<endl;
-      assert(false);
-    }
+      igl::squared_edge_lengths(V,F,L);
+      L=L.array().sqrt().eval();
   }
-}
+  
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization

include/igl/internal_angles.cpp

@@ -7,7 +7,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 "internal_angles.h"
-#include "edge_lengths.h"
+#include "squared_edge_lengths.h"
 #include "parallel_for.h"
 #include "get_seconds.h"
 
@@ -26,11 +26,11 @@ IGL_INLINE void igl::internal_angles(
     Matrix<
       Scalar,
       DerivedF::RowsAtCompileTime,
-      DerivedF::ColsAtCompileTime> L;
-    edge_lengths(V,F,L);
+      DerivedF::ColsAtCompileTime> L_sq;
+      igl::squared_edge_lengths(V,F,L_sq);
 
     assert(F.cols() == 3 && "F should contain triangles");
-    internal_angles(L,K);
+    igl::internal_angles_using_squared_edge_lengths(L_sq,K);
   }else
   {
     assert(V.cols() == 3 && "If F contains non-triangle facets, V must be 3D");
@@ -63,10 +63,37 @@ IGL_INLINE void igl::internal_angles(
 }
 
 template <typename DerivedL, typename DerivedK>
-IGL_INLINE void igl::internal_angles(
+IGL_INLINE void igl::internal_angles_using_squared_edge_lengths(
+  const Eigen::PlainObjectBase<DerivedL>& L_sq,
+  Eigen::PlainObjectBase<DerivedK> & K)
+{
+  typedef typename DerivedL::Index Index;
+  assert(L_sq.cols() == 3 && "Edge-lengths should come from triangles");
+  const Index m = L_sq.rows();
+  K.resize(m,3);
+  parallel_for(
+    m,
+    [&L_sq,&K](const Index f)
+    {
+      for(size_t d = 0;d<3;d++)
+      {
+        const auto & s1 = L_sq(f,d);
+        const auto & s2 = L_sq(f,(d+1)%3);
+        const auto & s3 = L_sq(f,(d+2)%3);
+        K(f,d) = acos((s3 + s2 - s1)/(2.*sqrt(s3*s2)));
+      }
+    },
+    1000l);
+}
+
+template <typename DerivedL, typename DerivedK>
+IGL_INLINE void igl::internal_angles_using_edge_lengths(
   const Eigen::PlainObjectBase<DerivedL>& L,
   Eigen::PlainObjectBase<DerivedK> & K)
 {
+  // Note:
+  //   Usage of internal_angles_using_squared_edge_lengths() is preferred to internal_angles_using_squared_edge_lengths()
+  //   This function is deprecated and probably will be removed in future versions
   typedef typename DerivedL::Index Index;
   assert(L.cols() == 3 && "Edge-lengths should come from triangles");
   const Index m = L.rows();

include/igl/internal_angles.h

@@ -28,12 +28,30 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedF>& F,
     Eigen::PlainObjectBase<DerivedK> & K);
   // Inputs:
+  //   L_sq  #F by 3 list of squared edge lengths
+  // Output:
+  //   K  #F by poly-size eigen Matrix of internal angles
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  //
+  // Note:
+  //   Usage of internal_angles_using_squared_edge_lengths is preferred to internal_angles_using_squared_edge_lengths
+  template <typename DerivedL, typename DerivedK>
+  IGL_INLINE void internal_angles_using_squared_edge_lengths(
+    const Eigen::PlainObjectBase<DerivedL>& L_sq,
+    Eigen::PlainObjectBase<DerivedK> & K);
+  // Inputs:
   //   L  #F by 3 list of edge lengths
+  // Output:
+  //   K  #F by poly-size eigen Matrix of internal angles
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  //
+  // Note:
+  //   Usage of internal_angles_using_squared_edge_lengths is preferred to internal_angles_using_squared_edge_lengths
+  //   This function is deprecated and probably will be removed in future versions
   template <typename DerivedL, typename DerivedK>
-  IGL_INLINE void internal_angles(
+  IGL_INLINE void internal_angles_using_edge_lengths(
     const Eigen::PlainObjectBase<DerivedL>& L,
-    Eigen::PlainObjectBase<DerivedK> & K);
-}
+    Eigen::PlainObjectBase<DerivedK> & K);}
 
 #ifndef IGL_STATIC_LIBRARY
 #  include "internal_angles.cpp"

include/igl/squared_edge_lengths.cpp

@@ -0,0 +1,105 @@
+// 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
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "squared_edge_lengths.h"
+#include "parallel_for.h"
+#include <iostream>
+  
+template <typename DerivedV, typename DerivedF, typename DerivedL>
+IGL_INLINE void igl::squared_edge_lengths(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedL>& L)
+{
+  using namespace std;
+  const int m = F.rows();
+  switch(F.cols())
+  {
+    case 2:
+    {
+      L.resize(F.rows(),1);
+      for(int i = 0;i<F.rows();i++)
+      {
+        L(i,0) = (V.row(F(i,1))-V.row(F(i,0))).squaredNorm();
+      }
+      break;
+    }
+    case 3:
+    {
+      L.resize(m,3);
+      // loop over faces
+      parallel_for(
+        m,
+        [&V,&F,&L](const int i)
+        {
+          L(i,0) = (V.row(F(i,1))-V.row(F(i,2))).squaredNorm();
+          L(i,1) = (V.row(F(i,2))-V.row(F(i,0))).squaredNorm();
+          L(i,2) = (V.row(F(i,0))-V.row(F(i,1))).squaredNorm();
+        },
+        1000);
+      break;
+    }
+    case 4:
+    {
+      L.resize(m,6);
+      // loop over faces
+      parallel_for(
+        m,
+        [&V,&F,&L](const int i)
+        {
+          L(i,0) = (V.row(F(i,3))-V.row(F(i,0))).squaredNorm();
+          L(i,1) = (V.row(F(i,3))-V.row(F(i,1))).squaredNorm();
+          L(i,2) = (V.row(F(i,3))-V.row(F(i,2))).squaredNorm();
+          L(i,3) = (V.row(F(i,1))-V.row(F(i,2))).squaredNorm();
+          L(i,4) = (V.row(F(i,2))-V.row(F(i,0))).squaredNorm();
+          L(i,5) = (V.row(F(i,0))-V.row(F(i,1))).squaredNorm();
+        },
+        1000);
+      break;
+    }
+    default:
+    {
+      cerr<< "squared_edge_lengths.h: Error: Simplex size ("<<F.cols()<<
+        ") not supported"<<endl;
+      assert(false);
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::squared_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::squared_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::squared_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::squared_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::squared_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::squared_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::squared_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::squared_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> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -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<double, -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::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, 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, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -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<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(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, 6, 0, -1, 6> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, 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, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 3, 0, -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, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::squared_edge_lengths<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+#endif

include/igl/squared_edge_lengths.h

@@ -0,0 +1,49 @@
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_SQUARED_EDGE_LENGTHS_H
+#define IGL_SQUARED_EDGE_LENGTHS_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+
+namespace igl
+{
+  // Constructs a list of squared lengths of edges opposite each index in a face
+  // (triangle/tet) list
+  //
+  // Templates:
+  //   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
+  //   DerivedF 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
+  //   F  #F by 2 list of mesh edges
+  //    or
+  //   F  #F by 3 list of mesh faces (must be triangles)
+  //    or
+  //   T  #T by 4 list of mesh elements (must be tets)
+  // Outputs:
+  //   L  #F by {1|3|6} list of edge lengths squared
+  //     for edges, column of lengths
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  //     for tets, columns correspond to edges
+  //     [3 0],[3 1],[3 2],[1 2],[2 0],[0 1]
+  //
+  template <typename DerivedV, typename DerivedF, typename DerivedL>
+  IGL_INLINE void squared_edge_lengths(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedL>& L);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "squared_edge_lengths.cpp"
+#endif
+
+#endif
+