manage meshes with holes

src/LocalGlobalSolver.cpp

@@ -8,14 +8,20 @@
 #include <igl/project_isometrically_to_plane.h>
 #include <igl/repdiag.h>
 
-LocalGlobalSolver::LocalGlobalSolver(const Eigen::Ref<Eigen::MatrixX3d> V, const Eigen::Ref<Eigen::MatrixX3i> F) : _F(F)
+LocalGlobalSolver::LocalGlobalSolver(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F)
+{
+  init(V, F);
+}
+
+void LocalGlobalSolver::init(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F)
 {
   using namespace Eigen;
 
   // number of vertices
   nV = V.rows();
   nF = F.rows();
 
+  _F = F;
   s1.resize(nF);
   s2.resize(nF);
   stressX.resize(nF, 2);
@@ -60,7 +66,7 @@ void LocalGlobalSolver::solveOneStep(Eigen::Ref<Eigen::MatrixX2d> U, double sMin
 
   Matrix<double, 2, -1> R(2, 2 * nF);
 #pragma omp parallel for schedule(static) num_threads(omp_get_max_threads() - 1)
-  for(int i = 0; i < _F.rows(); ++i)
+  for(int i = 0; i < nF; ++i)
   {
     Matrix2d B;
     B.col(0) = U.row(_F(i, 1)) - U.row(_F(i, 0));

src/LocalGlobalSolver.h

@@ -7,14 +7,17 @@
 class LocalGlobalSolver
 {
 public:
-  LocalGlobalSolver(const Eigen::Ref<Eigen::MatrixX3d> V, const Eigen::Ref<Eigen::MatrixX3i> F);
+  LocalGlobalSolver(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F);
+  void init(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F);
 
   void solveOneStep(Eigen::Ref<Eigen::MatrixX2d> U, double sMin, double sMax);
 
   void solve(Eigen::Ref<Eigen::MatrixX2d> U, double sMin, double sMax, int nbIter = -1);
 
   Eigen::VectorXd stretchAngles();
 
+  void conservativeResize(int _nF);
+
   Eigen::VectorXd s1;
   Eigen::VectorXd s2;
   Eigen::MatrixX2d stressX;

src/cli.cpp

@@ -16,14 +16,15 @@
 #include <geometrycentral/surface/vertex_position_geometry.h>
 #include <igl/loop.h>
 #include <igl/readOBJ.h>
+
 #include <thread>
 
 int main(int argc, char* argv[])
 {
   using namespace geometrycentral;
   using namespace geometrycentral::surface;
 
-  std::string filename = "beetle";
+  std::string filename = "beetle.obj";
   double wD = 0;
   double width = 200;
   double lim = 1e-6;
@@ -54,9 +55,9 @@ int main(int argc, char* argv[])
   // Load a mesh in OBJ format
   Eigen::MatrixXd V;
   Eigen::MatrixXi F;
-  if(!igl::readOBJ(std::string(DATA_PATH_STR) + filename + ".obj", V, F))
+  if(!igl::readOBJ(filename, V, F))
   {
-    std::cout << "File " << DATA_PATH_STR << filename << ".obj not found\n";
+    std::cout << "File " << filename << " not found\n";
     return 0;
   }
   // resize mesh to a 100mm width
@@ -77,8 +78,8 @@ int main(int argc, char* argv[])
 
   // Run local-global parameterization algorithm
   Timer paramTimer("Parameterization");
-  LocalGlobalSolver LGsolver(V, F);
-  Eigen::MatrixXd P = localGlobal(V, F, lambda1, lambda2, LGsolver);
+  Eigen::MatrixXd P;
+  LocalGlobalSolver LGsolver = localGlobal(V, F, P, lambda1, lambda2);
   paramTimer.stop();
 
   if(wD > 0) // smoothing
@@ -115,7 +116,6 @@ int main(int argc, char* argv[])
     for(int j = 0; j < 3; ++j)
       xTarget(3 * i + j) = V(i, j);
 
-
   // Simulation
   std::cout << "**********\nSIMULATION\n**********\n";
 
@@ -132,11 +132,9 @@ int main(int argc, char* argv[])
   Eigen::VectorXd d = (V - VTarget).cwiseProduct((V - VTarget)).rowwise().sum();
   d = d.array().sqrt();
   std::cout << "Avg distance = "
-            << 100 * d.sum() / d.size() / (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm()
-            << "\n";
+            << 100 * d.sum() / d.size() / (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm() << "\n";
   std::cout << "Max distance = "
-            << 100 * d.lpNorm<Eigen::Infinity>() /
-                  (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm()
+            << 100 * d.lpNorm<Eigen::Infinity>() / (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm()
             << "\n";
 
   // Directions optimizatiom
@@ -154,11 +152,9 @@ int main(int argc, char* argv[])
   d = (V - VTarget).cwiseProduct((V - VTarget)).rowwise().sum();
   d = d.array().sqrt();
   std::cout << "Avg distance = "
-            << 100 * d.sum() / d.size() / (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm()
-            << "\n";
+            << 100 * d.sum() / d.size() / (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm() << "\n";
   std::cout << "Max distance = "
-            << 100 * d.lpNorm<Eigen::Infinity>() /
-                  (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm()
+            << 100 * d.lpNorm<Eigen::Infinity>() / (VTarget.colwise().maxCoeff() - VTarget.colwise().minCoeff()).norm()
             << "\n";
 
   // Generate trajectories
@@ -212,9 +208,6 @@ int main(int argc, char* argv[])
   std::vector<std::vector<std::vector<Vector3>>> paths =
       generatePaths(geometryUV, th1, th2, layerHeight, nLayers, spacing, timeLimit);
 
-  std::ofstream s(std::string(DATA_PATH_STR) + filename + ".path");
   for(int i = 0; i < nLayers; ++i)
-  {
-    writePaths(std::string(DATA_PATH_STR) + filename + ".path", paths[i], (i + 1) * layerHeight);
-  }
+    writePaths(filename + ".path", paths[i], (i + 1) * layerHeight);
 }

src/main.cpp

@@ -34,7 +34,7 @@ int main(int argc, char* argv[])
   polyscope::view::style = polyscope::view::NavigateStyle::Planar;
   polyscope::options::groundPlaneMode = polyscope::GroundPlaneMode::ShadowOnly;
   polyscope::options::openImGuiWindowForUserCallback = false;
-  polyscope::loadColorMap("twilight", DATA_PATH_STR "twilight_colormap.png");
+  polyscope::loadColorMap("twilight", "twilight_colormap.png");
 
   // Some state about imgui windows to stack them
   float imguiStackMargin = 10;
@@ -51,7 +51,7 @@ int main(int argc, char* argv[])
   if(argc < 2)
     filename = igl::file_dialog_open();
   else
-    filename = std::string(DATA_PATH_STR) + std::string(argv[1]) + ".obj";
+    filename = argv[1];
   Eigen::MatrixXd V;
   Eigen::MatrixXi F;
   if(!igl::readOBJ(filename, V, F))
@@ -80,16 +80,17 @@ int main(int argc, char* argv[])
 
   // Run local-global parameterization algorithm
   Timer paramTimer("Parameterization");
-  LocalGlobalSolver LGsolver(V, F);
-  Eigen::MatrixXd P = localGlobal(V, F, lambda1, lambda2, LGsolver);
+
+  Eigen::MatrixXd P;
+  LocalGlobalSolver LGsolver = localGlobal(V, F, P, lambda1, lambda2);
   paramTimer.stop();
 
   // Add mesh and initial param to polyscope viewer
   polyscope::registerSurfaceMesh2D("Param", P, F)->setEdgeWidth(0.0);
-  polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma1", LGsolver.s1);
-  polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma2", LGsolver.s2);
+  polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma1", LGsolver.s1.head(F.rows()));
+  polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma2", LGsolver.s2.head(F.rows()));
   polyscope::getSurfaceMesh("Param")
-      ->addFaceScalarQuantity("stretch orientation", LGsolver.stretchAngles())
+      ->addFaceScalarQuantity("stretch orientation", LGsolver.stretchAngles().head(F.rows()))
       ->setColorMap("twilight")
       ->setMapRange({-PI / 2, PI / 2})
       ->setEnabled(true);
@@ -129,10 +130,10 @@ int main(int argc, char* argv[])
       // center P
       P.rowwise() -= P.colwise().sum() / P.rows();
 
-      polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma1", LGsolver.s1);
-      polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma2", LGsolver.s2);
+      polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma1", LGsolver.s1.head(F.rows()));
+      polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma2", LGsolver.s2.head(F.rows()));
       polyscope::getSurfaceMesh("Param")
-          ->addFaceScalarQuantity("stretch orientation", LGsolver.stretchAngles())
+          ->addFaceScalarQuantity("stretch orientation", LGsolver.stretchAngles().head(F.rows()))
           ->setColorMap("twilight")
           ->setMapRange({-PI / 2, PI / 2})
           ->setEnabled(true);
@@ -147,7 +148,7 @@ int main(int argc, char* argv[])
       ImGui::InputDouble("lambda1", &lambda1, 0, 0, "%.1e");
       ImGui::InputDouble("lambda2", &lambda2, 0, 0, "%.1e");
       ImGui::InputDouble("Thickness", &thickness, 0, 0, "%.1e");
-      if (ImGui::InputDouble("Curvature", &curvature, 0, 0, "%.1e"))
+      if(ImGui::InputDouble("Curvature", &curvature, 0, 0, "%.1e"))
       {
         deltaLambda = thickness * lambda1 * curvature;
       }
@@ -336,9 +337,8 @@ int main(int argc, char* argv[])
         auto adjointFunc = adjointFunction(geometry, F, MrInv, theta1, E1, lambda1, lambda2, deltaLambda, thickness);
 
         // Optimize this energy function using SGN [Zehnder et al. 2021]
-        sparse_gauss_newton(
-            geometry, V, xTarget, MrInv, theta1, theta2, adjointFunc, fixedIdx, n_iter, lim, wM, wL, E1, lambda1,
-            lambda2, deltaLambda, thickness, [&](const Eigen::VectorXd& X) {
+        sparse_gauss_newton(geometry, V, xTarget, MrInv, theta1, theta2, adjointFunc, fixedIdx, n_iter, lim, wM, wL, E1,
+                            lambda1, lambda2, deltaLambda, thickness, [&](const Eigen::VectorXd& X) {
                               // convert X to V for visualizing individual iterations
                               for(int i = 0; i < V.rows(); ++i)
                                 for(int j = 0; j < 3; ++j)

src/parameterization.cpp

@@ -8,16 +8,42 @@
 #include <igl/loop.h>
 #include <igl/map_vertices_to_circle.h>
 
-Eigen::MatrixXd tutte_embedding(const Eigen::MatrixXd &_V, const Eigen::MatrixXi &_F)
+Eigen::MatrixXd tutteEmbedding(const Eigen::MatrixXd &_V, const Eigen::MatrixXi &_F, const std::vector<int> &boundary_indices)
 {
   Eigen::VectorXi b;  // #constr boundary constraint indices
   Eigen::MatrixXd bc; // #constr-by-2 2D boundary constraint positions
   Eigen::MatrixXd P;  // #V-by-2 2D vertex positions
   Eigen::MatrixXi F = _F;
 
+  // Set boundary vertex positions
+  b.resize(boundary_indices.size(), 1);
+  for(size_t i = 0; i < boundary_indices.size(); ++i)
+    b(i) = boundary_indices[i];
+  igl::map_vertices_to_circle(_V, b, bc);
+
+  // make sure normals are consistent
+  Eigen::VectorXi C;
+  igl::bfs_orient(F, F, C);
+
+  // Compute interior vertex positions
+  igl::harmonic(F, b, bc, 1, P);
+
+  return P;
+}
+
+LocalGlobalSolver localGlobal(const Eigen::MatrixXd &V,
+                            const Eigen::MatrixXi &F,
+                            Eigen::MatrixXd &P,
+                            double lambda1,
+                            double lambda2)
+{
+  using namespace Eigen;
+
+  Eigen::MatrixXi _F = F;
+
   // Identify boundary loops
   std::vector<std::vector<int>> loops;
-  igl::boundary_loop(F, loops);
+  igl::boundary_loop(_F, loops);
 
   // find longest boundary loop
   size_t idxMax = -1;
@@ -26,7 +52,7 @@ Eigen::MatrixXd tutte_embedding(const Eigen::MatrixXd &_V, const Eigen::MatrixXi
   {
     double len = 0;
     for(size_t j = 0; j < loops[i].size(); ++j)
-      len += (_V.row(loops[i][(j + 1) % loops[i].size()]) - _V.row(loops[i][j])).norm();
+      len += (V.row(loops[i][(j + 1) % loops[i].size()]) - V.row(loops[i][j])).norm();
 
     if(len > maxLen)
     {
@@ -35,42 +61,20 @@ Eigen::MatrixXd tutte_embedding(const Eigen::MatrixXd &_V, const Eigen::MatrixXi
     }
   }
 
-  // Set boundary vertex positions
-  b.resize(loops[idxMax].size(), 1);
-  for(size_t i = 0; i < loops[idxMax].size(); ++i)
-    b(i) = loops[idxMax][i];
-  igl::map_vertices_to_circle(_V, b, bc);
-
   // Fill-in other holes
   for(size_t i = 0; i < loops.size(); ++i)
     if(i != idxMax)
     {
       int nB = loops[i].size();
-      int nF = F.rows();
-      F.conservativeResize(nF + nB - 2, 3);
+      int nF = _F.rows();
+      _F.conservativeResize(nF + nB - 2, 3);
       for(int j = 0; j < nB - 2; ++j)
-        F.row(nF + j) << loops[i][0], loops[i][j + 1], loops[i][j + 2];
+        _F.row(nF + j) << loops[i][0], loops[i][j + 1], loops[i][j + 2];
     }
 
-  // make sure normals are consistent
-  Eigen::VectorXi C;
-  igl::bfs_orient(F, F, C);
-
-  // Compute interior vertex positions
-  igl::harmonic(F, b, bc, 1, P);
-
-  return P;
-}
+  LocalGlobalSolver solver(V, _F);
 
-Eigen::MatrixXd localGlobal(const Eigen::MatrixXd &V,
-                            const Eigen::MatrixXi &F,
-                            double lambda1,
-                            double lambda2,
-                            LocalGlobalSolver &solver)
-{
-  using namespace Eigen;
-
-  MatrixXd P = tutte_embedding(V, F);
+  P = tutteEmbedding(V, F, loops[idxMax]);
 
   // run ARAP
   solver.solve(P, 1., 1.);
@@ -116,7 +120,10 @@ Eigen::MatrixXd localGlobal(const Eigen::MatrixXd &V,
   R(0, 1) = -R(1, 0);
   P = P * R;
 
-  return P;
+  // reinit with original topology
+  solver.init(V, F);
+  solver.solveOneStep(P, 1 / lambda2, 1 / lambda1);
+  return solver;
 }
 
 geometrycentral::surface::FaceData<Eigen::Matrix2d>
@@ -240,16 +247,13 @@ void subdivideMesh(geometrycentral::surface::VertexPositionGeometry &geometry,
 
   SurfaceMesh &mesh = geometry.mesh;
 
-  double avgEdgeLength = 0;
+  double maxEdgeLength = 0;
   geometry.requireVertexPositions();
   for(Edge e: mesh.edges())
-    //     avgEdgeLength += norm(geometry.vertexPositions[e.firstVertex()] -
-    //     geometry.vertexPositions[e.secondVertex()]);
-    // avgEdgeLength /= mesh.nEdges();
-    if(norm(geometry.vertexPositions[e.firstVertex()] - geometry.vertexPositions[e.secondVertex()]) > avgEdgeLength)
-      avgEdgeLength = norm(geometry.vertexPositions[e.firstVertex()] - geometry.vertexPositions[e.secondVertex()]);
+    if(norm(geometry.vertexPositions[e.firstVertex()] - geometry.vertexPositions[e.secondVertex()]) > maxEdgeLength)
+      maxEdgeLength = norm(geometry.vertexPositions[e.firstVertex()] - geometry.vertexPositions[e.secondVertex()]);
 
-  while(avgEdgeLength > threshold)
+  while(maxEdgeLength > threshold)
   {
     Eigen::MatrixX3i tempF = F;
     Eigen::SparseMatrix<double> S;
@@ -260,7 +264,7 @@ void subdivideMesh(geometrycentral::surface::VertexPositionGeometry &geometry,
 
     subdivMat.push_back(S);
 
-    avgEdgeLength /= 2;
+    maxEdgeLength /= 2;
   }
 }
 

src/parameterization.h

@@ -11,11 +11,11 @@
 #include <geometrycentral/surface/manifold_surface_mesh.h>
 #include <geometrycentral/surface/vertex_position_geometry.h>
 
-Eigen::MatrixXd localGlobal(const Eigen::MatrixXd& V,
-                            const Eigen::MatrixXi& F,
+LocalGlobalSolver localGlobal(const Eigen::MatrixXd &V,
+                            const Eigen::MatrixXi &F,
+                            Eigen::MatrixXd &P,
                             double lambda1,
-                            double lambda2,
-                            LocalGlobalSolver& solver);
+                            double lambda2);
 
 geometrycentral::surface::FaceData<Eigen::Matrix2d>
 precomputeParamData(geometrycentral::surface::VertexPositionGeometry& geometry);