API cleanup

DavidJourdan · Jul 4, 2024 · 8de1cd2 · 8de1cd2
1 parent 5543cae
commit 8de1cd2
Show file tree

Hide file tree

Showing 11 changed files with 107 additions and 68 deletions.
diff --git a/src/LocalGlobalSolver.h b/src/LocalGlobalSolver.h
@@ -3,7 +3,7 @@
 #include "solvers.h"
 
 #include <Eigen/Core>
-#include <Eigen/Sparse>
+#include <Eigen/SparseCore>
 #include <memory>
 
 class LocalGlobalSolver
@@ -18,8 +18,6 @@ class LocalGlobalSolver
 
   Eigen::VectorXd stretchAngles();
 
-  void conservativeResize(int _nF);
-
   Eigen::VectorXd s1;
   Eigen::VectorXd s2;
   Eigen::MatrixX2d stressX;

diff --git a/src/cli.cpp b/src/cli.cpp
@@ -78,8 +78,7 @@ int main(int argc, char* argv[])
 
   // Run local-global parameterization algorithm
   Timer paramTimer("Parameterization");
-  Eigen::MatrixXd P;
-  LocalGlobalSolver LGsolver = localGlobal(V, F, P, lambda1, lambda2);
+  Eigen::MatrixXd P = localGlobal(V, F, lambda1, lambda2);
   paramTimer.stop();
 
   if(wD > 0) // smoothing

diff --git a/src/main.cpp b/src/main.cpp
@@ -8,10 +8,10 @@
 #include "newton.h"
 #include "parameterization.h"
 #include "path_extraction.h"
+#include "save.h"
 #include "simulation_utils.h"
 #include "stretch_angles.h"
 #include "timer.h"
-#include "save.h"
 
 #include <geometrycentral/surface/manifold_surface_mesh.h>
 #include <geometrycentral/surface/vertex_position_geometry.h>
@@ -24,8 +24,6 @@
 
 #include <thread>
 
-
-
 int main(int argc, char* argv[])
 {
   using namespace geometrycentral;
@@ -84,16 +82,17 @@ int main(int argc, char* argv[])
   // Run local-global parameterization algorithm
   Timer paramTimer("Parameterization");
 
-  Eigen::MatrixXd P;
-  LocalGlobalSolver LGsolver = localGlobal(V, F, P, lambda1, lambda2);
+  Eigen::MatrixXd P = localGlobal(V, F, lambda1, lambda2);
   paramTimer.stop();
 
+  const auto &[sigma1, sigma2, angles] = computeSVDdata(V, P, F);
+
   // Add mesh and initial param to polyscope viewer
   polyscope::registerSurfaceMesh2D("Param", P, F)->setEdgeWidth(0.0);
-  polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma1", LGsolver.s1.head(F.rows()));
-  polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma2", LGsolver.s2.head(F.rows()));
+  polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma1", sigma1);
+  polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma2", sigma2);
   polyscope::getSurfaceMesh("Param")
-      ->addFaceScalarQuantity("stretch orientation", LGsolver.stretchAngles().head(F.rows()))
+      ->addFaceScalarQuantity("stretch orientation", angles)
       ->setColorMap("twilight")
       ->setMapRange({-PI / 2, PI / 2})
       ->setEnabled(true);
@@ -126,14 +125,12 @@ int main(int argc, char* argv[])
       polyscope::getSurfaceMesh("Param")->updateVertexPositions2D(P);
 
       // compute data from SVDs and display them
-      LGsolver.solveOneStep(P, 1 / lambda2, 1 / lambda1);
-      // center P
-      P.rowwise() -= P.colwise().sum() / P.rows();
+      const auto &[sigma1, sigma2, angles] = computeSVDdata(V, P, F);
 
-      polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma1", LGsolver.s1.head(F.rows()));
-      polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma2", LGsolver.s2.head(F.rows()));
+      polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma1", sigma1);
+      polyscope::getSurfaceMesh("Param")->addFaceScalarQuantity("sigma2", sigma2);
       polyscope::getSurfaceMesh("Param")
-          ->addFaceScalarQuantity("stretch orientation", LGsolver.stretchAngles().head(F.rows()))
+          ->addFaceScalarQuantity("stretch orientation", angles)
           ->setColorMap("twilight")
           ->setMapRange({-PI / 2, PI / 2})
           ->setEnabled(true);
@@ -145,32 +142,39 @@ int main(int argc, char* argv[])
 
     if(ImGui::TreeNode("Advanced"))
     {
-      if(ImGui::InputDouble("lambda1", &lambda1, 0, 0, "%.1e")){
+      if(ImGui::InputDouble("lambda1", &lambda1, 0, 0, "%.1e"))
+      {
         updateAndSaveAdvancedParams(lambda1, lambda2, thickness, deltaLambda, nLayers, nIter, lim);
       }
-      if(ImGui::InputDouble("lambda2", &lambda2, 0, 0, "%.1e")){
+      if(ImGui::InputDouble("lambda2", &lambda2, 0, 0, "%.1e"))
+      {
         updateAndSaveAdvancedParams(lambda1, lambda2, thickness, deltaLambda, nLayers, nIter, lim);
       }
-      if(ImGui::InputDouble("Thickness", &thickness, 0, 0, "%.1e")){
+      if(ImGui::InputDouble("Thickness", &thickness, 0, 0, "%.1e"))
+      {
         updateAndSaveAdvancedParams(lambda1, lambda2, thickness, deltaLambda, nLayers, nIter, lim);
       }
-      //if (ImGui::InputDouble("Curvature", &curvature, 0, 0, "%.1e"))
+      // if (ImGui::InputDouble("Curvature", &curvature, 0, 0, "%.1e"))
       //{
-      //  deltaLambda = thickness * lambda1 * curvature;
-      //  updateAndSaveAdvancedParams(lambda1, lambda2, thickness, deltaLambda, n_layers, nIter, lim);
-      //}
-      if(ImGui::InputInt("Nb of layers",&nLayers)){
+      //   deltaLambda = thickness * lambda1 * curvature;
+      //   updateAndSaveAdvancedParams(lambda1, lambda2, thickness, deltaLambda, n_layers, nIter, lim);
+      // }
+      if(ImGui::InputInt("Nb of layers", &nLayers))
+      {
         updateAndSaveAdvancedParams(lambda1, lambda2, thickness, deltaLambda, nLayers, nIter, lim);
       }
-      if(ImGui::InputInt("Iterations", &nIter)){
+      if(ImGui::InputInt("Iterations", &nIter))
+      {
         updateAndSaveAdvancedParams(lambda1, lambda2, thickness, deltaLambda, nLayers, nIter, lim);
       }
-      if(ImGui::InputDouble("Limit", &lim, 0, 0, "%.1e")){
+      if(ImGui::InputDouble("Limit", &lim, 0, 0, "%.1e"))
+      {
         updateAndSaveAdvancedParams(lambda1, lambda2, thickness, deltaLambda, nLayers, nIter, lim);
       }
       if(ImGui::Button("Rotate"))
       {
         // run ARAP
+        LocalGlobalSolver LGsolver(V, F);
         LGsolver.solve(P, 1., 1.);
 
         // center and rotate

diff --git a/src/newton.h b/src/newton.h
@@ -1,6 +1,6 @@
 #pragma once
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 #include <TinyAD/Support/GeometryCentral.hh>
 #include <TinyAD/ScalarFunction.hh>
 #include <geometrycentral/surface/intrinsic_geometry_interface.h>

diff --git a/src/parameterization.cpp b/src/parameterization.cpp
@@ -7,8 +7,10 @@
 #include <igl/harmonic.h>
 #include <igl/loop.h>
 #include <igl/map_vertices_to_circle.h>
+#include <igl/project_isometrically_to_plane.h>
 
-Eigen::MatrixXd tutteEmbedding(const Eigen::MatrixXd &_V, const Eigen::MatrixXi &_F, const std::vector<int> &boundary_indices)
+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
@@ -31,11 +33,7 @@ Eigen::MatrixXd tutteEmbedding(const Eigen::MatrixXd &_V, const Eigen::MatrixXi
   return P;
 }
 
-LocalGlobalSolver localGlobal(const Eigen::MatrixXd &V,
-                            const Eigen::MatrixXi &F,
-                            Eigen::MatrixXd &P,
-                            double lambda1,
-                            double lambda2)
+Eigen::MatrixXd localGlobal(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, double lambda1, double lambda2)
 {
   using namespace Eigen;
 
@@ -74,7 +72,7 @@ LocalGlobalSolver localGlobal(const Eigen::MatrixXd &V,
 
   LocalGlobalSolver solver(V, _F);
 
-  P = tutteEmbedding(V, F, loops[idxMax]);
+  MatrixXd P = tutteEmbedding(V, F, loops[idxMax]);
 
   // run ARAP
   solver.solve(P, 1., 1.);
@@ -123,21 +121,21 @@ LocalGlobalSolver localGlobal(const Eigen::MatrixXd &V,
   // reinit with original topology
   solver.init(V, F);
   solver.solveOneStep(P, 1 / lambda2, 1 / lambda1);
-  return solver;
+  return P;
 }
 
 geometrycentral::surface::FaceData<Eigen::Matrix2d>
-precomputeParamData(geometrycentral::surface::VertexPositionGeometry &geometry)
+precomputeParamData(geometrycentral::surface::VertexPositionGeometry& geometry)
 {
   using namespace geometrycentral;
   using namespace geometrycentral::surface;
 
-  SurfaceMesh &mesh = geometry.mesh;
+  SurfaceMesh& mesh = geometry.mesh;
 
   FaceData<Eigen::Matrix2d> rest_shapes(mesh);
   geometry.requireVertexPositions();
 
-  auto toEigen = [](const geometrycentral::Vector3 &_v) { return Eigen::Vector3d(_v.x, _v.y, _v.z); };
+  auto toEigen = [](const geometrycentral::Vector3& _v) { return Eigen::Vector3d(_v.x, _v.y, _v.z); };
 
   for(Face f: mesh.faces())
   {
@@ -163,9 +161,9 @@ precomputeParamData(geometrycentral::surface::VertexPositionGeometry &geometry)
 }
 
 geometrycentral::surface::FaceData<Eigen::Matrix2d>
-precomputeSimData(geometrycentral::surface::ManifoldSurfaceMesh &mesh,
-                  const Eigen::MatrixXd &P,
-                  const Eigen::MatrixXi &F)
+precomputeSimData(geometrycentral::surface::ManifoldSurfaceMesh& mesh,
+                  const Eigen::MatrixXd& P,
+                  const Eigen::MatrixXi& F)
 {
   using namespace geometrycentral;
   using namespace geometrycentral::surface;
@@ -187,12 +185,12 @@ precomputeSimData(geometrycentral::surface::ManifoldSurfaceMesh &mesh,
 }
 
 geometrycentral::surface::EdgeData<double>
-computeDualCotanWeights(geometrycentral::surface::IntrinsicGeometryInterface &geometry)
+computeDualCotanWeights(geometrycentral::surface::IntrinsicGeometryInterface& geometry)
 {
   using namespace geometrycentral;
   using namespace geometrycentral::surface;
 
-  SurfaceMesh &mesh = geometry.mesh;
+  SurfaceMesh& mesh = geometry.mesh;
 
   geometry.requireEdgeLengths();
   geometry.requireHalfedgeCotanWeights();
@@ -235,17 +233,17 @@ computeDualCotanWeights(geometrycentral::surface::IntrinsicGeometryInterface &ge
   return dualCotanWeights;
 }
 
-void subdivideMesh(geometrycentral::surface::VertexPositionGeometry &geometry,
-                   Eigen::MatrixXd &V,
-                   Eigen::MatrixXd &P,
-                   Eigen::MatrixXi &F,
-                   std::vector<Eigen::SparseMatrix<double>> &subdivMat,
+void subdivideMesh(geometrycentral::surface::VertexPositionGeometry& geometry,
+                   Eigen::MatrixXd& V,
+                   Eigen::MatrixXd& P,
+                   Eigen::MatrixXi& F,
+                   std::vector<Eigen::SparseMatrix<double>>& subdivMat,
                    double threshold)
 {
   using namespace geometrycentral;
   using namespace geometrycentral::surface;
 
-  SurfaceMesh &mesh = geometry.mesh;
+  SurfaceMesh& mesh = geometry.mesh;
 
   double maxEdgeLength = 0;
   geometry.requireVertexPositions();
@@ -269,15 +267,15 @@ void subdivideMesh(geometrycentral::surface::VertexPositionGeometry &geometry,
 }
 
 TinyAD::ScalarFunction<2, double, geometrycentral::surface::Vertex>
-parameterizationFunction(geometrycentral::surface::VertexPositionGeometry &geometry,
+parameterizationFunction(geometrycentral::surface::VertexPositionGeometry& geometry,
                          double wPhi,
                          double lambda1,
                          double lambda2)
 {
   using namespace geometrycentral;
   using namespace geometrycentral::surface;
 
-  SurfaceMesh &mesh = geometry.mesh;
+  SurfaceMesh& mesh = geometry.mesh;
 
   // precompute inverse jacobians
   FaceData<Eigen::Matrix2d> restShapes = precomputeParamData(geometry);
@@ -286,7 +284,7 @@ parameterizationFunction(geometrycentral::surface::VertexPositionGeometry &geome
   TinyAD::ScalarFunction<2, double, Vertex> func = TinyAD::scalar_function<2>(mesh.vertices());
 
   // Main objective term
-  func.add_elements<3>(mesh.faces(), [&, lambda1, lambda2, restShapes](auto &element) -> TINYAD_SCALAR_TYPE(element) {
+  func.add_elements<3>(mesh.faces(), [&, lambda1, lambda2, restShapes](auto& element) -> TINYAD_SCALAR_TYPE(element) {
     // Evaluate element using either double or TinyAD::Double
     using T = TINYAD_SCALAR_TYPE(element);
 
@@ -320,7 +318,7 @@ parameterizationFunction(geometrycentral::surface::VertexPositionGeometry &geome
   // Smoothness regularization (Dirichlet energy)
   func.add_elements<4>(
       mesh.edges(),
-      [&, lambda1, lambda2, wPhi, dualCotanWeights, restShapes](auto &element) -> TINYAD_SCALAR_TYPE(element) {
+      [&, lambda1, lambda2, wPhi, dualCotanWeights, restShapes](auto& element) -> TINYAD_SCALAR_TYPE(element) {
         // Evaluate element using either double or TinyAD::Double
         using T = TINYAD_SCALAR_TYPE(element);
 
@@ -355,3 +353,44 @@ parameterizationFunction(geometrycentral::surface::VertexPositionGeometry &geome
 
   return func;
 }
+
+std::tuple<Eigen::VectorXd, Eigen::VectorXd, Eigen::VectorXd>
+computeSVDdata(const Eigen::MatrixXd& V, const Eigen::MatrixXd& P, const Eigen::MatrixXi& F)
+{
+  using namespace Eigen;
+  int nF = F.rows();
+
+  MatrixX2d plane_V;
+  MatrixX3i plane_F;
+  SparseMatrix<double> ref_map;
+  igl::project_isometrically_to_plane(V, F, plane_V, plane_F, ref_map);
+
+  VectorXd sigma1(nF), sigma2(nF), angles(nF);
+  #pragma omp parallel for schedule(static) num_threads(omp_get_max_threads() - 1)
+  for(int i = 0; i < nF; ++i)
+  {
+    Matrix2d A;
+    A.col(0) = plane_V.row(nF + i) - plane_V.row(i);
+    A.col(1) = plane_V.row(2 * nF + i) - plane_V.row(i);
+
+    Matrix2d B;
+    B.col(0) = P.row(F(i, 1)) - P.row(F(i, 0));
+    B.col(1) = P.row(F(i, 2)) - P.row(F(i, 0));
+
+    // compute SVD of the V -> P parameterization
+    JacobiSVD<Matrix2d> svd;
+    Matrix2d M = B * A.inverse();
+    svd.compute(M, ComputeFullU | ComputeFullV);
+    Matrix2d S = svd.singularValues().asDiagonal();
+
+    // save the singular values sigma1, sigma2
+    sigma1(i) = S(0, 0);
+    sigma2(i) = S(1, 1);
+
+    // compute angles from the U, V matrices
+    Vector2d stressX = svd.matrixU().col(0);
+    geometrycentral::Vector2 dir{stressX(0), stressX(1)};
+    angles(i) = dir.pow(2).arg() / 2;
+  }
+  return std::make_tuple(sigma1, sigma2, angles);
+}
diff --git a/src/parameterization.h b/src/parameterization.h
@@ -11,11 +11,7 @@
 #include <geometrycentral/surface/manifold_surface_mesh.h>
 #include <geometrycentral/surface/vertex_position_geometry.h>
 
-LocalGlobalSolver localGlobal(const Eigen::MatrixXd &V,
-                            const Eigen::MatrixXi &F,
-                            Eigen::MatrixXd &P,
-                            double lambda1,
-                            double lambda2);
+Eigen::MatrixXd localGlobal(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, double lambda1, double lambda2);
 
 geometrycentral::surface::FaceData<Eigen::Matrix2d>
 precomputeParamData(geometrycentral::surface::VertexPositionGeometry& geometry);
@@ -46,3 +42,6 @@ parameterizationFunction(geometrycentral::surface::VertexPositionGeometry& geome
                          double wPhi,
                          double lambda1,
                          double lambda2);
+
+std::tuple<Eigen::VectorXd, Eigen::VectorXd, Eigen::VectorXd>
+computeSVDdata(const Eigen::MatrixXd& V, const Eigen::MatrixXd& P, const Eigen::MatrixXi& F);
diff --git a/src/path_extraction.cpp b/src/path_extraction.cpp
@@ -6,9 +6,9 @@
 #include "stripe_patterns.h"
 #include "timer.h"
 
+#include <geometrycentral/numerical/linear_solvers.h>
 #include <geometrycentral/surface/direction_fields.h>
 #include <geometrycentral/surface/manifold_surface_mesh.h>
-#include <geometrycentral/numerical/linear_solvers.h>
 #include <igl/boundary_loop.h>
 #include <igl/ramer_douglas_peucker.h>
 #include <polyscope/curve_network.h>
@@ -320,12 +320,12 @@ void stripePattern(VertexPositionGeometry& geometry,
   std::vector<std::vector<std::vector<Vector3>>> paths =
       generatePaths(geometryUV, th1, th2, layerHeight, nLayers, spacing, timeLimit);
 
-  //Update height every layers
+  // Update height every layers
   static double height = 0;
   for(int i = 0; i < nLayers; ++i)
   {
     height += layerHeight + static_cast<float>(i) / (nLayers - 1) * 2 * (0.8 / nLayers - layerHeight);
-    writePaths(filename + ".path", paths[i],height);
+    writePaths(filename + ".path", paths[i], height);
     drawPathsAndTravels(paths[i], spacing, i + 1);
   }
 }

diff --git a/src/simulation_utils.cpp b/src/simulation_utils.cpp
@@ -1,6 +1,6 @@
 #include "simulation_utils.h"
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 #include <igl/colon.h>
 #include <igl/slice.h>
 #include <igl/slice_into.h>