WIP simplify boundaries

petrasvestartas · petrasvestartas · commit 8b4ee00c4084 · 2025-06-12T16:58:54.000+02:00
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -20,6 +20,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Changed
 
+- Fixed point welding bug in map_mesh function.
+
 ### Removed
 
 
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -4,7 +4,7 @@ project(compas_libigl LANGUAGES CXX)
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_EXTENSIONS OFF)
 
-option(ENABLE_PRECOMPILED_HEADERS "Enable precompiled headers" OFF)
+option(ENABLE_PRECOMPILED_HEADERS "Enable precompiled headers" ON)
 option(MULTITHREADED_COMPILATION "Enable multi-threaded compilation (Ninja only)" ON)
 
 if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
@@ -152,13 +152,13 @@ endfunction()
 # === Add your modules ===
 add_nanobind_module(_nanobind src/nanobind.cpp)
 add_nanobind_module(_types_std src/types_std.cpp)
-add_nanobind_module(_boundaries src/boundaries.cpp)
-add_nanobind_module(_curvature src/curvature.cpp)
-add_nanobind_module(_geodistance src/geodistance.cpp)
-add_nanobind_module(_intersections src/intersections.cpp)
-add_nanobind_module(_isolines src/isolines.cpp)
-add_nanobind_module(_massmatrix src/massmatrix.cpp)
-add_nanobind_module(_meshing src/meshing.cpp)
-add_nanobind_module(_parametrisation src/parametrisation.cpp)
-add_nanobind_module(_planarize src/planarize.cpp)
+# add_nanobind_module(_boundaries src/boundaries.cpp)
+# add_nanobind_module(_curvature src/curvature.cpp)
+# add_nanobind_module(_geodistance src/geodistance.cpp)
+# add_nanobind_module(_intersections src/intersections.cpp)
+# add_nanobind_module(_isolines src/isolines.cpp)
+# add_nanobind_module(_massmatrix src/massmatrix.cpp)
+# add_nanobind_module(_meshing src/meshing.cpp)
+# add_nanobind_module(_parametrisation src/parametrisation.cpp)
+# add_nanobind_module(_planarize src/planarize.cpp)
 add_nanobind_module(_mapping src/mapping.cpp)
diff --git a/docs/examples/example_mapping_patterns.py b/docs/examples/example_mapping_patterns.py
@@ -23,8 +23,8 @@
 # Mapping: 3D Mesh, 2D Pattern, UV
 # ==============================================================================
 
-mesh_mapped = map_pattern_to_mesh("ZigZag", mesh, clip_boundaries=True, tolerance=1e-6, pattern_u=16, pattern_v=16)
-
+mesh_mapped0 = map_pattern_to_mesh("ZigZag", mesh, clip_boundaries=True, tolerance=1e-6, pattern_u=16, pattern_v=16, simplify_borders=True)
+mesh_mapped1 = map_pattern_to_mesh("ZigZag", mesh, clip_boundaries=True, tolerance=1e-6, pattern_u=16, pattern_v=16, simplify_borders=False)
 # ==============================================================================
 # Viewer
 # ==============================================================================
@@ -36,6 +36,8 @@
 viewer = Viewer(config=config)
 
 # viewer.scene.add(mesh, name="mesh", show_faces=False, linecolor=Color.grey(), opacity=0.2)
-viewer.scene.add(mesh_mapped, name="mesh_mapped", facecolor=Color.red())
+viewer.scene.add(mesh_mapped0, name="mesh_mapped0", facecolor=Color.red(), show_points=True)
+viewer.scene.add(mesh_mapped1, name="mesh_mapped1", facecolor=Color.blue(), show_points=True, show_faces=False)
+
 
 viewer.show()
diff --git a/src/compas_libigl/mapping.py b/src/compas_libigl/mapping.py
@@ -26,7 +26,7 @@
 from compas_libigl._types_std import VectorVectorInt  # noqa: F401
 
 
-def map_mesh(target_mesh, pattern_mesh, clip_boundaries=True, tolerance=1e-6):
+def map_mesh(target_mesh, pattern_mesh, clip_boundaries=True, simplify_borders=True, tolerance=1e-6):
     """
     Map a 2D pattern mesh onto a 3D target.
 
@@ -38,6 +38,8 @@ def map_mesh(target_mesh, pattern_mesh, clip_boundaries=True, tolerance=1e-6):
         A tuple of (vertices, faces) representing the pattern mesh.
     clip_boundaries : bool
         Whether to clip the pattern mesh to the boundaries of the target mesh.
+    simplify_borders : bool
+        Whether to simplify the border of the pattern mesh.
     tolerance : float
         The tolerance for point comparison, to remove duplicates.
 
@@ -63,14 +65,14 @@ def map_mesh(target_mesh, pattern_mesh, clip_boundaries=True, tolerance=1e-6):
 
     # Perform the mapping
     pv_numpy_copy, pf_numpy_cleaned, p_normals, pattern_is_boundary, pattern_groups = _mapping.map_mesh_with_automatic_parameterization(
-        v_numpy, f_numpy, pattern_v_numpy, pf, clip_boundaries, tolerance
+        v_numpy, f_numpy, pattern_v_numpy, pf, clip_boundaries, simplify_borders, tolerance
     )
 
     # Return the result as a tuple
     return pv_numpy_copy, pf_numpy_cleaned, p_normals, pattern_is_boundary, pattern_groups
 
 
-def map_pattern_to_mesh(name, mesh, clip_boundaries=True, tolerance=1e-6, pattern_u=16, pattern_v=16):
+def map_pattern_to_mesh(name, mesh, clip_boundaries=True, tolerance=1e-6, pattern_u=16, pattern_v=16, simplify_borders=True):
     """
     Map a 2D pattern mesh onto a 3D target.
 
@@ -108,6 +110,8 @@ def map_pattern_to_mesh(name, mesh, clip_boundaries=True, tolerance=1e-6, patter
         The number of pattern vertices in the u direction.
     pattern_v : int
         The number of pattern vertices in the v direction.
+    simplify_borders : bool
+        Whether to simplify the border of the pattern mesh.
 
     Returns
     -------
@@ -165,6 +169,6 @@ def map_pattern_to_mesh(name, mesh, clip_boundaries=True, tolerance=1e-6, patter
     pf = tessagon_mesh["face_list"]
 
     v, f = mesh.to_vertices_and_faces()
-    mapped_vertices, mapped_faces, mapped_normals, mapped_is_boundary, mapped_groups = map_mesh((v, f), (pv, pf), clip_boundaries=clip_boundaries, tolerance=tolerance)
+    mapped_vertices, mapped_faces, mapped_normals, mapped_is_boundary, mapped_groups = map_mesh((v, f), (pv, pf), clip_boundaries=clip_boundaries, simplify_borders=simplify_borders, tolerance=tolerance)
 
     return Mesh.from_vertices_and_faces(mapped_vertices, mapped_faces)
diff --git a/src/mapping.cpp b/src/mapping.cpp
@@ -1,5 +1,9 @@
 #include "mapping.hpp"
-
+#include <algorithm>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <unordered_set>
 
 // Custom hash function for tuple
 struct TupleHash {
@@ -18,13 +22,10 @@ std::tuple<int64_t, int64_t> grid_key(double x, double y, double tolerance) {
     };
 }
 
-bool is_same_point(double x1, double y1, double z1,
-    double x2, double y2, double z2,
-    double tol) {
-        double dx = x1 - x2;
-        double dy = y1 - y2;
-        double dz = z1 - z2;
-        return (dx*dx + dy*dy + dz*dz) < (tol*tol);
+bool is_same_point(double x1, double y1, double x2, double y2, double tol) {
+    double dx = x1 - x2;
+    double dy = y1 - y2;
+    return (dx*dx + dy*dy) < (tol*tol);
 }
 
 std::vector<std::vector<int>> map_mesh_cropped(
@@ -98,7 +99,7 @@ std::vector<std::vector<int>> map_mesh_cropped(
         }
         
         // Update pattern vertex position through barycentric interpolation, comment this out if you keep 2D pattern
-        pattern_v.row(id) = A * L(0, 0) + B * L(0, 1) + C * L(0, 2);
+        // pattern_v.row(id) = A * L(0, 0) + B * L(0, 1) + C * L(0, 2);
         
         // If normals are requested, interpolate them using the same barycentric coordinates
         if (compute_normals) {
@@ -166,13 +167,16 @@ bool paths_intersect(const Clipper2Lib::PathsD& paths1, const Clipper2Lib::Paths
 }
 
 
+
+
 std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, std::vector<bool>, std::vector<int>> eigen_to_clipper (
     Eigen::Ref<const compas::RowMatrixXd> flattned_target_uv,
     Eigen::Ref<const compas::RowMatrixXi> target_f, 
     
     Eigen::Ref<const compas::RowMatrixXd> pattern_v, 
     const std::vector<std::vector<int>>& pattern_f,
     bool clip_boundaries,
+    bool simplify_borders,
     double tolerance
 )   
 {
@@ -217,21 +221,18 @@ std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, std::vector<bool>
             path.emplace_back(pattern_v(point_id, 0), pattern_v(point_id, 1));
         paths.emplace_back(path);
 
+        if (!clip_boundaries){
+            patterns_to_keep.push_back(paths);
+            continue;
+        }
+
         // Get bounds of current pattern polygon
         Clipper2Lib::RectD pattern_bounds = Clipper2Lib::GetBounds(paths);
         
         // Check if the pattern bounds intersect with boundary bounds
         if (!pattern_bounds.Intersects(boundary_bounds))
             continue;
-
-        // Check if Elements are inside the boundary
-        // Clipper2Lib::PointD corners[] = {
-        //     {pattern_bounds.left, pattern_bounds.top},      // Top-left
-        //     {pattern_bounds.right, pattern_bounds.top},     // Top-right
-        //     {pattern_bounds.right, pattern_bounds.bottom},  // Bottom-right
-        //     {pattern_bounds.left, pattern_bounds.bottom}    // Bottom-left
-        // };
-
+    
         bool is_rect_inside_polygon = false;
         size_t num_corners_in_polygon = 0;
         for (const auto& corner : paths[0]) {
@@ -241,8 +242,6 @@ std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, std::vector<bool>
                 is_rect_inside_polygon = true;
                 bool has_collision = paths_intersect(paths, boundary);
 
-
-                
                 // check if the point is not in a hole
                 bool is_rect_outside_polygon_hole = true;
                 for (size_t i = 1; i < boundary.size(); i++) {                    
@@ -274,12 +273,10 @@ std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, std::vector<bool>
         // Fully enclosed polygons are added to the keep list
         // Other polygons edges are intersected with the boundary
 
-        if (num_corners_in_polygon == paths[0].size() && !is_hole_in_polygon){
+        if (num_corners_in_polygon == paths[0].size() && !is_hole_in_polygon){ //  
             patterns_to_keep.push_back(paths);
-        }else if (clip_boundaries){
-            if (paths_intersect(paths, boundary) || is_hole_in_polygon){
-                patterns_to_cut.push_back(paths);
-            }
+        }else if(num_corners_in_polygon > 0 || paths_intersect(paths, boundary) || is_hole_in_polygon){
+            patterns_to_cut.push_back(paths);
         }
 
 
@@ -290,10 +287,56 @@ std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, std::vector<bool>
     std::vector<Clipper2Lib::PathsD> solutions;
     for (const auto &subject : patterns_to_cut){
         Clipper2Lib::PathsD solution = Clipper2Lib::Intersect(subject, boundary, Clipper2Lib::FillRule::NonZero, 8);
-        if (solution.size() == 0)
+
+        if (solution.size() == 0){
+            // If no intersection is found, skip this polygon
             continue;
+        }else if (solution.size() > 1 || !simplify_borders){
+            // Highly likely we do not consider polygons with holes, we do not simplify them.
+            solutions.push_back(solution);
+        }else{
+            // After boolean intersection the pattern is merged with boundary polygons. This is often not wanted.
+            // Simplification is made by comparing polygon after boolean intersection and before. If the points lie on the initial polygon edges we keep them.
+            // Additionally we add an option to keep arbitrary points from a user given points or boundary valence points.
+            Clipper2Lib::PathsD simplified_paths{Clipper2Lib::PathD()};
+            simplified_paths[0].reserve(subject[0].size());
+
+            for (const auto &p : solution[0]){ // iterate the points of the intersection polygon
+                // First check if point is close to any vertex (corner)
+                bool point_added = false;
+                for(size_t i = 0; i < subject[0].size(); i++){
+                    double dx = p.x - subject[0][i].x;
+                    double dy = p.y - subject[0][i].y;
+                    if (dx*dx + dy*dy < tolerance*tolerance){
+                        simplified_paths[0].push_back(p);
+                        point_added = true;
+                        break;
+                    }
+                }
+                
+                // If not close to a vertex, check if close to any edge
+                if (!point_added) {
+                    for(size_t i = 0; i < subject[0].size()-1; i++){
+                        if (Clipper2Lib::PerpendicDistFromLineSqrd(p, subject[0][i], subject[0][i+1]) < tolerance*tolerance){
+                            simplified_paths[0].push_back(p);
+                            break;
+                        }
+                    }
+                    // Also check the closing edge between last and first point
+                    if (subject[0].size() > 1){
+                        size_t last = subject[0].size()-1;
+                        if (Clipper2Lib::PerpendicDistFromLineSqrd(p, subject[0][last], subject[0][0]) < tolerance*tolerance){
+                            simplified_paths[0].push_back(p);
+                        }
+                    }
+                }
+                }
+
+            if (simplified_paths[0].size() > 2)
+                solutions.push_back(simplified_paths);
+        }
         
-        solutions.push_back(solution);
+
     }
 
     std::vector<std::array<double, 3>> unique_points;
@@ -305,6 +348,7 @@ std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, std::vector<bool>
     groups.reserve(solutions.size() + patterns_to_keep.size());
     is_boundary.reserve(solutions.size() + patterns_to_keep.size());
     
+    // 9 bucket search to avoid rounding wrong rounded grid key, when the distance is close to the tolerance.
     auto find_or_add_point = [&](double x, double y) -> int {
         double z = 0.0;
         
@@ -323,7 +367,7 @@ std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, std::vector<bool>
                     // Check all points in this cell
                     for (int idx : grid_it->second) {
                         const auto& pt = unique_points[idx];
-                        if (is_same_point(pt[0], pt[1], pt[2], x, y, z, tolerance))
+                        if (is_same_point(pt[0], pt[1], x, y, tolerance))
                             return idx;
                     }
                 }
@@ -378,12 +422,15 @@ std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, std::vector<bool>
 
 }
 
+
+
 std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, compas::RowMatrixXd, std::vector<bool>, std::vector<int>> map_mesh_with_automatic_parameterization(
     Eigen::Ref<const compas::RowMatrixXd> target_v, 
     Eigen::Ref<const compas::RowMatrixXi> target_f, 
     Eigen::Ref<compas::RowMatrixXd> pattern_v, 
     const std::vector<std::vector<int>>& pattern_f,
     bool clip_boundaries,
+    bool simplify_borders,
     double tolerance)
 {
 
@@ -395,6 +442,35 @@ std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, compas::RowMatrix
     Eigen::VectorXi B;
     igl::boundary_loop(target_f, B);
 
+    // Get Boundary vertices
+    std::vector<Eigen::Vector3d> corner_vertex_coords;
+    
+    std::vector<std::vector<int>> VV;
+    igl::adjacency_list(target_f, VV);
+    
+    // Create a set of boundary vertices for faster lookup
+    std::unordered_set<int> boundary_set;
+    for (int i = 0; i < B.size(); i++) {
+        boundary_set.insert(B(i));
+    }
+    
+    for (int i = 0; i < B.size(); i++) {
+        int v = B(i);
+        int boundary_valence = 0;
+        for (int neighbor : VV[v])
+            if (boundary_set.count(neighbor) > 0)
+                boundary_valence++;
+        
+        if (boundary_valence < 2) {
+            Eigen::Vector3d vertex_pos(target_v(v, 0), target_v(v, 1), target_v(v, 2));
+            corner_vertex_coords.push_back(vertex_pos);
+            std::cout << "Corner vertex " << v << " at position (" 
+                      << target_v(v, 0) << ", " 
+                      << target_v(v, 1) << ", " 
+                      << target_v(v, 2) << ")" << std::endl;
+        }
+    }
+
     // Fix two points on the boundary
     Eigen::VectorXi fixed(2, 1);
     fixed(0) = B(0);
@@ -415,10 +491,10 @@ std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, compas::RowMatrix
     double scale_factor = 1.0 / std::max(size(0), size(1)); // Scale to fit in a 0-1 box while maintaining aspect ratio
     target_uv *= scale_factor;
 
-    // Clip the pattern
-    auto [clipped_pattern_v, clipped_pattern_f, clipped_pattern_is_boundary, clipped_pattern_groups] = eigen_to_clipper(target_uv, target_f, pattern_v, pattern_f, clip_boundaries, tolerance);
+    
 
-    // return std::make_tuple(clipped_pattern_v, clipped_pattern_f); // Comment this out to see 2d cropped pattern
+    // Clip the pattern
+    auto [clipped_pattern_v, clipped_pattern_f, clipped_pattern_is_boundary, clipped_pattern_groups] = eigen_to_clipper(target_uv, target_f, pattern_v, pattern_f, clip_boundaries, simplify_borders, tolerance);
 
     Eigen::MatrixXd clipped_pattern_uv;
     clipped_pattern_uv.setZero();
@@ -438,15 +514,13 @@ std::tuple<compas::RowMatrixXd, std::vector<std::vector<int>>, compas::RowMatrix
         clipped_pattern_uv,
         pattern_normals);
 
-
     return std::make_tuple(clipped_pattern_v, result, pattern_normals, clipped_pattern_is_boundary, clipped_pattern_groups);
 
 }
 
-
+// Define the nanobind module at global scope
 NB_MODULE(_mapping, m)
 {
-        
     m.def(
         "map_mesh_with_automatic_parameterization",
         &map_mesh_with_automatic_parameterization,
@@ -456,5 +530,7 @@ NB_MODULE(_mapping, m)
         "pattern_v"_a,
         "pattern_f"_a,
         "clip_boundaries"_a,
-        "tolerance"_a);
+        "simplify_borders"_a,
+        "tolerance"_a
+    );
 }
diff --git a/src/mapping.hpp b/src/mapping.hpp
