Merge pull request #105 from MarkGillespie/BFF

Add implementation of BFF

docs/docs/surface/algorithms/parameterization.md

@@ -0,0 +1,104 @@
+This section describes the algorithms in geometry-central for surface parameterization, which compute maps from surfaces meshes to the plane.
+
+Note that these procedures all depend on the _intrinsic_ geometry of a surface (via the `IntrinsicGeometryInterface`). Therefore, these routines can be run on abstract geometric domains as well as traditional surfaces in 3D.
+
+## Boundary First Flattening
+![parameterized face](/media/bff.png){: style="max-width: 15em; display: block; margin-left: auto; margin-right: auto;"}
+This algorithm is described in the paper [Boundary First Flattening](http://www.cs.cmu.edu/~kmcrane/Projects/BoundaryFirstFlattening/paper.pdf). It computes a _conformal_ parameterization of a surface mesh, allowing the user to specify target angles or scale factors along the boundary of the mesh. The input mesh must be a topological disk.
+
+`#include "geometrycentral/surface/boundary_first_flattening.h"`
+
+### Single Parameterizations
+
+A one-off utility function is provided to compute single parameterizations. Repeated parameterizations of the same mesh can be computed more efficiently using the utility class `BFF` below.
+
+Example:
+```cpp
+#include "geometrycentral/surface/boundary_first_flattening.h"
+#include "geometrycentral/surface/meshio.h"
+
+// Load a mesh
+std::unique_ptr<ManifoldSurfaceMesh> mesh;
+std::unique_ptr<VertexPositionGeometry> geometry;
+std::tie(mesh, geometry) = readManifoldSurfaceMesh(filename);
+
+VertexData<Vector2> parameterization = parameterizeBFF(*mesh, *geometry);
+```
+
+??? func "`#!cpp VertexData<Vector2> parameterizeBFF(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom)`"
+    Conformally parameterize the input mesh. Picks boundary conditions to minimize area distortion (i.e. sets the conformal scale factor to 0 along the boundary).
+
+??? func "`#!cpp VertexData<Vector2> parameterizeBFFfromScaleFactors(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, const VertexData<double>& boundaryScaleFactors)`"
+    Conformally parameterize the input mesh, setting the scale factors to the given values along the boundary.
+    Although `boundaryScaleFactors` is a `VertexData` object which stores values at all vertices, only the values at boundary vertices are used by the algorithm. All other values are ignored.
+    
+??? func "`#!cpp VertexData<Vector2> parameterizeBFFfromExteriorAngles(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom, const VertexData<double>& exteriorAngles)`"
+    Conformally parameterize the input mesh so that the boundary vertices of the parameterized mesh have the given exterior angles. The exterior angles must sum up to $2\pi$ along the boundary.
+    
+    Although `exteriorAngles` is a `VertexData` object which stores values at all vertices, only the values at boundary vertices are used by the algorithm. All other values are ignored.
+
+### Repeated Parameterization
+
+The stateful class `BFF` does precomputation when constructed to efficiently compute many parameterizations of the same mesh.
+
+Example:
+```cpp
+#include "geometrycentral/surface/boundary_first_flattening.h"
+#include "geometrycentral/surface/meshio.h"
+
+// Load a mesh
+std::unique_ptr<ManifoldSurfaceMesh> mesh;
+std::unique_ptr<VertexPositionGeometry> geometry;
+std::tie(mesh, geometry) = readManifoldSurfaceMesh(filename);
+
+VertexData<double> boundaryScaleFactors = /* target scale factors */;
+VertexData<double> exteriorAngles = /* target exteriorAngles */;
+
+BFF bff(*mesh, *geometry);
+VertexData<Vector2> parameterization1 = bff.flattenFromScaleFactors(boundaryScaleFactors);
+VertexData<Vector2> parameterization2 = bff.flattenFromExteriorAngles(exteriorAngles);
+```
+
+??? func "`#!cpp BFF::BFF(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geom)`"
+
+    Create a new solver for boundary first flattening. Most precomputation is done immediately when the object is constructed, although some additional precomputation may be done lazily later on.
+
+??? func "`#!cpp VertexData<Vector2> BFF::flatten()`"
+
+    Compute a conformal parameterization which minimizes area distortion (i.e. sets the scale factor to 0 along the boundary).
+
+??? func "`#!cpp VertexData<Vector2> BFF::flattenFromScaleFactors(const VertexData<double>& boundaryScaleFactors)`"
+
+    Compute a conformal parameterization with the given scale factor along the boundary.
+    Although `boundaryScaleFactors` is a `VertexData` object which stores values at all vertices, only the values at boundary vertices are used by the algorithm. All other values are ignored.
+
+??? func "`#!cpp VertexData<Vector2> BFF::flattenFromExteriorAngles(const VertexData<double>& exteriorAngles)`"
+
+    Compute a conformal parameterization with the given exterior angles along the boundary. The exterior angles must sum up to $2\pi$ along the boundary.
+
+    Although `exteriorAngles` is a `VertexData` object which stores values at all vertices, only the values at boundary vertices are used by the algorithm. All other values are ignored.
+
+### Citation
+
+If this algorithm contributes to academic work, please cite the following paper:
+
+```bib
+@article{Sawhney:2017:BFF,
+author = {Sawhney, Rohan and Crane, Keenan},
+title = {Boundary First Flattening},
+journal = {ACM Transactions on Graphics (TOG)},
+volume = {37},
+number = {1},
+month = dec,
+year = {2017},
+issn = {0730-0301},
+pages = {5:1--5:14},
+articleno = {5},
+numpages = {14},
+url = {http://doi.acm.org/10.1145/3132705},
+doi = {10.1145/3132705},
+acmid = {3132705},
+publisher = {ACM},
+address = {New York, NY, USA}
+}
+```

docs/mkdocs.yml

@@ -81,6 +81,7 @@ nav:
       - 'Mesh Graph Algorithms' : 'surface/algorithms/mesh_graph_algorithms.md'
       - 'Robust Geometry' : 'surface/algorithms/robust_geometry.md'
       - 'Flip Geodesics' : 'surface/algorithms/flip_geodesics.md'
+      - 'Parameterization' : 'surface/algorithms/parameterization.md'
     - Intrinsic Triangulations:
       - 'Basics' : 'surface/intrinsic_triangulations/basics.md'
       - 'Common Subdivision' : 'surface/intrinsic_triangulations/common_subdivision.md'

include/geometrycentral/surface/boundary_first_flattening.h

@@ -0,0 +1,52 @@
+#pragma once
+
+#include "geometrycentral/numerical/linear_algebra_utilities.h"
+#include "geometrycentral/numerical/linear_solvers.h"
+#include "geometrycentral/surface/intrinsic_geometry_interface.h"
+#include "geometrycentral/surface/manifold_surface_mesh.h"
+
+namespace geometrycentral {
+namespace surface {
+
+VertexData<Vector2> parameterizeBFF(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geo);
+VertexData<Vector2> parameterizeBFFfromScaleFactors(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geo,
+                                                    const VertexData<double>& boundaryScaleFactors);
+VertexData<Vector2> parameterizeBFFfromExteriorAngles(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geo,
+                                                      const VertexData<double>& exteriorAngles);
+
+class BFF {
+public:
+  BFF(ManifoldSurfaceMesh& mesh_, IntrinsicGeometryInterface& geo_);
+
+  VertexData<Vector2> flatten();
+  VertexData<Vector2> flattenFromScaleFactors(const VertexData<double>& uBdy);
+  VertexData<Vector2> flattenFromExteriorAngles(const VertexData<double>& kBdy);
+  VertexData<Vector2> flattenFromBoth(const Vector<double>& uBdy, const Vector<double>& kBdy);
+
+  Vector<double> dirichletToNeumann(const Vector<double>& uBdy);
+  Vector<double> neumannToDirichlet(const Vector<double>& kBdy);
+
+  std::array<Vector<double>, 2> computeBoundaryPositions(const Vector<double>& uBdy, const Vector<double>& kBdy);
+
+protected:
+  ManifoldSurfaceMesh& mesh;
+  IntrinsicGeometryInterface& geo;
+  VertexData<size_t> vIdx;
+  VertexData<int> iIdx, bIdx;
+
+
+  SparseMatrix<double> L, Lii, Lib, Lbb;
+  std::unique_ptr<PositiveDefiniteSolver<double>> Liisolver;
+  std::unique_ptr<PositiveDefiniteSolver<double>> Lsolver;
+  Vector<double> Omegai, Omegab;
+
+  Vector<bool> isInterior;
+  size_t nInterior, nBoundary, nVertices;
+
+  BlockDecompositionResult<double> Ldecomp;
+
+  void ensureHaveLSolver();
+};
+
+} // namespace surface
+} // namespace geometrycentral

src/CMakeLists.txt

@@ -17,6 +17,7 @@ SET(SRCS
   surface/rich_surface_mesh_data.cpp
 
   surface/base_geometry_interface.cpp
+  surface/boundary_first_flattening.cpp
   surface/intrinsic_geometry_interface.cpp
   surface/extrinsic_geometry_interface.cpp
   surface/exact_geodesic_helpers.cpp
@@ -89,6 +90,7 @@ SET(HEADERS
   ${INCLUDE_ROOT}/surface/barycentric_coordinate_helpers.h
   ${INCLUDE_ROOT}/surface/barycentric_coordinate_helpers.ipp
   ${INCLUDE_ROOT}/surface/base_geometry_interface.h
+  ${INCLUDE_ROOT}/surface/boundary_first_flattening.h
   ${INCLUDE_ROOT}/surface/detect_symmetry.h
   ${INCLUDE_ROOT}/surface/direction_fields.h
   ${INCLUDE_ROOT}/surface/edge_length_geometry.h

src/surface/boundary_first_flattening.cpp

@@ -0,0 +1,205 @@
+#include "geometrycentral/surface/boundary_first_flattening.h"
+
+namespace geometrycentral {
+namespace surface {
+
+VertexData<Vector2> parameterizeBFF(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geo) {
+  BFF bff(mesh, geo);
+  return bff.flatten();
+}
+
+VertexData<Vector2> parameterizeBFFfromScaleFactors(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geo,
+                                                    const VertexData<double>& boundaryScaleFactors) {
+  BFF bff(mesh, geo);
+  return bff.flattenFromScaleFactors(boundaryScaleFactors);
+}
+
+VertexData<Vector2> parameterizeBFFfromExteriorAngles(ManifoldSurfaceMesh& mesh, IntrinsicGeometryInterface& geo,
+                                                      const VertexData<double>& exteriorAngles) {
+  BFF bff(mesh, geo);
+  return bff.flattenFromExteriorAngles(exteriorAngles);
+}
+
+BFF::BFF(ManifoldSurfaceMesh& mesh_, IntrinsicGeometryInterface& geo_) : mesh(mesh_), geo(geo_) {
+
+  GC_SAFETY_ASSERT(mesh.eulerCharacteristic() == 2 && mesh.nBoundaryLoops() == 1,
+                   "Input to BFF must be a topological disk");
+
+  vIdx = mesh.getVertexIndices();
+  iIdx = VertexData<int>(mesh, -1);
+  bIdx = VertexData<int>(mesh, -1);
+
+  geo.requireCotanLaplacian();
+
+  L = geo.cotanLaplacian;
+
+  isInterior = Vector<bool>(mesh.nVertices());
+  size_t iB = 0, iI = 0;
+  for (Vertex v : mesh.vertices()) {
+    if (v.isBoundary()) {
+      bIdx[v] = iB++;
+      isInterior(vIdx[v]) = false;
+    } else {
+      iIdx[v] = iI++;
+      isInterior(vIdx[v]) = true;
+    }
+  }
+  nVertices = mesh.nVertices();
+  nInterior = mesh.nInteriorVertices();
+  nBoundary = nVertices - nInterior;
+
+  geo.requireCotanLaplacian();
+  SparseMatrix<double> L = geo.cotanLaplacian;
+  shiftDiagonal(L, 1e-12);
+
+  Ldecomp = blockDecomposeSquare(L, isInterior);
+
+  Lii = Ldecomp.AA;
+  Lib = Ldecomp.AB;
+  Lbb = Ldecomp.BB;
+
+  // TODO: extract this factorization from a full factorization of L
+  Liisolver.reset(new PositiveDefiniteSolver<double>(Lii));
+
+  geo.requireVertexAngleSums();
+  Omegai = Vector<double>(nInterior);
+  Omegab = Vector<double>(nBoundary);
+  for (Vertex v : mesh.vertices()) {
+    if (v.isBoundary()) {
+      Omegab(bIdx[v]) = M_PI - geo.vertexAngleSums[v];
+    } else {
+      Omegai(iIdx[v]) = 2 * M_PI - geo.vertexAngleSums[v];
+    }
+  }
+}
+
+VertexData<Vector2> BFF::flatten() {
+  VertexData<double> u(mesh, 0);
+  return flattenFromScaleFactors(u);
+}
+
+VertexData<Vector2> BFF::flattenFromScaleFactors(const VertexData<double>& uData) {
+  // Extract boundary values
+  Vector<double> uBdy, ignore;
+  decomposeVector(Ldecomp, uData.toVector(), ignore, uBdy);
+
+  // Compute complementary data
+  Vector<double> kBdy = dirichletToNeumann(uBdy);
+
+  // Flatten
+  return flattenFromBoth(uBdy, kBdy);
+}
+
+VertexData<Vector2> BFF::flattenFromExteriorAngles(const VertexData<double>& kData) {
+  // Extract boundary values
+  Vector<double> kBdy, ignore;
+  decomposeVector(Ldecomp, kData.toVector(), ignore, kBdy);
+
+  GC_SAFETY_ASSERT(abs(kBdy.sum() - 2 * M_PI) < 1e-3,
+                   "BFF error: target exterior angles must sum to 2 pi, but the input sums to " +
+                       std::to_string(kBdy.sum()));
+
+  // Compute complementary data
+  Vector<double> uBdy = neumannToDirichlet(kBdy);
+
+  // Flatten
+  return flattenFromBoth(uBdy, kBdy);
+}
+
+VertexData<Vector2> BFF::flattenFromBoth(const Vector<double>& uBdy, const Vector<double>& kBdy) {
+
+  Vector<double> boundaryX, boundaryY;
+  std::tie(boundaryX, boundaryY) = tuple_cat(computeBoundaryPositions(uBdy, kBdy));
+
+  Vector<double> interiorX = Liisolver->solve(-Lib * boundaryX);
+  Vector<double> interiorY = Liisolver->solve(-Lib * boundaryY);
+
+  VertexData<Vector2> parm(mesh);
+  for (Vertex v : mesh.vertices()) {
+    if (v.isBoundary()) {
+      size_t iV = bIdx[v];
+      parm[v] = Vector2{boundaryX(iV), boundaryY(iV)};
+    } else {
+      size_t iV = iIdx[v];
+      parm[v] = Vector2{interiorX(iV), interiorY(iV)};
+    }
+  }
+
+  return parm;
+}
+
+Vector<double> BFF::dirichletToNeumann(const Vector<double>& uBdy) {
+  return Omegab - (Lib.transpose() * Liisolver->solve(Omegai - Lib * uBdy)) - Lbb * uBdy;
+}
+
+Vector<double> BFF::neumannToDirichlet(const Vector<double>& kBdy) {
+  // Convert Neumann data to Dirichlet data by solving the Poisson equation and reading off values
+  ensureHaveLSolver();
+  Vector<double> rhs = reassembleVector(Ldecomp, Omegai, Vector<double>(Omegab - kBdy));
+  Vector<double> fullSolution = -Lsolver->solve(rhs);
+  Vector<double> uBdy, ignore;
+  decomposeVector(Ldecomp, fullSolution, ignore, uBdy);
+  double uMean = uBdy.mean(); // Ensure that u has mean 0
+  for (int i = 0; i < uBdy.size(); i++) uBdy(i) -= uMean;
+  return uBdy;
+}
+
+std::array<Vector<double>, 2> BFF::computeBoundaryPositions(const Vector<double>& uBdy, const Vector<double>& kBdy) {
+
+  auto src = [](Edge e) { return e.halfedge().vertex(); };
+  auto dst = [](Edge e) { return e.halfedge().next().vertex(); };
+
+  geo.requireEdgeLengths();
+
+  double phi = 0;
+
+  std::vector<Eigen::Triplet<double>> Ntriplets;
+  Vector<double> targetLength(nBoundary);
+  DenseMatrix<double> T(2, nBoundary);
+
+  for (BoundaryLoop b : mesh.boundaryLoops()) {
+    for (Edge e : b.adjacentEdges()) {
+      int iV = bIdx[src(e)];
+      GC_SAFETY_ASSERT(0 <= iV && iV < (int)nBoundary, "invalid boundary vertex index");
+
+      targetLength(iV) = geo.edgeLengths[e] * exp(0.5 * (uBdy(bIdx[src(e)]) + uBdy(bIdx[dst(e)])));
+
+      T(0, iV) = cos(phi);
+      T(1, iV) = sin(phi);
+
+      Ntriplets.emplace_back(iV, iV, geo.edgeLengths[e]);
+
+      phi += kBdy(bIdx[dst(e)]);
+    }
+  }
+
+  SparseMatrix<double> Ninv(nBoundary, nBoundary);
+  Ninv.setFromTriplets(std::begin(Ntriplets), std::end(Ntriplets));
+
+  Vector<double> roundedLength =
+      targetLength - Ninv * T.transpose() * (T * Ninv * T.transpose()).inverse() * T * targetLength;
+
+  std::array<Vector<double>, 2> bdyPositions{Vector<double>(nBoundary), Vector<double>(nBoundary)};
+  double x = 0, y = 0;
+  for (BoundaryLoop b : mesh.boundaryLoops()) {
+    for (Edge e : b.adjacentEdges()) {
+      size_t iV = bIdx[src(e)];
+
+      bdyPositions[0](iV) = x;
+      bdyPositions[1](iV) = y;
+
+      x += roundedLength(iV) * T(0, iV);
+      y += roundedLength(iV) * T(1, iV);
+    }
+  }
+
+  return bdyPositions;
+}
+
+void BFF::ensureHaveLSolver() {
+  if (!Lsolver) {
+    Lsolver.reset(new PositiveDefiniteSolver<double>(L));
+  }
+}
+} // namespace surface
+} // namespace geometrycentral