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.cpp b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp new file mode 100644 index 0000000000..71423f79bc --- /dev/null +++ b/components/image_simulation/wmtk/components/image_simulation/ConnectedComponentAnnotationHelper.cpp @@ -0,0 +1,718 @@ +namespace wmtk::components::image_simulation::tri { + +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: 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 { + const auto dists = comp_min_sq_dists(pos); + 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; + } + } + 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 + 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)); + } + } + } + } + + // 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 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 + // 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); + + 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; + } + } + + 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) { + 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()) { + 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()) { // equidistant vertex + continue; + } + 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)) { // will be erased + continue; + } + components[cidx].surrounding_comp_ids.insert(nbr_cidx); + components[nbr_cidx].surrounding_comp_ids.insert(cidx); + } + } + } + + // --- DEBUG: + { + 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.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); + 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::log_and_throw_error( + "DEBUG engulf_components: face row {} (verts [{},{},{}]) " + "has near-zero area = {} | pos [{},{}] [{},{}] [{},{}]", + i, + F.row(i).transpose(), + face_area[i], + V.row(F(i, 0)).transpose(), + V.row(F(i, 1)).transpose(), + V.row(F(i, 2)).transpose()); + } + } + 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; + } + 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(); + 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 = (abs(d1sq) > 1e-8) ? min_abs_sign / d1sq : min_abs_sign; + if (rel > kMinFrac) { + const size_t vid = row_vids[i]; + 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)); + } + } + } + // --- 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; + } +} + +} // namespace wmtk::components::image_simulation::tri diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.cpp index d94cb03ec9..a592d9c92d 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 @@ -9,7 +10,10 @@ #include #include #include +#include #include +#include +#include #include #include #include @@ -607,8 +611,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)) { @@ -616,6 +620,41 @@ 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 @@ -1583,9 +1622,127 @@ std::tuple 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& 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) { - 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> 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; + 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 << 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( diff --git a/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp b/components/image_simulation/wmtk/components/image_simulation/ImageSimulationMeshTri.hpp index 31923b05a7..0ce7c04d6c 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 @@ -204,7 +207,109 @@ class ImageSimulationMeshTri : public wmtk::TriMesh bool collapse_limit_length = true); std::tuple 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 faces; + std::unordered_set 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 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); + + /** + * @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, + 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. @@ -264,6 +369,11 @@ class ImageSimulationMeshTri : public wmtk::TriMesh std::vector 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 b12db9d87f..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,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 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; + 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> 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()); + } + 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::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"]; + 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; 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..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 @@ -26,7 +26,12 @@ 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", + "keep_largest_cc_tags" ] }, { @@ -50,7 +55,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", "keep_largest_cc"], "doc": "Image simulation contains multiple operations for modifying the image. Depending on the operation, more parameters might be required." }, { @@ -171,6 +176,56 @@ 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 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 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 b2847521f9..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 @@ -21,7 +21,12 @@ "log_file", "report", "DEBUG_output", - "DEBUG_sanity_checks" + "DEBUG_sanity_checks", + "fill_holes_tags", + "fill_holes_threshold", + "tight_seal_tag_sets", + "tight_seal_threshold", + "keep_largest_cc_tags" ] }, { @@ -45,7 +50,7 @@ "pointer": "/operation", "type": "string", "default": "remeshing", - "options": ["remeshing", "fill_holes_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." }, { @@ -166,5 +171,55 @@ "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 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 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." } ] diff --git a/tests/integration_tests.cpp b/tests/integration_tests.cpp index a6fe21369c..054a9c2fc9 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",