cotmatrix computed intrinsically (using edge lengths) + templates + tests

include/igl/cotmatrix_entries.cpp

@@ -135,6 +135,10 @@ IGL_INLINE void igl::cotmatrix_entries(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::cotmatrix_entries<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -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<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 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> >&);
 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> >&);

include/igl/cotmatrix_intrinsic.cpp

@@ -0,0 +1,67 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "cotmatrix_intrinsic.h"
+#include "cotmatrix_entries.h"
+#include <iostream>
+
+template <typename Derivedl, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::cotmatrix_intrinsic(
+  const Eigen::MatrixBase<Derivedl> & l, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<Scalar>& L)
+{
+  using namespace Eigen;
+  using namespace std;
+  // Cribbed from cotmatrix
+
+  const int nverts = F.maxCoeff()+1;
+  L.resize(nverts,nverts);
+  Matrix<int,Dynamic,2> edges;
+  int simplex_size = F.cols();
+  // 3 for triangles, 4 for tets
+  assert(simplex_size == 3);
+  // This is important! it could decrease the comptuation time by a factor of 2
+  // Laplacian for a closed 2d manifold mesh will have on average 7 entries per
+  // row
+  L.reserve(10*nverts);
+  edges.resize(3,2);
+  edges << 
+    1,2,
+    2,0,
+    0,1;
+  // Gather cotangents
+  Matrix<Scalar,Dynamic,Dynamic> C;
+  cotmatrix_entries(l,C);
+  
+  vector<Triplet<Scalar> > IJV;
+  IJV.reserve(F.rows()*edges.rows()*4);
+  // Loop over triangles
+  for(int i = 0; i < F.rows(); i++)
+  {
+    // loop over edges of element
+    for(int e = 0;e<edges.rows();e++)
+    {
+      int source = F(i,edges(e,0));
+      int dest = F(i,edges(e,1));
+      IJV.push_back(Triplet<Scalar>(source,dest,C(i,e)));
+      IJV.push_back(Triplet<Scalar>(dest,source,C(i,e)));
+      IJV.push_back(Triplet<Scalar>(source,source,-C(i,e)));
+      IJV.push_back(Triplet<Scalar>(dest,dest,-C(i,e)));
+    }
+  }
+  L.setFromTriplets(IJV.begin(),IJV.end());
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::cotmatrix_intrinsic<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cotmatrix_intrinsic<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cotmatrix_intrinsic<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cotmatrix_intrinsic<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>&);
+#endif
+

include/igl/cotmatrix_intrinsic.h

@@ -0,0 +1,40 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// 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 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_COTMATRIX_INTRINSIC_H
+#define IGL_COTMATRIX_INTRINSIC_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl 
+{
+  // Constructs the cotangent stiffness matrix (discrete laplacian) for a given
+  // mesh with faces F and edge lengths l.
+  //
+  // Inputs:
+  //   l  #F by 3 list of (half-)edge lengths
+  //   F  #F by 3 list of face indices into some (not necessarily
+  //     determined/embedable) list of vertex positions V. It is assumed #V ==
+  //     F.maxCoeff()+1
+  // Outputs:
+  //   L  #V by #V sparse Laplacian matrix
+  //
+  // See also: cotmatrix, intrinsic_delaunay_cotmatrix
+  template <typename Derivedl, typename DerivedF, typename Scalar>
+  IGL_INLINE void cotmatrix_intrinsic(
+    const Eigen::MatrixBase<Derivedl> & l, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    Eigen::SparseMatrix<Scalar>& L);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "cotmatrix.cpp"
+#endif
+
+#endif

include/igl/doublearea.cpp

@@ -232,6 +232,8 @@ IGL_INLINE void igl::doublearea_quad(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::doublearea<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 1, -1, -1>::Scalar, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<float, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 0, -1, 3> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<float, -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<float, -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> >&);

include/igl/sort.cpp

@@ -317,6 +317,8 @@ if(!ascending)
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<double, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, -1, 1, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 // generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<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> >&);

tests/include/igl/cotmatrix_intrinsic.cpp

@@ -0,0 +1,90 @@
+#include <test_common.h>
+#include <igl/cotmatrix_intrinsic.h>
+#include <igl/cotmatrix.h>
+#include <igl/edge_lengths.h>
+#include <igl/matlab_format.h>
+#include <igl/EPS.h>
+
+class cotmatrix_intrinsic : public ::testing::TestWithParam<std::string> {};
+
+TEST(cotmatrix_intrinsic, periodic)
+{
+  const Eigen::MatrixXi F = (Eigen::MatrixXi(18,3)<<
+    0,3,1,
+    3,4,1,
+    1,4,2,
+    4,5,2,
+    2,5,0,
+    5,3,0,
+    3,6,4,
+    6,7,4,
+    4,7,5,
+    7,8,5,
+    5,8,3,
+    8,6,3,
+    6,0,7,
+    0,1,7,
+    7,1,8,
+    1,2,8,
+    8,2,6,
+    2,0,6).finished();
+  const Eigen::MatrixXd l = (Eigen::MatrixXd(18,3)<<
+    0.47140452079103168,0.33333333333333331,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333331,
+    0.47140452079103168,0.33333333333333331,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333331,
+    0.47140452079103168,0.33333333333333337,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333337,
+    0.47140452079103168,0.33333333333333331,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333331,
+    0.47140452079103168,0.33333333333333331,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333331,
+    0.47140452079103168,0.33333333333333337,0.33333333333333331,
+    0.33333333333333331,0.47140452079103168,0.33333333333333337,
+    0.47140452079103168,0.33333333333333331,0.33333333333333337,
+    0.33333333333333337,0.47140452079103168,0.33333333333333331,
+    0.47140452079103168,0.33333333333333331,0.33333333333333337,
+    0.33333333333333337,0.47140452079103168,0.33333333333333331,
+    0.47140452079103173,0.33333333333333337,0.33333333333333337,
+    0.33333333333333337,0.47140452079103173,0.33333333333333337).finished();
+  Eigen::SparseMatrix<double> L;
+  igl::cotmatrix_intrinsic(l,F,L);
+  const Eigen::MatrixXd L_d = L;
+  const Eigen::MatrixXd L_gt = (Eigen::MatrixXd(9,9)<<
+    -4,1,1,1,0,0,1,0,0,
+    1,-4,1,0,1,0,0,1,0,
+    1,1,-4,0,0,1,0,0,1,
+    1,0,0,-4,1,1,1,0,0,
+    0,1,0,1,-4,1,0,1,0,
+    0,0,1,1,1,-4,0,0,1,
+    1,0,0,1,0,0,-4,1,1,
+    0,1,0,0,1,0,1,-4,1,
+    0,0,1,0,0,1,1,1,-4).finished();
+  test_common::assert_near(L_d,L_gt,igl::EPS<double>());
+}
+
+TEST_P(cotmatrix_intrinsic , manifold_meshes)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  test_common::load_mesh(GetParam(), V, F);
+  Eigen::MatrixXd l;
+  igl::edge_lengths(V,F,l);
+  Eigen::SparseMatrix<double> L,Li;
+  igl::cotmatrix(V,F,L);
+  igl::cotmatrix_intrinsic(l,F,Li);
+  // Augh, we don't have assert_near for sparse matrices...
+  // Instead test that bilinear form is near equal for 2 random vectors
+  const Eigen::VectorXd u = Eigen::VectorXd::Random(V.rows(),1);
+  const Eigen::VectorXd v = Eigen::VectorXd::Random(V.rows(),1);
+  const double uv = u.norm()*v.norm();
+  ASSERT_NEAR( u.dot(L*v)/uv, u.dot(Li*v)/uv, igl::EPS<double>());
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+ manifold_meshes,
+ cotmatrix_intrinsic,
+ ::testing::ValuesIn(test_common::manifold_meshes()),
+ test_common::string_test_name
+);

tests/include/igl/doublearea.cpp

@@ -32,7 +32,7 @@ TEST_P(doublearea, VF_vs_ABC )
 
 INSTANTIATE_TEST_CASE_P
 (
- all_meshes,
+ VF_vs_ABC,
  doublearea,
  ::testing::ValuesIn(test_common::all_meshes()),
  test_common::string_test_name