spectra module (#2216)

* lscm hessian and spectral

* spectra module + test

* try to use Eigen3_FOUND (not working locally)

* use fork

* Don't use size_t for small int (windows hell)

.gitignore

@@ -47,3 +47,4 @@ LibiglOptions.cmake
 
 # macos debris
 .DS_Store
+*~

CMakeLists.txt

@@ -93,6 +93,7 @@ option(LIBIGL_IMGUI            "Build target igl::imgui"            ${LIBIGL_TOP
 option(LIBIGL_OPENGL           "Build target igl::opengl"           ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_PNG              "Build target igl::png"              ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_PREDICATES       "Build target igl::predicates"       ${LIBIGL_TOPLEVEL_PROJECT})
+option(LIBIGL_SPECTRA          "Build target igl::spectra"          ${LIBIGL_TOPLEVEL_PROJECT})
 option(LIBIGL_XML              "Build target igl::xml"              ${LIBIGL_TOPLEVEL_PROJECT})
 
 # Copyleft modules. These modules are available under GPL license, and their dependencies are

cmake/igl/modules/spectra.cmake

@@ -0,0 +1,25 @@
+# 1. Define module
+igl_add_library(igl_spectra)
+
+# 2. Include headers
+include(GNUInstallDirs)
+target_include_directories(igl_spectra ${IGL_SCOPE}
+    $<BUILD_INTERFACE:${libigl_SOURCE_DIR}/include>
+    $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>
+)
+
+# 3. Target sources
+file(GLOB INC_FILES "${libigl_SOURCE_DIR}/include/igl/spectra/*.h")
+file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/include/igl/spectra/*.cpp")
+igl_target_sources(igl_spectra ${INC_FILES} ${SRC_FILES})
+
+# 4. Dependencies
+include(spectra)
+target_link_libraries(igl_spectra ${IGL_SCOPE}
+    igl::core
+    spectra::spectra
+)
+
+# 5. Unit tests
+file(GLOB SRC_FILES "${libigl_SOURCE_DIR}/tests/include/igl/spectra/*.cpp")
+igl_add_test(igl_spectra ${SRC_FILES})

cmake/libigl.cmake

@@ -19,6 +19,7 @@ igl_include_optional(glfw)
 igl_include_optional(imgui)
 igl_include_optional(predicates)
 igl_include_optional(png)
+igl_include_optional(spectra)
 igl_include_optional(xml)
 
 # Libigl copyleft modules

cmake/recipes/external/spectra.cmake

@@ -0,0 +1,16 @@
+if(TARGET spectra::spectra)
+  return()
+endif()
+include(FetchContent)
+
+message(STATUS "Third-party: creating target 'spectra::spectra'")
+
+# Use fork because yixuan/spectra struggles to find Eigen3
+FetchContent_Declare(
+  Spectra
+  GIT_REPOSITORY https://github.com/alecjacobson/spectra/
+  GIT_TAG bbdc521b70a733c52ebfc0ac1484c82e13c3d140
+)
+FetchContent_MakeAvailable(Spectra)
+
+add_library(spectra::spectra ALIAS Spectra)

include/igl/spectra/eigs.cpp

@@ -0,0 +1,139 @@
+#include "eigs.h"
+#include <igl/sort.h>
+#include <igl/slice.h>
+#include <Spectra/SymGEigsShiftSolver.h>
+#include <cassert>
+
+template <
+  typename EigsScalar,
+  typename DerivedU,
+  typename DerivedS,
+  typename Solver>
+IGL_INLINE bool igl::spectra::eigs(
+  const Eigen::SparseMatrix<EigsScalar> & A,
+  const Eigen::SparseMatrix<EigsScalar> & B,
+  const int k,
+  const igl::EigsType type,
+  Eigen::PlainObjectBase<DerivedU> & U,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  assert(k > 0 && "k should be positive");
+  assert(k < A.rows() && "k should be less than size of A");
+  assert(type == igl::EIGS_TYPE_SM && "Only SM supported");
+  // This seems like a hack. For the "eigs: grid" test this is necessary to get
+  // at least 1e-4 error for the first 5 eigen values. It's annoying that this
+  // means that the zero modes become O(sigma) and this is now rather large.
+  //
+  // I wonder if this is an issue with SparseLU and if UMFPACK would be better.
+  //
+  // Ideally this value would be 0.
+  const EigsScalar sigma = 1e-8;
+  return igl::spectra::eigs(A,B,k,sigma,U,S);
+}
+template <
+  typename EigsScalar,
+  typename DerivedU,
+  typename DerivedS,
+  typename Solver>
+IGL_INLINE bool igl::spectra::eigs(
+  const Eigen::SparseMatrix<EigsScalar> & A,
+  const Eigen::SparseMatrix<EigsScalar> & B,
+  const int k,
+  const EigsScalar sigma,
+  Eigen::PlainObjectBase<DerivedU> & U,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+
+  assert(k > 0 && "k should be positive");
+  assert(k < A.rows() && "k should be less than size of A");
+
+  class SparseMatProd
+  {
+    public:
+      using Scalar = EigsScalar;
+      const Eigen::SparseMatrix<Scalar> & m_B;
+      SparseMatProd(const Eigen::SparseMatrix<Scalar> & B) : m_B(B) {}
+      int rows() const { return m_B.rows(); }
+      int cols() const { return m_B.cols(); }
+      void perform_op(const Scalar *x_in, Scalar *y_out) const
+      {
+        typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> VectorXS;
+        Eigen::Map<const VectorXS> x(x_in, m_B.cols());
+        Eigen::Map<      VectorXS> y(y_out, m_B.rows());
+        y = m_B * x;
+      }
+  };
+
+  // Solver must expose .compute(A) and .solve(x)
+  class ShiftInvert
+  {
+  public:
+    using Scalar = EigsScalar;
+  private:
+    const Eigen::SparseMatrix<Scalar> & m_A;
+    const Eigen::SparseMatrix<Scalar> & m_B;
+    Scalar m_sigma;
+    Solver m_solver;
+  public:
+    bool m_solver_is_successfully_factorized;
+    ShiftInvert(
+        const Eigen::SparseMatrix<Scalar>& A, 
+        const Eigen::SparseMatrix<Scalar>& B, 
+        const Scalar sigma):
+        m_A(A), m_B(B)
+    {
+      assert(m_A.rows() == m_A.cols() && "A must be square");
+      assert(m_B.rows() == m_B.cols() && "B must be square");
+      assert(m_A.rows() == m_B.cols() && "A and B must have the same size");
+      set_shift(sigma, true);
+    }
+    void set_shift(const Scalar & sigma, const bool force = false)
+    {
+      if(sigma == m_sigma && !force)
+      {
+        return;
+      }
+      m_sigma = sigma;
+      const Eigen::SparseMatrix<Scalar> C = m_A + m_sigma * m_B;
+      m_solver.compute(C);
+      m_solver_is_successfully_factorized = (m_solver.info() == Eigen::Success);
+    }
+    int rows() const { return m_A.rows(); }
+    int cols() const { return m_A.cols(); }
+    void perform_op(const Scalar* x_in,Scalar* y_out) const
+    {
+      typedef Eigen::Matrix<Scalar, Eigen::Dynamic, 1> VectorXS;
+      Eigen::Map<const VectorXS>x(x_in, m_A.cols());
+      Eigen::Map<VectorXS>y(y_out, m_A.rows());
+      y = m_solver.solve(x);
+    }
+  };
+
+  SparseMatProd Bop(B);
+  ShiftInvert op(A, B, sigma);
+  if(!op.m_solver_is_successfully_factorized)
+  {
+    return false;
+  }
+  Spectra::SymGEigsShiftSolver<ShiftInvert, SparseMatProd, Spectra::GEigsMode::ShiftInvert> geigs(op, Bop, k, 2*k, sigma);
+
+  geigs.init();
+  int nconv = geigs.compute(Spectra::SortRule::LargestMagn);
+  if (geigs.info() != Spectra::CompInfo::Successful)
+  {
+    return false;
+  }
+  U = geigs.eigenvectors().template cast<typename DerivedU::Scalar>();
+  S = geigs.eigenvalues().template cast<typename DerivedS::Scalar>();
+
+  Eigen::VectorXi I;
+  igl::sort( Eigen::VectorXd(S), 1, false, S, I);
+  igl::slice(Eigen::MatrixXd(U),I,2,U);
+  return true;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template bool igl::spectra::eigs<double, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::SparseLU<Eigen::SparseMatrix<double, 0, int>, Eigen::COLAMDOrdering<int> > >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int> const&, int , igl::EigsType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/spectra/eigs.h

@@ -0,0 +1,72 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2023 Alec Jacobson <[email protected]>
+// 
+// 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_SPECTRA_EIGS_H
+#define IGL_SPECTRA_EIGS_H
+#include "../igl_inline.h"
+#include "../eigs.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+#include <Eigen/SparseLU>
+
+namespace igl
+{
+  namespace spectra
+  {
+    // Act like MATLAB's eigs function. Compute the first/last k eigen pairs of
+    // the generalized eigen value problem:
+    //
+    //     A u = s B u
+    //
+    // Solutions are approximate and sorted. 
+    //
+    // Ideally one should use ARPACK and the Eigen unsupported ARPACK module.
+    // This implementation does simple, naive power iterations.
+    //
+    // Inputs:
+    //   A  #A by #A symmetric matrix
+    //   B  #A by #A symmetric positive-definite matrix
+    //   k  number of eigen pairs to compute
+    //   type  whether to extract from the high or low end
+    // Outputs:
+    //   sU  #A by k list of sorted eigen vectors (descending)
+    //   sS  k list of sorted eigen values (descending)
+    //
+    // Known issues:
+    //   - only the 'sm' small magnitude eigen values are well supported
+    //   
+    template <
+      typename EigsScalar,
+      typename DerivedU,
+      typename DerivedS,
+      typename Solver = Eigen::SparseLU<Eigen::SparseMatrix<EigsScalar>> >
+    IGL_INLINE bool eigs(
+      const Eigen::SparseMatrix<EigsScalar> & A,
+      const Eigen::SparseMatrix<EigsScalar> & B,
+      const int k,
+      const igl::EigsType type,
+      Eigen::PlainObjectBase<DerivedU> & U,
+      Eigen::PlainObjectBase<DerivedS> & S);
+    template <
+      typename EigsScalar,
+      typename DerivedU,
+      typename DerivedS,
+      typename Solver = Eigen::SparseLU<Eigen::SparseMatrix<EigsScalar>> >
+    IGL_INLINE bool eigs(
+      const Eigen::SparseMatrix<EigsScalar> & A,
+      const Eigen::SparseMatrix<EigsScalar> & B,
+      const int k,
+      const EigsScalar sigma,
+      Eigen::PlainObjectBase<DerivedU> & U,
+      Eigen::PlainObjectBase<DerivedS> & S);
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "eigs.cpp"
+#endif
+#endif

include/igl/spectra/lscm.cpp

@@ -0,0 +1,29 @@
+#include "lscm.h"
+#include "../lscm_hessian.h"
+#include "../massmatrix.h"
+#include "../repdiag.h"
+#include "eigs.h"
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedV_uv>
+IGL_INLINE bool igl::spectra::lscm(
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  Eigen::PlainObjectBase<DerivedV_uv> & V_uv)
+{
+  using Scalar = typename DerivedV_uv::Scalar;
+  Eigen::SparseMatrix<Scalar> Q;
+  igl::lscm_hessian(V,F,Q);
+  Eigen::SparseMatrix<Scalar> M;
+  igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,M);
+  Eigen::SparseMatrix<Scalar> M2;
+  igl::repdiag(M,2,M2);
+
+  Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>  U;
+  Eigen::Matrix<Scalar, Eigen::Dynamic, 1>  S;
+  igl::spectra::eigs(Q,M2,3,igl::EIGS_TYPE_SM,U,S);
+
+  V_uv.resize(V.rows(),2);
+  V_uv<< U.col(0).head(V.rows()),U.col(0).tail(V.rows());
+}

include/igl/spectra/lscm.h

@@ -0,0 +1,52 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2023 Alec Jacobson
+//
+// 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_SPECTRA_LSCM_H
+#define IGL_SPECTRA_LSCM_H
+#include "../igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+namespace igl 
+{
+namespace spectra
+{
+  // Compute a Least-squares conformal map parametrization (equivalently
+  // derived in "Intrinsic Parameterizations of Surface Meshes" [Desbrun et al.
+  // 2002] and "Least Squares Conformal Maps for Automatic Texture Atlas
+  // Generation" [Lévy et al. 2002]), though this implementation follows the
+  // derivation in: "Spectral Conformal Parameterization" [Mullen et al. 2008]
+  // Free boundary version. "Spectral Conformal Parameterization" using Eigen
+  // decomposition. Assumes mesh is a single connected component topologically
+  // equivalent to a chunk of the plane.
+  //
+  //
+  // Inputs:
+  //   V  #V by 3 list of mesh vertex positions
+  //   F  #F by 3 list of mesh faces (must be triangles)
+  // Outputs:
+  //   UV #V by 2 list of 2D mesh vertex positions in UV space
+  // Returns true only on solver success.
+  //
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedV_uv>
+  IGL_INLINE bool lscm(
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    Eigen::PlainObjectBase<DerivedV_uv> & V_uv);
+}
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "lscm.cpp"
+#endif
+
+#endif
+

tests/include/igl/spectra/eigs.cpp

@@ -0,0 +1,34 @@
+#include <test_common.h>
+#include <igl/spectra/eigs.h>
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/triangulated_grid.h>
+
+#include <igl/matlab_format.h>
+#include <iostream>
+
+TEST_CASE("eigs: grid", "[igl/spectra]")
+{
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi F;
+  igl::triangulated_grid(10,10,V,F);
+  Eigen::SparseMatrix<double> L;
+  igl::cotmatrix(V,F,L);
+  Eigen::SparseMatrix<double> M;
+  igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,M);
+  Eigen::MatrixXd U;
+  Eigen::VectorXd S;
+  Eigen::VectorXd S_matlab(5,1);
+  S_matlab<<
+    -0.00000000000000600,
+    -9.76979543268284445,
+    -9.76979543268286221,
+    -19.53959086536570311,
+    -37.90080021472548566;
+  const bool success = igl::spectra::eigs(L,M,S_matlab.size(),igl::EIGS_TYPE_SM,U,S);
+  REQUIRE(success);
+
+
+  test_common::assert_near(S,S_matlab,1e-4);
+}
+