Merged feature into main

Author: rjc8237

CMakeLists.txt

@@ -15,7 +15,6 @@ if (PROJECT_IS_TOP_LEVEL)
 
   # Alternative library options
   option(USE_MULTIPRECISION "Use high precision floating point" OFF)
-  option(USE_EMBREE "Use Embree for rendering" OFF)
   option(USE_SUITESPARSE "Use suite sparse methods for matrix inversion" ON)
   option(USE_COMISO "Use Comiso for field generation" ON)
 
@@ -35,11 +34,11 @@ endif()
 
 list(PREPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
 
-#include(boost)
+include(boost)
 #include(eigen)
 
 # Optionally get multiprecision libraries
-if(USE_MULTIPRECISION)
+if(USE_MULTIPRECISION OR USE_MULTIPRECISION_OVERLAY)
   # Turning this on activates mpfr in the conformal code
   option(LIBIGL_COPYLEFT_CGAL "Use CGAL" ON)
   add_compile_definitions(WITH_MPFR)
@@ -48,6 +47,9 @@ if(USE_MULTIPRECISION)
   include_directories(${GMP_INCLUDE_DIR})
   link_directories(${MPFR_LIBRARIES_DIR})
 endif()
+if(USE_MULTIPRECISION) 
+  add_compile_definitions(MULTIPRECISION)
+endif()
 
 # Optionally get frame field libraries
 if (USE_COMISO)
@@ -64,9 +66,6 @@ if (RENDER_TEXTURE)
   option(LIBIGL_GLFW "Use GLFW" ON)
   option(LIBIGL_PNG "Use PNG" ON)
 endif()
-if (USE_EMBREE)
-  option(LIBIGL_EMBREE "Use EMBREE" ON)
-endif()
 
 # Set compile definitions
 add_compile_definitions(_USE_MATH_DEFINES)
@@ -79,9 +78,6 @@ endif()
 if (RENDER_TEXTURE)
   add_compile_definitions(RENDER_TEXTURE)
 endif()
-if (USE_EMBREE)
-  add_compile_definitions(USE_EMBREE)
-endif()
 if(USE_PYBIND)
   add_compile_definitions(PYBIND)
 endif()
@@ -97,12 +93,20 @@ if(USE_SUITESPARSE)
 endif()
 
 # Get external libraries
+include(boost)
 include(conformal_ideal_delaunay)
 include(libigl)
 include(spectra)
 include(cli11)
 include(json)
 include(geometry-central)
+#include(symmetric-dirichlet)
+
+include(directional)
+add_library(DirectionalLib INTERFACE)
+target_include_directories(DirectionalLib SYSTEM INTERFACE 
+  ${directional_SOURCE_DIR}/include
+)
 
 # Optionally create rendering library
 if(RENDER_TEXTURE)
@@ -119,12 +123,6 @@ if(RENDER_TEXTURE)
     plot
   )
 
-  # Link in embree if enabled
-  if(USE_EMBREE)
-    target_link_libraries(rendering PUBLIC
-      igl::embree
-    )
-  endif()
 
   # Change rendering libraries from null to singleton
   set(RENDER_LIBRARIES
@@ -144,6 +142,9 @@ endif()
 # Optionally build pybind 
 if(USE_PYBIND)
   include(pybind11)
+  set(PYBIND_LIBRARIES
+    pybind11::module
+  )
 endif()
 
 # Install executables to bin directory

cmake/directional.cmake

@@ -0,0 +1,8 @@
+include(FetchContent)
+FetchContent_Declare(
+    directional
+    SYSTEM
+    GIT_REPOSITORY https://github.com/avaxman/Directional.git
+)
+FetchContent_MakeAvailable(directional)
+

cmake/symmetric-dirichlet.cmake

@@ -0,0 +1,9 @@
+include(FetchContent)
+FetchContent_Declare(
+    symmetric-dirichlet
+    SYSTEM
+    GIT_REPOSITORY https://github.com/rjc8237/symmetric-dirichlet.git
+    GIT_TAG energy_dev
+)
+FetchContent_MakeAvailable(symmetric-dirichlet)
+

include/feature/feature/core/boundary_path.h

@@ -0,0 +1,152 @@
+#pragma once
+
+#include "feature/core/common.h"
+
+namespace Penner {
+namespace Feature {
+
+/**
+ * @brief Representation of a boundary edge path. If the boundary edge is unflipped and lies along
+ * the symmetry line, then this is just a single halfedge in the primal mesh. If the edge is
+ * flipped, then this is a path of halfedges in the mesh copy (type 2) homotopic to the symmetry
+ * line segment between the two vertices on the boundary.
+ *
+ */
+class BoundaryPath
+{
+public:
+    /**
+     * @brief Construct a new Boundary Path object from the given boundary vertex
+     * to the next ccw vertex on the boundary.
+     *
+     * @param m: underlying mesh
+     * @param vertex_index: starting vertex index (must be on the boundary)
+     */
+    BoundaryPath(const Mesh<Scalar>& m, int vertex_index);
+
+    /**
+     * @brief Compute the number of edges in the boundary path
+     * 
+     * @return number of edges
+     */
+    int size() const { return m_halfedge_path.size(); }
+
+    /**
+     * @brief Compute the length of the boundary path in the given mesh
+     *
+     * @param m: underlying mesh
+     * @return length of the current path in the mesh
+     */
+    Scalar compute_length(const Mesh<Scalar>& m) const;
+
+    /**
+     * @brief Compute the (doubled) log length of the boundary path in the given mesh
+     *
+     * @param m: underlying mesh
+     * @return log length of the current path in the mesh
+     */
+    Scalar compute_log_length(const Mesh<Scalar>& m) const;
+
+    /**
+     * @brief Compute the Jacobian of the length of the boundary path in the given mesh with
+     * respect to the halfedge lengths
+     *
+     * The Jacobian is represented as a sparse vector with index, value pairs.
+     *
+     * @param m: underlying mesh
+     * @return jacobian of length
+     */
+    std::vector<std::pair<int, Scalar>> compute_length_jacobian(const Mesh<Scalar>& m) const;
+
+    /**
+     * @brief Compute the Jacobian of the log length of the boundary path in the given mesh with
+     * respect to the halfedge lengths
+     *
+     * The Jacobian is represented as a sparse vector with index, value pairs.
+     *
+     * @param m: underlying mesh
+     * @return jacobian of log length
+     */
+    std::vector<std::pair<int, Scalar>> compute_log_length_jacobian(const Mesh<Scalar>& m) const;
+
+    /**
+     * @brief Check if the boundary path is valid
+     *
+     * @param m: underlying mesh
+     * @return true if the path is valid
+     * @return false otherwise
+     */
+    bool is_valid_boundary_path(const Mesh<Scalar>& m) const;
+
+    /**
+     * @brief Get the start vertex of the path
+     * 
+     * @return start vertex index
+     */
+    int get_start_vertex() const { return m_start_vertex; }
+
+    /**
+     * @brief Get the ending vertex of the path, which is the next original boundary vertex
+     * counterclockwise on the symmetry line.
+     * 
+     * @return end vertex index
+     */
+    int get_end_vertex() const { return m_end_vertex; }
+
+    /**
+     * @brief Get the representative path of halfedges from the start to end vertex;
+     * 
+     * This is either a single primal halfedge or a path of doubled halfedges bounding 
+     * the facets crossing the symmetry line between the vertices.
+     * 
+     * @return path from the start to end vertex
+     */
+    const std::vector<int>& get_halfedge_path() const { return m_halfedge_path; }
+
+    /**
+     * @brief Get transverse halfedges crossing the line of symmetry.
+     * 
+     * One halfedge per transverse edge is returned.
+     * 
+     * @return transverse halfedges.
+     */
+    const std::vector<int>& get_transverse_edges() const { return m_transverse_edges; }
+
+private:
+    int m_start_vertex;
+    int m_end_vertex;
+
+    Scalar compute_tri_length(Scalar side_length, Scalar base_length) const;
+    Scalar compute_quad_length(
+        Scalar side_length,
+        Scalar first_base_length,
+        Scalar second_base_length) const;
+    Scalar compute_tri_side_derivative(Scalar side_length, Scalar base_length) const;
+    Scalar compute_tri_base_derivative(Scalar side_length, Scalar base_length) const;
+    Scalar compute_quad_side_derivative(
+        Scalar side_length,
+        Scalar first_base_length,
+        Scalar second_base_length) const;
+    Scalar compute_quad_base_derivative(
+        Scalar side_length,
+        Scalar first_base_length,
+        Scalar second_base_length) const;
+
+    // Paths between boundary vertices contained in the primal mesh
+    std::vector<int> m_halfedge_path;
+
+    // Transverse edges crossing the symmetry line
+    std::vector<int> m_transverse_edges;
+};
+
+/**
+ * @brief Build all boundary paths for a mesh and a map to the corresponding halfedges.
+ * 
+ * @param m: underlying mesh
+ * @return list of boundary paths
+ * @return map from boundary path indices to halfedge indices
+ */
+std::tuple<std::vector<BoundaryPath>, MatrixX> build_boundary_paths(const Mesh<Scalar>& m);
+
+} // namespace Feature
+} // namespace Penner

include/feature/feature/core/common.h

@@ -0,0 +1,93 @@
+#pragma once
+
+#include "holonomy/interface.h"
+#include "holonomy/core/common.h"
+#include "holonomy/core/dual_loop.h"
+#include "holonomy/holonomy/marked_penner_cone_metric.h"
+#include "holonomy/holonomy/newton.h"
+
+namespace Penner {
+namespace Feature {
+
+using Holonomy::DualLoop;
+using Holonomy::MarkedPennerConeMetric;
+using Holonomy::NewtonParameters;
+using Holonomy::MarkedMetricParameters;
+typedef std::array<int, 2> VertexEdge;
+
+/**
+ * @brief Compute the square of a scalar.
+ *
+ * @param x: value to square
+ * @return squared value
+ */
+inline Scalar square(Scalar x)
+{
+    return x * x;
+}
+
+/**
+ * @brief Compute the one ring of halfedges emenating from a mesh vertex.
+ *
+ * @param m: mesh
+ * @param vertex_index: index of a vertex in the mesh
+ * @return halfedges emenating from a given vertex
+ */
+std::vector<int> generate_vertex_one_ring(const Mesh<Scalar>& m, int vertex_index);
+
+/**
+ * @brief Determine if a VF mesh is manifold
+ * @param F: mesh faces
+ * @return true iff the mesh is edge and vertex manifold
+ */
+bool is_manifold(const Eigen::MatrixXi& F);
+
+/**
+ * @brief Compute salient geometry aligned field directions for a mesh.
+ * 
+ * The parabolic anisotropy used for the relative threshold is ||k2| - |k1|| / max(|k1|, |k2|)
+ * This measurement is near 0 for parabolic regions and near 1 for highly anisotropic regions
+ * 
+ * @param V: mesh vertices
+ * @param F: mesh faces
+ * @param radius: (optional) vertex radius for fitting a smooth surface for field estimation
+ * @param abs_threshold: (optional) minimum threshold for mean anisotropy of principal curvatures
+ * @param rel_threshold: (optional) minimum threshold for parabolic anisotropy of principal curvatures
+ * @return |F|x3 matrix of per face directions
+ * @return per face mask indicating whether a direction is salient or not
+ */
+std::tuple<Eigen::MatrixXd, std::vector<bool>> compute_field_direction(
+    const Eigen::MatrixXd& V,
+    const Eigen::MatrixXi& F,
+    int radius=5,
+    Scalar abs_threshold=1.,
+    Scalar rel_threshold=0.9);
+
+/**
+ * @brief Reindex VF mesh.
+ * 
+ * @param V: mesh vertices
+ * @param F: mesh faces
+ * @param vtx_reindex: map from old to new vertex indices
+ * @return reindexed mesh vertices
+ * @return reindexed mesh faces
+ */
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXi> reindex_mesh(
+    const Eigen::MatrixXd& V,
+    const Eigen::MatrixXi& F,
+    const std::vector<int>& vtx_reindex);
+
+/**
+ * @brief Reindex list of edge endpoints under vertex reindexing.
+ * 
+ * @param endpoints: list of vertex endpoints before reindexing
+ * @param vtx_reindex: map from old to new vertex indices
+ * @return list of reindexed vertex endpoints
+ */
+std::vector<std::pair<int, int>> reindex_endpoints(
+    const std::vector<std::pair<int, int>>& endpoints,
+    const std::vector<int>& vtx_reindex);
+
+
+} // namespace Feature
+} // namespace Penner