Sync with f7df06a (#91)

Co-authored-by: Qingnan Zhou <qzhou@Qingnans-MacBook-Pro.local>

Author: qnzhou

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()

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
 )
 

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 <lagrange/SurfaceMesh.h>
+#include <lagrange/bvh/api.h>
+#include <lagrange/bvh/resolve_tjunctions.h>
+
+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 <typename Scalar, typename Index>
+void resolve_tjunctions(SurfaceMesh<Scalar, Index>& mesh, ResolveTJunctionsOptions options = {});
+
+} // namespace lagrange::bvh::internal

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 <lagrange/SurfaceMesh.h>
+#include <lagrange/bvh/api.h>
+
+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 <typename Scalar, typename Index>
+void resolve_tjunctions(SurfaceMesh<Scalar, Index>& mesh, ResolveTJunctionsOptions options = {});
+
+/// @}
+
+} // namespace lagrange::bvh

modules/bvh/python/src/bvh.cpp

@@ -16,6 +16,7 @@
 #include <lagrange/bvh/compute_mesh_distances.h>
 #include <lagrange/bvh/compute_uv_overlap.h>
 #include <lagrange/bvh/remove_interior_shells.h>
+#include <lagrange/bvh/resolve_tjunctions.h>
 #include <lagrange/bvh/weld_vertices.h>
 #include <lagrange/python/binding.h>
 #include <lagrange/python/bvh.h>
@@ -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

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 <lagrange/bvh/internal/resolve_tjunctions.h>
+
+#include <lagrange/SurfaceMeshTypes.h>
+#include <lagrange/bvh/AABB.h>
+#include <lagrange/bvh/api.h>
+#include <lagrange/internal/split_edges.h>
+#include <lagrange/triangulate_polygonal_facets.h>
+#include <lagrange/utils/assert.h>
+#include <lagrange/utils/function_ref.h>
+#include <lagrange/utils/span.h>
+#include <lagrange/views.h>
+
+// clang-format off
+#include <lagrange/utils/warnoff.h>
+#include <tbb/parallel_for.h>
+#include <lagrange/utils/warnon.h>
+// clang-format on
+
+#include <Eigen/Geometry>
+
+#include <algorithm>
+#include <utility>
+#include <vector>
+
+namespace lagrange::bvh::internal {
+
+namespace {
+
+template <typename Scalar, typename Index, int Dim>
+void resolve_tjunctions_impl(
+    SurfaceMesh<Scalar, Index>& mesh,
+    const ResolveTJunctionsOptions& options)
+{
+    using Tree = AABB<Scalar, Dim>;
+    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<Scalar>(options.tolerance);
+    if (options.tolerance < 0) {
+        Scalar diag = (vertices.colwise().maxCoeff() - vertices.colwise().minCoeff()).norm();
+        tol = static_cast<Scalar>(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<Box> boxes(static_cast<size_t>(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<std::vector<std::pair<Scalar, Index>>> 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<bool(typename Tree::Index)>([&](typename Tree::Index candidate) {
+                const Index v = static_cast<Index>(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<Index> edge_split_offsets(num_edges + 1, 0);
+    std::vector<Index> 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<Index>(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<span<Index>(Index)>([&](Index e) -> span<Index> {
+            const Index n = edge_split_offsets[e + 1] - edge_split_offsets[e];
+            return span<Index>(split_pts.data() + edge_split_offsets[e], n);
+        }),
+        function_ref<bool(Index)>([](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<bool(Index)>(is_new_facet));
+    }
+    mesh.remove_facets(facets_to_remove);
+}
+
+} // namespace
+
+template <typename Scalar, typename Index>
+void resolve_tjunctions(SurfaceMesh<Scalar, Index>& mesh, ResolveTJunctionsOptions options)
+{
+    const Index dim = mesh.get_dimension();
+    if (dim == 2) {
+        resolve_tjunctions_impl<Scalar, Index, 2>(mesh, options);
+    } else if (dim == 3) {
+        resolve_tjunctions_impl<Scalar, Index, 3>(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<Scalar, Index>( \
+        SurfaceMesh<Scalar, Index>&,                            \
+        ResolveTJunctionsOptions);
+LA_SURFACE_MESH_X(resolve_tjunctions, 0)
+
+} // namespace lagrange::bvh::internal