intrinsic delaunay cotmatrix + test

include/igl/intrinsic_delaunay_cotmatrix.cpp

@@ -0,0 +1,50 @@
+// 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 "intrinsic_delaunay_cotmatrix.h"
+#include "edge_lengths.h"
+#include "intrinsic_delaunay_triangulation.h"
+#include "cotmatrix_intrinsic.h"
+#include <cassert>
+
+template <typename DerivedV, typename DerivedF, typename Scalar>
+IGL_INLINE void igl::intrinsic_delaunay_cotmatrix(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<Scalar>& L)
+{
+  Eigen::Matrix<Scalar, Eigen::Dynamic, 3> l_intrinsic;
+  DerivedF F_intrinsic;
+  return igl::intrinsic_delaunay_cotmatrix(V,F,L,l_intrinsic,F_intrinsic);
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename Scalar,
+  typename Derivedl_intrinsic,
+  typename DerivedF_intrinsic>
+IGL_INLINE void igl::intrinsic_delaunay_cotmatrix(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedF> & F, 
+  Eigen::SparseMatrix<Scalar>& L,
+  Eigen::PlainObjectBase<Derivedl_intrinsic> & l_intrinsic,
+  Eigen::PlainObjectBase<DerivedF_intrinsic> & F_intrinsic)
+{
+  assert(F.cols() == 3 && "Only triangles are supported");
+  Eigen::Matrix<Scalar, Eigen::Dynamic, 3> l;
+  igl::edge_lengths(V,F,l);
+  igl::intrinsic_delaunay_triangulation(l,F,l_intrinsic,F_intrinsic);
+  igl::cotmatrix_intrinsic(l_intrinsic,F_intrinsic,L);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::intrinsic_delaunay_cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+#endif

include/igl/intrinsic_delaunay_cotmatrix.h

@@ -0,0 +1,51 @@
+// 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_INTRINSIC_DELAUNAY_COTMATRIX_H
+#define IGL_INTRINSIC_DELAUNAY_COTMATRIX_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // INTRINSIC_DELAUNAY_COTMATRIX Computes the discrete cotangent Laplacian of a
+  // mesh after converting it into its intrinsic Delaunay triangulation (see,
+  // e.g., [Fisher et al. 2007].
+  //
+  // Inputs:
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by 3 list of mesh elements (triangles or tetrahedra)
+  // Outputs: 
+  //   L  #V by #V cotangent matrix, each row i corresponding to V(i,:)
+  //   l_intrinsic  #F by 3 list of intrinsic edge-lengths used to compute L
+  //   F_intrinsic  #F by 3 list of intrinsic face indices used to compute L
+  //
+  // See also: intrinsic_delaunay_triangulation, cotmatrix, cotmatrix_intrinsic
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename Scalar,
+    typename Derivedl_intrinsic,
+    typename DerivedF_intrinsic>
+  IGL_INLINE void intrinsic_delaunay_cotmatrix(
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    Eigen::SparseMatrix<Scalar>& L,
+    Eigen::PlainObjectBase<Derivedl_intrinsic> & l_intrinsic,
+    Eigen::PlainObjectBase<DerivedF_intrinsic> & F_intrinsic);
+  template <typename DerivedV, typename DerivedF, typename Scalar>
+  IGL_INLINE void intrinsic_delaunay_cotmatrix(
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedF> & F, 
+    Eigen::SparseMatrix<Scalar>& L);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "intrinsic_delaunay_cotmatrix.cpp"
+#endif
+
+#endif

tests/include/igl/intrinsic_delaunay_cotmatrix.cpp

@@ -0,0 +1,73 @@
+#include <test_common.h>
+#include <igl/intrinsic_delaunay_cotmatrix.h>
+#include <igl/EPS.h>
+#include <igl/triangulated_grid.h>
+#include <igl/is_border_vertex.h>
+
+class intrinsic_delaunay_cotmatrix : public ::testing::TestWithParam<std::string> {};
+
+TEST(intrinsic_delaunay_cotmatrix,skewed_grid)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  const int s = 7;
+  igl::triangulated_grid(s,s,V,F);
+  // Skew against diagonal direction
+  V.col(0) -= 1.1*V.col(1);
+  Eigen::SparseMatrix<double> L;
+  Eigen::MatrixXi F_intrinsic;
+  Eigen::MatrixXd l_intrinsic;
+  igl::intrinsic_delaunay_cotmatrix(V,F,L,l_intrinsic,F_intrinsic);
+  Eigen::VectorXi LI,LJ;
+  Eigen::VectorXd LV;
+  igl::find(L,LI,LJ,LV);
+  const auto is_boundary_edge = [](const int i, const int j, const int s)->bool
+  {
+    const auto is_boundary_vertex = [](const int i, const int s)->bool
+    {
+      return (i<s) || (i>=s*s-s) || (i%s == 0) || (i%s == s-1);
+    };
+    return is_boundary_vertex(i,s) && is_boundary_vertex(j,s);
+  };
+  // Off diagonals should be all non-positive
+  for(int k = 0;k<LI.size();k++)
+  {
+    if(LI(k) != LJ(k) && !is_boundary_edge(LI(k),LJ(k),s))
+    {
+      ASSERT_GT(LV(k),-igl::EPS<double>());
+    }
+  }
+}
+
+TEST_P(intrinsic_delaunay_cotmatrix,manifold_meshes)
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  test_common::load_mesh(GetParam(), V, F);
+  Eigen::SparseMatrix<double> L;
+  Eigen::MatrixXi F_intrinsic;
+  Eigen::MatrixXd l_intrinsic;
+  igl::intrinsic_delaunay_cotmatrix(V,F,L,l_intrinsic,F_intrinsic);
+  Eigen::VectorXi LI,LJ;
+  Eigen::VectorXd LV;
+  igl::find(L,LI,LJ,LV);
+
+  const std::vector<bool> is_boundary_vertex = igl::is_border_vertex(F);
+  // Off diagonals should be all non-positive
+  for(int k = 0;k<LI.size();k++)
+  {
+    if(LI(k) != LJ(k) && 
+      !(is_boundary_vertex[LI(k)] && is_boundary_vertex[LJ(k)]))
+    {
+      ASSERT_GT(LV(k),-igl::EPS<double>());
+    }
+  }
+}
+
+INSTANTIATE_TEST_CASE_P
+(
+ manifold_meshes,
+ intrinsic_delaunay_cotmatrix,
+ ::testing::ValuesIn(test_common::manifold_meshes()),
+ test_common::string_test_name
+);