Adding 2D version of SCAF for Bijective Maps and a tutorial (#752)

* Add scaf code

Add the tutorial entry

Tutorial entry for SCAF. still need polishing

Put -4 in the energy computation of scaf:SymmD

* :school: Insert SCAF to 711

Fixes in tutorial.md

updating tutorial.html

Fix for gcc 4.8

Remove 4.1

Refactoring scaf WIP

Further refactoring WIP

Code Cleaning

* 🔨 Improve SLIM interface

* :recycle: refactor SCAF with less lines

* :recycle: refactor harmonic with boundary

* :sparkles: Add direct signature for scaf

* :memo: Change tutorial order wip

* :school: Timer in SLIM tutorial makes more sense

* :sparkles: Add direct interface for SLIM

* :art: format SCAF code.

* :snake: Python binding for slim and scaf

* :snake: expose SLIMData (maybe elevate its status more in the future)

* :alien: :snake: deprecate PYBIND11_PLUGIN

* :sparkles: harmonic with holes

* :recycle: SCAF, Decouple Initialization to tutorial

* :school: SCAF tutorial improvement

* :school: Tutorial order; insert scaf

* :art: Removing Adjusted_grad

*     :school: Inserted SCAF, not messing around with tutorials order until the new page is online

* :green_heart: Fix for static build

* :recycle: Uniform SLIM and SCAF and  python binding canonical

* :bug: Change in ScafData

* :apple: Fix namespace for build

* :sparkles: Introduced Mapping Energy Type to replace adhoc SLIMEnergy

* :recycle: Refactor subroutines from SLIM

* Revert python changes

* leave fill hole outside harmonic

* Revert Python changes

* hole fill

* Minor Fixes for SCAF

SCAF->IGL_SCAF
Topological_hole_fill 
"<igl/XXX>" in slim tutorial


Former-commit-id: 4edf5dede5a4430b86d4dc65705c59e643d955df

include/igl/MappingEnergyType.h

@@ -0,0 +1,26 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Zhongshi Jiang <jiangzs@nyu.edu>
+// 
+// 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_MAPPINGENERGYTYPE_H
+#define IGL_MAPPINGENERGYTYPE_H
+namespace igl
+{
+  // Energy Types used for Parameterization/Mapping. 
+  // Refer to SLIM [Rabinovich et al. 2017] for more details
+  // Todo: Integrate with ARAPEnergyType
+  
+  enum MappingEnergyType
+  {
+    ARAP,
+    LOG_ARAP,
+    SYMMETRIC_DIRICHLET,
+    CONFORMAL,
+    EXP_CONFORMAL,
+    EXP_SYMMETRIC_DIRICHLET
+  };
+}
+#endif

include/igl/mapping_energy_with_jacobians.cpp

@@ -0,0 +1,141 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Zhongshi Jiang <jiangzs@nyu.edu>
+// 
+// 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 "mapping_energy_with_jacobians.h"
+#include "polar_svd.h"
+
+IGL_INLINE double igl::mapping_energy_with_jacobians(
+  const Eigen::MatrixXd &Ji, 
+  const Eigen::VectorXd &areas, 
+  igl::MappingEnergyType slim_energy, 
+  double exp_factor){
+
+  double energy = 0;
+  if (Ji.cols() == 4)
+  {
+    Eigen::Matrix<double, 2, 2> ji;
+    for (int i = 0; i < Ji.rows(); i++)
+    {
+      ji(0, 0) = Ji(i, 0);
+      ji(0, 1) = Ji(i, 1);
+      ji(1, 0) = Ji(i, 2);
+      ji(1, 1) = Ji(i, 3);
+
+      typedef Eigen::Matrix<double, 2, 2> Mat2;
+      typedef Eigen::Matrix<double, 2, 1> Vec2;
+      Mat2 ri, ti, ui, vi;
+      Vec2 sing;
+      igl::polar_svd(ji, ri, ti, ui, sing, vi);
+      double s1 = sing(0);
+      double s2 = sing(1);
+
+      switch (slim_energy)
+      {
+        case igl::MappingEnergyType::ARAP:
+        {
+          energy += areas(i) * (pow(s1 - 1, 2) + pow(s2 - 1, 2));
+          break;
+        }
+        case igl::MappingEnergyType::SYMMETRIC_DIRICHLET:
+        {
+          energy += areas(i) * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2));
+          break;
+        }
+        case igl::MappingEnergyType::EXP_SYMMETRIC_DIRICHLET:
+        {
+          energy += areas(i) * exp(exp_factor * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2)));
+          break;
+        }
+        case igl::MappingEnergyType::LOG_ARAP:
+        {
+          energy += areas(i) * (pow(log(s1), 2) + pow(log(s2), 2));
+          break;
+        }
+        case igl::MappingEnergyType::CONFORMAL:
+        {
+          energy += areas(i) * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2));
+          break;
+        }
+        case igl::MappingEnergyType::EXP_CONFORMAL:
+        {
+          energy += areas(i) * exp(exp_factor * ((pow(s1, 2) + pow(s2, 2)) / (2 * s1 * s2)));
+          break;
+        }
+
+      }
+
+    }
+  }
+  else
+  {
+    Eigen::Matrix<double, 3, 3> ji;
+    for (int i = 0; i < Ji.rows(); i++)
+    {
+      ji(0, 0) = Ji(i, 0);
+      ji(0, 1) = Ji(i, 1);
+      ji(0, 2) = Ji(i, 2);
+      ji(1, 0) = Ji(i, 3);
+      ji(1, 1) = Ji(i, 4);
+      ji(1, 2) = Ji(i, 5);
+      ji(2, 0) = Ji(i, 6);
+      ji(2, 1) = Ji(i, 7);
+      ji(2, 2) = Ji(i, 8);
+
+      typedef Eigen::Matrix<double, 3, 3> Mat3;
+      typedef Eigen::Matrix<double, 3, 1> Vec3;
+      Mat3 ri, ti, ui, vi;
+      Vec3 sing;
+      igl::polar_svd(ji, ri, ti, ui, sing, vi);
+      double s1 = sing(0);
+      double s2 = sing(1);
+      double s3 = sing(2);
+
+      switch (slim_energy)
+      {
+        case igl::MappingEnergyType::ARAP:
+        {
+          energy += areas(i) * (pow(s1 - 1, 2) + pow(s2 - 1, 2) + pow(s3 - 1, 2));
+          break;
+        }
+        case igl::MappingEnergyType::SYMMETRIC_DIRICHLET:
+        {
+          energy += areas(i) * (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2) + pow(s3, 2) + pow(s3, -2));
+          break;
+        }
+        case igl::MappingEnergyType::EXP_SYMMETRIC_DIRICHLET:
+        {
+          energy += areas(i) * exp(exp_factor *
+                                    (pow(s1, 2) + pow(s1, -2) + pow(s2, 2) + pow(s2, -2) + pow(s3, 2) + pow(s3, -2)));
+          break;
+        }
+        case igl::MappingEnergyType::LOG_ARAP:
+        {
+          energy += areas(i) * (pow(log(s1), 2) + pow(log(std::abs(s2)), 2) + pow(log(std::abs(s3)), 2));
+          break;
+        }
+        case igl::MappingEnergyType::CONFORMAL:
+        {
+          energy += areas(i) * ((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow(s1 * s2 * s3, 2. / 3.)));
+          break;
+        }
+        case igl::MappingEnergyType::EXP_CONFORMAL:
+        {
+          energy += areas(i) * exp((pow(s1, 2) + pow(s2, 2) + pow(s3, 2)) / (3 * pow(s1 * s2 * s3, 2. / 3.)));
+          break;
+        }
+      }
+    }
+  }
+
+  return energy;
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+#endif

include/igl/mapping_energy_with_jacobians.h

@@ -0,0 +1,36 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2018 Zhongshi Jiang <jiangzs@nyu.edu>
+// 
+// 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_MAPPING_ENERGY_WITH_JACOBIANS_H
+#define IGL_MAPPING_ENERGY_WITH_JACOBIANS_H
+
+#include "igl_inline.h"
+#include <Eigen/Dense>
+#include "MappingEnergyType.h"
+
+namespace igl
+{
+   // compute the rotation-invariant energy of a mapping (represented in Jacobians and areas)
+   // Input:
+   // Ji: #F by 4 (9 if 3D) entries of jacobians
+   // areas: #F by 1 face areas
+   // slim_energy: energy type as in igl::MappingEnergyType
+   // exp_factor: see igl::MappingEnergyType
+   //
+   // Output:
+   // energy value
+   IGL_INLINE double mapping_energy_with_jacobians(const Eigen::MatrixXd &Ji, 
+                                                  const Eigen::VectorXd &areas, 
+                                                  igl::MappingEnergyType slim_energy, 
+                                                  double exp_factor);
+  
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "mapping_energy_with_jacobians.cpp"
+#endif
+
+#endif

include/igl/scaf.cpp

@@ -0,0 +1,707 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Zhongshi Jiang <jiangzs@nyu.edu>
+//
+// 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 "scaf.h"
+
+#include <Eigen/Dense>
+#include <Eigen/IterativeLinearSolvers>
+#include <Eigen/Sparse>
+#include <Eigen/SparseCholesky>
+#include <Eigen/SparseQR>
+#include <igl/PI.h>
+#include <igl/Timer.h>
+#include <igl/boundary_loop.h>
+#include <igl/cat.h>
+#include <igl/doublearea.h>
+#include <igl/flip_avoiding_line_search.h>
+#include <igl/flipped_triangles.h>
+#include <igl/grad.h>
+#include <igl/harmonic.h>
+#include <igl/local_basis.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/polar_svd.h>
+#include <igl/readOBJ.h>
+#include <igl/slice.h>
+#include <igl/slice_into.h>
+#include <igl/slim.h>
+#include <igl/triangle/triangulate.h>
+#include "mapping_energy_with_jacobians.h"
+
+#include <iostream>
+#include <map>
+#include <algorithm>
+#include <set>
+#include <vector>
+namespace igl
+{
+
+namespace scaf
+{
+void update_scaffold(igl::SCAFData &s)
+{
+  s.mv_num = s.m_V.rows();
+  s.mf_num = s.m_T.rows();
+
+  s.v_num = s.w_uv.rows();
+  s.sf_num = s.s_T.rows();
+
+  s.sv_num = s.v_num - s.mv_num;
+  s.f_num = s.sf_num + s.mf_num;
+
+  s.s_M = Eigen::VectorXd::Constant(s.sf_num, s.scaffold_factor);
+}
+
+void adjusted_grad(Eigen::MatrixXd &V,
+                   Eigen::MatrixXi &F,
+                   double area_threshold,
+                   Eigen::SparseMatrix<double> &Dx,
+                   Eigen::SparseMatrix<double> &Dy,
+                   Eigen::SparseMatrix<double> &Dz)
+{
+  Eigen::VectorXd M;
+  igl::doublearea(V, F, M);
+  std::vector<int> degen;
+  for (int i = 0; i < M.size(); i++)
+    if (M(i) < area_threshold)
+      degen.push_back(i);
+
+  Eigen::SparseMatrix<double> G;
+  igl::grad(V, F, G);
+
+  Dx = G.topRows(F.rows());
+  Dy = G.block(F.rows(), 0, F.rows(), V.rows());
+  Dz = G.bottomRows(F.rows());
+
+  // handcraft uniform gradient for faces area falling below threshold.
+  double sin60 = std::sin(M_PI / 3);
+  double cos60 = std::cos(M_PI / 3);
+  double deno = std::sqrt(sin60 * area_threshold);
+  Eigen::MatrixXd standard_grad(3, 3);
+  standard_grad << -sin60 / deno, sin60 / deno, 0,
+      -cos60 / deno, -cos60 / deno, 1 / deno,
+      0, 0, 0;
+
+  for (auto k : degen)
+    for (int j = 0; j < 3; j++)
+    {
+      Dx.coeffRef(k, F(k, j)) = standard_grad(0, j);
+      Dy.coeffRef(k, F(k, j)) = standard_grad(1, j);
+      Dz.coeffRef(k, F(k, j)) = standard_grad(2, j);
+    }
+}
+
+void compute_scaffold_gradient_matrix(SCAFData &s,
+                                      Eigen::SparseMatrix<double> &D1,
+                                      Eigen::SparseMatrix<double> &D2)
+{
+  using namespace Eigen;
+  Eigen::SparseMatrix<double> G;
+  MatrixXi F_s = s.s_T;
+  int vn = s.v_num;
+  MatrixXd V = MatrixXd::Zero(vn, 3);
+  V.leftCols(2) = s.w_uv;
+
+  double min_bnd_edge_len = INFINITY;
+  int acc_bnd = 0;
+  for (int i = 0; i < s.bnd_sizes.size(); i++)
+  {
+    int current_size = s.bnd_sizes[i];
+
+    for (int e = acc_bnd; e < acc_bnd + current_size - 1; e++)
+    {
+      min_bnd_edge_len = (std::min)(min_bnd_edge_len,
+                                    (s.w_uv.row(s.internal_bnd(e)) -
+                                     s.w_uv.row(s.internal_bnd(e + 1)))
+                                        .squaredNorm());
+    }
+    min_bnd_edge_len = (std::min)(min_bnd_edge_len,
+                                  (s.w_uv.row(s.internal_bnd(acc_bnd)) -
+                                   s.w_uv.row(s.internal_bnd(acc_bnd + current_size - 1)))
+                                      .squaredNorm());
+    acc_bnd += current_size;
+  }
+
+  double area_threshold = min_bnd_edge_len / 4.0;
+  Eigen::SparseMatrix<double> Dx, Dy, Dz;
+  adjusted_grad(V, F_s, area_threshold, Dx, Dy, Dz);
+
+  MatrixXd F1, F2, F3;
+  igl::local_basis(V, F_s, F1, F2, F3);
+  D1 = F1.col(0).asDiagonal() * Dx + F1.col(1).asDiagonal() * Dy +
+       F1.col(2).asDiagonal() * Dz;
+  D2 = F2.col(0).asDiagonal() * Dx + F2.col(1).asDiagonal() * Dy +
+       F2.col(2).asDiagonal() * Dz;
+}
+
+void mesh_improve(igl::SCAFData &s)
+{
+  using namespace Eigen;
+  MatrixXd m_uv = s.w_uv.topRows(s.mv_num);
+  MatrixXd V_bnd;
+  V_bnd.resize(s.internal_bnd.size(), 2);
+  for (int i = 0; i < s.internal_bnd.size(); i++) // redoing step 1.
+  {
+    V_bnd.row(i) = m_uv.row(s.internal_bnd(i));
+  }
+
+  if (s.rect_frame_V.size() == 0)
+  {
+    Matrix2d ob; // = rect_corners;
+    {
+      VectorXd uv_max = m_uv.colwise().maxCoeff();
+      VectorXd uv_min = m_uv.colwise().minCoeff();
+      VectorXd uv_mid = (uv_max + uv_min) / 2.;
+
+      Eigen::Array2d scaf_range(3, 3);
+      ob.row(0) = uv_mid.array() + scaf_range * ((uv_min - uv_mid).array());
+      ob.row(1) = uv_mid.array() + scaf_range * ((uv_max - uv_mid).array());
+    }
+    Vector2d rect_len;
+    rect_len << ob(1, 0) - ob(0, 0), ob(1, 1) - ob(0, 1);
+    int frame_points = 5;
+
+    s.rect_frame_V.resize(4 * frame_points, 2);
+    for (int i = 0; i < frame_points; i++)
+    {
+      // 0,0;0,1
+      s.rect_frame_V.row(i) << ob(0, 0), ob(0, 1) + i * rect_len(1) / frame_points;
+      // 0,0;1,1
+      s.rect_frame_V.row(i + frame_points)
+          << ob(0, 0) + i * rect_len(0) / frame_points,
+          ob(1, 1);
+      // 1,0;1,1
+      s.rect_frame_V.row(i + 2 * frame_points) << ob(1, 0), ob(1, 1) - i * rect_len(1) / frame_points;
+      // 1,0;0,1
+      s.rect_frame_V.row(i + 3 * frame_points)
+          << ob(1, 0) - i * rect_len(0) / frame_points,
+          ob(0, 1);
+      // 0,0;0,1
+    }
+    s.frame_ids = Eigen::VectorXi::LinSpaced(s.rect_frame_V.rows(), s.mv_num, s.mv_num + s.rect_frame_V.rows());
+  }
+
+  // Concatenate Vert and Edge
+  MatrixXd V;
+  MatrixXi E;
+  igl::cat(1, V_bnd, s.rect_frame_V, V);
+  E.resize(V.rows(), 2);
+  for (int i = 0; i < E.rows(); i++)
+    E.row(i) << i, i + 1;
+  int acc_bs = 0;
+  for (auto bs : s.bnd_sizes)
+  {
+    E(acc_bs + bs - 1, 1) = acc_bs;
+    acc_bs += bs;
+  }
+  E(V.rows() - 1, 1) = acc_bs;
+  assert(acc_bs == s.internal_bnd.size());
+
+  MatrixXd H = MatrixXd::Zero(s.component_sizes.size(), 2);
+  {
+    int hole_f = 0;
+    int hole_i = 0;
+    for (auto cs : s.component_sizes)
+    {
+      for (int i = 0; i < 3; i++)
+        H.row(hole_i) += m_uv.row(s.m_T(hole_f, i)); // redoing step 2
+      hole_f += cs;
+      hole_i++;
+    }
+  }
+  H /= 3.;
+
+  MatrixXd uv2;
+  igl::triangle::triangulate(V, E, H, std::basic_string<char>("qYYQ"), uv2, s.s_T);
+  auto bnd_n = s.internal_bnd.size();
+
+  for (auto i = 0; i < s.s_T.rows(); i++)
+    for (auto j = 0; j < s.s_T.cols(); j++)
+    {
+      auto &x = s.s_T(i, j);
+      if (x < bnd_n)
+        x = s.internal_bnd(x);
+      else
+        x += m_uv.rows() - bnd_n;
+    }
+
+  igl::cat(1, s.m_T, s.s_T, s.w_T);
+  s.w_uv.conservativeResize(m_uv.rows() - bnd_n + uv2.rows(), 2);
+  s.w_uv.bottomRows(uv2.rows() - bnd_n) = uv2.bottomRows(-bnd_n + uv2.rows());
+
+  update_scaffold(s);
+
+  // after_mesh_improve
+  compute_scaffold_gradient_matrix(s, s.Dx_s, s.Dy_s);
+
+  s.Dx_s.makeCompressed();
+  s.Dy_s.makeCompressed();
+  s.Dz_s.makeCompressed();
+  s.Ri_s = MatrixXd::Zero(s.Dx_s.rows(), s.dim * s.dim);
+  s.Ji_s.resize(s.Dx_s.rows(), s.dim * s.dim);
+  s.W_s.resize(s.Dx_s.rows(), s.dim * s.dim);
+}
+
+void add_new_patch(igl::SCAFData &s, const Eigen::MatrixXd &V_ref,
+                   const Eigen::MatrixXi &F_ref,
+                   const Eigen::RowVectorXd &center,
+                   const Eigen::MatrixXd &uv_init)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  assert(uv_init.rows() != 0);
+  Eigen::VectorXd M;
+  igl::doublearea(V_ref, F_ref, M);
+  s.mesh_measure += M.sum() / 2;
+
+  Eigen::VectorXi bnd;
+  Eigen::MatrixXd bnd_uv;
+
+  std::vector<std::vector<int>> all_bnds;
+  igl::boundary_loop(F_ref, all_bnds);
+  int num_holes = all_bnds.size() - 1;
+
+  s.component_sizes.push_back(F_ref.rows());
+
+  MatrixXd m_uv = s.w_uv.topRows(s.mv_num);
+  igl::cat(1, m_uv, uv_init, s.w_uv);
+
+  s.m_M.conservativeResize(s.mf_num + M.size());
+  s.m_M.bottomRows(M.size()) = M / 2;
+
+  for (auto cur_bnd : all_bnds)
+  {
+    s.internal_bnd.conservativeResize(s.internal_bnd.size() + cur_bnd.size());
+    s.internal_bnd.bottomRows(cur_bnd.size()) = Map<ArrayXi>(cur_bnd.data(), cur_bnd.size()) + s.mv_num;
+    s.bnd_sizes.push_back(cur_bnd.size());
+  }
+
+  s.m_T.conservativeResize(s.mf_num + F_ref.rows(), 3);
+  s.m_T.bottomRows(F_ref.rows()) = F_ref.array() + s.mv_num;
+  s.mf_num += F_ref.rows();
+
+  s.m_V.conservativeResize(s.mv_num + V_ref.rows(), 3);
+  s.m_V.bottomRows(V_ref.rows()) = V_ref;
+  s.mv_num += V_ref.rows();
+
+  s.rect_frame_V = MatrixXd();
+
+  mesh_improve(s);
+}
+
+void compute_jacobians(SCAFData &s, const Eigen::MatrixXd &V_new, bool whole)
+{
+  auto comp_J2 = [](const Eigen::MatrixXd &uv,
+                    const Eigen::SparseMatrix<double> &Dx,
+                    const Eigen::SparseMatrix<double> &Dy,
+                    Eigen::MatrixXd &Ji) {
+    // Ji=[D1*u,D2*u,D1*v,D2*v];
+    Ji.resize(Dx.rows(), 4);
+    Ji.col(0) = Dx * uv.col(0);
+    Ji.col(1) = Dy * uv.col(0);
+    Ji.col(2) = Dx * uv.col(1);
+    Ji.col(3) = Dy * uv.col(1);
+  };
+
+  Eigen::MatrixXd m_V_new = V_new.topRows(s.mv_num);
+  comp_J2(m_V_new, s.Dx_m, s.Dy_m, s.Ji_m);
+  if (whole)
+    comp_J2(V_new, s.Dx_s, s.Dy_s, s.Ji_s);
+}
+
+double compute_energy_from_jacobians(const Eigen::MatrixXd &Ji,
+                                     const Eigen::VectorXd &areas,
+                                     igl::MappingEnergyType energy_type)
+{
+  double energy = 0;
+  if (energy_type == igl::MappingEnergyType::SYMMETRIC_DIRICHLET)
+    energy = -4; // comply with paper description
+  return energy + igl::mapping_energy_with_jacobians(Ji, areas, energy_type, 0);
+}
+
+double compute_soft_constraint_energy(const SCAFData &s)
+{
+  double e = 0;
+  for (auto const &x : s.soft_cons)
+    e += s.soft_const_p * (x.second - s.w_uv.row(x.first)).squaredNorm();
+
+  return e;
+}
+
+double compute_energy(SCAFData &s, Eigen::MatrixXd &w_uv, bool whole)
+{
+  if (w_uv.rows() != s.v_num)
+    assert(!whole);
+  compute_jacobians(s, w_uv, whole);
+  double energy = compute_energy_from_jacobians(s.Ji_m, s.m_M, s.slim_energy);
+
+  if (whole)
+    energy += compute_energy_from_jacobians(s.Ji_s, s.s_M, s.scaf_energy);
+  energy += compute_soft_constraint_energy(s);
+  return energy;
+}
+
+void buildAm(const Eigen::VectorXd &sqrt_M,
+             const Eigen::SparseMatrix<double> &Dx,
+             const Eigen::SparseMatrix<double> &Dy,
+             const Eigen::MatrixXd &W,
+             Eigen::SparseMatrix<double> &Am)
+{
+  std::vector<Eigen::Triplet<double>> IJV;
+  Eigen::SparseMatrix<double> Dz;
+
+  Eigen::SparseMatrix<double> MDx = sqrt_M.asDiagonal() * Dx;
+  Eigen::SparseMatrix<double> MDy = sqrt_M.asDiagonal() * Dy;
+  igl::slim_buildA(MDx, MDy, Dz, W, IJV);
+
+  Am.setFromTriplets(IJV.begin(), IJV.end());
+  Am.makeCompressed();
+}
+
+void buildRhs(const Eigen::VectorXd &sqrt_M,
+              const Eigen::MatrixXd &W,
+              const Eigen::MatrixXd &Ri,
+              Eigen::VectorXd &f_rhs)
+{
+  const int dim = (W.cols() == 4) ? 2 : 3;
+  const int f_n = W.rows();
+  f_rhs.resize(dim * dim * f_n);
+
+  for (int i = 0; i < f_n; i++)
+  {
+    auto sqrt_area = sqrt_M(i);
+    f_rhs(i + 0 * f_n) = sqrt_area * (W(i, 0) * Ri(i, 0) + W(i, 1) * Ri(i, 1));
+    f_rhs(i + 1 * f_n) = sqrt_area * (W(i, 0) * Ri(i, 2) + W(i, 1) * Ri(i, 3));
+    f_rhs(i + 2 * f_n) = sqrt_area * (W(i, 2) * Ri(i, 0) + W(i, 3) * Ri(i, 1));
+    f_rhs(i + 3 * f_n) = sqrt_area * (W(i, 2) * Ri(i, 2) + W(i, 3) * Ri(i, 3));
+  }
+}
+
+void get_complement(const Eigen::VectorXi &bnd_ids, int v_n, Eigen::ArrayXi &unknown_ids)
+{ // get the complement of bnd_ids.
+  int assign = 0, i = 0;
+  for (int get = 0; i < v_n && get < bnd_ids.size(); i++)
+  {
+    if (bnd_ids(get) == i)
+      get++;
+    else
+      unknown_ids(assign++) = i;
+  }
+  while (i < v_n)
+    unknown_ids(assign++) = i++;
+  assert(assign + bnd_ids.size() == v_n);
+}
+
+void build_surface_linear_system(const SCAFData &s, Eigen::SparseMatrix<double> &L, Eigen::VectorXd &rhs)
+{
+  using namespace Eigen;
+  using namespace std;
+
+  const int v_n = s.v_num - (s.frame_ids.size());
+  const int dim = s.dim;
+  const int f_n = s.mf_num;
+
+  // to get the  complete A
+  Eigen::VectorXd sqrtM = s.m_M.array().sqrt();
+  Eigen::SparseMatrix<double> A(dim * dim * f_n, dim * v_n);
+  auto decoy_Dx_m = s.Dx_m;
+  decoy_Dx_m.conservativeResize(s.W_m.rows(), v_n);
+  auto decoy_Dy_m = s.Dy_m;
+  decoy_Dy_m.conservativeResize(s.W_m.rows(), v_n);
+  buildAm(sqrtM, decoy_Dx_m, decoy_Dy_m, s.W_m, A);
+
+  const VectorXi &bnd_ids = s.fixed_ids;
+  auto bnd_n = bnd_ids.size();
+  if (bnd_n == 0)
+  {
+
+    Eigen::SparseMatrix<double> At = A.transpose();
+    At.makeCompressed();
+
+    Eigen::SparseMatrix<double> id_m(At.rows(), At.rows());
+    id_m.setIdentity();
+
+    L = At * A;
+
+    Eigen::VectorXd frhs;
+    buildRhs(sqrtM, s.W_m, s.Ri_m, frhs);
+    rhs = At * frhs;
+  }
+  else
+  {
+    MatrixXd bnd_pos;
+    igl::slice(s.w_uv, bnd_ids, 1, bnd_pos);
+    ArrayXi known_ids(bnd_ids.size() * dim);
+    ArrayXi unknown_ids((v_n - bnd_ids.rows()) * dim);
+    get_complement(bnd_ids, v_n, unknown_ids);
+    VectorXd known_pos(bnd_ids.size() * dim);
+    for (int d = 0; d < dim; d++)
+    {
+      auto n_b = bnd_ids.rows();
+      known_ids.segment(d * n_b, n_b) = bnd_ids.array() + d * v_n;
+      known_pos.segment(d * n_b, n_b) = bnd_pos.col(d);
+      unknown_ids.block(d * (v_n - n_b), 0, v_n - n_b, unknown_ids.cols()) =
+          unknown_ids.topRows(v_n - n_b) + d * v_n;
+    }
+
+    Eigen::SparseMatrix<double> Au, Ae;
+    igl::slice(A, unknown_ids, 2, Au);
+    igl::slice(A, known_ids, 2, Ae);
+
+    Eigen::SparseMatrix<double> Aut = Au.transpose();
+    Aut.makeCompressed();
+
+    L = Aut * Au;
+
+    Eigen::VectorXd frhs;
+    buildRhs(sqrtM, s.W_m, s.Ri_m, frhs);
+
+    rhs = Aut * (frhs - Ae * known_pos);
+  }
+
+  // add soft constraints.
+  for (auto const &x : s.soft_cons)
+  {
+    int v_idx = x.first;
+
+    for (int d = 0; d < dim; d++)
+    {
+      rhs(d * (v_n) + v_idx) += s.soft_const_p * x.second(d); // rhs
+      L.coeffRef(d * v_n + v_idx,
+                 d * v_n + v_idx) += s.soft_const_p; // diagonal
+    }
+  }
+}
+
+void build_scaffold_linear_system(const SCAFData &s, Eigen::SparseMatrix<double> &L, Eigen::VectorXd &rhs)
+{
+  using namespace Eigen;
+
+  const int f_n = s.W_s.rows();
+  const int v_n = s.Dx_s.cols();
+  const int dim = s.dim;
+
+  Eigen::VectorXd sqrtM = s.s_M.array().sqrt();
+  Eigen::SparseMatrix<double> A(dim * dim * f_n, dim * v_n);
+  buildAm(sqrtM, s.Dx_s, s.Dy_s, s.W_s, A);
+
+  VectorXi bnd_ids;
+  igl::cat(1, s.fixed_ids, s.frame_ids, bnd_ids);
+
+  auto bnd_n = bnd_ids.size();
+  assert(bnd_n > 0);
+  MatrixXd bnd_pos;
+  igl::slice(s.w_uv, bnd_ids, 1, bnd_pos);
+
+  ArrayXi known_ids(bnd_ids.size() * dim);
+  ArrayXi unknown_ids((v_n - bnd_ids.rows()) * dim);
+
+  get_complement(bnd_ids, v_n, unknown_ids);
+
+  VectorXd known_pos(bnd_ids.size() * dim);
+  for (int d = 0; d < dim; d++)
+  {
+    auto n_b = bnd_ids.rows();
+    known_ids.segment(d * n_b, n_b) = bnd_ids.array() + d * v_n;
+    known_pos.segment(d * n_b, n_b) = bnd_pos.col(d);
+    unknown_ids.block(d * (v_n - n_b), 0, v_n - n_b, unknown_ids.cols()) =
+        unknown_ids.topRows(v_n - n_b) + d * v_n;
+  }
+  Eigen::VectorXd sqrt_M = s.s_M.array().sqrt();
+
+  // manual slicing for A(:, unknown/known)'
+  Eigen::SparseMatrix<double> Au, Ae;
+  igl::slice(A, unknown_ids, 2, Au);
+  igl::slice(A, known_ids, 2, Ae);
+
+  Eigen::SparseMatrix<double> Aut = Au.transpose();
+  Aut.makeCompressed();
+
+  L = Aut * Au;
+
+  Eigen::VectorXd frhs;
+  buildRhs(sqrtM, s.W_s, s.Ri_s, frhs);
+
+  rhs = Aut * (frhs - Ae * known_pos);
+}
+
+void solve_weighted_arap(SCAFData &s, Eigen::MatrixXd &uv)
+{
+  using namespace Eigen;
+  using namespace std;
+  int dim = s.dim;
+  igl::Timer timer;
+  timer.start();
+
+  VectorXi bnd_ids;
+  igl::cat(1, s.fixed_ids, s.frame_ids, bnd_ids);
+  const auto v_n = s.v_num;
+  const auto bnd_n = bnd_ids.size();
+  assert(bnd_n > 0);
+  MatrixXd bnd_pos;
+  igl::slice(s.w_uv, bnd_ids, 1, bnd_pos);
+
+  ArrayXi known_ids(bnd_n * dim);
+  ArrayXi unknown_ids((v_n - bnd_n) * dim);
+
+  get_complement(bnd_ids, v_n, unknown_ids);
+
+  VectorXd known_pos(bnd_ids.size() * dim);
+  for (int d = 0; d < dim; d++)
+  {
+    auto n_b = bnd_ids.rows();
+    known_ids.segment(d * n_b, n_b) = bnd_ids.array() + d * v_n;
+    known_pos.segment(d * n_b, n_b) = bnd_pos.col(d);
+    unknown_ids.block(d * (v_n - n_b), 0, v_n - n_b, unknown_ids.cols()) =
+        unknown_ids.topRows(v_n - n_b) + d * v_n;
+  }
+
+  Eigen::SparseMatrix<double> L;
+  Eigen::VectorXd rhs;
+
+  // fixed frame solving:
+  // x_e as the fixed frame, x_u for unknowns (mesh + unknown scaffold)
+  // min ||(A_u*x_u + A_e*x_e) - b||^2
+  // => A_u'*A_u*x_u  = Au'* (b - A_e*x_e) := Au'* b_u
+  //
+  // separate matrix build:
+  // min ||A_m x_m - b_m||^2 + ||A_s x_all - b_s||^2 + soft + proximal
+  // First change dimension of A_m to fit for x_all
+  // (Not just at the end, since x_all is flattened along dimensions)
+  // L = A_m'*A_m + A_s'*A_s + soft + proximal
+  // rhs = A_m'* b_m + A_s' * b_s + soft + proximal
+  //
+  Eigen::SparseMatrix<double> L_m, L_s;
+  Eigen::VectorXd rhs_m, rhs_s;
+  build_surface_linear_system(s, L_m, rhs_m);  // complete Am, with soft
+  build_scaffold_linear_system(s, L_s, rhs_s); // complete As, without proximal
+
+  L = L_m + L_s;
+  rhs = rhs_m + rhs_s;
+  L.makeCompressed();
+
+  Eigen::VectorXd unknown_Uc((v_n - s.frame_ids.size() - s.fixed_ids.size()) * dim), Uc(dim * v_n);
+
+  SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
+  unknown_Uc = solver.compute(L).solve(rhs);
+  igl::slice_into(unknown_Uc, unknown_ids.matrix(), 1, Uc);
+  igl::slice_into(known_pos, known_ids.matrix(), 1, Uc);
+
+  uv = Map<Matrix<double, -1, -1, Eigen::ColMajor>>(Uc.data(), v_n, dim);
+}
+
+double perform_iteration(SCAFData &s)
+{
+  Eigen::MatrixXd V_out = s.w_uv;
+  compute_jacobians(s, V_out, true);
+  igl::slim_update_weights_and_closest_rotations_with_jacobians(s.Ji_m, s.slim_energy, 0, s.W_m, s.Ri_m);
+  igl::slim_update_weights_and_closest_rotations_with_jacobians(s.Ji_s, s.scaf_energy, 0, s.W_s, s.Ri_s);
+  solve_weighted_arap(s, V_out);
+  auto whole_E = [&s](Eigen::MatrixXd &uv) { return compute_energy(s, uv, true); };
+
+  Eigen::MatrixXi w_T;
+  if (s.m_T.cols() == s.s_T.cols())
+    igl::cat(1, s.m_T, s.s_T, w_T);
+  else
+    w_T = s.s_T;
+  return igl::flip_avoiding_line_search(w_T, s.w_uv, V_out,
+                                        whole_E, -1) /
+         s.mesh_measure;
+}
+}
+}
+
+IGL_INLINE void igl::scaf_precompute(
+    const Eigen::MatrixXd &V,
+    const Eigen::MatrixXi &F,
+    const Eigen::MatrixXd &V_init,
+    igl::SCAFData &data,
+    igl::MappingEnergyType slim_energy,
+    Eigen::VectorXi &b,
+    Eigen::MatrixXd &bc,
+    double soft_p)
+{
+  Eigen::MatrixXd CN;
+  Eigen::MatrixXi FN;
+  igl::scaf::add_new_patch(data, V, F, Eigen::RowVector2d(0, 0), V_init);
+  data.soft_const_p = soft_p;
+  for (int i = 0; i < b.rows(); i++)
+    data.soft_cons[b(i)] = bc.row(i);
+  data.slim_energy = slim_energy;
+
+  auto &s = data;
+
+  if (!data.has_pre_calc)
+  {
+    int v_n = s.mv_num + s.sv_num;
+    int f_n = s.mf_num + s.sf_num;
+    int dim = s.dim;
+    Eigen::MatrixXd F1, F2, F3;
+    igl::local_basis(s.m_V, s.m_T, F1, F2, F3);
+    auto face_proj = [](Eigen::MatrixXd& F){
+      std::vector<Eigen::Triplet<double> >IJV;
+      int f_num = F.rows();
+      for(int i=0; i<F.rows(); i++) {
+        IJV.push_back(Eigen::Triplet<double>(i, i, F(i,0)));
+        IJV.push_back(Eigen::Triplet<double>(i, i+f_num, F(i,1)));
+        IJV.push_back(Eigen::Triplet<double>(i, i+2*f_num, F(i,2)));
+      }
+      Eigen::SparseMatrix<double> P(f_num, 3*f_num);
+      P.setFromTriplets(IJV.begin(), IJV.end());
+      return P;
+    };
+    Eigen::SparseMatrix<double> G;
+    igl::grad(s.m_V, s.m_T, G);
+    s.Dx_m = face_proj(F1) * G;
+    s.Dy_m = face_proj(F2) * G;
+
+    igl::scaf::compute_scaffold_gradient_matrix(s, s.Dx_s, s.Dy_s);
+
+    s.Dx_m.makeCompressed();
+    s.Dy_m.makeCompressed();
+    s.Ri_m = Eigen::MatrixXd::Zero(s.Dx_m.rows(), dim * dim);
+    s.Ji_m.resize(s.Dx_m.rows(), dim * dim);
+    s.W_m.resize(s.Dx_m.rows(), dim * dim);
+
+    s.Dx_s.makeCompressed();
+    s.Dy_s.makeCompressed();
+    s.Ri_s = Eigen::MatrixXd::Zero(s.Dx_s.rows(), dim * dim);
+    s.Ji_s.resize(s.Dx_s.rows(), dim * dim);
+    s.W_s.resize(s.Dx_s.rows(), dim * dim);
+
+    data.has_pre_calc = true;
+  }
+}
+
+IGL_INLINE Eigen::MatrixXd igl::scaf_solve(SCAFData &s, int iter_num)
+{
+  using namespace std;
+  using namespace Eigen;
+  s.energy = igl::scaf::compute_energy(s, s.w_uv, false) / s.mesh_measure;
+
+  for (int it = 0; it < iter_num; it++)
+  {
+    s.total_energy = igl::scaf::compute_energy(s, s.w_uv, true) / s.mesh_measure;
+    s.rect_frame_V = Eigen::MatrixXd();
+    igl::scaf::mesh_improve(s);
+
+    double new_weight = s.mesh_measure * s.energy / (s.sf_num * 100);
+    s.scaffold_factor = new_weight;
+    igl::scaf::update_scaffold(s);
+
+    s.total_energy = igl::scaf::perform_iteration(s);
+
+    s.energy =
+        igl::scaf::compute_energy(s, s.w_uv, false) / s.mesh_measure;
+  }
+
+  return s.w_uv.topRows(s.mv_num);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+#endif

include/igl/scaf.h

@@ -0,0 +1,101 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Zhongshi Jiang <jiangzs@nyu.edu>
+//
+// 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_SCAF_H
+#define IGL_SCAF_H
+
+#include "slim.h"
+#include "igl_inline.h"
+#include "MappingEnergyType.h"
+
+namespace igl
+{
+    // Use a similar interface to igl::slim
+    // Implement ready-to-use 2D version of the algorithm described in 
+    // SCAF: Simplicial Complex Augmentation Framework for Bijective Maps
+    // Zhongshi Jiang, Scott Schaefer, Daniele Panozzo, ACM Trancaction on Graphics (Proc. SIGGRAPH Asia 2017)
+    // For a complete implementation and customized UI, please refer to https://github.com/jiangzhongshi/scaffold-map
+
+    struct SCAFData
+    {
+      double scaffold_factor = 10;
+      igl::MappingEnergyType scaf_energy = igl::MappingEnergyType::SYMMETRIC_DIRICHLET;
+      igl::MappingEnergyType slim_energy = igl::MappingEnergyType::SYMMETRIC_DIRICHLET;
+
+      // Output
+      int dim = 2;
+      double total_energy; // scaffold + isometric 
+      double energy; // objective value
+
+      long mv_num = 0, mf_num = 0;
+      long sv_num = 0, sf_num = 0;
+      long v_num{}, f_num = 0;
+      Eigen::MatrixXd m_V; // input initial mesh V
+      Eigen::MatrixXi m_T; // input initial mesh F/T
+      // INTERNAL
+      Eigen::MatrixXd w_uv; // whole domain uv: mesh + free vertices
+      Eigen::MatrixXi s_T; // scaffold domain tets: scaffold tets
+      Eigen::MatrixXi w_T;
+
+      Eigen::VectorXd m_M; // mesh area or volume
+      Eigen::VectorXd s_M; // scaffold area or volume
+      Eigen::VectorXd w_M; // area/volume weights for whole
+      double mesh_measure; // area or volume
+      double proximal_p = 0;
+
+      Eigen::VectorXi frame_ids;
+      Eigen::VectorXi fixed_ids;
+
+      std::map<int, Eigen::RowVectorXd> soft_cons;
+      double soft_const_p = 1e4;
+
+      Eigen::VectorXi internal_bnd;
+      Eigen::MatrixXd rect_frame_V;
+      // multi-chart support
+      std::vector<int> component_sizes;
+      std::vector<int> bnd_sizes;
+    
+        // reweightedARAP interior variables.
+        bool has_pre_calc = false;
+        Eigen::SparseMatrix<double> Dx_s, Dy_s, Dz_s;
+        Eigen::SparseMatrix<double> Dx_m, Dy_m, Dz_m;
+        Eigen::MatrixXd Ri_m, Ji_m, Ri_s, Ji_s;
+        Eigen::MatrixXd W_m, W_s;
+    };
+
+
+// Compute necessary information to start using SCAF
+// Inputs:
+//		V           #V by 3 list of mesh vertex positions
+//		F           #F by 3/3 list of mesh faces (triangles/tets)
+//    data          igl::SCAFData
+//    slim_energy Energy type to minimize
+//    b           list of boundary indices into V (soft constraint)
+//    bc          #b by dim list of boundary conditions (soft constraint)
+//    soft_p      Soft penalty factor (can be zero)
+    IGL_INLINE void scaf_precompute(
+        const Eigen::MatrixXd &V,
+        const Eigen::MatrixXi &F,
+        const Eigen::MatrixXd &V_init,
+        SCAFData &data,
+        MappingEnergyType slim_energy,
+        Eigen::VectorXi& b,
+        Eigen::MatrixXd& bc,
+        double soft_p);
+
+
+// Run iter_num iterations of SCAF, with precomputed data
+// Outputs:
+//    V_o (in SLIMData): #V by dim list of mesh vertex positions
+  IGL_INLINE Eigen::MatrixXd scaf_solve(SCAFData &data, int iter_num);
+  }
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "scaf.cpp"
+#endif
+
+#endif //IGL_SCAF_H

include/igl/slice.cpp

@@ -355,6 +355,7 @@ template Eigen::Matrix<double, -1, -1, 0, -1, -1> igl::slice<Eigen::Matrix<doubl
 template void igl::slice<Eigen::PlainObjectBase<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::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, 3, 0, -1, 3>&);
 template void igl::slice<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
 template void igl::slice<unsigned int, unsigned int>(Eigen::SparseMatrix<unsigned int, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<unsigned int, 0, int>&);
+template void igl::slice<Eigen::SparseMatrix<double, 0, int>, Eigen::Array<int, -1, 1, 0, -1, 1>, Eigen::SparseMatrix<double, 0, int> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::DenseBase<Eigen::Array<int, -1, 1, 0, -1, 1> > const&, int, Eigen::SparseMatrix<double, 0, int>&);
 #ifdef WIN32
 template void igl::slice<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>>>(class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> const &, class Eigen::DenseBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> &);
 template void igl::slice<class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>>>(class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>> const &, class Eigen::DenseBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>> &);

include/igl/slim.h

@@ -9,6 +9,7 @@
 #define SLIM_H
 
 #include "igl_inline.h"
+#include "MappingEnergyType.h"
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 
@@ -29,16 +30,7 @@ struct SLIMData
   // Input
   Eigen::MatrixXd V; // #V by 3 list of mesh vertex positions
   Eigen::MatrixXi F; // #F by 3/3 list of mesh faces (triangles/tets)
-  enum SLIM_ENERGY
-  {
-    ARAP,
-    LOG_ARAP,
-    SYMMETRIC_DIRICHLET,
-    CONFORMAL,
-    EXP_CONFORMAL,
-    EXP_SYMMETRIC_DIRICHLET
-  };
-  SLIM_ENERGY slim_energy;
+  MappingEnergyType slim_energy;
 
   // Optional Input
   // soft constraints
@@ -64,9 +56,7 @@ struct SLIMData
   Eigen::VectorXd WGL_M;
   Eigen::VectorXd rhs;
   Eigen::MatrixXd Ri,Ji;
-  Eigen::VectorXd W_11; Eigen::VectorXd W_12; Eigen::VectorXd W_13;
-  Eigen::VectorXd W_21; Eigen::VectorXd W_22; Eigen::VectorXd W_23;
-  Eigen::VectorXd W_31; Eigen::VectorXd W_32; Eigen::VectorXd W_33;
+  Eigen::MatrixXd W;
   Eigen::SparseMatrix<double> Dx,Dy,Dz;
   int f_n,v_n;
   bool first_solve;
@@ -94,7 +84,7 @@ IGL_INLINE void slim_precompute(
   const Eigen::MatrixXi& F,
   const Eigen::MatrixXd& V_init,
   SLIMData& data,
-  SLIMData::SLIM_ENERGY slim_energy,
+  MappingEnergyType slim_energy,
   Eigen::VectorXi& b,
   Eigen::MatrixXd& bc,
   double soft_p);
@@ -104,6 +94,18 @@ IGL_INLINE void slim_precompute(
 //    V_o (in SLIMData): #V by dim list of mesh vertex positions
 IGL_INLINE Eigen::MatrixXd slim_solve(SLIMData& data, int iter_num);
 
+// Internal Routine. Exposed for Integration with SCAF
+IGL_INLINE void slim_update_weights_and_closest_rotations_with_jacobians(const Eigen::MatrixXd &Ji,
+                                                                    igl::MappingEnergyType slim_energy,
+                                                                    double exp_factor,
+                                                                    Eigen::MatrixXd &W,
+                                                                    Eigen::MatrixXd &Ri);
+
+IGL_INLINE void slim_buildA(const Eigen::SparseMatrix<double> &Dx,
+                        const Eigen::SparseMatrix<double> &Dy,
+                        const Eigen::SparseMatrix<double> &Dz,
+                        const Eigen::MatrixXd &W,
+                        std::vector<Eigen::Triplet<double> > & IJV);
 } // END NAMESPACE
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/sparse_cached.h

@@ -39,7 +39,7 @@ namespace igl
   // Example:
   //   Eigen::SparseMatrix<double> A;
   //   std::vector<Eigen::Triplet<double> > IJV;
-  //   buildA(IJV);
+  //   slim_buildA(IJV);
   //   if (A.rows() == 0)
   //   {
   //     A = Eigen::SparseMatrix<double>(rows,cols);

include/igl/topological_hole_fill.cpp

@@ -0,0 +1,55 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Zhongshi Jiang <jiangzs@nyu.edu>
+//
+// 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 "topological_hole_fill.h"
+  template <
+  typename DerivedF,
+  typename Derivedb,
+  typename VectorIndex,
+  typename DerivedF_filled>
+IGL_INLINE void igl::topological_hole_fill(
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<Derivedb> & b,
+  const std::vector<VectorIndex> & holes,
+  Eigen::PlainObjectBase<DerivedF_filled> &F_filled)
+{
+  int n_filled_faces = 0;
+  int num_holes = holes.size();
+  int real_F_num = F.rows();
+  const int V_rows = F.maxCoeff()+1;
+
+  for (int i = 0; i < num_holes; i++)
+    n_filled_faces += holes[i].size();
+  F_filled.resize(n_filled_faces + real_F_num, 3);
+  F_filled.topRows(real_F_num) = F;
+
+  int new_vert_id = V_rows;
+  int new_face_id = real_F_num;
+
+  for (int i = 0; i < num_holes; i++, new_vert_id++)
+  {
+    int cur_bnd_size = holes[i].size();
+    int it = 0;
+    int back = holes[i].size() - 1;
+    F_filled.row(new_face_id++) << holes[i][it], holes[i][back], new_vert_id;
+    while (it != back)
+    {
+      F_filled.row(new_face_id++)
+          << holes[i][(it + 1)],
+          holes[i][(it)], new_vert_id;
+      it++;
+    }
+  }
+  assert(new_face_id == F_filled.rows());
+  assert(new_vert_id == V_rows + num_holes);
+
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::topological_hole_fill<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, std::vector<int, std::allocator<int> > >(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1,0, -1, 1> > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/topological_hole_fill.h

@@ -0,0 +1,44 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Zhongshi Jiang <jiangzs@nyu.edu>
+//
+// 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_TOPOLOGICAL_HOLE_FILL_H
+#define IGL_TOPOLOGICAL_HOLE_FILL_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+    
+  // Topological fill hole on a mesh, with one additional vertex each hole
+  // Index of new abstract vertices will be F.maxCoeff() + (index of hole)
+  //
+  // Inputs:
+  //   F  #F by simplex-size list of element indices
+  //   b  #b boundary indices to preserve
+  //   holes vector of hole loops to fill
+  // Outputs:
+  //   F_filled  input F stacked with filled triangles.
+  //
+  template <
+  typename DerivedF,
+  typename Derivedb,
+  typename VectorIndex,
+  typename DerivedF_filled>
+IGL_INLINE void topological_hole_fill(
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<Derivedb> & b,
+  const std::vector<VectorIndex> & holes,
+  Eigen::PlainObjectBase<DerivedF_filled> &F_filled);
+
+}
+
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "topological_hole_fill.cpp"
+#endif
+
+#endif

tutorial/710_SLIM/main.cpp → tutorial/709_SLIM/main.cpp

@@ -1,15 +1,16 @@
 #include <iostream>
 
-#include "igl/slim.h"
+#include <igl/slim.h>
 
-#include "igl/vertex_components.h"
-#include "igl/readOBJ.h"
-#include "igl/writeOBJ.h"
-#include "igl/Timer.h"
+#include <igl/vertex_components.h>
+#include <igl/readOBJ.h>
+#include <igl/writeOBJ.h>
+#include <igl/Timer.h>
 
-#include "igl/boundary_loop.h"
-#include "igl/map_vertices_to_circle.h"
-#include "igl/harmonic.h"
+#include <igl/boundary_loop.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/harmonic.h>
+#include <igl/MappingEnergyType.h>
 #include <igl/serialize.h>
 #include <igl/read_triangle_mesh.h>
 #include <igl/opengl/glfw/Viewer.h>
@@ -96,9 +97,9 @@ void param_2d_demo_iter(igl::opengl::glfw::Viewer& viewer) {
 
     cout << "initialized parametrization" << endl;
 
-    sData.slim_energy = igl::SLIMData::SYMMETRIC_DIRICHLET;
+    sData.slim_energy = igl::MappingEnergyType::SYMMETRIC_DIRICHLET;
     Eigen::VectorXi b; Eigen::MatrixXd bc;
-    slim_precompute(V,F,uv_init,sData, igl::SLIMData::SYMMETRIC_DIRICHLET, b,bc,0);
+    slim_precompute(V,F,uv_init,sData, igl::MappingEnergyType::SYMMETRIC_DIRICHLET, b,bc,0);
 
     uv_scale_param = 15 * (1./sqrt(sData.mesh_area));
     viewer.data().set_mesh(V, F);
@@ -109,10 +110,12 @@ void param_2d_demo_iter(igl::opengl::glfw::Viewer& viewer) {
 
     first_iter = false;
   } else {
+    timer.start();
     slim_solve(sData,1); // 1 iter
     viewer.data().set_uv(sData.V_o*uv_scale_param);
-    cout << "time = " << timer.getElapsedTime() << endl;
   }
+  cout << "time = " << timer.getElapsedTime() << endl;
+  cout << "energy = " << sData.energy << endl;
 }
 
 void soft_const_demo_iter(igl::opengl::glfw::Viewer& viewer) {
@@ -127,7 +130,7 @@ void soft_const_demo_iter(igl::opengl::glfw::Viewer& viewer) {
     Eigen::VectorXi b; Eigen::MatrixXd bc;
     get_soft_constraint_for_circle(V_0,F,b,bc);
     double soft_const_p = 1e5;
-    slim_precompute(V,F,V_0,sData,igl::SLIMData::SYMMETRIC_DIRICHLET,b,bc,soft_const_p);
+    slim_precompute(V,F,V_0,sData,igl::MappingEnergyType::SYMMETRIC_DIRICHLET,b,bc,soft_const_p);
 
     viewer.data().set_mesh(V, F);
     viewer.core.align_camera_center(V,F);
@@ -144,6 +147,7 @@ void soft_const_demo_iter(igl::opengl::glfw::Viewer& viewer) {
 
 void deform_3d_demo_iter(igl::opengl::glfw::Viewer& viewer) {
   if (first_iter) {
+    timer.start();
     igl::readOBJ(TUTORIAL_SHARED_PATH "/cube_40k.obj", V, F);
 
     Eigen::MatrixXd V_0 = V;
@@ -152,16 +156,19 @@ void deform_3d_demo_iter(igl::opengl::glfw::Viewer& viewer) {
 
     double soft_const_p = 1e5;
     sData.exp_factor = 5.0;
-    slim_precompute(V,F,V_0,sData,igl::SLIMData::EXP_CONFORMAL,b,bc,soft_const_p);
+    slim_precompute(V,F,V_0,sData,igl::MappingEnergyType::EXP_CONFORMAL,b,bc,soft_const_p);
     //cout << "precomputed" << endl;
 
     first_iter = false;
     display_3d_mesh(viewer);
 
   } else {
+    timer.start();
     slim_solve(sData,1); // 1 iter
     display_3d_mesh(viewer);
   }
+  cout << "time = " << timer.getElapsedTime() << endl;
+  cout << "energy = " << sData.energy << endl;
 }
 
 void display_3d_mesh(igl::opengl::glfw::Viewer& viewer) {

tutorial/710_SCAF/CMakeLists.txt

@@ -0,0 +1,5 @@
+get_filename_component(PROJECT_NAME ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+project(${PROJECT_NAME})
+
+add_executable(${PROJECT_NAME}_bin main.cpp)
+target_link_libraries(${PROJECT_NAME}_bin igl::core igl::opengl igl::opengl_glfw igl::triangle tutorials)

tutorial/710_SCAF/main.cpp

@@ -0,0 +1,127 @@
+#include <igl/scaf.h>
+#include <igl/arap.h>
+#include <igl/boundary_loop.h>
+#include <igl/harmonic.h>
+#include <igl/map_vertices_to_circle.h>
+#include <igl/readOBJ.h>
+#include <igl/Timer.h>
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/MappingEnergyType.h>
+#include <igl/doublearea.h>
+#include <igl/PI.h>
+#include <igl/flipped_triangles.h>
+#include <igl/topological_hole_fill.h>
+
+#include "tutorial_shared_path.h"
+
+Eigen::MatrixXd V;
+Eigen::MatrixXi F;
+Eigen::MatrixXd V_uv;
+igl::Timer timer;
+igl::SCAFData scaf_data;
+
+bool show_uv = false;
+float uv_scale = 0.2;
+
+bool key_down(igl::opengl::glfw::Viewer& viewer, unsigned char key, int modifier)
+{
+  if (key == '1')
+    show_uv = false;
+  else if (key == '2')
+    show_uv = true;
+
+  if (key == ' ')
+  {
+    timer.start();
+    igl::scaf_solve(scaf_data, 1);
+    std::cout << "time = " << timer.getElapsedTime() << std::endl;
+  }
+
+  const auto& V_uv = uv_scale * scaf_data.w_uv.topRows(V.rows());
+  if (show_uv)
+  {
+    viewer.data().clear();
+    viewer.data().set_mesh(V_uv,F);
+    viewer.data().set_uv(V_uv);
+    viewer.core.align_camera_center(V_uv,F);
+  }
+  else
+  {
+    viewer.data().set_mesh(V,F);
+    viewer.data().set_uv(V_uv);
+    viewer.core.align_camera_center(V,F);
+  }
+
+  viewer.data().compute_normals();
+
+  return false;
+}
+
+int main(int argc, char *argv[])
+{
+  using namespace std;
+  // Load a mesh in OFF format
+  igl::readOBJ(TUTORIAL_SHARED_PATH "/camel_b.obj", V, F);
+
+  Eigen::MatrixXd bnd_uv, uv_init;
+
+  Eigen::VectorXd M;
+  igl::doublearea(V, F, M);
+  std::vector<std::vector<int>> all_bnds;
+  igl::boundary_loop(F, all_bnds);
+
+  // Heuristic primary boundary choice: longest
+  auto primary_bnd = std::max_element(all_bnds.begin(), all_bnds.end(), [](const std::vector<int> &a, const std::vector<int> &b) { return a.size()<b.size(); });
+
+  Eigen::VectorXi bnd = Eigen::Map<Eigen::VectorXi>(primary_bnd->data(), primary_bnd->size());
+
+  igl::map_vertices_to_circle(V, bnd, bnd_uv);
+  bnd_uv *= sqrt(M.sum() / (2 * igl::PI));
+  if (all_bnds.size() == 1)
+  {
+    if (bnd.rows() == V.rows()) // case: all vertex on boundary
+    {
+      uv_init.resize(V.rows(), 2);
+      for (int i = 0; i < bnd.rows(); i++)
+        uv_init.row(bnd(i)) = bnd_uv.row(i);
+    }
+    else
+    {
+      igl::harmonic(V, F, bnd, bnd_uv, 1, uv_init);
+      if (igl::flipped_triangles(uv_init, F).size() != 0)
+        igl::harmonic(F, bnd, bnd_uv, 1, uv_init); // fallback uniform laplacian
+    }
+  }
+  else
+  {
+    // if there is a hole, fill it and erase additional vertices.
+    all_bnds.erase(primary_bnd);
+    Eigen::MatrixXi F_filled;
+    igl::topological_hole_fill(F, bnd, all_bnds, F_filled);
+    igl::harmonic(F_filled, bnd, bnd_uv ,1, uv_init);
+    uv_init = uv_init.topRows(V.rows());
+  }
+
+  Eigen::VectorXi b; Eigen::MatrixXd bc;
+  igl::scaf_precompute(V, F, uv_init, scaf_data, igl::MappingEnergyType::SYMMETRIC_DIRICHLET, b, bc, 0);
+
+  // Plot the mesh
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(V, F);
+  const auto& V_uv = uv_scale * scaf_data.w_uv.topRows(V.rows());
+  viewer.data().set_uv(V_uv);
+  viewer.callback_key_down = &key_down;
+
+  // Enable wireframe
+  viewer.data().show_lines = true;
+
+  // Draw checkerboard texture
+  viewer.data().show_texture = true;
+
+
+  std::cerr << "Press space for running an iteration." << std::endl;
+  std::cerr << "Press 1 for Mesh 2 for UV" << std::endl;
+
+  // Launch the viewer
+  viewer.launch();
+}

tutorial/CMakeLists.txt

@@ -144,7 +144,8 @@ if(TUTORIALS_CHAPTER7)
   endif()
   add_subdirectory("707_SweptVolume")
   add_subdirectory("708_Picking")
-  add_subdirectory("710_SLIM")
+  add_subdirectory("709_SLIM")
+  add_subdirectory("710_SCAF")
   add_subdirectory("711_Subdivision")
   add_subdirectory("712_DataSmoothing")
   add_subdirectory("713_ShapeUp")

tutorial/images/710_SCAF.png.REMOVED.git-id

@@ -0,0 +1 @@
+0d5f21bca3dff436750d4e1500fad9344c6f36c1

tutorial/shared/camel_b.obj.REMOVED.git-id

@@ -0,0 +1 @@
+ad82107330848e2926775c0a8234ad6e301e03ff