diff --git a/CMakeLists.txt b/CMakeLists.txt index 193f26f4..6a1faa68 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -234,12 +234,13 @@ endif() # When using vcpkg, use find_package() to find external libraries if(DEFINED VCPKG_INSTALLED_DIR) set(LAGRANGE_USE_FIND_PACKAGE_DEFAULT ON) - option(LAGRANGE_WITH_EMBREE_3 "Use Embree3 with Lagrange" ON) else() set(LAGRANGE_USE_FIND_PACKAGE_DEFAULT OFF) - option(LAGRANGE_WITH_EMBREE_3 "Use Embree3 with Lagrange" OFF) endif() +# Default to Embree 4 +option(LAGRANGE_WITH_EMBREE_3 "Use Embree3 with Lagrange" OFF) + if(EMSCRIPTEN) set(LAGRANGE_DISABLE_FPE_DEFAULT ON) else() diff --git a/cmake/recipes/external/nanobind.cmake b/cmake/recipes/external/nanobind.cmake index 2df9c5e6..fbb2207e 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.9.2 + GIT_TAG v2.13.0 DOWNLOAD_ONLY ON ) diff --git a/modules/bvh/include/lagrange/bvh/internal/resolve_tjunctions.h b/modules/bvh/include/lagrange/bvh/internal/resolve_tjunctions.h new file mode 100644 index 00000000..b6fd1b21 --- /dev/null +++ b/modules/bvh/include/lagrange/bvh/internal/resolve_tjunctions.h @@ -0,0 +1,39 @@ +/* + * Copyright 2026 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 + +namespace lagrange::bvh::internal { + +/// +/// Resolve T-junctions formed by collinear, overlapping edges (BVH-accelerated variant). +/// +/// This behaves identically to @ref lagrange::bvh::resolve_tjunctions, but locates the candidate +/// vertices near each edge using an AABB tree built over the mesh vertices, rather than an +/// axis-aligned sort-and-sweep. It shares @ref lagrange::bvh::ResolveTJunctionsOptions and produces +/// the same output. Kept as an internal reference/benchmark alternative to the public variant. +/// +/// @param[in,out] mesh Input mesh (triangle or polygonal). Modified in place. +/// @param[in] options Optional settings. +/// +/// @tparam Scalar Mesh scalar type. +/// @tparam Index Mesh index type. +/// +/// @see lagrange::bvh::resolve_tjunctions +/// +template +void resolve_tjunctions(SurfaceMesh& mesh, ResolveTJunctionsOptions options = {}); + +} // namespace lagrange::bvh::internal diff --git a/modules/bvh/include/lagrange/bvh/resolve_tjunctions.h b/modules/bvh/include/lagrange/bvh/resolve_tjunctions.h new file mode 100644 index 00000000..0e48332a --- /dev/null +++ b/modules/bvh/include/lagrange/bvh/resolve_tjunctions.h @@ -0,0 +1,72 @@ +/* + * Copyright 2026 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::bvh { + +/// @addtogroup module-bvh +/// @{ + +/// +/// Option settings for @ref resolve_tjunctions. +/// +struct ResolveTJunctionsOptions +{ + /// Absolute distance tolerance used to detect T-junctions. A vertex whose distance to a + /// non-incident edge is within this tolerance, and whose projection lies strictly between the + /// edge endpoints, is treated as lying on that edge. If negative, the tolerance defaults to + /// `1e-6` times the bounding box diagonal of the mesh. + double tolerance = -1; + + /// If true, only boundary edges (edges adjacent to a single facet) are checked for + /// T-junctions. Set to false to also resolve T-junctions on interior edges. + bool boundary_only = true; + + /// If true (default), triangulate the facets affected by edge splitting, so a triangle-mesh + /// input yields a triangle-mesh output. If false, leave those facets as polygons. + bool triangulate_affected = true; +}; + +/// +/// Resolve T-junctions formed by collinear, overlapping edges. +/// +/// A T-junction occurs when a vertex lies on an edge that it is not topologically connected to, +/// causing the edge to overlap with the (unconnected) sub-edges incident to that vertex. This +/// function splits every such edge at the vertices lying on it, and splits the adjacent facets +/// accordingly, so the output mesh contains no overlapping collinear edges. +/// +/// Vertices are not moved: only edges and facets are subdivided so the topology conforms to the +/// existing vertex positions. The `tolerance` only controls detection, not geometric snapping. +/// +/// Both triangle and polygonal meshes are supported. By default (`triangulate_affected == true`) +/// the facets touched by a split are triangulated, so a triangle-mesh input yields a triangle-mesh +/// output. Set `triangulate_affected` to false to instead keep those facets as polygons, with the +/// split points inserted as additional (collinear) boundary vertices. +/// +/// @param[in,out] mesh Input mesh (triangle or polygonal). Modified in place. +/// @param[in] options Optional settings. +/// +/// @tparam Scalar Mesh scalar type. +/// @tparam Index Mesh index type. +/// +/// @note Consider running @ref remove_duplicate_vertices and @ref remove_degenerate_facets +/// beforehand to avoid propagating pre-existing degeneracies. +/// +template +void resolve_tjunctions(SurfaceMesh& mesh, ResolveTJunctionsOptions options = {}); + +/// @} + +} // namespace lagrange::bvh diff --git a/modules/bvh/python/src/bvh.cpp b/modules/bvh/python/src/bvh.cpp index b4eb1fdf..7ff707d0 100644 --- a/modules/bvh/python/src/bvh.cpp +++ b/modules/bvh/python/src/bvh.cpp @@ -16,6 +16,7 @@ #include #include #include +#include #include #include #include @@ -600,6 +601,36 @@ that do not share vertices are tested (vertex-adjacent facets are skipped). Vertex-adjacent facets are filtered before testing. The geometric test uses include_boundary=false (interior intersection only). )"); + + m.def( + "resolve_tjunctions", + [](MeshType& mesh, double tolerance, bool boundary_only, bool triangulate_affected) { + bvh::ResolveTJunctionsOptions opts; + opts.tolerance = tolerance; + opts.boundary_only = boundary_only; + opts.triangulate_affected = triangulate_affected; + bvh::resolve_tjunctions(mesh, std::move(opts)); + }, + "mesh"_a, + nb::kw_only(), + "tolerance"_a = bvh::ResolveTJunctionsOptions().tolerance, + "boundary_only"_a = bvh::ResolveTJunctionsOptions().boundary_only, + "triangulate_affected"_a = bvh::ResolveTJunctionsOptions().triangulate_affected, + R"(Resolve T-junctions formed by collinear, overlapping edges. + +A T-junction occurs when a vertex lies on an edge it is not topologically connected to. This +splits every such edge at the vertices lying on it, and splits the adjacent facets accordingly, +so the output mesh has no overlapping collinear edges. Vertices are not moved; only edges and +facets are subdivided. Both triangle and polygonal meshes are supported. + +:param mesh: Input mesh, triangle or polygonal (modified in place). +:param tolerance: Absolute distance tolerance for detecting vertices on an edge. If negative, + defaults to 1e-6 times the bounding box diagonal. +:param boundary_only: If True, only boundary edges are checked for T-junctions. Set to False to + also resolve T-junctions on interior edges. +:param triangulate_affected: If True (default), triangulate the facets affected by a split, so a + triangle mesh stays triangular. If False, keep those facets as polygons. +)"); } } // namespace lagrange::python diff --git a/modules/bvh/src/internal/resolve_tjunctions.cpp b/modules/bvh/src/internal/resolve_tjunctions.cpp new file mode 100644 index 00000000..f7d6b894 --- /dev/null +++ b/modules/bvh/src/internal/resolve_tjunctions.cpp @@ -0,0 +1,166 @@ +/* + * Copyright 2026 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 + +// clang-format off +#include +#include +#include +// clang-format on + +#include + +#include +#include +#include + +namespace lagrange::bvh::internal { + +namespace { + +template +void resolve_tjunctions_impl( + SurfaceMesh& mesh, + const ResolveTJunctionsOptions& options) +{ + using Tree = AABB; + using Box = typename Tree::Box; + using Point = typename Tree::Point; + + mesh.initialize_edges(); + const Index num_vertices = mesh.get_num_vertices(); + const Index num_edges = mesh.get_num_edges(); + if (num_vertices == 0 || num_edges == 0) return; + + auto vertices = vertex_view(mesh); + auto to_point = [&](Index v) { + Point p; + for (int d = 0; d < Dim; ++d) p[d] = vertices(v, d); + return p; + }; + + // Resolve the detection tolerance (default is relative to the bounding box diagonal). + Scalar tol = static_cast(options.tolerance); + if (options.tolerance < 0) { + Scalar diag = (vertices.colwise().maxCoeff() - vertices.colwise().minCoeff()).norm(); + tol = static_cast(1e-6) * diag; + } + if (tol < 0) tol = 0; + const Scalar tol_sq = tol * tol; + + // Build an AABB tree over the vertices (each stored as a degenerate point box). + std::vector boxes(static_cast(num_vertices)); + for (Index v = 0; v < num_vertices; ++v) { + const Point p = to_point(v); + boxes[v] = Box(p, p); + } + Tree tree; + tree.build({boxes.data(), boxes.size()}); + + // For each edge, query the tree with the edge's tolerance-expanded box for candidate vertices. + std::vector>> edge_splits(num_edges); + tbb::parallel_for(Index(0), num_edges, [&](Index e) { + if (options.boundary_only && !mesh.is_boundary_edge(e)) return; + auto ev = mesh.get_edge_vertices(e); + const Index v0 = ev[0]; + const Index v1 = ev[1]; + const Point p0 = to_point(v0); + const Point p1 = to_point(v1); + const Point edge_dir = p1 - p0; + const Scalar len_sq = edge_dir.squaredNorm(); + if (len_sq <= 0) return; // degenerate edge + + Box query(p0, p0); + query.extend(p1); + query = + Box((query.min() - Point::Constant(tol)).eval(), + (query.max() + Point::Constant(tol)).eval()); + + tree.intersect( + query, + function_ref([&](typename Tree::Index candidate) { + const Index v = static_cast(candidate); + if (v == v0 || v == v1) return true; + const Point pv = to_point(v); + const Scalar t = (pv - p0).dot(edge_dir) / len_sq; + if (t <= 0 || t >= 1) return true; // must lie strictly between endpoints + const Scalar dist_sq = (pv - (p0 + t * edge_dir)).squaredNorm(); + if (dist_sq > tol_sq) return true; + edge_splits[e].emplace_back(t, v); + return true; + })); + std::sort(edge_splits[e].begin(), edge_splits[e].end(), [](const auto& a, const auto& b) { + return a.first < b.first; + }); + }); + + // Build CSR split lists, ordered from get_edge_vertices(e)[0] to [1] (ascending t). + std::vector edge_split_offsets(num_edges + 1, 0); + std::vector split_pts; + for (Index e = 0; e < num_edges; e++) { + for (const auto& entry : edge_splits[e]) split_pts.push_back(entry.second); + edge_split_offsets[e + 1] = static_cast(split_pts.size()); + } + + if (split_pts.empty()) return; + + // Split edges without retriangulating: each affected facet gets a polygonal copy appended at + // id >= old_num_facets, leaving the originals (to be removed) in place. + const Index old_num_facets = mesh.get_num_facets(); + auto facets_to_remove = lagrange::internal::split_edges_only( + mesh, + function_ref(Index)>([&](Index e) -> span { + const Index n = edge_split_offsets[e + 1] - edge_split_offsets[e]; + return span(split_pts.data() + edge_split_offsets[e], n); + }), + function_ref([](Index) { return true; })); + + // Optionally triangulate only the new facets, then drop the original split facets. + if (options.triangulate_affected) { + auto is_new_facet = [old_num_facets](Index f) { return f >= old_num_facets; }; + triangulate_polygonal_facets(mesh, function_ref(is_new_facet)); + } + mesh.remove_facets(facets_to_remove); +} + +} // namespace + +template +void resolve_tjunctions(SurfaceMesh& mesh, ResolveTJunctionsOptions options) +{ + const Index dim = mesh.get_dimension(); + if (dim == 2) { + resolve_tjunctions_impl(mesh, options); + } else if (dim == 3) { + resolve_tjunctions_impl(mesh, options); + } else { + la_runtime_assert(false, "resolve_tjunctions: only 2D and 3D meshes are supported."); + } +} + +#define LA_X_resolve_tjunctions(_, Scalar, Index) \ + template LA_BVH_API void resolve_tjunctions( \ + SurfaceMesh&, \ + ResolveTJunctionsOptions); +LA_SURFACE_MESH_X(resolve_tjunctions, 0) + +} // namespace lagrange::bvh::internal diff --git a/modules/bvh/src/resolve_tjunctions.cpp b/modules/bvh/src/resolve_tjunctions.cpp new file mode 100644 index 00000000..a105ba98 --- /dev/null +++ b/modules/bvh/src/resolve_tjunctions.cpp @@ -0,0 +1,159 @@ +/* + * Copyright 2026 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 + +// clang-format off +#include +#include +#include +#include +// clang-format on + +#include +#include +#include +#include + +namespace lagrange::bvh { + +template +void resolve_tjunctions(SurfaceMesh& mesh, ResolveTJunctionsOptions options) +{ + mesh.initialize_edges(); + + const Index num_vertices = mesh.get_num_vertices(); + const Index num_edges = mesh.get_num_edges(); + const Index dim = mesh.get_dimension(); + if (num_vertices == 0 || num_edges == 0) return; + + auto vertices = vertex_view(mesh); + auto extent = (vertices.colwise().maxCoeff() - vertices.colwise().minCoeff()).eval(); + + // Resolve the detection tolerance (default is relative to the bounding box diagonal). + Scalar tol = static_cast(options.tolerance); + if (options.tolerance < 0) { + tol = static_cast(1e-6) * extent.norm(); + } + if (tol < 0) tol = 0; + const Scalar tol_sq = tol * tol; + + // Pick the axis with the largest extent to drive the sort-and-sweep candidate search. + Index axis = 0; + { + Scalar best = extent(0); + for (Index d = 1; d < dim; d++) { + if (extent(d) > best) { + best = extent(d); + axis = d; + } + } + } + + // Sort vertex ids by their coordinate along the dominant axis. + std::vector order(num_vertices); + std::iota(order.begin(), order.end(), Index(0)); + tbb::parallel_sort(order.begin(), order.end(), [&](Index a, Index b) { + return vertices(a, axis) < vertices(b, axis); + }); + std::vector sorted_axis(num_vertices); + for (Index i = 0; i < num_vertices; i++) sorted_axis[i] = vertices(order[i], axis); + + // Report each vertex on edge `e` (within tolerance, strictly between endpoints) to `visit`. + auto for_each_split_on_edge = [&](Index e, auto&& visit) { + if (options.boundary_only && !mesh.is_boundary_edge(e)) return; + auto ev = mesh.get_edge_vertices(e); + const Index v0 = ev[0]; + const Index v1 = ev[1]; + auto p0 = vertices.row(v0); + auto p1 = vertices.row(v1); + auto edge_dir = (p1 - p0).eval(); + const Scalar len_sq = edge_dir.squaredNorm(); + if (len_sq <= 0) return; // degenerate edge + + const Scalar lo = std::min(p0(axis), p1(axis)) - tol; + const Scalar hi = std::max(p0(axis), p1(axis)) + tol; + auto it_begin = std::lower_bound(sorted_axis.begin(), sorted_axis.end(), lo); + auto it_end = std::upper_bound(sorted_axis.begin(), sorted_axis.end(), hi); + for (auto it = it_begin; it != it_end; ++it) { + const Index v = order[static_cast(std::distance(sorted_axis.begin(), it))]; + if (v == v0 || v == v1) continue; + auto pv = vertices.row(v); + const Scalar t = (pv - p0).dot(edge_dir) / len_sq; + if (t <= 0 || t >= 1) continue; // must lie strictly between endpoints + const Scalar dist_sq = (pv - (p0 + t * edge_dir)).squaredNorm(); + if (dist_sq > tol_sq) continue; + visit(t, v); + } + }; + + // Pass 1: count split points per edge, then prefix-sum into CSR offsets. + std::vector edge_split_offsets(num_edges + 1, 0); + tbb::parallel_for(Index(0), num_edges, [&](Index e) { + Index count = 0; + for_each_split_on_edge(e, [&](Scalar, Index) { ++count; }); + edge_split_offsets[e + 1] = count; + }); + for (Index e = 0; e < num_edges; e++) edge_split_offsets[e + 1] += edge_split_offsets[e]; + + const Index num_split_pts = edge_split_offsets[num_edges]; + if (num_split_pts == 0) return; + + // Pass 2: fill each edge's range, ordered from get_edge_vertices(e)[0] to [1] (ascending t). + std::vector> split_scratch(num_split_pts); + std::vector split_pts(num_split_pts); + tbb::parallel_for(Index(0), num_edges, [&](Index e) { + const Index begin = edge_split_offsets[e]; + const Index end = edge_split_offsets[e + 1]; + Index cursor = begin; + for_each_split_on_edge(e, [&](Scalar t, Index v) { split_scratch[cursor++] = {t, v}; }); + std::sort( + split_scratch.begin() + begin, + split_scratch.begin() + end, + [](const auto& a, const auto& b) { return a.first < b.first; }); + for (Index i = begin; i < end; i++) split_pts[i] = split_scratch[i].second; + }); + + // Split edges without retriangulating: each affected facet gets a polygonal copy appended at + // id >= old_num_facets, leaving the originals (to be removed) in place. + const Index old_num_facets = mesh.get_num_facets(); + auto facets_to_remove = lagrange::internal::split_edges_only( + mesh, + function_ref(Index)>([&](Index e) -> span { + const Index n = edge_split_offsets[e + 1] - edge_split_offsets[e]; + return span(split_pts.data() + edge_split_offsets[e], n); + }), + function_ref([](Index) { return true; })); + + // Optionally triangulate only the new facets, then drop the original split facets. + if (options.triangulate_affected) { + auto is_new_facet = [old_num_facets](Index f) { return f >= old_num_facets; }; + triangulate_polygonal_facets(mesh, function_ref(is_new_facet)); + } + mesh.remove_facets(facets_to_remove); +} + +#define LA_X_resolve_tjunctions(_, Scalar, Index) \ + template LA_BVH_API void resolve_tjunctions( \ + SurfaceMesh&, \ + ResolveTJunctionsOptions); +LA_SURFACE_MESH_X(resolve_tjunctions, 0) + +} // namespace lagrange::bvh diff --git a/modules/bvh/tests/CMakeLists.txt b/modules/bvh/tests/CMakeLists.txt index 66bfdf19..ae1d4ef5 100644 --- a/modules/bvh/tests/CMakeLists.txt +++ b/modules/bvh/tests/CMakeLists.txt @@ -11,5 +11,5 @@ # lagrange_add_test() -lagrange_include_modules(primitive) -target_link_libraries(test_lagrange_bvh PRIVATE lagrange::primitive) +lagrange_include_modules(primitive subdivision) +target_link_libraries(test_lagrange_bvh PRIVATE lagrange::primitive lagrange::subdivision) diff --git a/modules/bvh/tests/test_resolve_tjunctions.cpp b/modules/bvh/tests/test_resolve_tjunctions.cpp new file mode 100644 index 00000000..0ffdca47 --- /dev/null +++ b/modules/bvh/tests/test_resolve_tjunctions.cpp @@ -0,0 +1,402 @@ +/* + * Copyright 2026 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 +#include +#include +#include + +#include +#include +#include +#include + +namespace { + +// Runs the shared correctness suite against `resolve`, a callable with the signature +// void(SurfaceMesh&, bvh::ResolveTJunctionsOptions). Both the public +// (sort-and-sweep) and internal (AABB) variants must produce identical results. +template +void run_resolve_tjunctions_tests(Resolve&& resolve) +{ + using namespace lagrange; + using Scalar = double; + using Index = uint32_t; + using Options = bvh::ResolveTJunctionsOptions; + + SECTION("basic t-junction 3d") + { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); // 0: a + mesh.add_vertex({2, 0, 0}); // 1: b + mesh.add_vertex({1, 1, 0}); // 2: top + mesh.add_vertex({1, 0, 0}); // 3: c (on edge a-b) + mesh.add_vertex({1, -1, 0}); // 4: bottom + mesh.add_triangle(0, 1, 2); // spans full edge (a, b) + mesh.add_triangle(0, 4, 3); // (a, bottom, c) + mesh.add_triangle(3, 4, 1); // (c, bottom, b) + + const Scalar area_before = compute_mesh_area(mesh); + resolve(mesh, Options{}); + + REQUIRE(mesh.get_num_vertices() == 5); + REQUIRE(mesh.get_num_facets() == 4); + REQUIRE_THAT(compute_mesh_area(mesh), Catch::Matchers::WithinRel(area_before, 1e-12)); + + mesh.initialize_edges(); + REQUIRE(mesh.find_edge_from_vertices(0, 1) == invalid()); // (a, b) is gone + REQUIRE(mesh.find_edge_from_vertices(0, 3) != invalid()); // (a, c) + REQUIRE(mesh.find_edge_from_vertices(3, 1) != invalid()); // (c, b) + } + + SECTION("basic t-junction 2d") + { + SurfaceMesh mesh(2); + mesh.add_vertex({0, 0}); + mesh.add_vertex({2, 0}); + mesh.add_vertex({1, 1}); + mesh.add_vertex({1, 0}); + mesh.add_vertex({1, -1}); + mesh.add_triangle(0, 1, 2); + mesh.add_triangle(0, 4, 3); + mesh.add_triangle(3, 4, 1); + + resolve(mesh, Options{}); + + REQUIRE(mesh.get_num_vertices() == 5); + REQUIRE(mesh.get_num_facets() == 4); + mesh.initialize_edges(); + REQUIRE(mesh.find_edge_from_vertices(0, 1) == invalid()); + } + + SECTION("multiple vertices on one edge") + { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); // 0: a + mesh.add_vertex({3, 0, 0}); // 1: b + mesh.add_vertex({1.5, 1, 0}); // 2: top + mesh.add_vertex({1, 0, 0}); // 3: on edge + mesh.add_vertex({2, 0, 0}); // 4: on edge + mesh.add_triangle(0, 1, 2); + + resolve(mesh, Options{}); + + REQUIRE(mesh.get_num_facets() == 3); // fan of 3 across the two split points + // The split facet (a pentagon before triangulation) must be fully triangulated. + for (Index f = 0; f < mesh.get_num_facets(); ++f) REQUIRE(mesh.get_facet_size(f) == 3); + mesh.initialize_edges(); + REQUIRE(mesh.find_edge_from_vertices(0, 1) == invalid()); + REQUIRE(mesh.find_edge_from_vertices(0, 3) != invalid()); + REQUIRE(mesh.find_edge_from_vertices(3, 4) != invalid()); + REQUIRE(mesh.find_edge_from_vertices(4, 1) != invalid()); + } + + SECTION("triangulate_affected option controls affected facets") + { + // A quad whose bottom edge (0-1) is spanned by a fan of three triangles meeting at + // vertices 4 and 5, which lie on that edge: the quad edge is the T-junction to split. + auto make_mesh = [] { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); // 0 + mesh.add_vertex({3, 0, 0}); // 1 + mesh.add_vertex({3, 1, 0}); // 2 + mesh.add_vertex({0, 1, 0}); // 3 + mesh.add_vertex({1, 0, 0}); // 4: on edge 0-1 + mesh.add_vertex({2, 0, 0}); // 5: on edge 0-1 + mesh.add_vertex({Scalar(1.5), -1, 0}); // 6: apex of the triangle fan + mesh.add_quad(0, 1, 2, 3); + mesh.add_triangle(0, 4, 6); + mesh.add_triangle(4, 5, 6); + mesh.add_triangle(5, 1, 6); + return mesh; + }; + + SECTION("triangulate_affected=false keeps the split facet as a polygon") + { + auto mesh = make_mesh(); + Options options; + options.triangulate_affected = false; + resolve(mesh, options); + + // The quad becomes a hexagon; the three triangles are untouched. + REQUIRE(mesh.get_num_facets() == 4); + REQUIRE_FALSE(mesh.is_triangle_mesh()); + Index num_tris = 0; + Index num_hexagons = 0; + for (Index f = 0; f < mesh.get_num_facets(); ++f) { + if (mesh.get_facet_size(f) == 3) ++num_tris; + if (mesh.get_facet_size(f) == 6) ++num_hexagons; + } + REQUIRE(num_tris == 3); + REQUIRE(num_hexagons == 1); + mesh.initialize_edges(); + REQUIRE(mesh.find_edge_from_vertices(0, 1) == invalid()); + REQUIRE(mesh.find_edge_from_vertices(0, 4) != invalid()); + REQUIRE(mesh.find_edge_from_vertices(4, 5) != invalid()); + REQUIRE(mesh.find_edge_from_vertices(5, 1) != invalid()); + } + + SECTION("triangulate_affected=true triangulates the split facet (the default)") + { + auto mesh = make_mesh(); + Options options; + options.triangulate_affected = true; + resolve(mesh, options); + + // The hexagon is triangulated (into 4 triangles) alongside the 3 original triangles. + for (Index f = 0; f < mesh.get_num_facets(); ++f) REQUIRE(mesh.get_facet_size(f) == 3); + REQUIRE(mesh.get_num_facets() == 7); + mesh.initialize_edges(); + REQUIRE(mesh.find_edge_from_vertices(0, 1) == invalid()); + } + } + + SECTION("no t-junction is a no-op") + { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); + mesh.add_vertex({1, 0, 0}); + mesh.add_vertex({1, 1, 0}); + mesh.add_vertex({0, 1, 0}); + mesh.add_triangle(0, 1, 2); + mesh.add_triangle(0, 2, 3); + + resolve(mesh, Options{}); + + REQUIRE(mesh.get_num_vertices() == 4); + REQUIRE(mesh.get_num_facets() == 2); + } + + SECTION("boundary_only option") + { + // Edge (a, b) is shared by two triangles (interior edge); vertex c lies on it. + auto make_mesh = [] { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); // 0: a + mesh.add_vertex({2, 0, 0}); // 1: b + mesh.add_vertex({1, 1, 0}); // 2: top + mesh.add_vertex({1, -1, 0}); // 3: bottom + mesh.add_vertex({1, 0, 0}); // 4: c on edge (a, b) + mesh.add_triangle(0, 1, 2); // (a, b, top) + mesh.add_triangle(1, 0, 3); // (b, a, bottom) -> edge (a, b) is interior + return mesh; + }; + + SECTION("default only checks boundary edges") + { + auto mesh = make_mesh(); + resolve(mesh, Options{}); + REQUIRE(mesh.get_num_facets() == 2); // interior edge (a, b) left untouched + } + + SECTION("boundary_only=false resolves interior edge") + { + auto mesh = make_mesh(); + Options options; + options.boundary_only = false; + resolve(mesh, options); + REQUIRE(mesh.get_num_facets() == 4); + mesh.initialize_edges(); + REQUIRE(mesh.find_edge_from_vertices(0, 1) == invalid()); + } + } + + SECTION("tolerance controls detection") + { + auto make_mesh = [] { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); + mesh.add_vertex({2, 0, 0}); + mesh.add_vertex({1, 1, 0}); + mesh.add_vertex({1, 1e-3, 0}); // near edge (a, b), offset by 1e-3 + mesh.add_triangle(0, 1, 2); + return mesh; + }; + + SECTION("default tolerance leaves near-vertex untouched") + { + auto mesh = make_mesh(); + resolve(mesh, Options{}); + REQUIRE(mesh.get_num_facets() == 1); + } + + SECTION("loose tolerance splits") + { + auto mesh = make_mesh(); + Options options; + options.tolerance = 1e-2; + resolve(mesh, options); + REQUIRE(mesh.get_num_facets() == 2); + } + } + + SECTION("indexed attribute is preserved") + { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); + mesh.add_vertex({2, 0, 0}); + mesh.add_vertex({1, 1, 0}); + mesh.add_vertex({1, 0, 0}); + mesh.add_vertex({1, -1, 0}); + mesh.add_triangle(0, 1, 2); + mesh.add_triangle(0, 4, 3); + mesh.add_triangle(3, 4, 1); + + // Indexed UV mirroring each corner's (x, y) position. + std::vector uv_values = {0, 0, 2, 0, 1, 1, 1, 0, 1, -1}; + std::vector uv_indices = {0, 1, 2, 0, 4, 3, 3, 4, 1}; + mesh.template create_attribute( + "uv", + AttributeElement::Indexed, + AttributeUsage::UV, + 2, + uv_values, + uv_indices); + + resolve(mesh, Options{}); + + REQUIRE(mesh.get_num_facets() == 4); + REQUIRE(mesh.has_attribute("uv")); + auto& uv_attr = mesh.template get_indexed_attribute("uv"); + auto uv = matrix_view(uv_attr.values()); + // Every UV must equal its corresponding vertex (x, y) since UV mirrors position. + auto positions = vertex_view(mesh); + auto facets = facet_view(mesh); + auto uv_idx = matrix_view(uv_attr.indices()); + for (Index f = 0; f < mesh.get_num_facets(); ++f) { + for (Index k = 0; k < 3; ++k) { + const Index vid = facets(f, k); + const Index cid = f * 3 + k; + const Index uid = uv_idx(cid, 0); + REQUIRE_THAT(uv(uid, 0), Catch::Matchers::WithinAbs(positions(vid, 0), 1e-12)); + REQUIRE_THAT(uv(uid, 1), Catch::Matchers::WithinAbs(positions(vid, 1), 1e-12)); + } + } + } +} + +} // namespace + +TEST_CASE("bvh::resolve_tjunctions", "[bvh][cleanup]") +{ + run_resolve_tjunctions_tests( + [](auto& mesh, auto options) { lagrange::bvh::resolve_tjunctions(mesh, options); }); +} + +TEST_CASE("bvh::internal::resolve_tjunctions", "[bvh][cleanup]") +{ + run_resolve_tjunctions_tests([](auto& mesh, auto options) { + lagrange::bvh::internal::resolve_tjunctions(mesh, options); + }); +} + +TEST_CASE("bvh::resolve_tjunctions benchmark", "[bvh][cleanup][!benchmark]") +{ + using namespace lagrange; + using Scalar = double; + using Index = uint32_t; + + // Build a T-junction-heavy mesh: randomly select a fraction of the dragon's facets, subdivide + // that subset to insert edge-midpoint vertices, then recombine with the untouched remainder. + const auto full = testing::load_surface_mesh("open/core/dragon.obj"); + const Index num_facets = full.get_num_facets(); + + auto make_input = [&](double fraction) { + std::vector shuffled(num_facets); + std::iota(shuffled.begin(), shuffled.end(), Index(0)); + std::shuffle(shuffled.begin(), shuffled.end(), std::mt19937(42)); + + const auto num_selected = static_cast(fraction * num_facets); + std::vector selected(shuffled.begin(), shuffled.begin() + num_selected); + std::vector unselected(shuffled.begin() + num_selected, shuffled.end()); + + auto selected_mesh = extract_submesh(full, {selected.data(), selected.size()}); + const auto unselected_mesh = extract_submesh(full, {unselected.data(), unselected.size()}); + selected_mesh = subdivision::midpoint_subdivision(selected_mesh); + return combine_meshes({&selected_mesh, &unselected_mesh}); + }; + + const auto input_01 = make_input(0.01); + const auto input_10 = make_input(0.10); + const auto input_50 = make_input(0.50); + + // sort-and-sweep is the public variant; AABB is the internal reference alternative. + BENCHMARK_ADVANCED("sort-and-sweep (1%)")(Catch::Benchmark::Chronometer meter) + { + // resolve_tjunctions mutates in place, so give each measured run its own fresh copy. + std::vector> meshes(static_cast(meter.runs()), input_01); + meter.measure([&](int i) { + auto& m = meshes[static_cast(i)]; + bvh::resolve_tjunctions(m); + return m.get_num_facets(); + }); + }; + + BENCHMARK_ADVANCED("AABB (1%)")(Catch::Benchmark::Chronometer meter) + { + std::vector> meshes(static_cast(meter.runs()), input_01); + meter.measure([&](int i) { + auto& m = meshes[static_cast(i)]; + bvh::internal::resolve_tjunctions(m); + return m.get_num_facets(); + }); + }; + + BENCHMARK_ADVANCED("sort-and-sweep (10%)")(Catch::Benchmark::Chronometer meter) + { + std::vector> meshes(static_cast(meter.runs()), input_10); + meter.measure([&](int i) { + auto& m = meshes[static_cast(i)]; + bvh::resolve_tjunctions(m); + return m.get_num_facets(); + }); + }; + + BENCHMARK_ADVANCED("AABB (10%)")(Catch::Benchmark::Chronometer meter) + { + std::vector> meshes(static_cast(meter.runs()), input_10); + meter.measure([&](int i) { + auto& m = meshes[static_cast(i)]; + bvh::internal::resolve_tjunctions(m); + return m.get_num_facets(); + }); + }; + + BENCHMARK_ADVANCED("sort-and-sweep (50%)")(Catch::Benchmark::Chronometer meter) + { + std::vector> meshes(static_cast(meter.runs()), input_50); + meter.measure([&](int i) { + auto& m = meshes[static_cast(i)]; + bvh::resolve_tjunctions(m); + return m.get_num_facets(); + }); + }; + + BENCHMARK_ADVANCED("AABB (50%)")(Catch::Benchmark::Chronometer meter) + { + std::vector> meshes(static_cast(meter.runs()), input_50); + meter.measure([&](int i) { + auto& m = meshes[static_cast(i)]; + bvh::internal::resolve_tjunctions(m); + return m.get_num_facets(); + }); + }; +} diff --git a/modules/core/include/lagrange/internal/attribute_string_utils.h b/modules/core/include/lagrange/internal/attribute_string_utils.h index c73ce36f..1a2d82d6 100644 --- a/modules/core/include/lagrange/internal/attribute_string_utils.h +++ b/modules/core/include/lagrange/internal/attribute_string_utils.h @@ -12,6 +12,7 @@ #pragma once #include +#include #include #include @@ -69,4 +70,13 @@ std::string_view value_type_name(const lagrange::Attribute& attr); template std::string_view value_type_name(); +/// +/// Returns a string representation of the attribute value type. +/// +/// @param[in] value_type Attribute value type. +/// +/// @return String representation. +/// +LA_CORE_API std::string_view value_type_name(AttributeValueType value_type); + } // namespace lagrange::internal diff --git a/modules/core/include/lagrange/internal/split_edges.h b/modules/core/include/lagrange/internal/split_edges.h index 8fe9a715..5034fe42 100644 --- a/modules/core/include/lagrange/internal/split_edges.h +++ b/modules/core/include/lagrange/internal/split_edges.h @@ -17,24 +17,65 @@ namespace lagrange::internal { -/** - * Split edges based on the input split points. - * - * @param mesh Input mesh which will be updated in place. All potential splits should - * be represented as vertices in the mesh. The mesh will be updated in - * place. - * @param get_edge_split_pts A function that returns the split points for each edge. The function - * takes an edge ID as input and returns a span of vertex indices - * representing the split pts. - * @param active_facet A function that takes a facet ID as input and returns true if the facet - * is active. Only active facets will be split. - * - * @return A vector of original facet IDs that were split. - */ +/// +/// Split edges based on the input split points, retriangulating the affected facets accordingly. +/// +/// @param mesh Input mesh which will be updated in place. All potential splits should +/// be represented as vertices in the mesh. The mesh will be updated in +/// place. +/// @param get_edge_split_pts A function that returns the split points for each edge. The function +/// takes an edge ID as input and returns a span of vertex indices +/// representing the split pts. +/// @param active_facet A function that takes a facet ID as input and returns true if the facet +/// is active. Only active facets will be split. +/// +/// @note The original facets that got split will remain in the mesh. The callers are responsible +/// for what to do with them (e.g., remove them if needed). +/// +/// @note Facet, corner and indexed attributes are propagated to the newly created facets/corners. +/// +/// @return A vector of original facet IDs that were split. +/// template std::vector split_edges( SurfaceMesh& mesh, function_ref(Index)> get_edge_split_pts, function_ref active_facet); + +/// +/// Split edges based on the input split points *without* retriangulating the affected facets (in +/// contrast to `split_edges`, which does retriangulate). +/// +/// Each original facet is left untouched in place. For every affected facet, a single new facet is +/// appended whose corner chain has the edge split points inserted into it, so the facet grows by +/// one corner per split point (e.g. a triangle with one split edge becomes a quad). The new facets +/// are appended after all existing facets, so their facet IDs are contiguous and all `>=` the facet +/// count at call time; callers can use this to identify the newly created facets (e.g. to +/// selectively triangulate only them). The mesh may become hybrid as a result. +/// +/// @param mesh Input mesh which will be updated in place. All potential splits should +/// be represented as vertices in the mesh. The mesh will be updated in +/// place. +/// @param get_edge_split_pts A function that returns the split points for each edge. The function +/// takes an edge ID as input and returns a span of vertex indices +/// representing the split pts. +/// @param active_facet A function that takes a facet ID as input and returns true if the facet +/// is active. Only active facets will be affected. +/// +/// @note The original facets that got split remain in the mesh unmodified. The callers are +/// responsible for what to do with them (e.g., remove them if needed). +/// +/// @note Facet, corner and indexed attributes are propagated to the newly created facets/corners. +/// Values at inserted split points are linearly interpolated along the split edge, so the split +/// edge must not be geometrically degenerate (its endpoints must have distinct positions). +/// +/// @return A vector of original facet IDs that were split. +/// +template +std::vector split_edges_only( + SurfaceMesh& mesh, + function_ref(Index)> get_edge_split_pts, + function_ref active_facet); + } // namespace lagrange::internal diff --git a/modules/core/include/lagrange/triangulate_polygonal_facets.h b/modules/core/include/lagrange/triangulate_polygonal_facets.h index 75c93f41..07fd5636 100644 --- a/modules/core/include/lagrange/triangulate_polygonal_facets.h +++ b/modules/core/include/lagrange/triangulate_polygonal_facets.h @@ -12,6 +12,7 @@ #pragma once #include +#include #include @@ -57,6 +58,25 @@ void triangulate_polygonal_facets( SurfaceMesh& mesh, const TriangulationOptions& options = {}); +/// +/// Triangulate polygonal facets of a mesh using a prescribed set of rules. +/// +/// @param[in, out] mesh Polygonal mesh to triangulate in place. +/// @param[in] should_triangulate Predicate determining whether a facet with more than 3 +/// vertices should be triangulated. Facets for which it +/// returns false are left untouched. This applies to both the +/// earcut and centroid-fan schemes. +/// @param[in] options Options for triangulation. +/// +/// @tparam Scalar Mesh scalar type. +/// @tparam Index Mesh index type. +/// +template +void triangulate_polygonal_facets( + SurfaceMesh& mesh, + function_ref should_triangulate, + const TriangulationOptions& options = {}); + /// @} } // namespace lagrange diff --git a/modules/core/python/include/lagrange/python/utils/StubType.h b/modules/core/python/include/lagrange/python/utils/StubType.h index b680b2be..a1a210e5 100644 --- a/modules/core/python/include/lagrange/python/utils/StubType.h +++ b/modules/core/python/include/lagrange/python/utils/StubType.h @@ -1,5 +1,5 @@ /* - * Copyright 2025 Adobe. All rights reserved. + * Copyright 2026 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/python/scripts/meshstat.py b/modules/core/python/scripts/meshstat.py index 7aae0fe7..ba471547 100755 --- a/modules/core/python/scripts/meshstat.py +++ b/modules/core/python/scripts/meshstat.py @@ -17,8 +17,6 @@ :func:`lagrange.scripts.meshstat.main`. """ -from __future__ import annotations - import logging import platform import sys diff --git a/modules/core/python/src/bind_utilities.h b/modules/core/python/src/bind_utilities.h index cbb44e7f..8f7c5b75 100644 --- a/modules/core/python/src/bind_utilities.h +++ b/modules/core/python/src/bind_utilities.h @@ -607,7 +607,9 @@ Vertices listed in `cone_vertices` are considered as cone vertices, which is alw m.def( "triangulate_polygonal_facets", - [](MeshType& mesh, std::string_view scheme) { + [](MeshType& mesh, + std::string_view scheme, + std::optional, Tensor, nb::list>> selected_facets) { lagrange::TriangulationOptions opt; if (scheme == "earcut") { opt.scheme = lagrange::TriangulationOptions::Scheme::Earcut; @@ -616,14 +618,71 @@ Vertices listed in `cone_vertices` are considered as cone vertices, which is alw } else { throw Error(lagrange::format("Unsupported triangulation scheme {}", scheme)); } - lagrange::triangulate_polygonal_facets(mesh, opt); + + if (!selected_facets.has_value()) { + // By default, triangulate every polygonal facet. + lagrange::triangulate_polygonal_facets(mesh, opt); + return; + } + + // Build a per-facet mask from either a list/tensor of facet ids or a boolean mask (a + // length-num_facets tensor whose `True` entries mark facets to triangulate). + const Index num_facets = mesh.get_num_facets(); + std::vector should_triangulate(static_cast(num_facets), 0); + auto mark_ids = [&](span ids) { + for (Index f : ids) { + if (f >= num_facets) { + throw Error( + lagrange::format( + "Facet index {} is out of range (mesh has {} facets)", + f, + num_facets)); + } + should_triangulate[f] = 1; + } + }; + auto& selected = selected_facets.value(); + if (const auto* list_ptr = std::get_if(&selected)) { + auto ids = nb::cast>(*list_ptr); + mark_ids({ids.data(), ids.size()}); + } else if (auto* mask_ptr = std::get_if>(&selected)) { + // Boolean per-facet mask: entry `f` is true iff facet `f` should be triangulated. + if (mask_ptr->ndim() != 1 || + mask_ptr->shape(0) != static_cast(num_facets)) { + throw Error( + lagrange::format( + "Facet mask must be a 1D array of length {} (the number of facets)", + num_facets)); + } + // Access the 1D buffer through a typed view (same pattern as + // `bind_surface_mesh.h`); Tensor<> enforces C-contiguity. + auto mask_view = mask_ptr->template view>(); + for (Index f = 0; f < num_facets; ++f) { + should_triangulate[f] = mask_view(f) ? 1 : 0; + } + } else { + auto [data, shape, stride] = tensor_to_span(std::get>(selected)); + la_runtime_assert(is_dense(shape, stride)); + mark_ids(data); + } + + lagrange::triangulate_polygonal_facets( + mesh, + lagrange::function_ref( + [&](Index f) { return should_triangulate[f] != 0; }), + opt); }, "mesh"_a, "scheme"_a = "earcut", + "selected_facets"_a = nb::none(), R"(Triangulate polygonal facets of the mesh. :param mesh: The input mesh to be triangulated in place. -:param scheme: The triangulation scheme (options are 'earcut' and 'centroid_fan'))"); +:param scheme: The triangulation scheme (options are 'earcut' and 'centroid_fan'). +:param selected_facets: Optional subset of facets to triangulate. Either a list/array of facet ids, + or a boolean per-facet mask (a length ``num_facets`` array whose ``True`` entries mark facets to + triangulate). Honored by both schemes; facets not selected are left untouched. If omitted, all + polygonal facets are triangulated.)"); nb::enum_(m, "ConnectivityType", "Mesh connectivity type") .value( diff --git a/modules/core/python/tests/test_triangulate_polygonal_facets.py b/modules/core/python/tests/test_triangulate_polygonal_facets.py index 61fc3c7d..6811b05c 100644 --- a/modules/core/python/tests/test_triangulate_polygonal_facets.py +++ b/modules/core/python/tests/test_triangulate_polygonal_facets.py @@ -10,6 +10,7 @@ # governing permissions and limitations under the License. # import lagrange +import numpy as np import pytest @@ -47,6 +48,72 @@ def test_cube_with_attribute(self, cube): normal_indices = normal_attr.indices assert normal_indices.num_elements == mesh.num_corners + @pytest.mark.parametrize("scheme", ["earcut", "centroid_fan"]) + def test_selected_facets_bool_mask(self, cube, scheme): + # Both schemes honor `selected_facets`, including for quads. Triangulate only two of the + # six cube (quad) facets and check the other four survive as quads. + mesh = cube + mask = np.zeros(mesh.num_facets, dtype=bool) + mask[0] = True + mask[2] = True + + lagrange.triangulate_polygonal_facets(mesh, scheme, mask) + + sizes = sorted(mesh.get_facet_size(f) for f in range(mesh.num_facets)) + assert sizes.count(4) == 4 # four untouched quads + # A quad becomes 2 triangles (earcut) or 4 triangles via a centroid fan (centroid_fan). + expected_triangles = {"earcut": 4, "centroid_fan": 8}[scheme] + assert sizes.count(3) == expected_triangles + + def test_selected_facets_inputs_are_equivalent(self, cube): + vertices = cube.vertices.copy() + facets = cube.facets.copy() + + def triangulate(selected): + mesh = lagrange.SurfaceMesh() + mesh.vertices = vertices + mesh.facets = facets + lagrange.triangulate_polygonal_facets(mesh, "centroid_fan", selected) + return sorted(mesh.get_facet_size(f) for f in range(mesh.num_facets)) + + mask = np.zeros(cube.num_facets, dtype=bool) + mask[0] = True + mask[2] = True + + from_mask = triangulate(mask) + from_list = triangulate([0, 2]) + from_array = triangulate(np.array([0, 2], dtype=np.uint32)) + + assert from_mask == from_list == from_array + + def test_selected_facets_default_triangulates_all(self, cube): + mesh = cube + lagrange.triangulate_polygonal_facets(mesh, "centroid_fan") + assert all(mesh.get_facet_size(f) == 3 for f in range(mesh.num_facets)) + + @pytest.mark.parametrize("scheme", ["earcut", "centroid_fan"]) + def test_selected_facets_empty_is_noop_with_edges(self, cube, scheme): + # An empty selection must be a no-op for both schemes, even when the mesh has edge + # connectivity (guards against add_polygons asserting on empty facet buffers). + mesh = cube + mesh.initialize_edges() + old_num_facets = mesh.num_facets + empty_mask = np.zeros(mesh.num_facets, dtype=bool) + lagrange.triangulate_polygonal_facets(mesh, scheme, empty_mask) + assert mesh.num_facets == old_num_facets + + def test_selected_facets_bad_mask_size(self, cube): + mesh = cube + with pytest.raises(RuntimeError): + lagrange.triangulate_polygonal_facets( + mesh, "centroid_fan", np.zeros(mesh.num_facets + 1, dtype=bool) + ) + + def test_selected_facets_out_of_range(self, cube): + mesh = cube + with pytest.raises(RuntimeError): + lagrange.triangulate_polygonal_facets(mesh, "centroid_fan", [mesh.num_facets]) + def test_cube_with_attribute_centroid_fan(self, cube): mesh = cube attr_id = lagrange.compute_normal(mesh) diff --git a/modules/core/src/internal/attribute_string_utils.cpp b/modules/core/src/internal/attribute_string_utils.cpp index 20a46608..7a375f26 100644 --- a/modules/core/src/internal/attribute_string_utils.cpp +++ b/modules/core/src/internal/attribute_string_utils.cpp @@ -63,6 +63,17 @@ std::string_view to_string(AttributeUsage usage) } } +std::string_view value_type_name(AttributeValueType value_type) +{ + switch (value_type) { +#define LA_X_value_type_name(_, ValueType) \ + case AttributeValueType::e_##ValueType: return #ValueType; + LA_ATTRIBUTE_X(value_type_name, 0) +#undef LA_X_value_type_name + default: la_debug_assert(false, "Unsupported enum type"); return ""; + } +} + #define LA_X_type_name(_, ValueType) \ template <> \ LA_CORE_API std::string_view value_type_name(const lagrange::Attribute&) \ diff --git a/modules/core/src/internal/split_edges.cpp b/modules/core/src/internal/split_edges.cpp index e176b5de..ea0e6614 100644 --- a/modules/core/src/internal/split_edges.cpp +++ b/modules/core/src/internal/split_edges.cpp @@ -21,6 +21,7 @@ #include #include +#include #include #include @@ -64,6 +65,31 @@ void interpolate_row( } } +template +void interpolate_row( + Eigen::MatrixBase& data, + Index row_to, + Index row_from0, + Index row_from1, + Scalar t) +{ + la_debug_assert(row_to < static_cast(data.rows())); + la_debug_assert(row_from0 < static_cast(data.rows())); + la_debug_assert(row_from1 < static_cast(data.rows())); + using ValueType = typename Derived::Scalar; + + if constexpr (std::is_integral_v) { + data.row(row_to) = (data.row(row_from0).template cast() * t + + data.row(row_from1).template cast() * (1 - t)) + .array() + .round() + .template cast() + .eval(); + } else { + data.row(row_to) = data.row(row_from0) * t + data.row(row_from1) * (1 - t); + } +} + template struct LocalBuffersT { @@ -84,7 +110,6 @@ std::vector split_edges( mesh.initialize_edges(); const Index num_output_vertices = mesh.get_num_vertices(); const Index num_input_facets = mesh.get_num_facets(); - const Index num_edges = mesh.get_num_edges(); std::vector original_triangle_index; std::vector split_triangles; @@ -94,8 +119,12 @@ std::vector split_edges( split_triangles_offsets.reserve(num_input_facets + 1); split_triangles_offsets.push_back(0); - // Compute the number of additional facets + // Compute the number of additional facets, and the CSR layout of split points per split facet. + std::vector parent_offsets; + parent_offsets.reserve(num_input_facets + 1); + parent_offsets.push_back(0); Index num_additional_facets = 0; + Index num_split_points = 0; for (Index fid = 0; fid < num_input_facets; fid++) { if (!active_facet(fid)) continue; Index chain_size = 3; @@ -107,23 +136,16 @@ std::vector split_edges( } if (chain_size == 3) continue; num_additional_facets += chain_size - 2; + num_split_points += chain_size - 3; split_triangles_offsets.push_back(num_additional_facets * 3); original_triangle_index.push_back(fid); + parent_offsets.push_back(num_split_points); } split_triangles.resize(num_additional_facets * 3, invalid()); - // Compute parent edge for each split point. - std::vector> parent_edge( - num_output_vertices, - {invalid(), invalid()}); - for (Index eid = 0; eid < num_edges; eid++) { - auto split_pts = get_edge_split_pts(eid); - for (auto vid : split_pts) { - if (parent_edge[vid][0] == invalid()) { - parent_edge[vid] = mesh.get_edge_vertices(eid); - } - } - } + // Parent edge {vertex, end0, end1} per split point, grouped per facet via `parent_offsets` and + // sorted by vertex so each facet resolves a shared vertex's edge by binary search, not a scan. + std::vector> parent_edges(num_split_points); using LocalBuffers = LocalBuffersT; tbb::enumerable_thread_specific local_buffers; @@ -138,6 +160,7 @@ std::vector split_edges( auto corners = mesh.get_facet_vertices(fid); Index corner_index[3] = {invalid(), invalid(), invalid()}; + Index parent_cursor = parent_offsets[idx]; for (Index i = 0; i < 3; i++) { Index v0 = corners[i]; chain.push_back(v0); @@ -146,15 +169,24 @@ std::vector split_edges( Index eid = mesh.get_edge(fid, i); la_debug_assert(eid != invalid()); auto split_pts = get_edge_split_pts(eid); - if (mesh.get_edge_vertices(eid)[0] == v0) { + auto edge_verts = mesh.get_edge_vertices(eid); + for (auto vid : split_pts) { + parent_edges[parent_cursor] = {vid, edge_verts[0], edge_verts[1]}; + ++parent_cursor; + } + if (edge_verts[0] == v0) { // Edge is oriented in the same direction as the triangle. chain.insert(chain.end(), split_pts.begin(), split_pts.end()); } else { // Edge is oriented in the opposite direction. - la_debug_assert(mesh.get_edge_vertices(eid)[1] == v0); + la_debug_assert(edge_verts[1] == v0); chain.insert(chain.end(), split_pts.rbegin(), split_pts.rend()); } } + la_debug_assert(parent_cursor == parent_offsets[idx + 1]); + std::sort( + parent_edges.begin() + parent_offsets[idx], + parent_edges.begin() + parent_offsets[idx + 1]); la_debug_assert(chain.size() > 3); size_t num_sub_triangles = chain.size() - 2; size_t curr_size = split_triangles_offsets[idx]; @@ -212,7 +244,8 @@ std::vector split_edges( return t; }; - auto barycentric_coordinates = [&](Index ori_fid, Index vid) -> std::array { + auto barycentric_coordinates = + [&](Index ori_fid, Index vid, Index p_begin, Index p_end) -> std::array { auto vid_0 = facets(ori_fid, 0); auto vid_1 = facets(ori_fid, 1); auto vid_2 = facets(ori_fid, 2); @@ -220,9 +253,14 @@ std::vector split_edges( if (vid == vid_1) return {0, 1, 0}; if (vid == vid_2) return {0, 0, 1}; - auto [v0, v1] = parent_edge[vid]; - la_debug_assert(v0 != invalid()); - la_debug_assert(v1 != invalid()); + auto it = std::lower_bound( + parent_edges.begin() + p_begin, + parent_edges.begin() + p_end, + vid, + [](const std::array& e, Index v) { return e[0] < v; }); + la_debug_assert(it != parent_edges.begin() + p_end && (*it)[0] == vid); + const Index v0 = (*it)[1]; + const Index v1 = (*it)[2]; auto t = edge_barycentric(v0, v1, vid); if (vid_0 != v0 && vid_0 != v1) { if (v0 == vid_1) { @@ -246,6 +284,22 @@ std::vector split_edges( } }; + // Barycentric coordinates of each new corner within its original facet, computed once and + // reused across every corner/indexed attribute rather than recomputed per attribute. + std::vector> corner_bc(static_cast(num_additional_facets) * 3); + tbb::parallel_for(size_t(0), original_triangle_index.size(), [&](size_t i) { + Index ori_fid = original_triangle_index[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 vid = facets(fid, k); + corner_bc[j * 3 + k] = + barycentric_coordinates(ori_fid, vid, parent_offsets[i], parent_offsets[i + 1]); + } + } + }); + auto map_corner_attribute = [&](auto&& data, auto&& corner_to_index) { for (size_t i = 0; i < original_triangle_index.size(); i++) { Index ori_fid = original_triangle_index[i]; @@ -259,8 +313,7 @@ std::vector split_edges( for (Index k = 0; k < 3; k++) { Index curr_cid = fid * 3 + k; - Index vid = facets(fid, k); - auto bc = barycentric_coordinates(ori_fid, vid); + auto& bc = corner_bc[j * 3 + k]; interpolate_row( data, corner_to_index(curr_cid), @@ -300,10 +353,164 @@ std::vector split_edges( return original_triangle_index; } -#define LA_X_split_edges(_, Scalar, Index) \ - template LA_CORE_API std::vector split_edges( \ - SurfaceMesh&, \ - function_ref(Index)>, \ +template +std::vector split_edges_only( + SurfaceMesh& mesh, + function_ref(Index)> get_edge_split_pts, + function_ref active_facet) +{ + const Index dim = mesh.get_dimension(); + mesh.initialize_edges(); + const Index num_input_facets = mesh.get_num_facets(); + const Index num_input_corners = mesh.get_num_corners(); + + auto vertices = vertex_view(mesh); + auto edge_barycentric = [&](Index v0, Index v1, Index v) { + Scalar diff = 0; + Scalar t = 0; + for (Index i = 0; i < dim; i++) { + if (vertices(v0, i) == vertices(v1, i)) continue; + Scalar d = vertices(v0, i) - vertices(v1, i); + if (std::abs(d) > diff) { + t = (vertices(v, i) - vertices(v1, i)) / d; + diff = std::abs(d); + } + } + la_runtime_assert( + diff > 0, + "Degenerate edge: cannot place a split point on an edge whose endpoints have " + "identical coordinates"); + return t; + }; + + std::vector original_facet_index; + std::vector new_facet_sizes; + std::vector new_facet_indices; + std::vector corner_src0; + std::vector corner_src1; + std::vector corner_weight0; + + std::vector chain; + for (Index fid = 0; fid < num_input_facets; fid++) { + if (!active_facet(fid)) continue; + + const Index facet_size = mesh.get_facet_size(fid); + const Index corner_begin = mesh.get_facet_corner_begin(fid); + auto corners = mesh.get_facet_vertices(fid); + + chain.clear(); + const size_t chain_start = new_facet_indices.size(); + for (Index i = 0; i < facet_size; i++) { + const Index v0 = corners[i]; + const Index i_next = (i + 1 == facet_size) ? 0 : i + 1; + const Index v1 = corners[i_next]; + const Index corner0 = corner_begin + i; + const Index corner1 = corner_begin + i_next; + + chain.push_back(v0); + corner_src0.push_back(corner0); + corner_src1.push_back(invalid()); + corner_weight0.push_back(1); + + const Index eid = mesh.get_edge(fid, i); + la_debug_assert(eid != invalid()); + auto split_pts = get_edge_split_pts(eid); + const bool same_direction = (mesh.get_edge_vertices(eid)[0] == v0); + la_debug_assert(same_direction || mesh.get_edge_vertices(eid)[1] == v0); + for (size_t k = 0; k < split_pts.size(); k++) { + Index vid = same_direction ? split_pts[k] : split_pts[split_pts.size() - 1 - k]; + chain.push_back(vid); + corner_src0.push_back(corner0); + corner_src1.push_back(corner1); + corner_weight0.push_back(edge_barycentric(v0, v1, vid)); + } + } + + if (static_cast(chain.size()) == facet_size) { + // No split points were introduced on this facet, discard the buffered corner data. + corner_src0.resize(chain_start); + corner_src1.resize(chain_start); + corner_weight0.resize(chain_start); + continue; + } + + original_facet_index.push_back(fid); + new_facet_sizes.push_back(static_cast(chain.size())); + new_facet_indices.insert(new_facet_indices.end(), chain.begin(), chain.end()); + } + + if (original_facet_index.empty()) return {}; + + // Clear edge data structure since the edge data is costly to maintain, + // and we are about to change facet connectivity anyway. + mesh.clear_edges(); + + mesh.add_hybrid( + {new_facet_sizes.data(), new_facet_sizes.size()}, + {new_facet_indices.data(), new_facet_indices.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_facet_index.size(); i++) { + data.row(num_input_facets + static_cast(i)) = + data.row(original_facet_index[i]); + } + }); + + // 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); + for (size_t j = 0; j < corner_src0.size(); j++) { + Index row_to = num_input_corners + static_cast(j); + if (corner_src1[j] == invalid()) { + data.row(row_to) = data.row(corner_src0[j]); + } else { + interpolate_row( + data, + row_to, + corner_src0[j], + corner_src1[j], + corner_weight0[j]); + } + } + }); + + // Map back to indexed attributes + for (auto attr_id : indexed_attribute_ids) { + map_attribute_in_place(mesh, attr_id, AttributeElement::Indexed); + } + + return original_facet_index; +} + +#define LA_X_split_edges(_, Scalar, Index) \ + template LA_CORE_API std::vector split_edges( \ + SurfaceMesh&, \ + function_ref(Index)>, \ + function_ref); \ + template LA_CORE_API std::vector split_edges_only( \ + SurfaceMesh&, \ + function_ref(Index)>, \ function_ref); LA_SURFACE_MESH_X(split_edges, 0) + } // namespace lagrange::internal diff --git a/modules/core/src/transform_mesh.cpp b/modules/core/src/transform_mesh.cpp index 0f7d1f5a..11c195be 100644 --- a/modules/core/src/transform_mesh.cpp +++ b/modules/core/src/transform_mesh.cpp @@ -34,6 +34,55 @@ namespace lagrange { +// Apply the geometric transform to one attribute's value matrix in place. Templated on +// only -- deliberately NOT on the mesh Index type. The heavy Eigen +// transform/cast machinery here is the dominant contributor to this TU's object size; keeping it +// Index-free means it is instantiated once per (Scalar, Dimension, ValueType) instead of being +// duplicated across every mesh index type. The row counter uses Eigen::Index for the same reason. +template +void transform_attribute_values( + RowMatrixView values, + AttributeUsage usage, + const Eigen::Transform& transform, + const Eigen::Matrix& cotransform, + const TransformOptions& options, + bool is_reflection) +{ + // Select higher-precision type between Scalar and ValueType + constexpr bool is_value_type_better = sizeof(ValueType) > sizeof(Scalar); + using HigherPrecisionType = std::conditional_t; + + auto A = transform.template cast(); + auto L = transform.linear().template cast(); + auto coL = cotransform.template cast(); + auto X = values.template cast().template leftCols().transpose(); + auto set = [&](auto&& Y) { + values.template leftCols() = Y.transpose().template cast(); + }; + HigherPrecisionType sign(options.reorient && is_reflection ? -1 : 1); + switch (usage) { + case AttributeUsage::Position: set(A * X); break; + case AttributeUsage::Normal: + set(sign * coL * X); + if (options.normalize_normals) { + tbb::parallel_for(Eigen::Index(0), Eigen::Index(values.rows()), [&](Eigen::Index c) { + values.row(c).template head().stableNormalize(); + }); + } + break; + case AttributeUsage::Tangent: [[fallthrough]]; + case AttributeUsage::Bitangent: + set(sign * L * X); + if (options.normalize_tangents_bitangents) { + tbb::parallel_for(Eigen::Index(0), Eigen::Index(values.rows()), [&](Eigen::Index c) { + values.row(c).template head().stableNormalize(); + }); + } + break; + default: break; + } +} + template void transform_mesh_internal( SurfaceMesh& mesh, @@ -51,8 +100,10 @@ void transform_mesh_internal( using AttributeType = std::decay_t; using ValueType = typename AttributeType::ValueType; + const AttributeUsage usage = attr_read.get_usage(); + // Skip if we don't need to modify the attribute (to avoid triggering copy-on-write) - switch (attr_read.get_usage()) { + switch (usage) { case AttributeUsage::Position: case AttributeUsage::Normal: case AttributeUsage::Tangent: @@ -60,59 +111,35 @@ void transform_mesh_internal( default: return; } - // Select higher-precision type between Scalar and ValueType - constexpr bool is_value_type_better = sizeof(ValueType) > sizeof(Scalar); - using HigherPrecisionType = std::conditional_t; - - // Apply geometric transform. - auto transform_values = [&](auto&& values) { - auto A = transform.template cast(); - auto L = transform.linear().template cast(); - auto coL = cotransform.template cast(); - auto X = values.template cast() - .template leftCols() - .transpose(); - auto set = [&](auto&& Y) { - values.template leftCols() = Y.transpose().template cast(); - }; - HigherPrecisionType sign(options.reorient && is_reflection ? -1 : 1); - if (!included_usages.test(attr_read.get_usage())) { + // Filter by value type/indexed. The heavy Eigen transform lives in the Index-free + // transform_attribute_values helper so it is not re-instantiated per mesh index type. + if constexpr (std::is_floating_point_v) { + // The included-usages check stays inside the floating-point branch so that a + // non-floating attribute of a transformable usage still triggers the type error below, + // matching the original behavior (the check used to live inside transform_values). + if (!included_usages.test(usage)) { logger().debug("Skipping transform for attribute: {}", name); return; } - switch (attr_read.get_usage()) { - case AttributeUsage::Position: set(A * X); break; - case AttributeUsage::Normal: - 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(); - }); - } - break; - case AttributeUsage::Tangent: [[fallthrough]]; - case AttributeUsage::Bitangent: - 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(); - }); - } - break; - default: break; - } - }; - - // Filter by value type/indexed - if constexpr (std::is_floating_point_v) { if constexpr (AttributeType::IsIndexed) { auto& attr = mesh.template ref_indexed_attribute(name); - transform_values(matrix_ref(attr.values())); + transform_attribute_values( + matrix_ref(attr.values()), + usage, + transform, + cotransform, + options, + is_reflection); } else { - transform_values(attribute_matrix_ref(mesh, name)); + transform_attribute_values( + attribute_matrix_ref(mesh, name), + usage, + transform, + cotransform, + options, + is_reflection); } } else { - LA_IGNORE(transform_values); std::string_view type_name; if constexpr (AttributeType::IsIndexed) { type_name = internal::value_type_name(attr_read.values()); @@ -122,7 +149,7 @@ void transform_mesh_internal( throw Error(format( "Invalid attribute value type ({}) for attribute usage: {}", type_name, - internal::to_string(attr_read.get_usage()))); + internal::to_string(usage))); } }); diff --git a/modules/core/src/triangulate_polygonal_facets.cpp b/modules/core/src/triangulate_polygonal_facets.cpp index 4fce7398..efd1d6a7 100644 --- a/modules/core/src/triangulate_polygonal_facets.cpp +++ b/modules/core/src/triangulate_polygonal_facets.cpp @@ -164,7 +164,8 @@ template void triangulate_polygonal_facets_earcut( SurfaceMesh& mesh, bool preserve_edges, - bool preserve_points) + bool preserve_points, + function_ref* should_triangulate = nullptr) { LAGRANGE_ZONE_SCOPED; @@ -198,11 +199,11 @@ void triangulate_polygonal_facets_earcut( if (!preserve_edges) { to_remove[f] = true; } - } else if (facet_size == 4) { + } else if (facet_size == 4 && (!should_triangulate || (*should_triangulate)(f))) { // Triangulate quad to_remove[f] = true; append_triangles_from_quad(mesh, f, new_to_old_corners, new_to_old_facets); - } else if (facet_size > 4) { + } else if (facet_size > 4 && (!should_triangulate || (*should_triangulate)(f))) { // Triangulate polygon via ear cutting to_remove[f] = true; append_triangles_from_polygon( @@ -311,7 +312,8 @@ template void triangulate_polygonal_facets_centroid_fan( SurfaceMesh& mesh, bool preserve_edges, - bool preserve_points) + bool preserve_points, + function_ref* should_triangulate = nullptr) { if (mesh.is_triangle_mesh()) { return; @@ -334,7 +336,8 @@ void triangulate_polygonal_facets_centroid_fan( for (Index fid = 0; fid < old_num_facets; ++fid) { const auto facet_size = mesh.get_facet_size(fid); if (facet_size != 3 && !(preserve_edges && facet_size == 2) && - !(preserve_points && facet_size == 1)) { + !(preserve_points && facet_size == 1) && + (!should_triangulate || (*should_triangulate)(fid))) { auto f = mesh.get_facet_vertices(fid); facets_to_remove.push_back(fid); @@ -370,6 +373,12 @@ void triangulate_polygonal_facets_centroid_fan( } } + if (facets_to_remove.empty()) { + // Nothing selected: leave the mesh untouched and avoid add_triangles() with empty data, + // which asserts on meshes that carry edge/connectivity information. + return; + } + mesh.add_vertices(new_vertex_count, {centroids.data(), centroids.size()}); mesh.add_triangles(new_triangle_count, {new_triangles.data(), new_triangles.size()}); @@ -547,10 +556,38 @@ void triangulate_polygonal_facets( } } +template +void triangulate_polygonal_facets( + SurfaceMesh& mesh, + function_ref should_triangulate, + const TriangulationOptions& options) +{ + switch (options.scheme) { + case TriangulationOptions::Scheme::Earcut: + triangulate_polygonal_facets_earcut( + mesh, + options.preserve_edges, + options.preserve_points, + &should_triangulate); + break; + case TriangulationOptions::Scheme::CentroidFan: + triangulate_polygonal_facets_centroid_fan( + mesh, + options.preserve_edges, + options.preserve_points, + &should_triangulate); + break; + } +} + // Iterate over mesh (scalar, index) types #define LA_X_triangulate_polygonal_facets(_, Scalar, Index) \ template LA_CORE_API void triangulate_polygonal_facets( \ SurfaceMesh& mesh, \ + const TriangulationOptions& options); \ + template LA_CORE_API void triangulate_polygonal_facets( \ + SurfaceMesh& mesh, \ + function_ref should_triangulate, \ const TriangulationOptions& options); LA_SURFACE_MESH_X(triangulate_polygonal_facets, 0) diff --git a/modules/core/src/weld_indexed_attribute.cpp b/modules/core/src/weld_indexed_attribute.cpp index d206546f..c99a670e 100644 --- a/modules/core/src/weld_indexed_attribute.cpp +++ b/modules/core/src/weld_indexed_attribute.cpp @@ -18,9 +18,8 @@ #include #include #include +#include #include -#include -#include #include // clang-format off @@ -31,51 +30,43 @@ // clang-format on #include +#include +#include +#include #include +#include namespace lagrange { namespace { -// Unsigned integer type using the most significant bit as a flag. -template -struct IndexWithFlagT +// Maps a value type to its corresponding real (non-integer) type for tolerance computations. +template +struct NonIntegerT { - explicit IndexWithFlagT(Index i = 0) { set_index(i); } - - void set_index(Index i) - { - i &= ~Mask; // Clear the flag - i |= (m_value & Mask); // Set the flag if it was set - m_value = i; - } - - Index index() const - { - return m_value & ~Mask; // Clear the flag - } - - void set_flag(bool enabled) - { - if (enabled) { - m_value |= Mask; // Set the flag - } else { - m_value &= ~Mask; // Clear the flag - } - } + using type = double; +}; +template <> +struct NonIntegerT +{ + using type = float; +}; +template <> +struct NonIntegerT +{ + using type = double; +}; - bool flag() const - { - return (m_value & Mask) != 0; // Check the flag +// Returns the default relative tolerance (equivalent to Eigen::NumTraits::dummy_precision()). +template +constexpr T default_rel_tolerance() +{ + if constexpr (std::is_same_v) { + return T(1e-5); + } else { + return T(1e-12); } - - Index m_value = 0; - - static constexpr Index Mask = Index(1) << (8 * sizeof(Index) - 1); - - static_assert(std::is_unsigned_v, "Index must be unsigned"); - static_assert(Mask != 0 && (Mask << 1 == 0)); -}; +} template struct IndexAndCornerT @@ -84,56 +75,117 @@ struct IndexAndCornerT Index corner; }; -template -void weld_indexed_attribute( - SurfaceMesh& mesh, - IndexedAttribute& attr, - span exclude_vertices, - bool merge_across_vertices, - Func equal) +// Check whether two rows of values are approximately equal. +// Returns true if all channels satisfy |a[c] - b[c]| <= atol + rtol * |b[c]|, +// and optionally if the angle between the two vectors is within the threshold. +template +bool rows_are_close( + const ValueType* row_a, + const ValueType* row_b, + size_t num_channels, + double eps_rel, + double eps_abs, + double cos_angle_abs) { - std::vector exclude_vertices_mask(mesh.get_num_vertices(), false); - for (auto vi : exclude_vertices) { - la_debug_assert(vi < mesh.get_num_vertices()); - exclude_vertices_mask[vi] = true; + using RealType = typename NonIntegerT::type; + + // Element-wise tolerance check. + for (size_t c = 0; c < num_channels; ++c) { + RealType diff = std::abs(static_cast(row_a[c]) - static_cast(row_b[c])); + RealType tol = static_cast(eps_abs) + + static_cast(eps_rel) * std::abs(static_cast(row_b[c])); + // Use !(diff <= tol) so that NaN causes the comparison to fail (returns false). + if (!(diff <= tol)) return false; } - const bool had_edges = mesh.has_edges(); - mesh.initialize_edges(); - const auto _ = make_scope_guard([&]() { - if (!had_edges) { - mesh.clear_edges(); + // Angle check (only if cos_angle_abs <= 1). + if (cos_angle_abs <= 1.0) { + RealType dot = 0; + RealType norm_a_sq = 0; + RealType norm_b_sq = 0; + for (size_t c = 0; c < num_channels; ++c) { + RealType a = static_cast(row_a[c]); + RealType b = static_cast(row_b[c]); + dot += a * b; + norm_a_sq += a * a; + norm_b_sq += b * b; } - }); - auto& attr_values = attr.values(); - auto values = matrix_view(attr_values); - auto corner_to_value = vector_ref(attr.indices()); + RealType threshold = + static_cast(cos_angle_abs) * std::sqrt(norm_a_sq) * std::sqrt(norm_b_sq); + // Use !(dot >= threshold) so that NaN causes the comparison to fail (returns false). + if (!(dot >= threshold)) { + return false; + } + } - const Index num_vertices = mesh.get_num_vertices(); - const Index num_values = static_cast(values.rows()); - const Index num_corners = mesh.get_num_corners(); + return true; +} - using IndexWithFlag = IndexWithFlagT; - std::vector corner_map(num_corners); - for (Index c = 0; c < num_corners; ++c) { - corner_map[c].set_index(c); +// Check whether two rows of values are exactly equal (per-element ==). +template +bool rows_are_equal(const ValueType* row_a, const ValueType* row_b, size_t num_channels) +{ + return std::equal(row_a, row_a + num_channels, row_b); +} + +// Check whether any channel in a row equals the invalid sentinel value. +template +bool row_has_invalid(const ValueType* row, size_t num_channels) +{ + const ValueType inv = lagrange::invalid(); + for (size_t c = 0; c < num_channels; ++c) { + if (row[c] == inv) return true; } + return false; +} - auto find_and_compress = [&](Index c) { - while (corner_map[c].index() != c) { - corner_map[c].set_index(corner_map[corner_map[c].index()].index()); - c = corner_map[c].index(); +// Assign a reduced index to the root corner of a group. Updates num_reduced in place. +template +void process_root( + Index c, + Index& num_reduced, + const std::vector& group_flagged, + span corner_to_value, + std::vector& corner_to_reduced, + std::vector& index_to_reduced) +{ + if (corner_to_reduced[c] != invalid()) { + // If the root corner has already been processed, we can skip it. + return; + } + if (group_flagged[c]) { + // If the group is flagged, it means a merge happened, and we assign a new index to + // the corner group. + corner_to_reduced[c] = num_reduced++; + } else { + // If the group is not flagged, we can preserve the original index. In other words, + // we assign a reduced index based on the original index associated to the corner, + // not based on the corner group. + Index i = corner_to_value[c]; + if (index_to_reduced[i] == invalid()) { + index_to_reduced[i] = num_reduced++; } - return c; - }; - - auto merge_groups = [&](Index c1, Index c2, bool flag = false) { - auto r1 = find_and_compress(c1); - auto r2 = find_and_compress(c2); - corner_map[r2].set_index(r1); - corner_map[r1].set_flag(corner_map[r1].flag() || corner_map[r2].flag() || flag); - }; + corner_to_reduced[c] = index_to_reduced[i]; + } + la_debug_assert(corner_to_reduced[c] != invalid()); +} +// Merge corners around each vertex that share the same index or have similar values. +// +// Templated on `Index` only: the per-vertex corner adjacency is supplied as an index-only CSR +// (`vertex_to_corners`) and the value comparison is type-erased through `values_close`. This keeps +// the heavy parallel machinery (TBB, SmallVector, DisjointSets) compiled once per index type +// rather than once per (ValueType, Scalar, Index) combination. +template +void merge_corners_per_vertex( + Index num_vertices, + const internal::InverseMapping& vertex_to_corners, + const std::vector& exclude_vertices_mask, + span corner_to_value, + function_ref values_close, + DisjointSets& corner_map, + std::vector& corner_flagged) +{ // Sort and find duplicate values shared by corners around the same vertex. using IndexAndCorner = IndexAndCornerT; tbb::parallel_for( @@ -143,8 +195,8 @@ void weld_indexed_attribute( if (exclude_vertices_mask[vi]) continue; SmallVector involved_indices_and_corners; - mesh.foreach_corner_around_vertex(vi, [&](Index ci) { - involved_indices_and_corners.push_back({corner_to_value(ci), ci}); + vertex_to_corners.foreach_mapped_to(vi, [&](Index ci) { + involved_indices_and_corners.push_back({corner_to_value[ci], ci}); }); la_debug_assert(involved_indices_and_corners.size() > 0); @@ -160,7 +212,7 @@ void weld_indexed_attribute( return (x.index != it_begin->index); }); for (auto it = it_begin; it != it_end; ++it) { - corner_map[it->corner].set_index(it_begin->corner); + corner_map.merge(it_begin->corner, it->corner); } it_begin = it_end; } @@ -174,7 +226,7 @@ void weld_indexed_attribute( // Update corner associated to the uniqued index to be the root of the group for (auto itr = first; itr != last; itr++) { Index& c = itr->corner; - c = corner_map[c].index(); + c = corner_map.find(c); } for (auto itr = first; itr != last; itr++) { @@ -183,95 +235,130 @@ void weld_indexed_attribute( // If the corner is not the root of the group, it means it has been merged with // another corner in this inner loop, and we don't need to compare against all // other uniqued indices again. - if (corner_map[c1].index() != c1) continue; + if (corner_map.find(c1) != c1) continue; // Quadratic loop to search for corners with similar values. for (auto itr2 = std::next(itr); itr2 != last; itr2++) { const auto& [i2, c2] = *itr2; - if (equal(i1, i2)) { + if (values_close(i1, i2)) { // Flag any corner group containing merged values. - merge_groups(c1, c2, true); + Index root = corner_map.merge(c1, c2); + corner_flagged[root] = 1; } } } } }); +} - if (merge_across_vertices) { - // Merge corner groups that share indices - auto index_to_corner = internal::invert_mapping( - {corner_to_value.data(), static_cast(num_corners)}, - num_values); - for (Index i = 0; i < num_values; i++) { - auto it_begin = index_to_corner.data.begin() + index_to_corner.offsets[i]; - auto it_end = index_to_corner.data.begin() + index_to_corner.offsets[i + 1]; - if (it_begin == it_end) continue; - for (auto it = it_begin + 1; it != it_end; ++it) { - merge_groups(*it_begin, *it); - } - } - } - +// Assign reduced indices to all corners based on the disjoint sets and flags. +template +Index assign_reduced_indices( + Index num_corners, + Index num_values, + DisjointSets& corner_map, + const std::vector& corner_flagged, + span corner_to_value, + std::vector& corner_to_reduced) +{ // Propagate flags to roots after all merges are done std::vector group_flagged(num_corners, false); for (Index c = 0; c < num_corners; ++c) { - if (corner_map[c].flag()) { - Index rc = find_and_compress(c); + if (corner_flagged[c]) { + Index rc = corner_map.find(c); group_flagged[rc] = true; } } Index num_reduced = 0; - std::vector corner_to_reduced(num_corners, invalid()); std::vector index_to_reduced(num_values, invalid()); - auto process_root = [&](Index c) { - if (corner_to_reduced[c] != invalid()) { - // If the root corner has already been processed, we can skip it. - return; - } - if (group_flagged[c]) { - // If the group is flagged, it means a merge happened, and we assign a new index to - // the corner group. - corner_to_reduced[c] = num_reduced++; - } else { - // If the group is not flagged, we can preserve the original index. In other words, - // we assign a reduced index based on the original index associated to the corner, - // not based on the corner group. - Index i = corner_to_value[c]; - if (index_to_reduced[i] == invalid()) { - index_to_reduced[corner_to_value[c]] = num_reduced++; - } - corner_to_reduced[c] = index_to_reduced[corner_to_value[c]]; - } - la_debug_assert(corner_to_reduced[c] != invalid()); - }; - // Assign reduced indices to corners. for (Index c = 0; c < num_corners; ++c) { - Index rc = find_and_compress(c); - process_root(rc); + Index rc = corner_map.find(c); + process_root( + rc, + num_reduced, + group_flagged, + corner_to_value, + corner_to_reduced, + index_to_reduced); if (rc != c) { corner_to_reduced[c] = corner_to_reduced[rc]; la_debug_assert(corner_to_reduced[c] != invalid()); } } - if (num_reduced == num_values) { - // Nothing to weld. - return; + return num_reduced; +} + +// Run the index-only welding pipeline: merge corners around each vertex, optionally merge corner +// groups that share the same value index across vertices, then assign reduced indices. +// +// Templated on `Index` only (not ValueType/Scalar): the value comparison is type-erased through +// `values_close`, so the heavy disjoint-set machinery is compiled once per index type. Fills +// `corner_to_reduced` and returns the number of reduced (welded) values. +template +Index weld_core( + Index num_corners, + Index num_vertices, + Index num_values, + const internal::InverseMapping& vertex_to_corners, + const std::vector& exclude_vertices_mask, + span corner_to_value, + function_ref values_close, + bool merge_across_vertices, + std::vector& corner_to_reduced) +{ + DisjointSets corner_map(num_corners); + std::vector corner_flagged(num_corners, 0); + + merge_corners_per_vertex( + num_vertices, + vertex_to_corners, + exclude_vertices_mask, + corner_to_value, + values_close, + corner_map, + corner_flagged); + + if (merge_across_vertices) { + // Merge corner groups that share indices + auto index_to_corner = internal::invert_mapping( + span(corner_to_value.data(), static_cast(num_corners)), + num_values); + for (Index i = 0; i < num_values; i++) { + auto it_begin = index_to_corner.data.begin() + index_to_corner.offsets[i]; + auto it_end = index_to_corner.data.begin() + index_to_corner.offsets[i + 1]; + if (it_begin == it_end) continue; + for (auto it = it_begin + 1; it != it_end; ++it) { + corner_map.merge(*it_begin, *it); + } + } } - auto reduced_to_corner = - internal::invert_mapping({corner_to_reduced.data(), corner_to_reduced.size()}, num_reduced); + corner_to_reduced.assign(num_corners, invalid()); + return assign_reduced_indices( + num_corners, + num_values, + corner_map, + corner_flagged, + corner_to_value, + corner_to_reduced); +} - Attribute attr_welded_values( - attr_values.get_element_type(), - attr_values.get_usage(), - attr_values.get_num_channels()); - attr_welded_values.resize_elements(num_reduced); - auto welded_values = matrix_ref(attr_welded_values); - welded_values.setZero(); +// Within each reduced group, sort the member corners by value index and move the unique ones to the +// front of the group's slice (via std::unique). The group offsets are NOT updated and no elements +// are erased; instead the per-group count of unique members is returned, and callers must read only +// the first `unique_count[ri]` entries of each group. Templated on `Index` only, so the heavy +// parallel sort is compiled once per index type rather than once per value type. +template +std::vector dedup_groups_by_value( + internal::InverseMapping& reduced_to_corner, + Index num_reduced, + span corner_to_value) +{ + std::vector unique_count(static_cast(num_reduced)); tbb::parallel_for(Index(0), num_reduced, [&](Index ri) { auto it_begin = reduced_to_corner.data.begin() + reduced_to_corner.offsets[ri]; auto it_end = reduced_to_corner.data.begin() + reduced_to_corner.offsets[ri + 1]; @@ -280,38 +367,131 @@ void weld_indexed_attribute( tbb::parallel_sort(it_begin, it_end, [&](Index ci, Index cj) { return corner_to_value[ci] < corner_to_value[cj]; }); - it_end = std::unique(it_begin, it_end, [&](Index ci, Index cj) { + auto new_end = std::unique(it_begin, it_end, [&](Index ci, Index cj) { return corner_to_value[ci] == corner_to_value[cj]; }); - for (auto it = it_begin; it != it_end; ++it) { - welded_values.row(ri) += values.row(corner_to_value[*it]); + unique_count[ri] = static_cast(std::distance(it_begin, new_end)); + }); + return unique_count; +} + +// Compute the welded attribute values by averaging merged groups. Only the final accumulation is +// ValueType-dependent; the index-only grouping/dedup is handled by `dedup_groups_by_value`. +template +void compute_welded_values( + Index num_reduced, + span values_data, + span corner_to_value, + const std::vector& corner_to_reduced, + size_t num_channels, + Attribute& attr_values) +{ + auto reduced_to_corner = + internal::invert_mapping({corner_to_reduced.data(), corner_to_reduced.size()}, num_reduced); + std::vector unique_count = + dedup_groups_by_value(reduced_to_corner, num_reduced, corner_to_value); + + Attribute attr_welded_values( + attr_values.get_element_type(), + attr_values.get_usage(), + attr_values.get_num_channels()); + attr_welded_values.resize_elements(num_reduced); + auto welded_data = attr_welded_values.ref_all(); + std::fill(welded_data.begin(), welded_data.end(), ValueType(0)); + + tbb::parallel_for(Index(0), num_reduced, [&](Index ri) { + const Index* group = reduced_to_corner.data.data() + reduced_to_corner.offsets[ri]; + const Index num = unique_count[ri]; + ValueType* dst = welded_data.data() + static_cast(ri) * num_channels; + for (Index k = 0; k < num; ++k) { + const ValueType* src = + values_data.data() + static_cast(corner_to_value[group[k]]) * num_channels; + for (size_t ch = 0; ch < num_channels; ++ch) { + dst[ch] += src[ch]; + } } - auto num = std::distance(it_begin, it_end); if (num > 1) { - welded_values.row(ri) /= static_cast(num); + for (size_t ch = 0; ch < num_channels; ++ch) { + dst[ch] /= static_cast(num); + } } }); attr_values = std::move(attr_welded_values); - tbb::parallel_for(Index(0), num_corners, [&](auto c) { - corner_to_value[c] = corner_to_reduced[c]; - }); } -template -bool allclose( - const Eigen::DenseBase& a, - const Eigen::DenseBase& b, - const typename DerivedA::RealScalar& rtol = - Eigen::NumTraits::dummy_precision(), - const typename DerivedA::RealScalar& atol = - Eigen::NumTraits::epsilon(), - const typename DerivedA::RealScalar& cos_angle_abs = 1) +// Thin per-(ValueType, Index) shell: builds the type-erased value comparator and forwards the +// index-only topology to `weld_core`, then computes the averaged welded values. All the heavy +// algorithmic code lives in `weld_core` (Index-only) and `compute_welded_values`. +template +void weld_indexed_attribute_impl( + IndexedAttribute& attr, + Index num_corners, + Index num_vertices, + const internal::InverseMapping& vertex_to_corners, + const std::vector& exclude_vertices_mask, + bool merge_across_vertices, + double eps_rel, + double eps_abs, + double cos_angle_abs) { - // TODO: Use two different checks for absolute and relative tolerances. - return ((a.derived() - b.derived()).array().abs() <= (atol + rtol * b.derived().array().abs())) - .all() && - (cos_angle_abs > 1 || (a.derived().dot(b.derived()) >= - cos_angle_abs * a.derived().norm() * b.derived().norm())); + auto& attr_values = attr.values(); + auto& attr_indices = attr.indices(); + const size_t num_channels = attr_values.get_num_channels(); + span values_data = attr_values.get_all(); + span corner_to_value = attr_indices.ref_all(); + + const Index num_values = static_cast(attr_values.get_num_elements()); + + // The only ValueType-dependent part of the merge pipeline: comparison of two value rows. + auto values_close = [&](Index i, Index j) -> bool { + const ValueType* row_i = values_data.data() + static_cast(i) * num_channels; + const ValueType* row_j = values_data.data() + static_cast(j) * num_channels; + if (rows_are_equal(row_i, row_j, num_channels)) { + return true; + } + if (row_has_invalid(row_i, num_channels) || row_has_invalid(row_j, num_channels)) { + return false; + } + // Debug-only finiteness check (equivalent to Eigen's allFinite()). + if constexpr (std::is_floating_point_v) { + la_debug_assert(std::all_of(row_i, row_i + num_channels, [](ValueType v) { + return std::isfinite(v); + })); + la_debug_assert(std::all_of(row_j, row_j + num_channels, [](ValueType v) { + return std::isfinite(v); + })); + } + return rows_are_close(row_i, row_j, num_channels, eps_rel, eps_abs, cos_angle_abs); + }; + + std::vector corner_to_reduced; + Index num_reduced = weld_core( + num_corners, + num_vertices, + num_values, + vertex_to_corners, + exclude_vertices_mask, + span(corner_to_value.data(), corner_to_value.size()), + values_close, + merge_across_vertices, + corner_to_reduced); + + if (num_reduced == num_values) { + // Nothing to weld. + return; + } + + compute_welded_values( + num_reduced, + values_data, + span(corner_to_value.data(), corner_to_value.size()), + corner_to_reduced, + num_channels, + attr_values); + + tbb::parallel_for(Index(0), num_corners, [&](Index c) { + corner_to_value[c] = corner_to_reduced[c]; + }); } } // namespace @@ -322,51 +502,56 @@ void weld_indexed_attribute( AttributeId attr_id, const WeldOptions& options) { + // Extract index-only topology once, independent of the attribute's value type: the set of + // corners incident to each vertex. Using the intrinsic corner-to-vertex map (rather than edge + // connectivity) keeps the heavy welding pipeline free of the Scalar type and avoids the need to + // initialize/clear mesh edges. + const Index num_vertices = mesh.get_num_vertices(); + const Index num_corners = mesh.get_num_corners(); + span corner_to_vertex = mesh.get_corner_to_vertex().get_all(); + const internal::InverseMapping vertex_to_corners = + internal::invert_mapping(corner_to_vertex, num_vertices); + + std::vector exclude_vertices_mask(static_cast(num_vertices), false); + for (auto vi : options.exclude_vertices) { + la_debug_assert(vi < static_cast(num_vertices)); + exclude_vertices_mask[vi] = true; + } + lagrange::internal::visit_attribute_write(mesh, attr_id, [&](auto&& attr) { using AttributeType = std::decay_t; if constexpr (AttributeType::IsIndexed) { using ValueType = typename AttributeType::ValueType; - using RealType = typename Eigen::NumTraits::NonInteger; - auto values = matrix_view(attr.values()); - - const RealType eps_rel = options.epsilon_rel.has_value() - ? safe_cast(options.epsilon_rel.value()) - : Eigen::NumTraits::dummy_precision(); - const RealType eps_abs = options.epsilon_abs.has_value() - ? safe_cast(options.epsilon_abs.value()) - : Eigen::NumTraits::epsilon(); - - constexpr RealType INVALID_COS_ANGLE_ABS = 2; // Out of the valid range. - const RealType cos_angle_abs = options.angle_abs.has_value() - ? std::cos(options.angle_abs.value()) - : INVALID_COS_ANGLE_ABS; - weld_indexed_attribute( - mesh, + using RealType = typename NonIntegerT::type; + + // safe_cast validates that a caller-provided tolerance is representable in + // the + // attribute's real type, throwing (rather than silently saturating to +/-inf inside + // rows_are_close) if it overflows a lower-precision ValueType such as float. Done once + // here rather than per element comparison. + const double eps_rel = + options.epsilon_rel.has_value() + ? static_cast(safe_cast(options.epsilon_rel.value())) + : static_cast(default_rel_tolerance()); + const double eps_abs = + options.epsilon_abs.has_value() + ? static_cast(safe_cast(options.epsilon_abs.value())) + : static_cast(std::numeric_limits::epsilon()); + + constexpr double INVALID_COS_ANGLE_ABS = 2.0; // Out of the valid range. + const double cos_angle_abs = options.angle_abs.has_value() + ? std::cos(options.angle_abs.value()) + : INVALID_COS_ANGLE_ABS; + weld_indexed_attribute_impl( attr, - options.exclude_vertices, + num_corners, + num_vertices, + vertex_to_corners, + exclude_vertices_mask, options.merge_across_vertices, - [&, eps_rel, eps_abs, cos_angle_abs](Index i, Index j) -> bool { - if (values.row(i) == values.row(j)) { - return true; - } - const bool invalid_i = - (values.row(i).array() == lagrange::invalid()).any(); - const bool invalid_j = - (values.row(j).array() == lagrange::invalid()).any(); - if (invalid_i || invalid_j) { - // Along with the equality check above, this ensures that we only merge - // invalid values with other invalid values, and we don't merge valid values - // with invalid values. - return false; - } - la_debug_assert(values.row(i).allFinite() && values.row(j).allFinite()); - return allclose( - values.row(i).template cast(), - values.row(j).template cast(), - eps_rel, - eps_abs, - cos_angle_abs); - }); + eps_rel, + eps_abs, + cos_angle_abs); } }); } diff --git a/modules/core/tests/test_split_edges.cpp b/modules/core/tests/test_split_edges.cpp new file mode 100644 index 00000000..943d442b --- /dev/null +++ b/modules/core/tests/test_split_edges.cpp @@ -0,0 +1,334 @@ +/* + * Copyright 2026 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 + +TEST_CASE("internal::split_edges_only", "[core]") +{ + using Scalar = double; + using Index = uint32_t; + using namespace lagrange; + + // Kept as plain lambdas (not pre-wrapped in function_ref) so they stay alive for the + // duration of the calls below; function_ref is non-owning and would otherwise dangle. + auto always_active = [](Index) { return true; }; + + SECTION("Single triangle, split one edge into a quad") + { + 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); + + // Split point in the middle of edge (0, 1). + Index split_vertex = mesh.get_num_vertices(); + mesh.add_vertex({0.5, 0, 0}); + + mesh.initialize_edges(); + Index split_eid = mesh.find_edge_from_vertices(0, 1); + std::vector split_pts = {split_vertex}; + auto get_edge_split_pts = [&](Index eid) -> span { + if (eid == split_eid) return span(split_pts.data(), split_pts.size()); + return span(split_pts.data(), 0); + }; + + auto updated_facets = + internal::split_edges_only(mesh, get_edge_split_pts, always_active); + + REQUIRE(updated_facets.size() == 1); + CHECK(updated_facets[0] == 0); + + REQUIRE(mesh.get_num_facets() == 2); + CHECK(mesh.get_facet_size(0) == 3); // Original triangle is left untouched. + + REQUIRE(mesh.get_facet_size(1) == 4); + auto quad = mesh.get_facet_vertices(1); + CHECK(quad[0] == 0); + CHECK(quad[1] == split_vertex); + CHECK(quad[2] == 1); + CHECK(quad[3] == 2); + } + + SECTION("Triangulated square, split diagonal into two quads") + { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); + mesh.add_vertex({1, 0, 0}); + mesh.add_vertex({1, 1, 0}); + mesh.add_vertex({0, 1, 0}); + mesh.add_triangle(0, 1, 2); + mesh.add_triangle(0, 2, 3); + + // Split point in the middle of the diagonal (0, 2). + Index split_vertex = mesh.get_num_vertices(); + mesh.add_vertex({0.5, 0.5, 0}); + + mesh.initialize_edges(); + Index split_eid = mesh.find_edge_from_vertices(0, 2); + std::vector split_pts = {split_vertex}; + auto get_edge_split_pts = [&](Index eid) -> span { + if (eid == split_eid) return span(split_pts.data(), split_pts.size()); + return span(split_pts.data(), 0); + }; + + auto updated_facets = + internal::split_edges_only(mesh, get_edge_split_pts, always_active); + + REQUIRE(updated_facets.size() == 2); + CHECK(updated_facets[0] == 0); + CHECK(updated_facets[1] == 1); + + REQUIRE(mesh.get_num_facets() == 4); + CHECK(mesh.get_facet_size(0) == 3); // Original triangles are left untouched. + CHECK(mesh.get_facet_size(1) == 3); + + REQUIRE(mesh.get_facet_size(2) == 4); + auto quad0 = mesh.get_facet_vertices(2); + CHECK(quad0[0] == 0); + CHECK(quad0[1] == 1); + CHECK(quad0[2] == 2); + CHECK(quad0[3] == split_vertex); + + REQUIRE(mesh.get_facet_size(3) == 4); + auto quad1 = mesh.get_facet_vertices(3); + CHECK(quad1[0] == 0); + CHECK(quad1[1] == split_vertex); + CHECK(quad1[2] == 2); + CHECK(quad1[3] == 3); + } + + SECTION("Attribute propagation") + { + 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); + + mesh.create_attribute( + "facet_value", + AttributeElement::Facet, + AttributeUsage::Scalar); + attribute_vector_ref(mesh, "facet_value") << 42; + + mesh.create_attribute( + "corner_value", + AttributeElement::Corner, + AttributeUsage::Scalar); + // Corner 0 -> vertex 0, corner 1 -> vertex 1, corner 2 -> vertex 2. + attribute_vector_ref(mesh, "corner_value") << 10, 20, 40; + + auto indexed_id = mesh.create_attribute( + "indexed_value", + AttributeElement::Indexed, + AttributeUsage::Scalar, + 1); + { + auto& iattr = mesh.ref_indexed_attribute(indexed_id); + iattr.values().resize_elements(3); + auto values = iattr.values().ref_all(); + values[0] = 10; + values[1] = 20; + values[2] = 40; + auto indices = iattr.indices().ref_all(); + indices[0] = 0; + indices[1] = 1; + indices[2] = 2; + } + + // Split point in the middle of edge (0, 1). + Index split_vertex = mesh.get_num_vertices(); + mesh.add_vertex({0.5, 0, 0}); + + mesh.initialize_edges(); + Index split_eid = mesh.find_edge_from_vertices(0, 1); + std::vector split_pts = {split_vertex}; + auto get_edge_split_pts = [&](Index eid) -> span { + if (eid == split_eid) return span(split_pts.data(), split_pts.size()); + return span(split_pts.data(), 0); + }; + + auto updated_facets = + internal::split_edges_only(mesh, get_edge_split_pts, always_active); + REQUIRE(updated_facets.size() == 1); + REQUIRE(mesh.get_num_facets() == 2); + REQUIRE(mesh.get_facet_size(1) == 4); + + // New quad facet inherits the original facet's attribute value. + CHECK(attribute_vector_ref(mesh, "facet_value")[1] == 42); + + // New quad corners are [v0, split_vertex, v1, v2]; split_vertex is the midpoint of + // edge (v0, v1), so its value interpolates corner_value[0] and corner_value[1]. + auto corner_value = attribute_vector_ref(mesh, "corner_value"); + Index quad_corner_begin = mesh.get_facet_corner_begin(1); + CHECK(corner_value[quad_corner_begin + 0] == 10); + CHECK(corner_value[quad_corner_begin + 1] == 15); + CHECK(corner_value[quad_corner_begin + 2] == 20); + CHECK(corner_value[quad_corner_begin + 3] == 40); + + // Indexed attribute round-trips through the same interpolation. + auto& iattr = mesh.get_indexed_attribute("indexed_value"); + auto indexed_values = iattr.values().get_all(); + auto indexed_indices = iattr.indices().get_all(); + CHECK(indexed_values[indexed_indices[quad_corner_begin + 0]] == 10); + CHECK(indexed_values[indexed_indices[quad_corner_begin + 1]] == 15); + CHECK(indexed_values[indexed_indices[quad_corner_begin + 2]] == 20); + CHECK(indexed_values[indexed_indices[quad_corner_begin + 3]] == 40); + } + + SECTION("Splitting a geometrically degenerate edge throws") + { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); + mesh.add_vertex({0, 0, 0}); // Coincident with vertex 0, so edge (0, 1) is degenerate. + mesh.add_vertex({0, 1, 0}); + mesh.add_triangle(0, 1, 2); + + // Split point on the degenerate edge (0, 1): interpolation weight is undefined. + Index split_vertex = mesh.get_num_vertices(); + mesh.add_vertex({0, 0, 0}); + + mesh.initialize_edges(); + Index split_eid = mesh.find_edge_from_vertices(0, 1); + std::vector split_pts = {split_vertex}; + auto get_edge_split_pts = [&](Index eid) -> span { + if (eid == split_eid) return span(split_pts.data(), split_pts.size()); + return span(split_pts.data(), 0); + }; + + LA_REQUIRE_THROWS( + internal::split_edges_only(mesh, get_edge_split_pts, always_active)); + } +} + +TEST_CASE("internal::split_edges shared split vertex", "[core]") +{ + using Scalar = double; + using Index = uint32_t; + using namespace lagrange; + + auto always_active = [](Index) { return true; }; + + // Two triangles with overlapping collinear edges on the x-axis, both carrying the same split + // vertex `c`; each triangle must interpolate `c` along its own edge, not a shared parent edge. + const Index a0 = 0, a1 = 1, a_apex = 2; // triangle A: edge (a0, a1) spans x = 0..4 + const Index b0 = 3, b1 = 4, b_apex = 5; // triangle B: edge (b0, b1) spans x = 1..3 + const Index c = 6; // midpoint of both edges + + auto make_mesh = [&]() { + SurfaceMesh mesh; + mesh.add_vertex({0, 0, 0}); + mesh.add_vertex({4, 0, 0}); + mesh.add_vertex({2, 2, 0}); + mesh.add_vertex({1, 0, 0}); + mesh.add_vertex({3, 0, 0}); + mesh.add_vertex({2, -2, 0}); + mesh.add_vertex({2, 0, 0}); + mesh.add_triangle(a0, a1, a_apex); + mesh.add_triangle(b0, b1, b_apex); + mesh.create_attribute( + "corner_value", + AttributeElement::Corner, + AttributeUsage::Scalar); + attribute_vector_ref(mesh, "corner_value") << 10, 40, 100, 0, 60, 200; + + auto indexed_id = mesh.create_attribute( + "indexed_value", + AttributeElement::Indexed, + AttributeUsage::Scalar, + 1); + { + auto& iattr = mesh.ref_indexed_attribute(indexed_id); + iattr.values().resize_elements(6); + auto values = iattr.values().ref_all(); + values[0] = 10, values[1] = 40, values[2] = 100, values[3] = 0, values[4] = 60, + values[5] = 200; + auto indices = iattr.indices().ref_all(); + for (Index i = 0; i < 6; ++i) indices[i] = i; + } + return mesh; + }; + + // Run `split_fn` on a fresh mesh with `c` declared as a split point on both collinear edges. + auto split = [&](auto&& split_fn) { + auto mesh = make_mesh(); + mesh.initialize_edges(); + const Index ea = mesh.find_edge_from_vertices(a0, a1); + const Index eb = mesh.find_edge_from_vertices(b0, b1); + std::vector split_pts = {c}; + auto get_edge_split_pts = [&](Index eid) -> span { + if (eid == ea || eid == eb) return span(split_pts.data(), split_pts.size()); + return span(split_pts.data(), 0); + }; + split_fn(mesh, get_edge_split_pts, always_active); + return mesh; + }; + + // In triangle A, c is the midpoint of (a0, a1) -> 0.5 * 10 + 0.5 * 40 = 25. + // In triangle B, c is the midpoint of (b0, b1) -> 0.5 * 0 + 0.5 * 60 = 30. + auto verify = [&](SurfaceMesh& mesh) { + const Scalar expected_a = 25; + const Scalar expected_b = 30; + auto corner_value = attribute_vector_ref(mesh, "corner_value"); + auto& iattr = mesh.get_indexed_attribute("indexed_value"); + auto indexed_values = iattr.values().get_all(); + auto indexed_indices = iattr.indices().get_all(); + Index checked_a = 0; + Index checked_b = 0; + for (Index f = 2; f < mesh.get_num_facets(); ++f) { + auto fv = mesh.get_facet_vertices(f); + const Index cb = mesh.get_facet_corner_begin(f); + bool from_a = false; + bool from_b = false; + for (Index k = 0; k < mesh.get_facet_size(f); ++k) { + if (fv[k] == a0 || fv[k] == a1 || fv[k] == a_apex) from_a = true; + if (fv[k] == b0 || fv[k] == b1 || fv[k] == b_apex) from_b = true; + } + REQUIRE(from_a != from_b); + for (Index k = 0; k < mesh.get_facet_size(f); ++k) { + if (fv[k] != c) continue; + CAPTURE(f, from_a); + const Scalar expected = from_a ? expected_a : expected_b; + CHECK(corner_value[cb + k] == expected); + CHECK(indexed_values[indexed_indices[cb + k]] == expected); + (from_a ? checked_a : checked_b) += 1; + } + } + REQUIRE(checked_a > 0); + REQUIRE(checked_b > 0); + }; + + SECTION("split_edges") + { + auto mesh = split([](SurfaceMesh& m, auto&& get_pts, auto&& active) { + internal::split_edges(m, get_pts, active); + }); + verify(mesh); + } + + SECTION("split_edges_only") + { + auto mesh = split([](SurfaceMesh& m, auto&& get_pts, auto&& active) { + internal::split_edges_only(m, get_pts, active); + }); + verify(mesh); + } +} diff --git a/modules/core/tests/test_transform_mesh.cpp b/modules/core/tests/test_transform_mesh.cpp index 5344cb75..d6bb4308 100644 --- a/modules/core/tests/test_transform_mesh.cpp +++ b/modules/core/tests/test_transform_mesh.cpp @@ -281,6 +281,91 @@ void test_transform_mesh_3d(bool pad_with_sign, TestCase test_case) } } +// Regression test for normalizing Normal/Tangent attributes on a 2D mesh. The transform used to +// normalize a hardcoded head<3>, which (a) reads past the end of a 2-channel (== dim) row on a 2D +// mesh (out-of-bounds / Eigen assert), and (b) would normalize the extra channel of a 3-channel +// (== dim+1) attribute. Normalization must operate on exactly the Dimension geometric components. +void test_transform_mesh_2d_attribute_normalization() +{ + using Scalar = double; + using Index = uint32_t; + + lagrange::SurfaceMesh mesh(2); + mesh.add_vertex({0, 0}); + mesh.add_vertex({1, 0}); + mesh.add_vertex({0, 1}); + mesh.add_triangle(0, 1, 2); + + // dim (2-channel) Normal/Tangent attributes with deliberately non-unit values so normalization + // is observable: each normal row is (3, 4) (norm 5), each tangent row is (0, 2). + auto id_nrm = mesh.create_attribute( + "normal", + lagrange::AttributeElement::Vertex, + lagrange::AttributeUsage::Normal, + 2, + std::array{3., 4., 3., 4., 3., 4.}); + auto id_tan = mesh.create_attribute( + "tangent", + lagrange::AttributeElement::Vertex, + lagrange::AttributeUsage::Tangent, + 2, + std::array{0., 2., 0., 2., 0., 2.}); + + // dim+1 (3-channel) Normal/Tangent attributes: the first 2 channels are the geometric vector + // (3, 4); the 3rd is an extra channel (e.g. a sign/padding) set to 7 that must be preserved + // verbatim by both the transform (leftCols) and the normalization (head). + auto id_nrm3 = mesh.create_attribute( + "normal3", + lagrange::AttributeElement::Vertex, + lagrange::AttributeUsage::Normal, + 3, + std::array{3., 4., 7., 3., 4., 7., 3., 4., 7.}); + auto id_tan3 = mesh.create_attribute( + "tangent3", + lagrange::AttributeElement::Vertex, + lagrange::AttributeUsage::Tangent, + 3, + std::array{0., 2., 7., 0., 2., 7., 0., 2., 7.}); + + lagrange::TransformOptions opt; + opt.normalize_normals = true; + opt.normalize_tangents_bitangents = true; + + // Identity transform: leaves directions unchanged, so we exercise (and isolate) the + // normalization path. With the previous head<3> code the 2-channel case would read out of + // bounds and the 3-channel case would normalize across the extra channel. + lagrange::transform_mesh(mesh, Eigen::Affine2d::Identity(), opt); + + auto nrm = lagrange::matrix_view(mesh.get_attribute(id_nrm)); + auto tan = lagrange::matrix_view(mesh.get_attribute(id_tan)); + auto nrm3 = lagrange::matrix_view(mesh.get_attribute(id_nrm3)); + auto tan3 = lagrange::matrix_view(mesh.get_attribute(id_tan3)); + REQUIRE(nrm.cols() == 2); + REQUIRE(tan.cols() == 2); + REQUIRE(nrm3.cols() == 3); + REQUIRE(tan3.cols() == 3); + for (Index v = 0; v < 3; ++v) { + // dim case: the 2-vector is normalized. + REQUIRE_THAT(nrm.row(v).norm(), Catch::Matchers::WithinAbs(1.0, 1e-12)); + REQUIRE_THAT(nrm.row(v)(0), Catch::Matchers::WithinAbs(0.6, 1e-12)); + REQUIRE_THAT(nrm.row(v)(1), Catch::Matchers::WithinAbs(0.8, 1e-12)); + REQUIRE_THAT(tan.row(v).norm(), Catch::Matchers::WithinAbs(1.0, 1e-12)); + REQUIRE_THAT(tan.row(v)(0), Catch::Matchers::WithinAbs(0.0, 1e-12)); + REQUIRE_THAT(tan.row(v)(1), Catch::Matchers::WithinAbs(1.0, 1e-12)); + + // dim+1 case: only the first 2 (geometric) channels are normalized; the extra 3rd channel + // is left untouched at 7. + REQUIRE_THAT(nrm3.row(v).head<2>().norm(), Catch::Matchers::WithinAbs(1.0, 1e-12)); + REQUIRE_THAT(nrm3.row(v)(0), Catch::Matchers::WithinAbs(0.6, 1e-12)); + REQUIRE_THAT(nrm3.row(v)(1), Catch::Matchers::WithinAbs(0.8, 1e-12)); + REQUIRE_THAT(nrm3.row(v)(2), Catch::Matchers::WithinAbs(7.0, 1e-12)); + REQUIRE_THAT(tan3.row(v).head<2>().norm(), Catch::Matchers::WithinAbs(1.0, 1e-12)); + REQUIRE_THAT(tan3.row(v)(0), Catch::Matchers::WithinAbs(0.0, 1e-12)); + REQUIRE_THAT(tan3.row(v)(1), Catch::Matchers::WithinAbs(1.0, 1e-12)); + REQUIRE_THAT(tan3.row(v)(2), Catch::Matchers::WithinAbs(7.0, 1e-12)); + } +} + } // namespace TEST_CASE("transform_mesh_2d", "[next]") @@ -290,6 +375,11 @@ TEST_CASE("transform_mesh_2d", "[next]") } } +TEST_CASE("transform_mesh_2d_attribute_normalization", "[next]") +{ + test_transform_mesh_2d_attribute_normalization(); +} + TEST_CASE("transform_mesh_3d", "[next]") { for (int i = 0; i < static_cast(TestCase::NumTestCases); ++i) { diff --git a/modules/core/tests/test_triangulate_polygonal_facets.cpp b/modules/core/tests/test_triangulate_polygonal_facets.cpp index b4ef025e..7672a900 100644 --- a/modules/core/tests/test_triangulate_polygonal_facets.cpp +++ b/modules/core/tests/test_triangulate_polygonal_facets.cpp @@ -360,6 +360,132 @@ void test_centroid_fan() } } +template +void test_should_triangulate() +{ + using namespace lagrange; + + // Two disjoint (convex, planar) pentagons: both schemes honor `should_triangulate` for every + // facet with more than 3 vertices, so selecting one pentagon exercises the predicate. + auto make_mesh = [] { + SurfaceMesh mesh(3); + // Pentagon 0. + mesh.add_vertex({0, 0, 0}); // 0 + mesh.add_vertex({2, 0, 0}); // 1 + mesh.add_vertex({Scalar(2.5), Scalar(1.5), 0}); // 2 + mesh.add_vertex({1, Scalar(2.5), 0}); // 3 + mesh.add_vertex({Scalar(-0.5), Scalar(1.5), 0}); // 4 + // Pentagon 1. + mesh.add_vertex({5, 0, 0}); // 5 + mesh.add_vertex({7, 0, 0}); // 6 + mesh.add_vertex({Scalar(7.5), Scalar(1.5), 0}); // 7 + mesh.add_vertex({6, Scalar(2.5), 0}); // 8 + mesh.add_vertex({Scalar(4.5), Scalar(1.5), 0}); // 9 + mesh.add_polygon({0, 1, 2, 3, 4}); // facet 0 (pentagon) + mesh.add_polygon({5, 6, 7, 8, 9}); // facet 1 (pentagon) + return mesh; + }; + + for (auto scheme : + {TriangulationOptions::Scheme::Earcut, TriangulationOptions::Scheme::CentroidFan}) { + TriangulationOptions options; + options.scheme = scheme; + + // Only triangulate the second facet. The first pentagon is left untouched. + auto mesh = make_mesh(); + const Index target_facet = 1; + triangulate_polygonal_facets( + mesh, + function_ref([&](Index f) { return f == target_facet; }), + options); + + // Mesh is not fully triangulated since we skipped one facet. + REQUIRE(!mesh.is_triangle_mesh()); + + // Count facets by size: the untouched pentagon (5) must survive, and the triangulated + // pentagon must have become triangles. + Index num_tris = 0; + Index num_pentagons = 0; + for (Index f = 0; f < mesh.get_num_facets(); ++f) { + switch (mesh.get_facet_size(f)) { + case 3: ++num_tris; break; + case 5: ++num_pentagons; break; + default: break; + } + } + REQUIRE(num_pentagons == 1); + // A pentagon triangulates into 3 triangles (earcut) or 5 triangles (centroid fan). + REQUIRE(num_tris >= 3); + } + + // Quads must also honor the predicate. Regression test: the earcut scheme previously + // triangulated every quad unconditionally, ignoring `should_triangulate`. + for (auto scheme : + {TriangulationOptions::Scheme::Earcut, TriangulationOptions::Scheme::CentroidFan}) { + TriangulationOptions options; + options.scheme = scheme; + + SurfaceMesh mesh(3); + mesh.add_vertex({0, 0, 0}); + mesh.add_vertex({1, 0, 0}); + mesh.add_vertex({1, 1, 0}); + mesh.add_vertex({0, 1, 0}); + mesh.add_vertex({2, 0, 0}); + mesh.add_vertex({2, 1, 0}); + mesh.add_quad(0, 1, 2, 3); // facet 0 + mesh.add_quad(1, 4, 5, 2); // facet 1 + + const Index target_facet = 1; + triangulate_polygonal_facets( + mesh, + function_ref([&](Index f) { return f == target_facet; }), + options); + + // Facet 0 stays a quad; facet 1 is triangulated. + Index num_tris = 0; + Index num_quads = 0; + for (Index f = 0; f < mesh.get_num_facets(); ++f) { + switch (mesh.get_facet_size(f)) { + case 3: ++num_tris; break; + case 4: ++num_quads; break; + default: break; + } + } + REQUIRE(num_quads == 1); + REQUIRE(num_tris >= 2); + } + + // A predicate that always returns false is a no-op for both schemes, even with edge + // connectivity (an empty selection must not add empty facet buffers, which asserts). + for (auto scheme : + {TriangulationOptions::Scheme::Earcut, TriangulationOptions::Scheme::CentroidFan}) { + TriangulationOptions options; + options.scheme = scheme; + + auto mesh = make_mesh(); + mesh.initialize_edges(); + const Index old_num_facets = mesh.get_num_facets(); + triangulate_polygonal_facets( + mesh, + function_ref([](Index) { return false; }), + options); + REQUIRE(!mesh.is_triangle_mesh()); + REQUIRE(mesh.get_num_facets() == old_num_facets); + } + + // A predicate that always returns true is equivalent to the default behavior. + { + TriangulationOptions options; + auto mesh = make_mesh(); + triangulate_polygonal_facets( + mesh, + function_ref([](Index) { return true; }), + options); + mesh.compress_if_regular(); + REQUIRE(mesh.is_triangle_mesh()); + } +} + } // namespace TEST_CASE("earcut", "[core]") @@ -413,5 +539,11 @@ TEST_CASE("triangulate_polygonal_facets: centroid fan", "[core]") LA_SURFACE_MESH_X(centroid_fan, 0) } +TEST_CASE("triangulate_polygonal_facets: should_triangulate", "[core]") +{ +#define LA_X_should_triangulate(_, Scalar, Index) test_should_triangulate(); + LA_SURFACE_MESH_X(should_triangulate, 0) +} + // TODO: Test removal degenerate facets, once we allow sizes <= 2 // TODO: Test with 2d meshes diff --git a/modules/image_io/include/lagrange/image_io/load_image.h b/modules/image_io/include/lagrange/image_io/load_image.h index 60a35190..b69d6338 100644 --- a/modules/image_io/include/lagrange/image_io/load_image.h +++ b/modules/image_io/include/lagrange/image_io/load_image.h @@ -38,6 +38,12 @@ struct LoadImageResult LA_IMAGE_IO_API LoadImageResult load_image(const fs::path& path, spdlog::level::level_enum error_lvl = spdlog::level::err); +// Load image from a memory buffer (PNG, JPEG, or other stb-supported formats). +LA_IMAGE_IO_API LoadImageResult load_image_from_buffer( + const void* buffer, + size_t size, + spdlog::level::level_enum error_lvl = spdlog::level::err); + // Load png or jpg image using stb library. Produces uint8 data. LA_IMAGE_IO_API LoadImageResult load_image_stb(const fs::path& path, spdlog::level::level_enum error_lvl = spdlog::level::err); diff --git a/modules/image_io/src/load_image.cpp b/modules/image_io/src/load_image.cpp index 9304f0ca..afa43b9f 100644 --- a/modules/image_io/src/load_image.cpp +++ b/modules/image_io/src/load_image.cpp @@ -19,6 +19,8 @@ #include +#include + namespace lagrange { namespace image_io { @@ -247,5 +249,84 @@ LoadImageResult load_image_bin(const fs::path& path, spdlog::level::level_enum e return rtn; } +LoadImageResult +load_image_from_buffer(const void* buffer, size_t size, spdlog::level::level_enum error_lvl) +{ + LoadImageResult rtn; + if (buffer == nullptr || size == 0) { + logger().log(error_lvl, "load_image_from_buffer error: empty buffer"); + return rtn; + } + if (size > static_cast(std::numeric_limits::max())) { + logger().log(error_lvl, "load_image_from_buffer error: buffer too large"); + return rtn; + } + + const auto* buf = reinterpret_cast(buffer); + const int len = static_cast(size); + int w, h, ch; + + if (stbi_is_16_bit_from_memory(buf, len)) { + rtn.precision = image::ImagePrecision::uint16; + uint16_t* data = stbi_load_16_from_memory(buf, len, &w, &h, &ch, STBI_default); + if (data == nullptr) { + logger().log(error_lvl, "load_image_from_buffer error: stbi failed to decode image"); + return rtn; + } + if (ch != 1 && ch != 3 && ch != 4) { + logger().log( + error_lvl, + "load_image_from_buffer error: unsupported channel count {}", + ch); + stbi_image_free(data); + return rtn; + } + if (w <= 0 || h <= 0) { + stbi_image_free(data); + return rtn; + } + size_t _w = static_cast(w); + size_t _h = static_cast(h); + size_t _ch = static_cast(ch); + rtn.valid = true; + rtn.width = _w; + rtn.height = _h; + rtn.channel = static_cast(ch); + rtn.storage = std::make_shared(sizeof(uint16_t) * _ch * _w, _h, 1); + std::copy_n(data, _ch * _w * _h, reinterpret_cast(rtn.storage->data())); + stbi_image_free(data); + } else { + rtn.precision = image::ImagePrecision::uint8; + unsigned char* data = stbi_load_from_memory(buf, len, &w, &h, &ch, STBI_default); + if (data == nullptr) { + logger().log(error_lvl, "load_image_from_buffer error: stbi failed to decode image"); + return rtn; + } + if (ch != 1 && ch != 3 && ch != 4) { + logger().log( + error_lvl, + "load_image_from_buffer error: unsupported channel count {}", + ch); + stbi_image_free(data); + return rtn; + } + if (w <= 0 || h <= 0) { + stbi_image_free(data); + return rtn; + } + size_t _w = static_cast(w); + size_t _h = static_cast(h); + size_t _ch = static_cast(ch); + rtn.valid = true; + rtn.width = _w; + rtn.height = _h; + rtn.channel = static_cast(ch); + rtn.storage = std::make_shared(_ch * _w, _h, 1); + std::copy_n(data, _ch * _w * _h, rtn.storage->data()); + stbi_image_free(data); + } + return rtn; +} + } // namespace image_io } // namespace lagrange diff --git a/modules/io/include/lagrange/io/internal/scene_utils.h b/modules/io/include/lagrange/io/internal/scene_utils.h index c1c822ac..b5eee4c0 100644 --- a/modules/io/include/lagrange/io/internal/scene_utils.h +++ b/modules/io/include/lagrange/io/internal/scene_utils.h @@ -30,4 +30,20 @@ bool try_load_image( const LoadOptions& options, scene::ImageExperimental& image); +/** + * Load an image from a memory buffer (PNG, JPEG, or other stb-supported formats). + * + * @param[in] buffer Pointer to raw image file bytes. + * @param[in] size Number of bytes in the buffer. + * @param[in] options Load options. + * @param[out] image This will be filled with the loaded data. + * + * @return true if successful. + */ +bool try_load_image_from_buffer( + const void* buffer, + size_t size, + const LoadOptions& options, + scene::ImageExperimental& image); + } // namespace lagrange::io::internal diff --git a/modules/io/src/internal/scene_utils.cpp b/modules/io/src/internal/scene_utils.cpp index a554723d..174dcd34 100644 --- a/modules/io/src/internal/scene_utils.cpp +++ b/modules/io/src/internal/scene_utils.cpp @@ -16,6 +16,35 @@ namespace lagrange::io::internal { +namespace { + +void populate_image_buffer( + const image_io::LoadImageResult& result, + scene::ImageBufferExperimental& buf) +{ + buf.width = result.width; + buf.height = result.height; + buf.num_channels = static_cast(result.channel); + switch (result.precision) { + case image::ImagePrecision::uint8: buf.element_type = AttributeValueType::e_uint8_t; break; + case image::ImagePrecision::int8: buf.element_type = AttributeValueType::e_int8_t; break; + case image::ImagePrecision::uint16: buf.element_type = AttributeValueType::e_uint16_t; break; + case image::ImagePrecision::uint32: buf.element_type = AttributeValueType::e_uint32_t; break; + case image::ImagePrecision::int32: buf.element_type = AttributeValueType::e_int32_t; break; + case image::ImagePrecision::float32: buf.element_type = AttributeValueType::e_float; break; + case image::ImagePrecision::float64: buf.element_type = AttributeValueType::e_double; break; + case image::ImagePrecision::float16: [[fallthrough]]; + default: throw std::runtime_error("Unsupported image precision"); + } + + la_runtime_assert(result.storage != nullptr); + const size_t num_bytes = + buf.width * buf.height * buf.num_channels * buf.get_bits_per_element() / 8; + buf.data.assign(result.storage->data(), result.storage->data() + num_bytes); +} + +} // namespace + bool try_load_image( const std::string& name, const LoadOptions& options, @@ -31,30 +60,23 @@ bool try_load_image( image_io::LoadImageResult result = image_io::load_image(path, error_lvl); if (!result.valid) return false; - scene::ImageBufferExperimental& buffer = image.image; - buffer.width = result.width; - buffer.height = result.height; - buffer.num_channels = static_cast(result.channel); - switch (result.precision) { - case image::ImagePrecision::uint8: buffer.element_type = AttributeValueType::e_uint8_t; break; - case image::ImagePrecision::int8: buffer.element_type = AttributeValueType::e_int8_t; break; - case image::ImagePrecision::uint32: buffer.element_type = AttributeValueType::e_uint32_t; break; - case image::ImagePrecision::int32: buffer.element_type = AttributeValueType::e_int32_t; break; - case image::ImagePrecision::float32: buffer.element_type = AttributeValueType::e_float; break; - case image::ImagePrecision::float64: buffer.element_type = AttributeValueType::e_double; break; - case image::ImagePrecision::float16: [[fallthrough]]; - default: throw std::runtime_error("Unsupported image precision"); - } + populate_image_buffer(result, image.image); + return true; +} - la_runtime_assert(result.storage != nullptr); - const size_t num_bytes = - buffer.width * buffer.height * buffer.num_channels * buffer.get_bits_per_element() / 8; - buffer.data.reserve(num_bytes); - std::copy( - result.storage->data(), - result.storage->data() + num_bytes, - std::back_inserter(buffer.data)); +bool try_load_image_from_buffer( + const void* buffer, + size_t size, + const LoadOptions& options, + scene::ImageExperimental& image) +{ + spdlog::level::level_enum error_lvl = spdlog::level::err; + if (options.quiet) error_lvl = spdlog::level::off; + + image_io::LoadImageResult result = image_io::load_image_from_buffer(buffer, size, error_lvl); + if (!result.valid) return false; + populate_image_buffer(result, image.image); return true; } diff --git a/modules/io/src/load_fbx.cpp b/modules/io/src/load_fbx.cpp index 9c8d7685..6e95fe7c 100644 --- a/modules/io/src/load_fbx.cpp +++ b/modules/io/src/load_fbx.cpp @@ -247,6 +247,25 @@ MeshType convert_mesh_ufbx_to_lagrange(const ufbx_mesh* mesh, const LoadOptions& attr.values().ref_all().begin()); } + // Per-facet material slot index from ufbx_mesh::face_material (UFBX_NO_INDEX → -1). + // The values index into node->materials[] for the primary instancing node, which is + // dense (no null slots) since create_load_opts() leaves connect_broken_elements off. + // Written before triangulation so triangulate_polygonal_facets replicates it correctly. + if (opt.load_materials && mesh->face_material.count == mesh->num_faces && mesh->num_faces > 0) { + auto attr_id = lmesh.template create_attribute( + AttributeName::material_id, + AttributeElement::Facet, + AttributeUsage::Scalar, + 1); + auto& attr = lmesh.template ref_attribute(attr_id); + attr.resize_elements(mesh->num_faces); + auto vals = attr.ref_all(); + for (size_t f = 0; f < mesh->num_faces; ++f) { + const uint32_t raw = mesh->face_material.data[f]; + vals[f] = (raw == UFBX_NO_INDEX) ? invalid() : safe_cast(raw); + } + } + if (opt.stitch_vertices) { stitch_mesh(lmesh); } @@ -482,12 +501,17 @@ SceneType load_scene_fbx(const ufbx_scene* scene, const LoadOptions& opt) limage.uri = texture->relative_filename.data; // note: there is no width/height anywhere. read the image from disk, or read png from texture->content. bool loaded = false; - if (texture->content.size > 0) { - // TODO: read image from embedded data. - // But our image_io module does not support loading from buffer - logger().warn( - "Loading fbx embedded textures is currently unsupported, missing data for {}", - limage.name); + if (texture->content.size > 0 && opt.load_images) { + loaded = internal::try_load_image_from_buffer( + texture->content.data, + texture->content.size, + opt, + limage); + if (!loaded && !opt.quiet) { + logger().warn( + "Failed to load embedded texture image for texture '{}'", + limage.name); + } } else if (opt.load_images) { loaded |= internal::try_load_image(texture->filename.data, opt, limage); loaded |= internal::try_load_image(texture->relative_filename.data, opt, limage); @@ -555,10 +579,12 @@ SceneType load_scene_fbx(const ufbx_scene* scene, const LoadOptions& opt) la_runtime_assert(mesh_idx != lagrange::invalid()); std::vector material_idxs; + // node->materials[] is dense (no null slots) unless ufbx is loaded with + // connect_broken_elements, which create_load_opts() never enables. This keeps + // the per-instance list parallel to mesh->face_material indices. for (const ufbx_material* material : node->materials) { - if (material) { - material_idxs.push_back(element_index[material->element_id]); - } + la_runtime_assert(material != nullptr); + material_idxs.push_back(element_index[material->element_id]); } lnode.meshes.push_back({mesh_idx, material_idxs}); } diff --git a/modules/io/tests/test_fbx.cpp b/modules/io/tests/test_fbx.cpp index 5d8346fd..2b1ea2ba 100644 --- a/modules/io/tests/test_fbx.cpp +++ b/modules/io/tests/test_fbx.cpp @@ -9,17 +9,21 @@ * 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 + TEST_CASE("load_fbx", "[io][fbx]" LA_CORP_FLAG) { lagrange::io::LoadOptions options; @@ -64,6 +68,86 @@ TEST_CASE("load_fbx with duplicate attr", "[io][fbx]") options)); } +TEST_CASE("load_fbx face_material preserved", "[io][fbx][scene]" LA_CORP_FLAG) +{ + // Verify that ufbx's per-face material assignment is preserved as a material_id + // facet attribute on the loaded mesh, not silently dropped. + lagrange::io::LoadOptions options; + options.quiet = true; + auto scene = lagrange::io::load_scene_fbx( + lagrange::testing::get_data_path( + "corp/io/human_alloy_failed_00036_Esmee003_Basics_001.fbx"), + options); + + // At least one raw mesh must carry the per-facet material_id attribute written + // by the loader from ufbx::face_material. + bool found_attr = false; + for (const auto& mesh : scene.meshes) { + if (mesh.has_attribute(lagrange::AttributeName::material_id)) { + found_attr = true; + break; + } + } + REQUIRE(found_attr); + + // After scene_to_meshes_and_materials, at least one output mesh must have more than + // one distinct per-instance material index — confirming the per-face data was loaded + // instead of a flat materials.front() assignment. + auto result = lagrange::scene::scene_to_meshes_and_materials(scene); + bool found_multi = false; + for (size_t mi = 0; mi < result.meshes.size(); ++mi) { + const auto& mesh = result.meshes[mi]; + if (!mesh.has_attribute(lagrange::AttributeName::material_id)) continue; + const auto& attr = mesh.template get_attribute( + lagrange::AttributeName::material_id); + const auto vals = attr.get_all(); + if (vals.size() < 2) continue; + // Count how many per-instance material indices this mesh references. + // invalid() means "no material". + // Also verify all valid IDs are in-range for this mesh's material list. + const auto num_mats = + static_cast(result.material_ids[mi].size()); + const auto inv = lagrange::invalid(); + bool all_in_range = true; + lagrange::SurfaceMesh32d::Index first_valid = inv; + bool has_second = false; + for (auto id : vals) { + if (id == inv) continue; + if (id >= num_mats) { + all_in_range = false; + break; + } + if (first_valid == inv) { + first_valid = id; + } else if (id != first_valid) { + has_second = true; + } + } + if (all_in_range && has_second) { + found_multi = true; + break; + } + } + REQUIRE(found_multi); +} + +TEST_CASE("load_fbx embedded textures", "[io][fbx]") +{ + lagrange::io::LoadOptions options; + options.load_images = true; + options.quiet = true; + + auto scene = lagrange::io::load_scene_fbx( + lagrange::testing::get_data_path("open/io/Walking.fbx"), + options); + + // Walking.fbx has embedded textures — at least one image must have decoded pixel data. + bool any_loaded = std::any_of(scene.images.begin(), scene.images.end(), [](const auto& img) { + return img.image.width > 0 && img.image.height > 0 && !img.image.data.empty(); + }); + REQUIRE(any_loaded); +} + TEST_CASE("load_fbx with geometric transform", "[io][fbx]" LA_CORP_FLAG) { lagrange::io::LoadOptions options; @@ -85,3 +169,54 @@ TEST_CASE("load_fbx with geometric transform", "[io][fbx]" LA_CORP_FLAG) lagrange::testing::ensure_approx_equivalent_mesh(mesh, expected); lagrange::testing::ensure_approx_equivalent_mesh(flat, expected); } + +TEST_CASE( + "scene_to_meshes_and_materials mesh/material count match", + "[io][fbx][scene]" LA_CORP_FLAG) +{ + lagrange::io::LoadOptions options; + options.quiet = true; + auto scene = lagrange::io::load_scene_fbx( + lagrange::testing::get_data_path( + "corp/io/human_alloy_failed_00036_Esmee003_Basics_001.fbx"), + options); + + auto result = lagrange::scene::scene_to_meshes_and_materials(scene); + + // meshes[i] and material_ids[i] must be in 1-to-1 correspondence. + REQUIRE(result.meshes.size() == result.material_ids.size()); + + for (size_t i = 0; i < result.meshes.size(); ++i) { + const auto& mesh = result.meshes[i]; + const auto& mat_ids = result.material_ids[i]; + + // Every material index in the per-mesh list must be a valid scene material index. + for (auto mat_id : mat_ids) { + REQUIRE(mat_id < scene.materials.size()); + } + + // scene_to_meshes_and_materials must always produce a material_id facet attribute. + REQUIRE(mesh.has_attribute(lagrange::AttributeName::material_id)); + const auto& attr = mesh.template get_attribute( + lagrange::AttributeName::material_id); + REQUIRE(attr.get_num_elements() == mesh.get_num_facets()); + + for (auto facet_mat_id : attr.get_all()) { + const bool is_invalid = + facet_mat_id == lagrange::invalid(); + // Per-facet value is either a valid scene material index or the no-material sentinel. + REQUIRE((is_invalid || facet_mat_id < scene.materials.size())); + if (!is_invalid) { + // It must also appear in the per-mesh material_ids list. + bool found = false; + for (auto id : mat_ids) { + if (id == facet_mat_id) { + found = true; + break; + } + } + REQUIRE(found); + } + } + } +} diff --git a/modules/polyddg/include/lagrange/polyddg/compute_smooth_direction_field.h b/modules/polyddg/include/lagrange/polyddg/compute_smooth_direction_field.h index b2254316..a416f6d3 100644 --- a/modules/polyddg/include/lagrange/polyddg/compute_smooth_direction_field.h +++ b/modules/polyddg/include/lagrange/polyddg/compute_smooth_direction_field.h @@ -63,6 +63,17 @@ struct SmoothDirectionFieldOptions /// - @c AttributeElement::Vertex: @c \@smooth_direction_field /// - @c AttributeElement::Facet: @c \@smooth_direction_field_facets std::optional direction_field_attribute; + + /// Alignment tradeoff parameter λ_t balancing smoothness against alignment strength in the + /// constrained solve (L_reg − λ_t M) u = M q (Knöppel et al. 2013, Eq. 16 / Algorithm 3). + /// Only used when @c alignment_attribute is set. + /// + /// - @c 0 (default): the paper's recommended balanced value. + /// - Negative (toward −∞): stronger alignment, weaker smoothness. + /// - Positive (toward λ₁, the smallest generalized eigenvalue): weaker alignment, smoother + /// field. Must stay strictly below λ₁; values at or above it make the system indefinite + /// and the solve fails. + double alignment_lambda = 0.0; }; /// diff --git a/modules/polyddg/python/src/polyddg.cpp b/modules/polyddg/python/src/polyddg.cpp index 28c16ab4..f98aba86 100644 --- a/modules/polyddg/python/src/polyddg.cpp +++ b/modules/polyddg/python/src/polyddg.cpp @@ -788,7 +788,8 @@ are the principal directions. All four quantities are stored as vertex attribute double beta, lagrange::AttributeElement output_element_type, std::string_view alignment_attribute, - std::string_view direction_field_attribute) { + std::string_view direction_field_attribute, + double alignment_lambda) { polyddg::SmoothDirectionFieldOptions opts; opts.nrosy = nrosy; opts.lambda = beta; @@ -796,6 +797,7 @@ are the principal directions. All four quantities are stored as vertex attribute opts.alignment_attribute = alignment_attribute; if (!direction_field_attribute.empty()) opts.direction_field_attribute = direction_field_attribute; + opts.alignment_lambda = alignment_lambda; return polyddg::compute_smooth_direction_field(mesh, ops, opts); }, "mesh"_a, @@ -806,6 +808,7 @@ are the principal directions. All four quantities are stored as vertex attribute "output_element_type"_a = default_sdf_opts.output_element_type, "alignment_attribute"_a = default_sdf_opts.alignment_attribute, "direction_field_attribute"_a = "", + "alignment_lambda"_a = default_sdf_opts.alignment_lambda, R"(Compute the globally smoothest n-direction field on a surface mesh. Dispatches to a vertex-based or facet-based implementation depending on @@ -839,6 +842,11 @@ per-vertex for ``Vertex``, per-facet for ``Facet``. :param direction_field_attribute: Output attribute name. Pass ``""`` (the default) to use the canonical name (``"@smooth_direction_field"`` for Vertex, ``"@smooth_direction_field_facets"`` for Facet). +:param alignment_lambda: Tradeoff parameter :math:`\lambda_t` balancing smoothness against + alignment strength in the constrained solve (Knöppel et al. 2013, Eq. 16). Only used when + ``alignment_attribute`` is set. ``0`` (default) is the recommended balanced value; negative + values strengthen alignment; positive values (below the smallest generalized eigenvalue) + yield a smoother, less-aligned field. :return: Attribute ID of the output direction field.)"); @@ -849,7 +857,8 @@ per-vertex for ``Vertex``, per-facet for ``Facet``. double beta, lagrange::AttributeElement output_element_type, std::string_view alignment_attribute, - std::string_view direction_field_attribute) { + std::string_view direction_field_attribute, + double alignment_lambda) { polyddg::SmoothDirectionFieldOptions opts; opts.nrosy = nrosy; opts.lambda = beta; @@ -857,6 +866,7 @@ per-vertex for ``Vertex``, per-facet for ``Facet``. opts.alignment_attribute = alignment_attribute; if (!direction_field_attribute.empty()) opts.direction_field_attribute = direction_field_attribute; + opts.alignment_lambda = alignment_lambda; return polyddg::compute_smooth_direction_field(mesh, opts); }, "mesh"_a, @@ -866,6 +876,7 @@ per-vertex for ``Vertex``, per-facet for ``Facet``. "output_element_type"_a = default_sdf_opts.output_element_type, "alignment_attribute"_a = default_sdf_opts.alignment_attribute, "direction_field_attribute"_a = "", + "alignment_lambda"_a = default_sdf_opts.alignment_lambda, R"(Compute the globally smoothest n-direction field on a surface mesh. Convenience overload that constructs a :class:`DifferentialOperators` instance internally. @@ -887,6 +898,11 @@ Dispatches to a vertex-based or facet-based implementation depending on (default) for per-vertex output, or ``AttributeElement.Facet`` for per-facet output. :param alignment_attribute: Name of an alignment vector attribute (zero = unconstrained). Must match ``output_element_type``. If empty, the unconstrained smoothest field is computed. +:param alignment_lambda: Tradeoff parameter :math:`\lambda_t` balancing smoothness against + alignment strength in the constrained solve (Knöppel et al. 2013, Eq. 16). Only used when + ``alignment_attribute`` is set. ``0`` (default) is the recommended balanced value; negative + values strengthen alignment; positive values (below the smallest generalized eigenvalue) + yield a smoother, less-aligned field. :param direction_field_attribute: Output attribute name. Pass ``""`` (the default) to use the canonical name. diff --git a/modules/polyddg/python/tests/test_polyddg.py b/modules/polyddg/python/tests/test_polyddg.py index c7ea3e97..74471edb 100644 --- a/modules/polyddg/python/tests/test_polyddg.py +++ b/modules/polyddg/python/tests/test_polyddg.py @@ -380,3 +380,55 @@ def test_determinism(self, octahedron): d2 = np.array(octahedron.attribute(id2).data) # Fields may differ by global sign flip; compare abs dot products. assert np.allclose(np.abs(d1), np.abs(d2), atol=1e-10) or np.allclose(d1, -d2, atol=1e-10) + + def test_alignment_lambda(self): + """alignment_lambda is accepted and modulates the smoothness/alignment tradeoff. + + On a flat grid (λ_1 ≈ 0), the λ_t = 0 solve is smoothness-dominated and ignores a + high-frequency guidance field, while a strongly negative λ_t recovers it. + """ + n = 12 + grid = lagrange.SurfaceMesh() + for i in range(n): + for j in range(n): + grid.add_vertex([float(i), float(j), 0.0]) + for i in range(n - 1): + for j in range(n - 1): + v00, v10 = i * n + j, (i + 1) * n + j + v01, v11 = i * n + (j + 1), (i + 1) * n + (j + 1) + grid.add_triangle(v00, v10, v11) + grid.add_triangle(v00, v11, v01) + ops = lagrange.polyddg.DifferentialOperators(grid) + nv = grid.num_vertices + + # Dense, high-frequency guidance: alternate the target angle by ±22.5° per vertex. + bases = np.stack([np.array(ops.vertex_basis(v)) for v in range(nv)]) # (nv, 3, 2) + guide_angle = np.where(np.arange(nv) % 2 == 0, 1.0, -1.0) * (np.pi / 8.0) + local = np.stack([np.cos(guide_angle), np.sin(guide_angle)], axis=1) # (nv, 2) + align = np.einsum("vij,vj->vi", bases, local) + grid.create_attribute( + "@align", + element=lagrange.AttributeElement.Vertex, + usage=lagrange.AttributeUsage.Vector, + initial_values=align, + ) + + def mean_align_error(lambda_t, out_name): + attr_id = lagrange.polyddg.compute_smooth_direction_field( + grid, + ops, + nrosy=4, + alignment_attribute="@align", + alignment_lambda=lambda_t, + direction_field_attribute=out_name, + ) + data = np.array(grid.attribute(attr_id).data).reshape(-1, 3) + assert np.allclose(np.linalg.norm(data, axis=1), 1.0, atol=1e-10) + out2 = np.einsum("vij,vi->vj", bases, data) # project to local frames + out_angle = np.arctan2(out2[:, 1], out2[:, 0]) + return np.mean(1.0 - np.cos(4.0 * (out_angle - guide_angle))) + + err_default = mean_align_error(0.0, "@sdf_lt0") + err_strong = mean_align_error(-100.0, "@sdf_lt_strong") + assert err_default > 0.1 + assert err_strong < 0.1 * err_default diff --git a/modules/polyddg/src/compute_smooth_direction_field.cpp b/modules/polyddg/src/compute_smooth_direction_field.cpp index ff654af4..6ea21690 100644 --- a/modules/polyddg/src/compute_smooth_direction_field.cpp +++ b/modules/polyddg/src/compute_smooth_direction_field.cpp @@ -35,29 +35,41 @@ using solver::SolverLDLT; // solve_connection_laplacian — shared solve kernel for both vertex and facet paths // ============================================================================= -// Solves for the smoothest (or alignment-constrained) n-rosy field given a -// regularized Laplacian L_reg = L + εM and mass matrix M. -// -// When has_constraints is true, solves (L_reg) x = M q and normalizes in M-norm. -// When has_constraints is false, finds the smallest generalized eigenvector of -// L_reg x = σ M x via Spectra, falling back to inverse power iteration. -// -// fn_name is used only in runtime-assert/warning messages. +/** + * Solves for the smoothest (or alignment-constrained) n-rosy field given a regularized + * Laplacian L_reg = L + εM and mass matrix M. + * + * has_constraints=true: solve (L_reg − λ_t M) x = M q (Knöppel et al. 2013, Eq. 16) and + * normalize in M-norm; lambda_t is the alignment tradeoff parameter. + * + * has_constraints=false: smallest generalized eigenvector of L_reg x = σ M x via Spectra, + * falling back to inverse power iteration (lambda_t is unused). + * + * fn_name is used only in runtime-assert/warning messages. + */ template static void solve_connection_laplacian( const Eigen::SparseMatrix& L_reg, const Eigen::SparseMatrix& M, const Eigen::Matrix& q, bool has_constraints, + Scalar lambda_t, Eigen::Index size, std::string_view fn_name, Eigen::Matrix& x) { if (has_constraints) { - SolverLDLT> slv(L_reg); + // Eq. 16 system matrix: (L + εM) − λ_t M. Only materialize the shift when λ_t ≠ 0. + Eigen::SparseMatrix A_shifted; + if (lambda_t != Scalar(0)) A_shifted = L_reg - lambda_t * M; + const Eigen::SparseMatrix& A = (lambda_t != Scalar(0)) ? A_shifted : L_reg; + SolverLDLT> slv(A); la_runtime_assert( slv.info() == Eigen::Success, - lagrange::format("{}: factorization of L + eps*M failed", fn_name)); + lagrange::format( + "{}: factorization of (L + eps*M - lambda_t*M) failed; alignment_lambda may be " + "too large (must stay below the smallest generalized eigenvalue)", + fn_name)); x = slv.solve((M * q).eval()); la_runtime_assert( @@ -119,6 +131,7 @@ AttributeId compute_smooth_direction_field_on_facets( const Index num_edges = mesh.get_num_edges(); const Index n = static_cast(options.nrosy); const int n_int = static_cast(options.nrosy); + const Scalar lambda_t = static_cast(options.alignment_lambda); // ---- 1. Edge-to-face mapping ---- // For each edge store the local vertex index of edge.v0 in each adjacent face. @@ -277,6 +290,7 @@ AttributeId compute_smooth_direction_field_on_facets( M, q, has_constraints, + lambda_t, static_cast(2 * num_facets), "compute_smooth_direction_field_on_facets", x); @@ -335,6 +349,7 @@ AttributeId compute_smooth_direction_field( const Index n = static_cast(options.nrosy); const int n_int = static_cast(options.nrosy); const Scalar lambda = static_cast(options.lambda); + const Scalar lambda_t = static_cast(options.alignment_lambda); // Build the connection Laplacian L of size (#V*2) x (#V*2). auto L = ops.connection_laplacian_nrosy(n, lambda); @@ -411,6 +426,7 @@ AttributeId compute_smooth_direction_field( M, q, has_constraints, + lambda_t, static_cast(num_vertices * 2), "compute_smooth_direction_field", x); diff --git a/modules/polyddg/tests/test_compute_smooth_direction_field.cpp b/modules/polyddg/tests/test_compute_smooth_direction_field.cpp index 964f4304..60a42a09 100644 --- a/modules/polyddg/tests/test_compute_smooth_direction_field.cpp +++ b/modules/polyddg/tests/test_compute_smooth_direction_field.cpp @@ -27,6 +27,8 @@ #include #include +#include +#include TEST_CASE("compute_smooth_direction_field", "[polyddg]") { @@ -342,6 +344,75 @@ TEST_CASE("compute_smooth_direction_field", "[polyddg]") } } + SECTION("flat grid: alignment_lambda modulates alignment strength") + { + // A flat grid has a zero-energy smooth mode (λ_1 ≈ 0), so λ_t = 0 is smoothness-dominated + // and ignores high-frequency guidance while a strongly negative λ_t recovers it. + constexpr Index N = 12; + SurfaceMesh grid; + for (Index i = 0; i < N; ++i) { + for (Index j = 0; j < N; ++j) { + grid.add_vertex({static_cast(i), static_cast(j), 0.0}); + } + } + for (Index i = 0; i + 1 < N; ++i) { + for (Index j = 0; j + 1 < N; ++j) { + const Index v00 = i * N + j, v10 = (i + 1) * N + j; + const Index v01 = i * N + (j + 1), v11 = (i + 1) * N + (j + 1); + grid.add_triangle(v00, v10, v11); + grid.add_triangle(v00, v11, v01); + } + } + polyddg::DifferentialOperators ops(grid); + + const Index nv = grid.get_num_vertices(); + + // Dense, high-frequency guidance: alternate the target angle by ±22.5° (in each local + // frame) with vertex parity. A smooth field cannot follow this alternation. + std::vector guide_angle(nv); + auto align_id = internal::find_or_create_attribute( + grid, + "@lambda_t_alignment", + AttributeElement::Vertex, + AttributeUsage::Vector, + 3, + internal::ResetToDefault::Yes); + auto align_data = attribute_matrix_ref(grid, align_id); + for (Index vid = 0; vid < nv; ++vid) { + const Scalar a = (vid % 2 == 0 ? 1.0 : -1.0) * (internal::pi / 8.0); + guide_angle[vid] = a; + align_data.row(vid) = + (ops.vertex_basis(vid) * Eigen::Matrix(std::cos(a), std::sin(a))) + .transpose(); + } + + auto mean_align_error = [&](Scalar lambda_t, std::string_view out_name) { + polyddg::SmoothDirectionFieldOptions opts; + opts.nrosy = 4; + opts.alignment_attribute = "@lambda_t_alignment"; + opts.alignment_lambda = lambda_t; + opts.direction_field_attribute = out_name; + auto result = polyddg::compute_smooth_direction_field(grid, ops, opts); + auto data = attribute_matrix_view(grid, result); + Scalar err = 0; + for (Index vid = 0; vid < nv; ++vid) { + REQUIRE_THAT(data.row(vid).norm(), Catch::Matchers::WithinAbs(1.0, 1e-10)); + Eigen::Matrix out_2d = + ops.vertex_basis(vid).transpose() * data.row(vid).transpose(); + Scalar out_angle = std::atan2(out_2d(1), out_2d(0)); + err += 1.0 - std::cos(4.0 * (out_angle - guide_angle[vid])); + } + return err / static_cast(nv); + }; + + // λ_t = 0 favors smoothness and cannot follow the alternating guidance; a strongly + // negative λ_t strengthens alignment and drives the mean error toward zero. + const Scalar err_default = mean_align_error(0.0, "@sdf_lt_default"); + const Scalar err_strong = mean_align_error(-100.0, "@sdf_lt_strong"); + REQUIRE(err_default > 0.1); + REQUIRE(err_strong < 0.1 * err_default); + } + SECTION("torus: per-face zero-energy condition") { primitive::TorusOptions torus_opts; @@ -591,6 +662,90 @@ TEST_CASE("compute_smooth_direction_field_on_facets", "[polyddg]") REQUIRE(cos4_diff > 0.85); } + SECTION("flat grid: alignment_lambda modulates alignment strength") + { + // Facet-path counterpart of the vertex-path test in the "compute_smooth_direction_field" + // TEST_CASE above: a flat grid has a zero-energy smooth mode (λ_1 ≈ 0), so λ_t = 0 is + // smoothness-dominated and ignores high-frequency guidance while a strongly negative + // λ_t recovers it. + constexpr Index N = 12; + SurfaceMesh grid; + for (Index i = 0; i < N; ++i) { + for (Index j = 0; j < N; ++j) { + grid.add_vertex({static_cast(i), static_cast(j), 0.0}); + } + } + for (Index i = 0; i + 1 < N; ++i) { + for (Index j = 0; j + 1 < N; ++j) { + const Index v00 = i * N + j, v10 = (i + 1) * N + j; + const Index v01 = i * N + (j + 1), v11 = (i + 1) * N + (j + 1); + grid.add_triangle(v00, v10, v11); + grid.add_triangle(v00, v11, v01); + } + } + polyddg::DifferentialOperators ops(grid); + + const Index nf = grid.get_num_facets(); + + // Dense, high-frequency guidance: alternate the target angle by ±22.5° (in each local + // frame) using (grid-cell parity) XOR (triangle type). This grid is split into two + // triangles per cell ("A" = (v00,v10,v11), even fid; "B" = (v00,v11,v01), odd fid), and + // each triangle's own local frame (facet_basis's first-edge convention) differs by a + // fixed 45° between the two types; the connection Laplacian's n=4 transport therefore + // rotates by a fixed 180° across every A/B edge, making "same target for A and B" (or a + // plain per-cell checkerboard) a zero-energy pattern the smoothest field reproduces for + // free. XOR-ing in the triangle type breaks that degeneracy so the pattern is genuinely + // high-frequency with respect to this connection. + const auto centroid_id = ops.get_centroid_attribute_id(); + const auto centroid_view = attribute_matrix_view(grid, centroid_id); + std::vector guide_angle(nf); + auto align_id = internal::find_or_create_attribute( + grid, + "@lambda_t_alignment_facets", + AttributeElement::Facet, + AttributeUsage::Vector, + 3, + internal::ResetToDefault::Yes); + auto align_data = attribute_matrix_ref(grid, align_id); + for (Index fid = 0; fid < nf; ++fid) { + const Index cell_i = static_cast(std::floor(centroid_view(fid, 0))); + const Index cell_j = static_cast(std::floor(centroid_view(fid, 1))); + const Index parity = (cell_i + cell_j + (fid % 2)) % 2; + const Scalar a = (parity == 0 ? 1.0 : -1.0) * (internal::pi / 8.0); + guide_angle[fid] = a; + align_data.row(fid) = + (ops.facet_basis(fid) * Eigen::Matrix(std::cos(a), std::sin(a))) + .transpose(); + } + + auto mean_align_error = [&](Scalar lambda_t, std::string_view out_name) { + polyddg::SmoothDirectionFieldOptions opts; + opts.nrosy = 4; + opts.output_element_type = AttributeElement::Facet; + opts.alignment_attribute = "@lambda_t_alignment_facets"; + opts.alignment_lambda = lambda_t; + opts.direction_field_attribute = out_name; + auto result = polyddg::compute_smooth_direction_field(grid, ops, opts); + auto data = attribute_matrix_view(grid, result); + Scalar err = 0; + for (Index fid = 0; fid < nf; ++fid) { + REQUIRE_THAT(data.row(fid).norm(), Catch::Matchers::WithinAbs(1.0, 1e-10)); + Eigen::Matrix out_2d = + ops.facet_basis(fid).transpose() * data.row(fid).transpose(); + Scalar out_angle = std::atan2(out_2d(1), out_2d(0)); + err += 1.0 - std::cos(4.0 * (out_angle - guide_angle[fid])); + } + return err / static_cast(nf); + }; + + // λ_t = 0 favors smoothness and cannot follow the alternating guidance; a strongly + // negative λ_t strengthens alignment and drives the mean error toward zero. + const Scalar err_default = mean_align_error(0.0, "@sdf_facets_lt_default"); + const Scalar err_strong = mean_align_error(-100.0, "@sdf_facets_lt_strong"); + REQUIRE(err_default > 0.5); + REQUIRE(err_strong < 0.1); + } + SECTION("convenience overload (no ops argument)") { SurfaceMesh mesh; diff --git a/modules/python/lagrange/scripts/meshstat.py b/modules/python/lagrange/scripts/meshstat.py index ce339c81..28541867 100644 --- a/modules/python/lagrange/scripts/meshstat.py +++ b/modules/python/lagrange/scripts/meshstat.py @@ -11,8 +11,6 @@ # """Print basic information about a mesh file.""" -from __future__ import annotations - import argparse import json import logging diff --git a/modules/scene/include/lagrange/scene/scene_convert.h b/modules/scene/include/lagrange/scene/scene_convert.h index 85934f7f..8a318820 100644 --- a/modules/scene/include/lagrange/scene/scene_convert.h +++ b/modules/scene/include/lagrange/scene/scene_convert.h @@ -74,7 +74,10 @@ SurfaceMesh scene_to_mesh( /// @tparam Scalar Input scene scalar type. /// @tparam Index Input scene index type. /// -/// @return List of meshes with transforms applied. +/// @return List of meshes with transforms applied. Each mesh carries a `material_id` facet +/// attribute (see @ref AttributeName::material_id) whose values are global indices into +/// `Scene::materials`. Facets with no associated material are marked with +/// `invalid()`. /// template std::vector> scene_to_meshes( @@ -125,15 +128,30 @@ Scene simple_scene_to_scene(const SimpleScene& sim template struct MeshesAndMaterialsResult { - /// List of meshes with transforms applied. + /// List of meshes with transforms applied. Each mesh carries a `material_id` facet attribute + /// (see @ref AttributeName::material_id) whose values are global indices into the scene's + /// material list (i.e. the same indexing as `Scene::materials` and `material_ids` below). + /// Facets with no associated material are marked with `invalid()`. std::vector> meshes; - /// List of material IDs for each mesh. + /// List of material IDs used by each mesh (indices into `Scene::materials`). This is the set of + /// materials referenced by the corresponding mesh instance; use the per-facet `material_id` + /// attribute on each mesh to recover the exact facet-to-material assignment. std::vector> material_ids; }; /// -/// Converts a scene into a list of meshes with all the transforms applied and a list of material IDs. +/// Converts a scene into a list of meshes with all the transforms applied and a list of material +/// IDs. +/// +/// In addition to the per-mesh `material_ids` list, every output mesh is given a `material_id` +/// facet attribute holding the global scene material index used by each facet, so that the +/// facet-to-material assignment is preserved even when a single mesh uses multiple materials. +/// +/// **Convention for existing `material_id` attributes:** If an input mesh already carries a +/// `material_id` facet attribute, its values are interpreted as *instance-local* indices into the +/// corresponding `SceneMeshInstance::materials` list. Because instances are flattened and meshes +/// are duplicated, this function remaps them to global `Scene::materials` indices on output. /// /// @param[in] scene Scene to convert. /// @param[in] transform_options Options to use when applying mesh transformations. diff --git a/modules/scene/python/src/bind_scene.h b/modules/scene/python/src/bind_scene.h index 5d325110..7bc0f35a 100644 --- a/modules/scene/python/src/bind_scene.h +++ b/modules/scene/python/src/bind_scene.h @@ -783,6 +783,10 @@ void bind_scene(nb::module_& m) "reorient"_a = TransformOptions{}.reorient, R"(Converts a scene into a list of meshes with all the transforms applied and a list of material IDs. +Each output mesh also carries a ``material_id`` facet attribute holding the global scene material +index used by each facet, so the facet-to-material assignment is preserved even when a single mesh +uses multiple materials. Facets with no associated material are marked with ``invalid_index``. + :param scene: Scene to convert. :param normalize_normals: If enabled, normals are normalized after transformation. :param normalize_tangents_bitangents: If enabled, tangents and bitangents are normalized after transformation. diff --git a/modules/scene/src/scene_convert.cpp b/modules/scene/src/scene_convert.cpp index 4ac5d998..dc0ad1c2 100644 --- a/modules/scene/src/scene_convert.cpp +++ b/modules/scene/src/scene_convert.cpp @@ -11,12 +11,17 @@ */ #include +#include +#include #include +#include +#include #include #include #include #include #include +#include namespace lagrange::scene { @@ -47,6 +52,7 @@ MeshesAndMaterialsResult scene_to_meshes_and_materials( const TransformOptions& transform_options) { MeshesAndMaterialsResult ret; + ElementId num_materials = scene.materials.size(); for (ElementId node_id = 0; node_id < scene.nodes.size(); ++node_id) { const auto& node = scene.nodes[node_id]; @@ -57,12 +63,100 @@ MeshesAndMaterialsResult scene_to_meshes_and_materials( utils::compute_global_node_transform(scene, node_id).template cast(); for (const SceneMeshInstance& mesh_instance : node.meshes) { const auto mesh_id = mesh_instance.mesh; - ret.meshes.emplace_back( - transformed_mesh( - scene.meshes.at(mesh_id), - world_from_mesh, - transform_options)); - ret.material_ids.push_back(mesh_instance.materials); + auto mesh = transformed_mesh( + scene.meshes.at(mesh_id), + world_from_mesh, + transform_options); + + // Write a self-contained, global material_id facet attribute on the output mesh so + // that each facet directly identifies the scene material it uses. Values index into + // scene.materials (and the parallel material_ids[i] list); facets with no material are + // marked with invalid(). + const auto& materials = mesh_instance.materials; + if (mesh.has_attribute(AttributeName::material_id)) { + // Per the SceneMeshInstance convention, an existing material_id facet attribute + // holds instance-local indices into mesh_instance.materials. Remap those to global + // scene material indices. + const auto& attr_base = mesh.get_attribute_base(AttributeName::material_id); + la_runtime_assert( + attr_base.get_element_type() == AttributeElement::Facet && + attr_base.get_num_channels() == 1, + "existing material_id attribute must be a 1-channel Facet attribute"); + if (!mesh.template is_attribute_type(AttributeName::material_id)) { + cast_attribute_in_place(mesh, AttributeName::material_id); + } + auto values = + mesh.template ref_attribute(AttributeName::material_id).ref_all(); + size_t num_invalid_slots = 0; + size_t num_invalid_materials = 0; + for (auto& value : values) { + if (value < materials.size()) { + const auto global_id = materials[value]; + if (global_id < num_materials) { + value = static_cast(global_id); + } else { + value = invalid(); + ++num_invalid_materials; + } + } else { + value = invalid(); + ++num_invalid_slots; + } + } + if (num_invalid_slots || num_invalid_materials) { + logger().warn( + "Scene node {} mesh {} has a material_id facet attribute with {} invalid " + "slots and {} invalid materials. Invalid slots/materials will be marked as " + "invalid().", + node_id, + mesh_id, + num_invalid_slots, + num_invalid_materials); + } + } else { + // No per-facet attribute: assign all facets to the first listed material (or + // invalid if none). When the instance lists multiple materials, picking the first + // is the best we can do without facet-level data. + const Index global_id = + materials.empty() ? invalid() : static_cast(materials.front()); + const Index value = global_id != invalid() && global_id < num_materials + ? static_cast(global_id) + : invalid(); + std::vector values(mesh.get_num_facets(), value); + mesh.template create_attribute( + AttributeName::material_id, + AttributeElement::Facet, + AttributeUsage::Scalar, + 1, + values); + if (materials.size() > 1) { + logger().warn( + "Scene node {} mesh {} has no material_id facet attribute, but the " + "instance lists {} materials. All facets will be assigned to the first " + "material ({}).", + node_id, + mesh_id, + materials.size(), + global_id); + } + if (global_id == invalid()) { + logger().warn( + "Scene node {} mesh {} has no material. All facets will be marked as " + "invalid().", + node_id, + mesh_id); + } else if (global_id >= num_materials) { + logger().warn( + "Scene node {} mesh {} references material {} which is out of range. All " + "facets will be marked as invalid().", + node_id, + mesh_id, + global_id); + } + } + + ret.meshes.emplace_back(std::move(mesh)); + ret.material_ids.push_back(materials); } } diff --git a/modules/scene/tests/test_scene.cpp b/modules/scene/tests/test_scene.cpp index 101c247e..0ff2a82e 100644 --- a/modules/scene/tests/test_scene.cpp +++ b/modules/scene/tests/test_scene.cpp @@ -11,10 +11,12 @@ */ #include +#include #include #include #include #include +#include #include #include @@ -119,6 +121,95 @@ TEST_CASE("Scene: convert", "[scene]") REQUIRE(facet_view(mesh) == facet_view(mesh2)); } +TEST_CASE("Scene: scene_to_meshes_and_materials material_id", "[scene]") +{ + using namespace lagrange; + using Scalar = double; + using Index = uint32_t; + + auto make_mesh = []() { + SurfaceMesh mesh; + mesh.add_vertices(4); + vertex_ref(mesh).setRandom(); + mesh.add_triangle(0, 1, 2); + mesh.add_triangle(1, 2, 3); + return mesh; + }; + + SECTION("single material per instance") + { + scene::Scene scene; + scene::ElementId mesh0 = scene.add(make_mesh()); + scene::ElementId mesh1 = scene.add(make_mesh()); + scene::ElementId mat0 = scene.add(scene::MaterialExperimental{}); + scene::ElementId mat1 = scene.add(scene::MaterialExperimental{}); + + scene::Node node; + node.meshes.push_back({mesh0, {mat1}}); + node.meshes.push_back({mesh1, {mat0}}); + scene.root_nodes.push_back(scene.add(std::move(node))); + + auto result = scene::scene_to_meshes_and_materials(scene); + REQUIRE(result.meshes.size() == 2); + REQUIRE(result.material_ids.size() == 2); + + // Every facet of the first output mesh uses the instance's single (global) material. + REQUIRE(result.meshes[0].has_attribute(AttributeName::material_id)); + auto ids0 = attribute_vector_view(result.meshes[0], AttributeName::material_id); + REQUIRE((ids0.array() == static_cast(mat1)).all()); + + auto ids1 = attribute_vector_view(result.meshes[1], AttributeName::material_id); + REQUIRE((ids1.array() == static_cast(mat0)).all()); + } + + SECTION("no material") + { + scene::Scene scene; + scene::ElementId mesh0 = scene.add(make_mesh()); + + scene::Node node; + node.meshes.push_back({mesh0, {}}); + scene.root_nodes.push_back(scene.add(std::move(node))); + + auto result = scene::scene_to_meshes_and_materials(scene); + REQUIRE(result.meshes.size() == 1); + + auto ids = attribute_vector_view(result.meshes[0], AttributeName::material_id); + REQUIRE((ids.array() == invalid()).all()); + } + + SECTION("multiple materials remapped from local indices") + { + // A mesh carrying a per-facet material_id attribute holds instance-local indices, which + // must be remapped to global scene material indices on output. + auto mesh = make_mesh(); + std::vector local_ids = {1, 0}; // facet 0 -> materials[1], facet 1 -> materials[0] + mesh.template create_attribute( + AttributeName::material_id, + AttributeElement::Facet, + AttributeUsage::Scalar, + 1, + local_ids); + + scene::Scene scene; + scene::ElementId mesh0 = scene.add(std::move(mesh)); + scene::ElementId mat_a = scene.add(scene::MaterialExperimental{}); + scene.add(scene::MaterialExperimental{}); // unreferenced material + scene::ElementId mat_b = scene.add(scene::MaterialExperimental{}); + + scene::Node node; + node.meshes.push_back({mesh0, {mat_a, mat_b}}); // local 0 -> mat_a, local 1 -> mat_b + scene.root_nodes.push_back(scene.add(std::move(node))); + + auto result = scene::scene_to_meshes_and_materials(scene); + REQUIRE(result.meshes.size() == 1); + + auto ids = attribute_vector_view(result.meshes[0], AttributeName::material_id); + REQUIRE(ids(0) == static_cast(mat_b)); // local 1 + REQUIRE(ids(1) == static_cast(mat_a)); // local 0 + } +} + TEST_CASE("Scene: scene_to_simple_scene empty", "[scene]") { using Scalar = double; diff --git a/modules/testing/include/lagrange/testing/common.h b/modules/testing/include/lagrange/testing/common.h index 238ab632..04095d53 100644 --- a/modules/testing/include/lagrange/testing/common.h +++ b/modules/testing/include/lagrange/testing/common.h @@ -180,5 +180,14 @@ SurfaceMesh load_surface_mesh(const fs::path& relative_path) /// LA_TESTING_API void setup_mkl_reproducibility(); +/// +/// Disable interactive Windows error dialogs (CRT assert/error report windows, the abort() message +/// box, and Windows Error Reporting popups) and route their output to stderr instead, so that a +/// failed assert terminates the process immediately rather than blocking a headless CI agent until +/// the pipeline timeout. This function has no effect on other platforms, and is a no-op when a +/// debugger is attached (to preserve interactive debugging behavior). +/// +LA_TESTING_API void disable_windows_error_dialogs(); + } // namespace testing } // namespace lagrange diff --git a/modules/testing/main/main.cpp b/modules/testing/main/main.cpp index 3f24b8d9..aab25344 100644 --- a/modules/testing/main/main.cpp +++ b/modules/testing/main/main.cpp @@ -10,6 +10,7 @@ * governing permissions and limitations under the License. */ #include +#include #include #include @@ -74,6 +75,8 @@ static_assert(false, "Emscripten must be compiled with pthreads support"); int main(int argc, char* argv[]) { + lagrange::testing::disable_windows_error_dialogs(); + Catch::Session session; int log_level = spdlog::level::warn; bool fpe_flag = false; diff --git a/modules/testing/src/common.cpp b/modules/testing/src/common.cpp index 7c2db2cf..6987811f 100644 --- a/modules/testing/src/common.cpp +++ b/modules/testing/src/common.cpp @@ -25,9 +25,40 @@ #include #endif +#ifdef _WIN32 + #ifndef WIN32_LEAN_AND_MEAN + #define WIN32_LEAN_AND_MEAN + #endif + #ifndef NOMINMAX + #define NOMINMAX + #endif + #include + + #include + #include + #include +#endif + namespace lagrange { namespace testing { +#if defined(_WIN32) && defined(_DEBUG) +namespace { + +// Terminate on Debug CRT assert/error reports instead of continuing with "Ignore" semantics. +int crt_report_hook(int report_type, char* message, int* /*return_value*/) +{ + if (report_type == _CRT_ERROR || report_type == _CRT_ASSERT) { + fputs(message, stderr); + fflush(stderr); + abort(); + } + return FALSE; +} + +} // namespace +#endif + fs::path get_data_dir() { #ifdef TEST_DATA_DIR @@ -153,5 +184,36 @@ void setup_mkl_reproducibility() #endif } +void disable_windows_error_dialogs() +{ +#ifdef _WIN32 + // Keep interactive dialogs and debug breaks when running under a debugger. + if (IsDebuggerPresent()) { + return; + } + + // Route assert() failure messages to stderr instead of a message box. + _set_error_mode(_OUT_TO_STDERR); + + // Write the abort() message to stderr, but don't invoke Windows Error Reporting. + _set_abort_behavior(_WRITE_ABORT_MSG, _WRITE_ABORT_MSG | _CALL_REPORTFAULT); + + // Suppress OS-level error popups for hard crashes (e.g. access violations). + SetErrorMode( + GetErrorMode() | SEM_FAILCRITICALERRORS | SEM_NOGPFAULTERRORBOX | SEM_NOOPENFILEERRORBOX); + + #ifdef _DEBUG + // Route Debug CRT reports to stderr instead of the "Debug Assertion Failed" dialog. + for (int report_type : {_CRT_WARN, _CRT_ERROR, _CRT_ASSERT}) { + _CrtSetReportMode(report_type, _CRTDBG_MODE_FILE); + _CrtSetReportFile(report_type, _CRTDBG_FILE_STDERR); + } + + // Ensure assert/error reports terminate the process instead of merely logging. + _CrtSetReportHook(crt_report_hook); + #endif +#endif +} + } // namespace testing } // namespace lagrange diff --git a/modules/ui/src/types/Camera_xcode264_workaround.cpp b/modules/ui/src/types/Camera_xcode264_workaround.cpp index 76864022..534f68ed 100644 --- a/modules/ui/src/types/Camera_xcode264_workaround.cpp +++ b/modules/ui/src/types/Camera_xcode264_workaround.cpp @@ -16,7 +16,7 @@ // prevents the crash. This works across all build systems (CMake, Make, etc.) // // Bug: VectorCombine::foldSelectShuffle() segfaults on (r_0 * r_arc).inverse() -// Affects: Apple Clang (Xcode 26.4, __apple_build_version__ 21000000-21999999), x86_64 only +// Affects: Apple Clang (Xcode 26.4, __apple_build_version__ 21000099), x86_64 only // // Why this function is in a separate file with file-level #pragma: // - Localized #pragma around the function (function-level) does NOT work @@ -27,15 +27,25 @@ // - Therefore, the #pragma must be at FILE level to disable the pass for the // entire compilation unit containing the problematic template instantiations // -// TODO: Remove this file when Apple fixes the bug in future Xcode release +// TODO: Remove this file once we stop supporting Xcode 26... #include // Only disable optimizations for the specific problematic configuration -// Note: __apple_build_version__ is only defined by Apple Clang, not Homebrew LLVM +// +// Note: +// - __apple_build_version__ is only defined by Apple Clang, not Homebrew LLVM +// - Bug appears in Xcode 26.4 (clang-2100.0.123.102) +// - Got fixed in Xcode-27.0.0-Beta.3 (clang-2100.3.25.1) +// +// Version numbers: +// - 21000099 = Xcode 26.4.0 +// - 21000323 = Xcode 27.0.0 Beta 2 +// - 21000325 = Xcode 27.0.0 Beta 3 + #if LAGRANGE_TARGET_OS(APPLE) && LAGRANGE_TARGET_PLATFORM(x86_64) && \ defined(__apple_build_version__) && __apple_build_version__ >= 21000000 && \ - __apple_build_version__ < 22000000 + __apple_build_version__ < 21000325 #pragma clang optimize off #endif @@ -119,6 +129,6 @@ void Camera::rotate_arcball( // Re-enable optimizations if they were disabled #if LAGRANGE_TARGET_OS(APPLE) && LAGRANGE_TARGET_PLATFORM(x86_64) && \ defined(__apple_build_version__) && __apple_build_version__ >= 21000000 && \ - __apple_build_version__ < 22000000 + __apple_build_version__ < 21000325 #pragma clang optimize on #endif diff --git a/modules/ui/src/utils/math_xcode264_workaround.cpp b/modules/ui/src/utils/math_xcode264_workaround.cpp index 6f4557db..e9538fe1 100644 --- a/modules/ui/src/utils/math_xcode264_workaround.cpp +++ b/modules/ui/src/utils/math_xcode264_workaround.cpp @@ -16,7 +16,7 @@ // prevents the crash. This works across all build systems (CMake, Make, etc.) // // Bug: VectorCombine::foldSelectShuffle() segfaults on (perspective * view).inverse() -// Affects: Apple Clang (Xcode 26.4, __apple_build_version__ 21000000-21999999), x86_64 only +// Affects: Apple Clang (Xcode 26.4, __apple_build_version__ 21000099), x86_64 only // // Why this function is in a separate file with file-level #pragma: // - Localized #pragma around the function (function-level) does NOT work @@ -27,15 +27,25 @@ // - Therefore, the #pragma must be at FILE level to disable the pass for the // entire compilation unit containing the problematic template instantiations // -// TODO: Remove this file when Apple fixes the bug in future Xcode release +// TODO: Remove this file once we stop supporting Xcode 26... #include // Only disable optimizations for the specific problematic configuration -// Note: __apple_build_version__ is only defined by Apple Clang, not Homebrew LLVM +// +// Note: +// - __apple_build_version__ is only defined by Apple Clang, not Homebrew LLVM +// - Bug appears in Xcode 26.4 (clang-2100.0.123.102) +// - Got fixed in Xcode-27.0.0-Beta.3 (clang-2100.3.25.1) +// +// Version numbers: +// - 21000099 = Xcode 26.4.0 +// - 21000323 = Xcode 27.0.0 Beta 2 +// - 21000325 = Xcode 27.0.0 Beta 3 + #if LAGRANGE_TARGET_OS(APPLE) && LAGRANGE_TARGET_PLATFORM(x86_64) && \ defined(__apple_build_version__) && __apple_build_version__ >= 21000000 && \ - __apple_build_version__ < 22000000 + __apple_build_version__ < 21000325 #pragma clang optimize off #endif @@ -67,6 +77,6 @@ Eigen::Vector3f unproject_point( // Re-enable optimizations if they were disabled #if LAGRANGE_TARGET_OS(APPLE) && LAGRANGE_TARGET_PLATFORM(x86_64) && \ defined(__apple_build_version__) && __apple_build_version__ >= 21000000 && \ - __apple_build_version__ < 22000000 + __apple_build_version__ < 21000325 #pragma clang optimize on #endif diff --git a/modules/ui/tests/main.cpp b/modules/ui/tests/main.cpp index 6c490475..8e23c24d 100644 --- a/modules/ui/tests/main.cpp +++ b/modules/ui/tests/main.cpp @@ -30,6 +30,8 @@ struct MiniGLContext int main(int argc, char* argv[]) { + lagrange::testing::disable_windows_error_dialogs(); + #ifdef LAGRANGE_UI_OPENGL_TESTS MiniGLContext opengl_context; #endif diff --git a/modules/xatlas/examples/unwrap_xatlas.cpp b/modules/xatlas/examples/unwrap_xatlas.cpp index 4639ee99..1978fedd 100644 --- a/modules/xatlas/examples/unwrap_xatlas.cpp +++ b/modules/xatlas/examples/unwrap_xatlas.cpp @@ -35,7 +35,7 @@ int main(int argc, char** argv) struct { fs::path input; - fs::path output = "output.ply"; + fs::path output = "output.obj"; fs::path uv_mesh_output; std::string uv_attribute_name = "texcoord"; std::string atlas_attribute_name;