From 7c1b8a31f2d2b53e0a9720708b8764d9c28ae8bd Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Wed, 4 Mar 2026 22:46:22 -0500 Subject: [PATCH 01/10] compiles, write unit tests and check every case --- .../ImageSimulationMeshTri.cpp | 376 +++++++++++++++++- .../ImageSimulationMeshTri.hpp | 35 +- .../image_simulation/image_simulation.cpp | 26 +- .../image_simulation_spec.hpp | 42 +- .../image_simulation_spec.json | 42 +- 5 files changed, 513 insertions(+), 8 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index 498167c05a..64630ee4e4 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -1523,9 +1524,380 @@ std::tuple ImageSimulationMeshTri::get_max_avg_energy() return std::make_tuple(max_energy, avg_energy); } -void ImageSimulationMeshTri::fill_holes_topo() +std::vector +ImageSimulationMeshTri::compute_connected_components() const { - log_and_throw_error("fill_holes_topo() is not implemented"); + const size_t n_faces = tri_capacity(); + std::vector comp_id(n_faces, -1); + std::vector components; + + // Pass 1: BFS to assign component ids and collect faces per component + for (size_t fid = 0; fid < n_faces; ++fid) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + if (comp_id[fid] != -1) continue; + + const int64_t tag = m_face_attribute[fid].tags[0]; + const int comp_idx = (int)components.size(); + ConnectedComponent& comp = components.emplace_back(); + if (tag == -1) log_and_throw_error(fmt::format("Face {} has invalid tag -1", fid)); + comp.tag = tag; + + comp.faces.push_back(fid); + comp_id[fid] = comp_idx; + + for (size_t i = 0; i < comp.faces.size(); ++i) { // BFS loop with comp.faces growing + const size_t cur = comp.faces[i]; + const auto fvs = oriented_tri_vids(cur); + const Vector2d& fp0 = m_vertex_attribute[fvs[0]].m_pos; + const Vector2d& fp1 = m_vertex_attribute[fvs[1]].m_pos; + const Vector2d& fp2 = m_vertex_attribute[fvs[2]].m_pos; + comp.area += 0.5 * std::abs( + (fp1[0] - fp0[0]) * (fp2[1] - fp0[1]) - + (fp1[1] - fp0[1]) * (fp2[0] - fp0[0])); + for (int j = 0; j < 3; ++j) { + const Tuple edge_tup = tuple_from_edge(cur, j); + const auto t_opp = edge_tup.switch_face(*this); + if (!t_opp) continue; + const size_t nbr = t_opp->fid(*this); + if (m_face_attribute[nbr].tags[0] == tag && comp_id[nbr] == -1) { + comp_id[nbr] = comp_idx; + comp.faces.push_back(nbr); + } + } + } + } + + // Pass 2: populate surrounding_comp_ids by scanning all boundary edges + for (size_t fid = 0; fid < n_faces; ++fid) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + const size_t cidx = comp_id[fid]; + for (int j = 0; j < 3; ++j) { + const Tuple edge_tup = tuple_from_edge(fid, j); + const auto t_opp = edge_tup.switch_face(*this); + if (!t_opp) { + components[cidx].touches_boundary = true; + continue; + } + const size_t nbr_cidx = comp_id[t_opp->fid(*this)]; + if (nbr_cidx != cidx) { + components[cidx].surrounding_comp_ids.insert(nbr_cidx); + } + } + } + + return components; +} + +void ImageSimulationMeshTri::engulf_component( + std::vector& components, + const size_t hole_comp_id, const size_t engulfing_comp_id) +{ if (hole_comp_id >= components.size() || engulfing_comp_id >= components.size()) { + log_and_throw_error("Invalid component ids for engulfing: hole_comp_id = {}, engulfing_comp_id = {}, total components = {}", + hole_comp_id, engulfing_comp_id, components.size()); + } + auto& hole_comp = components[hole_comp_id]; + auto& engulfing_comp = components[engulfing_comp_id]; + + for (const size_t fid : hole_comp.faces) { + m_face_attribute[fid].tags[0] = engulfing_comp.tag; + } + engulfing_comp.area += hole_comp.area; + engulfing_comp.faces.insert( + engulfing_comp.faces.end(), hole_comp.faces.begin(), hole_comp.faces.end()); + engulfing_comp.surrounding_comp_ids.erase(hole_comp_id); + hole_comp = ConnectedComponent(); +} + +void ImageSimulationMeshTri::engulf_components( + std::vector& components, + const std::vector& hole_comp_ids, const std::unordered_set& engulfing_comp_ids) +{ + + + if (engulfing_comp_ids.empty()) { + log_and_throw_error("No engulfing components found for hole component(s) {}", fmt::join(hole_comp_ids, ",")); + } + else if (engulfing_comp_ids.size() == 1) { + for (const size_t comp_id : hole_comp_ids) { + engulf_component(components, comp_id, *engulfing_comp_ids.begin()); + } + return; + } + else { + // Build face->component mapping from current components state + const size_t n_faces = tri_capacity(); + std::vector comp_id_map(n_faces, std::numeric_limits::max()); + for (size_t i = 0; i < components.size(); ++i) { + for (const size_t fid : components[i].faces) { + comp_id_map[fid] = i; + } + } + + // Build boundary between hole components and engulfing components + std::vector> boundary_edge_centroids; + for (const size_t comp_id : hole_comp_ids) { + auto& comp = components[comp_id]; + for (const size_t fid : comp.faces) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + for (int j = 0; j < 3; ++j) { + const Tuple edge_tup = tuple_from_edge(fid, j); + const auto t_opp = edge_tup.switch_face(*this); + if (!t_opp) continue; + const size_t nbr_cidx = comp_id_map[t_opp->fid(*this)]; + if (engulfing_comp_ids.count(nbr_cidx)) { // neighbor is one of the engulfing components + const auto vs = get_edge_vids(edge_tup); + const Vector2d& p0 = m_vertex_attribute[vs[0]].m_pos; + const Vector2d& p1 = m_vertex_attribute[vs[1]].m_pos; + boundary_edge_centroids.emplace_back((p0 + p1) / 2.0, nbr_cidx); + } + } + } + } + + // For each hole face find nearest boundary edge centroid and assign to corresponding engulfing component + for (const size_t comp_id : hole_comp_ids) { + auto& comp = components[comp_id]; + for (const size_t fid : comp.faces) { + const auto fvs = oriented_tri_vids(fid); + const Vector2d& fp0 = m_vertex_attribute[fvs[0]].m_pos; + const Vector2d& fp1 = m_vertex_attribute[fvs[1]].m_pos; + const Vector2d& fp2 = m_vertex_attribute[fvs[2]].m_pos; + const Vector2d face_centroid = (fp0 + fp1 + fp2) / 3.0; + + double min_dist = std::numeric_limits::max(); + size_t closest_comp_id = std::numeric_limits::max(); + for (const auto& [edge_centroid, nbr_cidx] : boundary_edge_centroids) { + const double dist = (face_centroid - edge_centroid).norm(); + if (dist < min_dist) { + min_dist = dist; + closest_comp_id = nbr_cidx; + } + } + if (closest_comp_id == std::numeric_limits::max()) { + log_and_throw_error("Failed to find closest component for face {}", fid); + } + m_face_attribute[fid].tags[0] = components[closest_comp_id].tag; + components[closest_comp_id].area += 0.5 * std::abs( + (fp1[0] - fp0[0]) * (fp2[1] - fp0[1]) - + (fp1[1] - fp0[1]) * (fp2[0] - fp0[0])); + components[closest_comp_id].faces.push_back(fid); + comp_id_map[fid] = closest_comp_id; + } + } + // Update surrounding_comp_ids for newly merged components + for (const size_t comp_id : hole_comp_ids) { + auto& comp = components[comp_id]; + for (const size_t fid : comp.faces) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + const size_t cidx = comp_id_map[fid]; + for (int j = 0; j < 3; ++j) { + const Tuple edge_tup = tuple_from_edge(fid, j); + const auto t_opp = edge_tup.switch_face(*this); + if (!t_opp) continue; + const size_t nbr_cidx = comp_id_map[t_opp->fid(*this)]; + if (nbr_cidx == std::numeric_limits::max()) continue; + if (nbr_cidx == comp_id) continue; + if (nbr_cidx != cidx) { + components[cidx].surrounding_comp_ids.insert(nbr_cidx); + components[nbr_cidx].surrounding_comp_ids.insert(cidx); + } + } + } + } + + for (const size_t comp_id : hole_comp_ids) { + for (auto& c : components) { + c.surrounding_comp_ids.erase(comp_id); + } + components[comp_id] = ConnectedComponent(); + } + } +} + +void ImageSimulationMeshTri::extract_hole_clusters(std::vector& components, std::unordered_set& tags, std::vector>& hole_clusters, double threshold) +{ + // BFS-cluster all non-fill-tag components by adjacency. + // Each cluster is a maximal group of connected non-fill-tag components. + // A cluster is enclosed if every surrounding component outside it is in tags. + std::vector cluster_id(components.size(), -1); + std::vector component_clusters; + for (size_t i = 0; i < components.size(); ++i) { + if (components[i].faces.empty() || tags.count(components[i].tag) || components[i].touches_boundary) continue; + if (cluster_id[i] != -1) continue; + const size_t cid = component_clusters.size(); + ComponentCluster& cluster = component_clusters.emplace_back(); + cluster.comp_ids.push_back(i); + cluster_id[i] = cid; + for (size_t qi = 0; qi < cluster.comp_ids.size(); ++qi) { // BFS loop with cluster.comp_ids growing + const size_t cur = cluster.comp_ids[qi]; + cluster.area += components[cur].area; + for (const size_t sid : components[cur].surrounding_comp_ids) { + if (components[sid].faces.empty()) continue; + if (tags.count(components[sid].tag)) continue; + if (components[sid].touches_boundary) { + cluster.enclosed = false; // neighbour escapes to mesh boundary + continue; + } + if (cluster_id[sid] == -1) { + cluster_id[sid] = cid; + cluster.comp_ids.push_back(sid); + } + } + } + // also mark not enclosed if any member touches a non-fill-tag component outside this cluster + for (const size_t ci : cluster.comp_ids) { + for (const size_t sid : components[ci].surrounding_comp_ids) { + if (!components[sid].faces.empty() && + !tags.count(components[sid].tag) && + cluster_id[sid] != cid) { + cluster.enclosed = false; + } + } + } + } + + // fill enclosed clusters whose total area is below threshold + for (const ComponentCluster& cluster : component_clusters) { + if (!cluster.enclosed) continue; + if (cluster.area >= threshold) continue; + hole_clusters.push_back({}); + for (const size_t ci : cluster.comp_ids) { + hole_clusters.back().push_back(ci); + } + } +} + +void ImageSimulationMeshTri::recompute_surface_info() +{ + for (const Tuple& e : get_edges()) { + SmartTuple ee(*this, e); + m_edge_attribute[ee.eid()].m_is_surface_fs = 0; + } + for (size_t vid = 0; vid < vert_capacity(); ++vid) { + m_vertex_attribute[vid].m_is_on_surface = false; + } + for (const Tuple& e : get_edges()) { + SmartTuple ee(*this, e); + const auto t_opp = ee.switch_face(); + if (!t_opp) continue; + bool has_diff_tag = false; + for (size_t j = 0; j < m_tags_count; ++j) { + if (m_face_attribute[ee.fid()].tags[j] != m_face_attribute[t_opp->fid()].tags[j]) { + has_diff_tag = true; + break; + } + } + if (!has_diff_tag) continue; + m_edge_attribute[ee.eid()].m_is_surface_fs = 1; + const size_t v1 = ee.vid(); + const size_t v2 = ee.switch_vertex().vid(); + m_vertex_attribute[v1].m_is_on_surface = true; + m_vertex_attribute[v2].m_is_on_surface = true; + } +} +void ImageSimulationMeshTri::fill_holes_topo( + const std::vector& fill_holes_tags, + double threshold) +{ + if (m_tags_count == 0) { + logger().warn("fill_holes_topo: no tags, skipping"); + return; + } + + // -- Step 1: compute all connected components with their surrounding component ids + auto components = compute_connected_components(); + + // -- Step 2: for each fill_tag, fill all enclosed components + bool any_filled = false; + for (const int64_t fill_tag : fill_holes_tags) { + + std::vector> hole_clusters; + std::unordered_set tags{fill_tag}; + extract_hole_clusters(components, tags, hole_clusters, threshold); + + for (const auto& hole_cluster : hole_clusters) { + // Collect the fill-tag component indices that surround this cluster + std::unordered_set engulfing_comp_ids; + for (const size_t ci : hole_cluster) { + for (const size_t sid : components[ci].surrounding_comp_ids) { + if (!components[sid].faces.empty() && components[sid].tag == fill_tag) { + engulfing_comp_ids.insert(sid); + } + } + } + engulf_components(components, hole_cluster, engulfing_comp_ids); + any_filled = true; + } + if (hole_clusters.empty()) { + logger().info( + "fill_holes_topo: no enclosed components found for fill_tag {}", fill_tag); + } + } + + if (!any_filled) return; + + // -- Step 3: recompute m_is_surface_fs and m_is_on_surface + recompute_surface_info(); +} +void ImageSimulationMeshTri::keep_largest_connected_component(const std::vector& lcc_tags) +{ + auto components = compute_connected_components(); + if (components.empty()) { + logger().warn("keep_largest_connected_component: no components found, skipping"); + return; + } + // For each tag in lcc_tags, find the largest component with that tag and mark all other components with that tag for removal + for (const int64_t lcc_tag : lcc_tags) { + int largest_comp_id = -1; + double largest_area = -1.0; + for (size_t i = 0; i < components.size(); ++i) { + if (components[i].faces.empty() || components[i].tag != lcc_tag) continue; + if (components[i].area > largest_area) { + largest_area = components[i].area; + largest_comp_id = i; + } + } + if (largest_comp_id == -1) { + logger().warn( + "keep_largest_connected_component: no component found for tag {}, skipping", + lcc_tag); + continue; + } + for (size_t i = 0; i < components.size(); ++i) { + if (components[i].faces.empty() || components[i].tag != lcc_tag || (int)i == largest_comp_id) + continue; + // Copy surrounding_comp_ids — engulf_components will reset components[i] + std::unordered_set surr = components[i].surrounding_comp_ids; + engulf_components(components, std::vector{i}, surr); + } + } + + recompute_surface_info(); +} + +void ImageSimulationMeshTri::tight_seal_topo(const std::vector>& tight_seal_tag_sets, double threshold) +{ + auto components = compute_connected_components(); + + for (const auto& tag_set : tight_seal_tag_sets) { + std::vector> hole_clusters; + std::unordered_set tags_copy(tag_set.begin(), tag_set.end()); + extract_hole_clusters(components, tags_copy, hole_clusters, threshold); + for (const auto& hole_cluster : hole_clusters) { + // Collect component indices in tag_set that surround this cluster + std::unordered_set engulfing_comp_ids; + for (const size_t ci : hole_cluster) { + for (const size_t sid : components[ci].surrounding_comp_ids) { + if (!components[sid].faces.empty() && tag_set.count(components[sid].tag)) { + engulfing_comp_ids.insert(sid); + } + } + } + engulf_components(components, hole_cluster, engulfing_comp_ids); + } + } + + recompute_surface_info(); } simplex::RawSimplexCollection ImageSimulationMeshTri::get_order1_edges_for_vertex( diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp index 0a1caf4bc5..2ae0cb7d5b 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp @@ -1,5 +1,8 @@ #pragma once +#include +#include + #include #include #include @@ -196,7 +199,37 @@ class ImageSimulationMeshTri : public wmtk::TriMesh bool collapse_limit_length = true); std::tuple get_max_avg_energy(); - void fill_holes_topo(); + // Returns: (face ids in component, set of distinct tags of all neighboring faces across tag boundaries) - defined according to tag_0 + struct ConnectedComponent + { + int64_t tag = -1; + std::vector faces; + std::unordered_set surrounding_comp_ids; + double area = 0.0; + bool touches_boundary = false; + }; + std::vector compute_connected_components() const; + void engulf_component(std::vector& components, const size_t hole_comp_id, const size_t engulfing_comp_id); + void engulf_components(std::vector& components, const std::vector& hole_comp_ids, const std::unordered_set& engulfing_comp_ids); + void keep_largest_connected_component(const std::vector& lcc_tags); + + // A cluster is a maximal group of connected non-fill-tag ConnectedComponents. + // enclosed = true iff every surrounding component outside the cluster has fill_tag. + struct ComponentCluster + { + std::vector comp_ids; // indices into components vector + double area = 0.0; + bool enclosed = true; + }; + void extract_hole_clusters(std::vector& components, std::unordered_set& tags, std::vector>& hole_clusters, double threshold); + + void recompute_surface_info(); + void fill_holes_topo( + const std::vector& fill_holes_tags, + double threshold = std::numeric_limits::infinity()); + + void tight_seal_topo(const std::vector>& tight_seal_tag_sets, double threshold = std::numeric_limits::infinity()); + /** * @brief Get all edges on the surface that are incident to vid. diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp index b12db9d87f..02ec049927 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp @@ -174,7 +174,31 @@ void run_2D(const nlohmann::json& json_params, const InputData& input_data) if (operation == "remeshing") { mesh.mesh_improvement(max_its); // <-- tetwild } else if (operation == "fill_holes_topo") { - mesh.fill_holes_topo(); + const std::vector fill_holes_tags = + json_params["fill_holes_tags"].get>(); + const double raw_threshold = json_params["fill_holes_threshold"].get(); + const double threshold = + raw_threshold < 0 ? std::numeric_limits::infinity() : raw_threshold; + mesh.fill_holes_topo(fill_holes_tags, threshold); + } else if (operation == "tight_seal_topo") { + // Accept either [[t1,t2],[t3,...]] (list of sets) or flat [t1,t2] (one set) + std::vector> tag_sets; + const auto& raw = json_params["tight_seal_tag_sets"]; + if (!raw.empty() && raw[0].is_array()) { + for (const auto& s : raw) { + std::unordered_set ts; + for (const auto& v : s) ts.insert(v.get()); + tag_sets.push_back(std::move(ts)); + } + } else { + std::unordered_set ts; + for (const auto& v : raw) ts.insert(v.get()); + tag_sets.push_back(std::move(ts)); + } + const double raw_threshold = json_params["tight_seal_threshold"].get(); + const double threshold = + raw_threshold < 0 ? std::numeric_limits::infinity() : raw_threshold; + mesh.tight_seal_topo(tag_sets, threshold); } else { log_and_throw_error("Unknown image simulation operation"); } diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp index f72379d359..9c9e193817 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp @@ -26,7 +26,11 @@ nlohmann::json image_simulation_spec = R"( "log_file", "report", "DEBUG_output", - "DEBUG_sanity_checks" + "DEBUG_sanity_checks", + "fill_holes_tags", + "fill_holes_threshold", + "tight_seal_tag_sets", + "tight_seal_threshold" ] }, { @@ -50,7 +54,7 @@ nlohmann::json image_simulation_spec = R"( "pointer": "/operation", "type": "string", "default": "remeshing", - "options": ["remeshing", "fill_holes_topo"], + "options": ["remeshing", "fill_holes_topo", "tight_seal_topo"], "doc": "Image simulation contains multiple operations for modifying the image. Depending on the operation, more parameters might be required." }, { @@ -171,6 +175,40 @@ nlohmann::json image_simulation_spec = R"( "type": "bool", "default": false, "doc": "Perform sanity checks after every operation. This can be very slow and should only be used for debugging." + }, + { + "pointer": "/fill_holes_tags", + "type": "list", + "default": [0], + "doc": "For fill_holes_topo: list of tag values used to fill enclosed connected components (processed in order)." + }, + { + "pointer": "/fill_holes_tags/*", + "type": "int", + "doc": "A tag value to fill enclosed connected components of other tags." + }, + { + "pointer": "/fill_holes_threshold", + "type": "float", + "default": -1, + "doc": "For fill_holes_topo: only fill a connected component if its area is less than this threshold. Negative value means no threshold (fill all enclosed components)." + }, + { + "pointer": "/tight_seal_tag_sets", + "type": "list", + "default": [[0]], + "doc": "For tight_seal_topo: list of tag sets. Each set defines a group of tags whose enclosed holes are filled. Can be a flat list [t1,t2] (treated as one set) or a list of lists [[t1,t2],[t3]]." + }, + { + "pointer": "/tight_seal_tag_sets/*", + "type": "int", + "doc": "A tag value in the tight seal tag set (when tight_seal_tag_sets is a flat list of ints)." + }, + { + "pointer": "/tight_seal_threshold", + "type": "float", + "default": -1, + "doc": "For tight_seal_topo: only fill a hole cluster if its total area is less than this threshold. Negative value means no threshold." } ] )"_json; diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json index b2847521f9..2adc919a1b 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json @@ -21,7 +21,11 @@ "log_file", "report", "DEBUG_output", - "DEBUG_sanity_checks" + "DEBUG_sanity_checks", + "fill_holes_tags", + "fill_holes_threshold", + "tight_seal_tag_sets", + "tight_seal_threshold" ] }, { @@ -45,7 +49,7 @@ "pointer": "/operation", "type": "string", "default": "remeshing", - "options": ["remeshing", "fill_holes_topo"], + "options": ["remeshing", "fill_holes_topo", "tight_seal_topo"], "doc": "Image simulation contains multiple operations for modifying the image. Depending on the operation, more parameters might be required." }, { @@ -166,5 +170,39 @@ "type": "bool", "default": false, "doc": "Perform sanity checks after every operation. This can be very slow and should only be used for debugging." + }, + { + "pointer": "/fill_holes_tags", + "type": "list", + "default": [0], + "doc": "For fill_holes_topo: list of tag values used to fill enclosed connected components (processed in order)." + }, + { + "pointer": "/fill_holes_tags/*", + "type": "int", + "doc": "A tag value to fill enclosed connected components of other tags." + }, + { + "pointer": "/fill_holes_threshold", + "type": "float", + "default": -1, + "doc": "For fill_holes_topo: only fill a connected component if its area is less than this threshold. Negative value means no threshold (fill all enclosed components)." + }, + { + "pointer": "/tight_seal_tag_sets", + "type": "list", + "default": [[0]], + "doc": "For tight_seal_topo: list of tag sets. Each set defines a group of tags whose enclosed holes are filled. Can be a flat list [t1,t2] (treated as one set) or a list of lists [[t1,t2],[t3]]." + }, + { + "pointer": "/tight_seal_tag_sets/*", + "type": "int", + "doc": "A tag value in the tight seal tag set (when tight_seal_tag_sets is a flat list of ints)." + }, + { + "pointer": "/tight_seal_threshold", + "type": "float", + "default": -1, + "doc": "For tight_seal_topo: only fill a hole cluster if its total area is less than this threshold. Negative value means no threshold." } ] From da5fed50a2f4b28b801da8f18dd764260d145fa2 Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Fri, 6 Mar 2026 15:40:42 -0500 Subject: [PATCH 02/10] tight seal remeshing workin. write tests and setup CI --- .../ConnectedComponentAnnotationHelper.inl | 559 ++++++++++++++++++ .../ImageSimulationMeshTri.cpp | 320 ++-------- .../ImageSimulationMeshTri.hpp | 5 + .../image_simulation/image_simulation.cpp | 18 +- .../image_simulation_spec.hpp | 25 +- .../image_simulation_spec.json | 9 +- 6 files changed, 645 insertions(+), 291 deletions(-) create mode 100644 components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl diff --git a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl new file mode 100644 index 0000000000..c7f3a84941 --- /dev/null +++ b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl @@ -0,0 +1,559 @@ +std::vector +ImageSimulationMeshTri::compute_connected_components() const +{ + const size_t n_faces = tri_capacity(); + std::vector comp_id(n_faces, -1); + std::vector components; + + // Pass 1: BFS to assign component ids and collect faces per component + for (size_t fid = 0; fid < n_faces; ++fid) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + if (comp_id[fid] != -1) continue; + + const int64_t tag = m_face_attribute[fid].tags[0]; + const int comp_idx = (int)components.size(); + ConnectedComponent& comp = components.emplace_back(); + if (tag == -1) log_and_throw_error(fmt::format("Face {} has invalid tag -1", fid)); + comp.tag = tag; + + comp.faces.push_back(fid); + comp_id[fid] = comp_idx; + + for (size_t i = 0; i < comp.faces.size(); ++i) { // BFS loop with comp.faces growing + const size_t cur = comp.faces[i]; + const auto fvs = oriented_tri_vids(cur); + const Vector2d& fp0 = m_vertex_attribute[fvs[0]].m_pos; + const Vector2d& fp1 = m_vertex_attribute[fvs[1]].m_pos; + const Vector2d& fp2 = m_vertex_attribute[fvs[2]].m_pos; + comp.area += 0.5 * std::abs( + (fp1[0] - fp0[0]) * (fp2[1] - fp0[1]) - + (fp1[1] - fp0[1]) * (fp2[0] - fp0[0])); + for (int j = 0; j < 3; ++j) { + const Tuple edge_tup = tuple_from_edge(cur, j); + const auto t_opp = edge_tup.switch_face(*this); + if (!t_opp) continue; + const size_t nbr = t_opp->fid(*this); + if (m_face_attribute[nbr].tags[0] == tag && comp_id[nbr] == -1) { + comp_id[nbr] = comp_idx; + comp.faces.push_back(nbr); + } + } + } + } + + // Pass 2: populate surrounding_comp_ids by scanning all boundary edges + for (size_t fid = 0; fid < n_faces; ++fid) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + const size_t cidx = comp_id[fid]; + for (int j = 0; j < 3; ++j) { + const Tuple edge_tup = tuple_from_edge(fid, j); + const auto t_opp = edge_tup.switch_face(*this); + if (!t_opp) { + components[cidx].touches_boundary = true; + continue; + } + const size_t nbr_cidx = comp_id[t_opp->fid(*this)]; + if (nbr_cidx != cidx) { + components[cidx].surrounding_comp_ids.insert(nbr_cidx); + } + } + } + + return components; +} + +void ImageSimulationMeshTri::engulf_component( + std::vector& components, + const size_t hole_comp_id, const size_t engulfing_comp_id) +{ if (hole_comp_id >= components.size() || engulfing_comp_id >= components.size()) { + log_and_throw_error("Invalid component ids for engulfing: hole_comp_id = {}, engulfing_comp_id = {}, total components = {}", + hole_comp_id, engulfing_comp_id, components.size()); + } + auto& hole_comp = components[hole_comp_id]; + auto& engulfing_comp = components[engulfing_comp_id]; + + for (const size_t fid : hole_comp.faces) { + m_face_attribute[fid].tags[0] = engulfing_comp.tag; + } + engulfing_comp.area += hole_comp.area; + engulfing_comp.faces.insert( + engulfing_comp.faces.end(), hole_comp.faces.begin(), hole_comp.faces.end()); + engulfing_comp.surrounding_comp_ids.erase(hole_comp_id); + hole_comp = ConnectedComponent(); +} + +void ImageSimulationMeshTri::engulf_components( + std::vector& components, + const std::vector& hole_comp_ids, const std::unordered_set& engulfing_comp_ids) +{ + + + if (engulfing_comp_ids.empty()) { + log_and_throw_error("No engulfing components found for hole component(s) {}", fmt::join(hole_comp_ids, ",")); + } + else if (engulfing_comp_ids.size() == 1) { + for (const size_t comp_id : hole_comp_ids) { + engulf_component(components, comp_id, *engulfing_comp_ids.begin()); + } + return; + } + else { + // Build face->component mapping from current components state + std::vector comp_id_map(tri_capacity(), std::numeric_limits::max()); + for (size_t i = 0; i < components.size(); ++i) { + for (const size_t fid : components[i].faces) { + comp_id_map[fid] = i; + } + } + + // Build boundary edge centroids between hole and engulfing components + std::vector> boundary_edge_centroids; + for (const size_t comp_id : hole_comp_ids) { + for (const size_t fid : components[comp_id].faces) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + for (int j = 0; j < 3; ++j) { + const Tuple edge_tup = tuple_from_edge(fid, j); + const auto t_opp = edge_tup.switch_face(*this); + if (!t_opp) continue; + const size_t nbr_cidx = comp_id_map[t_opp->fid(*this)]; + if (engulfing_comp_ids.count(nbr_cidx)) { // neighbor is one of the engulfing components + const auto vs = get_edge_vids(edge_tup); + const Vector2d& p0 = m_vertex_attribute[vs[0]].m_pos; + const Vector2d& p1 = m_vertex_attribute[vs[1]].m_pos; + boundary_edge_centroids.emplace_back((p0 + p1) / 2.0, nbr_cidx); + } + } + } + } + + // Helper: find closest engulfing comp for a given vertex position + const auto compute_vertex_closest_comp = [&](const Vector2d& pos) -> size_t { + double min_dist = std::numeric_limits::max(); + size_t closest = std::numeric_limits::max(); + for (const auto& [centroid, nbr_cidx] : boundary_edge_centroids) { + const double dist = (pos - centroid).squaredNorm(); + if (dist < min_dist) { + min_dist = dist; + closest = nbr_cidx; + } + } + return closest; + }; + + // Compute vertex_closest_comp for all vertices of hole faces + const std::unordered_set hole_comp_ids_set(hole_comp_ids.begin(), hole_comp_ids.end()); + std::unordered_map vertex_closest_comp; // vid -> closest engulfing comp_id + for (const size_t comp_id : hole_comp_ids) { + for (const size_t fid : components[comp_id].faces) { + const auto fvs = oriented_tri_vids(fid); + for (const size_t vid : fvs) { + if (vertex_closest_comp.count(vid)) continue; + vertex_closest_comp[vid] = compute_vertex_closest_comp(m_vertex_attribute[vid].m_pos); + } + } + } + + // Find hole edges where the two endpoint vertices disagree on closest comp -> split. + // Deduplicate by sorted vertex pair to avoid splitting same edge twice. + std::set> edges_to_split_set; + for (const size_t comp_id : hole_comp_ids) { + for (const size_t fid : components[comp_id].faces) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + for (int j = 0; j < 3; ++j) { + const Tuple edge_tup = tuple_from_edge(fid, j); + const auto vs = get_edge_vids(edge_tup); + const size_t va = vs[0], vb = vs[1]; + if (vertex_closest_comp.count(va) && vertex_closest_comp.count(vb)) { + const size_t ca = vertex_closest_comp.at(va); + const size_t cb = vertex_closest_comp.at(vb); + // Skip if either endpoint is equidistant + // would have no centroids for that comp and produce garbage values. + if (ca != std::numeric_limits::max() && + cb != std::numeric_limits::max() && ca != cb) { + edges_to_split_set.emplace(std::min(va, vb), std::max(va, vb)); + } + } + } + } + } + + // Helper: signed distance between two engulfing comps at position p. + // Negative = closer to comp_a centroids, positive = closer to comp_b centroids. + const auto voronoi_sign = [&](const Vector2d& p, const size_t comp_a, const size_t comp_b) -> double { + double min_a = std::numeric_limits::max(); + double min_b = std::numeric_limits::max(); + for (const auto& [centroid, cidx] : boundary_edge_centroids) { + const double d = (p - centroid).squaredNorm(); + if (cidx == comp_a) min_a = std::min(min_a, d); + if (cidx == comp_b) min_b = std::min(min_b, d); + } + return min_a - min_b; + }; + + // Split disagreement edges; split_edge_after binary-searches vmid onto the + // equidistant surface between the two competing engulfing comps via m_voronoi_split_fn. + for (const auto& [va, vb] : edges_to_split_set) { + auto [edge_tup, _eid] = tuple_from_edge({{va, vb}}); + if (!edge_tup.is_valid(*this)) continue; + const size_t comp_a = vertex_closest_comp.at(va); + const size_t comp_b = vertex_closest_comp.at(vb); + // Skip if either endpoint was relabeled equidistant by a previous iteration's + // near-endpoint check — voronoi_sign with equidistant comp returns garbage. + if (comp_a == std::numeric_limits::max() || + comp_b == std::numeric_limits::max()) + continue; + + // If the Voronoi boundary is so close to va or vb that vmid would land on top of + // it, skip the split — it would create a degenerate face. + // Use linear interpolation: t = |sign_a| / (|sign_a| + |sign_b|) is the fraction + // along the edge where the zero-crossing lies. + // Also relabel the near-boundary vertex as equidistant so assignment is correct. + { + const double sign_a = voronoi_sign(m_vertex_attribute[va].m_pos, comp_a, comp_b); + const double sign_b = voronoi_sign(m_vertex_attribute[vb].m_pos, comp_a, comp_b); + const double total = std::abs(sign_a) + std::abs(sign_b); + constexpr double kMinFrac = 0.01; + if (total == 0.0) continue; // both equidistant + const double t = std::abs(sign_a) / total; + if (t < kMinFrac) { vertex_closest_comp[va] = std::numeric_limits::max(); continue; } + if (t > 1.0 - kMinFrac) { vertex_closest_comp[vb] = std::numeric_limits::max(); continue; } + } + + m_voronoi_split_fn = [&, comp_a, comp_b](const Vector2d& p) -> double { + return voronoi_sign(p, comp_a, comp_b); + }; + std::vector new_tris; + const bool split_success = split_edge(edge_tup, new_tris); + m_voronoi_split_fn = nullptr; + if (!split_success || new_tris.empty()) { + if (!split_success) logger().warn( + "engulf_components: split_edge({},{}) failed; assignment will use centroid fallback", + va, vb); + } else { + const size_t vmid = new_tris[0].vid(*this); + vertex_closest_comp[vmid] = std::numeric_limits::max(); // equidistant + + // Extend comp_id_map and register genuinely new hole faces + if (tri_capacity() > comp_id_map.size()) + comp_id_map.resize(tri_capacity(), std::numeric_limits::max()); + for (const Tuple& nt : new_tris) { + const size_t nfid = nt.fid(*this); + if (comp_id_map[nfid] != std::numeric_limits::max()) continue; + const int64_t ftag = m_face_attribute[nfid].tags[0]; + for (const size_t hcid : hole_comp_ids) { + if (components[hcid].tag == ftag) { + comp_id_map[nfid] = hcid; + components[hcid].faces.push_back(nfid); + break; + } + } + } + } + } + + // Assign all hole faces to engulfing comp. + // Boundary (equidistant) vertices have vertex_closest_comp == max — skip them. + // All definite vertices of a face must agree on the comp; if all are boundary, skip face. + std::unordered_map> newly_assigned; // engulfing cidx -> fids + for (const size_t comp_id : hole_comp_ids) { + for (const size_t fid : components[comp_id].faces) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + const auto fvs = oriented_tri_vids(fid); + // Ensure all vertices have an entry (new faces may have new verts) + for (const size_t vid : fvs) { + if (!vertex_closest_comp.count(vid)) + vertex_closest_comp[vid] = + compute_vertex_closest_comp(m_vertex_attribute[vid].m_pos); + } + // Collect definite (non-boundary) vertex assignments + size_t closest = std::numeric_limits::max(); + bool conflict = false; + for (const size_t vid : fvs) { + const size_t c = vertex_closest_comp.at(vid); + if (c == std::numeric_limits::max()) continue; // equidistant vertex + if (closest == std::numeric_limits::max()) { + closest = c; + } else if (closest != c) { + conflict = true; + break; + } + } + if (conflict || closest == std::numeric_limits::max()) { + const Vector2d centroid = + (m_vertex_attribute[fvs[0]].m_pos + + m_vertex_attribute[fvs[1]].m_pos + + m_vertex_attribute[fvs[2]].m_pos) / 3.0; + closest = compute_vertex_closest_comp(centroid); + if (conflict) logger().warn("Face {} has conflicting vertex assignments after split; using centroid fallback for assignment",fid); + if (closest == std::numeric_limits::max()) logger().warn("Face {} has no definite vertex assignments after split; using centroid fallback for assignment", fid); + } + + const auto& fp0 = m_vertex_attribute[fvs[0]].m_pos; + const auto& fp1 = m_vertex_attribute[fvs[1]].m_pos; + const auto& fp2 = m_vertex_attribute[fvs[2]].m_pos; + m_face_attribute[fid].tags[0] = components[closest].tag; + components[closest].area += 0.5 * std::abs( + (fp1[0] - fp0[0]) * (fp2[1] - fp0[1]) - + (fp1[1] - fp0[1]) * (fp2[0] - fp0[0])); + components[closest].faces.push_back(fid); + comp_id_map[fid] = closest; + newly_assigned[closest].push_back(fid); + } + } + + // Update surrounding_comp_ids using newly assigned faces (not stale comp.faces) + for (const auto& [cidx, fids] : newly_assigned) { + for (const size_t fid : fids) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + for (int j = 0; j < 3; ++j) { + const Tuple edge_tup = tuple_from_edge(fid, j); + const auto t_opp = edge_tup.switch_face(*this); + if (!t_opp) continue; + const size_t nbr_cidx = comp_id_map[t_opp->fid(*this)]; + if (nbr_cidx == std::numeric_limits::max()) continue; + if (nbr_cidx == cidx) continue; + if (hole_comp_ids_set.count(nbr_cidx)) continue; // will be erased + components[cidx].surrounding_comp_ids.insert(nbr_cidx); + components[nbr_cidx].surrounding_comp_ids.insert(cidx); + } + } + } + + // --- DEBUG: export mesh with per-vertex distance and label data --- + { + static int debug_call = 0; + const std::string debug_path = fmt::format("/tmp/engulf_debug_{}.vtu", debug_call++); + + // Collect all valid vertices referenced by hole faces after splitting + std::map vid_to_row; + std::vector row_vids; + for (const size_t comp_id : hole_comp_ids) { + for (const size_t fid : components[comp_id].faces) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + for (const size_t v : oriented_tri_vids(fid)) { + if (!vid_to_row.count(v)) { + vid_to_row[v] = (int)row_vids.size(); + row_vids.push_back(v); + } + } + } + } + + const int nv = (int)row_vids.size(); + Eigen::MatrixXd V(nv, 3); + Eigen::VectorXd closest_label(nv), voronoi_sdf(nv); + V.setZero(); closest_label.setZero(); voronoi_sdf.setZero(); + + for (int i = 0; i < nv; ++i) { + const size_t vid = row_vids[i]; + const Vector2d& p = m_vertex_attribute[vid].m_pos; + V(i, 0) = p[0]; V(i, 1) = p[1]; V(i, 2) = 0.0; + + // Find the two closest comps by per-comp minimum squared distance — + // same computation as closest_comp / voronoi_sign uses. + std::unordered_map comp_min_sq; + for (const auto& [centroid, cidx] : boundary_edge_centroids) { + const double d = (p - centroid).squaredNorm(); + auto [it, ins] = comp_min_sq.emplace(cidx, d); + if (!ins) it->second = std::min(it->second, d); + } + size_t c1 = std::numeric_limits::max(); + size_t c2 = std::numeric_limits::max(); + double d1sq = std::numeric_limits::max(); + double d2sq = std::numeric_limits::max(); + for (const auto& [cidx, dsq] : comp_min_sq) { + if (dsq < d1sq) { d2sq = d1sq; c2 = c1; d1sq = dsq; c1 = cidx; } + else if (dsq < d2sq) { d2sq = dsq; c2 = cidx; } + } + // voronoi_sign(p, c1, c2) = min_sq_dist_c1 - min_sq_dist_c2 + // This is exactly what is used in the splitting λ, so values here are + // directly comparable to the binary-search target (zero crossing). + voronoi_sdf[i] = (c2 != std::numeric_limits::max()) + ? voronoi_sign(p, c1, c2) + : 0.0; + + // Label: closest comp from vertex_closest_comp (-1 = equidistant) + const size_t c = vertex_closest_comp.count(vid) + ? vertex_closest_comp.at(vid) + : std::numeric_limits::max(); + closest_label[i] = (c == std::numeric_limits::max()) ? -1.0 : (double)c; + } + + // Build face connectivity + std::vector> face_list; + for (const size_t comp_id : hole_comp_ids) { + for (const size_t fid : components[comp_id].faces) { + if (!tuple_from_tri(fid).is_valid(*this)) continue; + const auto fvs = oriented_tri_vids(fid); + face_list.push_back({(size_t)vid_to_row.at(fvs[0]), + (size_t)vid_to_row.at(fvs[1]), + (size_t)vid_to_row.at(fvs[2])}); + } + } + Eigen::MatrixXi F((int)face_list.size(), 3); + Eigen::VectorXd face_area((int)face_list.size()); + for (int i = 0; i < (int)face_list.size(); ++i) { + F(i, 0) = (int)face_list[i][0]; + F(i, 1) = (int)face_list[i][1]; + F(i, 2) = (int)face_list[i][2]; + const double ax = V(F(i,0),0), ay = V(F(i,0),1); + const double bx = V(F(i,1),0), by = V(F(i,1),1); + const double cx = V(F(i,2),0), cy = V(F(i,2),1); + face_area[i] = 0.5 * std::abs((bx-ax)*(cy-ay) - (by-ay)*(cx-ax)); + } + // Check for zero/near-zero area faces + constexpr double kAreaTol = 1e-10; + int zero_area_count = 0; + for (int i = 0; i < (int)face_list.size(); ++i) { + if (face_area[i] < kAreaTol) { + ++zero_area_count; + wmtk::logger().warn( + "DEBUG engulf_components: face row {} (verts [{},{},{}]) " + "has near-zero area = {} | pos [{},{}] [{},{}] [{},{}]", + i, F(i,0), F(i,1), F(i,2), face_area[i], + V(F(i,0),0), V(F(i,0),1), + V(F(i,1),0), V(F(i,1),1), + V(F(i,2),0), V(F(i,2),1)); + } + } + if (zero_area_count > 0) + wmtk::logger().warn( + "DEBUG engulf_components: {} / {} hole faces have near-zero area", + zero_area_count, face_list.size()); + + // Verify: equidistant-labeled vertices should have voronoi_sdf ≈ 0 for SOME + // component pair — the pair used when splitting may differ from the globally + // nearest two, so we check the minimum |voronoi_sign| over all pairs. + for (int i = 0; i < nv; ++i) { + if (closest_label[i] != -1.0) continue; + // Build per-component minimum squared distances + std::unordered_map comp_min_sq; + for (const auto& [centroid, cidx] : boundary_edge_centroids) { + const double d = (m_vertex_attribute[row_vids[i]].m_pos - centroid).squaredNorm(); + auto [it, ins] = comp_min_sq.emplace(cidx, d); + if (!ins) it->second = std::min(it->second, d); + } + if (comp_min_sq.size() < 2) continue; + // Find the nearest centroid distance (denominator for rel) + double d1sq = std::numeric_limits::max(); + for (const auto& [cidx, dsq] : comp_min_sq) + d1sq = std::min(d1sq, dsq); + // Find the pair (ca, cb) that minimizes |voronoi_sign| = |d_ca - d_cb| + double min_abs_sign = std::numeric_limits::max(); + size_t best_ca = 0, best_cb = 0; + const std::vector> comp_vec(comp_min_sq.begin(), comp_min_sq.end()); + for (size_t ai = 0; ai < comp_vec.size(); ++ai) + for (size_t bi = ai + 1; bi < comp_vec.size(); ++bi) { + const double s = std::abs(comp_vec[ai].second - comp_vec[bi].second); + if (s < min_abs_sign) { min_abs_sign = s; best_ca = comp_vec[ai].first; best_cb = comp_vec[bi].first; } + } + constexpr double kMinFrac = 0.01; + const double rel = (d1sq > 0) ? min_abs_sign / d1sq : min_abs_sign; + if (rel > kMinFrac) { + const size_t vid = row_vids[i]; + wmtk::logger().warn( + "DEBUG engulf_components: vertex {} labeled equidistant " + "but best pair voronoi_sign(ca={},cb={})={:.4f} rel={:.4f} | pos [{}, {}]", + vid, best_ca, best_cb, min_abs_sign, rel, V(i, 0), V(i, 1)); + } + } + + paraviewo::VTUWriter dbg_writer; + dbg_writer.add_field("closest_comp", closest_label); + dbg_writer.add_field("voronoi_sdf", voronoi_sdf); + dbg_writer.add_cell_field("face_area", face_area); + dbg_writer.write_mesh(debug_path, V, F); + wmtk::logger().info("engulf_components debug mesh written to {}", debug_path); + } + // --- END DEBUG --- + + for (const size_t comp_id : hole_comp_ids) { + for (auto& c : components) { + c.surrounding_comp_ids.erase(comp_id); + } + components[comp_id] = ConnectedComponent(); + } + } +} + +void ImageSimulationMeshTri::extract_hole_clusters(std::vector& components, std::unordered_set& tags, std::vector>& hole_clusters, double threshold) +{ + // BFS-cluster all non-fill-tag components by adjacency. + // Each cluster is a maximal group of connected non-fill-tag components. + // A cluster is enclosed if every surrounding component outside it is in tags. + std::vector cluster_id(components.size(), -1); + std::vector component_clusters; + for (size_t i = 0; i < components.size(); ++i) { + if (components[i].faces.empty() || tags.count(components[i].tag) || components[i].touches_boundary) continue; + if (cluster_id[i] != -1) continue; + const size_t cid = component_clusters.size(); + ComponentCluster& cluster = component_clusters.emplace_back(); + cluster.comp_ids.push_back(i); + cluster_id[i] = cid; + for (size_t qi = 0; qi < cluster.comp_ids.size(); ++qi) { // BFS loop with cluster.comp_ids growing + const size_t cur = cluster.comp_ids[qi]; + cluster.area += components[cur].area; + for (const size_t sid : components[cur].surrounding_comp_ids) { + if (components[sid].faces.empty()) continue; + if (tags.count(components[sid].tag)) continue; + if (components[sid].touches_boundary) { + cluster.enclosed = false; // neighbour escapes to mesh boundary + continue; + } + if (cluster_id[sid] == -1) { + cluster_id[sid] = cid; + cluster.comp_ids.push_back(sid); + } + } + } + // also mark not enclosed if any member touches a non-fill-tag component outside this cluster + for (const size_t ci : cluster.comp_ids) { + for (const size_t sid : components[ci].surrounding_comp_ids) { + if (!components[sid].faces.empty() && + !tags.count(components[sid].tag) && + cluster_id[sid] != cid) { + cluster.enclosed = false; + } + } + } + } + + // fill enclosed clusters whose total area is below threshold + for (const ComponentCluster& cluster : component_clusters) { + if (!cluster.enclosed) continue; + if (cluster.area >= threshold) continue; + hole_clusters.push_back({}); + for (const size_t ci : cluster.comp_ids) { + hole_clusters.back().push_back(ci); + } + } +} + +void ImageSimulationMeshTri::recompute_surface_info() +{ + for (const Tuple& e : get_edges()) { + SmartTuple ee(*this, e); + m_edge_attribute[ee.eid()].m_is_surface_fs = 0; + } + for (size_t vid = 0; vid < vert_capacity(); ++vid) { + m_vertex_attribute[vid].m_is_on_surface = false; + } + for (const Tuple& e : get_edges()) { + SmartTuple ee(*this, e); + const auto t_opp = ee.switch_face(); + if (!t_opp) continue; + bool has_diff_tag = false; + for (size_t j = 0; j < m_tags_count; ++j) { + if (m_face_attribute[ee.fid()].tags[j] != m_face_attribute[t_opp->fid()].tags[j]) { + has_diff_tag = true; + break; + } + } + if (!has_diff_tag) continue; + m_edge_attribute[ee.eid()].m_is_surface_fs = 1; + const size_t v1 = ee.vid(); + const size_t v2 = ee.switch_vertex().vid(); + m_vertex_attribute[v1].m_is_on_surface = true; + m_vertex_attribute[v2].m_is_on_surface = true; + } +} diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index 64630ee4e4..a353f27c81 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -9,6 +9,8 @@ #include #include #include +#include +#include #include #include #include @@ -603,8 +605,8 @@ bool ImageSimulationMeshTri::split_edge_after(const Tuple& loc) const size_t v2_id = cache.v2_id; /// check inversion & rounding - m_vertex_attribute[v_id].m_pos = - (m_vertex_attribute[v1_id].m_pos + m_vertex_attribute[v2_id].m_pos) / 2; + auto& p = m_vertex_attribute[v_id].m_pos; + p = (m_vertex_attribute[v1_id].m_pos + m_vertex_attribute[v2_id].m_pos) / 2; for (const Tuple& t : locs) { if (is_inverted(t)) { @@ -612,6 +614,33 @@ bool ImageSimulationMeshTri::split_edge_after(const Tuple& loc) } } + // If a Voronoi split function is set, binary-search vmid onto its zero-crossing. + // p0 stays on the negative side, p1 on the positive side. + if (m_voronoi_split_fn) { + Vector2d p0 = m_vertex_attribute[v1_id].m_pos; + Vector2d p1 = m_vertex_attribute[v2_id].m_pos; + if (m_voronoi_split_fn(p0) >= 0) std::swap(p0, p1); // ensure p0 is negative side + for (int i = 0; i < 20; ++i) { + p = 0.5 * (p0 + p1); + bool inv = false; + for (const Tuple& t : locs) { + if (is_inverted(t)) { inv = true; break; } + } + if (inv || (p1 - p0).squaredNorm() < 1e-20) break; + if (m_voronoi_split_fn(p) < 0) p0 = p; + else p1 = p; + } + // final inversion guard: revert to midpoint if needed + bool inv = false; + for (const Tuple& t : locs) { + if (is_inverted(t)) { + inv = true; + logger().warn("Voronoi split resulted in inversion, reverting to midpoint. Iteration: {}", debug_print_counter++); + break; } + } + if (inv) p = (m_vertex_attribute[v1_id].m_pos + m_vertex_attribute[v2_id].m_pos) / 2; + } + // update face attributes { // v1 - v_new @@ -1524,280 +1553,9 @@ std::tuple ImageSimulationMeshTri::get_max_avg_energy() return std::make_tuple(max_energy, avg_energy); } -std::vector -ImageSimulationMeshTri::compute_connected_components() const -{ - const size_t n_faces = tri_capacity(); - std::vector comp_id(n_faces, -1); - std::vector components; - - // Pass 1: BFS to assign component ids and collect faces per component - for (size_t fid = 0; fid < n_faces; ++fid) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; - if (comp_id[fid] != -1) continue; - - const int64_t tag = m_face_attribute[fid].tags[0]; - const int comp_idx = (int)components.size(); - ConnectedComponent& comp = components.emplace_back(); - if (tag == -1) log_and_throw_error(fmt::format("Face {} has invalid tag -1", fid)); - comp.tag = tag; - - comp.faces.push_back(fid); - comp_id[fid] = comp_idx; - - for (size_t i = 0; i < comp.faces.size(); ++i) { // BFS loop with comp.faces growing - const size_t cur = comp.faces[i]; - const auto fvs = oriented_tri_vids(cur); - const Vector2d& fp0 = m_vertex_attribute[fvs[0]].m_pos; - const Vector2d& fp1 = m_vertex_attribute[fvs[1]].m_pos; - const Vector2d& fp2 = m_vertex_attribute[fvs[2]].m_pos; - comp.area += 0.5 * std::abs( - (fp1[0] - fp0[0]) * (fp2[1] - fp0[1]) - - (fp1[1] - fp0[1]) * (fp2[0] - fp0[0])); - for (int j = 0; j < 3; ++j) { - const Tuple edge_tup = tuple_from_edge(cur, j); - const auto t_opp = edge_tup.switch_face(*this); - if (!t_opp) continue; - const size_t nbr = t_opp->fid(*this); - if (m_face_attribute[nbr].tags[0] == tag && comp_id[nbr] == -1) { - comp_id[nbr] = comp_idx; - comp.faces.push_back(nbr); - } - } - } - } - - // Pass 2: populate surrounding_comp_ids by scanning all boundary edges - for (size_t fid = 0; fid < n_faces; ++fid) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; - const size_t cidx = comp_id[fid]; - for (int j = 0; j < 3; ++j) { - const Tuple edge_tup = tuple_from_edge(fid, j); - const auto t_opp = edge_tup.switch_face(*this); - if (!t_opp) { - components[cidx].touches_boundary = true; - continue; - } - const size_t nbr_cidx = comp_id[t_opp->fid(*this)]; - if (nbr_cidx != cidx) { - components[cidx].surrounding_comp_ids.insert(nbr_cidx); - } - } - } - - return components; -} - -void ImageSimulationMeshTri::engulf_component( - std::vector& components, - const size_t hole_comp_id, const size_t engulfing_comp_id) -{ if (hole_comp_id >= components.size() || engulfing_comp_id >= components.size()) { - log_and_throw_error("Invalid component ids for engulfing: hole_comp_id = {}, engulfing_comp_id = {}, total components = {}", - hole_comp_id, engulfing_comp_id, components.size()); - } - auto& hole_comp = components[hole_comp_id]; - auto& engulfing_comp = components[engulfing_comp_id]; - - for (const size_t fid : hole_comp.faces) { - m_face_attribute[fid].tags[0] = engulfing_comp.tag; - } - engulfing_comp.area += hole_comp.area; - engulfing_comp.faces.insert( - engulfing_comp.faces.end(), hole_comp.faces.begin(), hole_comp.faces.end()); - engulfing_comp.surrounding_comp_ids.erase(hole_comp_id); - hole_comp = ConnectedComponent(); -} - -void ImageSimulationMeshTri::engulf_components( - std::vector& components, - const std::vector& hole_comp_ids, const std::unordered_set& engulfing_comp_ids) -{ - - - if (engulfing_comp_ids.empty()) { - log_and_throw_error("No engulfing components found for hole component(s) {}", fmt::join(hole_comp_ids, ",")); - } - else if (engulfing_comp_ids.size() == 1) { - for (const size_t comp_id : hole_comp_ids) { - engulf_component(components, comp_id, *engulfing_comp_ids.begin()); - } - return; - } - else { - // Build face->component mapping from current components state - const size_t n_faces = tri_capacity(); - std::vector comp_id_map(n_faces, std::numeric_limits::max()); - for (size_t i = 0; i < components.size(); ++i) { - for (const size_t fid : components[i].faces) { - comp_id_map[fid] = i; - } - } - - // Build boundary between hole components and engulfing components - std::vector> boundary_edge_centroids; - for (const size_t comp_id : hole_comp_ids) { - auto& comp = components[comp_id]; - for (const size_t fid : comp.faces) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; - for (int j = 0; j < 3; ++j) { - const Tuple edge_tup = tuple_from_edge(fid, j); - const auto t_opp = edge_tup.switch_face(*this); - if (!t_opp) continue; - const size_t nbr_cidx = comp_id_map[t_opp->fid(*this)]; - if (engulfing_comp_ids.count(nbr_cidx)) { // neighbor is one of the engulfing components - const auto vs = get_edge_vids(edge_tup); - const Vector2d& p0 = m_vertex_attribute[vs[0]].m_pos; - const Vector2d& p1 = m_vertex_attribute[vs[1]].m_pos; - boundary_edge_centroids.emplace_back((p0 + p1) / 2.0, nbr_cidx); - } - } - } - } - - // For each hole face find nearest boundary edge centroid and assign to corresponding engulfing component - for (const size_t comp_id : hole_comp_ids) { - auto& comp = components[comp_id]; - for (const size_t fid : comp.faces) { - const auto fvs = oriented_tri_vids(fid); - const Vector2d& fp0 = m_vertex_attribute[fvs[0]].m_pos; - const Vector2d& fp1 = m_vertex_attribute[fvs[1]].m_pos; - const Vector2d& fp2 = m_vertex_attribute[fvs[2]].m_pos; - const Vector2d face_centroid = (fp0 + fp1 + fp2) / 3.0; - - double min_dist = std::numeric_limits::max(); - size_t closest_comp_id = std::numeric_limits::max(); - for (const auto& [edge_centroid, nbr_cidx] : boundary_edge_centroids) { - const double dist = (face_centroid - edge_centroid).norm(); - if (dist < min_dist) { - min_dist = dist; - closest_comp_id = nbr_cidx; - } - } - if (closest_comp_id == std::numeric_limits::max()) { - log_and_throw_error("Failed to find closest component for face {}", fid); - } - m_face_attribute[fid].tags[0] = components[closest_comp_id].tag; - components[closest_comp_id].area += 0.5 * std::abs( - (fp1[0] - fp0[0]) * (fp2[1] - fp0[1]) - - (fp1[1] - fp0[1]) * (fp2[0] - fp0[0])); - components[closest_comp_id].faces.push_back(fid); - comp_id_map[fid] = closest_comp_id; - } - } - // Update surrounding_comp_ids for newly merged components - for (const size_t comp_id : hole_comp_ids) { - auto& comp = components[comp_id]; - for (const size_t fid : comp.faces) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; - const size_t cidx = comp_id_map[fid]; - for (int j = 0; j < 3; ++j) { - const Tuple edge_tup = tuple_from_edge(fid, j); - const auto t_opp = edge_tup.switch_face(*this); - if (!t_opp) continue; - const size_t nbr_cidx = comp_id_map[t_opp->fid(*this)]; - if (nbr_cidx == std::numeric_limits::max()) continue; - if (nbr_cidx == comp_id) continue; - if (nbr_cidx != cidx) { - components[cidx].surrounding_comp_ids.insert(nbr_cidx); - components[nbr_cidx].surrounding_comp_ids.insert(cidx); - } - } - } - } - - for (const size_t comp_id : hole_comp_ids) { - for (auto& c : components) { - c.surrounding_comp_ids.erase(comp_id); - } - components[comp_id] = ConnectedComponent(); - } - } -} - -void ImageSimulationMeshTri::extract_hole_clusters(std::vector& components, std::unordered_set& tags, std::vector>& hole_clusters, double threshold) -{ - // BFS-cluster all non-fill-tag components by adjacency. - // Each cluster is a maximal group of connected non-fill-tag components. - // A cluster is enclosed if every surrounding component outside it is in tags. - std::vector cluster_id(components.size(), -1); - std::vector component_clusters; - for (size_t i = 0; i < components.size(); ++i) { - if (components[i].faces.empty() || tags.count(components[i].tag) || components[i].touches_boundary) continue; - if (cluster_id[i] != -1) continue; - const size_t cid = component_clusters.size(); - ComponentCluster& cluster = component_clusters.emplace_back(); - cluster.comp_ids.push_back(i); - cluster_id[i] = cid; - for (size_t qi = 0; qi < cluster.comp_ids.size(); ++qi) { // BFS loop with cluster.comp_ids growing - const size_t cur = cluster.comp_ids[qi]; - cluster.area += components[cur].area; - for (const size_t sid : components[cur].surrounding_comp_ids) { - if (components[sid].faces.empty()) continue; - if (tags.count(components[sid].tag)) continue; - if (components[sid].touches_boundary) { - cluster.enclosed = false; // neighbour escapes to mesh boundary - continue; - } - if (cluster_id[sid] == -1) { - cluster_id[sid] = cid; - cluster.comp_ids.push_back(sid); - } - } - } - // also mark not enclosed if any member touches a non-fill-tag component outside this cluster - for (const size_t ci : cluster.comp_ids) { - for (const size_t sid : components[ci].surrounding_comp_ids) { - if (!components[sid].faces.empty() && - !tags.count(components[sid].tag) && - cluster_id[sid] != cid) { - cluster.enclosed = false; - } - } - } - } +#include "ConnectedComponentAnnotationHelper.inl" - // fill enclosed clusters whose total area is below threshold - for (const ComponentCluster& cluster : component_clusters) { - if (!cluster.enclosed) continue; - if (cluster.area >= threshold) continue; - hole_clusters.push_back({}); - for (const size_t ci : cluster.comp_ids) { - hole_clusters.back().push_back(ci); - } - } -} - -void ImageSimulationMeshTri::recompute_surface_info() -{ - for (const Tuple& e : get_edges()) { - SmartTuple ee(*this, e); - m_edge_attribute[ee.eid()].m_is_surface_fs = 0; - } - for (size_t vid = 0; vid < vert_capacity(); ++vid) { - m_vertex_attribute[vid].m_is_on_surface = false; - } - for (const Tuple& e : get_edges()) { - SmartTuple ee(*this, e); - const auto t_opp = ee.switch_face(); - if (!t_opp) continue; - bool has_diff_tag = false; - for (size_t j = 0; j < m_tags_count; ++j) { - if (m_face_attribute[ee.fid()].tags[j] != m_face_attribute[t_opp->fid()].tags[j]) { - has_diff_tag = true; - break; - } - } - if (!has_diff_tag) continue; - m_edge_attribute[ee.eid()].m_is_surface_fs = 1; - const size_t v1 = ee.vid(); - const size_t v2 = ee.switch_vertex().vid(); - m_vertex_attribute[v1].m_is_on_surface = true; - m_vertex_attribute[v2].m_is_on_surface = true; - } -} -void ImageSimulationMeshTri::fill_holes_topo( - const std::vector& fill_holes_tags, - double threshold) +void ImageSimulationMeshTri::fill_holes_topo(const std::vector& fill_holes_tags, double threshold) { if (m_tags_count == 0) { logger().warn("fill_holes_topo: no tags, skipping"); @@ -1886,13 +1644,25 @@ void ImageSimulationMeshTri::tight_seal_topo(const std::vector engulfing_comp_ids; + std::unordered_set surrounding_tags; for (const size_t ci : hole_cluster) { for (const size_t sid : components[ci].surrounding_comp_ids) { if (!components[sid].faces.empty() && tag_set.count(components[sid].tag)) { engulfing_comp_ids.insert(sid); + surrounding_tags.insert(components[sid].tag); } } } + // check surrounding tags is the same as tag_set + if (surrounding_tags != tag_set) continue; // harmless hole, skip + std::cout<<"Engulfing_comp_ids: "; + for (const auto& id : engulfing_comp_ids) { + std::cout< face_tags; }; tbb::enumerable_thread_specific swap_cache; + + // When set, split_edge_after binary-searches vmid onto the zero-crossing of this function. + // Negative = stays on v1 side, positive = stays on v2 side. + // Set before split_edge(), cleared immediately after. + std::function m_voronoi_split_fn = nullptr; }; } // namespace wmtk::components::image_simulation::tri diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp index 02ec049927..1aa69e0a25 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp @@ -181,24 +181,21 @@ void run_2D(const nlohmann::json& json_params, const InputData& input_data) raw_threshold < 0 ? std::numeric_limits::infinity() : raw_threshold; mesh.fill_holes_topo(fill_holes_tags, threshold); } else if (operation == "tight_seal_topo") { - // Accept either [[t1,t2],[t3,...]] (list of sets) or flat [t1,t2] (one set) + // tight_seal_tag_sets is a list of lists: [[t1,t2],[t3,...]] std::vector> tag_sets; - const auto& raw = json_params["tight_seal_tag_sets"]; - if (!raw.empty() && raw[0].is_array()) { - for (const auto& s : raw) { - std::unordered_set ts; - for (const auto& v : s) ts.insert(v.get()); - tag_sets.push_back(std::move(ts)); - } - } else { + for (const auto& s : json_params["tight_seal_tag_sets"]) { std::unordered_set ts; - for (const auto& v : raw) ts.insert(v.get()); + for (const auto& v : s) ts.insert(v.get()); tag_sets.push_back(std::move(ts)); } const double raw_threshold = json_params["tight_seal_threshold"].get(); const double threshold = raw_threshold < 0 ? std::numeric_limits::infinity() : raw_threshold; mesh.tight_seal_topo(tag_sets, threshold); + } else if (operation == "keep_largest_cc") { + const std::vector lcc_tags = + json_params["keep_largest_cc_tags"].get>(); + mesh.keep_largest_connected_component(lcc_tags); } else { log_and_throw_error("Unknown image simulation operation"); } @@ -266,6 +263,7 @@ void image_simulation(nlohmann::json json_params) output_filename.string()); } output_filename.replace_extension(""); // extension is added back later + json_params["output"] = output_filename.string(); // propagate resolved path to run_2D/run_3D auto get_unique_vtu_name = [&output_filename]() -> std::string { static size_t vtu_counter = 0; diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp index 9c9e193817..c2833c4f2e 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp @@ -30,7 +30,8 @@ nlohmann::json image_simulation_spec = R"( "fill_holes_tags", "fill_holes_threshold", "tight_seal_tag_sets", - "tight_seal_threshold" + "tight_seal_threshold", + "keep_largest_cc_tags" ] }, { @@ -54,7 +55,7 @@ nlohmann::json image_simulation_spec = R"( "pointer": "/operation", "type": "string", "default": "remeshing", - "options": ["remeshing", "fill_holes_topo", "tight_seal_topo"], + "options": ["remeshing", "fill_holes_topo", "tight_seal_topo", "keep_largest_cc"], "doc": "Image simulation contains multiple operations for modifying the image. Depending on the operation, more parameters might be required." }, { @@ -197,18 +198,34 @@ nlohmann::json image_simulation_spec = R"( "pointer": "/tight_seal_tag_sets", "type": "list", "default": [[0]], - "doc": "For tight_seal_topo: list of tag sets. Each set defines a group of tags whose enclosed holes are filled. Can be a flat list [t1,t2] (treated as one set) or a list of lists [[t1,t2],[t3]]." + "doc": "For tight_seal_topo: list of tag sets. Each inner list defines a group of tags whose enclosed holes are filled, e.g. [[1,2],[3]]." }, { "pointer": "/tight_seal_tag_sets/*", + "type": "list", + "doc": "One tag set: a list of tag values that are treated as a group for hole filling." + }, + { + "pointer": "/tight_seal_tag_sets/*/*", "type": "int", - "doc": "A tag value in the tight seal tag set (when tight_seal_tag_sets is a flat list of ints)." + "doc": "A tag value in this tag set." }, { "pointer": "/tight_seal_threshold", "type": "float", "default": -1, "doc": "For tight_seal_topo: only fill a hole cluster if its total area is less than this threshold. Negative value means no threshold." + }, + { + "pointer": "/keep_largest_cc_tags", + "type": "list", + "default": [0], + "doc": "For keep_largest_cc: list of tag values for which only the largest connected component is kept; all smaller components are merged into their neighbours." + }, + { + "pointer": "/keep_largest_cc_tags/*", + "type": "int", + "doc": "A tag value whose smaller connected components will be removed." } ] )"_json; diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json index 2adc919a1b..c17a454a1c 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json @@ -192,12 +192,17 @@ "pointer": "/tight_seal_tag_sets", "type": "list", "default": [[0]], - "doc": "For tight_seal_topo: list of tag sets. Each set defines a group of tags whose enclosed holes are filled. Can be a flat list [t1,t2] (treated as one set) or a list of lists [[t1,t2],[t3]]." + "doc": "For tight_seal_topo: list of tag sets. Each inner list defines a group of tags whose enclosed holes are filled, e.g. [[1,2],[3]]." }, { "pointer": "/tight_seal_tag_sets/*", + "type": "list", + "doc": "One tag set: a list of tag values that are treated as a group for hole filling." + }, + { + "pointer": "/tight_seal_tag_sets/*/*", "type": "int", - "doc": "A tag value in the tight seal tag set (when tight_seal_tag_sets is a flat list of ints)." + "doc": "A tag value in this tag set." }, { "pointer": "/tight_seal_threshold", From bfac9225f8960127766aa50dd7395877ae849d64 Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Fri, 13 Mar 2026 15:45:27 -0400 Subject: [PATCH 03/10] image sim topo annot 2d CI, pending data upload --- .../ConnectedComponentAnnotationHelper.inl | 90 +++++++++---------- .../image_simulation_spec.json | 16 +++- 2 files changed, 56 insertions(+), 50 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl index c7f3a84941..766406ac06 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl +++ b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl @@ -126,16 +126,24 @@ void ImageSimulationMeshTri::engulf_components( } } - // Helper: find closest engulfing comp for a given vertex position + // Helper: per-comp minimum squared distance from pos to any centroid of that comp. + const auto comp_min_sq_dists = [&](const Vector2d& pos) { + std::unordered_map result; + for (const auto& [centroid, cidx] : boundary_edge_centroids) { + const double d = (pos - centroid).squaredNorm(); + auto [it, ins] = result.emplace(cidx, d); + if (!ins) it->second = std::min(it->second, d); + } + return result; + }; + + // Helper: find closest engulfing comp for a given vertex position. const auto compute_vertex_closest_comp = [&](const Vector2d& pos) -> size_t { - double min_dist = std::numeric_limits::max(); + const auto dists = comp_min_sq_dists(pos); size_t closest = std::numeric_limits::max(); - for (const auto& [centroid, nbr_cidx] : boundary_edge_centroids) { - const double dist = (pos - centroid).squaredNorm(); - if (dist < min_dist) { - min_dist = dist; - closest = nbr_cidx; - } + double min_dist = std::numeric_limits::max(); + for (const auto& [cidx, d] : dists) { + if (d < min_dist) { min_dist = d; closest = cidx; } } return closest; }; @@ -166,13 +174,13 @@ void ImageSimulationMeshTri::engulf_components( if (vertex_closest_comp.count(va) && vertex_closest_comp.count(vb)) { const size_t ca = vertex_closest_comp.at(va); const size_t cb = vertex_closest_comp.at(vb); - // Skip if either endpoint is equidistant - // would have no centroids for that comp and produce garbage values. + // Skip if either endpoint is equidistant if (ca != std::numeric_limits::max() && cb != std::numeric_limits::max() && ca != cb) { edges_to_split_set.emplace(std::min(va, vb), std::max(va, vb)); } } + else log_and_throw_error(fmt::format("Vertex {} or {} has no closest component; this should not happen", va, vb)); } } } @@ -180,14 +188,10 @@ void ImageSimulationMeshTri::engulf_components( // Helper: signed distance between two engulfing comps at position p. // Negative = closer to comp_a centroids, positive = closer to comp_b centroids. const auto voronoi_sign = [&](const Vector2d& p, const size_t comp_a, const size_t comp_b) -> double { - double min_a = std::numeric_limits::max(); - double min_b = std::numeric_limits::max(); - for (const auto& [centroid, cidx] : boundary_edge_centroids) { - const double d = (p - centroid).squaredNorm(); - if (cidx == comp_a) min_a = std::min(min_a, d); - if (cidx == comp_b) min_b = std::min(min_b, d); - } - return min_a - min_b; + const auto dists = comp_min_sq_dists(p); + const double da = dists.count(comp_a) ? dists.at(comp_a) : std::numeric_limits::max(); + const double db = dists.count(comp_b) ? dists.at(comp_b) : std::numeric_limits::max(); + return da - db; }; // Split disagreement edges; split_edge_after binary-searches vmid onto the @@ -211,9 +215,18 @@ void ImageSimulationMeshTri::engulf_components( { const double sign_a = voronoi_sign(m_vertex_attribute[va].m_pos, comp_a, comp_b); const double sign_b = voronoi_sign(m_vertex_attribute[vb].m_pos, comp_a, comp_b); - const double total = std::abs(sign_a) + std::abs(sign_b); - constexpr double kMinFrac = 0.01; - if (total == 0.0) continue; // both equidistant + + if (sign_a * sign_b > 0.0) log_and_throw_error(fmt::format("Endpoints {} and {} of edge to split have same sign {}, {}: this should not happen", va, vb, sign_a, sign_b)); + if (sign_a > 0.0 || sign_b < 0.0) log_and_throw_error(fmt::format("Endpoint signs for edge to split are sign_a = {}, sign_b = {}: this should not happen", sign_a, sign_b)); + + const double total = sign_b - sign_a; + if (total < 1e-10) // both comps are very close to equidistant — skip split and mark both endpoints as equidistant to avoid numerical issues with voronoi_sign near the boundary + { + vertex_closest_comp[va] = std::numeric_limits::max(); + vertex_closest_comp[vb] = std::numeric_limits::max(); + continue; + } + constexpr double kMinFrac = 1e-3; // manual tolerance to avoid degenerate splits; if the boundary is within 0.1% of the edge length from an endpoint, skip the split and mark that endpoint as equidistant. const double t = std::abs(sign_a) / total; if (t < kMinFrac) { vertex_closest_comp[va] = std::numeric_limits::max(); continue; } if (t > 1.0 - kMinFrac) { vertex_closest_comp[vb] = std::numeric_limits::max(); continue; } @@ -319,7 +332,7 @@ void ImageSimulationMeshTri::engulf_components( } } - // --- DEBUG: export mesh with per-vertex distance and label data --- + // --- DEBUG: { static int debug_call = 0; const std::string debug_path = fmt::format("/tmp/engulf_debug_{}.vtu", debug_call++); @@ -349,14 +362,8 @@ void ImageSimulationMeshTri::engulf_components( const Vector2d& p = m_vertex_attribute[vid].m_pos; V(i, 0) = p[0]; V(i, 1) = p[1]; V(i, 2) = 0.0; - // Find the two closest comps by per-comp minimum squared distance — - // same computation as closest_comp / voronoi_sign uses. - std::unordered_map comp_min_sq; - for (const auto& [centroid, cidx] : boundary_edge_centroids) { - const double d = (p - centroid).squaredNorm(); - auto [it, ins] = comp_min_sq.emplace(cidx, d); - if (!ins) it->second = std::min(it->second, d); - } + // Find the two closest comps using the shared helper. + const auto comp_min_sq = comp_min_sq_dists(p); size_t c1 = std::numeric_limits::max(); size_t c2 = std::numeric_limits::max(); double d1sq = std::numeric_limits::max(); @@ -407,7 +414,7 @@ void ImageSimulationMeshTri::engulf_components( for (int i = 0; i < (int)face_list.size(); ++i) { if (face_area[i] < kAreaTol) { ++zero_area_count; - wmtk::logger().warn( + wmtk::log_and_throw_error( "DEBUG engulf_components: face row {} (verts [{},{},{}]) " "has near-zero area = {} | pos [{},{}] [{},{}] [{},{}]", i, F(i,0), F(i,1), F(i,2), face_area[i], @@ -417,7 +424,7 @@ void ImageSimulationMeshTri::engulf_components( } } if (zero_area_count > 0) - wmtk::logger().warn( + wmtk::log_and_throw_error( "DEBUG engulf_components: {} / {} hole faces have near-zero area", zero_area_count, face_list.size()); @@ -426,13 +433,7 @@ void ImageSimulationMeshTri::engulf_components( // nearest two, so we check the minimum |voronoi_sign| over all pairs. for (int i = 0; i < nv; ++i) { if (closest_label[i] != -1.0) continue; - // Build per-component minimum squared distances - std::unordered_map comp_min_sq; - for (const auto& [centroid, cidx] : boundary_edge_centroids) { - const double d = (m_vertex_attribute[row_vids[i]].m_pos - centroid).squaredNorm(); - auto [it, ins] = comp_min_sq.emplace(cidx, d); - if (!ins) it->second = std::min(it->second, d); - } + const auto comp_min_sq = comp_min_sq_dists(m_vertex_attribute[row_vids[i]].m_pos); if (comp_min_sq.size() < 2) continue; // Find the nearest centroid distance (denominator for rel) double d1sq = std::numeric_limits::max(); @@ -448,22 +449,15 @@ void ImageSimulationMeshTri::engulf_components( if (s < min_abs_sign) { min_abs_sign = s; best_ca = comp_vec[ai].first; best_cb = comp_vec[bi].first; } } constexpr double kMinFrac = 0.01; - const double rel = (d1sq > 0) ? min_abs_sign / d1sq : min_abs_sign; + const double rel = (abs(d1sq) > 1e-8) ? min_abs_sign / d1sq : min_abs_sign; if (rel > kMinFrac) { const size_t vid = row_vids[i]; - wmtk::logger().warn( + wmtk::log_and_throw_error( "DEBUG engulf_components: vertex {} labeled equidistant " "but best pair voronoi_sign(ca={},cb={})={:.4f} rel={:.4f} | pos [{}, {}]", vid, best_ca, best_cb, min_abs_sign, rel, V(i, 0), V(i, 1)); } } - - paraviewo::VTUWriter dbg_writer; - dbg_writer.add_field("closest_comp", closest_label); - dbg_writer.add_field("voronoi_sdf", voronoi_sdf); - dbg_writer.add_cell_field("face_area", face_area); - dbg_writer.write_mesh(debug_path, V, F); - wmtk::logger().info("engulf_components debug mesh written to {}", debug_path); } // --- END DEBUG --- diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json index c17a454a1c..aab00f1ba6 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json @@ -25,7 +25,8 @@ "fill_holes_tags", "fill_holes_threshold", "tight_seal_tag_sets", - "tight_seal_threshold" + "tight_seal_threshold", + "keep_largest_cc_tags" ] }, { @@ -49,7 +50,7 @@ "pointer": "/operation", "type": "string", "default": "remeshing", - "options": ["remeshing", "fill_holes_topo", "tight_seal_topo"], + "options": ["remeshing", "fill_holes_topo", "tight_seal_topo", "keep_largest_cc"], "doc": "Image simulation contains multiple operations for modifying the image. Depending on the operation, more parameters might be required." }, { @@ -209,5 +210,16 @@ "type": "float", "default": -1, "doc": "For tight_seal_topo: only fill a hole cluster if its total area is less than this threshold. Negative value means no threshold." + }, + { + "pointer": "/keep_largest_cc_tags", + "type": "list", + "default": [0], + "doc": "For keep_largest_cc: list of tag values for which only the largest connected component is kept; all smaller components are merged into their neighbours." + }, + { + "pointer": "/keep_largest_cc_tags/*", + "type": "int", + "doc": "A tag value whose smaller connected components will be removed." } ] From 49ff7c4088750a2a02de192adb6e23662e615bc9 Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Fri, 13 Mar 2026 16:49:59 -0400 Subject: [PATCH 04/10] pre-commit fix --- .../ConnectedComponentAnnotationHelper.inl | 248 ++++++++++++------ 1 file changed, 171 insertions(+), 77 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl index 766406ac06..3db84dabd9 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl +++ b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl @@ -26,8 +26,8 @@ ImageSimulationMeshTri::compute_connected_components() const const Vector2d& fp1 = m_vertex_attribute[fvs[1]].m_pos; const Vector2d& fp2 = m_vertex_attribute[fvs[2]].m_pos; comp.area += 0.5 * std::abs( - (fp1[0] - fp0[0]) * (fp2[1] - fp0[1]) - - (fp1[1] - fp0[1]) * (fp2[0] - fp0[0])); + (fp1[0] - fp0[0]) * (fp2[1] - fp0[1]) - + (fp1[1] - fp0[1]) * (fp2[0] - fp0[0])); for (int j = 0; j < 3; ++j) { const Tuple edge_tup = tuple_from_edge(cur, j); const auto t_opp = edge_tup.switch_face(*this); @@ -64,10 +64,16 @@ ImageSimulationMeshTri::compute_connected_components() const void ImageSimulationMeshTri::engulf_component( std::vector& components, - const size_t hole_comp_id, const size_t engulfing_comp_id) -{ if (hole_comp_id >= components.size() || engulfing_comp_id >= components.size()) { - log_and_throw_error("Invalid component ids for engulfing: hole_comp_id = {}, engulfing_comp_id = {}, total components = {}", - hole_comp_id, engulfing_comp_id, components.size()); + const size_t hole_comp_id, + const size_t engulfing_comp_id) +{ + if (hole_comp_id >= components.size() || engulfing_comp_id >= components.size()) { + log_and_throw_error( + "Invalid component ids for engulfing: hole_comp_id = {}, engulfing_comp_id = {}, total " + "components = {}", + hole_comp_id, + engulfing_comp_id, + components.size()); } auto& hole_comp = components[hole_comp_id]; auto& engulfing_comp = components[engulfing_comp_id]; @@ -77,27 +83,28 @@ void ImageSimulationMeshTri::engulf_component( } engulfing_comp.area += hole_comp.area; engulfing_comp.faces.insert( - engulfing_comp.faces.end(), hole_comp.faces.begin(), hole_comp.faces.end()); - engulfing_comp.surrounding_comp_ids.erase(hole_comp_id); + engulfing_comp.faces.end(), + hole_comp.faces.begin(), + hole_comp.faces.end()); + engulfing_comp.surrounding_comp_ids.erase(hole_comp_id); hole_comp = ConnectedComponent(); } void ImageSimulationMeshTri::engulf_components( std::vector& components, - const std::vector& hole_comp_ids, const std::unordered_set& engulfing_comp_ids) + const std::vector& hole_comp_ids, + const std::unordered_set& engulfing_comp_ids) { - - if (engulfing_comp_ids.empty()) { - log_and_throw_error("No engulfing components found for hole component(s) {}", fmt::join(hole_comp_ids, ",")); - } - else if (engulfing_comp_ids.size() == 1) { + log_and_throw_error( + "No engulfing components found for hole component(s) {}", + fmt::join(hole_comp_ids, ",")); + } else if (engulfing_comp_ids.size() == 1) { for (const size_t comp_id : hole_comp_ids) { engulf_component(components, comp_id, *engulfing_comp_ids.begin()); } return; - } - else { + } else { // Build face->component mapping from current components state std::vector comp_id_map(tri_capacity(), std::numeric_limits::max()); for (size_t i = 0; i < components.size(); ++i) { @@ -116,7 +123,8 @@ void ImageSimulationMeshTri::engulf_components( const auto t_opp = edge_tup.switch_face(*this); if (!t_opp) continue; const size_t nbr_cidx = comp_id_map[t_opp->fid(*this)]; - if (engulfing_comp_ids.count(nbr_cidx)) { // neighbor is one of the engulfing components + if (engulfing_comp_ids.count( + nbr_cidx)) { // neighbor is one of the engulfing components const auto vs = get_edge_vids(edge_tup); const Vector2d& p0 = m_vertex_attribute[vs[0]].m_pos; const Vector2d& p1 = m_vertex_attribute[vs[1]].m_pos; @@ -143,20 +151,26 @@ void ImageSimulationMeshTri::engulf_components( size_t closest = std::numeric_limits::max(); double min_dist = std::numeric_limits::max(); for (const auto& [cidx, d] : dists) { - if (d < min_dist) { min_dist = d; closest = cidx; } + if (d < min_dist) { + min_dist = d; + closest = cidx; + } } return closest; }; // Compute vertex_closest_comp for all vertices of hole faces - const std::unordered_set hole_comp_ids_set(hole_comp_ids.begin(), hole_comp_ids.end()); + const std::unordered_set hole_comp_ids_set( + hole_comp_ids.begin(), + hole_comp_ids.end()); std::unordered_map vertex_closest_comp; // vid -> closest engulfing comp_id for (const size_t comp_id : hole_comp_ids) { for (const size_t fid : components[comp_id].faces) { const auto fvs = oriented_tri_vids(fid); for (const size_t vid : fvs) { if (vertex_closest_comp.count(vid)) continue; - vertex_closest_comp[vid] = compute_vertex_closest_comp(m_vertex_attribute[vid].m_pos); + vertex_closest_comp[vid] = + compute_vertex_closest_comp(m_vertex_attribute[vid].m_pos); } } } @@ -179,18 +193,25 @@ void ImageSimulationMeshTri::engulf_components( cb != std::numeric_limits::max() && ca != cb) { edges_to_split_set.emplace(std::min(va, vb), std::max(va, vb)); } - } - else log_and_throw_error(fmt::format("Vertex {} or {} has no closest component; this should not happen", va, vb)); + } else + log_and_throw_error( + fmt::format( + "Vertex {} or {} has no closest component; this should not happen", + va, + vb)); } } } // Helper: signed distance between two engulfing comps at position p. // Negative = closer to comp_a centroids, positive = closer to comp_b centroids. - const auto voronoi_sign = [&](const Vector2d& p, const size_t comp_a, const size_t comp_b) -> double { + const auto voronoi_sign = + [&](const Vector2d& p, const size_t comp_a, const size_t comp_b) -> double { const auto dists = comp_min_sq_dists(p); - const double da = dists.count(comp_a) ? dists.at(comp_a) : std::numeric_limits::max(); - const double db = dists.count(comp_b) ? dists.at(comp_b) : std::numeric_limits::max(); + const double da = + dists.count(comp_a) ? dists.at(comp_a) : std::numeric_limits::max(); + const double db = + dists.count(comp_b) ? dists.at(comp_b) : std::numeric_limits::max(); return da - db; }; @@ -215,21 +236,46 @@ void ImageSimulationMeshTri::engulf_components( { const double sign_a = voronoi_sign(m_vertex_attribute[va].m_pos, comp_a, comp_b); const double sign_b = voronoi_sign(m_vertex_attribute[vb].m_pos, comp_a, comp_b); - - if (sign_a * sign_b > 0.0) log_and_throw_error(fmt::format("Endpoints {} and {} of edge to split have same sign {}, {}: this should not happen", va, vb, sign_a, sign_b)); - if (sign_a > 0.0 || sign_b < 0.0) log_and_throw_error(fmt::format("Endpoint signs for edge to split are sign_a = {}, sign_b = {}: this should not happen", sign_a, sign_b)); - const double total = sign_b - sign_a; - if (total < 1e-10) // both comps are very close to equidistant — skip split and mark both endpoints as equidistant to avoid numerical issues with voronoi_sign near the boundary + if (sign_a * sign_b > 0.0) + log_and_throw_error( + fmt::format( + "Endpoints {} and {} of edge to split have same sign {}, {}: this " + "should not happen", + va, + vb, + sign_a, + sign_b)); + if (sign_a > 0.0 || sign_b < 0.0) + log_and_throw_error( + fmt::format( + "Endpoint signs for edge to split are sign_a = {}, sign_b = {}: this " + "should not happen", + sign_a, + sign_b)); + + const double total = sign_b - sign_a; + if (total < 1e-10) // both comps are very close to equidistant — skip split and mark + // both endpoints as equidistant to avoid numerical issues with + // voronoi_sign near the boundary { vertex_closest_comp[va] = std::numeric_limits::max(); vertex_closest_comp[vb] = std::numeric_limits::max(); continue; } - constexpr double kMinFrac = 1e-3; // manual tolerance to avoid degenerate splits; if the boundary is within 0.1% of the edge length from an endpoint, skip the split and mark that endpoint as equidistant. + constexpr double kMinFrac = + 1e-3; // manual tolerance to avoid degenerate splits; if the boundary is within + // 0.1% of the edge length from an endpoint, skip the split and mark that + // endpoint as equidistant. const double t = std::abs(sign_a) / total; - if (t < kMinFrac) { vertex_closest_comp[va] = std::numeric_limits::max(); continue; } - if (t > 1.0 - kMinFrac) { vertex_closest_comp[vb] = std::numeric_limits::max(); continue; } + if (t < kMinFrac) { + vertex_closest_comp[va] = std::numeric_limits::max(); + continue; + } + if (t > 1.0 - kMinFrac) { + vertex_closest_comp[vb] = std::numeric_limits::max(); + continue; + } } m_voronoi_split_fn = [&, comp_a, comp_b](const Vector2d& p) -> double { @@ -239,9 +285,12 @@ void ImageSimulationMeshTri::engulf_components( const bool split_success = split_edge(edge_tup, new_tris); m_voronoi_split_fn = nullptr; if (!split_success || new_tris.empty()) { - if (!split_success) logger().warn( - "engulf_components: split_edge({},{}) failed; assignment will use centroid fallback", - va, vb); + if (!split_success) + logger().warn( + "engulf_components: split_edge({},{}) failed; assignment will use centroid " + "fallback", + va, + vb); } else { const size_t vmid = new_tris[0].vid(*this); vertex_closest_comp[vmid] = std::numeric_limits::max(); // equidistant @@ -293,21 +342,29 @@ void ImageSimulationMeshTri::engulf_components( } if (conflict || closest == std::numeric_limits::max()) { const Vector2d centroid = - (m_vertex_attribute[fvs[0]].m_pos + - m_vertex_attribute[fvs[1]].m_pos + - m_vertex_attribute[fvs[2]].m_pos) / 3.0; + (m_vertex_attribute[fvs[0]].m_pos + m_vertex_attribute[fvs[1]].m_pos + + m_vertex_attribute[fvs[2]].m_pos) / + 3.0; closest = compute_vertex_closest_comp(centroid); - if (conflict) logger().warn("Face {} has conflicting vertex assignments after split; using centroid fallback for assignment",fid); - if (closest == std::numeric_limits::max()) logger().warn("Face {} has no definite vertex assignments after split; using centroid fallback for assignment", fid); + if (conflict) + logger().warn( + "Face {} has conflicting vertex assignments after split; using " + "centroid fallback for assignment", + fid); + if (closest == std::numeric_limits::max()) + logger().warn( + "Face {} has no definite vertex assignments after split; using " + "centroid fallback for assignment", + fid); } - + const auto& fp0 = m_vertex_attribute[fvs[0]].m_pos; const auto& fp1 = m_vertex_attribute[fvs[1]].m_pos; const auto& fp2 = m_vertex_attribute[fvs[2]].m_pos; m_face_attribute[fid].tags[0] = components[closest].tag; components[closest].area += 0.5 * std::abs( - (fp1[0] - fp0[0]) * (fp2[1] - fp0[1]) - - (fp1[1] - fp0[1]) * (fp2[0] - fp0[0])); + (fp1[0] - fp0[0]) * (fp2[1] - fp0[1]) - + (fp1[1] - fp0[1]) * (fp2[0] - fp0[0])); components[closest].faces.push_back(fid); comp_id_map[fid] = closest; newly_assigned[closest].push_back(fid); @@ -332,7 +389,7 @@ void ImageSimulationMeshTri::engulf_components( } } - // --- DEBUG: + // --- DEBUG: { static int debug_call = 0; const std::string debug_path = fmt::format("/tmp/engulf_debug_{}.vtu", debug_call++); @@ -355,12 +412,16 @@ void ImageSimulationMeshTri::engulf_components( const int nv = (int)row_vids.size(); Eigen::MatrixXd V(nv, 3); Eigen::VectorXd closest_label(nv), voronoi_sdf(nv); - V.setZero(); closest_label.setZero(); voronoi_sdf.setZero(); + V.setZero(); + closest_label.setZero(); + voronoi_sdf.setZero(); for (int i = 0; i < nv; ++i) { const size_t vid = row_vids[i]; const Vector2d& p = m_vertex_attribute[vid].m_pos; - V(i, 0) = p[0]; V(i, 1) = p[1]; V(i, 2) = 0.0; + V(i, 0) = p[0]; + V(i, 1) = p[1]; + V(i, 2) = 0.0; // Find the two closest comps using the shared helper. const auto comp_min_sq = comp_min_sq_dists(p); @@ -369,20 +430,26 @@ void ImageSimulationMeshTri::engulf_components( double d1sq = std::numeric_limits::max(); double d2sq = std::numeric_limits::max(); for (const auto& [cidx, dsq] : comp_min_sq) { - if (dsq < d1sq) { d2sq = d1sq; c2 = c1; d1sq = dsq; c1 = cidx; } - else if (dsq < d2sq) { d2sq = dsq; c2 = cidx; } + if (dsq < d1sq) { + d2sq = d1sq; + c2 = c1; + d1sq = dsq; + c1 = cidx; + } else if (dsq < d2sq) { + d2sq = dsq; + c2 = cidx; + } } // voronoi_sign(p, c1, c2) = min_sq_dist_c1 - min_sq_dist_c2 // This is exactly what is used in the splitting λ, so values here are // directly comparable to the binary-search target (zero crossing). - voronoi_sdf[i] = (c2 != std::numeric_limits::max()) - ? voronoi_sign(p, c1, c2) - : 0.0; + voronoi_sdf[i] = + (c2 != std::numeric_limits::max()) ? voronoi_sign(p, c1, c2) : 0.0; // Label: closest comp from vertex_closest_comp (-1 = equidistant) const size_t c = vertex_closest_comp.count(vid) - ? vertex_closest_comp.at(vid) - : std::numeric_limits::max(); + ? vertex_closest_comp.at(vid) + : std::numeric_limits::max(); closest_label[i] = (c == std::numeric_limits::max()) ? -1.0 : (double)c; } @@ -392,9 +459,10 @@ void ImageSimulationMeshTri::engulf_components( for (const size_t fid : components[comp_id].faces) { if (!tuple_from_tri(fid).is_valid(*this)) continue; const auto fvs = oriented_tri_vids(fid); - face_list.push_back({(size_t)vid_to_row.at(fvs[0]), - (size_t)vid_to_row.at(fvs[1]), - (size_t)vid_to_row.at(fvs[2])}); + face_list.push_back( + {(size_t)vid_to_row.at(fvs[0]), + (size_t)vid_to_row.at(fvs[1]), + (size_t)vid_to_row.at(fvs[2])}); } } Eigen::MatrixXi F((int)face_list.size(), 3); @@ -403,10 +471,10 @@ void ImageSimulationMeshTri::engulf_components( F(i, 0) = (int)face_list[i][0]; F(i, 1) = (int)face_list[i][1]; F(i, 2) = (int)face_list[i][2]; - const double ax = V(F(i,0),0), ay = V(F(i,0),1); - const double bx = V(F(i,1),0), by = V(F(i,1),1); - const double cx = V(F(i,2),0), cy = V(F(i,2),1); - face_area[i] = 0.5 * std::abs((bx-ax)*(cy-ay) - (by-ay)*(cx-ax)); + const double ax = V(F(i, 0), 0), ay = V(F(i, 0), 1); + const double bx = V(F(i, 1), 0), by = V(F(i, 1), 1); + const double cx = V(F(i, 2), 0), cy = V(F(i, 2), 1); + face_area[i] = 0.5 * std::abs((bx - ax) * (cy - ay) - (by - ay) * (cx - ax)); } // Check for zero/near-zero area faces constexpr double kAreaTol = 1e-10; @@ -417,16 +485,24 @@ void ImageSimulationMeshTri::engulf_components( wmtk::log_and_throw_error( "DEBUG engulf_components: face row {} (verts [{},{},{}]) " "has near-zero area = {} | pos [{},{}] [{},{}] [{},{}]", - i, F(i,0), F(i,1), F(i,2), face_area[i], - V(F(i,0),0), V(F(i,0),1), - V(F(i,1),0), V(F(i,1),1), - V(F(i,2),0), V(F(i,2),1)); + i, + F(i, 0), + F(i, 1), + F(i, 2), + face_area[i], + V(F(i, 0), 0), + V(F(i, 0), 1), + V(F(i, 1), 0), + V(F(i, 1), 1), + V(F(i, 2), 0), + V(F(i, 2), 1)); } } if (zero_area_count > 0) wmtk::log_and_throw_error( "DEBUG engulf_components: {} / {} hole faces have near-zero area", - zero_area_count, face_list.size()); + zero_area_count, + face_list.size()); // Verify: equidistant-labeled vertices should have voronoi_sdf ≈ 0 for SOME // component pair — the pair used when splitting may differ from the globally @@ -437,16 +513,21 @@ void ImageSimulationMeshTri::engulf_components( if (comp_min_sq.size() < 2) continue; // Find the nearest centroid distance (denominator for rel) double d1sq = std::numeric_limits::max(); - for (const auto& [cidx, dsq] : comp_min_sq) - d1sq = std::min(d1sq, dsq); + for (const auto& [cidx, dsq] : comp_min_sq) d1sq = std::min(d1sq, dsq); // Find the pair (ca, cb) that minimizes |voronoi_sign| = |d_ca - d_cb| double min_abs_sign = std::numeric_limits::max(); size_t best_ca = 0, best_cb = 0; - const std::vector> comp_vec(comp_min_sq.begin(), comp_min_sq.end()); + const std::vector> comp_vec( + comp_min_sq.begin(), + comp_min_sq.end()); for (size_t ai = 0; ai < comp_vec.size(); ++ai) for (size_t bi = ai + 1; bi < comp_vec.size(); ++bi) { const double s = std::abs(comp_vec[ai].second - comp_vec[bi].second); - if (s < min_abs_sign) { min_abs_sign = s; best_ca = comp_vec[ai].first; best_cb = comp_vec[bi].first; } + if (s < min_abs_sign) { + min_abs_sign = s; + best_ca = comp_vec[ai].first; + best_cb = comp_vec[bi].first; + } } constexpr double kMinFrac = 0.01; const double rel = (abs(d1sq) > 1e-8) ? min_abs_sign / d1sq : min_abs_sign; @@ -455,7 +536,13 @@ void ImageSimulationMeshTri::engulf_components( wmtk::log_and_throw_error( "DEBUG engulf_components: vertex {} labeled equidistant " "but best pair voronoi_sign(ca={},cb={})={:.4f} rel={:.4f} | pos [{}, {}]", - vid, best_ca, best_cb, min_abs_sign, rel, V(i, 0), V(i, 1)); + vid, + best_ca, + best_cb, + min_abs_sign, + rel, + V(i, 0), + V(i, 1)); } } } @@ -470,7 +557,11 @@ void ImageSimulationMeshTri::engulf_components( } } -void ImageSimulationMeshTri::extract_hole_clusters(std::vector& components, std::unordered_set& tags, std::vector>& hole_clusters, double threshold) +void ImageSimulationMeshTri::extract_hole_clusters( + std::vector& components, + std::unordered_set& tags, + std::vector>& hole_clusters, + double threshold) { // BFS-cluster all non-fill-tag components by adjacency. // Each cluster is a maximal group of connected non-fill-tag components. @@ -478,13 +569,16 @@ void ImageSimulationMeshTri::extract_hole_clusters(std::vector cluster_id(components.size(), -1); std::vector component_clusters; for (size_t i = 0; i < components.size(); ++i) { - if (components[i].faces.empty() || tags.count(components[i].tag) || components[i].touches_boundary) continue; + if (components[i].faces.empty() || tags.count(components[i].tag) || + components[i].touches_boundary) + continue; if (cluster_id[i] != -1) continue; const size_t cid = component_clusters.size(); ComponentCluster& cluster = component_clusters.emplace_back(); cluster.comp_ids.push_back(i); cluster_id[i] = cid; - for (size_t qi = 0; qi < cluster.comp_ids.size(); ++qi) { // BFS loop with cluster.comp_ids growing + for (size_t qi = 0; qi < cluster.comp_ids.size(); + ++qi) { // BFS loop with cluster.comp_ids growing const size_t cur = cluster.comp_ids[qi]; cluster.area += components[cur].area; for (const size_t sid : components[cur].surrounding_comp_ids) { @@ -500,11 +594,11 @@ void ImageSimulationMeshTri::extract_hole_clusters(std::vector Date: Fri, 13 Mar 2026 16:52:37 -0400 Subject: [PATCH 05/10] precommit fix --- .../ImageSimulationMeshTri.cpp | 60 ++++++++++++------- .../ImageSimulationMeshTri.hpp | 25 ++++++-- 2 files changed, 57 insertions(+), 28 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index a353f27c81..701129bf83 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -9,10 +9,10 @@ #include #include #include -#include #include -#include #include +#include +#include #include #include #include @@ -624,19 +624,27 @@ bool ImageSimulationMeshTri::split_edge_after(const Tuple& loc) p = 0.5 * (p0 + p1); bool inv = false; for (const Tuple& t : locs) { - if (is_inverted(t)) { inv = true; break; } + if (is_inverted(t)) { + inv = true; + break; + } } if (inv || (p1 - p0).squaredNorm() < 1e-20) break; - if (m_voronoi_split_fn(p) < 0) p0 = p; - else p1 = p; + if (m_voronoi_split_fn(p) < 0) + p0 = p; + else + p1 = p; } // final inversion guard: revert to midpoint if needed bool inv = false; for (const Tuple& t : locs) { - if (is_inverted(t)) { - inv = true; - logger().warn("Voronoi split resulted in inversion, reverting to midpoint. Iteration: {}", debug_print_counter++); - break; } + if (is_inverted(t)) { + inv = true; + logger().warn( + "Voronoi split resulted in inversion, reverting to midpoint. Iteration: {}", + debug_print_counter++); + break; + } } if (inv) p = (m_vertex_attribute[v1_id].m_pos + m_vertex_attribute[v2_id].m_pos) / 2; } @@ -1555,7 +1563,9 @@ std::tuple ImageSimulationMeshTri::get_max_avg_energy() #include "ConnectedComponentAnnotationHelper.inl" -void ImageSimulationMeshTri::fill_holes_topo(const std::vector& fill_holes_tags, double threshold) +void ImageSimulationMeshTri::fill_holes_topo( + const std::vector& fill_holes_tags, + double threshold) { if (m_tags_count == 0) { logger().warn("fill_holes_topo: no tags, skipping"); @@ -1568,11 +1578,10 @@ void ImageSimulationMeshTri::fill_holes_topo(const std::vector& fill_ho // -- Step 2: for each fill_tag, fill all enclosed components bool any_filled = false; for (const int64_t fill_tag : fill_holes_tags) { - std::vector> hole_clusters; std::unordered_set tags{fill_tag}; extract_hole_clusters(components, tags, hole_clusters, threshold); - + for (const auto& hole_cluster : hole_clusters) { // Collect the fill-tag component indices that surround this cluster std::unordered_set engulfing_comp_ids; @@ -1588,7 +1597,8 @@ void ImageSimulationMeshTri::fill_holes_topo(const std::vector& fill_ho } if (hole_clusters.empty()) { logger().info( - "fill_holes_topo: no enclosed components found for fill_tag {}", fill_tag); + "fill_holes_topo: no enclosed components found for fill_tag {}", + fill_tag); } } @@ -1604,7 +1614,8 @@ void ImageSimulationMeshTri::keep_largest_connected_component(const std::vector< logger().warn("keep_largest_connected_component: no components found, skipping"); return; } - // For each tag in lcc_tags, find the largest component with that tag and mark all other components with that tag for removal + // For each tag in lcc_tags, find the largest component with that tag and mark all other + // components with that tag for removal for (const int64_t lcc_tag : lcc_tags) { int largest_comp_id = -1; double largest_area = -1.0; @@ -1622,7 +1633,8 @@ void ImageSimulationMeshTri::keep_largest_connected_component(const std::vector< continue; } for (size_t i = 0; i < components.size(); ++i) { - if (components[i].faces.empty() || components[i].tag != lcc_tag || (int)i == largest_comp_id) + if (components[i].faces.empty() || components[i].tag != lcc_tag || + (int)i == largest_comp_id) continue; // Copy surrounding_comp_ids — engulf_components will reset components[i] std::unordered_set surr = components[i].surrounding_comp_ids; @@ -1633,7 +1645,9 @@ void ImageSimulationMeshTri::keep_largest_connected_component(const std::vector< recompute_surface_info(); } -void ImageSimulationMeshTri::tight_seal_topo(const std::vector>& tight_seal_tag_sets, double threshold) +void ImageSimulationMeshTri::tight_seal_topo( + const std::vector>& tight_seal_tag_sets, + double threshold) { auto components = compute_connected_components(); @@ -1655,14 +1669,16 @@ void ImageSimulationMeshTri::tight_seal_topo(const std::vector get_max_avg_energy(); - // Returns: (face ids in component, set of distinct tags of all neighboring faces across tag boundaries) - defined according to tag_0 + // Returns: (face ids in component, set of distinct tags of all neighboring faces across tag + // boundaries) - defined according to tag_0 struct ConnectedComponent { int64_t tag = -1; @@ -209,8 +210,14 @@ class ImageSimulationMeshTri : public wmtk::TriMesh bool touches_boundary = false; }; std::vector compute_connected_components() const; - void engulf_component(std::vector& components, const size_t hole_comp_id, const size_t engulfing_comp_id); - void engulf_components(std::vector& components, const std::vector& hole_comp_ids, const std::unordered_set& engulfing_comp_ids); + void engulf_component( + std::vector& components, + const size_t hole_comp_id, + const size_t engulfing_comp_id); + void engulf_components( + std::vector& components, + const std::vector& hole_comp_ids, + const std::unordered_set& engulfing_comp_ids); void keep_largest_connected_component(const std::vector& lcc_tags); // A cluster is a maximal group of connected non-fill-tag ConnectedComponents. @@ -221,14 +228,20 @@ class ImageSimulationMeshTri : public wmtk::TriMesh double area = 0.0; bool enclosed = true; }; - void extract_hole_clusters(std::vector& components, std::unordered_set& tags, std::vector>& hole_clusters, double threshold); + void extract_hole_clusters( + std::vector& components, + std::unordered_set& tags, + std::vector>& hole_clusters, + double threshold); void recompute_surface_info(); void fill_holes_topo( const std::vector& fill_holes_tags, double threshold = std::numeric_limits::infinity()); - - void tight_seal_topo(const std::vector>& tight_seal_tag_sets, double threshold = std::numeric_limits::infinity()); + + void tight_seal_topo( + const std::vector>& tight_seal_tag_sets, + double threshold = std::numeric_limits::infinity()); /** From dc6001a3ca207ef95e7e000ce766a999d68d7382 Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Tue, 17 Mar 2026 11:30:51 +0100 Subject: [PATCH 06/10] PR fixes --- cmake/wmtk_data.cmake | 2 +- ...=> ConnectedComponentAnnotationHelper.cpp} | 197 ++++++++++++------ .../ImageSimulationMeshTri.cpp | 3 +- .../image_simulation/image_simulation.cpp | 12 +- 4 files changed, 141 insertions(+), 73 deletions(-) rename components/image_simulation/wmtk/components/image_simulation/{ConnectedComponentAnnotationHelper.inl => ConnectedComponentAnnotationHelper.cpp} (85%) diff --git a/cmake/wmtk_data.cmake b/cmake/wmtk_data.cmake index 8fd854759c..88a40a38d6 100644 --- a/cmake/wmtk_data.cmake +++ b/cmake/wmtk_data.cmake @@ -16,7 +16,7 @@ ExternalProject_Add( SOURCE_DIR ${WMTK_DATA_ROOT} GIT_REPOSITORY https://github.com/wildmeshing/data2.git - GIT_TAG 6d9e2b2290fccde62abecab0e44f45133d0a2c70 + GIT_TAG 6975028477e56db021b09ed27a6fff644b7956b8 CONFIGURE_COMMAND "" BUILD_COMMAND "" diff --git a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp similarity index 85% rename from components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl rename to components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp index 3db84dabd9..989c34ad44 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.inl +++ b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp @@ -1,5 +1,5 @@ -std::vector -ImageSimulationMeshTri::compute_connected_components() const +std::vector +wmtk::components::image_simulation::tri::ImageSimulationMeshTri::compute_connected_components() const { const size_t n_faces = tri_capacity(); std::vector comp_id(n_faces, -1); @@ -7,13 +7,19 @@ ImageSimulationMeshTri::compute_connected_components() const // Pass 1: BFS to assign component ids and collect faces per component for (size_t fid = 0; fid < n_faces; ++fid) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; - if (comp_id[fid] != -1) continue; + if (!tuple_from_tri(fid).is_valid(*this)) { + continue; + } + if (comp_id[fid] != -1) { + continue; + } const int64_t tag = m_face_attribute[fid].tags[0]; const int comp_idx = (int)components.size(); ConnectedComponent& comp = components.emplace_back(); - if (tag == -1) log_and_throw_error(fmt::format("Face {} has invalid tag -1", fid)); + if (tag == -1) { + log_and_throw_error(fmt::format("Face {} has invalid tag -1", fid)); + } comp.tag = tag; comp.faces.push_back(fid); @@ -31,7 +37,9 @@ ImageSimulationMeshTri::compute_connected_components() const for (int j = 0; j < 3; ++j) { const Tuple edge_tup = tuple_from_edge(cur, j); const auto t_opp = edge_tup.switch_face(*this); - if (!t_opp) continue; + if (!t_opp) { + continue; + } const size_t nbr = t_opp->fid(*this); if (m_face_attribute[nbr].tags[0] == tag && comp_id[nbr] == -1) { comp_id[nbr] = comp_idx; @@ -43,7 +51,9 @@ ImageSimulationMeshTri::compute_connected_components() const // Pass 2: populate surrounding_comp_ids by scanning all boundary edges for (size_t fid = 0; fid < n_faces; ++fid) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; + if (!tuple_from_tri(fid).is_valid(*this)) { + continue; + } const size_t cidx = comp_id[fid]; for (int j = 0; j < 3; ++j) { const Tuple edge_tup = tuple_from_edge(fid, j); @@ -62,8 +72,8 @@ ImageSimulationMeshTri::compute_connected_components() const return components; } -void ImageSimulationMeshTri::engulf_component( - std::vector& components, +void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::engulf_component( + std::vector& components, const size_t hole_comp_id, const size_t engulfing_comp_id) { @@ -90,8 +100,8 @@ void ImageSimulationMeshTri::engulf_component( hole_comp = ConnectedComponent(); } -void ImageSimulationMeshTri::engulf_components( - std::vector& components, +void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::engulf_components( + std::vector& components, const std::vector& hole_comp_ids, const std::unordered_set& engulfing_comp_ids) { @@ -117,11 +127,15 @@ void ImageSimulationMeshTri::engulf_components( std::vector> boundary_edge_centroids; for (const size_t comp_id : hole_comp_ids) { for (const size_t fid : components[comp_id].faces) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; + if (!tuple_from_tri(fid).is_valid(*this)) { + continue; + } for (int j = 0; j < 3; ++j) { const Tuple edge_tup = tuple_from_edge(fid, j); const auto t_opp = edge_tup.switch_face(*this); - if (!t_opp) continue; + if (!t_opp) { + continue; + } const size_t nbr_cidx = comp_id_map[t_opp->fid(*this)]; if (engulfing_comp_ids.count( nbr_cidx)) { // neighbor is one of the engulfing components @@ -140,7 +154,9 @@ void ImageSimulationMeshTri::engulf_components( for (const auto& [centroid, cidx] : boundary_edge_centroids) { const double d = (pos - centroid).squaredNorm(); auto [it, ins] = result.emplace(cidx, d); - if (!ins) it->second = std::min(it->second, d); + if (!ins) { + it->second = std::min(it->second, d); + } } return result; }; @@ -168,7 +184,9 @@ void ImageSimulationMeshTri::engulf_components( for (const size_t fid : components[comp_id].faces) { const auto fvs = oriented_tri_vids(fid); for (const size_t vid : fvs) { - if (vertex_closest_comp.count(vid)) continue; + if (vertex_closest_comp.count(vid)) { + continue; + } vertex_closest_comp[vid] = compute_vertex_closest_comp(m_vertex_attribute[vid].m_pos); } @@ -180,7 +198,9 @@ void ImageSimulationMeshTri::engulf_components( std::set> edges_to_split_set; for (const size_t comp_id : hole_comp_ids) { for (const size_t fid : components[comp_id].faces) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; + if (!tuple_from_tri(fid).is_valid(*this)) { + continue; + } for (int j = 0; j < 3; ++j) { const Tuple edge_tup = tuple_from_edge(fid, j); const auto vs = get_edge_vids(edge_tup); @@ -193,12 +213,13 @@ void ImageSimulationMeshTri::engulf_components( cb != std::numeric_limits::max() && ca != cb) { edges_to_split_set.emplace(std::min(va, vb), std::max(va, vb)); } - } else + } else { log_and_throw_error( fmt::format( "Vertex {} or {} has no closest component; this should not happen", va, vb)); + } } } } @@ -219,14 +240,17 @@ void ImageSimulationMeshTri::engulf_components( // equidistant surface between the two competing engulfing comps via m_voronoi_split_fn. for (const auto& [va, vb] : edges_to_split_set) { auto [edge_tup, _eid] = tuple_from_edge({{va, vb}}); - if (!edge_tup.is_valid(*this)) continue; + if (!edge_tup.is_valid(*this)) { + continue; + } const size_t comp_a = vertex_closest_comp.at(va); const size_t comp_b = vertex_closest_comp.at(vb); // Skip if either endpoint was relabeled equidistant by a previous iteration's // near-endpoint check — voronoi_sign with equidistant comp returns garbage. if (comp_a == std::numeric_limits::max() || - comp_b == std::numeric_limits::max()) + comp_b == std::numeric_limits::max()) { continue; + } // If the Voronoi boundary is so close to va or vb that vmid would land on top of // it, skip the split — it would create a degenerate face. @@ -237,7 +261,7 @@ void ImageSimulationMeshTri::engulf_components( const double sign_a = voronoi_sign(m_vertex_attribute[va].m_pos, comp_a, comp_b); const double sign_b = voronoi_sign(m_vertex_attribute[vb].m_pos, comp_a, comp_b); - if (sign_a * sign_b > 0.0) + if (sign_a * sign_b > 0.0) { log_and_throw_error( fmt::format( "Endpoints {} and {} of edge to split have same sign {}, {}: this " @@ -246,13 +270,15 @@ void ImageSimulationMeshTri::engulf_components( vb, sign_a, sign_b)); - if (sign_a > 0.0 || sign_b < 0.0) + } + if (sign_a > 0.0 || sign_b < 0.0) { log_and_throw_error( fmt::format( "Endpoint signs for edge to split are sign_a = {}, sign_b = {}: this " "should not happen", sign_a, sign_b)); + } const double total = sign_b - sign_a; if (total < 1e-10) // both comps are very close to equidistant — skip split and mark @@ -285,22 +311,26 @@ void ImageSimulationMeshTri::engulf_components( const bool split_success = split_edge(edge_tup, new_tris); m_voronoi_split_fn = nullptr; if (!split_success || new_tris.empty()) { - if (!split_success) - logger().warn( + if (!split_success) { + log_and_throw_error( "engulf_components: split_edge({},{}) failed; assignment will use centroid " "fallback", va, vb); + } } else { const size_t vmid = new_tris[0].vid(*this); vertex_closest_comp[vmid] = std::numeric_limits::max(); // equidistant // Extend comp_id_map and register genuinely new hole faces - if (tri_capacity() > comp_id_map.size()) + if (tri_capacity() > comp_id_map.size()) { comp_id_map.resize(tri_capacity(), std::numeric_limits::max()); + } for (const Tuple& nt : new_tris) { const size_t nfid = nt.fid(*this); - if (comp_id_map[nfid] != std::numeric_limits::max()) continue; + if (comp_id_map[nfid] != std::numeric_limits::max()) { + continue; + } const int64_t ftag = m_face_attribute[nfid].tags[0]; for (const size_t hcid : hole_comp_ids) { if (components[hcid].tag == ftag) { @@ -319,20 +349,25 @@ void ImageSimulationMeshTri::engulf_components( std::unordered_map> newly_assigned; // engulfing cidx -> fids for (const size_t comp_id : hole_comp_ids) { for (const size_t fid : components[comp_id].faces) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; + if (!tuple_from_tri(fid).is_valid(*this)) { + continue; + } const auto fvs = oriented_tri_vids(fid); // Ensure all vertices have an entry (new faces may have new verts) for (const size_t vid : fvs) { - if (!vertex_closest_comp.count(vid)) + if (!vertex_closest_comp.count(vid)) { vertex_closest_comp[vid] = compute_vertex_closest_comp(m_vertex_attribute[vid].m_pos); + } } // Collect definite (non-boundary) vertex assignments size_t closest = std::numeric_limits::max(); bool conflict = false; for (const size_t vid : fvs) { const size_t c = vertex_closest_comp.at(vid); - if (c == std::numeric_limits::max()) continue; // equidistant vertex + if (c == std::numeric_limits::max()) { // equidistant vertex + continue; + } if (closest == std::numeric_limits::max()) { closest = c; } else if (closest != c) { @@ -346,16 +381,18 @@ void ImageSimulationMeshTri::engulf_components( m_vertex_attribute[fvs[2]].m_pos) / 3.0; closest = compute_vertex_closest_comp(centroid); - if (conflict) + if (conflict) { logger().warn( "Face {} has conflicting vertex assignments after split; using " "centroid fallback for assignment", fid); - if (closest == std::numeric_limits::max()) + } + if (closest == std::numeric_limits::max()) { logger().warn( "Face {} has no definite vertex assignments after split; using " "centroid fallback for assignment", fid); + } } const auto& fp0 = m_vertex_attribute[fvs[0]].m_pos; @@ -374,15 +411,25 @@ void ImageSimulationMeshTri::engulf_components( // Update surrounding_comp_ids using newly assigned faces (not stale comp.faces) for (const auto& [cidx, fids] : newly_assigned) { for (const size_t fid : fids) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; + if (!tuple_from_tri(fid).is_valid(*this)) { + continue; + } for (int j = 0; j < 3; ++j) { const Tuple edge_tup = tuple_from_edge(fid, j); const auto t_opp = edge_tup.switch_face(*this); - if (!t_opp) continue; + if (!t_opp) { + continue; + } const size_t nbr_cidx = comp_id_map[t_opp->fid(*this)]; - if (nbr_cidx == std::numeric_limits::max()) continue; - if (nbr_cidx == cidx) continue; - if (hole_comp_ids_set.count(nbr_cidx)) continue; // will be erased + if (nbr_cidx == std::numeric_limits::max()) { + continue; + } + if (nbr_cidx == cidx) { + continue; + } + if (hole_comp_ids_set.count(nbr_cidx)) { // will be erased + continue; + } components[cidx].surrounding_comp_ids.insert(nbr_cidx); components[nbr_cidx].surrounding_comp_ids.insert(cidx); } @@ -399,7 +446,9 @@ void ImageSimulationMeshTri::engulf_components( std::vector row_vids; for (const size_t comp_id : hole_comp_ids) { for (const size_t fid : components[comp_id].faces) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; + if (!tuple_from_tri(fid).is_valid(*this)) { + continue; + } for (const size_t v : oriented_tri_vids(fid)) { if (!vid_to_row.count(v)) { vid_to_row[v] = (int)row_vids.size(); @@ -419,9 +468,7 @@ void ImageSimulationMeshTri::engulf_components( for (int i = 0; i < nv; ++i) { const size_t vid = row_vids[i]; const Vector2d& p = m_vertex_attribute[vid].m_pos; - V(i, 0) = p[0]; - V(i, 1) = p[1]; - V(i, 2) = 0.0; + V.row(i) = Vector3d(p[0], p[1], 0.0); // Find the two closest comps using the shared helper. const auto comp_min_sq = comp_min_sq_dists(p); @@ -457,7 +504,9 @@ void ImageSimulationMeshTri::engulf_components( std::vector> face_list; for (const size_t comp_id : hole_comp_ids) { for (const size_t fid : components[comp_id].faces) { - if (!tuple_from_tri(fid).is_valid(*this)) continue; + if (!tuple_from_tri(fid).is_valid(*this)) { + continue; + } const auto fvs = oriented_tri_vids(fid); face_list.push_back( {(size_t)vid_to_row.at(fvs[0]), @@ -486,41 +535,43 @@ void ImageSimulationMeshTri::engulf_components( "DEBUG engulf_components: face row {} (verts [{},{},{}]) " "has near-zero area = {} | pos [{},{}] [{},{}] [{},{}]", i, - F(i, 0), - F(i, 1), - F(i, 2), + F.row(i).transpose(), face_area[i], - V(F(i, 0), 0), - V(F(i, 0), 1), - V(F(i, 1), 0), - V(F(i, 1), 1), - V(F(i, 2), 0), - V(F(i, 2), 1)); + V.row(F(i, 0)).transpose(), + V.row(F(i, 1)).transpose(), + V.row(F(i, 2)).transpose()); } } - if (zero_area_count > 0) + if (zero_area_count > 0) { wmtk::log_and_throw_error( "DEBUG engulf_components: {} / {} hole faces have near-zero area", zero_area_count, face_list.size()); + } // Verify: equidistant-labeled vertices should have voronoi_sdf ≈ 0 for SOME // component pair — the pair used when splitting may differ from the globally // nearest two, so we check the minimum |voronoi_sign| over all pairs. for (int i = 0; i < nv; ++i) { - if (closest_label[i] != -1.0) continue; + if (closest_label[i] != -1.0) { + continue; + } const auto comp_min_sq = comp_min_sq_dists(m_vertex_attribute[row_vids[i]].m_pos); - if (comp_min_sq.size() < 2) continue; + if (comp_min_sq.size() < 2) { + continue; + } // Find the nearest centroid distance (denominator for rel) double d1sq = std::numeric_limits::max(); - for (const auto& [cidx, dsq] : comp_min_sq) d1sq = std::min(d1sq, dsq); + for (const auto& [cidx, dsq] : comp_min_sq) { + d1sq = std::min(d1sq, dsq); + } // Find the pair (ca, cb) that minimizes |voronoi_sign| = |d_ca - d_cb| double min_abs_sign = std::numeric_limits::max(); size_t best_ca = 0, best_cb = 0; const std::vector> comp_vec( comp_min_sq.begin(), comp_min_sq.end()); - for (size_t ai = 0; ai < comp_vec.size(); ++ai) + for (size_t ai = 0; ai < comp_vec.size(); ++ai) { for (size_t bi = ai + 1; bi < comp_vec.size(); ++bi) { const double s = std::abs(comp_vec[ai].second - comp_vec[bi].second); if (s < min_abs_sign) { @@ -529,6 +580,7 @@ void ImageSimulationMeshTri::engulf_components( best_cb = comp_vec[bi].first; } } + } constexpr double kMinFrac = 0.01; const double rel = (abs(d1sq) > 1e-8) ? min_abs_sign / d1sq : min_abs_sign; if (rel > kMinFrac) { @@ -557,8 +609,8 @@ void ImageSimulationMeshTri::engulf_components( } } -void ImageSimulationMeshTri::extract_hole_clusters( - std::vector& components, +void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::extract_hole_clusters( + std::vector& components, std::unordered_set& tags, std::vector>& hole_clusters, double threshold) @@ -570,9 +622,12 @@ void ImageSimulationMeshTri::extract_hole_clusters( std::vector component_clusters; for (size_t i = 0; i < components.size(); ++i) { if (components[i].faces.empty() || tags.count(components[i].tag) || - components[i].touches_boundary) + components[i].touches_boundary) { + continue; + } + if (cluster_id[i] != -1) { continue; - if (cluster_id[i] != -1) continue; + } const size_t cid = component_clusters.size(); ComponentCluster& cluster = component_clusters.emplace_back(); cluster.comp_ids.push_back(i); @@ -582,8 +637,12 @@ void ImageSimulationMeshTri::extract_hole_clusters( const size_t cur = cluster.comp_ids[qi]; cluster.area += components[cur].area; for (const size_t sid : components[cur].surrounding_comp_ids) { - if (components[sid].faces.empty()) continue; - if (tags.count(components[sid].tag)) continue; + if (components[sid].faces.empty()) { + continue; + } + if (tags.count(components[sid].tag)) { + continue; + } if (components[sid].touches_boundary) { cluster.enclosed = false; // neighbour escapes to mesh boundary continue; @@ -608,8 +667,12 @@ void ImageSimulationMeshTri::extract_hole_clusters( // fill enclosed clusters whose total area is below threshold for (const ComponentCluster& cluster : component_clusters) { - if (!cluster.enclosed) continue; - if (cluster.area >= threshold) continue; + if (!cluster.enclosed) { + continue; + } + if (cluster.area >= threshold) { + continue; + } hole_clusters.push_back({}); for (const size_t ci : cluster.comp_ids) { hole_clusters.back().push_back(ci); @@ -617,7 +680,7 @@ void ImageSimulationMeshTri::extract_hole_clusters( } } -void ImageSimulationMeshTri::recompute_surface_info() +void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::recompute_surface_info() { for (const Tuple& e : get_edges()) { SmartTuple ee(*this, e); @@ -629,7 +692,9 @@ void ImageSimulationMeshTri::recompute_surface_info() for (const Tuple& e : get_edges()) { SmartTuple ee(*this, e); const auto t_opp = ee.switch_face(); - if (!t_opp) continue; + if (!t_opp) { + continue; + } bool has_diff_tag = false; for (size_t j = 0; j < m_tags_count; ++j) { if (m_face_attribute[ee.fid()].tags[j] != m_face_attribute[t_opp->fid()].tags[j]) { @@ -637,7 +702,9 @@ void ImageSimulationMeshTri::recompute_surface_info() break; } } - if (!has_diff_tag) continue; + if (!has_diff_tag) { + continue; + } m_edge_attribute[ee.eid()].m_is_surface_fs = 1; const size_t v1 = ee.vid(); const size_t v2 = ee.switch_vertex().vid(); diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index 701129bf83..c2f2b7eb91 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp @@ -1,4 +1,5 @@ #include "ImageSimulationMeshTri.hpp" +#include "ConnectedComponentAnnotationHelper.cpp" #include #include @@ -1561,8 +1562,6 @@ std::tuple ImageSimulationMeshTri::get_max_avg_energy() return std::make_tuple(max_energy, avg_energy); } -#include "ConnectedComponentAnnotationHelper.inl" - void ImageSimulationMeshTri::fill_holes_topo( const std::vector& fill_holes_tags, double threshold) diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp index 1aa69e0a25..76b710db1a 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp @@ -175,8 +175,8 @@ void run_2D(const nlohmann::json& json_params, const InputData& input_data) mesh.mesh_improvement(max_its); // <-- tetwild } else if (operation == "fill_holes_topo") { const std::vector fill_holes_tags = - json_params["fill_holes_tags"].get>(); - const double raw_threshold = json_params["fill_holes_threshold"].get(); + json_params["fill_holes_tags"]; + const double raw_threshold = json_params["fill_holes_threshold"]; const double threshold = raw_threshold < 0 ? std::numeric_limits::infinity() : raw_threshold; mesh.fill_holes_topo(fill_holes_tags, threshold); @@ -185,16 +185,18 @@ void run_2D(const nlohmann::json& json_params, const InputData& input_data) std::vector> tag_sets; for (const auto& s : json_params["tight_seal_tag_sets"]) { std::unordered_set ts; - for (const auto& v : s) ts.insert(v.get()); + for (const auto& v : s) { + ts.insert(v.get()); + } tag_sets.push_back(std::move(ts)); } - const double raw_threshold = json_params["tight_seal_threshold"].get(); + const double raw_threshold = json_params["tight_seal_threshold"]; const double threshold = raw_threshold < 0 ? std::numeric_limits::infinity() : raw_threshold; mesh.tight_seal_topo(tag_sets, threshold); } else if (operation == "keep_largest_cc") { const std::vector lcc_tags = - json_params["keep_largest_cc_tags"].get>(); + json_params["keep_largest_cc_tags"]; mesh.keep_largest_connected_component(lcc_tags); } else { log_and_throw_error("Unknown image simulation operation"); From 7905d7e7b944433c76ad2c9224051f6ee297350e Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Tue, 17 Mar 2026 14:09:40 +0100 Subject: [PATCH 07/10] PR fixes --- .../ConnectedComponentAnnotationHelper.cpp | 15 +++++++++++---- .../image_simulation/image_simulation.cpp | 6 ++---- 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp index 989c34ad44..5f53d745b2 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp @@ -1,5 +1,6 @@ std::vector -wmtk::components::image_simulation::tri::ImageSimulationMeshTri::compute_connected_components() const +wmtk::components::image_simulation::tri::ImageSimulationMeshTri::compute_connected_components() + const { const size_t n_faces = tri_capacity(); std::vector comp_id(n_faces, -1); @@ -73,7 +74,9 @@ wmtk::components::image_simulation::tri::ImageSimulationMeshTri::compute_connect } void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::engulf_component( - std::vector& components, + std::vector< + wmtk::components::image_simulation::tri::ImageSimulationMeshTri::ConnectedComponent>& + components, const size_t hole_comp_id, const size_t engulfing_comp_id) { @@ -101,7 +104,9 @@ void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::engulf_com } void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::engulf_components( - std::vector& components, + std::vector< + wmtk::components::image_simulation::tri::ImageSimulationMeshTri::ConnectedComponent>& + components, const std::vector& hole_comp_ids, const std::unordered_set& engulfing_comp_ids) { @@ -610,7 +615,9 @@ void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::engulf_com } void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::extract_hole_clusters( - std::vector& components, + std::vector< + wmtk::components::image_simulation::tri::ImageSimulationMeshTri::ConnectedComponent>& + components, std::unordered_set& tags, std::vector>& hole_clusters, double threshold) diff --git a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp index 76b710db1a..e67052cc73 100644 --- a/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp @@ -174,8 +174,7 @@ void run_2D(const nlohmann::json& json_params, const InputData& input_data) if (operation == "remeshing") { mesh.mesh_improvement(max_its); // <-- tetwild } else if (operation == "fill_holes_topo") { - const std::vector fill_holes_tags = - json_params["fill_holes_tags"]; + const std::vector fill_holes_tags = json_params["fill_holes_tags"]; const double raw_threshold = json_params["fill_holes_threshold"]; const double threshold = raw_threshold < 0 ? std::numeric_limits::infinity() : raw_threshold; @@ -195,8 +194,7 @@ void run_2D(const nlohmann::json& json_params, const InputData& input_data) raw_threshold < 0 ? std::numeric_limits::infinity() : raw_threshold; mesh.tight_seal_topo(tag_sets, threshold); } else if (operation == "keep_largest_cc") { - const std::vector lcc_tags = - json_params["keep_largest_cc_tags"]; + const std::vector lcc_tags = json_params["keep_largest_cc_tags"]; mesh.keep_largest_connected_component(lcc_tags); } else { log_and_throw_error("Unknown image simulation operation"); From f0b74821c765b3cdc13e116d79aba7dbd75c2eb9 Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Tue, 17 Mar 2026 15:01:09 +0100 Subject: [PATCH 08/10] add doc, and integration tests --- .../ImageSimulationMeshTri.hpp | 67 +++++++++++++++++-- tests/integration_tests.cpp | 3 + 2 files changed, 66 insertions(+), 4 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp index 587e4153aa..6055ecb79e 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp +++ b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp @@ -199,8 +199,16 @@ class ImageSimulationMeshTri : public wmtk::TriMesh bool collapse_limit_length = true); std::tuple get_max_avg_energy(); - // Returns: (face ids in component, set of distinct tags of all neighboring faces across tag - // boundaries) - defined according to tag_0 + /** + * @brief Connected Component according to tag_0 + * + * @param tag value of tag_0 + * @param faces vector of face ids in the connected component + * @param surrounding_comp_ids set of distinct connected component ids (in an external list) of + * all neighboring faces across tag_0 boundaries + * @param area total area of the connected component + * @param touches_boundary whether the connected component touches the boundary of the mesh + */ struct ConnectedComponent { int64_t tag = -1; @@ -209,25 +217,75 @@ class ImageSimulationMeshTri : public wmtk::TriMesh double area = 0.0; bool touches_boundary = false; }; + + /** + * @brief Compute connected components of the mesh according to tag_0, and return a vector of + * ConnectedComponent structs. + */ std::vector compute_connected_components() const; + + /** + * @brief Replace a component with id hole_comp_id with another component with id + * engulfing_comp_id, and keep the components vector consistent. This is used to fill holes in + * the mesh. + * @param components vector of ConnectedComponent structs, representing the connected components + * of the mesh according to tag_0 + * @param hole_comp_id the id of the component to be engulfed (i.e. the hole) + * @param engulfing_comp_id the id of the component that engulfs the hole + */ void engulf_component( std::vector& components, const size_t hole_comp_id, const size_t engulfing_comp_id); + + /** + * @brief Replace each component in hole_comp_ids with nearest component from + * engulfing_comp_ids, and keep the components vector consistent. + * @param components vector of ConnectedComponent structs, representing the connected components + * of the mesh according to tag_0 + * @param hole_comp_ids vector of ids of components to be engulfed (i.e. the holes) + * @param engulfing_comp_ids set of ids of components that can engulf the holes + */ void engulf_components( std::vector& components, const std::vector& hole_comp_ids, const std::unordered_set& engulfing_comp_ids); + + /** + * @brief Keep only the largest connected component for each of the distinct tag_0 values, and + * engulf all other components. + * + * @param lcc_tags + */ void keep_largest_connected_component(const std::vector& lcc_tags); - // A cluster is a maximal group of connected non-fill-tag ConnectedComponents. - // enclosed = true iff every surrounding component outside the cluster has fill_tag. + /** + * @brief A cluster of connected components. Used to store a maximal group of connected + * non-fill-tag ConnectedComponents for filling holes. + * @param comp_ids vector of indices into the components vector, representing the connected + * components in the cluster + * @param area total area of the cluster + * @param enclosed whether the cluster is enclosed by fill_tag components (i.e. whether the + * cluster is a hole that can be filled) + */ struct ComponentCluster { std::vector comp_ids; // indices into components vector double area = 0.0; bool enclosed = true; }; + + /** + * @brief Extract clusters of connected components that represent holes in the mesh. + * + * @param components vector of ConnectedComponent structs, representing the connected components + * of the mesh according to tag_0 + * @param tags set of tag_0 values which define the holes to be filled + * @param hole_clusters vector of clusters of connected components that represent holes in the + * mesh. + * @param threshold area threshold for the hole clusters. Clusters with area below the threshold + * will not be filled. + */ void extract_hole_clusters( std::vector& components, std::unordered_set& tags, @@ -235,6 +293,7 @@ class ImageSimulationMeshTri : public wmtk::TriMesh double threshold); void recompute_surface_info(); + void fill_holes_topo( const std::vector& fill_holes_tags, double threshold = std::numeric_limits::infinity()); diff --git a/tests/integration_tests.cpp b/tests/integration_tests.cpp index 9201c9ad83..78f47f8ec2 100644 --- a/tests/integration_tests.cpp +++ b/tests/integration_tests.cpp @@ -21,6 +21,9 @@ const path integration_tests_dir = data_dir / "integration_tests"; */ std::vector input_files{ "image_simulation__remeshing_2d.json", + "image_simulation_fill_holes_2d.json", + "image_simulation_keep_lcc_2d.json", + "image_simulation_tight_seal_2d.json", "isotropic_remeshing_bunny.json", "isotropic_remeshing_piece.json", "manifold_extraction_randmesh.json", From 2c259814d18881efdd70acd19bbf878a7fe3558b Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Wed, 18 Mar 2026 17:48:05 +0100 Subject: [PATCH 09/10] namespace --- .../ConnectedComponentAnnotationHelper.cpp | 28 +++++++++---------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp index 5f53d745b2..139e6536dd 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp @@ -1,5 +1,7 @@ -std::vector -wmtk::components::image_simulation::tri::ImageSimulationMeshTri::compute_connected_components() +namespace wmtk::components::image_simulation::tri { + +std::vector +ImageSimulationMeshTri::compute_connected_components() const { const size_t n_faces = tri_capacity(); @@ -73,10 +75,8 @@ wmtk::components::image_simulation::tri::ImageSimulationMeshTri::compute_connect return components; } -void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::engulf_component( - std::vector< - wmtk::components::image_simulation::tri::ImageSimulationMeshTri::ConnectedComponent>& - components, +void ImageSimulationMeshTri::engulf_component( + std::vector& components, const size_t hole_comp_id, const size_t engulfing_comp_id) { @@ -103,10 +103,8 @@ void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::engulf_com hole_comp = ConnectedComponent(); } -void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::engulf_components( - std::vector< - wmtk::components::image_simulation::tri::ImageSimulationMeshTri::ConnectedComponent>& - components, +void ImageSimulationMeshTri::engulf_components( + std::vector& components, const std::vector& hole_comp_ids, const std::unordered_set& engulfing_comp_ids) { @@ -614,10 +612,8 @@ void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::engulf_com } } -void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::extract_hole_clusters( - std::vector< - wmtk::components::image_simulation::tri::ImageSimulationMeshTri::ConnectedComponent>& - components, +void ImageSimulationMeshTri::extract_hole_clusters( + std::vector& components, std::unordered_set& tags, std::vector>& hole_clusters, double threshold) @@ -687,7 +683,7 @@ void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::extract_ho } } -void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::recompute_surface_info() +void ImageSimulationMeshTri::recompute_surface_info() { for (const Tuple& e : get_edges()) { SmartTuple ee(*this, e); @@ -719,3 +715,5 @@ void wmtk::components::image_simulation::tri::ImageSimulationMeshTri::recompute_ m_vertex_attribute[v2].m_is_on_surface = true; } } + +} // namespace wmtk::components::image_simulation::tri From 6ea4a04ecdea60b7cfa0ae360606c837c90bf9ab Mon Sep 17 00:00:00 2001 From: Uday Kusupati Date: Wed, 18 Mar 2026 18:01:17 +0100 Subject: [PATCH 10/10] namespace --- .../image_simulation/ConnectedComponentAnnotationHelper.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp index 139e6536dd..71423f79bc 100644 --- a/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp +++ b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp @@ -1,8 +1,7 @@ namespace wmtk::components::image_simulation::tri { std::vector -ImageSimulationMeshTri::compute_connected_components() - const +ImageSimulationMeshTri::compute_connected_components() const { const size_t n_faces = tri_capacity(); std::vector comp_id(n_faces, -1);