diff --git a/.gitignore b/.gitignore index 782a0d4..c7dc9a9 100644 --- a/.gitignore +++ b/.gitignore @@ -80,3 +80,6 @@ LagrangeOptions.cmake # clangd compile_commands.json keyfile.txt + +# ctags +tags diff --git a/VERSION b/VERSION index 9cd1a39..2ece8e1 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -6.27.0 +6.28.0 diff --git a/cmake/recipes/external/assimp.cmake b/cmake/recipes/external/assimp.cmake index 6dd35d1..303942b 100644 --- a/cmake/recipes/external/assimp.cmake +++ b/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) diff --git a/cmake/recipes/external/assimp.patch b/cmake/recipes/external/assimp.patch new file mode 100644 index 0000000..649bb75 --- /dev/null +++ b/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" ) diff --git a/cmake/recipes/external/nanobind.cmake b/cmake/recipes/external/nanobind.cmake index ea0f2e8..9f64f53 100644 --- a/cmake/recipes/external/nanobind.cmake +++ b/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 ) diff --git a/cmake/recipes/external/nasoq.cmake b/cmake/recipes/external/nasoq.cmake index 330efef..de020c7 100644 --- a/cmake/recipes/external/nasoq.cmake +++ b/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) diff --git a/cmake/recipes/external/nlohmann_json.cmake b/cmake/recipes/external/nlohmann_json.cmake index 75f5fe5..43969d5 100644 --- a/cmake/recipes/external/nlohmann_json.cmake +++ b/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) diff --git a/modules/core/include/lagrange/SurfaceMesh.h b/modules/core/include/lagrange/SurfaceMesh.h index 11307a6..07af12d 100644 --- a/modules/core/include/lagrange/SurfaceMesh.h +++ b/modules/core/include/lagrange/SurfaceMesh.h @@ -551,6 +551,28 @@ class SurfaceMesh /// void remove_facets(function_ref 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 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 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 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 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 func) const; + /// /// Applies a function to each facet around a prescribed vertex. /// @@ -2470,6 +2500,28 @@ class SurfaceMesh /// std::pair reindex_facets_internal(span 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 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. diff --git a/modules/core/include/lagrange/compute_pointcloud_pca.h b/modules/core/include/lagrange/compute_pointcloud_pca.h index 812f1c5..6e1cc23 100644 --- a/modules/core/include/lagrange/compute_pointcloud_pca.h +++ b/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 -#include +#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS + #include +#endif -#include +#include 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 -struct ComputePointcloudPCAOutput +struct ComputePointcloudPCAOptions { - // The point around which the covariance matrix is evaluated. - // Col vector to be consistent with `components` and weights. - Eigen::Matrix center; - // Each column is a component, sorted by weight magnitudes. - // n_rows == n_cols == space dimension - Eigen::Matrix components; - // Each entry is a weight for the corresponding principal component - // n_rows == space dimension - // Col vector to be consistent with `components`. - Eigen::Matrix 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 -ComputePointcloudPCAOutput compute_pointcloud_pca( - const Eigen::MatrixBase& points, - const bool should_shift_centeroid, - const bool should_normalize) +template +struct PointcloudPCAOutput { - // - using Scalar = typename Derived1::Scalar; - using MatrixXS = - Eigen::Matrix; - using RowVectorXS = Eigen::Matrix; - - // 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 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 center; - ComputePointcloudPCAOutput 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, 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 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 +PointcloudPCAOutput compute_pointcloud_pca( + span points, + ComputePointcloudPCAOptions options = {}); } // namespace lagrange diff --git a/modules/core/include/lagrange/compute_triangle_normal.h b/modules/core/include/lagrange/compute_triangle_normal.h index b0a0428..2f94043 100644 --- a/modules/core/include/lagrange/compute_triangle_normal.h +++ b/modules/core/include/lagrange/compute_triangle_normal.h @@ -15,3 +15,5 @@ #include #endif +// For SurfaceMesh class, please refer to compute_facet_normal.h + diff --git a/modules/core/include/lagrange/legacy/compute_pointcloud_pca.h b/modules/core/include/lagrange/legacy/compute_pointcloud_pca.h new file mode 100644 index 0000000..1d73c43 --- /dev/null +++ b/modules/core/include/lagrange/legacy/compute_pointcloud_pca.h @@ -0,0 +1,123 @@ +/* + * Copyright 2019 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ +#pragma once + +// Should we move the definition to a .cpp file? Any .cpp file +// seeing Eigen/Eigenvalues would take ages to compile. +#include +#include + +#include +#include + +namespace lagrange { +LAGRANGE_LEGACY_INLINE +namespace legacy { + +// 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 +struct ComputePointcloudPCAOutput +{ + // The point around which the covariance matrix is evaluated. + // Col vector to be consistent with `components` and weights. + Eigen::Matrix center; + // Each column is a component, sorted by weight magnitudes. + // n_rows == n_cols == space dimension + Eigen::Matrix components; + // Each entry is a weight for the corresponding principal component + // n_rows == space dimension + // Col vector to be consistent with `components`. + Eigen::Matrix weights; +}; + +// 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 +ComputePointcloudPCAOutput compute_pointcloud_pca( + const Eigen::MatrixBase& points, + const bool should_shift_centeroid, + const bool should_normalize) +{ + // + using Scalar = typename Derived1::Scalar; + using MatrixXS = + Eigen::Matrix; + using RowVectorXS = Eigen::Matrix; + + // 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 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"); + + ComputePointcloudPCAOutput output; + output.center = center.transpose(); + output.components = eigs.eigenvectors(); // Columns are eigenvectors + output.weights = eigs.eigenvalues(); + + // make sure components follow the right hand rule + if (output.components.determinant() < 0.) { + output.components.col(0) *= -1; + } + + return output; +} + +} // namespace legacy +} // namespace lagrange diff --git a/modules/core/include/lagrange/legacy/orient_outward.h b/modules/core/include/lagrange/legacy/orient_outward.h new file mode 100644 index 0000000..0beb4ba --- /dev/null +++ b/modules/core/include/lagrange/legacy/orient_outward.h @@ -0,0 +1,96 @@ +/* + * Copyright 2020 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ +#pragma once + +#include +#include +#include +#include +#include + +#include + +namespace lagrange { +LAGRANGE_LEGACY_INLINE +namespace legacy { + +/// +/// Orient the facets of a mesh so that the signed volume of each connected component is positive or negative. +/// +/// @param[in,out] mesh Mesh whose facets needs to be re-oriented. +/// @param[in] positive Orient to have positive signed volume. +/// +/// @tparam VertexArray Type of vertex array. +/// @tparam FacetArray Type of facet array. +/// +template +void orient_outward(lagrange::Mesh& mesh, bool positive = true) +{ + using Scalar = typename VertexArray::Scalar; + using Index = typename FacetArray::Scalar; + using RowVector3r = Eigen::Matrix; + + if (mesh.get_num_facets() == 0) { + // TODO: Fix igl::bfs_orient to work on empty meshes. + return; + } + + auto tetra_signed_volume = [](const RowVector3r& p1, + const RowVector3r& p2, + const RowVector3r& p3, + const RowVector3r& p4) { + return (p2 - p1).dot((p3 - p1).cross(p4 - p1)) / 6.0; + }; + + auto component_signed_volumes = [&](const auto& vertices, + const auto& facets, + const auto& components, + auto& signed_volumes) { + la_runtime_assert(vertices.cols() == 3); + la_runtime_assert(facets.cols() == 3); + std::array t; + t[3] = RowVector3r::Zero(vertices.cols()); + signed_volumes.resize(components.maxCoeff() + 1); + signed_volumes.setZero(); + for (Eigen::Index f = 0; f < facets.rows(); ++f) { + for (Eigen::Index lv = 0; lv < facets.cols(); ++lv) { + t[lv] = vertices.row(facets(f, lv)); + } + double vol = tetra_signed_volume(t[0], t[1], t[2], t[3]); + signed_volumes(components(f)) -= vol; + } + }; + + const auto& vertices = mesh.get_vertices(); + FacetArray facets; + mesh.export_facets(facets); + + // Orient connected components on the surface + Eigen::Matrix components; + lagrange::internal::bfs_orient(facets, facets, components); + + // Signed volumes per components + Eigen::Matrix signed_volumes; + component_signed_volumes(vertices, facets, components, signed_volumes); + + // Flip facets in each compoenent with incorrect volume sign + for (Eigen::Index f = 0; f < facets.rows(); ++f) { + if ((positive ? 1 : -1) * signed_volumes(components(f)) < 0) { + facets.row(f) = facets.row(f).reverse().eval(); + } + } + + mesh.import_facets(facets); +} + +} // namespace legacy +} // namespace lagrange diff --git a/modules/core/include/lagrange/legacy/select_facets_by_normal_similarity.h b/modules/core/include/lagrange/legacy/select_facets_by_normal_similarity.h new file mode 100644 index 0000000..ba5a3d7 --- /dev/null +++ b/modules/core/include/lagrange/legacy/select_facets_by_normal_similarity.h @@ -0,0 +1,269 @@ +/* + * Copyright 2019 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ + +#pragma once + +#include +#include +#include +#include + +#include + +namespace lagrange { +LAGRANGE_LEGACY_INLINE +namespace legacy { + +/** + * Given a seed vertex, selects faces around it based on the change + * in triangle normals. + */ + +// +// The parameters for the function +// +template +struct SelectFacetsByNormalSimilarityParameters +{ + using Scalar = typename MeshType::Scalar; + using Index = typename MeshType::Index; + + // + // This are input and have to be set + // + // Increasing this would select a larger region + Scalar flood_error_limit = std::numeric_limits::max(); + // This tries to smooth the selection boundary (reduce ears) + bool should_smooth_boundary = true; + + // + // These parameters have default values + // Default value comes from segmentation/Constants.h + // + // If this is not empty, then only faces are considered which are marked as true + std::vector is_facet_selectable = std::vector(); + // Internal param... + Scalar flood_second_to_first_order_limit_ratio = 1.0 / 6.0; + // Number of iterations for smoothing the boundary + Index num_smooth_iterations = 3; + // Use DFS or BFS search + enum class SearchType { BFS = 0, DFS = 1 }; + SearchType search_type = SearchType::DFS; + + // Set selectable facets using a vector of indices + void set_selectable_facets(MeshType& mesh_ref, const std::vector& selectable_facets) + { + if (selectable_facets.empty()) { + is_facet_selectable.clear(); + } else { + is_facet_selectable.resize(mesh_ref.get_num_facets(), false); + for (auto facet_id : range_sparse(mesh_ref.get_num_facets(), selectable_facets)) { + is_facet_selectable[facet_id] = true; + } + } + } + +}; // SelectFaceByNormalSimilarityParameters + +// +// Perform the selection +// mesh, input: the mesh +// seed_facet_id, input: the index of the seed +// parameters, input: the parameters above +// returns a vector that contains true for facets that are selected +// +template +std::vector select_facets_by_normal_similarity( + MeshType& mesh, + const typename MeshType::Index seed_facet_id, + const SelectFacetsByNormalSimilarityParameters& parameters) +{ + // =================================================== + // Helper typedef and lambdas + // =================================================== + + using VertexType = typename MeshType::VertexType; + using Index = typename MeshType::Index; + using Scalar = typename MeshType::Scalar; + using IndexList = typename MeshType::IndexList; + using AttributeArray = typename MeshType::AttributeArray; + using Parameters = SelectFacetsByNormalSimilarityParameters; + + // Check if a face is selectable + auto is_facet_selectable = [&mesh, ¶meters](const Index facet_id) -> bool { + if (parameters.is_facet_selectable.empty()) { + return true; + } else { + la_runtime_assert(facet_id < mesh.get_num_facets()); + return parameters.is_facet_selectable[facet_id]; + } + }; + + // The energy between two normals + auto normal_direction_error = [](const VertexType& n1, const VertexType& n2) -> Scalar { + return static_cast(1.0) - std::abs(n1.dot(n2)); + }; + + // =================================================== + // Let's begin + // current code just copied from + // Segmentation/SuperTriangulation/SingleFloodFromSeed() + // =================================================== + + if (mesh.get_vertex_per_facet() != 3) { + throw std::runtime_error("Input mesh must be a triangle mesh."); + } + + // We need this things + if (!mesh.is_connectivity_initialized()) { + mesh.initialize_connectivity(); + } + if (!mesh.has_facet_attribute("normal")) { + compute_triangle_normal(mesh); + } + + // This would be the final value that will be returned + std::vector is_facet_selected(mesh.get_num_facets(), false); + + // Normals (and the seed normal) + const AttributeArray& facet_normals = mesh.get_facet_attribute("normal"); + const VertexType seed_n = facet_normals.row(seed_facet_id); + + // DFS search utility + std::vector is_facet_processed(mesh.get_num_facets(), false); + std::deque search_queue; + + // + // (0): Add the seed neighbors to the queue to initialize the queue + // + is_facet_processed[seed_facet_id] = true; + is_facet_selected[seed_facet_id] = true; + { + const IndexList& facets_adjacent_to_seed = mesh.get_facets_adjacent_to_facet(seed_facet_id); + for (Index j = 0; j < static_cast(facets_adjacent_to_seed.size()); j++) { + const Index ne_fid = facets_adjacent_to_seed[j]; + const VertexType ne_n = facet_normals.row(ne_fid); + + if ((!is_facet_processed[ne_fid]) && is_facet_selectable(ne_fid)) { + Scalar error = normal_direction_error(seed_n, ne_n); + if (error < parameters.flood_error_limit) { + is_facet_selected[ne_fid] = true; + search_queue.push_back(ne_fid); + } + } + } // End of for over neighbours of seed + } // end of step 0 + + // + // (1): Run through all neighbors, process them, and push them into a queue + // for further flood + // + while (!search_queue.empty()) { + Index fid = invalid(); + if (parameters.search_type == Parameters::SearchType::DFS) { + fid = search_queue.front(); + search_queue.pop_front(); + } else { + fid = search_queue.back(); + search_queue.pop_back(); + } + + const VertexType center_n = facet_normals.row(fid); + + + const IndexList& facets_adjacent_to_candidate = mesh.get_facets_adjacent_to_facet(fid); + for (Index j = 0; j < static_cast(facets_adjacent_to_candidate.size()); j++) { + const Index ne_fid = facets_adjacent_to_candidate[j]; + const VertexType ne_n = facet_normals.row(ne_fid); + + if ((!is_facet_processed[ne_fid]) && is_facet_selectable(ne_fid)) { + const Scalar error_1 = normal_direction_error(seed_n, ne_n); + const Scalar error_2 = normal_direction_error(center_n, ne_n); + is_facet_processed[ne_fid] = true; + + if ((error_1 < parameters.flood_error_limit) && + (error_2 < parameters.flood_error_limit)) { + is_facet_selected[ne_fid] = true; + search_queue.push_back(ne_fid); + } else if ( + error_2 < parameters.flood_error_limit * + parameters.flood_second_to_first_order_limit_ratio) { + // Second order approximation + is_facet_selected[ne_fid] = true; + search_queue.push_back(ne_fid); + } + } // end of is neighbour selectable + } // end of loops over neighbours + } // end of dfs stack + + // + // Smooth the selection + // + if (parameters.should_smooth_boundary) { + for (Index st_iter_2 = 0; st_iter_2 < parameters.num_smooth_iterations; st_iter_2++) { + // Go through boundary faces - check spikes + for (Index i = 0; i < mesh.get_num_facets(); i++) { + const IndexList& facets_adjacent_to_candidate = + mesh.get_facets_adjacent_to_facet(i); + + if (is_facet_selectable(i) && (facets_adjacent_to_candidate.size() > 2)) { + const bool select_flag = is_facet_selected[i]; + Index neighbor_diff_flag = 0; + + const Index ne_fid_0 = facets_adjacent_to_candidate[0]; + const bool ne_select_0 = is_facet_selected[ne_fid_0]; + if (select_flag != ne_select_0) neighbor_diff_flag++; + + const Index ne_fid_1 = facets_adjacent_to_candidate[1]; + const bool ne_select_1 = is_facet_selected[ne_fid_1]; + if (select_flag != ne_select_1) neighbor_diff_flag++; + + const Index ne_fid_2 = facets_adjacent_to_candidate[2]; + const bool ne_select_2 = is_facet_selected[ne_fid_2]; + if (select_flag != ne_select_2) neighbor_diff_flag++; + + if (neighbor_diff_flag > 1) { + const VertexType self_n = facet_normals.row(i); + const VertexType ne_n_0 = facet_normals.row(ne_fid_0); + const VertexType ne_n_1 = facet_normals.row(ne_fid_1); + const VertexType ne_n_2 = facet_normals.row(ne_fid_2); + + const Scalar e_0 = normal_direction_error(self_n, ne_n_0); + const Scalar e_1 = normal_direction_error(self_n, ne_n_1); + const Scalar e_2 = normal_direction_error(self_n, ne_n_2); + + const bool e_0_ok = e_0 < parameters.flood_error_limit; + const bool e_1_ok = e_1 < parameters.flood_error_limit; + const bool e_2_ok = e_2 < parameters.flood_error_limit; + + if (ne_select_0 && ne_select_1 && (e_0_ok || e_1_ok)) { + is_facet_selected[i] = true; + } else if (ne_select_1 && ne_select_2 && (e_1_ok || e_2_ok)) { + is_facet_selected[i] = true; + } else if (ne_select_2 && ne_select_0 && (e_2_ok || e_0_ok)) { + is_facet_selected[i] = true; + } + } + } // end of non selected facet that has two selectable neighbours + } // end of loop over faces + } // end of for(num_smooth_iterations) + } // end of should smooth boundary + + + return is_facet_selected; + +} // select_faces_by_normal_similarity + + + +} // namespace legacy +} // namespace lagrange diff --git a/modules/core/include/lagrange/marching_triangles.h b/modules/core/include/lagrange/marching_triangles.h index 9861351..ebc8cf5 100644 --- a/modules/core/include/lagrange/marching_triangles.h +++ b/modules/core/include/lagrange/marching_triangles.h @@ -1,5 +1,5 @@ /* - * Copyright 2024 Adobe. All rights reserved. + * Copyright 2019 Adobe. All rights reserved. * This file is licensed to you under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. You may obtain a copy * of the License at http://www.apache.org/licenses/LICENSE-2.0 diff --git a/modules/core/include/lagrange/mesh_cleanup/legacy/split_long_edges.h b/modules/core/include/lagrange/mesh_cleanup/legacy/split_long_edges.h new file mode 100644 index 0000000..f154084 --- /dev/null +++ b/modules/core/include/lagrange/mesh_cleanup/legacy/split_long_edges.h @@ -0,0 +1,280 @@ +/* + * Copyright 2016 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ +#pragma once + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace lagrange { +LAGRANGE_LEGACY_INLINE +namespace legacy { + +// TODO: Should accept const ref, needs refactor (non-moving export facet attribute) +template +std::unique_ptr +split_long_edges(MeshType& mesh, typename MeshType::Scalar sq_tol, bool recursive = false) +{ + using Index = typename MeshType::Index; + using Scalar = typename MeshType::Scalar; + using Edge = typename MeshType::Edge; + if (mesh.get_vertex_per_facet() != 3) { + throw "Only triangle is supported"; + } + const Index dim = mesh.get_dim(); + const Index num_vertices = mesh.get_num_vertices(); + const Index num_facets = mesh.get_num_facets(); + const auto& vertices = mesh.get_vertices(); + const auto& facets = mesh.get_facets(); + std::vector additional_vertices; + EdgeMap> splitting_pts; + EdgeSet visited; + std::vector> vertex_mapping; + + typename MeshType::AttributeArray active_facets; + const bool has_active_region = mesh.has_facet_attribute("__is_active"); + if (has_active_region) { + mesh.export_facet_attribute("__is_active", active_facets); + } else { + active_facets.resize(num_facets, 1); + active_facets.setZero(); + } + auto is_active = [&active_facets](Index fid) { return active_facets(fid, 0) != 0; }; + + auto split_edge = [&vertices, + &additional_vertices, + &splitting_pts, + &visited, + sq_tol, + num_vertices, + &vertex_mapping](const Edge& edge) { + if (visited.find(edge) != visited.end()) return; + visited.insert(edge); + if (splitting_pts.find(edge) != splitting_pts.end()) return; + Scalar sq_length = (vertices.row(edge[0]) - vertices.row(edge[1])).squaredNorm(); + if (sq_length <= sq_tol) return; + + const Scalar num_segments = std::ceil(sqrt(sq_length / sq_tol)); + const auto v0 = vertices.row(edge[0]).eval(); + const auto v1 = vertices.row(edge[1]).eval(); + std::vector split_pts_indices; + split_pts_indices.push_back(edge[0]); + const Index base = safe_cast(additional_vertices.size()) + num_vertices; + for (Index i = 1; i < num_segments; i++) { + Scalar ratio = i / num_segments; + additional_vertices.push_back(v0 * (1.0f - ratio) + v1 * ratio); + vertex_mapping.emplace_back(edge[0], edge[1], 1.0f - ratio); + split_pts_indices.push_back(base + i - 1); + } + split_pts_indices.push_back(edge[1]); + splitting_pts.insert({edge, split_pts_indices}); + }; + + for (Index i = 0; i < num_facets; i++) { + if (has_active_region && !is_active(i)) continue; + split_edge({{facets(i, 0), facets(i, 1)}}); + split_edge({{facets(i, 1), facets(i, 2)}}); + split_edge({{facets(i, 2), facets(i, 0)}}); + } + + // Concatenate original vertices and additional vertices. + const Index total_num_vertices = num_vertices + safe_cast(additional_vertices.size()); + la_runtime_assert(vertex_mapping.size() == additional_vertices.size()); + typename MeshType::VertexArray all_vertices(total_num_vertices, dim); + all_vertices.block(0, 0, num_vertices, dim) = vertices; + for (Index i = num_vertices; i < total_num_vertices; i++) { + all_vertices.row(i) = additional_vertices[i - num_vertices]; + } + + // Re-triangulate facets. + std::vector> out_facets; + std::vector facet_map; + for (Index i = 0; i < num_facets; i++) { + Index corners[3]; + std::vector chain; + + // Copy over inactive facets. + if (has_active_region && !is_active(i)) { + la_runtime_assert(facets.row(i).maxCoeff() < total_num_vertices); + la_runtime_assert(facets.row(i).minCoeff() >= 0); + out_facets.push_back(facets.row(i)); + facet_map.push_back(i); + continue; + } + + for (Index j = 0; j < 3; j++) { + Edge e{{facets(i, j), facets(i, (j + 1) % 3)}}; + corners[j] = safe_cast(chain.size()); + chain.push_back(e[0]); + auto itr = splitting_pts.find(e); + if (itr != splitting_pts.end()) { + const auto& pts = itr->second; + la_runtime_assert(pts.size() >= 3); + if (pts[0] == e[0]) { + auto start = std::next(pts.cbegin()); + auto end = std::prev(pts.cend()); + chain.insert(chain.end(), start, end); + } else { + la_runtime_assert(pts[0] == e[1]); + auto start = std::next(pts.crbegin()); + auto end = std::prev(pts.crend()); + chain.insert(chain.end(), start, end); + } + } + } + if (chain.size() == 3) { + // No split. + la_runtime_assert(facets.row(i).maxCoeff() < total_num_vertices); + la_runtime_assert(facets.row(i).minCoeff() >= 0); + out_facets.push_back(facets.row(i)); + facet_map.push_back(i); + active_facets(i, 0) = 0; + } else { + auto additional_facets = + split_triangle(all_vertices, chain, corners[0], corners[1], corners[2]); + std::copy( + additional_facets.begin(), + additional_facets.end(), + std::back_inserter(out_facets)); + facet_map.insert(facet_map.end(), additional_facets.size(), i); + active_facets(i, 0) = 1; + } + } + + const Index num_out_facets = safe_cast(out_facets.size()); + typename MeshType::FacetArray all_facets(num_out_facets, 3); + for (Index i = 0; i < num_out_facets; i++) { + la_runtime_assert(out_facets[i].minCoeff() >= 0); + la_runtime_assert(out_facets[i].maxCoeff() < total_num_vertices); + all_facets.row(i) = out_facets[i]; + } + + // Mark active facets (i.e. facets that are split). + if (has_active_region) { + mesh.import_facet_attribute("__is_active", active_facets); + } else { + mesh.add_facet_attribute("__is_active"); + mesh.import_facet_attribute("__is_active", active_facets); + } + + auto out_mesh = create_mesh(all_vertices, all_facets); + + // Port vertex attributes + auto vertex_attributes = mesh.get_vertex_attribute_names(); + auto map_vertex_fn = [&](Eigen::Index i, + std::vector>& weights) { + weights.clear(); + if (i < num_vertices) { + weights.emplace_back(i, 1.0); + } else { + const auto& record = vertex_mapping[i - num_vertices]; + Index v0 = std::get<0>(record); + Index v1 = std::get<1>(record); + Scalar ratio = std::get<2>(record); + weights.emplace_back(v0, ratio); + weights.emplace_back(v1, 1.0 - ratio); + } + }; + + for (const auto& name : vertex_attributes) { + auto attr = mesh.get_vertex_attribute_array(name); + auto attr2 = to_shared_ptr(attr->row_slice(total_num_vertices, map_vertex_fn)); + out_mesh->add_vertex_attribute(name); + out_mesh->set_vertex_attribute_array(name, std::move(attr2)); + } + + // Port facet attributes + auto facet_attributes = mesh.get_facet_attribute_names(); + for (const auto& name : facet_attributes) { + auto attr = mesh.get_facet_attribute_array(name); + auto attr2 = to_shared_ptr(attr->row_slice(facet_map)); + out_mesh->add_facet_attribute(name); + out_mesh->set_facet_attribute_array(name, std::move(attr2)); + } + + // Port corner attributes + map_corner_attributes(mesh, *out_mesh, facet_map); + + // Port indexed attributes. + const auto& indexed_attr_names = mesh.get_indexed_attribute_names(); + typename MeshType::AttributeArray attr; + typename MeshType::IndexArray indices; + for (const auto& attr_name : indexed_attr_names) { + std::tie(attr, indices) = mesh.get_indexed_attribute(attr_name); + assert(indices.rows() == facets.rows()); + + auto get_vertex_index_in_facet = [&facets](Index fid, Index vid) -> Index { + if (vid == facets(fid, 0)) + return 0; + else if (vid == facets(fid, 1)) + return 1; + else { + assert(vid == facets(fid, 2)); + return 2; + } + }; + + typename MeshType::AttributeArray out_attr(num_out_facets * 3, attr.cols()); + for (auto i : range(num_out_facets)) { + const auto old_fid = facet_map[i]; + assert(old_fid < facets.rows()); + + for (Index j = 0; j < 3; j++) { + const auto vid = all_facets(i, j); + if (vid < num_vertices) { + const auto old_j = get_vertex_index_in_facet(old_fid, vid); + // Vertex is part of the input vertices. + out_attr.row(i * 3 + j) = attr.row(indices(old_fid, old_j)); + } else { + const auto& record = vertex_mapping[vid - num_vertices]; + Index v0 = std::get<0>(record); + Index v1 = std::get<1>(record); + Scalar ratio = std::get<2>(record); + Index old_j0 = get_vertex_index_in_facet(old_fid, v0); + Index old_j1 = get_vertex_index_in_facet(old_fid, v1); + + out_attr.row(i * 3 + j) = attr.row(indices(old_fid, old_j0)) * ratio + + attr.row(indices(old_fid, old_j1)) * (1.0f - ratio); + } + } + } + + if (attr_name == "normal") { + for (auto i : range(num_out_facets * 3)) { + out_attr.row(i).stableNormalize(); + } + } + + const std::string tmp_name = "__" + attr_name; + out_mesh->add_corner_attribute(tmp_name); + out_mesh->import_corner_attribute(tmp_name, out_attr); + map_corner_attribute_to_indexed_attribute(*out_mesh, tmp_name); + rename_indexed_attribute(*out_mesh, tmp_name, attr_name); + out_mesh->remove_corner_attribute(tmp_name); + } + + if (recursive && total_num_vertices > num_vertices) { + return split_long_edges(*out_mesh, sq_tol, recursive); + } else + return out_mesh; +} +} // namespace legacy +} // namespace lagrange diff --git a/modules/core/include/lagrange/mesh_cleanup/legacy/split_triangle.h b/modules/core/include/lagrange/mesh_cleanup/legacy/split_triangle.h new file mode 100644 index 0000000..d6dfdfa --- /dev/null +++ b/modules/core/include/lagrange/mesh_cleanup/legacy/split_triangle.h @@ -0,0 +1,172 @@ +/* + * Copyright 2016 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include + +namespace lagrange { +LAGRANGE_LEGACY_INLINE +namespace legacy { +/** + * Split triangle into smaller triangles based on the chain of spliting + * points. + * + * Arguments: + * vertices: list of vertices that contains all 3 corner of the triangle + * and all splitting points. + * chain: A chain of indices into vertices that iterates all + * splitting points and corners in counterclockwise order. + * v0: index into chain that represents corner 0. + * v1: index into chain that represents corner 1. + * v2: index into chain that represents corner 2. + * + * Returns: + * facets: output facets. + */ +template +std::vector> split_triangle( + const VertexArray& vertices, + const std::vector& chain, + const Index v0, + const Index v1, + const Index v2) +{ + const Index chain_size = safe_cast(chain.size()); + auto next = [&chain_size](Index i) { return (i + 1) % chain_size; }; + auto prev = [&chain_size](Index i) { return (i + chain_size - 1) % chain_size; }; + auto sq_length = [&vertices, &chain](Index vi, Index vj) { + return (vertices.row(chain[vi]) - vertices.row(chain[vj])).squaredNorm(); + }; + auto is_corner = [v0, v1, v2](Index v) { return v == v0 || v == v1 || v == v2; }; + std::vector> facets; + + const double NOT_USED = std::numeric_limits::max(); + Eigen::Matrix candidates; + candidates << v0, next(v0), prev(v0), 0, 0, 0, v1, next(v1), prev(v1), 0, 0, 0, v2, next(v2), + prev(v2), 0, 0, 0; + + Eigen::Matrix candidate_lengths; + candidate_lengths << sq_length(candidates(0, 1), candidates(0, 2)), NOT_USED, + sq_length(candidates(1, 1), candidates(1, 2)), NOT_USED, + sq_length(candidates(2, 1), candidates(2, 2)), NOT_USED; + + auto index_comp = [&candidate_lengths](Index i, Index j) { + // Use greater than operator so the heap is min heap. + return candidate_lengths.row(i).minCoeff() > candidate_lengths.row(j).minCoeff(); + }; + std::priority_queue, decltype(index_comp)> Q(index_comp); + Q.push(0); + Q.push(1); + Q.push(2); + + Eigen::Matrix visited(chain_size, 3); + visited.setZero(); + visited(v0, 0) = 1; + visited(v1, 1) = 1; + visited(v2, 2) = 1; + + while (!Q.empty()) { + Index idx = Q.top(); + Q.pop(); + Index selection = 0; + if (candidate_lengths(idx, 0) != NOT_USED && + candidate_lengths(idx, 0) <= candidate_lengths(idx, 1)) { + selection = 0; + } else if ( + candidate_lengths(idx, 1) != NOT_USED && + candidate_lengths(idx, 1) <= candidate_lengths(idx, 0)) { + selection = 1; + } else { + // Select neither candidate from this corner. + continue; + } + la_runtime_assert(candidate_lengths.row(idx).minCoeff() <= candidate_lengths.minCoeff()); + + Index base_v = candidates(idx, selection * 3); + Index right_v = candidates(idx, selection * 3 + 1); + Index left_v = candidates(idx, selection * 3 + 2); + la_runtime_assert(base_v >= 0 && base_v < chain_size); + la_runtime_assert(right_v >= 0 && right_v < chain_size); + la_runtime_assert(left_v >= 0 && left_v < chain_size); + la_runtime_assert(visited(base_v, idx) >= 1); + + // A special case. + if (is_corner(base_v) && is_corner(right_v) && is_corner(left_v)) { + if (chain_size != 3) { + candidate_lengths(idx, selection) = NOT_USED; + Q.push(idx); + continue; + } + } + + if (visited.row(base_v).sum() > 1 || visited(right_v, idx) > 1 || + visited(left_v, idx) > 1) { + // This canadidate is invalid. + candidate_lengths(idx, selection) = NOT_USED; + Q.push(idx); + continue; + } + + visited(base_v, idx) = 2; + visited(right_v, idx) = 1; + visited(left_v, idx) = 1; + facets.push_back({chain[base_v], chain[right_v], chain[left_v]}); + + // Update candidate from this corner. + if (visited.row(right_v).sum() == 1) { + Index right_to_right = next(right_v); + candidate_lengths(idx, 0) = sq_length(left_v, right_to_right); + candidates(idx, 0) = right_v; + candidates(idx, 1) = right_to_right; + candidates(idx, 2) = left_v; + } else { + candidate_lengths(idx, 0) = NOT_USED; + } + + if (visited.row(left_v).sum() == 1) { + Index left_to_left = prev(left_v); + candidate_lengths(idx, 1) = sq_length(right_v, left_to_left); + candidates(idx, 3) = left_v; + candidates(idx, 4) = right_v; + candidates(idx, 5) = left_to_left; + } else { + candidate_lengths(idx, 1) = NOT_USED; + } + Q.push(idx); + } + + Eigen::Matrix visited_sum = + (visited.array() > 0).rowwise().count(); + if ((visited_sum.array() > 1).count() == 3) { + Index f[3]; + Index count = 0; + for (Index i = 0; i < chain_size; i++) { + if (visited_sum[i] > 1) { + la_runtime_assert(count < 3); + f[count] = chain[i]; + count++; + } + } + la_runtime_assert(count == 3); + facets.push_back({f[0], f[1], f[2]}); + } + return facets; +} +} // namespace legacy +} // namespace lagrange diff --git a/modules/core/include/lagrange/mesh_cleanup/resolve_nonmanifoldness.h b/modules/core/include/lagrange/mesh_cleanup/resolve_nonmanifoldness.h index 8913126..03813d6 100644 --- a/modules/core/include/lagrange/mesh_cleanup/resolve_nonmanifoldness.h +++ b/modules/core/include/lagrange/mesh_cleanup/resolve_nonmanifoldness.h @@ -1,5 +1,5 @@ /* - * Copyright 2024 Adobe. All rights reserved. + * Copyright 2018 Adobe. All rights reserved. * This file is licensed to you under the Apache License, Version 2.0 (the "License"); * you may not use this file except in compliance with the License. You may obtain a copy * of the License at http://www.apache.org/licenses/LICENSE-2.0 diff --git a/modules/core/include/lagrange/mesh_cleanup/split_long_edges.h b/modules/core/include/lagrange/mesh_cleanup/split_long_edges.h index 0ccd484..a393073 100644 --- a/modules/core/include/lagrange/mesh_cleanup/split_long_edges.h +++ b/modules/core/include/lagrange/mesh_cleanup/split_long_edges.h @@ -11,266 +11,44 @@ */ #pragma once -#include +#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS + #include +#endif -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include -namespace lagrange { - -// TODO: Should accept const ref, needs refactor (non-moving export facet attribute) -template -std::unique_ptr -split_long_edges(MeshType& mesh, typename MeshType::Scalar sq_tol, bool recursive = false) -{ - using Index = typename MeshType::Index; - using Scalar = typename MeshType::Scalar; - using Edge = typename MeshType::Edge; - if (mesh.get_vertex_per_facet() != 3) { - throw "Only triangle is supported"; - } - const Index dim = mesh.get_dim(); - const Index num_vertices = mesh.get_num_vertices(); - const Index num_facets = mesh.get_num_facets(); - const auto& vertices = mesh.get_vertices(); - const auto& facets = mesh.get_facets(); - std::vector additional_vertices; - EdgeMap> splitting_pts; - EdgeSet visited; - std::vector> vertex_mapping; - - typename MeshType::AttributeArray active_facets; - const bool has_active_region = mesh.has_facet_attribute("__is_active"); - if (has_active_region) { - mesh.export_facet_attribute("__is_active", active_facets); - } else { - active_facets.resize(num_facets, 1); - active_facets.setZero(); - } - auto is_active = [&active_facets](Index fid) { return active_facets(fid, 0) != 0; }; - - auto split_edge = [&vertices, - &additional_vertices, - &splitting_pts, - &visited, - sq_tol, - num_vertices, - &vertex_mapping](const Edge& edge) { - if (visited.find(edge) != visited.end()) return; - visited.insert(edge); - if (splitting_pts.find(edge) != splitting_pts.end()) return; - Scalar sq_length = (vertices.row(edge[0]) - vertices.row(edge[1])).squaredNorm(); - if (sq_length <= sq_tol) return; - - const Scalar num_segments = std::ceil(sqrt(sq_length / sq_tol)); - const auto v0 = vertices.row(edge[0]).eval(); - const auto v1 = vertices.row(edge[1]).eval(); - std::vector split_pts_indices; - split_pts_indices.push_back(edge[0]); - const Index base = safe_cast(additional_vertices.size()) + num_vertices; - for (Index i = 1; i < num_segments; i++) { - Scalar ratio = i / num_segments; - additional_vertices.push_back(v0 * (1.0f - ratio) + v1 * ratio); - vertex_mapping.emplace_back(edge[0], edge[1], 1.0f - ratio); - split_pts_indices.push_back(base + i - 1); - } - split_pts_indices.push_back(edge[1]); - splitting_pts.insert({edge, split_pts_indices}); - }; - - for (Index i = 0; i < num_facets; i++) { - if (has_active_region && !is_active(i)) continue; - split_edge({{facets(i, 0), facets(i, 1)}}); - split_edge({{facets(i, 1), facets(i, 2)}}); - split_edge({{facets(i, 2), facets(i, 0)}}); - } - - // Concatenate original vertices and additional vertices. - const Index total_num_vertices = num_vertices + safe_cast(additional_vertices.size()); - la_runtime_assert(vertex_mapping.size() == additional_vertices.size()); - typename MeshType::VertexArray all_vertices(total_num_vertices, dim); - all_vertices.block(0, 0, num_vertices, dim) = vertices; - for (Index i = num_vertices; i < total_num_vertices; i++) { - all_vertices.row(i) = additional_vertices[i - num_vertices]; - } - - // Re-triangulate facets. - std::vector> out_facets; - std::vector facet_map; - for (Index i = 0; i < num_facets; i++) { - Index corners[3]; - std::vector chain; - - // Copy over inactive facets. - if (has_active_region && !is_active(i)) { - la_runtime_assert(facets.row(i).maxCoeff() < total_num_vertices); - la_runtime_assert(facets.row(i).minCoeff() >= 0); - out_facets.push_back(facets.row(i)); - facet_map.push_back(i); - continue; - } - - for (Index j = 0; j < 3; j++) { - Edge e{{facets(i, j), facets(i, (j + 1) % 3)}}; - corners[j] = safe_cast(chain.size()); - chain.push_back(e[0]); - auto itr = splitting_pts.find(e); - if (itr != splitting_pts.end()) { - const auto& pts = itr->second; - la_runtime_assert(pts.size() >= 3); - if (pts[0] == e[0]) { - auto start = std::next(pts.cbegin()); - auto end = std::prev(pts.cend()); - chain.insert(chain.end(), start, end); - } else { - la_runtime_assert(pts[0] == e[1]); - auto start = std::next(pts.crbegin()); - auto end = std::prev(pts.crend()); - chain.insert(chain.end(), start, end); - } - } - } - if (chain.size() == 3) { - // No split. - la_runtime_assert(facets.row(i).maxCoeff() < total_num_vertices); - la_runtime_assert(facets.row(i).minCoeff() >= 0); - out_facets.push_back(facets.row(i)); - facet_map.push_back(i); - active_facets(i, 0) = 0; - } else { - auto additional_facets = - split_triangle(all_vertices, chain, corners[0], corners[1], corners[2]); - std::copy( - additional_facets.begin(), - additional_facets.end(), - std::back_inserter(out_facets)); - facet_map.insert(facet_map.end(), additional_facets.size(), i); - active_facets(i, 0) = 1; - } - } - - const Index num_out_facets = safe_cast(out_facets.size()); - typename MeshType::FacetArray all_facets(num_out_facets, 3); - for (Index i = 0; i < num_out_facets; i++) { - la_runtime_assert(out_facets[i].minCoeff() >= 0); - la_runtime_assert(out_facets[i].maxCoeff() < total_num_vertices); - all_facets.row(i) = out_facets[i]; - } +#include - // Mark active facets (i.e. facets that are split). - if (has_active_region) { - mesh.import_facet_attribute("__is_active", active_facets); - } else { - mesh.add_facet_attribute("__is_active"); - mesh.import_facet_attribute("__is_active", active_facets); - } - auto out_mesh = create_mesh(all_vertices, all_facets); - - // Port vertex attributes - auto vertex_attributes = mesh.get_vertex_attribute_names(); - auto map_vertex_fn = [&](Eigen::Index i, - std::vector>& weights) { - weights.clear(); - if (i < num_vertices) { - weights.emplace_back(i, 1.0); - } else { - const auto& record = vertex_mapping[i - num_vertices]; - Index v0 = std::get<0>(record); - Index v1 = std::get<1>(record); - Scalar ratio = std::get<2>(record); - weights.emplace_back(v0, ratio); - weights.emplace_back(v1, 1.0 - ratio); - } - }; - - for (const auto& name : vertex_attributes) { - auto attr = mesh.get_vertex_attribute_array(name); - auto attr2 = to_shared_ptr(attr->row_slice(total_num_vertices, map_vertex_fn)); - out_mesh->add_vertex_attribute(name); - out_mesh->set_vertex_attribute_array(name, std::move(attr2)); - } - - // Port facet attributes - auto facet_attributes = mesh.get_facet_attribute_names(); - for (const auto& name : facet_attributes) { - auto attr = mesh.get_facet_attribute_array(name); - auto attr2 = to_shared_ptr(attr->row_slice(facet_map)); - out_mesh->add_facet_attribute(name); - out_mesh->set_facet_attribute_array(name, std::move(attr2)); - } - - // Port corner attributes - map_corner_attributes(mesh, *out_mesh, facet_map); - - // Port indexed attributes. - const auto& indexed_attr_names = mesh.get_indexed_attribute_names(); - typename MeshType::AttributeArray attr; - typename MeshType::IndexArray indices; - for (const auto& attr_name : indexed_attr_names) { - std::tie(attr, indices) = mesh.get_indexed_attribute(attr_name); - assert(indices.rows() == facets.rows()); - - auto get_vertex_index_in_facet = [&facets](Index fid, Index vid) -> Index { - if (vid == facets(fid, 0)) - return 0; - else if (vid == facets(fid, 1)) - return 1; - else { - assert(vid == facets(fid, 2)); - return 2; - } - }; - - typename MeshType::AttributeArray out_attr(num_out_facets * 3, attr.cols()); - for (auto i : range(num_out_facets)) { - const auto old_fid = facet_map[i]; - assert(old_fid < facets.rows()); - - for (Index j = 0; j < 3; j++) { - const auto vid = all_facets(i, j); - if (vid < num_vertices) { - const auto old_j = get_vertex_index_in_facet(old_fid, vid); - // Vertex is part of the input vertices. - out_attr.row(i * 3 + j) = attr.row(indices(old_fid, old_j)); - } else { - const auto& record = vertex_mapping[vid - num_vertices]; - Index v0 = std::get<0>(record); - Index v1 = std::get<1>(record); - Scalar ratio = std::get<2>(record); - Index old_j0 = get_vertex_index_in_facet(old_fid, v0); - Index old_j1 = get_vertex_index_in_facet(old_fid, v1); - - out_attr.row(i * 3 + j) = attr.row(indices(old_fid, old_j0)) * ratio + - attr.row(indices(old_fid, old_j1)) * (1.0f - ratio); - } - } - } - - if (attr_name == "normal") { - for (auto i : range(num_out_facets * 3)) { - out_attr.row(i).stableNormalize(); - } - } +namespace lagrange { - const std::string tmp_name = "__" + attr_name; - out_mesh->add_corner_attribute(tmp_name); - out_mesh->import_corner_attribute(tmp_name, out_attr); - map_corner_attribute_to_indexed_attribute(*out_mesh, tmp_name); - rename_indexed_attribute(*out_mesh, tmp_name, attr_name); - out_mesh->remove_corner_attribute(tmp_name); - } +struct SplitLongEdgesOptions +{ + /// Maximum edge length. Edges longer than this value will be split. + float max_edge_length = 0.1f; + + /// If true, the operation will be applied recursively until no edge is longer than + /// `max_edge_length`. + bool recursive = true; + + /// The facet attribute name that will be used to determine the active region. + /// Only edges that are part of the active region will be split. + /// If empty, all edges will be considered. + /// The attribute must be of value type `uint8_t`. + std::string_view active_region_attribute = ""; + + /// The edge length attribute name that will be used to determine the edge length. + /// If this attribute does not exist, the function will compute the edge lengths. + std::string_view edge_length_attribute = "@edge_length"; +}; + +/// +/// Split edges that are longer than `options.max_edge_length`. +/// +/// @param[in,out] mesh Input mesh. +/// @param[in] options Optional settings. +/// +template +void split_long_edges(SurfaceMesh& mesh, SplitLongEdgesOptions options = {}); - if (recursive && total_num_vertices > num_vertices) { - return split_long_edges(*out_mesh, sq_tol, recursive); - } else - return out_mesh; -} } // namespace lagrange diff --git a/modules/core/include/lagrange/mesh_cleanup/split_triangle.h b/modules/core/include/lagrange/mesh_cleanup/split_triangle.h index 1c5b46d..b7e41cc 100644 --- a/modules/core/include/lagrange/mesh_cleanup/split_triangle.h +++ b/modules/core/include/lagrange/mesh_cleanup/split_triangle.h @@ -11,158 +11,7 @@ */ #pragma once -#include -#include -#include -#include +#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS + #include +#endif -#include -#include - -namespace lagrange { -/** - * Split triangle into smaller triangles based on the chain of spliting - * points. - * - * Arguments: - * vertices: list of vertices that contains all 3 corner of the triangle - * and all splitting points. - * chain: A chain of indices into vertices that iterates all - * splitting points and corners in counterclockwise order. - * v0: index into chain that represents corner 0. - * v1: index into chain that represents corner 1. - * v2: index into chain that represents corner 2. - * - * Returns: - * facets: output facets. - */ -template -std::vector> split_triangle( - const VertexArray& vertices, - const std::vector& chain, - const Index v0, - const Index v1, - const Index v2) -{ - const Index chain_size = safe_cast(chain.size()); - auto next = [&chain_size](Index i) { return (i + 1) % chain_size; }; - auto prev = [&chain_size](Index i) { return (i + chain_size - 1) % chain_size; }; - auto sq_length = [&vertices, &chain](Index vi, Index vj) { - return (vertices.row(chain[vi]) - vertices.row(chain[vj])).squaredNorm(); - }; - auto is_corner = [v0, v1, v2](Index v) { return v == v0 || v == v1 || v == v2; }; - std::vector> facets; - - const double NOT_USED = std::numeric_limits::max(); - Eigen::Matrix candidates; - candidates << v0, next(v0), prev(v0), 0, 0, 0, v1, next(v1), prev(v1), 0, 0, 0, v2, next(v2), - prev(v2), 0, 0, 0; - - Eigen::Matrix candidate_lengths; - candidate_lengths << sq_length(candidates(0, 1), candidates(0, 2)), NOT_USED, - sq_length(candidates(1, 1), candidates(1, 2)), NOT_USED, - sq_length(candidates(2, 1), candidates(2, 2)), NOT_USED; - - auto index_comp = [&candidate_lengths](Index i, Index j) { - // Use greater than operator so the heap is min heap. - return candidate_lengths.row(i).minCoeff() > candidate_lengths.row(j).minCoeff(); - }; - std::priority_queue, decltype(index_comp)> Q(index_comp); - Q.push(0); - Q.push(1); - Q.push(2); - - Eigen::Matrix visited(chain_size, 3); - visited.setZero(); - visited(v0, 0) = 1; - visited(v1, 1) = 1; - visited(v2, 2) = 1; - - while (!Q.empty()) { - Index idx = Q.top(); - Q.pop(); - Index selection = 0; - if (candidate_lengths(idx, 0) != NOT_USED && - candidate_lengths(idx, 0) <= candidate_lengths(idx, 1)) { - selection = 0; - } else if ( - candidate_lengths(idx, 1) != NOT_USED && - candidate_lengths(idx, 1) <= candidate_lengths(idx, 0)) { - selection = 1; - } else { - // Select neither candidate from this corner. - continue; - } - la_runtime_assert(candidate_lengths.row(idx).minCoeff() <= candidate_lengths.minCoeff()); - - Index base_v = candidates(idx, selection * 3); - Index right_v = candidates(idx, selection * 3 + 1); - Index left_v = candidates(idx, selection * 3 + 2); - la_runtime_assert(base_v >= 0 && base_v < chain_size); - la_runtime_assert(right_v >= 0 && right_v < chain_size); - la_runtime_assert(left_v >= 0 && left_v < chain_size); - la_runtime_assert(visited(base_v, idx) >= 1); - - // A special case. - if (is_corner(base_v) && is_corner(right_v) && is_corner(left_v)) { - if (chain_size != 3) { - candidate_lengths(idx, selection) = NOT_USED; - Q.push(idx); - continue; - } - } - - if (visited.row(base_v).sum() > 1 || visited(right_v, idx) > 1 || - visited(left_v, idx) > 1) { - // This canadidate is invalid. - candidate_lengths(idx, selection) = NOT_USED; - Q.push(idx); - continue; - } - - visited(base_v, idx) = 2; - visited(right_v, idx) = 1; - visited(left_v, idx) = 1; - facets.push_back({chain[base_v], chain[right_v], chain[left_v]}); - - // Update candidate from this corner. - if (visited.row(right_v).sum() == 1) { - Index right_to_right = next(right_v); - candidate_lengths(idx, 0) = sq_length(left_v, right_to_right); - candidates(idx, 0) = right_v; - candidates(idx, 1) = right_to_right; - candidates(idx, 2) = left_v; - } else { - candidate_lengths(idx, 0) = NOT_USED; - } - - if (visited.row(left_v).sum() == 1) { - Index left_to_left = prev(left_v); - candidate_lengths(idx, 1) = sq_length(right_v, left_to_left); - candidates(idx, 3) = left_v; - candidates(idx, 4) = right_v; - candidates(idx, 5) = left_to_left; - } else { - candidate_lengths(idx, 1) = NOT_USED; - } - Q.push(idx); - } - - Eigen::Matrix visited_sum = - (visited.array() > 0).rowwise().count(); - if ((visited_sum.array() > 1).count() == 3) { - Index f[3]; - Index count = 0; - for (Index i = 0; i < chain_size; i++) { - if (visited_sum[i] > 1) { - la_runtime_assert(count < 3); - f[count] = chain[i]; - count++; - } - } - la_runtime_assert(count == 3); - facets.push_back({f[0], f[1], f[2]}); - } - return facets; -} -} // namespace lagrange diff --git a/modules/core/include/lagrange/orient_outward.h b/modules/core/include/lagrange/orient_outward.h index 85f85e6..5f2fce9 100644 --- a/modules/core/include/lagrange/orient_outward.h +++ b/modules/core/include/lagrange/orient_outward.h @@ -11,82 +11,34 @@ */ #pragma once -#include -#include -#include -#include +#include -#include +#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS + #include +#endif namespace lagrange { /// -/// Orient the facets of a mesh so that the signed volume of each connected component is positive or negative. +/// Options for orienting the facets of a mesh. /// -/// @param[in,out] mesh Mesh whose facets needs to be re-oriented. -/// @param[in] positive Orient to have positive signed volume. -/// -/// @tparam VertexArray Type of vertex array. -/// @tparam FacetArray Type of facet array. -/// -template -void orient_outward(lagrange::Mesh& mesh, bool positive = true) +struct OrientOptions { - using Scalar = typename VertexArray::Scalar; - using Index = typename FacetArray::Scalar; - using RowVector3r = Eigen::Matrix; - - if (mesh.get_num_facets() == 0) { - // TODO: Fix igl::bfs_orient to work on empty meshes. - return; - } - - auto tetra_signed_volume = [](const RowVector3r& p1, - const RowVector3r& p2, - const RowVector3r& p3, - const RowVector3r& p4) { - return (p2 - p1).dot((p3 - p1).cross(p4 - p1)) / 6.0; - }; - - auto component_signed_volumes = [&](const auto& vertices, - const auto& facets, - const auto& components, - auto& signed_volumes) { - la_runtime_assert(vertices.cols() == 3); - la_runtime_assert(facets.cols() == 3); - std::array t; - t[3] = RowVector3r::Zero(vertices.cols()); - signed_volumes.resize(components.maxCoeff() + 1); - signed_volumes.setZero(); - for (Eigen::Index f = 0; f < facets.rows(); ++f) { - for (Eigen::Index lv = 0; lv < facets.cols(); ++lv) { - t[lv] = vertices.row(facets(f, lv)); - } - double vol = tetra_signed_volume(t[0], t[1], t[2], t[3]); - signed_volumes(components(f)) -= vol; - } - }; - - const auto& vertices = mesh.get_vertices(); - FacetArray facets; - mesh.export_facets(facets); - - // Orient connected components on the surface - Eigen::Matrix components; - lagrange::internal::bfs_orient(facets, facets, components); + /// Orient to have positive signed volume. + bool positive = true; +}; - // Signed volumes per components - Eigen::Matrix signed_volumes; - component_signed_volumes(vertices, facets, components, signed_volumes); - - // Flip facets in each compoenent with incorrect volume sign - for (Eigen::Index f = 0; f < facets.rows(); ++f) { - if ((positive ? 1 : -1) * signed_volumes(components(f)) < 0) { - facets.row(f) = facets.row(f).reverse().eval(); - } - } - - mesh.import_facets(facets); -} +/// +/// Orient the facets of a mesh so that the signed volume of each connected component is positive or +/// negative. +/// +/// @param[in,out] mesh Mesh whose facets needs to be re-oriented. +/// @param[in] options Orientation options. +/// +/// @tparam Scalar Mesh scalar type. +/// @tparam Index Mesh index type. +/// +template +void orient_outward(lagrange::SurfaceMesh& mesh, const OrientOptions& options = {}); } // namespace lagrange diff --git a/modules/core/include/lagrange/select_facets_by_normal_similarity.h b/modules/core/include/lagrange/select_facets_by_normal_similarity.h index e8ebe58..6485381 100644 --- a/modules/core/include/lagrange/select_facets_by_normal_similarity.h +++ b/modules/core/include/lagrange/select_facets_by_normal_similarity.h @@ -11,253 +11,95 @@ */ #pragma once -#include -#include -#include -#include - -#include +#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS + #include +#endif +#include // needed in options +#include // needed in options namespace lagrange { -/** - * Given a seed vertex, selects faces around it based on the change - * in triangle normals. - */ - -// -// The parameters for the function -// -template -struct SelectFacetsByNormalSimilarityParameters -{ - using Scalar = typename MeshType::Scalar; - using Index = typename MeshType::Index; - - // - // This are input and have to be set - // - // Increasing this would select a larger region - Scalar flood_error_limit = std::numeric_limits::max(); - // This tries to smooth the selection boundary (reduce ears) - bool should_smooth_boundary = true; - - // - // These parameters have default values - // Default value comes from segmentation/Constants.h - // - // If this is not empty, then only faces are considered which are marked as true - std::vector is_facet_selectable = std::vector(); - // Internal param... - Scalar flood_second_to_first_order_limit_ratio = 1.0 / 6.0; - // Number of iterations for smoothing the boundary - Index num_smooth_iterations = 3; - // Use DFS or BFS search - enum class SearchType { BFS = 0, DFS = 1 }; - SearchType search_type = SearchType::DFS; - // Set selectable facets using a vector of indices - void set_selectable_facets(MeshType& mesh_ref, const std::vector& selectable_facets) - { - if (selectable_facets.empty()) { - is_facet_selectable.clear(); - } else { - is_facet_selectable.resize(mesh_ref.get_num_facets(), false); - for (auto facet_id : range_sparse(mesh_ref.get_num_facets(), selectable_facets)) { - is_facet_selectable[facet_id] = true; - } - } - } +/// +/// @defgroup group-surfacemesh-utils Mesh utilities +/// @ingroup group-surfacemesh +/// +/// Various mesh processing utilities +/// +/// @{ -}; // SelectFaceByNormalSimilarityParameters - -// -// Perform the selection -// mesh, input: the mesh -// seed_facet_id, input: the index of the seed -// parameters, input: the parameters above -// returns a vector that contains true for facets that are selected -// -template -std::vector select_facets_by_normal_similarity( - MeshType& mesh, - const typename MeshType::Index seed_facet_id, - const SelectFacetsByNormalSimilarityParameters& parameters) +/// +/// Option struct for selecting facets based on normal similarity. +/// +struct SelectFacetsByNormalSimilarityOptions { - // =================================================== - // Helper typedef and lambdas - // =================================================== - - using VertexType = typename MeshType::VertexType; - using Index = typename MeshType::Index; - using Scalar = typename MeshType::Scalar; - using IndexList = typename MeshType::IndexList; - using AttributeArray = typename MeshType::AttributeArray; - using Parameters = SelectFacetsByNormalSimilarityParameters; - - // Check if a face is selectable - auto is_facet_selectable = [&mesh, ¶meters](const Index facet_id) -> bool { - if (parameters.is_facet_selectable.empty()) { - return true; - } else { - la_runtime_assert(facet_id < mesh.get_num_facets()); - return parameters.is_facet_selectable[facet_id]; - } - }; - - // The energy between two normals - auto normal_direction_error = [](const VertexType& n1, const VertexType& n2) -> Scalar { - return static_cast(1.0) - std::abs(n1.dot(n2)); + /// Increasing this would select a larger region + double flood_error_limit = std::numeric_limits::max(); + + /// + /// There are two types of error limits when flood-search goes from one facet to its neighboring + /// facet: + /// + /// 1. first-order-limit := difference(seed normal, neighboring facet normal) + /// 2. second-order-limit := difference(the one facet's normal, neighboring facet normal) + /// + /// Setting flood_second_to_first_order_limit_ratio > 0.0 allows the selected region to grow on + /// low-curvature areas even though the normals differ from the seed normal. + double flood_second_to_first_order_limit_ratio = 1.0 / 6.0; + + /// The attribute name for the facet normal + std::string_view facet_normal_attribute_name = "@facet_normal"; + + /// Users can specify whether a facet is selectable by an uint8 attribute. + /// e.g. "@is_facet_selectable" + std::optional is_facet_selectable_attribute_name; + + /// The attribute name for the selection output + std::string_view output_attribute_name = "@is_selected"; + + /// select_facets_by_normal_similarity uses either BFS or DFS in its flooding search + enum class SearchType { + BFS = 0, ///< Breadth-First Search + DFS = 1 ///< Depth-First Search }; - // =================================================== - // Let's begin - // current code just copied from - // Segmentation/SuperTriangulation/SingleFloodFromSeed() - // =================================================== - - if (mesh.get_vertex_per_facet() != 3) { - throw std::runtime_error("Input mesh must be a triangle mesh."); - } - - // We need this things - if (!mesh.is_connectivity_initialized()) { - mesh.initialize_connectivity(); - } - if (!mesh.has_facet_attribute("normal")) { - compute_triangle_normal(mesh); - } - - // This would be the final value that will be returned - std::vector is_facet_selected(mesh.get_num_facets(), false); - - // Normals (and the seed normal) - const AttributeArray& facet_normals = mesh.get_facet_attribute("normal"); - const VertexType seed_n = facet_normals.row(seed_facet_id); - - // DFS search utility - std::vector is_facet_processed(mesh.get_num_facets(), false); - std::deque search_queue; - - // - // (0): Add the seed neighbors to the queue to initialize the queue - // - is_facet_processed[seed_facet_id] = true; - is_facet_selected[seed_facet_id] = true; - { - const IndexList& facets_adjacent_to_seed = mesh.get_facets_adjacent_to_facet(seed_facet_id); - for (Index j = 0; j < static_cast(facets_adjacent_to_seed.size()); j++) { - const Index ne_fid = facets_adjacent_to_seed[j]; - const VertexType ne_n = facet_normals.row(ne_fid); - - if ((!is_facet_processed[ne_fid]) && is_facet_selectable(ne_fid)) { - Scalar error = normal_direction_error(seed_n, ne_n); - if (error < parameters.flood_error_limit) { - is_facet_selected[ne_fid] = true; - search_queue.push_back(ne_fid); - } - } - } // End of for over neighbours of seed - } // end of step 0 - - // - // (1): Run through all neighbors, process them, and push them into a queue - // for further flood - // - while (!search_queue.empty()) { - Index fid = invalid(); - if (parameters.search_type == Parameters::SearchType::DFS) { - fid = search_queue.front(); - search_queue.pop_front(); - } else { - fid = search_queue.back(); - search_queue.pop_back(); - } - - const VertexType center_n = facet_normals.row(fid); - - - const IndexList& facets_adjacent_to_candidate = mesh.get_facets_adjacent_to_facet(fid); - for (Index j = 0; j < static_cast(facets_adjacent_to_candidate.size()); j++) { - const Index ne_fid = facets_adjacent_to_candidate[j]; - const VertexType ne_n = facet_normals.row(ne_fid); - - if ((!is_facet_processed[ne_fid]) && is_facet_selectable(ne_fid)) { - const Scalar error_1 = normal_direction_error(seed_n, ne_n); - const Scalar error_2 = normal_direction_error(center_n, ne_n); - is_facet_processed[ne_fid] = true; - - if ((error_1 < parameters.flood_error_limit) && - (error_2 < parameters.flood_error_limit)) { - is_facet_selected[ne_fid] = true; - search_queue.push_back(ne_fid); - } else if ( - error_2 < parameters.flood_error_limit * - parameters.flood_second_to_first_order_limit_ratio) { - // Second order approximation - is_facet_selected[ne_fid] = true; - search_queue.push_back(ne_fid); - } - } // end of is neighbour selectable - } // end of loops over neighbours - } // end of dfs stack - - // - // Smooth the selection - // - if (parameters.should_smooth_boundary) { - for (Index st_iter_2 = 0; st_iter_2 < parameters.num_smooth_iterations; st_iter_2++) { - // Go through boundary faces - check spikes - for (Index i = 0; i < mesh.get_num_facets(); i++) { - const IndexList& facets_adjacent_to_candidate = - mesh.get_facets_adjacent_to_facet(i); - - if (is_facet_selectable(i) && (facets_adjacent_to_candidate.size() > 2)) { - const bool select_flag = is_facet_selected[i]; - Index neighbor_diff_flag = 0; - - const Index ne_fid_0 = facets_adjacent_to_candidate[0]; - const bool ne_select_0 = is_facet_selected[ne_fid_0]; - if (select_flag != ne_select_0) neighbor_diff_flag++; - - const Index ne_fid_1 = facets_adjacent_to_candidate[1]; - const bool ne_select_1 = is_facet_selected[ne_fid_1]; - if (select_flag != ne_select_1) neighbor_diff_flag++; - - const Index ne_fid_2 = facets_adjacent_to_candidate[2]; - const bool ne_select_2 = is_facet_selected[ne_fid_2]; - if (select_flag != ne_select_2) neighbor_diff_flag++; - - if (neighbor_diff_flag > 1) { - const VertexType self_n = facet_normals.row(i); - const VertexType ne_n_0 = facet_normals.row(ne_fid_0); - const VertexType ne_n_1 = facet_normals.row(ne_fid_1); - const VertexType ne_n_2 = facet_normals.row(ne_fid_2); + /// The search type (BFS or DFS) + SearchType search_type = SearchType::DFS; - const Scalar e_0 = normal_direction_error(self_n, ne_n_0); - const Scalar e_1 = normal_direction_error(self_n, ne_n_1); - const Scalar e_2 = normal_direction_error(self_n, ne_n_2); + /// Smoothing boundary of selected region (reduce ears) + int num_smooth_iterations = 3; - const bool e_0_ok = e_0 < parameters.flood_error_limit; - const bool e_1_ok = e_1 < parameters.flood_error_limit; - const bool e_2_ok = e_2 < parameters.flood_error_limit; - if (ne_select_0 && ne_select_1 && (e_0_ok || e_1_ok)) { - is_facet_selected[i] = true; - } else if (ne_select_1 && ne_select_2 && (e_1_ok || e_2_ok)) { - is_facet_selected[i] = true; - } else if (ne_select_2 && ne_select_0 && (e_2_ok || e_0_ok)) { - is_facet_selected[i] = true; - } - } - } // end of non selected facet that has two selectable neighbours - } // end of loop over faces - } // end of for(num_smooth_iterations) - } // end of should smooth boundary +}; // SelectFaceByNormalSimilarityOptions +/** + * Given a seed facet, selects facets around it based on the change in triangle normals. + * + * @param[in, out] mesh The input mesh. + * @param[in] seed_facet_id The index of the seed facet + * @param[in] options Optional arguments (greedy, output_attribute_name). + * + * @tparam Scalar Mesh scalar type. + * @tparam Index Mesh index type. + * + * @return AttributeId The attribute id of the selection, whose attribute name is given + * by options.output_attribute_name. + * + * @note Currently only support triangular mesh and will throw error if + * mesh.is_triangle_mesh() is false! The function will check if the mesh contains + * facet normal by looking for options.facet_normal_attribute_name, and if not + * found, will call compute_facet_normal(meshm, + * options.facet_normal_attribute_name). + * + * @see SelectFacetByNormalSimilarityOptions + */ - return is_facet_selected; +template +AttributeId select_facets_by_normal_similarity( + SurfaceMesh& mesh, + const Index seed_facet_id, + const SelectFacetsByNormalSimilarityOptions& options = {}); -} // select_faces_by_normal_similarity +/// @} } // namespace lagrange diff --git a/modules/core/include/lagrange/types/TransformOptions.h b/modules/core/include/lagrange/types/TransformOptions.h index d8962e7..491f044 100644 --- a/modules/core/include/lagrange/types/TransformOptions.h +++ b/modules/core/include/lagrange/types/TransformOptions.h @@ -32,6 +32,11 @@ struct TransformOptions /// If enabled, tangents and bitangents are normalized after transformation. bool normalize_tangents_bitangents = true; + + /// If enabled, when applying a transform with negative determinant: + /// 1. Normals, tangents and bitangents are flipped. + /// 2. Facet orientations are reversed. + bool reorient = false; }; /// diff --git a/modules/core/include/lagrange/utils/warning.h b/modules/core/include/lagrange/utils/warning.h index 0d0a883..57c9a9c 100644 --- a/modules/core/include/lagrange/utils/warning.h +++ b/modules/core/include/lagrange/utils/warning.h @@ -111,8 +111,16 @@ /// Ignore self move /// @hideinitializer +#if defined(__GNUC__) && __GNUC__ >= 13 +// -Wself-move is introduced in GCC 13. +#define LA_IGNORE_SELF_MOVE_GCC LA_DISABLE_WARNING_GCC(-Wself-move) +#else +#define LA_IGNORE_SELF_MOVE_GCC +#endif +/// @hideinitializer #define LA_IGNORE_SELF_MOVE_WARNING_BEGIN LA_DISABLE_WARNING_BEGIN \ LA_DISABLE_WARNING_CLANG(-Wself-move) \ + LA_IGNORE_SELF_MOVE_GCC \ LA_DISABLE_WARNING_MSVC(26800) /// @hideinitializer #define LA_IGNORE_SELF_MOVE_WARNING_END LA_DISABLE_WARNING_END diff --git a/modules/core/include/lagrange/utils/warnoff.h b/modules/core/include/lagrange/utils/warnoff.h index 53b423c..6dee24c 100644 --- a/modules/core/include/lagrange/utils/warnoff.h +++ b/modules/core/include/lagrange/utils/warnoff.h @@ -69,6 +69,7 @@ #pragma GCC diagnostic ignored "-Wstrict-aliasing" #pragma GCC diagnostic ignored "-Wstringop-overflow" #pragma GCC diagnostic ignored "-Wswitch-default" + #pragma GCC diagnostic ignored "-Wtautological-compare" #pragma GCC diagnostic ignored "-Wundef" #pragma GCC diagnostic ignored "-Wuninitialized" #pragma GCC diagnostic ignored "-Wunused-but-set-variable" diff --git a/modules/core/python/scripts/meshconvert.py b/modules/core/python/scripts/meshconvert.py index caec91e..6c99477 100755 --- a/modules/core/python/scripts/meshconvert.py +++ b/modules/core/python/scripts/meshconvert.py @@ -21,6 +21,7 @@ def parse_args(): parser.add_argument("input_mesh", help="input mesh file") parser.add_argument("output_mesh", help="output mesh file") parser.add_argument("--triangulate", "-t", action="store_true", help="triangulate the mesh") + parser.add_argument("--stitch-vertices", '-s', action="store_true", help="stitch vertices") parser.add_argument( "--logging-level", "-l", @@ -58,7 +59,7 @@ def main(): lagrange.io.save_mesh(args.output_mesh, mesh) else: lagrange.logger.info(f"Loading input mesh from {args.input_mesh}") - mesh = lagrange.io.load_mesh(args.input_mesh) + mesh = lagrange.io.load_mesh(args.input_mesh, stitch_vertices=args.stitch_vertices) if not mesh.is_regular and args.triangulate: lagrange.triangulate_polygonal_facets(mesh) lagrange.logger.info(f"Saving output mesh in {args.output_mesh}") diff --git a/modules/core/python/src/bind_mesh_cleanup.h b/modules/core/python/src/bind_mesh_cleanup.h index a06003b..ebc0d0a 100644 --- a/modules/core/python/src/bind_mesh_cleanup.h +++ b/modules/core/python/src/bind_mesh_cleanup.h @@ -20,6 +20,7 @@ #include #include #include +#include // clang-format off #include @@ -159,6 +160,40 @@ is not. R"(Resolve both vertex and edge nonmanifoldness in a mesh. :param mesh: The input mesh for inplace modification.)"); + + m.def( + "split_long_edges", + [](MeshType& mesh, + float max_edge_length, + bool recursive, + std::optional active_region_attribute, + std::optional edge_length_attribute) { + SplitLongEdgesOptions opts; + opts.max_edge_length = max_edge_length; + opts.recursive = recursive; + if (active_region_attribute.has_value()) + opts.active_region_attribute = active_region_attribute.value(); + if (edge_length_attribute.has_value()) + opts.edge_length_attribute = edge_length_attribute.value(); + split_long_edges(mesh, std::move(opts)); + }, + "mesh"_a, + "max_edge_length"_a = 0.1f, + "recursive"_a = true, + "active_region_attribute"_a = nb::none(), + "edge_length_attribute"_a = nb::none(), + R"(Split edges that are longer than `max_edge_length`. + +:param mesh: The input mesh for inplace modification. +:param max_edge_length: Maximum edge length. Edges longer than this value will be split. +:param recursive: If true, the operation will be applied recursively until no edge is longer than `max_edge_length`. +:param active_region_attribute: The facet attribute name that will be used to determine the active region. + If `None`, all edges will be considered. + Otherwise, only edges that are part of the active region will be split. + The attribute must be of value type `uint8_t`. +:param edge_length_attribute: The edge length attribute name that will be used to determine the edge length. + If `None`, the function will compute the edge lengths. +)"); } } // namespace lagrange::python diff --git a/modules/core/python/src/bind_surface_mesh.h b/modules/core/python/src/bind_surface_mesh.h index f18988a..93a233d 100644 --- a/modules/core/python/src/bind_surface_mesh.h +++ b/modules/core/python/src/bind_surface_mesh.h @@ -1210,11 +1210,18 @@ If not provided, the edges are initialized in an arbitrary order. surface_mesh_class.def( "clone", - [](MeshType& self) -> MeshType { - auto py_self = nb::find(self); - return nb::cast(py_self.attr("__deepcopy__")()); + [](MeshType& self, bool strip) -> MeshType { + if (strip) { + return MeshType::stripped_copy(self); + } else { + auto py_self = nb::find(self); + return nb::cast(py_self.attr("__deepcopy__")()); + } }, - R"(Create a deep copy of this mesh.)"); + "strip"_a = false, + R"(Create a deep copy of this mesh. + +:param strip: If True, strip the mesh of all attributes except for the reserved attributes.)"); } } // namespace lagrange::python diff --git a/modules/core/python/src/bind_utilities.h b/modules/core/python/src/bind_utilities.h index dff8f91..87cf747 100644 --- a/modules/core/python/src/bind_utilities.h +++ b/modules/core/python/src/bind_utilities.h @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -49,6 +50,8 @@ #include #include #include +#include +#include // clang-format off #include @@ -319,6 +322,28 @@ Vertices listed in `cone_vertices` are considered as cone vertices, which is alw :returns: the id of the indexed normal attribute.)"); + using ConstArray3d = nb::ndarray, nb::c_contig, nb::device::cpu>; + m.def( + "compute_pointcloud_pca", + [](ConstArray3d points, bool shift_centroid, bool normalize) { + ComputePointcloudPCAOptions options; + options.shift_centroid = shift_centroid; + options.normalize = normalize; + PointcloudPCAOutput output = + compute_pointcloud_pca({points.data(), points.size()}, options); + return std::make_tuple(output.center, output.eigenvectors, output.eigenvalues); + }, + "points"_a, + "shift_centroid"_a = ComputePointcloudPCAOptions().shift_centroid, + "normalize"_a = ComputePointcloudPCAOptions().normalize, + R"(Compute principal components of a point cloud. + +:param points: Input points. +:param shift_centroid: When true: covariance = (P-centroid)^T (P-centroid), when false: covariance = (P)^T (P). +:param normalize: Should we divide the result by number of points? + +:returns: tuple of (center, eigenvectors, eigenvalues).)"); + m.def( "compute_greedy_coloring", [](MeshType& mesh, @@ -1352,8 +1377,7 @@ A mesh considered as manifold if it is both vertex and edge manifold. if (std::holds_alternative(attribute)) { opt.attribute_id = std::get(attribute); } else { - opt.attribute_id = - mesh.get_attribute_id(std::get(attribute)); + opt.attribute_id = mesh.get_attribute_id(std::get(attribute)); } opt.isovalue = isovalue; opt.keep_below = keep_below; @@ -1535,6 +1559,112 @@ The input mesh must be a triangle mesh. :param output_attribute_name: The output attribute name. If none, cast will replace the input attribute. :returns: The id of the new attribute.)"); + + m.def( + "select_facets_by_normal_similarity", + [](SurfaceMesh& mesh, + Index seed_facet_id, + std::optional flood_error_limit, + std::optional flood_second_to_first_order_limit_ratio, + std::optional facet_normal_attribute_name, + std::optional is_facet_selectable_attribute_name, + std::optional output_attribute_name, + std::optional search_type, + std::optional num_smooth_iterations) { + // Set options in the C++ struct + SelectFacetsByNormalSimilarityOptions options; + if (flood_error_limit.has_value()) + options.flood_error_limit = flood_error_limit.value(); + if (flood_second_to_first_order_limit_ratio.has_value()) + options.flood_second_to_first_order_limit_ratio = flood_second_to_first_order_limit_ratio.value(); + if (facet_normal_attribute_name.has_value()) + options.facet_normal_attribute_name = facet_normal_attribute_name.value(); + if (is_facet_selectable_attribute_name.has_value()) { + options.is_facet_selectable_attribute_name = is_facet_selectable_attribute_name; + } + if (output_attribute_name.has_value()) + options.output_attribute_name = output_attribute_name.value(); + if (search_type.has_value()) { + if (search_type.value() == "BFS") + options.search_type = SelectFacetsByNormalSimilarityOptions::SearchType::BFS; + else if (search_type.value() == "DFS") + options.search_type = SelectFacetsByNormalSimilarityOptions::SearchType::DFS; + else + throw std::runtime_error(fmt::format("Invalid search type: {}", search_type.value())); + } + if (num_smooth_iterations.has_value()) + options.num_smooth_iterations = num_smooth_iterations.value(); + + return select_facets_by_normal_similarity(mesh, seed_facet_id, options); + }, + "mesh"_a, /* `_a` is a literal for nanobind to create nb::args, a required argument */ + "seed_facet_id"_a, + "flood_error_limit"_a = nb::none(), + "flood_second_to_first_order_limit_ratio"_a = nb::none(), + "facet_normal_attribute_name"_a = nb::none(), + "is_facet_selectable_attribute_name"_a = nb::none(), + "output_attribute_name"_a = nb::none(), + "search_type"_a = nb::none(), + "num_smooth_iterations"_a = nb::none(), + R"(Select facets by normal similarity (Pythonic API). + +:param mesh: Input mesh. +:param seed_facet_id: Index of the seed facet. +:param flood_error_limit: Tolerance for normals of the seed and the selected facets. Higher limit leads to larger selected region. +:param flood_second_to_first_order_limit_ratio: Ratio of the flood_error_limit and the tolerance for normals of neighboring selected facets. Higher ratio leads to more curvature in selected region. +:param facet_normal_attribute_name: Attribute name of the facets normal. If the mesh doesn't have this attribute, it will call compute_facet_normal to compute it. +:param is_facet_selectable_attribute_name: If provided, this function will look for this attribute to determine if a facet is selectable. +:param output_attribute_name: Attribute name of whether a facet is selected. +:param search_type: Use 'BFS' for breadth-first search or 'DFS' for depth-first search. +:param num_smooth_iterations: Number of iterations to smooth the boundary of the selected region. + +:returns: Id of the attribute on whether a facet is selected.)", + nb::sig("def select_facets_by_normal_similarity(mesh: SurfaceMesh, " + "seed_facet_id: int, " + "flood_error_limit: float | None = None, " + "flood_second_to_first_order_limit_ratio: float | None = None, " + "facet_normal_attribute_name: str | None = None, " + "is_facet_selectable_attribute_name: str | None = None, " + "output_attribute_name: str | None = None, " + "search_type: typing.Literal['BFS', 'DFS'] | None = None," + "num_smooth_iterations: int | None = None) -> int") + ); + + m.def( + "select_facets_in_frustum", + [](SurfaceMesh& mesh, + std::array, 4> frustum_plane_points, + std::array, 4> frustum_plane_normals, + std::optional greedy, + std::optional output_attribute_name) { + // Set options in the C++ struct + Frustum frustum; + for (size_t i=0; i < 4; ++i) { + frustum.planes[i].point = frustum_plane_points[i]; + frustum.planes[i].normal = frustum_plane_normals[i]; + } + FrustumSelectionOptions options; + if (greedy.has_value()) + options.greedy = greedy.value(); + if (output_attribute_name.has_value()) + options.output_attribute_name = output_attribute_name.value(); + + return select_facets_in_frustum(mesh, frustum, options); + }, + "mesh"_a, + "frustum_plane_points"_a, + "frustum_plane_normals"_a, + "greedy"_a = nb::none(), + "output_attribute_name"_a = nb::none(), + R"(Select facets in a frustum (Pythonic API). + +:param mesh: Input mesh. +:param frustum_plane_points: Four points on each of the frustum planes. +:param frustum_plane_normals: Four normals of each of the frustum planes. +:param greedy: If true, the function returns as soon as the first facet is found. +:param output_attribute_name: Attribute name of whether a facet is selected. + +:returns: Whether any facets got selected.)"); } } // namespace lagrange::python diff --git a/modules/core/python/tests/assets.py b/modules/core/python/tests/assets.py index 4af0cd6..c19dbf7 100644 --- a/modules/core/python/tests/assets.py +++ b/modules/core/python/tests/assets.py @@ -90,6 +90,44 @@ def cube(): return mesh +@pytest.fixture +def cube_triangular(): + vertices = np.array( + [ + [0, 0, 0], + [1, 0, 0], + [1, 1, 0], + [0, 1, 0], + [0, 0, 1], + [1, 0, 1], + [1, 1, 1], + [0, 1, 1], + ], + dtype=float, + ) + facets = np.array( + [ + [0, 3, 2], + [2, 1, 0], + [4, 5, 6], + [6, 7, 4], + [1, 2, 6], + [6, 5, 1], + [4, 7, 3], + [3, 0, 4], + [2, 3, 7], + [7, 6, 2], + [0, 1, 5], + [5, 4, 0], + ], + dtype=np.uint32, + ) + mesh = lagrange.SurfaceMesh() + mesh.vertices = vertices + mesh.facets = facets + return mesh + + @pytest.fixture def cube_with_uv(cube): mesh = cube diff --git a/modules/core/python/tests/test_compute_pointcloud_pca.py b/modules/core/python/tests/test_compute_pointcloud_pca.py new file mode 100644 index 0000000..f96624d --- /dev/null +++ b/modules/core/python/tests/test_compute_pointcloud_pca.py @@ -0,0 +1,39 @@ +# +# Copyright 2024 Adobe. All rights reserved. +# This file is licensed to you under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. You may obtain a copy +# of the License at http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software distributed under +# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS +# OF ANY KIND, either express or implied. See the License for the specific language +# governing permissions and limitations under the License. +# +import lagrange +import numpy as np +import pytest + +a = 0.1 +b = 0.4 +c = 1.2 +def make_points(): + return np.array([ + [a, 0, 0], + [-a, 0, 0], + [0, -b, 0], + [0, b, 0], + [0, 0, c], + [0, 0, -c], + ], dtype=np.float64) + +class TestComputePointcloudPCA: + def test_simple(self): + points = make_points() + print(points.shape) + center, eigenvectors, eigenvalues = lagrange.compute_pointcloud_pca(points) + assert center[0] == 0 + assert center[1] == 0 + assert center[2] == 0 + assert eigenvalues[0] == 2 * a * a + assert eigenvalues[1] == 2 * b * b + assert eigenvalues[2] == 2 * c * c diff --git a/modules/core/python/tests/test_select_facets_by_normal_similarity.py b/modules/core/python/tests/test_select_facets_by_normal_similarity.py new file mode 100644 index 0000000..6596cd4 --- /dev/null +++ b/modules/core/python/tests/test_select_facets_by_normal_similarity.py @@ -0,0 +1,24 @@ +# +# Copyright 2024 Adobe. All rights reserved. +# This file is licensed to you under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. You may obtain a copy +# of the License at http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software distributed under +# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS +# OF ANY KIND, either express or implied. See the License for the specific language +# governing permissions and limitations under the License. +# +import lagrange +import numpy as np +import pytest + +from .assets import * + + +class TestSelectFacetsByNormalSimilarity: + def test_cube(self, single_triangle): + mesh = single_triangle + lagrange.select_facets_by_normal_similarity(mesh, 0, output_attribute_name="is_selected") + is_selected = mesh.attribute("is_selected") + assert is_selected.data[0] == 1 diff --git a/modules/core/python/tests/test_select_facets_in_frustum.py b/modules/core/python/tests/test_select_facets_in_frustum.py new file mode 100644 index 0000000..2ae1043 --- /dev/null +++ b/modules/core/python/tests/test_select_facets_in_frustum.py @@ -0,0 +1,30 @@ +# +# Copyright 2024 Adobe. All rights reserved. +# This file is licensed to you under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. You may obtain a copy +# of the License at http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software distributed under +# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS +# OF ANY KIND, either express or implied. See the License for the specific language +# governing permissions and limitations under the License. +# +import lagrange +import numpy as np +import pytest + +from .assets import * + + +class TestSelectFacetsInFrustum: + def test_cube_big_frustum(self, cube): + mesh = cube + points = np.array([[0, 0, -5566], [0, 0, 5566], [0, -5566, 0], [0, 5566, 0]]) + normals = np.array([[0, 0, 1], [0, 0, -1], [0, 1, 0], [0, -1, 0]]) + assert lagrange.select_facets_in_frustum( + mesh, points, normals, output_attribute_name="is_selected" + ) + is_selected = mesh.attribute("is_selected") + assert ( + mesh.num_facets == np.asarray(is_selected.data).sum() + ) # all facets should be selected diff --git a/modules/core/python/tests/test_split_long_edges.py b/modules/core/python/tests/test_split_long_edges.py new file mode 100644 index 0000000..e74702d --- /dev/null +++ b/modules/core/python/tests/test_split_long_edges.py @@ -0,0 +1,77 @@ +# +# Copyright 2024 Adobe. All rights reserved. +# This file is licensed to you under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. You may obtain a copy +# of the License at http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software distributed under +# the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS +# OF ANY KIND, either express or implied. See the License for the specific language +# governing permissions and limitations under the License. +# +import lagrange +import numpy as np + +from .assets import single_triangle + + +class TestSplitLongEdges: + def test_triangle(self, single_triangle): + mesh = lagrange.SurfaceMesh() + mesh.add_vertex([0, 0, 0]) + mesh.add_vertex([1, 0, 0]) + mesh.add_vertex([0, 1, 0]) + mesh.add_triangle(0, 1, 2) + + # No split + lagrange.split_long_edges(mesh, 10) + assert mesh.num_vertices == 3 + assert mesh.num_facets == 1 + + # Split with 0.5 threshold + lagrange.split_long_edges(mesh, 0.5, recursive=False) + assert mesh.num_vertices == 7 + assert mesh.num_facets == 5 + + # Split with 0.5 threshold recursively + lagrange.split_long_edges(mesh, 0.5, recursive=True) + assert mesh.num_vertices == 9 + assert mesh.num_facets == 9 + + def test_two_triangles(self): + mesh = lagrange.SurfaceMesh() + mesh.add_vertex([0, 0, 0]) + mesh.add_vertex([1, 0, 0]) + mesh.add_vertex([0, 1, 0]) + mesh.add_vertex([1, 1, 0]) + mesh.add_triangle(0, 1, 2) + mesh.add_triangle(2, 1, 3) + + mesh.create_attribute("active", initial_values=[1, 0]) + lagrange.split_long_edges(mesh, 0.5, active_region_attribute="active", recursive=True) + assert mesh.num_facets == 12 + + def test_two_triangles_with_uv(self): + mesh = lagrange.SurfaceMesh() + mesh.add_vertex([0, 0, 0]) + mesh.add_vertex([1, 0, 0]) + mesh.add_vertex([0, 1, 0]) + mesh.add_vertex([1, 1, 0]) + mesh.add_triangle(0, 1, 2) + mesh.add_triangle(2, 1, 3) + + mesh.create_attribute( + "uv", + usage=lagrange.AttributeUsage.UV, + initial_values=np.array([[0, 0], [1, 0], [0, 1], [1, 1]], dtype=np.float32), + initial_indices=np.array([[0, 1, 2], [2, 1, 3]], dtype=np.uint32), + ) + lagrange.split_long_edges(mesh, 0.5, recursive=True) + assert mesh.num_facets == 18 + assert mesh.has_attribute("uv") + + uv_attr = mesh.indexed_attribute("uv") + uv_values = uv_attr.values.data + uv_indices = uv_attr.indices.data + + assert np.allclose(mesh.vertices[mesh.facets.ravel(order="C"), :2], uv_values[uv_indices]) diff --git a/modules/core/src/SurfaceMesh.cpp b/modules/core/src/SurfaceMesh.cpp index d8f6711..a1b8a96 100644 --- a/modules/core/src/SurfaceMesh.cpp +++ b/modules/core/src/SurfaceMesh.cpp @@ -1956,6 +1956,44 @@ void SurfaceMesh::remove_facets(function_ref should_ resize_edges_internal(num_edges_new); } +template +void SurfaceMesh::flip_facets(span facets_to_flip) +{ + std::vector old_to_new_corners(get_num_corners()); + std::iota(old_to_new_corners.begin(), old_to_new_corners.end(), 0); + for (Index f : facets_to_flip) { + const Index c_first = get_facet_corner_begin(f); + const Index c_last = get_facet_corner_end(f); + std::reverse(old_to_new_corners.begin() + c_first, old_to_new_corners.begin() + c_last); + } + reindex_corners_internal( + [&](Index c) { return old_to_new_corners[c]; }, + CornerMappingType::ReversingFacets); +} + +template +void SurfaceMesh::flip_facets(std::initializer_list facets_to_flip) +{ + flip_facets(span(facets_to_flip)); +} + +template +void SurfaceMesh::flip_facets(function_ref should_flip_func) +{ + std::vector old_to_new_corners(get_num_corners()); + std::iota(old_to_new_corners.begin(), old_to_new_corners.end(), 0); + for (Index f = 0; f < get_num_facets(); ++f) { + if (should_flip_func(f)) { + const Index c_first = get_facet_corner_begin(f); + const Index c_last = get_facet_corner_end(f); + std::reverse(old_to_new_corners.begin() + c_first, old_to_new_corners.begin() + c_last); + } + } + reindex_corners_internal( + [&](Index c) { return old_to_new_corners[c]; }, + CornerMappingType::ReversingFacets); +} + namespace { template @@ -2581,6 +2619,24 @@ bool SurfaceMesh::is_boundary_edge(Index e) const return (get_next_corner_around_edge(c) == invalid()); } +template +void SurfaceMesh::foreach_facet_around_facet(Index f, function_ref func) + const +{ + // Loop over edges + for (Index c = get_facet_corner_begin(f); c != get_facet_corner_end(f); ++c) { + Index e = get_corner_edge(c); + auto func_skip_self = [f, &func](Index fi) { + if (fi != f) { + func(fi); + } + }; + + // Loop over facets around edge + foreach_facet_around_edge(e, func_skip_self); + } +} + template void SurfaceMesh::foreach_facet_around_vertex( Index v, @@ -2765,101 +2821,6 @@ auto SurfaceMesh::reindex_facets_internal(span old_t const Index num_corners = get_num_corners(); const Index num_facets = get_num_facets(); - auto reindex_edges = [&](auto old_to_new_edges) { - // Update content of EdgeIndex attributes - seq_foreach_attribute_write(*this, [&](auto&& attr) { - update_element_index( - attr, - AttributeUsage::EdgeIndex, - num_edges, - old_to_new_edges); - }); - - // Move edge attributes - seq_foreach_attribute_write(*this, [&](auto&& attr) { - for (Index e = 0; e < num_edges; ++e) { - if (old_to_new_edges(e) != invalid()) { - if (old_to_new_edges(e) != e) { - auto from = attr.get_row(e); - auto to = attr.ref_row(old_to_new_edges(e)); - std::copy(from.begin(), from.end(), to.begin()); - } - } - } - }); - }; - - auto reindex_corners = [&](auto old_to_new_corners) { - // Repair connectivity chains and skip over deleted corners - std::array, 2> ids = { - {{{m_reserved_ids.vertex_to_first_corner(), - m_reserved_ids.next_corner_around_vertex()}}, - {{m_reserved_ids.edge_to_first_corner(), m_reserved_ids.next_corner_around_edge()}}}}; - for (auto [id_first, id_next] : ids) { - if (id_first != invalid_attribute_id()) { - auto first_corner = ref_attribute(id_first).ref_all(); - auto next_corner = ref_attribute(id_next).ref_all(); - auto next_valid_corner = [&](Index ci) { - // Find the next corner that is not discarded - Index cj = next_corner[ci]; - while (cj != invalid() && old_to_new_corners(cj) == invalid()) { - cj = next_corner[cj]; - } - return cj; - }; - for (Index& c0 : first_corner) { - if (id_first == m_reserved_ids.vertex_to_first_corner() && - c0 == invalid()) { - // Isolated vertex, skip - continue; - } - if (old_to_new_corners(c0) == invalid()) { - // Replace head by next valid corner - c0 = next_valid_corner(c0); - } - for (Index c = c0; c != invalid(); c = next_corner[c]) { - // Iterate over the singly linked list and replace invalid items - if (next_corner[c] != invalid()) { - next_corner[c] = next_valid_corner(c); - } - } - } - } - } - - // Update content of CornerIndex attributes - seq_foreach_attribute_write(*this, [&](auto&& attr) { - update_element_index( - attr, - AttributeUsage::CornerIndex, - num_corners, - old_to_new_corners); - }); - - // Move corner & indexed attributes - auto move_corner_attributes = [&](auto&& attr) { - for (Index c = 0; c < num_corners; ++c) { - if (old_to_new_corners(c) != invalid()) { - if (old_to_new_corners(c) != c) { - auto from = attr.get_row(c); - auto to = attr.ref_row(old_to_new_corners(c)); - std::copy(from.begin(), from.end(), to.begin()); - } - } - } - }; - seq_foreach_attribute_write( - *this, - [&](auto&& attr) { - using AttributeType = std::decay_t; - if constexpr (AttributeType::IsIndexed) { - move_corner_attributes(attr.indices()); - } else { - move_corner_attributes(attr); - } - }); - }; - // Compute corner remapping Index new_num_corners = 0; if (is_hybrid()) { @@ -2874,7 +2835,9 @@ auto SurfaceMesh::reindex_facets_internal(span old_t } } } - reindex_corners([&](Index c) noexcept { return old_to_new_corners[c]; }); + reindex_corners_internal( + [&](Index c) noexcept { return old_to_new_corners[c]; }, + CornerMappingType::RemovingFacets); } else { const Index nvpf = get_vertex_per_facet(); for (Index f = 0; f < num_facets; ++f) { @@ -2882,16 +2845,42 @@ auto SurfaceMesh::reindex_facets_internal(span old_t new_num_corners += nvpf; } } - reindex_corners([&](Index c) { - const Index new_f = old_to_new_facets[c / nvpf]; - const Index lv = c % nvpf; - if (new_f == invalid()) { - return invalid(); - } else { - return new_f * nvpf + lv; + reindex_corners_internal( + [&](Index c) { + const Index new_f = old_to_new_facets[c / nvpf]; + const Index lv = c % nvpf; + if (new_f == invalid()) { + return invalid(); + } else { + return new_f * nvpf + lv; + } + }, + CornerMappingType::RemovingFacets); + } + + auto reindex_edges = [&](auto old_to_new_edges) { + // Update content of EdgeIndex attributes + seq_foreach_attribute_write(*this, [&](auto&& attr) { + update_element_index( + attr, + AttributeUsage::EdgeIndex, + num_edges, + old_to_new_edges); + }); + + // Move edge attributes + seq_foreach_attribute_write(*this, [&](auto&& attr) { + for (Index e = 0; e < num_edges; ++e) { + if (old_to_new_edges(e) != invalid()) { + if (old_to_new_edges(e) != e) { + auto from = attr.get_row(e); + auto to = attr.ref_row(old_to_new_edges(e)); + std::copy(from.begin(), from.end(), to.begin()); + } + } } }); - } + }; // Compute edge remapping Index new_num_edges = 0; @@ -2910,6 +2899,7 @@ auto SurfaceMesh::reindex_facets_internal(span old_t reindex_edges([&](Index e) noexcept { return old_to_new_edges[e]; }); } + // Update content of FacetIndex attributes auto remap_f = [&](Index i) { return old_to_new_facets[i]; }; seq_foreach_attribute_write(*this, [&](auto&& attr) { @@ -2932,6 +2922,162 @@ auto SurfaceMesh::reindex_facets_internal(span old_t return std::make_pair(new_num_corners, new_num_edges); } +template +void SurfaceMesh::reindex_corners_internal( + function_ref old_to_new_corners, + CornerMappingType mapping_type) +{ + const Index num_facets = get_num_facets(); + const Index num_corners = get_num_corners(); + + // + // When reversing facets, we need to remap the following attributes differently: + // 1. $corner_to_edge + // 2. $edge_to_first_corner + // 3. $next_corner_around_edge + // 4. $facet_to_first_corner + // + // For attributes (1)-(3), the corners are used to identify an edge, when reversing facets we + // need to use the previous corner around the facet to denote the same edge instead. + // + // For attribute (4), we do not remap it at all, since the first corner of the facet remains + // the same. + // + // Ideally we'd have a different ElementType for FacetVertex and FacetEdge, but for now we're + // just hardcoding this behavior. + // + std::vector corner_to_prev_corner; + if (mapping_type == CornerMappingType::ReversingFacets && + m_reserved_ids.edge_to_first_corner() != invalid_attribute_id()) { + corner_to_prev_corner.resize(num_corners); + std::iota(corner_to_prev_corner.begin(), corner_to_prev_corner.end(), 0); + for (Index f = 0; f < num_facets; ++f) { + const Index c0 = get_facet_corner_begin(f); + if (old_to_new_corners(c0) == c0) { + // Facet not flipped, continue + continue; + } + const Index nv = get_facet_size(f); + for (Index lv = 0; lv < nv; ++lv) { + const Index c = c0 + lv; + corner_to_prev_corner[c] = c0 + (lv + nv - 1) % nv; + } + } + } + + auto old_to_new_corners_shifted = [&](Index c) { + if (mapping_type == CornerMappingType::ReversingFacets && !corner_to_prev_corner.empty()) { + return corner_to_prev_corner[old_to_new_corners(c)]; + } else { + return old_to_new_corners(c); + } + }; + + // Repair connectivity chains and skip over deleted corners. + std::array, 2> ids = { + {{{m_reserved_ids.vertex_to_first_corner(), m_reserved_ids.next_corner_around_vertex()}}, + {{m_reserved_ids.edge_to_first_corner(), m_reserved_ids.next_corner_around_edge()}}}}; + for (auto [id_first, id_next] : ids) { + if (id_first != invalid_attribute_id()) { + auto first_corner = ref_attribute(id_first).ref_all(); + auto next_corner = ref_attribute(id_next).ref_all(); + auto next_valid_corner = [&](Index ci) { + // Find the next corner that is not discarded + Index cj = next_corner[ci]; + while (cj != invalid() && old_to_new_corners(cj) == invalid()) { + cj = next_corner[cj]; + } + return cj; + }; + for (Index& c0 : first_corner) { + if (id_first == m_reserved_ids.vertex_to_first_corner() && c0 == invalid()) { + // Isolated vertex, skip + continue; + } + if (old_to_new_corners(c0) == invalid()) { + // Replace head by next valid corner + c0 = next_valid_corner(c0); + } + for (Index c = c0; c != invalid(); c = next_corner[c]) { + // Iterate over the singly linked list and replace invalid items + if (next_corner[c] != invalid()) { + next_corner[c] = next_valid_corner(c); + } + } + } + } + } + + seq_foreach_named_attribute_write(*this, [&](auto name, auto&& attr) { + if (mapping_type == CornerMappingType::ReversingFacets && + name == s_reserved_names.facet_to_first_corner()) { + // Don't remap the "first corner" of a facet when reversing facet indices + return; + } + if (mapping_type == CornerMappingType::ReversingFacets && + (name == s_reserved_names.edge_to_first_corner() || + name == s_reserved_names.next_corner_around_edge())) { + update_element_index( + attr, + AttributeUsage::CornerIndex, + num_corners, + old_to_new_corners_shifted); + } else { + update_element_index( + attr, + AttributeUsage::CornerIndex, + num_corners, + old_to_new_corners); + } + }); + + // Move corner & indexed attributes + auto move_corner_attributes = [&](auto&& attr, bool shift) { + using AttributeType = std::decay_t; + using ValueType = typename AttributeType::ValueType; + std::vector tmp; + if (mapping_type == CornerMappingType::ReversingFacets) { + tmp.resize(attr.get_num_channels()); + } + for (Index c = 0; c < num_corners; ++c) { + Index new_c = old_to_new_corners(c); + if (shift) { + new_c = corner_to_prev_corner[new_c]; + } + if (new_c != invalid()) { + if (new_c != c) { + auto to = attr.ref_row(new_c); + if (mapping_type == CornerMappingType::RemovingFacets) { + auto from = attr.get_row(c); + std::copy(from.begin(), from.end(), to.begin()); + } else { + if (new_c < c) { + auto from = attr.ref_row(c); + std::copy(to.begin(), to.end(), tmp.begin()); + std::copy(from.begin(), from.end(), to.begin()); + std::copy(tmp.begin(), tmp.end(), from.begin()); + } + } + } + } + } + }; + seq_foreach_named_attribute_write( + *this, + [&](auto name, auto&& attr) { + using AttributeType = std::decay_t; + bool shift = + (mapping_type == CornerMappingType::ReversingFacets && + (name == s_reserved_names.corner_to_edge() || + name == s_reserved_names.next_corner_around_edge())); + if constexpr (AttributeType::IsIndexed) { + move_corner_attributes(attr.indices(), shift); + } else { + move_corner_attributes(attr, shift); + } + }); +} + template template void SurfaceMesh::resize_elements_internal(Index num_elements) diff --git a/modules/core/src/compute_pointcloud_pca.cpp b/modules/core/src/compute_pointcloud_pca.cpp new file mode 100644 index 0000000..b523267 --- /dev/null +++ b/modules/core/src/compute_pointcloud_pca.cpp @@ -0,0 +1,89 @@ +/* + * Copyright 2024 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ + +#include +#include +#include + +#include +#include + +namespace lagrange { + +template +PointcloudPCAOutput compute_pointcloud_pca( + span points, + ComputePointcloudPCAOptions options) +{ + // Prechecks + constexpr uint32_t dimension = 3; + la_runtime_assert(points.size() % dimension == 0); + + const size_t num_points = points.size() / dimension; + la_runtime_assert(num_points >= 2, "There must be at least two points"); + + using MatrixXS = Eigen::Matrix; + using RowVector3S = Eigen::Matrix; + + Eigen::Map vertices(points.data(), num_points, dimension); + + MatrixXS covariance; + RowVector3S center; + + if (options.shift_centroid) { + center = vertices.colwise().mean(); + } else { + center = RowVector3S::Zero(); + } + + if (options.normalize) { + // We may instead divide by points.rows()-1 which applies the Bessel's correction. + // https://en.wikipedia.org/wiki/Bessel%27s_correction + covariance = ((vertices.rowwise() - center).transpose() * (vertices.rowwise() - center)) / + vertices.rows(); + } else { + covariance = ((vertices.rowwise() - center).transpose() * (vertices.rowwise() - center)); + } + + // SelfAdjointEigenSolver is guaranteed to return the eigenvalues sorted + // in increasing order. Note that this is not true for the general EigenValueSolver. + Eigen::SelfAdjointEigenSolver 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"); + + // make sure components follow the right hand rule + auto eigenvectors = eigs.eigenvectors(); + auto eigenvalues = eigs.eigenvalues(); + if (eigenvectors.determinant() < 0.) { + eigenvectors.col(0) *= -1; + } + + PointcloudPCAOutput output = PointcloudPCAOutput(); + for (size_t i = 0; i < dimension; ++i) { + output.center[i] = center(i); + for (size_t j = 0; j < dimension; ++j) { + output.eigenvectors[i][j] = eigenvectors(i, j); + } + output.eigenvalues[i] = eigenvalues(i); + } + + return output; +} + +#define LA_X_compute_pointcloud_pca(_, Scalar) \ + template LA_CORE_API PointcloudPCAOutput compute_pointcloud_pca( \ + span, \ + ComputePointcloudPCAOptions); +LA_SURFACE_MESH_SCALAR_X(compute_pointcloud_pca, 0) + +} // namespace lagrange diff --git a/modules/core/src/mesh_cleanup/split_long_edges.cpp b/modules/core/src/mesh_cleanup/split_long_edges.cpp new file mode 100644 index 0000000..24757f7 --- /dev/null +++ b/modules/core/src/mesh_cleanup/split_long_edges.cpp @@ -0,0 +1,361 @@ +/* + * Copyright 2024 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "split_triangle.h" + +// clang-format off +#include +#include +#include +#include +// clang-format on + +namespace lagrange { + +namespace { + +template +void interpolate_row( + Eigen::MatrixBase& data, + Index row_to, + Index row_from_1, + Index row_from_2, + Scalar t) +{ + la_debug_assert(row_to < static_cast(data.rows())); + la_debug_assert(row_from_1 < static_cast(data.rows())); + la_debug_assert(row_from_2 < static_cast(data.rows())); + using ValueType = typename Derived::Scalar; + + if constexpr (std::is_integral_v) { + data.row(row_to) = (data.row(row_from_1).template cast() * (1 - t) + + data.row(row_from_2).template cast() * t) + .array() + .round() + .template cast() + .eval(); + } else { + data.row(row_to) = data.row(row_from_1) * (1 - t) + data.row(row_from_2) * t; + } +} + +} // namespace + +template +void split_long_edges(SurfaceMesh& mesh, SplitLongEdgesOptions options) +{ + la_runtime_assert(mesh.is_triangle_mesh(), "Input mesh is not a triangle mesh."); + mesh.initialize_edges(); + + const Index dim = mesh.get_dimension(); + const Index num_input_vertices = mesh.get_num_vertices(); + const Index num_edges = mesh.get_num_edges(); + const Index num_input_facets = mesh.get_num_facets(); + + constexpr std::string_view internal_active_region_attribute_name = + "@__active_region_internal__"; + AttributeId internal_active_region_attribute_id = invalid(); + if (!options.active_region_attribute.empty()) { + // Cast attribute to a known type for easy access. + internal_active_region_attribute_id = cast_attribute( + mesh, + options.active_region_attribute, + internal_active_region_attribute_name); + } + + // An active facet means its edges are checked against the edge length threshold. + auto is_active = [&](Index fid) { + if (internal_active_region_attribute_id == invalid()) return true; + auto active_view = + attribute_vector_view(mesh, internal_active_region_attribute_id); + return static_cast(active_view[fid]); + }; + + // An edge is active if it is adjacent to an active facet. + auto is_edge_active = [&](Index eid) { + bool r = false; + mesh.foreach_facet_around_edge(eid, [&](auto fid) { r |= is_active(fid); }); + return r; + }; + + // A facet may be involved in splitting if one of its edges is active. + auto is_facet_involved = [&](Index fid) { + auto e0 = mesh.get_edge(fid, 0); + auto e1 = mesh.get_edge(fid, 1); + auto e2 = mesh.get_edge(fid, 2); + return is_edge_active(e0) || is_edge_active(e1) || is_edge_active(e2); + }; + + EdgeLengthOptions edge_length_options; + edge_length_options.output_attribute_name = options.edge_length_attribute; + auto edge_length_attr_id = compute_edge_lengths(mesh, edge_length_options); + auto edge_lengths = attribute_vector_view(mesh, edge_length_attr_id); + auto input_vertices = vertex_view(mesh); + + // Computer split locations + Index estimated_num_additional_vertices = + static_cast(std::ceil(edge_lengths.sum() / options.max_edge_length)); + std::vector additional_vertices; + additional_vertices.reserve(estimated_num_additional_vertices * dim); + std::vector> additional_vertex_sources; + additional_vertex_sources.reserve(estimated_num_additional_vertices); + std::vector edge_split_indices(num_edges + 1); + edge_split_indices[0] = 0; + for (Index eid = 0; eid < num_edges; eid++) { + auto l = edge_lengths[eid]; + if (!is_edge_active(eid) || l <= options.max_edge_length) { + edge_split_indices[eid + 1] = edge_split_indices[eid]; + continue; + } + + auto n = static_cast(std::ceil(l / options.max_edge_length)); + assert(n > 1); + + auto [v0, v1] = mesh.get_edge_vertices(eid); + auto p0 = input_vertices.row(v0).eval(); + auto p1 = input_vertices.row(v1).eval(); + for (Index i = 1; i < n; i++) { + Scalar t = static_cast(i) / n; + auto p = ((1 - t) * p0 + t * p1).eval(); + additional_vertices.insert(additional_vertices.end(), p.data(), p.data() + dim); + additional_vertex_sources.emplace_back(v0, v1, t); + } + + edge_split_indices[eid + 1] = edge_split_indices[eid] + n - 1; + }; + if (additional_vertices.empty()) return; + + // Update vertices + const Index num_additional_vertices = static_cast(additional_vertices.size()) / dim; + mesh.add_vertices( + num_additional_vertices, + {additional_vertices.data(), additional_vertices.size()}); + auto output_vertices = vertex_view(mesh); + + // Interpolate vertex attributes for newly added vertices + par_foreach_named_attribute_write( + mesh, + [&](std::string_view name, auto&& attr) { + using ValueType = typename std::decay_t::ValueType; + if (mesh.attr_name_is_reserved(name)) return; + auto data = matrix_ref(attr); + static_assert(std::is_same_v::Scalar>); + for (Index i = 0; i < num_additional_vertices; i++) { + Index v0, v1; + Scalar t; + std::tie(v0, v1, t) = additional_vertex_sources[i]; + interpolate_row(data, num_input_vertices + i, v0, v1, t); + } + }); + + // Split triangles + std::vector original_triangle_index; + std::vector split_triangles; + std::vector split_triangles_offsets; + original_triangle_index.reserve(num_input_facets); + split_triangles.reserve(num_input_facets); + split_triangles_offsets.reserve(num_input_facets + 1); + split_triangles_offsets.push_back(0); + + // Compute the number of additional facets + Index num_additional_facets = 0; + for (Index fid = 0; fid < num_input_facets; fid++) { + if (!is_facet_involved(fid)) continue; + Index chain_size = 3; + for (Index i = 0; i < 3; i++) { + Index eid = mesh.get_edge(fid, i); + la_debug_assert(eid != invalid()); + Index n = edge_split_indices[eid + 1] - edge_split_indices[eid]; + chain_size += n; + } + if (chain_size == 3) continue; + num_additional_facets += chain_size - 2; + split_triangles_offsets.push_back(num_additional_facets * 3); + original_triangle_index.push_back(fid); + } + split_triangles.resize(num_additional_facets * 3, invalid()); + + tbb::enumerable_thread_specific> chain_buffer; + tbb::enumerable_thread_specific> queue_buffer; + tbb::enumerable_thread_specific> visited_buffer; + tbb::parallel_for(size_t(0), original_triangle_index.size(), [&](size_t idx) { + auto fid = original_triangle_index[idx]; + auto& chain = chain_buffer.local(); + chain.clear(); + chain.reserve(64); + + auto corners = mesh.get_facet_vertices(fid); + Index corner_index[3] = {invalid(), invalid(), invalid()}; + for (Index i = 0; i < 3; i++) { + Index v0 = corners[i]; + chain.push_back(v0); + corner_index[i] = static_cast(chain.size() - 1); + + Index eid = mesh.get_edge(fid, i); + la_debug_assert(eid != invalid()); + Index n = edge_split_indices[eid + 1] - edge_split_indices[eid]; + if (mesh.get_edge_vertices(eid)[0] == v0) { + // Edge is oriented in the same direction as the triangle. + for (Index j = 0; j < n; j++) { + chain.push_back(num_input_vertices + edge_split_indices[eid] + j); + } + } else { + // Edge is oriented in the opposite direction. + la_debug_assert(mesh.get_edge_vertices(eid)[1] == v0); + for (Index j = 0; j < n; j++) { + chain.push_back(num_input_vertices + edge_split_indices[eid] + n - 1 - j); + } + } + } + la_debug_assert(chain.size() > 3); + size_t num_sub_triangles = chain.size() - 2; + size_t curr_size = split_triangles_offsets[idx]; + span sub_triangulation(split_triangles.data() + curr_size, num_sub_triangles * 3); + auto& visited = visited_buffer.local(); + visited.resize(chain.size() * 3); + split_triangle( + static_cast(output_vertices.rows()), + span(output_vertices.data(), output_vertices.size()), + chain, + visited, + queue_buffer.local(), + corner_index[0], + corner_index[1], + corner_index[2], + sub_triangulation); + }); + + // Clear edge data structure since the edge data is costly to maintain, + // and we will be updating triangulation anyway. + mesh.clear_edges(); + + la_debug_assert(split_triangles.size() % 3 == 0); + mesh.add_triangles(num_additional_facets, {split_triangles.data(), split_triangles.size()}); + + // Propagate facet attributes + par_foreach_named_attribute_write( + mesh, + [&](std::string_view name, auto&& attr) { + if (mesh.attr_name_is_reserved(name)) return; + auto data = matrix_ref(attr); + for (size_t i = 0; i < original_triangle_index.size(); i++) { + for (size_t j = split_triangles_offsets[i] / 3; + j < split_triangles_offsets[i + 1] / 3; + j++) { + data.row(num_input_facets + static_cast(j)) = + data.row(original_triangle_index[i]); + } + } + }); + + auto map_corner_attribute = [&](auto&& data, auto&& corner_to_index) { + auto facets = facet_view(mesh); + for (size_t i = 0; i < original_triangle_index.size(); i++) { + for (size_t j = split_triangles_offsets[i] / 3; j < split_triangles_offsets[i + 1] / 3; + j++) { + Index fid = num_input_facets + static_cast(j); + + for (Index k = 0; k < 3; k++) { + Index curr_cid = fid * 3 + k; + Index vid = facets(fid, k); + if (vid < num_input_vertices) { + // Corner comes from a corner in the original triangulation. + Index ori_cid = invalid(); + for (Index l = 0; l < 3; l++) { + if (vid == facets(original_triangle_index[i], l)) { + ori_cid = original_triangle_index[i] * 3 + l; + } + } + la_debug_assert(ori_cid != invalid()); + data.row(corner_to_index(curr_cid)) = data.row(corner_to_index(ori_cid)); + } else { + // Corner comes from a split edge. + auto [v0, v1, t] = additional_vertex_sources[vid - num_input_vertices]; + + Index ori_c0 = invalid(); + Index ori_c1 = invalid(); + for (Index l = 0; l < 3; l++) { + if (v0 == facets(original_triangle_index[i], l)) { + ori_c0 = original_triangle_index[i] * 3 + l; + } + if (v1 == facets(original_triangle_index[i], l)) { + ori_c1 = original_triangle_index[i] * 3 + l; + } + } + la_debug_assert(ori_c0 != invalid()); + la_debug_assert(ori_c1 != invalid()); + interpolate_row( + data, + corner_to_index(curr_cid), + corner_to_index(ori_c0), + corner_to_index(ori_c1), + t); + } + } + } + } + }; + + if (internal_active_region_attribute_id != invalid()) { + // Remove the internal active region attribute + mesh.delete_attribute(internal_active_region_attribute_name); + } + + // Convert indexed attributes to corner attributes + std::vector indexed_attribute_ids; + seq_foreach_named_attribute_read( + mesh, + [&](std::string_view name, [[maybe_unused]] auto&& attr) { + if (mesh.attr_name_is_reserved(name)) return; + indexed_attribute_ids.push_back(mesh.get_attribute_id(name)); + }); + for (auto& attr_id : indexed_attribute_ids) { + attr_id = map_attribute_in_place(mesh, attr_id, AttributeElement::Corner); + } + + // Propagate corner attributes + par_foreach_named_attribute_write( + mesh, + [&](std::string_view name, auto&& attr) { + if (mesh.attr_name_is_reserved(name)) return; + auto data = matrix_ref(attr); + map_corner_attribute(data, [&](Index cid) { return cid; }); + }); + + // Map back to indexed attributes + for (auto attr_id : indexed_attribute_ids) { + map_attribute_in_place(mesh, attr_id, AttributeElement::Indexed); + } + + mesh.remove_facets(original_triangle_index); + + if (options.recursive) { + split_long_edges(mesh, options); + } +} + +#define LA_X_split_long_edges(_, Scalar, Index) \ + template LA_CORE_API void split_long_edges( \ + SurfaceMesh&, \ + SplitLongEdgesOptions); +LA_SURFACE_MESH_X(split_long_edges, 0) +} // namespace lagrange diff --git a/modules/core/src/mesh_cleanup/split_triangle.cpp b/modules/core/src/mesh_cleanup/split_triangle.cpp new file mode 100644 index 0000000..13f811f --- /dev/null +++ b/modules/core/src/mesh_cleanup/split_triangle.cpp @@ -0,0 +1,204 @@ +/* + * Copyright 2024 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ +#include +#include +#include +#include + +#include "split_triangle.h" + +#include + +#include +#include + +namespace lagrange { + +template +void split_triangle( + size_t num_points, + span points, + span chain, + span visited_buffer, + std::vector& queue_buffer, + const Index v0, + const Index v1, + const Index v2, + span triangulation) +{ + la_runtime_assert(points.size() % num_points == 0); + la_runtime_assert(triangulation.size() % 3 == 0); + la_runtime_assert( + triangulation.size() / 3 + 2 == chain.size(), + "Inconsistent triangulation size vs chian size."); + la_runtime_assert(visited_buffer.size() == chain.size() * 3); + Eigen::Map> + vertices(points.data(), num_points, points.size() / num_points); + + const Index chain_size = static_cast(chain.size()); + auto next = [&chain_size](Index i) { return (i + 1) % chain_size; }; + auto prev = [&chain_size](Index i) { return (i + chain_size - 1) % chain_size; }; + auto sq_length = [&](Index vi, Index vj) { + return (vertices.row(chain[vi]) - vertices.row(chain[vj])).squaredNorm(); + }; + auto is_corner = [v0, v1, v2](Index v) { return v == v0 || v == v1 || v == v2; }; + + constexpr Scalar NOT_USED = invalid(); + Eigen::Matrix candidates; + // clang-format off + candidates << + v0, next(v0), prev(v0), 0, 0, 0, + v1, next(v1), prev(v1), 0, 0, 0, + v2, next(v2), prev(v2), 0, 0, 0; + // clang-format on + + Eigen::Matrix candidate_lengths; + // clang-format off + candidate_lengths << + sq_length(candidates(0, 1), candidates(0, 2)), NOT_USED, + sq_length(candidates(1, 1), candidates(1, 2)), NOT_USED, + sq_length(candidates(2, 1), candidates(2, 2)), NOT_USED; + // clang-format on + + auto index_comp = [&candidate_lengths](Index i, Index j) { + // Use greater than operator so the heap is min heap. + return candidate_lengths.row(i).minCoeff() > candidate_lengths.row(j).minCoeff(); + }; + auto& Q = queue_buffer; + Q.clear(); + Q.push_back(0); + Q.push_back(1); + Q.push_back(2); + std::make_heap(Q.begin(), Q.end(), index_comp); + + Eigen::Map> visited( + visited_buffer.data(), + chain_size, + 3); + visited.setZero(); + visited(v0, 0) = 1; + visited(v1, 1) = 1; + visited(v2, 2) = 1; + + size_t num_generated_triangles = 0; + while (!Q.empty()) { + Index idx = Q.front(); + std::pop_heap(Q.begin(), Q.end(), index_comp); + Q.pop_back(); + Index selection = 0; + if (candidate_lengths(idx, 0) != NOT_USED && + candidate_lengths(idx, 0) <= candidate_lengths(idx, 1)) { + selection = 0; + } else if ( + candidate_lengths(idx, 1) != NOT_USED && + candidate_lengths(idx, 1) <= candidate_lengths(idx, 0)) { + selection = 1; + } else { + // Select neither candidate from this corner. + continue; + } + la_runtime_assert(candidate_lengths.row(idx).minCoeff() <= candidate_lengths.minCoeff()); + + Index base_v = candidates(idx, selection * 3); + Index right_v = candidates(idx, selection * 3 + 1); + Index left_v = candidates(idx, selection * 3 + 2); + la_runtime_assert(base_v >= 0 && base_v < chain_size); + la_runtime_assert(right_v >= 0 && right_v < chain_size); + la_runtime_assert(left_v >= 0 && left_v < chain_size); + la_runtime_assert(visited(base_v, idx) >= 1); + + // A special case. + if (is_corner(base_v) && is_corner(right_v) && is_corner(left_v)) { + if (chain_size != 3) { + candidate_lengths(idx, selection) = NOT_USED; + Q.push_back(idx); + std::push_heap(Q.begin(), Q.end(), index_comp); + continue; + } + } + + if (visited.row(base_v).sum() > 1 || visited(right_v, idx) > 1 || + visited(left_v, idx) > 1) { + // This canadidate is invalid. + candidate_lengths(idx, selection) = NOT_USED; + Q.push_back(idx); + std::push_heap(Q.begin(), Q.end(), index_comp); + continue; + } + + visited(base_v, idx) = 2; + visited(right_v, idx) = 1; + visited(left_v, idx) = 1; + triangulation[num_generated_triangles * 3] = chain[base_v]; + triangulation[num_generated_triangles * 3 + 1] = chain[right_v]; + triangulation[num_generated_triangles * 3 + 2] = chain[left_v]; + num_generated_triangles++; + + // Update candidate from this corner. + if (visited.row(right_v).sum() == 1) { + Index right_to_right = next(right_v); + candidate_lengths(idx, 0) = sq_length(left_v, right_to_right); + candidates(idx, 0) = right_v; + candidates(idx, 1) = right_to_right; + candidates(idx, 2) = left_v; + } else { + candidate_lengths(idx, 0) = NOT_USED; + } + + if (visited.row(left_v).sum() == 1) { + Index left_to_left = prev(left_v); + candidate_lengths(idx, 1) = sq_length(right_v, left_to_left); + candidates(idx, 3) = left_v; + candidates(idx, 4) = right_v; + candidates(idx, 5) = left_to_left; + } else { + candidate_lengths(idx, 1) = NOT_USED; + } + Q.push_back(idx); + std::push_heap(Q.begin(), Q.end(), index_comp); + } + + Eigen::Matrix visited_sum = + (visited.array() > 0).rowwise().count(); + if ((visited_sum.array() > 1).count() == 3) { + Index f[3]; + Index count = 0; + for (Index i = 0; i < chain_size; i++) { + if (visited_sum[i] > 1) { + la_runtime_assert(count < 3); + f[count] = chain[i]; + count++; + } + } + la_runtime_assert(count == 3); + triangulation[num_generated_triangles * 3] = f[0]; + triangulation[num_generated_triangles * 3 + 1] = f[1]; + triangulation[num_generated_triangles * 3 + 2] = f[2]; + num_generated_triangles++; + } + la_debug_assert(num_generated_triangles + 2 == chain.size()); +} + +#define LA_X_split_triangle(_, Scalar, Index) \ + template LA_CORE_API void split_triangle( \ + size_t, \ + span, \ + span, \ + span, \ + std::vector&, \ + const Index, \ + const Index, \ + const Index, \ + span); +LA_SURFACE_MESH_X(split_triangle, 0) + +} // namespace lagrange diff --git a/modules/core/src/mesh_cleanup/split_triangle.h b/modules/core/src/mesh_cleanup/split_triangle.h new file mode 100644 index 0000000..44a991c --- /dev/null +++ b/modules/core/src/mesh_cleanup/split_triangle.h @@ -0,0 +1,49 @@ +/* + * Copyright 2024 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ +#pragma once + +#include + +#include + +namespace lagrange { + +/// +/// Split a triangle into smaller triangles based on the chain of splitting +/// points on the edges. +/// +/// @param num_points Total number of points. +/// @param points Coordinates of the points. +/// @param chain Chain of indices into points that iterates all vertices and +/// splitting points in counterclockwise order around the triangle. +/// @param visited_buffer Temporary buffer that should be of size 3 * chain.size(). +/// @param queue_buffer Temporary buffer used for priority queue. +/// @param v0 Index of the first corner of the input triangle in `chain`. +/// @param v1 Index of the second corner of the input triangle in `chain`. +/// @param v2 Index of the third corner of the input triangle in `chain`. +/// +/// @param[out] triangulation The (n-2) by 3 array of output triangles, +/// where n is the size of the chain. +/// +template +void split_triangle( + size_t num_points, + span points, + span chain, + span visited_buffer, + std::vector& queue_buffer, + const Index v0, + const Index v1, + const Index v2, + span triangulation); + +} // namespace lagrange diff --git a/modules/core/src/normalize_meshes.cpp b/modules/core/src/normalize_meshes.cpp index 7ebba83..d091aa2 100644 --- a/modules/core/src/normalize_meshes.cpp +++ b/modules/core/src/normalize_meshes.cpp @@ -57,9 +57,9 @@ void normalize_meshes(span*> meshes) } } -#define LA_X_normalize_meshes(_, Scalar, Index) \ - template LA_CORE_API void normalize_mesh(SurfaceMesh & mesh); \ - template LA_CORE_API void normalize_meshes(span*> meshes); +#define LA_X_normalize_meshes(_, Scalar, Index) \ + template LA_CORE_API void normalize_mesh(SurfaceMesh& mesh); \ + template LA_CORE_API void normalize_meshes(span*> meshes); LA_SURFACE_MESH_X(normalize_meshes, 0) } // namespace lagrange diff --git a/modules/core/src/orient_outward.cpp b/modules/core/src/orient_outward.cpp new file mode 100644 index 0000000..86fb399 --- /dev/null +++ b/modules/core/src/orient_outward.cpp @@ -0,0 +1,207 @@ +/* + * Copyright 2024 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ + +#include + +#include +#include +#include +#include + +#include + +namespace lagrange { + +namespace { + +template +bool contains(const std::vector& v, const T& x) +{ + return std::find(v.begin(), v.end(), x) != v.end(); +} + +template +std::vector bfs_orient(const lagrange::SurfaceMesh& mesh) +{ + // Mesh traversal to orient facets consistently. + + // 1. Find oriented patches, i.e. connected components obtained by traversing oriented edges. + std::vector facet_to_patch(mesh.get_num_facets(), 0); + std::vector> adj_patches; + Index num_patches = 0; + { + std::vector marked(mesh.get_num_facets(), false); + std::stack pending; + for (Index f = 0; f < mesh.get_num_facets(); ++f) { + if (marked[f]) { + continue; + } + marked[f] = true; + pending.push(f); + facet_to_patch[f] = num_patches; + adj_patches.push_back({}); + while (!pending.empty()) { + Index fi = pending.top(); + pending.pop(); + const Index pi = num_patches; + for (Index ci = mesh.get_facet_corner_begin(fi); + ci != mesh.get_facet_corner_end(fi); + ++ci) { + const Index vi = mesh.get_corner_vertex(ci); + const Index e = mesh.get_corner_edge(ci); + mesh.foreach_corner_around_edge(e, [&](Index cj) { + const Index fj = mesh.get_corner_facet(cj); + const Index vj = mesh.get_corner_vertex(cj); + if (!marked[fj] && vi != vj) { + // TODO: Assert that vertex(ci + 1) == vj and vertex(cj + i) == vi + la_debug_assert(fi != fj); + marked[fj] = true; + pending.push(fj); + facet_to_patch[fj] = pi; + } + // Adjacent facets with different patch numbers are connected through + // non-oriented edges. Make them adjacent in our "patch graph". + if (marked[fj] && facet_to_patch[fj] != pi) { + const Index pj = facet_to_patch[fj]; + if (!contains(adj_patches[pi], pj)) { + adj_patches[pi].push_back(pj); + adj_patches[pj].push_back(pi); + } + } + }); + } + } + ++num_patches; + } + } + + // 2. Traverse the patch graph and greedily flip adjacent oriented patches. If the mesh is + // orientable, this will result in an consistently oriented mesh. Otherwise (e.g. klein bottle + // or mobius strip), we just get... something. + std::vector should_flip_patch(num_patches, false); + { + std::vector marked(num_patches, false); + std::stack pending; + for (Index p = 0; p < num_patches; ++p) { + if (marked[p]) { + continue; + } + marked[p] = true; + pending.push(p); + while (!pending.empty()) { + Index pi = pending.top(); + pending.pop(); + for (Index pj : adj_patches[pi]) { + if (!marked[pj]) { + marked[pj] = true; + pending.push(pj); + should_flip_patch[pj] = !should_flip_patch[pi]; + } + } + } + } + } + + // 3. Mark facets as needing to be flipped based on patch information + std::vector should_flip(mesh.get_num_facets(), false); + for (Index f = 0; f < mesh.get_num_facets(); ++f) { + should_flip[f] = should_flip_patch[facet_to_patch[f]]; + } + + return should_flip; +} + +} // namespace + +template +void orient_outward(lagrange::SurfaceMesh& mesh, const OrientOptions& options) +{ + la_runtime_assert(mesh.get_dimension() == 3); + + bool had_edges = mesh.has_edges(); + mesh.initialize_edges(); + + // Orient facets consistently within a connected component (if possible) + auto should_flip = bfs_orient(mesh); + + // Compute connected components + ComponentOptions component_options; + component_options.output_attribute_name = "@__component_id__"; + auto num_components = compute_components(mesh, component_options); + auto components = + mesh.template get_attribute(component_options.output_attribute_name).get_all(); + + // Compute signed volumes per component + auto tetra_signed_volume = [](const Eigen::RowVector3d& p1, + const Eigen::RowVector3d& p2, + const Eigen::RowVector3d& p3, + const Eigen::RowVector3d& p4) { + return (p2 - p1).dot((p3 - p1).cross(p4 - p1)) / 6.0; + }; + + std::vector signed_volumes(num_components, 0); + { + auto vertices = vertex_view(mesh).template leftCols<3>().template cast(); + Eigen::RowVector3d zero = Eigen::RowVector3d::Zero(); + for (Index f = 0; f < mesh.get_num_facets(); ++f) { + auto facet = mesh.get_facet_vertices(f); + Index nv = mesh.get_facet_size(f); + Scalar sign = (should_flip[f] ? -1 : 1); + if (nv < 3) { + continue; + } else if (nv == 3) { + signed_volumes[components[f]] -= sign * tetra_signed_volume( + vertices.row(facet[0]), + vertices.row(facet[1]), + vertices.row(facet[2]), + zero); + } else { + Eigen::RowVector3d bary = zero; + for (Index lv = 0; lv < nv; ++lv) { + bary += vertices.row(facet[lv]); + } + bary /= nv; + for (Index lv = 0; lv < nv; ++lv) { + signed_volumes[components[f]] -= sign * tetra_signed_volume( + vertices.row(facet[lv]), + vertices.row(facet[(lv + 1) % nv]), + bary, + zero); + } + } + } + } + + // Reverse orientation of each facet whose component's signed volume doesn't match the desired + // orientation + for (Index f = 0; f < mesh.get_num_facets(); ++f) { + // signbit is true if positive, false if negative + if (std::signbit(signed_volumes[components[f]]) != !options.positive) { + should_flip[f] = !should_flip[f]; + } + } + + // Flip facets in each component with incorrect volume sign + if (!had_edges) { + mesh.clear_edges(); + } + mesh.delete_attribute(component_options.output_attribute_name); + mesh.flip_facets([&](Index f) { return should_flip[f]; }); +} + +#define LA_X_orient_outward(_, Scalar, Index) \ + template LA_CORE_API void orient_outward( \ + SurfaceMesh& mesh, \ + const OrientOptions& options); +LA_SURFACE_MESH_X(orient_outward, 0) + +} // namespace lagrange diff --git a/modules/core/src/orientation.cpp b/modules/core/src/orientation.cpp index 469e075..8f02505 100644 --- a/modules/core/src/orientation.cpp +++ b/modules/core/src/orientation.cpp @@ -99,10 +99,10 @@ AttributeId compute_edge_is_oriented( return id; } -#define LA_X_topology(_, Scalar, Index) \ - template LA_CORE_API bool is_oriented(const SurfaceMesh& mesh); \ - template LA_CORE_API AttributeId compute_edge_is_oriented( \ - SurfaceMesh& mesh, \ +#define LA_X_topology(_, Scalar, Index) \ + template LA_CORE_API bool is_oriented(const SurfaceMesh& mesh); \ + template LA_CORE_API AttributeId compute_edge_is_oriented( \ + SurfaceMesh& mesh, \ const OrientationOptions& options); LA_SURFACE_MESH_X(topology, 0) diff --git a/modules/core/src/select_facets_by_normal_similarity.cpp b/modules/core/src/select_facets_by_normal_similarity.cpp new file mode 100644 index 0000000..a829e92 --- /dev/null +++ b/modules/core/src/select_facets_by_normal_similarity.cpp @@ -0,0 +1,212 @@ +/* + * Copyright 2024 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace lagrange { + +template +AttributeId select_facets_by_normal_similarity( + SurfaceMesh& mesh, + const Index seed_facet_id, + const SelectFacetsByNormalSimilarityOptions& options) +{ + // =================================================== + // Helper typedef and lambdas + // =================================================== + + using Vector3D = typename Eigen::Vector3; + + const Index num_facets = mesh.get_num_facets(); + const Scalar flood_error_limit = static_cast(options.flood_error_limit); + const Scalar flood_second_to_first_order_limit_ratio = + static_cast(options.flood_second_to_first_order_limit_ratio); + + // Check if a face is selectable + std::function is_facet_selectable; + if (options.is_facet_selectable_attribute_name.has_value()) { + AttributeId id = internal::find_attribute( + mesh, + options.is_facet_selectable_attribute_name.value(), + Facet, + AttributeUsage::Scalar, + /* number of channels */ 1); + la_debug_assert(id != invalid_attribute_id()); + lagrange::span attr_view = + mesh.template get_attribute(id).get_all(); + is_facet_selectable = [attr_view, num_facets](const Index facet_id) -> bool { + la_runtime_assert(facet_id < num_facets); + return static_cast(attr_view[facet_id]); + }; + } else { + is_facet_selectable = [](const Index) -> bool { return true; }; + } + + // The energy between two normals + auto normal_direction_error = [](const Vector3D& n1, const Vector3D& n2) -> Scalar { + // We assume normals from compute_facet_normal are normalized + return static_cast(1.0) - std::abs(n1.dot(n2)); + }; + + // =================================================== + // Let's begin + // current code just copied from + // Segmentation/SuperTriangulation/SingleFloodFromSeed() + // =================================================== + + la_runtime_assert( + mesh.is_triangle_mesh(), + "At the moment, select_facets_by_normal_similarity only supports triangle meshes."); + + AttributeId id; + if (mesh.has_attribute(options.facet_normal_attribute_name)) { + id = mesh.get_attribute_id(options.facet_normal_attribute_name); + } else { + FacetNormalOptions fn_options; + fn_options.output_attribute_name = options.facet_normal_attribute_name; + id = compute_facet_normal(mesh, fn_options); + } + auto facet_normals = attribute_matrix_view(mesh, id); + const Vector3D seed_n = facet_normals.row(seed_facet_id); + + // This would be the final value that will be returned + AttributeId selected_facet_id = internal::find_or_create_attribute( + mesh, + options.output_attribute_name, + Facet, + AttributeUsage::Scalar, + /* number of channels */ 1, + internal::ResetToDefault::Yes); + lagrange::span is_facet_selected = + mesh.template ref_attribute(selected_facet_id).ref_all(); + + // DFS search utility + std::vector is_facet_processed(num_facets, false); + std::deque search_queue; + + // + // (0): Add the seed neighbors to the queue to initialize the queue + // + is_facet_processed[seed_facet_id] = true; + is_facet_selected[seed_facet_id] = true; + { + auto process_neighboring_facet = [&](const Index ne_fid) { + const Vector3D ne_n = facet_normals.row(ne_fid); + + if ((!is_facet_processed[ne_fid]) && is_facet_selectable(ne_fid)) { + Scalar error = normal_direction_error(seed_n, ne_n); + if (error < flood_error_limit) { + is_facet_selected[ne_fid] = true; + search_queue.push_back(ne_fid); + } + } + }; + mesh.initialize_edges(); // required to loop over edge neighbors + mesh.foreach_facet_around_facet(seed_facet_id, process_neighboring_facet); + } // end of step 0 + + // + // (1): Run through all neighbors, process them, and push them into a queue + // for further flood + // + while (!search_queue.empty()) { + Index fid = invalid(); // set to invalid before assignment + if (options.search_type == SelectFacetsByNormalSimilarityOptions::SearchType::DFS) { + // depth first: first in, first out + fid = search_queue.front(); + search_queue.pop_front(); + } else { + // breadth first: last in, first out + fid = search_queue.back(); + search_queue.pop_back(); + } + + const Vector3D center_n = + facet_normals.row(fid); // center_n must be similar to seed_n, fid must be selected + auto process_neighboring_facet = [&](const Index ne_fid) { + const Vector3D ne_n = facet_normals.row(ne_fid); + + if ((!is_facet_processed[ne_fid]) && is_facet_selectable(ne_fid)) { + const Scalar error_1 = normal_direction_error(seed_n, ne_n); + const Scalar error_2 = normal_direction_error(center_n, ne_n); + is_facet_processed[ne_fid] = true; + + if ((error_1 < flood_error_limit) && (error_2 < flood_error_limit)) { + // center_n ~ ne_n ~ seed_n + is_facet_selected[ne_fid] = true; + search_queue.push_back(ne_fid); + } else if (error_2 < flood_error_limit * flood_second_to_first_order_limit_ratio) { + // center_n ~ ne_n !~ seed_n + // Second order approximation + is_facet_selected[ne_fid] = true; + search_queue.push_back(ne_fid); + } + } // end of is neighbour selectable + }; // end of loops over neighbours + mesh.foreach_facet_around_facet(fid, process_neighboring_facet); + } // end of DFS/BFS stack + + // + // Smooth the selection by selecting more facets + // We select additional facets if more than half of its neighboring facets agree + // + for (Index st_iter_2 = 0; st_iter_2 < static_cast(options.num_smooth_iterations); + st_iter_2++) { + // Go through boundary faces - check spikes + for (Index i = 0; i < num_facets; i++) { + // We consider adding facet i if it wasn't already selected and is selectable + if (!is_facet_selected[i] && is_facet_selectable(i)) { + const Vector3D self_n = facet_normals.row(i); + Index agreeing_neighbors = 0; + Index num_neighbors = 0; + auto record_neighbor_votes = [&](const Index ne_fid) { + num_neighbors++; + if (is_facet_selected[ne_fid]) { // neighbor can only vote if its selected + const Vector3D ne_n = facet_normals.row(ne_fid); + if (normal_direction_error(self_n, ne_n) < flood_error_limit) { + agreeing_neighbors++; + } + } + }; + mesh.foreach_facet_around_facet(i, record_neighbor_votes); + // an n-polygon can have <= n neighboring facets (could be a boundary facet) + la_debug_assert(num_neighbors <= mesh.get_facet_size(i)); + + if (agreeing_neighbors > num_neighbors / 2) { + is_facet_selected[i] = true; + } + } + } // end of loop over faces + } // end of for(num_smooth_iterations) + + + return selected_facet_id; + +} // select_faces_by_normal_similarity + +#define LA_X_select_facets_by_normal_similarity(_, Scalar, Index) \ + template LA_CORE_API AttributeId select_facets_by_normal_similarity( \ + SurfaceMesh& mesh, \ + const Index seed_facet_id, \ + const SelectFacetsByNormalSimilarityOptions& options); +LA_SURFACE_MESH_X(select_facets_by_normal_similarity, 0) + +} // namespace lagrange diff --git a/modules/core/src/select_facets_in_frustum.cpp b/modules/core/src/select_facets_in_frustum.cpp index 148c469..9ad86bb 100644 --- a/modules/core/src/select_facets_in_frustum.cpp +++ b/modules/core/src/select_facets_in_frustum.cpp @@ -12,9 +12,9 @@ #include -#include -#include #include +#include +#include #include #include @@ -25,7 +25,6 @@ #include - namespace lagrange { template @@ -34,30 +33,30 @@ bool select_facets_in_frustum( const Frustum& frustum, const FrustumSelectionOptions& options) { - using Point3D = typename Eigen::Vector3; - const Point3D n0(frustum.planes[0].normal.data()); - const Point3D n1(frustum.planes[1].normal.data()); - const Point3D n2(frustum.planes[2].normal.data()); - const Point3D n3(frustum.planes[3].normal.data()); - const Point3D p0(frustum.planes[0].point.data()); - const Point3D p1(frustum.planes[1].point.data()); - const Point3D p2(frustum.planes[2].point.data()); - const Point3D p3(frustum.planes[3].point.data()); + using Vector3 = typename Eigen::Vector3; + const Vector3 n0(frustum.planes[0].normal.data()); + const Vector3 n1(frustum.planes[1].normal.data()); + const Vector3 n2(frustum.planes[2].normal.data()); + const Vector3 n3(frustum.planes[3].normal.data()); + const Vector3 p0(frustum.planes[0].point.data()); + const Vector3 p1(frustum.planes[1].point.data()); + const Vector3 p2(frustum.planes[2].point.data()); + const Vector3 p3(frustum.planes[3].point.data()); const Index num_facets = mesh.get_num_facets(); // Per-thread buffers to store intermediate results. struct LocalBuffers { - Point3D v0, v1, v2; // Triangle vertices. - Point3D q0, q1, q2, q3; // Intermediate tet vertices. + Vector3 v0, v1, v2; // Triangle vertices. + Vector3 q0, q1, q2, q3; // Intermediate tet vertices. std::array r0, r1, r2; // Triangle projections. }; tbb::enumerable_thread_specific temp_vars; - auto edge_overlap_with_negative_octant = [](const Point3D& q0, const Point3D& q1) -> bool { + auto edge_overlap_with_negative_octant = [](const Vector3& q0, const Vector3& q1) -> bool { Scalar t_min = 0, t_max = 1; - const Point3D e = {q0[0] - q1[0], q0[1] - q1[1], q0[2] - q1[2]}; + const Vector3 e = {q0[0] - q1[0], q0[1] - q1[1], q0[2] - q1[2]}; if (e[0] > 0) { t_max = std::min(t_max, -q1[0] / e[0]); } else if (e[0] < 0) { @@ -86,8 +85,8 @@ bool select_facets_in_frustum( }; auto compute_plane = - [](const Point3D& q0, const Point3D& q1, const Point3D& q2) -> std::pair { - const Point3D n = { + [](const Vector3& q0, const Vector3& q1, const Vector3& q2) -> std::pair { + const Vector3 n = { q0[2] * (q2[1] - q1[1]) + q1[2] * (q0[1] - q2[1]) + q2[2] * (q1[1] - q0[1]), q0[0] * (q2[2] - q1[2]) + q1[0] * (q0[2] - q2[2]) + q2[0] * (q1[2] - q0[2]), q0[1] * (q2[0] - q1[0]) + q1[1] * (q0[0] - q2[0]) + q2[1] * (q1[0] - q0[0])}; @@ -104,10 +103,10 @@ bool select_facets_in_frustum( return 0; }; - auto triangle_intersects_negative_axis = [&](const Point3D& q0, - const Point3D& q1, - const Point3D& q2, - const Point3D& n, + auto triangle_intersects_negative_axis = [&](const Vector3& q0, + const Vector3& q1, + const Vector3& q2, + const Vector3& n, const Scalar c, const int axis) -> bool { auto& r0 = temp_vars.local().r0; @@ -143,7 +142,7 @@ bool select_facets_in_frustum( }; auto triangle_intersects_negative_axes = - [&](const Point3D& q0, const Point3D& q1, const Point3D& q2) -> bool { + [&](const Vector3& q0, const Vector3& q1, const Vector3& q2) -> bool { const auto r = compute_plane(q0, q1, q2); if (triangle_intersects_negative_axis(q0, q1, q2, r.first, r.second, 0)) return true; if (triangle_intersects_negative_axis(q0, q1, q2, r.first, r.second, 1)) return true; @@ -152,7 +151,7 @@ bool select_facets_in_frustum( }; auto tet_overlap_with_negative_octant = - [&](const Point3D& q0, const Point3D& q1, const Point3D& q2, const Point3D& q3) -> bool { + [&](const Vector3& q0, const Vector3& q1, const Vector3& q2, const Vector3& q3) -> bool { // Check 1: Check if any tet vertices is in negative octant. if ((q0.array() < 0).all()) return true; if ((q1.array() < 0).all()) return true; diff --git a/modules/core/src/thicken_and_close_mesh.cpp b/modules/core/src/thicken_and_close_mesh.cpp index d7c970b..eefbefe 100644 --- a/modules/core/src/thicken_and_close_mesh.cpp +++ b/modules/core/src/thicken_and_close_mesh.cpp @@ -21,6 +21,8 @@ #include #include +#include + namespace lagrange { namespace { diff --git a/modules/core/src/transform_mesh.cpp b/modules/core/src/transform_mesh.cpp index bf0afaa..a08411d 100644 --- a/modules/core/src/transform_mesh.cpp +++ b/modules/core/src/transform_mesh.cpp @@ -105,6 +105,8 @@ void transform_mesh( auto cotransform = compute_cotransform(transform); + bool is_reflection = (transform.linear().determinant() < 0); + par_foreach_named_attribute_read(mesh, [&](auto&& name, auto&& attr_read) { using AttributeType = std::decay_t; using ValueType = typename AttributeType::ValueType; @@ -133,10 +135,11 @@ void transform_mesh( auto set = [&](auto&& Y) { values.template leftCols() = Y.transpose().template cast(); }; + HigherPrecisionType sign(options.reorient && is_reflection ? -1 : 1); switch (attr_read.get_usage()) { case AttributeUsage::Position: set(A * X); break; case AttributeUsage::Normal: - set(coL * X); + set(sign * coL * X); if (options.normalize_normals) { tbb::parallel_for(Index(0), Index(values.rows()), [&](Index c) { values.row(c).template head<3>().stableNormalize(); @@ -145,7 +148,7 @@ void transform_mesh( break; case AttributeUsage::Tangent: [[fallthrough]]; case AttributeUsage::Bitangent: - set(L * X); + set(sign * L * X); if (options.normalize_tangents_bitangents) { tbb::parallel_for(Index(0), Index(values.rows()), [&](Index c) { values.row(c).template head<3>().stableNormalize(); @@ -179,9 +182,10 @@ void transform_mesh( } }); - // TODO: We must flip facets due to transform with negative scale - // if (options.reorient_facets && transform.linear().determinant() < 0) { - // } + // We must flip facets due to transform with negative scale + if (options.reorient && is_reflection) { + mesh.flip_facets([](Index /*f*/) { return true; }); + } } template @@ -194,22 +198,22 @@ SurfaceMesh transformed_mesh( return mesh; } -#define LA_X_transform_mesh(_, Scalar, Index) \ +#define LA_X_transform_mesh(_, Scalar, Index) \ template LA_CORE_API void transform_mesh( \ - SurfaceMesh & mesh, \ - const Eigen::Transform& transform, \ - const TransformOptions& options); \ + SurfaceMesh & mesh, \ + const Eigen::Transform& transform, \ + const TransformOptions& options); \ template LA_CORE_API SurfaceMesh transformed_mesh( \ - SurfaceMesh mesh, \ - const Eigen::Transform& transform, \ - const TransformOptions& options); \ + SurfaceMesh mesh, \ + const Eigen::Transform& transform, \ + const TransformOptions& options); \ template LA_CORE_API void transform_mesh( \ - SurfaceMesh & mesh, \ - const Eigen::Transform& transform, \ - const TransformOptions& options); \ + SurfaceMesh & mesh, \ + const Eigen::Transform& transform, \ + const TransformOptions& options); \ template LA_CORE_API SurfaceMesh transformed_mesh( \ - SurfaceMesh mesh, \ - const Eigen::Transform& transform, \ + SurfaceMesh mesh, \ + const Eigen::Transform& transform, \ const TransformOptions& options); LA_SURFACE_MESH_X(transform_mesh, 0) diff --git a/modules/core/tests/mesh_cleanup/test_split_long_edges.cpp b/modules/core/tests/mesh_cleanup/test_split_long_edges.cpp index 714544b..f365818 100644 --- a/modules/core/tests/mesh_cleanup/test_split_long_edges.cpp +++ b/modules/core/tests/mesh_cleanup/test_split_long_edges.cpp @@ -9,15 +9,208 @@ * OF ANY KIND, either express or implied. See the License for the specific language * governing permissions and limitations under the License. */ -#include +#include #include +#include -#include -#include +#include +#include +#include +#include #include -#include +#include #include +#include +#include +#include + +TEST_CASE("split_long_edges", "[surface][cleanup]") +{ + using namespace lagrange; + using Scalar = double; + using Index = uint32_t; + + auto check_edge_length = + [](auto& mesh, Scalar max_edge_length, std::string_view active_attr_name = "") { + auto is_active = [&](Index fid) { + if (active_attr_name.empty()) return true; + auto active_view = attribute_vector_view(mesh, active_attr_name); + return static_cast(active_view[fid]); + }; + auto attr_id = compute_edge_lengths(mesh); + auto edge_lengths = attribute_vector_view(mesh, attr_id); + for (Index fid = 0; fid < mesh.get_num_facets(); fid++) { + if (!is_active(fid)) continue; + for (Index i = 0; i < 3; i++) { + Index eid = mesh.get_edge(fid, i); + REQUIRE(edge_lengths(eid) <= max_edge_length); + } + } + }; + + auto check_area = [](auto& mesh, Scalar total_area) { + REQUIRE_THAT(compute_mesh_area(mesh), Catch::Matchers::WithinRel(total_area, 1e-6)); + }; + + SECTION("Single triangle") + { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); + mesh.add_vertex({1, 0, 0}); + mesh.add_vertex({0, 1, 0}); + mesh.add_triangle(0, 1, 2); + + SplitLongEdgesOptions options; + options.max_edge_length = 0.5; + options.recursive = true; + split_long_edges(mesh, options); + REQUIRE(mesh.get_num_vertices() == 9); + REQUIRE(mesh.get_num_facets() == 9); + check_edge_length(mesh, options.max_edge_length); + check_area(mesh, 0.5); + } + + SECTION("Two triangles") + { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); + mesh.add_vertex({1, 0, 0}); + mesh.add_vertex({0, 1, 0}); + mesh.add_vertex({1, 1, 0}); + mesh.add_triangle(0, 1, 2); + mesh.add_triangle(2, 1, 3); + + SECTION("non-recursive") + { + SplitLongEdgesOptions options; + options.max_edge_length = 0.5; + options.recursive = false; + split_long_edges(mesh, options); + REQUIRE(mesh.get_num_facets() == 10); + check_area(mesh, 1); + } + + SECTION("recursive") + { + SplitLongEdgesOptions options; + options.max_edge_length = 0.5; + options.recursive = true; + split_long_edges(mesh, options); + REQUIRE(mesh.get_num_facets() == 18); + check_edge_length(mesh, options.max_edge_length); + check_area(mesh, 1); + } + + SECTION("with active region") + { + uint8_t active_buffer[2] = {1, 0}; + mesh.create_attribute( + "active", + AttributeElement::Facet, + AttributeUsage::Scalar, + 1, + active_buffer); + + SplitLongEdgesOptions options; + options.max_edge_length = 0.5; + options.active_region_attribute = "active"; + options.recursive = true; + split_long_edges(mesh, options); + REQUIRE(mesh.get_num_facets() == 12); + REQUIRE(mesh.has_attribute("active")); + check_edge_length(mesh, options.max_edge_length, "active"); + check_area(mesh, 1); + } + + SECTION("with uv") + { + { + Scalar uv[8] = {0, 0, 1, 0, 0, 1, 1, 1}; + Index uv_indices[6] = {0, 1, 2, 2, 1, 3}; + mesh.create_attribute( + "uv", + AttributeElement::Indexed, + AttributeUsage::Vector, + 2, + uv, + uv_indices); + } + SplitLongEdgesOptions options; + options.max_edge_length = 0.5; + options.recursive = true; + split_long_edges(mesh, options); + REQUIRE(mesh.get_num_facets() == 18); + check_edge_length(mesh, options.max_edge_length); + check_area(mesh, 1); + + REQUIRE(mesh.has_attribute("uv")); + + // Ensure uv values are correctly interpolated. + auto& attr = mesh.get_indexed_attribute("uv"); + auto& uv_values = matrix_view(attr.values()); + auto& uv_indices = vector_view(attr.indices()); + auto vertices = vertex_view(mesh); + auto facets = facet_view(mesh); + Index num_facets = mesh.get_num_facets(); + for (Index fid = 0; fid < num_facets; fid++) { + auto v0 = vertices.row(facets(fid, 0)).template segment<2>(0).eval(); + auto v1 = vertices.row(facets(fid, 1)).template segment<2>(0).eval(); + auto v2 = vertices.row(facets(fid, 2)).template segment<2>(0).eval(); + + auto uv_idx_0 = uv_indices(fid * 3); + auto uv_idx_1 = uv_indices(fid * 3 + 1); + auto uv_idx_2 = uv_indices(fid * 3 + 2); + + auto uv0 = uv_values.row(uv_idx_0).template segment<2>(0).eval(); + auto uv1 = uv_values.row(uv_idx_1).template segment<2>(0).eval(); + auto uv2 = uv_values.row(uv_idx_2).template segment<2>(0).eval(); + + REQUIRE_THAT((v0 - uv0).norm(), Catch::Matchers::WithinAbs(0, 1e-6)); + REQUIRE_THAT((v1 - uv1).norm(), Catch::Matchers::WithinAbs(0, 1e-6)); + REQUIRE_THAT((v2 - uv2).norm(), Catch::Matchers::WithinAbs(0, 1e-6)); + } + } + } +} + +TEST_CASE("split_long_edges benchmark", "[surface][mesh_cleanup][!benchmark]") +{ + using namespace lagrange; + using Scalar = double; + using Index = uint32_t; + + auto mesh = lagrange::testing::load_surface_mesh("open/core/ball.obj"); + + SplitLongEdgesOptions options; + options.max_edge_length = 0.1f; + options.recursive = true; + + BENCHMARK_ADVANCED("split_long_edges")(Catch::Benchmark::Chronometer meter) + { + auto mesh_copy = SurfaceMesh::stripped_copy(mesh); + meter.measure([&mesh_copy, options] { split_long_edges(mesh_copy, options); }); + REQUIRE(mesh_copy.get_num_facets() > mesh.get_num_facets()); + }; + +#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS + using MeshType = TriangleMesh3D; + auto legacy_mesh = to_legacy_mesh(mesh); + BENCHMARK("legacy::split_long_edges") + { + auto split_mesh = legacy::split_long_edges( + *legacy_mesh, + options.max_edge_length * options.max_edge_length, + options.recursive); + REQUIRE(split_mesh->get_num_facets() > legacy_mesh->get_num_facets()); + return split_mesh; + }; +#endif +} +#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS + #include + #include + #include TEST_CASE("SplitLongEdgesTest", "[split_long_edges][triangle_mesh][cleanup]" LA_SLOW_DEBUG_FLAG) { using namespace lagrange; @@ -149,3 +342,4 @@ TEST_CASE("SplitLongEdgesTest", "[split_long_edges][triangle_mesh][cleanup]" LA_ REQUIRE(mesh2->is_uv_initialized()); } } +#endif diff --git a/modules/core/tests/test_compute_pointcloud_pca.cpp b/modules/core/tests/test_compute_pointcloud_pca.cpp index fc2b302..7b56680 100644 --- a/modules/core/tests/test_compute_pointcloud_pca.cpp +++ b/modules/core/tests/test_compute_pointcloud_pca.cpp @@ -21,10 +21,139 @@ #include #include -TEST_CASE("ComputePointcloudPCA", "[compute_pointcloud_pca][symmetry]") +TEST_CASE("compute_pointcloud_pca", "[compute_pointcloud_pca][symmetry]") { using namespace lagrange; + using MatrixX3dRowMajor = Eigen::Matrix; + + const double eps = 1e-10; + + // An arbitrary rotation + const Eigen::Matrix3d rotation = + Eigen::AngleAxisd(M_PI * 0.2657, Eigen::Vector3d(-1, 4, -7).normalized()) + .toRotationMatrix(); + + // An arbitrary translation + const Eigen::Vector3d translation(1.34, -5.214, 0.35654); + + + // Some points on the x,y, and z axes + const double a = 0.1; + const double b = 0.4; + const double c = 1.2; + + MatrixX3dRowMajor points(6, 3); + points.row(0) << a, 0, 0; + points.row(1) << -a, 0, 0; + points.row(2) << 0, -b, 0; + points.row(3) << 0, b, 0; + points.row(4) << 0, 0, c; + points.row(5) << 0, 0, -c; + + // Make sure that the pca was correct + auto verify_pca = [eps, a, b, c]( + const double mass, + const Eigen::MatrixXd& pts, + const std::array& eigenvalues, + const std::array, 3>& eigenvectors, + const Eigen::Matrix3d& R, + const Eigen::Vector3d& t) { + const double a2 = a * a; + const double b2 = b * b; + const double c2 = c * c; + + auto approx0 = Catch::Approx(0).margin(eps); + + Eigen::Vector3d weights(eigenvalues.data()); + Eigen::Map> components( + (double*)(eigenvectors.data()), + 3, + 3); + REQUIRE(weights[0] == Catch::Approx(mass * 2 * a2)); + REQUIRE(weights[1] == Catch::Approx(mass * 2 * b2)); + REQUIRE(weights[2] == Catch::Approx(mass * 2 * c2)); + REQUIRE((components.col(0) - R * Eigen::Vector3d(1, 0, 0)).norm() == approx0); + REQUIRE((components.col(1) - R * Eigen::Vector3d(0, 1, 0)).norm() == approx0); + REQUIRE((components.col(2) - R * Eigen::Vector3d(0, 0, 1)).norm() == approx0); + + Eigen::MatrixXd pminustr = pts.rowwise() - t.transpose(); + REQUIRE( + (mass * pminustr.transpose() * pminustr - + components * weights.asDiagonal() * components.transpose()) + .norm() == approx0); + }; + + ComputePointcloudPCAOptions options; + + SECTION("Simple case") + { + MatrixX3dRowMajor points_tr = points; + span points_span(points_tr.data(), points_tr.size()); + options.shift_centroid = false; + options.normalize = false; + auto out = compute_pointcloud_pca(points_span, options); + verify_pca( + 1 /*mass*/, + points_tr, + out.eigenvalues, + out.eigenvectors, + Eigen::Matrix3d::Identity(), + Eigen::Vector3d::Zero()); + } + + SECTION("With rotation") + { + MatrixX3dRowMajor points_tr = points * rotation.transpose(); + span points_span(points_tr.data(), points_tr.size()); + options.shift_centroid = false; + options.normalize = false; + auto out = compute_pointcloud_pca(points_span, options); + verify_pca( + 1 /* mass */, + points_tr, + out.eigenvalues, + out.eigenvectors, + rotation, + Eigen::Vector3d::Zero()); + Eigen::Vector3d center(out.center.data()); + REQUIRE(center.norm() == Catch::Approx(0.).margin(eps)); + } + + SECTION("With rotation and translation") + { + MatrixX3dRowMajor points_tr = + (points * rotation.transpose()).rowwise() + translation.transpose(); + span points_span(points_tr.data(), points_tr.size()); + options.shift_centroid = true; + options.normalize = false; + auto out = compute_pointcloud_pca(points_span, options); + verify_pca(1 /* mass */, points_tr, out.eigenvalues, out.eigenvectors, rotation, translation); + Eigen::Vector3d center(out.center.data()); + REQUIRE((center - translation).norm() == Catch::Approx(0.).margin(eps)); + } + + SECTION("With rotation and translation, also scale the covariance matrix") + { + MatrixX3dRowMajor points_tr = + (points * rotation.transpose()).rowwise() + translation.transpose(); + span points_span(points_tr.data(), points_tr.size()); + const double mass = safe_cast(1.) / (points.rows()); + options.shift_centroid = true; + options.normalize = true; + auto out = compute_pointcloud_pca(points_span, options); + verify_pca(mass, points_tr, out.eigenvalues, out.eigenvectors, rotation, translation); + Eigen::Vector3d center(out.center.data()); + REQUIRE((center - translation).norm() == Catch::Approx(0.).margin(eps)); + } +} + + +#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS + +TEST_CASE("legacy::compute_pointcloud_pca", "[compute_pointcloud_pca][symmetry]") +{ + using namespace lagrange; using MatrixX3dRowMajor = Eigen::Matrix; using MatrixX3dColMajor = Eigen::Matrix; @@ -86,8 +215,10 @@ TEST_CASE("ComputePointcloudPCA", "[compute_pointcloud_pca][symmetry]") SECTION("Simple case") { MatrixType points_tr = points; - auto out = - compute_pointcloud_pca(points_tr, false /*shift_center*/, false /*normalize*/); + auto out = legacy::compute_pointcloud_pca( + points_tr, + false /*shift_center*/, + false /*normalize*/); verify_pca( 1 /*mass*/, points_tr, @@ -100,8 +231,10 @@ TEST_CASE("ComputePointcloudPCA", "[compute_pointcloud_pca][symmetry]") SECTION("With rotation") { MatrixType points_tr = points * rotation.transpose(); - auto out = - compute_pointcloud_pca(points_tr, false /*shift_center*/, false /*normalize*/); + auto out = legacy::compute_pointcloud_pca( + points_tr, + false /*shift_center*/, + false /*normalize*/); verify_pca( 1 /* mass */, points_tr, @@ -116,8 +249,10 @@ TEST_CASE("ComputePointcloudPCA", "[compute_pointcloud_pca][symmetry]") { MatrixType points_tr = (points * rotation.transpose()).rowwise() + translation.transpose(); - auto out = - compute_pointcloud_pca(points_tr, true /*shift_center*/, false /*normalize*/); + auto out = legacy::compute_pointcloud_pca( + points_tr, + true /*shift_center*/, + false /*normalize*/); verify_pca(1 /* mass */, points_tr, out.weights, out.components, rotation, translation); REQUIRE((out.center - translation).norm() == Catch::Approx(0.).margin(eps)); } @@ -127,7 +262,10 @@ TEST_CASE("ComputePointcloudPCA", "[compute_pointcloud_pca][symmetry]") MatrixType points_tr = (points * rotation.transpose()).rowwise() + translation.transpose(); const double mass = safe_cast(1.) / (points.rows()); - auto out = compute_pointcloud_pca(points_tr, true /*shift_center*/, true /*normalize*/); + auto out = legacy::compute_pointcloud_pca( + points_tr, + true /*shift_center*/, + true /*normalize*/); verify_pca(mass, points_tr, out.weights, out.components, rotation, translation); REQUIRE((out.center - translation).norm() == Catch::Approx(0.).margin(eps)); } @@ -139,8 +277,10 @@ TEST_CASE("ComputePointcloudPCA", "[compute_pointcloud_pca][symmetry]") SECTION("Simple case") { MatrixType points_tr = points; - auto out = - compute_pointcloud_pca(points_tr, false /*shift_center*/, false /*normalize*/); + auto out = legacy::compute_pointcloud_pca( + points_tr, + false /*shift_center*/, + false /*normalize*/); verify_pca( 1 /*mass*/, points_tr, @@ -153,8 +293,10 @@ TEST_CASE("ComputePointcloudPCA", "[compute_pointcloud_pca][symmetry]") SECTION("With rotation") { MatrixType points_tr = points * rotation.transpose(); - auto out = - compute_pointcloud_pca(points_tr, false /*shift_center*/, false /*normalize*/); + auto out = legacy::compute_pointcloud_pca( + points_tr, + false /*shift_center*/, + false /*normalize*/); verify_pca( 1 /* mass */, points_tr, @@ -169,8 +311,10 @@ TEST_CASE("ComputePointcloudPCA", "[compute_pointcloud_pca][symmetry]") { MatrixType points_tr = (points * rotation.transpose()).rowwise() + translation.transpose(); - auto out = - compute_pointcloud_pca(points_tr, true /*shift_center*/, false /*normalize*/); + auto out = legacy::compute_pointcloud_pca( + points_tr, + true /*shift_center*/, + false /*normalize*/); verify_pca(1 /* mass */, points_tr, out.weights, out.components, rotation, translation); REQUIRE((out.center - translation).norm() == Catch::Approx(0.).margin(eps)); } @@ -180,10 +324,15 @@ TEST_CASE("ComputePointcloudPCA", "[compute_pointcloud_pca][symmetry]") MatrixType points_tr = (points * rotation.transpose()).rowwise() + translation.transpose(); const double mass = safe_cast(1.) / (points.rows()); - auto out = compute_pointcloud_pca(points_tr, true /*shift_center*/, true /*normalize*/); + auto out = legacy::compute_pointcloud_pca( + points_tr, + true /*shift_center*/, + true /*normalize*/); verify_pca(mass, points_tr, out.weights, out.components, rotation, translation); REQUIRE((out.center - translation).norm() == Catch::Approx(0.).margin(eps)); } } } // end of TEST + +#endif diff --git a/modules/core/tests/test_legacy_select_facets_by_normal_similarity.cpp b/modules/core/tests/test_legacy_select_facets_by_normal_similarity.cpp new file mode 100644 index 0000000..c8afaa6 --- /dev/null +++ b/modules/core/tests/test_legacy_select_facets_by_normal_similarity.cpp @@ -0,0 +1,299 @@ +/* + * Copyright 2024 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ +#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS +#include +#include + +#include +#include + +// +// Compare the selection with some predefined expectations on a cylinder. +// +TEST_CASE("legacy::select_facets_by_normal_similarity", "[select_facets_by_normal_similarity]") +{ + // ========================= + // Set this to true, only if the results need to be visualized. + // ========================= + const bool should_dump_meshes = false; + + + // ========================= + // Helper typedefs + // ========================= + using namespace lagrange; + + using MeshType = TriangleMesh3D; + using Parameters = SelectFacetsByNormalSimilarityParameters; + using Index = MeshType::Index; + using AttributeArray = MeshType::AttributeArray; + using Scalar = MeshType::Scalar; + using VertexType = MeshType::VertexType; + + // ========================= + // Helper lambdas + // ========================= + + // return number of true indices in a bit vector + auto get_bitvector_num_true = [](const std::vector& is_selected) -> Index { + Index ans = 0; + for (Index i = 0; i < static_cast(is_selected.size()); ++i) { + if (is_selected[i]) { + ++ans; + } + } + return ans; + }; + + // convert bitvector to attribute array + auto get_bitvector_as_attribute = [](const std::vector& is_selected) -> AttributeArray { + AttributeArray ans(is_selected.size(), 1); + for (Index i = 0; i < static_cast(is_selected.size()); ++i) { + ans(i) = static_cast(is_selected[i]); + } + return ans; + }; + + // Get midpoint of a facet + auto get_facet_midpoint = [](std::shared_ptr mesh, + const Index facet_id) -> VertexType { + const auto face_vertex_ids = mesh->get_facets().row(facet_id).eval(); + const auto v0 = mesh->get_vertices().row(face_vertex_ids(0)).eval(); + const auto v1 = mesh->get_vertices().row(face_vertex_ids(1)).eval(); + const auto v2 = mesh->get_vertices().row(face_vertex_ids(2)).eval(); + const auto midpoint = ((v0 + v1 + v2) / 3.).eval(); + return midpoint; + }; + + // generate a cylinder -- copied from test_compute_curvature + auto generate_cylinder = [](const Scalar radius, + const Scalar height, + const Index n_radial_segments, + const Index n_vertical_segments) -> std::shared_ptr { + // 2 vertices for closing the top and bottom + const Index n_vertices = n_radial_segments * (n_vertical_segments + 1) + 2; + Vertices3D vertices(n_vertices, 3); + + for (const auto h : range(n_vertical_segments + 1)) { + for (const auto r : range(n_radial_segments)) { + const double angle = 2 * M_PI * r / n_radial_segments; + vertices.row(h * n_radial_segments + r) << radius * cos(angle), radius * sin(angle), + height * h / n_vertical_segments; + } + } + vertices.row(n_vertices - 2) << 0, 0, 0; + vertices.row(n_vertices - 1) << 0, 0, height; + + // set the triangles on the sides + const Index n_top_facets = n_radial_segments; + const Index n_bottom_facets = n_radial_segments; + const Index n_side_facets = 2 * n_radial_segments * n_vertical_segments; + const Index n_facets = n_top_facets + n_bottom_facets + n_side_facets; + + Triangles facets(n_facets, 3); + for (const auto h : range(n_vertical_segments)) { + for (const auto r : range(n_radial_segments)) { + const Index i0 = h * n_radial_segments + r; + const Index i1 = h * n_radial_segments + (r + 1) % n_radial_segments; + const Index j0 = (h + 1) * n_radial_segments + r; + const Index j1 = (h + 1) * n_radial_segments + (r + 1) % n_radial_segments; + // Alternate the diagonal edges + if (r % 2 == 0) { + facets.row(2 * h * n_radial_segments + 2 * r) << i0, i1, j1; + facets.row(2 * h * n_radial_segments + 2 * r + 1) << i0, j1, j0; + } else { + facets.row(2 * h * n_radial_segments + 2 * r) << i0, i1, j0; + facets.row(2 * h * n_radial_segments + 2 * r + 1) << j0, i1, j1; + } + } + } + + // triangles on the top + for (const auto r : range(n_radial_segments)) { + const Index center = n_vertices - 2; + const Index i0 = r; + const Index i1 = (r + 1) % n_radial_segments; + facets.row(n_side_facets + r) << center, i0, i1; + } + + // triangles on the bottom + for (const auto r : range(n_radial_segments)) { + const Index center = n_vertices - 1; + const Index i0 = n_vertical_segments * n_radial_segments + r; + const Index i1 = n_vertical_segments * n_radial_segments + (r + 1) % n_radial_segments; + facets.row(n_side_facets + n_top_facets + r) << center, i0, i1; + } + + return create_mesh(vertices, facets); + }; + + // ========================= + // A simple case, the cylinder mesh + // ========================= + SECTION("cylinder") + { + // Properties of the cylinder + const Index n_radial_segments = 25; + const Index n_vertical_segments = 15; + const Scalar height = 1.5; + const Scalar radius = 2.; + // + const Index n_side_facets = 2 * n_radial_segments * n_vertical_segments; + const Index face_on_bottom_id = n_side_facets + 1; + const Index face_on_side_id = 1; + // + const Index n_vertices = n_radial_segments * (n_vertical_segments + 1) + 2; + const Index vertex_on_bottom_mid_id = n_vertices - 2; + // + auto mesh = generate_cylinder(radius, height, n_radial_segments, n_vertical_segments); + + Parameters parameters; + parameters.flood_error_limit = 0.1; + + SECTION("select-on-bottom") + { + // Set the seed and run the selection + // but prevent one face to be selected + const Index seed_id = face_on_bottom_id; + parameters.should_smooth_boundary = false; + // + const Index dont_select_this = face_on_bottom_id + 1; + // parameters.is_facet_selectable.resize(mesh->get_num_facets(), true); + // parameters.is_facet_selectable[dont_select_this] = false; + std::vector selectables; + for (auto i : range(mesh->get_num_facets())) { + if (i != dont_select_this) { + selectables.push_back(i); + } + } + parameters.set_selectable_facets(*mesh, selectables); + + auto is_facet_selected = + lagrange::select_facets_by_normal_similarity(*mesh, seed_id, parameters); + + // Make sure only and only the bottom is selected + for (const auto facet_id : range(mesh->get_num_facets())) { + const auto midpoint = get_facet_midpoint(mesh, facet_id); + if (facet_id == dont_select_this) { + REQUIRE(is_facet_selected[facet_id] == false); + } else if (midpoint.z() == Catch::Approx(0.)) { + REQUIRE(is_facet_selected[facet_id] == true); + } else { + REQUIRE(is_facet_selected[facet_id] == false); + } + } // end of testing + + // Dump the mesh if needed + if (should_dump_meshes) { + if (!mesh->has_facet_attribute("is_selected")) { + mesh->add_facet_attribute("is_selected"); + } + auto attrib = get_bitvector_as_attribute(is_facet_selected); + mesh->import_facet_attribute("is_selected", attrib); + io::save_mesh("test_select_facets_by_normal_similarity__cylinder_0.vtk", *mesh); + io::save_mesh("test_select_facets_by_normal_similarity__cylinder_0.obj", *mesh); + } + } // end of select on bottom + SECTION("select-on-side") + { + // Set the seed and run the selection + const Index seed_id = face_on_side_id; + parameters.should_smooth_boundary = false; + // + std::vector is_facet_selected = + lagrange::select_facets_by_normal_similarity(*mesh, seed_id, parameters); + + // Make sure only the faces on the row of the seed + // and the adjacent rows are selected + for (const auto facet_id : range(mesh->get_num_facets())) { + const auto facet_midpoint = get_facet_midpoint(mesh, facet_id); + const double dtheta = 2 * M_PI / n_radial_segments; + const Scalar y_min_lim = -sin(dtheta) * radius; + const Scalar y_max_lim = sin(2 * dtheta) * radius; + const Scalar x_min_lim = 0; + const bool is_face_in_the_right_region = (facet_midpoint.y() < y_max_lim) && + (facet_midpoint.y() > y_min_lim) && + (facet_midpoint.x() > x_min_lim); + + if (facet_midpoint.z() == Catch::Approx(0.)) { // don't select bottom + REQUIRE(is_facet_selected[facet_id] == false); + } else if (facet_midpoint.z() == Catch::Approx(height)) { // don't select top + REQUIRE(is_facet_selected[facet_id] == false); + } else if (is_face_in_the_right_region) { // only select right part of the side + // std::cout << facet_id << ";" << facet_midpoint.y() << "; " << y_min_lim << + // ";" << y_max_lim << std::endl; + REQUIRE(is_facet_selected[facet_id] == true); + } else { // don't select other parts + REQUIRE(is_facet_selected[facet_id] == false); + } + } // end of testing + + if (should_dump_meshes) { + if (!mesh->has_facet_attribute("is_selected")) { + mesh->add_facet_attribute("is_selected"); + } + auto attrib = get_bitvector_as_attribute(is_facet_selected); + mesh->import_facet_attribute("is_selected", attrib); + io::save_mesh("test_select_facets_by_normal_similarity__cylinder_1.vtk", *mesh); + io::save_mesh("test_select_facets_by_normal_similarity__cylinder_1.obj", *mesh); + } + } // end of select on side + SECTION("smooth-selection-boundary") + { + // Move the bottom vertex of the mesh, so that some of the bottom triangles + // also get selected. I push it towards the seed. + { + auto vertices = mesh->get_vertices(); + vertices.row(vertex_on_bottom_mid_id) << 1.9, 0.3, -.3; + mesh->initialize(vertices, mesh->get_facets()); + } + + // Set the seed and run the selection + // once with and once without smoothing the selection boundary + const Index seed_id = face_on_side_id; + // + parameters.should_smooth_boundary = false; + std::vector is_facet_selected_no_smooth = + lagrange::select_facets_by_normal_similarity(*mesh, seed_id, parameters); + // + parameters.should_smooth_boundary = true; + std::vector is_facet_selected_with_smooth = + lagrange::select_facets_by_normal_similarity(*mesh, seed_id, parameters); + + + // Smoothing should allow one extra triangle to appear + const Index n_selected_no_smooth = get_bitvector_num_true(is_facet_selected_no_smooth); + const Index n_selected_with_smooth = + get_bitvector_num_true(is_facet_selected_with_smooth); + REQUIRE(n_selected_no_smooth + 1 == n_selected_with_smooth); + + if (should_dump_meshes) { + if (!mesh->has_facet_attribute("is_selected_no_smooth")) { + mesh->add_facet_attribute("is_selected_no_smooth"); + } + auto attrib = get_bitvector_as_attribute(is_facet_selected_no_smooth); + mesh->import_facet_attribute("is_selected_no_smooth", attrib); + // + if (!mesh->has_facet_attribute("is_selected_with_smooth")) { + mesh->add_facet_attribute("is_selected_with_smooth"); + } + attrib = get_bitvector_as_attribute(is_facet_selected_with_smooth); + mesh->import_facet_attribute("is_selected_with_smooth", attrib); + // + io::save_mesh("test_select_facets_by_normal_similarity__cylinder_2.vtk", *mesh); + io::save_mesh("test_select_facets_by_normal_similarity__cylinder_2.obj", *mesh); + } + } // end of smooth-boundary + } // end of cuboid + +} // end of the test section +#endif diff --git a/modules/core/tests/test_legacy_select_facets_in_frustum.cpp b/modules/core/tests/test_legacy_select_facets_in_frustum.cpp new file mode 100644 index 0000000..bd48ab2 --- /dev/null +++ b/modules/core/tests/test_legacy_select_facets_in_frustum.cpp @@ -0,0 +1,247 @@ +/* + * Copyright 2024 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ +#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS +#include +#include + +#include +#include +#include +#include + +template +void run_legacy() +{ + using namespace lagrange; + + SECTION("rectangle") + { + // 2 +-----+ 3 + // |\ | + // | \ | + // | \| + // 0 +-----+ 1 + Vertices3D vertices(4, 3); + vertices << 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0; + Triangles facets(2, 3); + facets << 0, 1, 2, 2, 1, 3; + + auto mesh = create_mesh(vertices, facets); + using MeshType = decltype(mesh)::element_type; + using VertexType = typename MeshType::VertexType; + + auto select_point = [&](Scalar x, Scalar y, Scalar margin = 0.1) { + select_facets_in_frustum( + *mesh, + VertexType(1, 0, 0), + VertexType(x - margin, 0, 0), + VertexType(-1, 0, 0), + VertexType(x + margin, 0, 0), + VertexType(0, 1, 0), + VertexType(0, y - margin, 0), + VertexType(0, -1, 0), + VertexType(0, y + margin, 0)); + }; + + SECTION("select all") + { + select_facets_in_frustum( + *mesh, + VertexType(1, 0, 0), + VertexType(-1, 0, 0), + VertexType(-1, 0, 0), + VertexType(2, 0, 0), + VertexType(0, 1, 0), + VertexType(0, -1, 0), + VertexType(0, -1, 0), + VertexType(0, 2, 0)); + + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) > 0); + REQUIRE(attr(1, 0) > 0); + } + + SECTION("select none") + { + select_facets_in_frustum( + *mesh, + VertexType(1, 0, 0), + VertexType(1.1, 0, 0), + VertexType(-1, 0, 0), + VertexType(2, 0, 0), + VertexType(0, 1, 0), + VertexType(0, -1, 0), + VertexType(0, -1, 0), + VertexType(0, 2, 0)); + + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) == 0); + REQUIRE(attr(1, 0) == 0); + } + + SECTION("select none again") + { + select_facets_in_frustum( + *mesh, + VertexType(1, 0, 0), + VertexType(2, 0, 0), + VertexType(-1, 0, 0), + VertexType(-1, 0, 0), + VertexType(0, 1, 0), + VertexType(0, 2, 0), + VertexType(0, -1, 0), + VertexType(0, -1, 0)); + + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) == 0); + REQUIRE(attr(1, 0) == 0); + } + + SECTION("select none 3") + { + select_facets_in_frustum( + *mesh, + VertexType(1, 0, 0), + VertexType(-1, 0, 0), + VertexType(-1, 0, 0), + VertexType(2, 0, 0), + VertexType(0, 0, 1), + VertexType(0, 0, 0.5), + VertexType(0, 0, -1), + VertexType(0, 0, 1.0)); + + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) == 0); + REQUIRE(attr(1, 0) == 0); + } + + SECTION("select all again") + { + select_facets_in_frustum( + *mesh, + VertexType(1, 0, 0), + VertexType(0.4, 0, 0), + VertexType(-1, 0, 0), + VertexType(0.6, 0, 0), + VertexType(0, 0, 1), + VertexType(0, 0, -0.1), + VertexType(0, 0, -1), + VertexType(0, 0, 0.1)); + + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) > 0); + REQUIRE(attr(1, 0) > 0); + } + + SECTION("select just the origin") + { + select_point(0, 0); + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) > 0); + REQUIRE(attr(1, 0) == 0); + } + + SECTION("select (1, 1)") + { + select_point(1, 1); + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) == 0); + REQUIRE(attr(1, 0) > 0); + } + + SECTION("select (0, 1)") + { + select_point(0, 1); + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) > 0); + REQUIRE(attr(1, 0) > 0); + } + + SECTION("select (1, 0)") + { + select_point(1, 0); + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) > 0); + REQUIRE(attr(1, 0) > 0); + } + + SECTION("select (0.5, 0.5)") + { + select_point(0.5, 0.5); + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) > 0); + REQUIRE(attr(1, 0) > 0); + } + + SECTION("select (0.25, 0.25)") + { + select_point(0.25, 0.25); + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) > 0); + REQUIRE(attr(1, 0) == 0); + } + + SECTION("select (0.75, 0.75)") + { + select_point(0.75, 0.75); + REQUIRE(mesh->has_facet_attribute("is_selected")); + + const auto& attr = mesh->get_facet_attribute("is_selected"); + REQUIRE(attr.rows() == 2); + REQUIRE(attr(0, 0) == 0); + REQUIRE(attr(1, 0) > 0); + } + } +} + +TEST_CASE("legacy::select_facets_in_frustum", "[select_facets_in_frustum]") +{ + run_legacy(); +} +TEST_CASE("legacy::select_facets_in_frustum", "[select_facets_in_frustum]") +{ + run_legacy(); +} +#endif diff --git a/modules/core/tests/test_orient_outward.cpp b/modules/core/tests/test_orient_outward.cpp index 653f990..77097e1 100644 --- a/modules/core/tests/test_orient_outward.cpp +++ b/modules/core/tests/test_orient_outward.cpp @@ -11,14 +11,135 @@ */ #include +#include #include +#include +#include +#include -TEST_CASE("orient_outward", "[Mesh]") +TEST_CASE("orient_outward", "[mesh][orient]") { - auto mesh_in = lagrange::testing::load_mesh("open/core/torus3_in.obj"); - auto mesh_out = lagrange::testing::load_mesh("open/core/torus3_out.obj"); + using Scalar = double; + using Index = uint32_t; + + auto mesh_in = lagrange::testing::load_surface_mesh("open/core/torus3_in.obj"); + auto mesh_out = lagrange::testing::load_surface_mesh("open/core/torus3_out.obj"); + + SECTION("without edges") + { + lagrange::orient_outward(mesh_in); + REQUIRE(vertex_view(mesh_in) == vertex_view(mesh_out)); + REQUIRE(facet_view(mesh_in) == facet_view(mesh_out)); + lagrange::testing::check_mesh(mesh_in); + } + + SECTION("with edges") + { + mesh_in.initialize_edges(); + lagrange::testing::check_mesh(mesh_in); + lagrange::orient_outward(mesh_in); + REQUIRE(vertex_view(mesh_in) == vertex_view(mesh_out)); + REQUIRE(facet_view(mesh_in) == facet_view(mesh_out)); + lagrange::testing::check_mesh(mesh_in); + } +} + +TEST_CASE("orient_outward: poly", "[mesh][orient]") +{ + using Scalar = double; + using Index = uint32_t; + + std::vector names = {"hexaSphere.obj", "noisy-sphere.obj"}; + std::vector is_positive = {true, false}; + + for (size_t i = 0; i < names.size(); ++i) { + for (auto with_edges : {false, true}) { + auto mesh = + lagrange::testing::load_surface_mesh("open/core/poly/" + names[i]); + auto expected = mesh; + if (with_edges) { + mesh.initialize_edges(); + } + + // Should be a no-op + lagrange::OrientOptions opt; + opt.positive = is_positive[i]; + lagrange::orient_outward(mesh, opt); + REQUIRE(vertex_view(mesh) == vertex_view(expected)); + REQUIRE( + vector_view(mesh.get_corner_to_vertex()) == + vector_view(expected.get_corner_to_vertex())); + lagrange::testing::check_mesh(mesh); + + // Randomly flip some facets + mesh.flip_facets([](Index f) { return f % 10 == 0; }); + REQUIRE(vertex_view(mesh) == vertex_view(expected)); + REQUIRE( + vector_view(mesh.get_corner_to_vertex()) != + vector_view(expected.get_corner_to_vertex())); + + // Orient should bring it back to its original state + lagrange::orient_outward(mesh, opt); + REQUIRE(vertex_view(mesh) == vertex_view(expected)); + REQUIRE( + vector_view(mesh.get_corner_to_vertex()) == + vector_view(expected.get_corner_to_vertex())); + lagrange::testing::check_mesh(mesh); + } + } +} + +TEST_CASE("flip_facets: two triangles", "[mesh][flip]") +{ + using Index = uint32_t; + + for (auto with_edges : {false, true}) { + lagrange::SurfaceMesh32f mesh; + mesh.add_vertices(4); + mesh.add_triangle(0, 1, 3); + mesh.add_triangle(2, 3, 1); + if (with_edges) { + mesh.initialize_edges(); + } + + lagrange::testing::check_mesh(mesh); + mesh.flip_facets([](Index f) { return f == 1; }); + lagrange::testing::check_mesh(mesh); + } +} + +TEST_CASE("flip_facets: poly", "[mesh][flip]") +{ + using Scalar = double; + using Index = uint32_t; + + std::vector names = {"poly/tetris.obj", "square.obj"}; + for (auto name : names) { + for (auto with_edges : {false, true}) { + auto mesh = lagrange::testing::load_surface_mesh("open/core/" + name); + if (with_edges) { + mesh.initialize_edges(); + } + + lagrange::testing::check_mesh(mesh); + mesh.flip_facets([](Index f) { return f % 5 == 0; }); + lagrange::testing::check_mesh(mesh); + } + } +} + +#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS + +TEST_CASE("orient_outward: legacy", "[mesh][orient]") +{ + auto mesh_in = + lagrange::testing::load_mesh("open/core/torus3_in.obj"); + auto mesh_out = + lagrange::testing::load_mesh("open/core/torus3_out.obj"); lagrange::orient_outward(*mesh_in); REQUIRE(mesh_in->get_vertices() == mesh_out->get_vertices()); REQUIRE(mesh_in->get_facets() == mesh_out->get_facets()); } + +#endif diff --git a/modules/core/tests/test_select_facets_by_normal_similarity.cpp b/modules/core/tests/test_select_facets_by_normal_similarity.cpp index 79187e2..4f8d40a 100644 --- a/modules/core/tests/test_select_facets_by_normal_similarity.cpp +++ b/modules/core/tests/test_select_facets_by_normal_similarity.cpp @@ -12,19 +12,25 @@ #include #include +#include +#include +#include #include #include +#include +namespace { // // Compare the selection with some predefined expectations on a cylinder. // -TEST_CASE("select_facets_by_normal_similarity", "[select_facets_by_normal_similarity]") +template +void run() { // ========================= // Set this to true, only if the results need to be visualized. // ========================= - const bool should_dump_meshes = false; + // const bool should_dump_meshes = false; // ========================= @@ -32,19 +38,16 @@ TEST_CASE("select_facets_by_normal_similarity", "[select_facets_by_normal_simila // ========================= using namespace lagrange; - using MeshType = TriangleMesh3D; - using Parameters = SelectFacetsByNormalSimilarityParameters; - using Index = MeshType::Index; - using AttributeArray = MeshType::AttributeArray; - using Scalar = MeshType::Scalar; - using VertexType = MeshType::VertexType; + using Index = uint32_t; + using MeshType = SurfaceMesh; + using VertexType = typename Eigen::Vector3; // ========================= // Helper lambdas // ========================= // return number of true indices in a bit vector - auto get_bitvector_num_true = [](const std::vector& is_selected) -> Index { + auto get_bitvector_num_true = [](const lagrange::span is_selected) -> Index { Index ans = 0; for (Index i = 0; i < static_cast(is_selected.size()); ++i) { if (is_selected[i]) { @@ -54,34 +57,29 @@ TEST_CASE("select_facets_by_normal_similarity", "[select_facets_by_normal_simila return ans; }; - // convert bitvector to attribute array - auto get_bitvector_as_attribute = [](const std::vector& is_selected) -> AttributeArray { - AttributeArray ans(is_selected.size(), 1); - for (Index i = 0; i < static_cast(is_selected.size()); ++i) { - ans(i) = static_cast(is_selected[i]); - } - return ans; - }; // Get midpoint of a facet - auto get_facet_midpoint = [](std::shared_ptr mesh, - const Index facet_id) -> VertexType { - const auto face_vertex_ids = mesh->get_facets().row(facet_id).eval(); - const auto v0 = mesh->get_vertices().row(face_vertex_ids(0)).eval(); - const auto v1 = mesh->get_vertices().row(face_vertex_ids(1)).eval(); - const auto v2 = mesh->get_vertices().row(face_vertex_ids(2)).eval(); - const auto midpoint = ((v0 + v1 + v2) / 3.).eval(); - return midpoint; + auto get_facet_midpoint = [](const MeshType& mesh, const Index facet_id) -> VertexType { + Index nvertices = mesh.get_facet_size(facet_id); + VertexType midpoint = VertexType::Zero(); + for (Index i = 0; i < nvertices; i++) { + const auto v = mesh.get_position(mesh.get_facet_vertex(facet_id, i)); + midpoint += VertexType(v.data()); + } + return static_cast(1.0 / double(nvertices)) * midpoint; }; // generate a cylinder -- copied from test_compute_curvature auto generate_cylinder = [](const Scalar radius, const Scalar height, const Index n_radial_segments, - const Index n_vertical_segments) -> std::shared_ptr { + const Index n_vertical_segments) -> MeshType { // 2 vertices for closing the top and bottom const Index n_vertices = n_radial_segments * (n_vertical_segments + 1) + 2; - Vertices3D vertices(n_vertices, 3); + MeshType mesh; + mesh.add_vertices(n_vertices); + + auto vertices = vertex_ref(mesh); // returns an Eigen::Map for (const auto h : range(n_vertical_segments + 1)) { for (const auto r : range(n_radial_segments)) { @@ -94,12 +92,6 @@ TEST_CASE("select_facets_by_normal_similarity", "[select_facets_by_normal_simila vertices.row(n_vertices - 1) << 0, 0, height; // set the triangles on the sides - const Index n_top_facets = n_radial_segments; - const Index n_bottom_facets = n_radial_segments; - const Index n_side_facets = 2 * n_radial_segments * n_vertical_segments; - const Index n_facets = n_top_facets + n_bottom_facets + n_side_facets; - - Triangles facets(n_facets, 3); for (const auto h : range(n_vertical_segments)) { for (const auto r : range(n_radial_segments)) { const Index i0 = h * n_radial_segments + r; @@ -108,11 +100,11 @@ TEST_CASE("select_facets_by_normal_similarity", "[select_facets_by_normal_simila const Index j1 = (h + 1) * n_radial_segments + (r + 1) % n_radial_segments; // Alternate the diagonal edges if (r % 2 == 0) { - facets.row(2 * h * n_radial_segments + 2 * r) << i0, i1, j1; - facets.row(2 * h * n_radial_segments + 2 * r + 1) << i0, j1, j0; + mesh.add_triangle(i0, i1, j1); + mesh.add_triangle(i0, j1, j0); } else { - facets.row(2 * h * n_radial_segments + 2 * r) << i0, i1, j0; - facets.row(2 * h * n_radial_segments + 2 * r + 1) << j0, i1, j1; + mesh.add_triangle(i0, i1, j0); + mesh.add_triangle(j0, i1, j1); } } } @@ -122,7 +114,7 @@ TEST_CASE("select_facets_by_normal_similarity", "[select_facets_by_normal_simila const Index center = n_vertices - 2; const Index i0 = r; const Index i1 = (r + 1) % n_radial_segments; - facets.row(n_side_facets + r) << center, i0, i1; + mesh.add_triangle(center, i0, i1); } // triangles on the bottom @@ -130,10 +122,10 @@ TEST_CASE("select_facets_by_normal_similarity", "[select_facets_by_normal_simila const Index center = n_vertices - 1; const Index i0 = n_vertical_segments * n_radial_segments + r; const Index i1 = n_vertical_segments * n_radial_segments + (r + 1) % n_radial_segments; - facets.row(n_side_facets + n_top_facets + r) << center, i0, i1; + mesh.add_triangle(center, i0, i1); } - return create_mesh(vertices, facets); + return mesh; }; // ========================= @@ -146,75 +138,85 @@ TEST_CASE("select_facets_by_normal_similarity", "[select_facets_by_normal_simila const Index n_vertical_segments = 15; const Scalar height = 1.5; const Scalar radius = 2.; - // + const Index n_side_facets = 2 * n_radial_segments * n_vertical_segments; const Index face_on_bottom_id = n_side_facets + 1; const Index face_on_side_id = 1; - // + const Index n_vertices = n_radial_segments * (n_vertical_segments + 1) + 2; const Index vertex_on_bottom_mid_id = n_vertices - 2; - // - auto mesh = generate_cylinder(radius, height, n_radial_segments, n_vertical_segments); - Parameters parameters; - parameters.flood_error_limit = 0.1; + auto mesh = generate_cylinder(radius, height, n_radial_segments, n_vertical_segments); + SelectFacetsByNormalSimilarityOptions options; + options.flood_error_limit = 0.1; SECTION("select-on-bottom") { + SelectFacetsByNormalSimilarityOptions options_bottom = options; // Set the seed and run the selection // but prevent one face to be selected const Index seed_id = face_on_bottom_id; - parameters.should_smooth_boundary = false; - // + options_bottom.num_smooth_iterations = 0; + options_bottom.is_facet_selectable_attribute_name = "@is_facet_selectable"; + const Index dont_select_this = face_on_bottom_id + 1; - // parameters.is_facet_selectable.resize(mesh->get_num_facets(), true); - // parameters.is_facet_selectable[dont_select_this] = false; - std::vector selectables; - for (auto i : range(mesh->get_num_facets())) { - if (i != dont_select_this) { - selectables.push_back(i); - } - } - parameters.set_selectable_facets(*mesh, selectables); + AttributeId selectability_id = internal::find_or_create_attribute( + mesh, + options_bottom.is_facet_selectable_attribute_name.value(), + Facet, + AttributeUsage::Scalar, + /* number of channels*/ 1, + internal::ResetToDefault::Yes); + lagrange::span attr_ref = + mesh.template ref_attribute(selectability_id).ref_all(); + std::fill(attr_ref.begin(), attr_ref.end(), 1); + attr_ref[dont_select_this] = 0; + - auto is_facet_selected = - lagrange::select_facets_by_normal_similarity(*mesh, seed_id, parameters); + lagrange::select_facets_by_normal_similarity(mesh, seed_id, options_bottom); // Make sure only and only the bottom is selected - for (const auto facet_id : range(mesh->get_num_facets())) { + AttributeId selected_id = internal::find_attribute( + mesh, + options_bottom.output_attribute_name, + Facet, + AttributeUsage::Scalar, + /* num channels */ 1); + const lagrange::span is_facet_selected = + mesh.template ref_attribute(selected_id).ref_all(); + for (const auto facet_id : range(mesh.get_num_facets())) { const auto midpoint = get_facet_midpoint(mesh, facet_id); if (facet_id == dont_select_this) { - REQUIRE(is_facet_selected[facet_id] == false); + REQUIRE(is_facet_selected[facet_id] == 0); } else if (midpoint.z() == Catch::Approx(0.)) { - REQUIRE(is_facet_selected[facet_id] == true); + REQUIRE(is_facet_selected[facet_id] != 0); } else { - REQUIRE(is_facet_selected[facet_id] == false); + REQUIRE(is_facet_selected[facet_id] == 0); } } // end of testing - // Dump the mesh if needed - if (should_dump_meshes) { - if (!mesh->has_facet_attribute("is_selected")) { - mesh->add_facet_attribute("is_selected"); - } - auto attrib = get_bitvector_as_attribute(is_facet_selected); - mesh->import_facet_attribute("is_selected", attrib); - io::save_mesh("test_select_facets_by_normal_similarity__cylinder_0.vtk", *mesh); - io::save_mesh("test_select_facets_by_normal_similarity__cylinder_0.obj", *mesh); - } } // end of select on bottom SECTION("select-on-side") { + SelectFacetsByNormalSimilarityOptions options_side = options; // Set the seed and run the selection const Index seed_id = face_on_side_id; - parameters.should_smooth_boundary = false; - // - std::vector is_facet_selected = - lagrange::select_facets_by_normal_similarity(*mesh, seed_id, parameters); + options_side.num_smooth_iterations = 0; + + lagrange::select_facets_by_normal_similarity(mesh, seed_id, options_side); + + AttributeId selected_id = internal::find_attribute( + mesh, + options_side.output_attribute_name, + Facet, + AttributeUsage::Scalar, + /* num channels */ 1); + const lagrange::span is_facet_selected = + mesh.template ref_attribute(selected_id).ref_all(); // Make sure only the faces on the row of the seed // and the adjacent rows are selected - for (const auto facet_id : range(mesh->get_num_facets())) { + for (const auto facet_id : range(mesh.get_num_facets())) { const auto facet_midpoint = get_facet_midpoint(mesh, facet_id); const double dtheta = 2 * M_PI / n_radial_segments; const Scalar y_min_lim = -sin(dtheta) * radius; @@ -225,74 +227,81 @@ TEST_CASE("select_facets_by_normal_similarity", "[select_facets_by_normal_simila (facet_midpoint.x() > x_min_lim); if (facet_midpoint.z() == Catch::Approx(0.)) { // don't select bottom - REQUIRE(is_facet_selected[facet_id] == false); + REQUIRE(is_facet_selected[facet_id] == 0); } else if (facet_midpoint.z() == Catch::Approx(height)) { // don't select top - REQUIRE(is_facet_selected[facet_id] == false); + REQUIRE(is_facet_selected[facet_id] == 0); } else if (is_face_in_the_right_region) { // only select right part of the side - // std::cout << facet_id << ";" << facet_midpoint.y() << "; " << y_min_lim << - // ";" << y_max_lim << std::endl; - REQUIRE(is_facet_selected[facet_id] == true); + REQUIRE(is_facet_selected[facet_id] != 0); } else { // don't select other parts - REQUIRE(is_facet_selected[facet_id] == false); + REQUIRE(is_facet_selected[facet_id] == 0); } } // end of testing - if (should_dump_meshes) { - if (!mesh->has_facet_attribute("is_selected")) { - mesh->add_facet_attribute("is_selected"); - } - auto attrib = get_bitvector_as_attribute(is_facet_selected); - mesh->import_facet_attribute("is_selected", attrib); - io::save_mesh("test_select_facets_by_normal_similarity__cylinder_1.vtk", *mesh); - io::save_mesh("test_select_facets_by_normal_similarity__cylinder_1.obj", *mesh); - } } // end of select on side SECTION("smooth-selection-boundary") { + SelectFacetsByNormalSimilarityOptions options_nosmoothbdry = options; + SelectFacetsByNormalSimilarityOptions options_smoothbdry = options; // Move the bottom vertex of the mesh, so that some of the bottom triangles // also get selected. I push it towards the seed. { - auto vertices = mesh->get_vertices(); - vertices.row(vertex_on_bottom_mid_id) << 1.9, 0.3, -.3; - mesh->initialize(vertices, mesh->get_facets()); + auto p = mesh.ref_position(vertex_on_bottom_mid_id); + p[0] = 1.9f; + p[1] = 0.3f; + p[2] = -.3f; + FacetNormalOptions fn_options; + fn_options.output_attribute_name = options_smoothbdry.facet_normal_attribute_name; + compute_facet_normal(mesh, fn_options); } // Set the seed and run the selection // once with and once without smoothing the selection boundary const Index seed_id = face_on_side_id; - // - parameters.should_smooth_boundary = false; - std::vector is_facet_selected_no_smooth = - lagrange::select_facets_by_normal_similarity(*mesh, seed_id, parameters); - // - parameters.should_smooth_boundary = true; - std::vector is_facet_selected_with_smooth = - lagrange::select_facets_by_normal_similarity(*mesh, seed_id, parameters); + + options_nosmoothbdry.num_smooth_iterations = 0; + options_nosmoothbdry.output_attribute_name = "@is_selected_nosmoothbdry"; + lagrange::select_facets_by_normal_similarity(mesh, seed_id, options_nosmoothbdry); + + AttributeId selected_id_nosmoothbdry = internal::find_attribute( + mesh, + options_nosmoothbdry.output_attribute_name, + Facet, + AttributeUsage::Scalar, + /* num channels */ 1); + const lagrange::span is_facet_selected_nosmoothbdry = + mesh.template ref_attribute(selected_id_nosmoothbdry).ref_all(); + + options_smoothbdry.num_smooth_iterations = 3; + options_smoothbdry.output_attribute_name = "@is_selected_smoothbdry"; + lagrange::select_facets_by_normal_similarity(mesh, seed_id, options_smoothbdry); + + AttributeId selected_id_smoothbdry = internal::find_attribute( + mesh, + options_smoothbdry.output_attribute_name, + Facet, + AttributeUsage::Scalar, + /* num channels */ 1); + const lagrange::span is_facet_selected_smoothbdry = + mesh.template ref_attribute(selected_id_smoothbdry).ref_all(); // Smoothing should allow one extra triangle to appear - const Index n_selected_no_smooth = get_bitvector_num_true(is_facet_selected_no_smooth); - const Index n_selected_with_smooth = - get_bitvector_num_true(is_facet_selected_with_smooth); - REQUIRE(n_selected_no_smooth + 1 == n_selected_with_smooth); - - if (should_dump_meshes) { - if (!mesh->has_facet_attribute("is_selected_no_smooth")) { - mesh->add_facet_attribute("is_selected_no_smooth"); - } - auto attrib = get_bitvector_as_attribute(is_facet_selected_no_smooth); - mesh->import_facet_attribute("is_selected_no_smooth", attrib); - // - if (!mesh->has_facet_attribute("is_selected_with_smooth")) { - mesh->add_facet_attribute("is_selected_with_smooth"); - } - attrib = get_bitvector_as_attribute(is_facet_selected_with_smooth); - mesh->import_facet_attribute("is_selected_with_smooth", attrib); - // - io::save_mesh("test_select_facets_by_normal_similarity__cylinder_2.vtk", *mesh); - io::save_mesh("test_select_facets_by_normal_similarity__cylinder_2.obj", *mesh); - } + const Index n_selected_nosmoothbdry = + get_bitvector_num_true(is_facet_selected_nosmoothbdry); + const Index n_selected_smoothbdry = + get_bitvector_num_true(is_facet_selected_smoothbdry); + REQUIRE(n_selected_nosmoothbdry + 1 == n_selected_smoothbdry); + } // end of smooth-boundary - } // end of cuboid + } // end of cylinder +} // end of run +} // end of namespace -} // end of the test section +TEST_CASE("select_facets_by_normal_similarity", "[select_facets_by_normal_similarity]") +{ + run(); +} +TEST_CASE("select_facets_by_normal_similarity", "[select_facets_by_normal_similarity]") +{ + run(); +} diff --git a/modules/core/tests/test_select_facets_in_frustum.cpp b/modules/core/tests/test_select_facets_in_frustum.cpp index 2de899e..c9d38b2 100644 --- a/modules/core/tests/test_select_facets_in_frustum.cpp +++ b/modules/core/tests/test_select_facets_in_frustum.cpp @@ -12,237 +12,12 @@ #include #include -#include -#include -#include -#include - #include +#include namespace { -#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS - -template -void run_legacy() -{ - using namespace lagrange; - - SECTION("rectangle") - { - // 2 +-----+ 3 - // |\ | - // | \ | - // | \| - // 0 +-----+ 1 - Vertices3D vertices(4, 3); - vertices << 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0; - Triangles facets(2, 3); - facets << 0, 1, 2, 2, 1, 3; - - auto mesh = create_mesh(vertices, facets); - using MeshType = decltype(mesh)::element_type; - using VertexType = typename MeshType::VertexType; - - auto select_point = [&](Scalar x, Scalar y, Scalar margin = 0.1) { - select_facets_in_frustum( - *mesh, - VertexType(1, 0, 0), - VertexType(x - margin, 0, 0), - VertexType(-1, 0, 0), - VertexType(x + margin, 0, 0), - VertexType(0, 1, 0), - VertexType(0, y - margin, 0), - VertexType(0, -1, 0), - VertexType(0, y + margin, 0)); - }; - - SECTION("select all") - { - select_facets_in_frustum( - *mesh, - VertexType(1, 0, 0), - VertexType(-1, 0, 0), - VertexType(-1, 0, 0), - VertexType(2, 0, 0), - VertexType(0, 1, 0), - VertexType(0, -1, 0), - VertexType(0, -1, 0), - VertexType(0, 2, 0)); - - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) > 0); - REQUIRE(attr(1, 0) > 0); - } - - SECTION("select none") - { - select_facets_in_frustum( - *mesh, - VertexType(1, 0, 0), - VertexType(1.1, 0, 0), - VertexType(-1, 0, 0), - VertexType(2, 0, 0), - VertexType(0, 1, 0), - VertexType(0, -1, 0), - VertexType(0, -1, 0), - VertexType(0, 2, 0)); - - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) == 0); - REQUIRE(attr(1, 0) == 0); - } - - SECTION("select none again") - { - select_facets_in_frustum( - *mesh, - VertexType(1, 0, 0), - VertexType(2, 0, 0), - VertexType(-1, 0, 0), - VertexType(-1, 0, 0), - VertexType(0, 1, 0), - VertexType(0, 2, 0), - VertexType(0, -1, 0), - VertexType(0, -1, 0)); - - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) == 0); - REQUIRE(attr(1, 0) == 0); - } - - SECTION("select none 3") - { - select_facets_in_frustum( - *mesh, - VertexType(1, 0, 0), - VertexType(-1, 0, 0), - VertexType(-1, 0, 0), - VertexType(2, 0, 0), - VertexType(0, 0, 1), - VertexType(0, 0, 0.5), - VertexType(0, 0, -1), - VertexType(0, 0, 1.0)); - - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) == 0); - REQUIRE(attr(1, 0) == 0); - } - - SECTION("select all again") - { - select_facets_in_frustum( - *mesh, - VertexType(1, 0, 0), - VertexType(0.4, 0, 0), - VertexType(-1, 0, 0), - VertexType(0.6, 0, 0), - VertexType(0, 0, 1), - VertexType(0, 0, -0.1), - VertexType(0, 0, -1), - VertexType(0, 0, 0.1)); - - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) > 0); - REQUIRE(attr(1, 0) > 0); - } - - SECTION("select just the origin") - { - select_point(0, 0); - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) > 0); - REQUIRE(attr(1, 0) == 0); - } - - SECTION("select (1, 1)") - { - select_point(1, 1); - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) == 0); - REQUIRE(attr(1, 0) > 0); - } - - SECTION("select (0, 1)") - { - select_point(0, 1); - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) > 0); - REQUIRE(attr(1, 0) > 0); - } - - SECTION("select (1, 0)") - { - select_point(1, 0); - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) > 0); - REQUIRE(attr(1, 0) > 0); - } - - SECTION("select (0.5, 0.5)") - { - select_point(0.5, 0.5); - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) > 0); - REQUIRE(attr(1, 0) > 0); - } - - SECTION("select (0.25, 0.25)") - { - select_point(0.25, 0.25); - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) > 0); - REQUIRE(attr(1, 0) == 0); - } - - SECTION("select (0.75, 0.75)") - { - select_point(0.75, 0.75); - REQUIRE(mesh->has_facet_attribute("is_selected")); - - const auto& attr = mesh->get_facet_attribute("is_selected"); - REQUIRE(attr.rows() == 2); - REQUIRE(attr(0, 0) == 0); - REQUIRE(attr(1, 0) > 0); - } - } -} -#endif - template void run() { @@ -462,18 +237,7 @@ void run() } } } -} // namespace - -#ifdef LAGRANGE_ENABLE_LEGACY_FUNCTIONS -TEST_CASE("legacy::select_facets_in_frustum", "[select_facets_in_frustum]") -{ - run_legacy(); -} -TEST_CASE("legacy::select_facets_in_frustum", "[select_facets_in_frustum]") -{ - run_legacy(); -} -#endif +} // end namespace TEST_CASE("select_facets_in_frustum", "[select_facets_in_frustum]") { @@ -482,4 +246,4 @@ TEST_CASE("select_facets_in_frustum", "[select_facets_in_frustum]") TEST_CASE("select_facets_in_frustum", "[select_facets_in_frustum]") { run(); -} \ No newline at end of file +} diff --git a/modules/core/tests/test_surface_mesh.cpp b/modules/core/tests/test_surface_mesh.cpp index 3e370de..7bda306 100644 --- a/modules/core/tests/test_surface_mesh.cpp +++ b/modules/core/tests/test_surface_mesh.cpp @@ -2876,13 +2876,200 @@ void test_metadata_attribute() lagrange::testing::check_mesh(mesh); } +template +void test_foreach_facet_around_facet() +{ + // STEF + + { + // manifold hybrid mesh + lagrange::SurfaceMesh mesh; + mesh.add_vertices(10); + mesh.add_triangle(0, 1, 2); + mesh.add_polygon({1, 3, 4, 2}); + mesh.add_polygon({0, 2, 6, 8, 5}); + mesh.add_polygon({2, 4, 7, 6}); + mesh.add_triangle(6, 9, 8); + mesh.add_triangle(6, 7, 9); + mesh.initialize_edges(); // This is needed to ensure m_reserved_ids.corner_to_edge != + // invalid_attribute_id() + lagrange::testing::check_mesh(mesh); + + + std::vector visited(mesh.get_num_facets(), false); + auto visit_neighbors = [&](Index f) { visited[f] = true; }; + + std::vector neighbor_of_facet_0 = {false, true, true, false, false, false}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(0, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_0[f]); + } + + std::vector neighbor_of_facet_1 = {true, false, false, true, false, false}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(1, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_1[f]); + } + + std::vector neighbor_of_facet_2 = {true, false, false, true, true, false}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(2, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_2[f]); + } + + std::vector neighbor_of_facet_3 = {false, true, true, false, false, true}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(3, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_3[f]); + } + + std::vector neighbor_of_facet_4 = {false, false, true, false, false, true}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(4, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_4[f]); + } + + std::vector neighbor_of_facet_5 = {false, false, false, true, true, false}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(5, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_5[f]); + } + } + + + { + // manifold hybrid mesh cylinder + lagrange::SurfaceMesh mesh; + mesh.add_vertices(6); + mesh.add_triangle(0, 2, 3); + mesh.add_triangle(2, 4, 3); + mesh.add_polygon({1, 5, 4, 2}); + mesh.add_polygon({0, 3, 5, 1}); + mesh.initialize_edges(); // TODO is this needed? + lagrange::testing::check_mesh(mesh); + + std::vector visited(mesh.get_num_facets(), false); + auto visit_neighbors = [&](Index f) { visited[f] = true; }; + + std::vector neighbor_of_facet_0 = {false, true, false, true}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(0, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_0[f]); + } + + std::vector neighbor_of_facet_1 = {true, false, true, false}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(1, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_1[f]); + } + + std::vector neighbor_of_facet_2 = {false, true, false, true}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(2, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_2[f]); + } + std::vector neighbor_of_facet_3 = {true, false, true, false}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(3, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_3[f]); + } + } + + { + // non-manifold triangular mesh + lagrange::SurfaceMesh mesh; + mesh.add_vertices(6); + mesh.add_triangle(0, 4, 5); + mesh.add_triangle(0, 1, 4); + mesh.add_triangle(1, 2, 4); + mesh.add_triangle(1, 3, 4); + mesh.initialize_edges(); // TODO is this needed? + /* lagrange::testing::check_mesh(mesh); */ + + std::vector visited(mesh.get_num_facets(), false); + auto visit_neighbors = [&](Index f) { visited[f] = true; }; + + std::vector neighbor_of_facet_0 = {false, true, false, false}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(0, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_0[f]); + } + + std::vector neighbor_of_facet_1 = {true, false, true, true}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(1, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_1[f]); + } + + std::vector neighbor_of_facet_2 = {false, true, false, true}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(2, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_2[f]); + } + std::vector neighbor_of_facet_3 = {false, true, true, false}; + visited.assign(mesh.get_num_facets(), false); + mesh.foreach_facet_around_facet(3, visit_neighbors); + for (Index f = 0; f < mesh.get_num_facets(); f++) { + REQUIRE(visited[f] == neighbor_of_facet_3[f]); + } + } + + { + // one triangle + lagrange::SurfaceMesh mesh; + mesh.add_vertices(3); + mesh.add_triangle(0, 1, 2); + mesh.initialize_edges(); // TODO is this needed? + lagrange::testing::check_mesh(mesh); + + bool has_neighbor = false; + auto visit_neighbor = [&has_neighbor](Index) { has_neighbor = true; }; + mesh.foreach_facet_around_facet(0, visit_neighbor); + REQUIRE(!has_neighbor); + } + + { + // two polygons sharing two edges + lagrange::SurfaceMesh mesh; + mesh.add_vertices(6); + mesh.add_polygon({0, 1, 5, 4}); + mesh.add_polygon({1, 2, 3, 4, 5}); + + mesh.initialize_edges(); // TODO is this needed? + lagrange::testing::check_mesh(mesh); + + std::vector visited(mesh.get_num_facets(), 0); + auto visit_neighbor = [&visited](Index f) { visited[f]++; }; + + mesh.foreach_facet_around_facet(0, visit_neighbor); + REQUIRE(visited[0] == 0); + REQUIRE(visited[1] == 2); + + visited.assign(mesh.get_num_facets(), 0); + mesh.foreach_facet_around_facet(1, visit_neighbor); + REQUIRE(visited[0] == 2); + REQUIRE(visited[1] == 0); + } +} + } // namespace -TEST_CASE("SurfaceMesh Construction", "[mesh]") -{ +TEST_CASE("SurfaceMesh Construction", "[mesh]"){ #define LA_X_test_mesh_construction(_, Scalar, Index) test_mesh_construction(); - LA_SURFACE_MESH_X(test_mesh_construction, 0) -} + LA_SURFACE_MESH_X(test_mesh_construction, 0)} TEST_CASE("SurfaceMesh: Remove Elements", "[mesh]") { @@ -2899,20 +3086,15 @@ TEST_CASE("SurfaceMesh: Remove Elements", "[mesh]") } } -TEST_CASE("SurfaceMesh: Storage", "[mesh]") -{ +TEST_CASE("SurfaceMesh: Storage", "[mesh]"){ #define LA_X_test_mesh_storage(_, Scalar, Index) test_mesh_storage(); - LA_SURFACE_MESH_X(test_mesh_storage, 0) -} + LA_SURFACE_MESH_X(test_mesh_storage, 0)} TEST_CASE("SurfaceMesh: Copy and Move", "[mesh]") { - SECTION("Without edges") - { + SECTION("Without edges"){ #define LA_X_test_copy_move_false(_, Scalar, Index) test_copy_move(false); - LA_SURFACE_MESH_X(test_copy_move_false, 0) - } - SECTION("With edges") + LA_SURFACE_MESH_X(test_copy_move_false, 0)} SECTION("With edges") { #define LA_X_test_copy_move_true(_, Scalar, Index) test_copy_move(true); LA_SURFACE_MESH_X(test_copy_move_true, 0) @@ -2984,11 +3166,9 @@ TEST_CASE("SurfaceMesh: Edit Facets With Edges", "[mesh]") LA_SURFACE_MESH_X(test_edit_facets_with_edges, 0) } -TEST_CASE("SurfaceMesh: User Edges", "[mesh]") -{ +TEST_CASE("SurfaceMesh: User Edges", "[mesh]"){ #define LA_X_test_user_edges(_, Scalar, Index) test_user_edges(); - LA_SURFACE_MESH_X(test_user_edges, 0) -} + LA_SURFACE_MESH_X(test_user_edges, 0)} TEST_CASE("SurfaceMesh: Reserved Attributes Basic", "[mesh]") { @@ -3015,11 +3195,9 @@ TEST_CASE("SurfaceMesh: Element Index Type", "[mesh]") LA_ATTRIBUTE_X(test_element_index_type_aux, 0) } -TEST_CASE("SurfaceMesh: Element Index Resize", "[mesh]") -{ +TEST_CASE("SurfaceMesh: Element Index Resize", "[mesh]"){ #define LA_X_test_element_index_resize(_, Scalar, Index) test_element_index_resize(); - LA_SURFACE_MESH_X(test_element_index_resize, 0) -} + LA_SURFACE_MESH_X(test_element_index_resize, 0)} TEST_CASE("SurfaceMesh: Resize Attribute Type", "[mesh]") { @@ -3038,23 +3216,17 @@ TEST_CASE("SurfaceMesh: Copy Attribute", "[mesh]") LA_ATTRIBUTE_X(test_copy_attribute_aux, 0) } -TEST_CASE("SurfaceMesh: Shrink To Fit", "[mesh]") -{ +TEST_CASE("SurfaceMesh: Shrink To Fit", "[mesh]"){ #define LA_X_test_shrink_to_fit(_, Scalar, Index) test_shrink_to_fit(); - LA_SURFACE_MESH_X(test_shrink_to_fit, 0) -} + LA_SURFACE_MESH_X(test_shrink_to_fit, 0)} -TEST_CASE("SurfaceMesh: Compress If Regular", "[mesh]") -{ +TEST_CASE("SurfaceMesh: Compress If Regular", "[mesh]"){ #define LA_X_test_compress_if_regular(_, Scalar, Index) test_compress_if_regular(); - LA_SURFACE_MESH_X(test_compress_if_regular, 0) -} + LA_SURFACE_MESH_X(test_compress_if_regular, 0)} -TEST_CASE("SurfaceMesh: Test facets of size 1 and 2", "[mesh]") -{ +TEST_CASE("SurfaceMesh: Test facets of size 1 and 2", "[mesh]"){ #define LA_X_test_1_and_2_facets(_, Scalar, Index) test_1_and_2_facets(); - LA_SURFACE_MESH_X(test_1_and_2_facets, 0) -} + LA_SURFACE_MESH_X(test_1_and_2_facets, 0)} TEST_CASE("SurfaceMesh: Value Attribute", "[mesh]") { @@ -3065,11 +3237,9 @@ TEST_CASE("SurfaceMesh: Value Attribute", "[mesh]") LA_ATTRIBUTE_X(test_value_attribute_aux, 0) } -TEST_CASE("SurfaceMesh: metadata", "[mesh]") -{ +TEST_CASE("SurfaceMesh: metadata", "[mesh]"){ #define LA_X_test_metadata_attribute(_, Scalar, Index) test_metadata_attribute(); - LA_SURFACE_MESH_X(test_metadata_attribute, 0) -} + LA_SURFACE_MESH_X(test_metadata_attribute, 0)} TEST_CASE("SurfaceMesh: sanity check", "[mesh]") { @@ -3092,3 +3262,10 @@ TEST_CASE("SurfaceMesh: sanity check", "[mesh]") lagrange::testing::check_mesh(mesh); } + +TEST_CASE("SurfaceMesh: foreach_facet_around_facet", "[mesh]") +{ +#define LA_X_test_foreach_facet_around_facet(_, Scalar, Index) \ + test_foreach_facet_around_facet(); + LA_SURFACE_MESH_X(test_foreach_facet_around_facet, 0) +} diff --git a/modules/core/tests/test_transform_mesh.cpp b/modules/core/tests/test_transform_mesh.cpp index e6f8b01..d656f78 100644 --- a/modules/core/tests/test_transform_mesh.cpp +++ b/modules/core/tests/test_transform_mesh.cpp @@ -35,13 +35,16 @@ enum class TestCase : int { Rotation = 3, SymmetryXY = 4, SymmetryXZ = 5, - NumTestCases = 6, + NegativeScaling = 6, + NegativeScalingReorient = 7, + NumTestCases = 8, }; void test_transform_mesh_2d(TestCase test_case) { using Scalar = double; using Index = uint32_t; + using RowVector3i = Eigen::Matrix; lagrange::SurfaceMesh mesh(2); mesh.add_vertex({0, 0}); @@ -67,6 +70,7 @@ void test_transform_mesh_2d(TestCase test_case) auto& uv_attr = mesh.get_indexed_attribute(id_uv); auto vertices = vertex_view(mesh); + auto facets = facet_view(mesh); auto uv = lagrange::matrix_view(uv_attr.values()); for (Index v = 0; v < 3; ++v) { @@ -83,6 +87,7 @@ void test_transform_mesh_2d(TestCase test_case) lagrange::transform_mesh(mesh, Eigen::Affine2d(Eigen::Scaling(Scalar(2)))); REQUIRE(vertices.row(1) == Eigen::RowVector2d(2, 0)); REQUIRE(vertices.row(2) == Eigen::RowVector2d(0, 2)); + REQUIRE(facets.row(0) == RowVector3i(0, 1, 2)); break; case TestCase::NonUniformScaling: lagrange::transform_mesh(mesh, Eigen::Affine2d(Eigen::Scaling(Scalar(2), Scalar(3)))); @@ -108,6 +113,24 @@ void test_transform_mesh_2d(TestCase test_case) REQUIRE(vertices.row(1) == Eigen::RowVector2d(1, 0)); REQUIRE(vertices.row(2) == Eigen::RowVector2d(0, -1)); break; + case TestCase::NegativeScaling: + lagrange::transform_mesh(mesh, Eigen::Affine2d(Eigen::Scaling(Scalar(-1), Scalar(2)))); + REQUIRE(vertices.row(1) == Eigen::RowVector2d(-1, 0)); + REQUIRE(vertices.row(2) == Eigen::RowVector2d(0, 2)); + REQUIRE(facets.row(0) == RowVector3i(0, 1, 2)); + break; + case TestCase::NegativeScalingReorient: { + lagrange::TransformOptions topt; + topt.reorient = true; + lagrange::transform_mesh( + mesh, + Eigen::Affine2d(Eigen::Scaling(Scalar(-1), Scalar(2))), + topt); + REQUIRE(vertices.row(1) == Eigen::RowVector2d(-1, 0)); + REQUIRE(vertices.row(2) == Eigen::RowVector2d(0, 2)); + REQUIRE(facets.row(0) == RowVector3i(2, 1, 0)); + break; + } default: break; } } @@ -117,6 +140,7 @@ void test_transform_mesh_3d(bool pad_with_sign, TestCase test_case) { using Scalar = double; using Index = uint32_t; + using RowVector3i = Eigen::Matrix; lagrange::SurfaceMesh mesh; mesh.add_vertex({0, 0, 0}); @@ -151,6 +175,7 @@ void test_transform_mesh_3d(bool pad_with_sign, TestCase test_case) auto& bitangent_attr = mesh.get_indexed_attribute(id_bitangent); auto vertices = vertex_view(mesh); + auto facets = facet_view(mesh); auto uv = lagrange::matrix_view(uv_attr.values()); auto nrm = lagrange::matrix_view(nrm_attr.values()); auto tangent = lagrange::matrix_view(tangent_attr.values()); @@ -228,6 +253,33 @@ void test_transform_mesh_3d(bool pad_with_sign, TestCase test_case) REQUIRE(tangent.row(0).head<3>() == Eigen::RowVector3d(1, 0, 0)); REQUIRE(bitangent.row(0).head<3>() == Eigen::RowVector3d(0, -1, 0)); break; + + case TestCase::NegativeScaling: + lagrange::transform_mesh( + mesh, + Eigen::Affine3d(Eigen::Scaling(Scalar(-1)))); + REQUIRE(vertices.row(1) == Eigen::RowVector3d(-1, 0, 0)); + REQUIRE(vertices.row(2) == Eigen::RowVector3d(0, -1, 0)); + REQUIRE(nrm.row(0).head<3>() == Eigen::RowVector3d(0, 0, 1)); + REQUIRE(tangent.row(0).head<3>() == Eigen::RowVector3d(-1, 0, 0)); + REQUIRE(bitangent.row(0).head<3>() == Eigen::RowVector3d(0, -1, 0)); + REQUIRE(facets.row(0) == RowVector3i(0, 1, 2)); + break; + case TestCase::NegativeScalingReorient: { + lagrange::TransformOptions topt; + topt.reorient = true; + lagrange::transform_mesh( + mesh, + Eigen::Affine3d(Eigen::Scaling(Scalar(-1))), + topt); + REQUIRE(vertices.row(1) == Eigen::RowVector3d(-1, 0, 0)); + REQUIRE(vertices.row(2) == Eigen::RowVector3d(0, -1, 0)); + REQUIRE(nrm.row(0).head<3>() == Eigen::RowVector3d(0, 0, -1)); + REQUIRE(tangent.row(0).head<3>() == Eigen::RowVector3d(1, 0, 0)); + REQUIRE(bitangent.row(0).head<3>() == Eigen::RowVector3d(0, 1, 0)); + REQUIRE(facets.row(0) == RowVector3i(2, 1, 0)); + break; + } default: break; } } @@ -245,7 +297,8 @@ TEST_CASE("transform_mesh_2d", "[next]") TEST_CASE("transform_mesh_3d", "[next]") { for (int i = 0; i < static_cast(TestCase::NumTestCases); ++i) { - test_transform_mesh_3d(false, static_cast(i)); - test_transform_mesh_3d(true, static_cast(i)); + for (bool pad_with_sign : {false, true}) { + test_transform_mesh_3d(pad_with_sign, static_cast(i)); + } } } diff --git a/modules/io/src/load_fbx.cpp b/modules/io/src/load_fbx.cpp index ce7b3fd..ed09ac5 100644 --- a/modules/io/src/load_fbx.cpp +++ b/modules/io/src/load_fbx.cpp @@ -34,6 +34,7 @@ #include #include +#include #include namespace lagrange::io { @@ -95,7 +96,7 @@ MeshType convert_mesh_ufbx_to_lagrange(const ufbx_mesh* mesh, const LoadOptions& AttributeUsage::Normal); auto& attr = lmesh.template ref_indexed_attribute(id); attr.indices().resize_elements(mesh->vertex_normal.indices.count); - attr.values().resize_elements(mesh->vertex_normal.values.count * dim); + attr.values().resize_elements(mesh->vertex_normal.values.count); std::copy_n( mesh->vertex_normal.indices.begin(), mesh->vertex_normal.indices.count, @@ -119,7 +120,7 @@ MeshType convert_mesh_ufbx_to_lagrange(const ufbx_mesh* mesh, const LoadOptions& uv_dim); auto& attr = lmesh.template ref_indexed_attribute(id); attr.indices().resize_elements(uv_set.vertex_uv.indices.count); - attr.values().resize_elements(uv_set.vertex_uv.values.count * uv_dim); + attr.values().resize_elements(uv_set.vertex_uv.values.count); std::copy_n( uv_set.vertex_uv.indices.begin(), uv_set.vertex_uv.indices.count, @@ -190,7 +191,7 @@ MeshType convert_mesh_ufbx_to_lagrange(const ufbx_mesh* mesh, const LoadOptions& dim); auto& attr = lmesh.template ref_indexed_attribute(id); attr.indices().resize_elements(mesh->vertex_tangent.indices.count); - attr.values().resize_elements(mesh->vertex_tangent.values.count * dim); + attr.values().resize_elements(mesh->vertex_tangent.values.count); std::copy_n( mesh->vertex_tangent.indices.begin(), mesh->vertex_tangent.indices.count, @@ -198,7 +199,7 @@ MeshType convert_mesh_ufbx_to_lagrange(const ufbx_mesh* mesh, const LoadOptions& std::copy_n( &mesh->vertex_tangent.values[0].v[0], mesh->vertex_tangent.values.count * dim, - attr.indices().ref_all().begin()); + attr.values().ref_all().begin()); } // bitangent @@ -210,7 +211,7 @@ MeshType convert_mesh_ufbx_to_lagrange(const ufbx_mesh* mesh, const LoadOptions& dim); auto& attr = lmesh.template ref_indexed_attribute(id); attr.indices().resize_elements(mesh->vertex_bitangent.indices.count); - attr.values().resize_elements(mesh->vertex_bitangent.values.count * dim); + attr.values().resize_elements(mesh->vertex_bitangent.values.count); std::copy_n( mesh->vertex_bitangent.indices.begin(), mesh->vertex_bitangent.indices.count, @@ -218,7 +219,7 @@ MeshType convert_mesh_ufbx_to_lagrange(const ufbx_mesh* mesh, const LoadOptions& std::copy_n( &mesh->vertex_bitangent.values[0].v[0], mesh->vertex_bitangent.values.count * dim, - attr.indices().ref_all().begin()); + attr.values().ref_all().begin()); } // vertex color @@ -230,7 +231,7 @@ MeshType convert_mesh_ufbx_to_lagrange(const ufbx_mesh* mesh, const LoadOptions& color_dim); auto& attr = lmesh.template ref_indexed_attribute(id); attr.indices().resize_elements(mesh->vertex_color.indices.count); - attr.values().resize_elements(mesh->vertex_color.values.count * color_dim); + attr.values().resize_elements(mesh->vertex_color.values.count); std::copy_n( mesh->vertex_color.indices.begin(), mesh->vertex_color.indices.count, @@ -238,7 +239,7 @@ MeshType convert_mesh_ufbx_to_lagrange(const ufbx_mesh* mesh, const LoadOptions& std::copy_n( &mesh->vertex_color.values[0].v[0], mesh->vertex_color.values.count * color_dim, - attr.indices().ref_all().begin()); + attr.values().ref_all().begin()); } if (opt.stitch_vertices) { @@ -299,12 +300,18 @@ SceneType load_simple_scene_fbx(const ufbx_scene* scene, const LoadOptions& opt) SceneType lscene; + lscene.reserve_meshes(static_cast(scene->meshes.count)); for (size_t i = 0; i < scene->meshes.count; ++i) { const ufbx_mesh* mesh = scene->meshes[i]; - auto lmesh = convert_mesh_ufbx_to_lagrange(mesh, opt); - element_index[mesh->element_id] = lscene.add_mesh(lmesh); + element_index[mesh->element_id] = lscene.add_mesh(MeshType()); } + tbb::parallel_for(size_t(0), scene->meshes.count, [&](size_t i) { + const ufbx_mesh* mesh = scene->meshes[i]; + lscene.ref_mesh(static_cast(element_index[mesh->element_id])) = + convert_mesh_ufbx_to_lagrange(mesh, opt); + }); + for (size_t i = 0; i < scene->nodes.count; ++i) { const ufbx_node* node = scene->nodes[i]; if (node->mesh) { diff --git a/modules/io/src/load_mesh_ply.cpp b/modules/io/src/load_mesh_ply.cpp index b3208bc..715e73b 100644 --- a/modules/io/src/load_mesh_ply.cpp +++ b/modules/io/src/load_mesh_ply.cpp @@ -57,6 +57,8 @@ void extract_normal( std::string attr_name = fmt::format("{}_{}{}", internal::to_string(element), internal::to_string(usage), suffix); + logger().debug("Reading normal attribute {} -> {}", name, attr_name); + auto id = mesh.template create_attribute(attr_name, element, usage, 3); auto attr = mesh.template ref_attribute(id).ref_all(); la_runtime_assert(static_cast(attr.size()) == num_entries * 3); @@ -83,6 +85,8 @@ void extract_vertex_uv( std::string attr_name = fmt::format("{}_{}{}", internal::to_string(element), internal::to_string(usage), suffix); + logger().debug("Reading uv attribute {} -> {}", name, attr_name); + auto id = mesh.template create_attribute(attr_name, element, usage, 2); auto attr = mesh.template ref_attribute(id).ref_all(); for (Index i = 0; i < num_vertices; ++i) { @@ -109,6 +113,8 @@ void extract_color( fmt::format("{}_{}{}", internal::to_string(element), internal::to_string(usage), suffix); Index num_channels = has_alpha ? 4 : 3; + logger().debug("Reading color attribute {} -> {}", name, attr_name); + auto id = mesh.template create_attribute(attr_name, element, usage, num_channels); auto attr = mesh.template ref_attribute(id).ref_all(); la_debug_assert(static_cast(attr.size()) == num_entries * num_channels); diff --git a/modules/subdivision/python/src/subdivision.cpp b/modules/subdivision/python/src/subdivision.cpp index 991056d..cd76702 100644 --- a/modules/subdivision/python/src/subdivision.cpp +++ b/modules/subdivision/python/src/subdivision.cpp @@ -136,6 +136,8 @@ void populate_subdivision_module(nb::module_& m) :param mesh: The source mesh. :param num_levels: The number of levels of subdivision to apply. :param scheme: The subdivision scheme to use. +:param adaptive: Whether to use adaptive subdivision. +:param max_edge_length: The maximum edge length for adaptive subdivision. :param vertex_boundary_interpolation: Vertex boundary interpolation rule. :param face_varying_interpolation: Face-varying interpolation rule. :param use_limit_surface: Interpolate all data to the limit surface. diff --git a/modules/testing/include/lagrange/testing/require_approx.h b/modules/testing/include/lagrange/testing/require_approx.h new file mode 100644 index 0000000..7d5d6dd --- /dev/null +++ b/modules/testing/include/lagrange/testing/require_approx.h @@ -0,0 +1,39 @@ +/* + * Copyright 2024 Adobe. All rights reserved. + * This file is licensed to you under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. You may obtain a copy + * of the License at http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software distributed under + * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS + * OF ANY KIND, either express or implied. See the License for the specific language + * governing permissions and limitations under the License. + */ + +#pragma once + +#include + +#include + +namespace lagrange::testing { + +template +void require_approx( + const Eigen::MatrixBase& A, + const Eigen::MatrixBase& B, + typename DerivedA::Scalar eps_rel, + typename DerivedA::Scalar eps_abs) +{ + REQUIRE(A.rows() == B.rows()); + REQUIRE(A.cols() == B.cols()); + for (Eigen::Index i = 0; i < A.size(); ++i) { + REQUIRE_THAT( + A.derived().data()[i], + Catch::Matchers::WithinRel(B.derived().data()[i], eps_rel) || + (Catch::Matchers::WithinAbs(0, eps_abs) && + Catch::Matchers::WithinAbs(B.derived().data()[i], eps_abs))); + } +} + +} // namespace lagrange::testing