Topological Annotations 2D

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 ""

components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp

@@ -1,4 +1,5 @@
 #include "ImageSimulationMeshTri.hpp"
+#include "ConnectedComponentAnnotationHelper.cpp"
 
 #include <igl/Timer.h>
 #include <igl/is_edge_manifold.h>
@@ -9,7 +10,10 @@
 #include <Eigen/Core>
 #include <Eigen/Geometry>
 #include <atomic>
+#include <map>
 #include <paraviewo/VTUWriter.hpp>
+#include <set>
+#include <unordered_map>
 #include <wmtk/ExecutionScheduler.hpp>
 #include <wmtk/optimization/AMIPSEnergy.hpp>
 #include <wmtk/optimization/DirichletEnergy.hpp>
@@ -607,15 +611,50 @@ 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)) {
             return false;
         }
     }
 
+    // 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
@@ -1583,9 +1622,127 @@ std::tuple<double, double> ImageSimulationMeshTri::get_max_avg_energy()
     return std::make_tuple(max_energy, avg_energy);
 }
 
-void ImageSimulationMeshTri::fill_holes_topo()
+void ImageSimulationMeshTri::fill_holes_topo(
+    const std::vector<int64_t>& 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<std::vector<size_t>> hole_clusters;
+        std::unordered_set<int64_t> 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<size_t> 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<int64_t>& 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<size_t> surr = components[i].surrounding_comp_ids;
+            engulf_components(components, std::vector<size_t>{i}, surr);
+        }
+    }
+
+    recompute_surface_info();
+}
+
+void ImageSimulationMeshTri::tight_seal_topo(
+    const std::vector<std::unordered_set<int64_t>>& tight_seal_tag_sets,
+    double threshold)
 {
-    log_and_throw_error("fill_holes_topo() is not implemented");
+    auto components = compute_connected_components();
+
+    for (const auto& tag_set : tight_seal_tag_sets) {
+        std::vector<std::vector<size_t>> hole_clusters;
+        std::unordered_set<int64_t> 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<size_t> engulfing_comp_ids;
+            std::unordered_set<int64_t> 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 << id << " ";
+            }
+            std::cout << std::endl;
+            std::cout << "Surrounding_tags: ";
+            for (const auto& tag : surrounding_tags) {
+                std::cout << tag << " ";
+            }
+            std::cout << std::endl;
+            engulf_components(components, hole_cluster, engulfing_comp_ids);
+        }
+    }
+
+    recompute_surface_info();
 }
 
 simplex::RawSimplexCollection ImageSimulationMeshTri::get_order1_edges_for_vertex(

components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp

@@ -1,5 +1,8 @@
 #pragma once
 
+#include <limits>
+#include <unordered_set>
+
 #include <wmtk/utils/PartitionMesh.h>
 #include <wmtk/utils/VectorUtils.h>
 #include <wmtk/AttributeCollection.hpp>
@@ -204,7 +207,109 @@ class ImageSimulationMeshTri : public wmtk::TriMesh
         bool collapse_limit_length = true);
     std::tuple<double, double> get_max_avg_energy();
 
-    void fill_holes_topo();
+    /**
+     * @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;
+        std::vector<size_t> faces;
+        std::unordered_set<size_t> surrounding_comp_ids;
+        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<ConnectedComponent> 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<ConnectedComponent>& 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<ConnectedComponent>& components,
+        const std::vector<size_t>& hole_comp_ids,
+        const std::unordered_set<size_t>& 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<int64_t>& lcc_tags);
+
+    /**
+     * @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<size_t> 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<ConnectedComponent>& components,
+        std::unordered_set<int64_t>& tags,
+        std::vector<std::vector<size_t>>& hole_clusters,
+        double threshold);
+
+    void recompute_surface_info();
+
+    void fill_holes_topo(
+        const std::vector<int64_t>& fill_holes_tags,
+        double threshold = std::numeric_limits<double>::infinity());
+
+    void tight_seal_topo(
+        const std::vector<std::unordered_set<int64_t>>& tight_seal_tag_sets,
+        double threshold = std::numeric_limits<double>::infinity());
+
 
     /**
      * @brief Get all edges on the surface that are incident to vid.
@@ -264,6 +369,11 @@ class ImageSimulationMeshTri : public wmtk::TriMesh
         std::vector<int64_t> face_tags;
     };
     tbb::enumerable_thread_specific<SwapInfoCache> 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<double(const Vector2d&)> m_voronoi_split_fn = nullptr;
 };
 
 } // namespace wmtk::components::image_simulation::tri

components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp

@@ -174,7 +174,28 @@ 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<int64_t> 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<double>::infinity() : raw_threshold;
+        mesh.fill_holes_topo(fill_holes_tags, threshold);
+    } else if (operation == "tight_seal_topo") {
+        // tight_seal_tag_sets is a list of lists: [[t1,t2],[t3,...]]
+        std::vector<std::unordered_set<int64_t>> tag_sets;
+        for (const auto& s : json_params["tight_seal_tag_sets"]) {
+            std::unordered_set<int64_t> ts;
+            for (const auto& v : s) {
+                ts.insert(v.get<int64_t>());
+            }
+            tag_sets.push_back(std::move(ts));
+        }
+        const double raw_threshold = json_params["tight_seal_threshold"];
+        const double threshold =
+            raw_threshold < 0 ? std::numeric_limits<double>::infinity() : raw_threshold;
+        mesh.tight_seal_topo(tag_sets, threshold);
+    } else if (operation == "keep_largest_cc") {
+        const std::vector<int64_t> lcc_tags = json_params["keep_largest_cc_tags"];
+        mesh.keep_largest_connected_component(lcc_tags);
     } else {
         log_and_throw_error("Unknown image simulation operation");
     }
@@ -242,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;