Sync with 3384bdc (#47)

.gitignore

@@ -80,3 +80,6 @@ LagrangeOptions.cmake
 # clangd
 compile_commands.json
 keyfile.txt
+
+# ctags
+tags

VERSION

@@ -1 +1 @@
-6.27.0
+6.28.0

cmake/recipes/external/assimp.cmake

@@ -43,6 +43,12 @@ CPMAddPackage(
     NAME assimp
     GITHUB_REPOSITORY assimp/assimp
     GIT_TAG c1deb808faadd85a7a007447b62ae238a4be2337
+
+    PATCHES
+        # Prevent Assimp from meddling with compiler flags in debug mode.
+        # See internal lagrange-lib/#1303 for more details.
+        # Remember to update this patch when updating Assimp.
+        assimp.patch
 )
 
 set_target_properties(assimp PROPERTIES FOLDER third_party/assimp)

cmake/recipes/external/assimp.patch

@@ -0,0 +1,12 @@
+diff --git i/CMakeLists.txt w/CMakeLists.txt
+index 88f69174a..7a3c580f7 100644
+--- i/CMakeLists.txt
++++ w/CMakeLists.txt
+@@ -291,7 +291,6 @@ ELSEIF(MSVC)
+   ENDIF()
+   # supress warning for double to float conversion if Double precission is activated
+   ADD_COMPILE_OPTIONS(/wd4244) 
+-  SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} /D_DEBUG /Zi /Od")
+   SET(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE}")
+   SET(CMAKE_SHARED_LINKER_FLAGS_RELEASE "${CMAKE_SHARED_LINKER_FLAGS_RELEASE} /DEBUG:FULL /PDBALTPATH:%_PDB% /OPT:REF /OPT:ICF")
+ ELSEIF (CMAKE_CXX_COMPILER_ID MATCHES "Clang" )

cmake/recipes/external/nanobind.cmake

@@ -19,7 +19,7 @@ include(CPM)
 CPMAddPackage(
     NAME nanobind
     GITHUB_REPOSITORY wjakob/nanobind
-    GIT_TAG v2.1.0
+    GIT_TAG v2.2.0
     DOWNLOAD_ONLY ON
 )
 

cmake/recipes/external/nasoq.cmake

@@ -48,7 +48,7 @@ include(CPM)
 CPMAddPackage(
     NAME nasoq
     GITHUB_REPOSITORY sympiler/nasoq
-    GIT_TAG fc2051dfa991160cd6dd326d0fb1580ffb77b93b
+    GIT_TAG 03bc29fc10fcc29b68e59fa18b081a1a45779e9b
 )
 
 target_link_libraries(nasoq PUBLIC BLAS::BLAS)

cmake/recipes/external/nlohmann_json.cmake

@@ -16,13 +16,13 @@ endif()
 message(STATUS "Third-party (external): creating target 'nlohmann_json::nlohmann_json'")
 
 # nlohmann_json is a big repo for a single header, so we just download the release archive
-set(NLOHMANNJSON_VERSION "v3.11.2")
+set(NLOHMANNJSON_VERSION "v3.11.3")
 
 include(CPM)
 CPMAddPackage(
     NAME nlohmann_json
     URL "https://github.com/nlohmann/json/releases/download/${NLOHMANNJSON_VERSION}/include.zip"
-    URL_HASH SHA256=e5c7a9f49a16814be27e4ed0ee900ecd0092bfb7dbfca65b5a421b774dccaaed
+    URL_HASH SHA256=a22461d13119ac5c78f205d3df1db13403e58ce1bb1794edc9313677313f4a9d
 )
 
 add_library(nlohmann_json INTERFACE)

modules/core/include/lagrange/SurfaceMesh.h

@@ -551,6 +551,28 @@ class SurfaceMesh
     ///
     void remove_facets(function_ref<bool(Index)> should_remove_func);
 
+    ///
+    /// Reverses the orientation of a list of facets.
+    ///
+    /// @param[in]  facets_to_flip  List of facets to be reversed.
+    ///
+    void flip_facets(span<const Index> facets_to_flip);
+
+    ///
+    /// Reverses the orientation of a list of facets.
+    ///
+    /// @param[in]  facets_to_flip  List of facets to be reversed.
+    ///
+    void flip_facets(std::initializer_list<const Index> facets_to_flip);
+
+    ///
+    /// Reverses the orientation of a list of facets.
+    ///
+    /// @param[in]  should_flip_func  Indicator function returning whether a facet should be
+    ///                               reversed.
+    ///
+    void flip_facets(function_ref<bool(Index)> should_flip_func);
+
     ///
     /// Clear buffer for mesh vertices and other vertex attributes. Since this function also removes
     /// any invalid facet, the entire mesh will be cleared.
@@ -2352,6 +2374,14 @@ class SurfaceMesh
     ///
     void foreach_facet_around_edge(Index e, function_ref<void(Index)> func) const;
 
+    ///
+    /// Applies a function to each of the facets around a prescribed facet.
+    ///
+    /// @param[in]  f     Queried facet index.
+    /// @param[in]  func  Callback to apply to each incident facet.
+    ///
+    void foreach_facet_around_facet(Index f, function_ref<void(Index)> func) const;
+
     ///
     /// Applies a function to each facet around a prescribed vertex.
     ///
@@ -2470,6 +2500,28 @@ class SurfaceMesh
     ///
     std::pair<Index, Index> reindex_facets_internal(span<const Index> old_to_new);
 
+    ///
+    /// Assumptions on the corner mapping received by the internal reindexing method.
+    ///
+    enum class CornerMappingType {
+        /// Buffer is being compressed, so it is safe to move data without an intermediate buffer.
+        RemovingFacets,
+
+        /// Facets are being flipped, so it is safe to simply swap rows i and j when i < j.
+        ReversingFacets,
+    };
+
+    ///
+    /// Reindex mesh corners according to the given mapping. This function is used internally when
+    /// reindexing facets and when reversing facet orientations.
+    ///
+    /// @param[in]  old_to_new_corners  Mapping old corner index -> new corner index.
+    /// @param[in]  mapping_type        Corner mapping type.
+    ///
+    void reindex_corners_internal(
+        function_ref<Index(Index)> old_to_new_corners,
+        CornerMappingType mapping_type);
+
     ///
     /// Resize the buffers associated to a specific element type in the mesh. Newly inserted
     /// elements will be default-initialized.

modules/core/include/lagrange/compute_pointcloud_pca.h

@@ -11,109 +11,68 @@
  */
 #pragma once
 
-// Should we move the definition to a .cpp file? Any .cpp file
-// seeing Eigen/Eigenvalues would take ages to compile.
-#include <Eigen/Core>
-#include <Eigen/Eigenvalues>
+#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS
+    #include <lagrange/legacy/compute_pointcloud_pca.h>
+#endif
 
-#include <lagrange/utils/assert.h>
+#include <lagrange/utils/span.h>
 
 namespace lagrange {
 
-// compute_pointcloud_pca()
-//
-// Finds the principle components for a pointcloud
-// Assumes that the points are supplied in a matrix where each
-// ``row'' is a point.
-//
-// This is closely related to the inertia tensor, principal directions
-// and principal moments. But it is not exactly the same.
-//
-// COVARIANCE_MATRIX (tensor) = (P-eC)^T (P-eC)
-// Where C is the centroid, and e is column vector of ones.
-// eigenvalues and eigenvectors of this matrix would be
-// the principal weights and components.
-//
-// MOMENT OF INERTIA (tensor) = trace(P^T P)I - P^T P
-// eigenvalues and eigenvectors of this matrix would be
-// the principal moments and directions.
-//
-// This refactors segmentation/Pca.h
-
-template <typename Scalar>
-struct ComputePointcloudPCAOutput
+struct ComputePointcloudPCAOptions
 {
-    // The point around which the covariance matrix is evaluated.
-    // Col vector to be consistent with `components` and weights.
-    Eigen::Matrix<Scalar, Eigen::Dynamic, 1> center;
-    // Each column is a component, sorted by weight magnitudes.
-    // n_rows == n_cols == space dimension
-    Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic> components;
-    // Each entry is a weight for the corresponding principal component
-    // n_rows == space dimension
-    // Col vector to be consistent with `components`.
-    Eigen::Matrix<Scalar, Eigen::Dynamic, 1> weights;
+    /**
+     * when true : covariance = (P - centroid) ^ T(P - centroid)
+     * when false : covariance = (P) ^ T(P)
+     */
+    bool shift_centroid = false;
+
+    /**
+     * should we divide the result by number of points?
+     */
+    bool normalize = false;
 };
 
-// Input:
-// points, Each row should be a point in the point cloud.
-//         To be consistent with the rest of Lagrange.
-// should_shift_centeroid,
-//   when true: covariance = (P-centroid)^T (P-centroid)
-//   when false: covariance = (P)^T (P)
-// should_normalize, should we divide the result by number of points?
-template <typename Derived1>
-ComputePointcloudPCAOutput<typename Derived1::Scalar> compute_pointcloud_pca(
-    const Eigen::MatrixBase<Derived1>& points,
-    const bool should_shift_centeroid,
-    const bool should_normalize)
+template <typename Scalar>
+struct PointcloudPCAOutput
 {
-    //
-    using Scalar = typename Derived1::Scalar;
-    using MatrixXS =
-        Eigen::Matrix<Scalar, Derived1::ColsAtCompileTime, Derived1::ColsAtCompileTime>;
-    using RowVectorXS = Eigen::Matrix<Scalar, 1, Derived1::ColsAtCompileTime>;
-
-    // Prechecks
-    la_runtime_assert(points.rows() >= 2, "There must be at least two points");
-
-    // Compute covariance
-    MatrixXS covariance;
-    RowVectorXS center;
-
-    if (should_shift_centeroid) {
-        center = points.colwise().mean();
-    } else {
-        center = RowVectorXS::Zero(points.cols());
-    } // end of should_shift_centroid
-
-    if (should_normalize) {
-        // We may instead divide by points.rows()-1 which applies the Bessel's correction.
-        // https://en.wikipedia.org/wiki/Bessel%27s_correction
-        covariance =
-            ((points.rowwise() - center).transpose() * (points.rowwise() - center)) / points.rows();
-    } else {
-        covariance = ((points.rowwise() - center).transpose() * (points.rowwise() - center));
-    } // end of should_normalize
-
-    // SelfAdjointEigenSolver is guaranteed to return the eigenvalues sorted
-    // in increasing order. Note that this is not true for the general EigenValueSolver.
-    Eigen::SelfAdjointEigenSolver<MatrixXS> eigs(covariance);
-    // Unless something is going wrong, the covariance matrix should
-    // be well behaved. Let's double check.
-    la_runtime_assert(eigs.info() == Eigen::Success, "Eigen decomposition failed");
+    /// The point around which the covariance matrix is evaluated.
+    std::array<Scalar, 3> center;
 
-    ComputePointcloudPCAOutput<Scalar> output;
-    output.center = center.transpose();
-    output.components = eigs.eigenvectors(); // Columns are eigenvectors
-    output.weights = eigs.eigenvalues();
+    /// The 3 components, sorted by weight magnitudes.
+    std::array<std::array<Scalar, 3>, 3> eigenvectors;
 
-    // make sure components follow the right hand rule
-    if (output.components.determinant() < 0.) {
-        output.components.col(0) *= -1;
-    }
+    /// Each entry is a weight for the corresponding principal component
+    std::array<Scalar, 3> eigenvalues;
+};
 
-    return output;
-}
+/**
+ * Finds the principal components for a pointcloud
+ *
+ * Assumes that the points are supplied in a matrix where each
+ * ``row'' is a point.
+ *
+ * This is closely related to the inertia tensor, principal directions
+ * and principal moments. But it is not exactly the same.
+ *
+ * COVARIANCE_MATRIX (tensor) = (P-eC)^T (P-eC)
+ * Where C is the centroid, and e is column vector of ones.
+ * eigenvalues and eigenvectors of this matrix would be
+ * the principal weights and components.
+ *
+ * MOMENT OF INERTIA (tensor) = trace(P^T P)I - P^T P
+ * eigenvalues and eigenvectors of this matrix would be
+ * the principal moments and directions.
+ *
+ * @param[in]  points               The input points.
+ * @param[in]  options              Options struct, see above.
+ *
+ * @return PointCloudPCAOutput      Contains center, eigenvectors (components) and eigenvalues
+ * (weights). See above for more details.
+ */
+template <typename Scalar>
+PointcloudPCAOutput<Scalar> compute_pointcloud_pca(
+    span<const Scalar> points,
+    ComputePointcloudPCAOptions options = {});
 
 } // namespace lagrange

modules/core/include/lagrange/compute_triangle_normal.h

@@ -15,3 +15,5 @@
     #include <lagrange/legacy/compute_triangle_normal.h>
 #endif
 
+// For SurfaceMesh class, please refer to compute_facet_normal.h
+