added intrinsic version of cotmatrix_entries + test

include/igl/cotmatrix_entries.cpp

@@ -1,6 +1,7 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2018 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 
@@ -50,6 +51,7 @@ IGL_INLINE void igl::cotmatrix_entries(
       C.resize(m,3);
       for(int i = 0;i<m;i++)
       {
+        // Alec: I'm doubtful that using l2 here is actually improving numerics.
         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;
@@ -97,14 +99,45 @@ IGL_INLINE void igl::cotmatrix_entries(
   }
 }
 
+template <typename Derivedl, typename DerivedC>
+IGL_INLINE void igl::cotmatrix_entries(
+  const Eigen::MatrixBase<Derivedl>& l,
+  Eigen::PlainObjectBase<DerivedC>& C)
+{
+  using namespace Eigen;
+  const int m = l.rows();
+  assert(l.cols() == 4 && "Only triangles accepted");
+  //Compute squared Edge lengths 
+  Matrix<typename DerivedC::Scalar,Dynamic,3> l2;
+  l2 = l.array().square();
+  // Alec: It's a little annoying that there's duplicate code here. The
+  // "extrinic" version above is first computing squared edge lengths, taking
+  // the square root and calling this. We can't have a cotmatrix_entries(l,l2,C)
+  // overload because it will confuse Eigen with the cotmatrix_entries(V,F,C)
+  // overload. In the end, I'd like to be convinced that using l2 directly above
+  // is actually better numerically (or significantly faster) than just calling
+  // edge_lengths and this cotmatrix_entries(l,C);
+  //
+  // double area
+  Matrix<typename DerivedC::Scalar,Dynamic,1> dblA;
+  doublearea(l,0.,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++)
+  {
+    // Alec: I'm doubtful that using l2 here is actually improving numerics.
+    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;
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
-// generated by autoexplicit.sh
+template void igl::cotmatrix_entries<Eigen::Matrix<double, -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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::cotmatrix_entries<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::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
-// generated by autoexplicit.sh
 template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/cotmatrix_entries.h

@@ -1,6 +1,7 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
 // Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2018 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 
@@ -28,6 +29,17 @@ namespace igl
     const Eigen::MatrixBase<DerivedV>& V,
     const Eigen::MatrixBase<DerivedF>& F,
     Eigen::PlainObjectBase<DerivedC>& C);
+  // Intrinsic version.
+  //
+  // Inputs:
+  //   l  #F by 3 list of triangle edge lengths (see edge_lengths)
+  // Outputs:
+  //   C  #F by 3 list of 1/2*cotangents corresponding angles
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  template <typename Derivedl, typename DerivedC>
+  IGL_INLINE void cotmatrix_entries(
+    const Eigen::MatrixBase<Derivedl>& l,
+    Eigen::PlainObjectBase<DerivedC>& C);
 }
 
 #ifndef IGL_STATIC_LIBRARY

tests/include/igl/cotmatrix_entries.cpp

@@ -1,5 +1,7 @@
 #include <test_common.h>
 #include <igl/cotmatrix_entries.h>
+#include <igl/edge_lengths.h>
+#include <igl/EPS.h>
 
 TEST(cotmatrix_entries, simple)
 {
@@ -30,7 +32,7 @@ TEST(cotmatrix_entries, simple)
   for(int f = 0;f<C1.rows();f++)
   {
 #ifdef IGL_EDGE_LENGTHS_SQUARED_H
-      //Hard assert if we have edge_lenght_squared
+      //Hard assert if we have edge_length_squared
       for(int v = 0;v<3;v++)
         if (C1(f,v) > 0.1)
             ASSERT_EQ(0.5, C1(f,v));
@@ -39,7 +41,7 @@ TEST(cotmatrix_entries, simple)
        //All cotangents sum 1.0 for those triangles
        ASSERT_EQ(1.0, C1.row(f).sum());
 #else
-      //Soft assert if we have not edge_lenght_squared
+      //Soft assert if we have not edge_length_squared
       for(int v = 0;v<3;v++)
         if (C1(f,v) > 0.1)
             ASSERT_NEAR(0.5, C1(f,v), epsilon);
@@ -72,7 +74,7 @@ TEST(cotmatrix_entries, simple)
   for(int f = 0;f<C1.rows();f++)
   {    
 #ifdef IGL_EDGE_LENGTHS_SQUARED_H
-      //Hard assert if we have edge_lenght_squared
+      //Hard assert if we have edge_length_squared
       for(int v = 0;v<3;v++)
         if (C1(f,v) > 0.1)
             ASSERT_EQ(0.5, C1(f,v));
@@ -81,7 +83,7 @@ TEST(cotmatrix_entries, simple)
        //All cotangents sum 1.0 for those triangles
        ASSERT_EQ(1.0, C1.row(f).sum());
 #else
-      //Soft assert if we have not edge_lenght_squared
+      //Soft assert if we have not edge_length_squared
       for(int v = 0;v<3;v++)
         if (C1(f,v) > 0.1)
             ASSERT_NEAR(0.5, C1(f,v), epsilon);
@@ -134,3 +136,19 @@ TEST(cotmatrix_entries, simple)
   }
 
 }//TEST(cotmatrix_entries, simple)
+
+TEST(cotmatrix_entries, intrinsic)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  //This is a cube of dimensions 1.0x1.0x1.0
+  test_common::load_mesh("cube.obj", V, F);
+  Eigen::MatrixXd Cext,Cint;
+  // compute C extrinsically
+  igl::cotmatrix_entries(V,F,Cext);
+  // compute C intrinsically
+  Eigen::MatrixXd l;
+  igl::edge_lengths(V,F,l);
+  igl::cotmatrix_entries(l,Cint);
+  test_common::assert_near(Cext,Cint,igl::EPS<double>());
+}