diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index c69e5890..3e99a20e 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -27,15 +27,19 @@ // #include #include #include +#include #include #include +#include #include #include #include #include #include +#include #include +#include namespace roofer::reconstruction { @@ -48,8 +52,14 @@ namespace roofer::reconstruction { typedef CGAL::Triangulation_vertex_base_2 VertexBase; typedef CGAL::Triangulation_vertex_base_with_info_2 VertexBaseWithInfo; + struct TriFaceInfo { + FaceInfo* label = nullptr; + }; typedef CGAL::Constrained_triangulation_face_base_2 FaceBase; - typedef CGAL::Triangulation_data_structure_2 + typedef CGAL::Triangulation_face_base_with_info_2 + FaceBaseWithInfo; + typedef CGAL::Triangulation_data_structure_2 TriangulationDataStructure; typedef CGAL::Constrained_Delaunay_triangulation_2< K, TriangulationDataStructure, Tag> @@ -62,7 +72,26 @@ namespace roofer::reconstruction { typedef T::Face_handle Face_handle; typedef std::pair Edge; - typedef std::unordered_set ConstraintSet; + struct ConstraintToRestore { + Vertex_handle other; + FaceInfo* left_label = nullptr; + FaceInfo* right_label = nullptr; + }; + + struct FaceLabelSeed { + Face_handle face; + FaceInfo* label = nullptr; + }; + + struct LabelledConstraint { + Vertex_handle first; + Vertex_handle second; + FaceInfo* left_label = nullptr; + FaceInfo* right_label = nullptr; + }; + + typedef std::unordered_set VertexSet; + typedef std::vector ConstraintsToRestore; struct ForcedRegionLabel { T::Point_2 seed; @@ -148,14 +177,7 @@ namespace roofer::reconstruction { using ForcedFaceLabels = std::vector>; using TriangleFaceLabels = std::unordered_map; - struct LabelledRegion { - T::Point_2 sample; - FaceInfo* label; - double area; - }; - struct RegionLabelling { - std::vector regions; TriangleFaceLabels face_labels; }; @@ -171,11 +193,28 @@ namespace roofer::reconstruction { return forced_face_labels; } - FaceInfo* select_region_label( + FaceInfo* select_label_by_weight( + const std::unordered_map& votes, + const SourceFaceIds& source_ids) { + FaceInfo* selected = nullptr; + double selected_count = -1; + std::size_t selected_id = std::numeric_limits::max(); + for (const auto& [source, count] : votes) { + auto id = source_ids.at(source); + if (count > selected_count || + (count == selected_count && id < selected_id)) { + selected = source; + selected_count = count; + selected_id = id; + } + } + return selected; + } + + FaceInfo* select_region_label_from_source( T& tri, const std::vector& region, CGAL::Arr_walk_along_line_point_location& walk_pl, - Arrangement_2& source_arrangement, const SourceFaceIds& source_ids, - const ForcedFaceLabels& forced_face_labels) { + Arrangement_2& source_arrangement, const SourceFaceIds& source_ids) { // For this triangle region collect total overlap area for each source // arrangement face. std::unordered_map overlap_area; @@ -188,61 +227,112 @@ namespace roofer::reconstruction { } // Select the arrangement face with the largest overlap area. - FaceInfo* selected = nullptr; - double selected_area = -1; - std::size_t selected_id = std::numeric_limits::max(); - for (const auto& [source, area] : overlap_area) { - auto id = source_ids.at(source); - if (area > selected_area || - (area == selected_area && id < selected_id)) { - selected = source; - selected_area = area; - selected_id = id; + return select_label_by_weight(overlap_area, source_ids); + } + + void seed_region_labels_from_source( + T& tri, + CGAL::Arr_walk_along_line_point_location& walk_pl, + Arrangement_2& source_arrangement, const SourceFaceIds& source_ids) { + auto regions = constrained_regions(tri); + for (const auto& region : regions) { + auto* label = select_region_label_from_source( + tri, region, walk_pl, source_arrangement, source_ids); + if (label == nullptr) continue; + + for (auto face : region) { + face->info().label = label; + } + } + for (auto face : tri.all_face_handles()) { + if (tri.is_infinite(face)) { + face->info().label = &source_arrangement.unbounded_face()->data(); } } + } - // Apply the forced regions determined by non-manifold repairs. + FaceInfo* forced_label_for_region( + const std::vector& region, + const ForcedFaceLabels& forced_face_labels) { for (const auto& [forced_face, forced_label] : forced_face_labels) { if (std::find(region.begin(), region.end(), forced_face) != region.end()) { - selected = forced_label; - break; + return forced_label; } } - return selected; + return nullptr; + } + + FaceInfo* seeded_label_for_region( + const std::vector& region, + const std::vector& label_seeds, + const SourceFaceIds& source_ids) { + std::unordered_map seed_votes; + for (const auto& seed : label_seeds) { + if (seed.label == nullptr) continue; + if (std::find(region.begin(), region.end(), seed.face) != + region.end()) { + seed_votes[seed.label] += 1; + } + } + return select_label_by_weight(seed_votes, source_ids); + } + + FaceInfo* existing_label_for_region(T& tri, + const std::vector& region, + const SourceFaceIds& source_ids) { + std::unordered_map overlap_area; + for (auto face : region) { + auto* label = face->info().label; + if (label == nullptr) continue; + overlap_area[label] += std::abs(tri.triangle(face).area()); + } + return select_label_by_weight(overlap_area, source_ids); } RegionLabelling compute_region_labelling( - T& tri, - CGAL::Arr_walk_along_line_point_location& walk_pl, - Arrangement_2& source_arrangement, const SourceFaceIds& source_ids, - const std::vector& forced_labels) { + T& tri, const SourceFaceIds& source_ids, + const std::vector& forced_labels, + FaceInfo* unbounded_label, + const std::vector& label_seeds = {}) { RegionLabelling labelling; auto forced_face_labels = locate_forced_labels(tri, forced_labels); + // Full relabelling pass: label_seeds only cover newly inserted + // constraints; regions without seeds still need forced/existing-label + // resolution, especially after constraint removal or region merges. + // TODO: only do full relabelling if absolutely necessary (ie. no forced + // labels) auto regions = constrained_regions(tri); for (const auto& region : regions) { - auto* label = - select_region_label(tri, region, walk_pl, source_arrangement, - source_ids, forced_face_labels); + auto* label = forced_label_for_region(region, forced_face_labels); + if (label == nullptr) { + label = seeded_label_for_region(region, label_seeds, source_ids); + } + if (label == nullptr) { + label = existing_label_for_region(tri, region, source_ids); + } if (label == nullptr) continue; - double area = 0; for (auto face : region) { + face->info().label = label; labelling.face_labels[face] = label; - area += std::abs(tri.triangle(face).area()); } - labelling.regions.push_back( - {CGAL::centroid(tri.triangle(region.front())), label, area}); } for (auto face : tri.all_face_handles()) { if (tri.is_infinite(face)) { - labelling.face_labels[face] = - &source_arrangement.unbounded_face()->data(); + face->info().label = unbounded_label; + labelling.face_labels[face] = unbounded_label; } } return labelling; } + void relabel_from_existing_region_labels(T& tri, + const SourceFaceIds& source_ids, + FaceInfo* unbounded_label) { + compute_region_labelling(tri, source_ids, {}, unbounded_label); + } + double squared_distance_to_face(T& tri, Face_handle face, const T::Point_2& point) { auto triangle = tri.triangle(face); @@ -349,11 +439,31 @@ namespace roofer::reconstruction { return count; } + FaceInfo* incident_region_label(T& tri, Vertex_handle vertex) { + Face_circulator face = tri.incident_faces(vertex), done(face); + if (face == nullptr) return nullptr; + + do { + if (!tri.is_infinite(face) && face->info().label != nullptr) { + return face->info().label; + } + } while (++face != done); + return nullptr; + } + + void seed_removed_vertex_region(T& tri, const T::Point_2& point, + FaceInfo* label) { + if (label == nullptr) return; + + auto face = tri.locate(point); + if (!tri.is_infinite(face)) face->info().label = label; + } + void remove_dangling_constraints_and_vertices(T& tri) { // Build the constraint graph and peel its degree-one vertices. Updating // the graph as edges are peeled also detects constraint chains that only // become dangling after their outermost edge has been removed. - std::unordered_map adjacency; + std::unordered_map adjacency; for (const auto& edge : tri.constrained_edges()) { auto first = edge.first->vertex(tri.cw(edge.second)); auto second = edge.first->vertex(tri.ccw(edge.second)); @@ -397,7 +507,12 @@ namespace roofer::reconstruction { vertices_to_remove.push_back(vertex); } } - for (auto vertex : vertices_to_remove) tri.remove(vertex); + for (auto vertex : vertices_to_remove) { + auto* label = incident_region_label(tri, vertex); + auto point = vertex->point(); + tri.remove(vertex); + seed_removed_vertex_region(tri, point, label); + } } // Calculate azimuth for each incident triangulation edge and sample the @@ -673,6 +788,18 @@ namespace roofer::reconstruction { tri.insert_constraint(first, second); } + LabelledConstraint labelled_constraint_from_side( + Vertex_handle first, Vertex_handle second, const T::Point_2& side_point, + FaceInfo* side_label, FaceInfo* other_label); + + void insert_labelled_constraints( + T& tri, const std::vector& constraints, + std::vector& label_seeds); + + void insert_labelled_constraints_and_relabel( + T& tri, const std::vector& constraints, + const SourceFaceIds& source_ids, FaceInfo* unbounded_label); + std::optional find_problematic_vertex( T& tri, const TriangleFaceLabels& face_labels, const SourceFaceIds& source_ids, double repair_radius, @@ -708,6 +835,7 @@ namespace roofer::reconstruction { ForcedRegionLabel repair_vertex(T& tri, const RepairCandidate& candidate, const SourceFaceIds& source_ids, + FaceInfo* unbounded_label, double maximum_radius, double height_tolerance) { const auto& sectors = candidate.sectors; @@ -773,6 +901,8 @@ namespace roofer::reconstruction { std::vector split_vertices; split_vertices.reserve(sectors.size()); + std::vector labelled_constraints; + labelled_constraints.reserve(sectors.size() * 2); const auto split_radii = repair_cell_radii(repair_rays, radius); for (std::size_t i = 0; i < sectors.size(); ++i) { T::Point_2 split_point( @@ -781,15 +911,16 @@ namespace roofer::reconstruction { auto split_vertex = tri.insert(split_point); split_vertex->info() = false; split_vertices.push_back(split_vertex); + + auto previous = (i + sectors.size() - 1) % sectors.size(); + labelled_constraints.push_back({split_vertex, sectors[i].neighbour, + sectors[i].face_info, + sectors[previous].face_info}); } // remove incident constraints and out main vertex tri.remove_incident_constraints(vertex); tri.remove(vertex); - for (std::size_t i = 0; i < sectors.size(); ++i) { - insert_constraint_if_distinct(tri, split_vertices[i], - sectors[i].neighbour); - } // std::optional forced_seed; @@ -809,10 +940,14 @@ namespace roofer::reconstruction { } continue; } - insert_constraint_if_distinct(tri, split_vertices[i], - split_vertices[next]); + labelled_constraints.push_back(labelled_constraint_from_side( + split_vertices[i], split_vertices[next], original_point, + selected_face, sectors[i].face_info)); } + insert_labelled_constraints_and_relabel(tri, labelled_constraints, + source_ids, unbounded_label); + if (!forced_seed) { throw roofer::rooferException( "Non-manifold junction repair has no merge sector"); @@ -820,10 +955,137 @@ namespace roofer::reconstruction { return {*forced_seed, selected_face}; } - void get_incident_constraints(T& tri, Vertex_handle vthis, - Vertex_handle vexcept1, - Vertex_handle vexcept2, - ConstraintSet& constraints_to_restore) { + std::optional face_is_left_of_segment(T& tri, Face_handle face, + int edge_index, + Vertex_handle from, + Vertex_handle to) { + if (tri.is_infinite(face)) return std::nullopt; + + auto opposite = face->vertex(edge_index); + auto orientation = + CGAL::orientation(from->point(), to->point(), opposite->point()); + if (orientation == CGAL::LEFT_TURN) return true; + if (orientation == CGAL::RIGHT_TURN) return false; + return std::nullopt; + } + + void set_side_label(ConstraintToRestore& constraint, bool is_left, + FaceInfo* label) { + if (label == nullptr) return; + if (is_left) { + constraint.left_label = label; + } else { + constraint.right_label = label; + } + } + + ConstraintToRestore get_labelled_constraint(T& tri, const Edge& edge, + Vertex_handle from, + Vertex_handle to) { + ConstraintToRestore constraint{to, nullptr, nullptr}; + + auto face = edge.first; + auto neighbour = face->neighbor(edge.second); + auto face_is_left_side = + face_is_left_of_segment(tri, face, edge.second, from, to); + auto neighbour_index = neighbour->index(face); + auto neighbour_is_left_side = + face_is_left_of_segment(tri, neighbour, neighbour_index, from, to); + + if (face_is_left_side) { + set_side_label(constraint, *face_is_left_side, face->info().label); + set_side_label(constraint, !*face_is_left_side, + neighbour->info().label); + } + + if (neighbour_is_left_side) { + set_side_label(constraint, *neighbour_is_left_side, + neighbour->info().label); + set_side_label(constraint, !*neighbour_is_left_side, + face->info().label); + } + + return constraint; + } + + void add_constraint_face_label_seed(T& tri, Face_handle face, + int edge_index, Vertex_handle from, + Vertex_handle to, FaceInfo* left_label, + FaceInfo* right_label, + std::vector& seeds) { + auto face_side = face_is_left_of_segment(tri, face, edge_index, from, to); + if (!face_side) return; + + auto* label = *face_side ? left_label : right_label; + if (label == nullptr) return; + + face->info().label = label; + seeds.push_back({face, label}); + } + + // TODO: do we really need to insert constraints twice? just one side should + // be enough since the incident edges are always circular (and this edge's + // left label is the neighbour's right label) + void add_constraint_label_seeds(T& tri, const Edge& edge, + Vertex_handle from, Vertex_handle to, + FaceInfo* left_label, FaceInfo* right_label, + std::vector& seeds) { + add_constraint_face_label_seed(tri, edge.first, edge.second, from, to, + left_label, right_label, seeds); + + auto neighbour = edge.first->neighbor(edge.second); + auto neighbour_index = neighbour->index(edge.first); + add_constraint_face_label_seed(tri, neighbour, neighbour_index, from, to, + left_label, right_label, seeds); + } + + LabelledConstraint labelled_constraint_from_side( + Vertex_handle first, Vertex_handle second, const T::Point_2& side_point, + FaceInfo* side_label, FaceInfo* other_label) { + auto orientation = + CGAL::orientation(first->point(), second->point(), side_point); + if (orientation == CGAL::RIGHT_TURN) { + return {first, second, other_label, side_label}; + } + return {first, second, side_label, other_label}; + } + + void insert_labelled_constraints( + T& tri, const std::vector& constraints, + std::vector& label_seeds) { + for (const auto& constraint : constraints) { + if (constraint.first == constraint.second) continue; + insert_constraint_if_distinct(tri, constraint.first, constraint.second); + } + + for (const auto& constraint : constraints) { + if (constraint.first == constraint.second) continue; + + Face_handle face; + int index; + if (tri.is_edge(constraint.first, constraint.second, face, index) && + tri.is_constrained({face, index})) { + add_constraint_label_seeds(tri, {face, index}, constraint.first, + constraint.second, constraint.left_label, + constraint.right_label, label_seeds); + } + } + } + + void insert_labelled_constraints_and_relabel( + T& tri, const std::vector& constraints, + const SourceFaceIds& source_ids, FaceInfo* unbounded_label) { + std::vector label_seeds; + insert_labelled_constraints(tri, constraints, label_seeds); + compute_region_labelling(tri, source_ids, {}, unbounded_label, + label_seeds); + } + + // collect the incident constraints on this vertex, excluding the ones that + // contain vexcept1 or vexcept2 + void get_incident_constraints( + T& tri, Vertex_handle vthis, Vertex_handle vexcept1, + Vertex_handle vexcept2, ConstraintsToRestore& constraints_to_restore) { // std::cout << "vthis degree=" << tri.degree(vthis) << std::endl; Edge_circulator ec = tri.incident_edges(vthis), done(ec); if (ec != nullptr) { @@ -839,26 +1101,33 @@ namespace roofer::reconstruction { }; if (vother != vexcept1 && vother != vexcept2) { - constraints_to_restore.insert(vother); + constraints_to_restore.push_back( + get_labelled_constraint(tri, *ec, vthis, vother)); } } } while (++ec != done); } } - void restore_constraints(T& tri, T::Point_2& pnew, bool boundary_vertex, - ConstraintSet& constraints_to_restore) { + void restore_constraints(T& tri, const T::Point_2& pnew, + bool is_boundary_vertex, + const ConstraintsToRestore& constraints_to_restore, + const SourceFaceIds& source_ids, + FaceInfo* unbounded_label) { auto vnew = tri.insert(pnew); - vnew->info() = boundary_vertex; - - // restore constraints - // std::cout << "restoring " << constraints_to_restore.size() << " - // constraints\n"; - for (auto& vh : constraints_to_restore) { - // std::cout << "reinsert constrained " << *vnew << " - " << *vh << - // std::endl; - tri.insert_constraint(vnew, vh); + vnew->info() = is_boundary_vertex; + std::vector labelled_constraints; + labelled_constraints.reserve(constraints_to_restore.size()); + + for (const auto& constraint : constraints_to_restore) { + if (vnew == constraint.other) continue; + + labelled_constraints.push_back({vnew, constraint.other, + constraint.left_label, + constraint.right_label}); } + insert_labelled_constraints_and_relabel(tri, labelled_constraints, + source_ids, unbounded_label); } class ArrangementSnapper : public ArrangementSnapperInterface { @@ -875,16 +1144,13 @@ namespace roofer::reconstruction { for (auto face : arr.face_handles()) { source_face_ids[&face->data()] = next_source_face_id++; } + auto* unbounded_label = &arr.unbounded_face()->data(); // map from arr vertices to tri vertices std::unordered_map vertex_map; - // Segment_list_2 seg_list; - // Polyline_list_2 output_list; for (auto arrVertex : arr.vertex_handles()) { - // auto& p = v->point(); - // check if this vertex is on the footprint Arrangement_2::Halfedge_around_vertex_circulator ec = arrVertex->incident_halfedges(), @@ -911,6 +1177,10 @@ namespace roofer::reconstruction { vertex_map[arrEdge->target()]); } } + // TODOX: use more of a combinatorics approach rather than the area + // voting. Ie. take one edge per arr face (and store face), link to + // matching constraint in cdt, and kick off face labeling from there. + seed_region_labels_from_source(tri, walk_pl, arr, source_face_ids); // Detect triangles with 3 short edges => collapse triangle to point // (remove 2 vertices) @@ -943,7 +1213,7 @@ namespace roofer::reconstruction { ) { // std::cout << "small triangle between " << p0 << " and " << p1 // << " and " << p2 << std::endl; - ConstraintSet constraints_to_restore; + ConstraintsToRestore constraints_to_restore; // collect incident constraint edges for (size_t i = 0; i < 3; ++i) { @@ -978,7 +1248,8 @@ namespace roofer::reconstruction { tri.remove(v2); restore_constraints(tri, pnew, boundary_vertex, - constraints_to_restore); + constraints_to_restore, source_face_ids, + unbounded_label); found_small_face = true; break; // we need to restart the loop because we may have @@ -1009,7 +1280,7 @@ namespace roofer::reconstruction { // << std::endl; // auto vi = ceit->second; - ConstraintSet constraints_to_restore; + ConstraintsToRestore constraints_to_restore; get_incident_constraints(tri, v1, v1, v2, constraints_to_restore); get_incident_constraints(tri, v2, v1, v2, constraints_to_restore); @@ -1034,7 +1305,8 @@ namespace roofer::reconstruction { tri.remove(v2); restore_constraints(tri, pnew, boundary_vertex, - constraints_to_restore); + constraints_to_restore, source_face_ids, + unbounded_label); found_short_edge = true; break; @@ -1068,8 +1340,12 @@ namespace roofer::reconstruction { // std::cout << "flat triangle between " << s0 << " and " << p0 << // std::endl; tri.remove_constrained_edge(fit, 0); - if (!tri.is_constrained(e2)) tri.insert_constraint(v0, v1); - if (!tri.is_constrained(e1)) tri.insert_constraint(v0, v2); + if (!tri.is_constrained(e2)) + insert_constraint_if_distinct(tri, v0, v1); + if (!tri.is_constrained(e1)) + insert_constraint_if_distinct(tri, v0, v2); + relabel_from_existing_region_labels(tri, source_face_ids, + unbounded_label); } } else if ((CGAL::squared_distance(s1, p1) < sq_dist_thres) && tri.is_constrained(e1)) { @@ -1077,8 +1353,12 @@ namespace roofer::reconstruction { // std::cout << "flat triangle between " << s1 << " and " << p1 << // std::endl; tri.remove_constrained_edge(fit, 1); - if (!tri.is_constrained(e2)) tri.insert_constraint(v1, v0); - if (!tri.is_constrained(e0)) tri.insert_constraint(v1, v2); + if (!tri.is_constrained(e2)) + insert_constraint_if_distinct(tri, v1, v0); + if (!tri.is_constrained(e0)) + insert_constraint_if_distinct(tri, v1, v2); + relabel_from_existing_region_labels(tri, source_face_ids, + unbounded_label); } } else if ((CGAL::squared_distance(s2, p2) < sq_dist_thres) && tri.is_constrained(e2)) { @@ -1086,8 +1366,12 @@ namespace roofer::reconstruction { // std::cout << "flat triangle between " << s2 << " and " << p2 << // std::endl; tri.remove_constrained_edge(fit, 2); - if (!tri.is_constrained(e1)) tri.insert_constraint(v2, v0); - if (!tri.is_constrained(e0)) tri.insert_constraint(v2, v1); + if (!tri.is_constrained(e1)) + insert_constraint_if_distinct(tri, v2, v0); + if (!tri.is_constrained(e0)) + insert_constraint_if_distinct(tri, v2, v1); + relabel_from_existing_region_labels(tri, source_face_ids, + unbounded_label); } } } @@ -1104,7 +1388,7 @@ namespace roofer::reconstruction { RegionLabelling intermediate_labelling; while (true) { intermediate_labelling = compute_region_labelling( - tri, walk_pl, arr, source_face_ids, forced_region_labels); + tri, source_face_ids, forced_region_labels, unbounded_label); if (!cfg.repair_non_manifold_vertices) break; auto candidate = find_problematic_vertex( @@ -1113,8 +1397,8 @@ namespace roofer::reconstruction { exterior_height_provider); if (!candidate) break; forced_region_labels.push_back(repair_vertex( - tri, *candidate, source_face_ids, cfg.manifold_repair_radius, - cfg.manifold_height_tolerance)); + tri, *candidate, source_face_ids, unbounded_label, + cfg.manifold_repair_radius, cfg.manifold_height_tolerance)); } // convert back from triangulation to arrangement @@ -1140,6 +1424,7 @@ namespace roofer::reconstruction { auto& p2_ = v2->point(); auto p1 = Arrangement_2::Point_2(p1_.x(), p1_.y()); auto p2 = Arrangement_2::Point_2(p2_.x(), p2_.y()); + if (vertex2arr_map[v1] == vertex2arr_map[v2]) continue; // std::cout << p1 << " -- " << p2 << std::endl; @@ -1158,19 +1443,77 @@ namespace roofer::reconstruction { std::unordered_map> output_face_labels; - for (const auto& region : intermediate_labelling.regions) { - if (!(region.area > 0)) continue; + for (auto face : tri.finite_face_handles()) { + auto* label = face->info().label; + if (label == nullptr) continue; + + const auto triangle = tri.triangle(face); + const auto area = std::abs(triangle.area()); + if (!(area > 0)) continue; + const auto sample = CGAL::centroid(triangle); auto object = snap_walk_pl.locate( - Arrangement_2::Point_2(region.sample.x(), region.sample.y())); + Arrangement_2::Point_2(sample.x(), sample.y())); if (auto located_face = std::get_if(&object)) { auto output_face = arr_snap.non_const_handle(*located_face); if (output_face->is_unbounded()) continue; - output_face_labels[output_face][region.label] += region.area; + output_face_labels[output_face][label] += area; } } arr_snap.unbounded_face()->data() = arr.unbounded_face()->data(); + + auto add_adjacent_label_votes = + [&](Arrangement_2::Face_handle output_face) { + auto& labels = output_face_labels[output_face]; + auto add_boundary_votes = [&](auto ccb) { + auto current = ccb; + do { + auto adjacent = current->twin()->face(); + if (adjacent != output_face) { + if (adjacent->is_unbounded()) { + labels[unbounded_label] += 1; + } else if (auto adjacent_labels = + output_face_labels.find(adjacent); + adjacent_labels != output_face_labels.end()) { + for (const auto& [label, weight] : + adjacent_labels->second) { + labels[label] += weight; + } + } + } + } while (++current != ccb); + }; + + for (auto ccb = output_face->outer_ccbs_begin(); + ccb != output_face->outer_ccbs_end(); ++ccb) { + add_boundary_votes(*ccb); + } + for (auto hole = output_face->holes_begin(); + hole != output_face->holes_end(); ++hole) { + add_boundary_votes(*hole); + } + if (labels.empty()) output_face_labels.erase(output_face); + }; + + bool propagated_label; + do { + propagated_label = false; + for (auto output_face : arr_snap.face_handles()) { + if (output_face->is_unbounded()) continue; + if (auto labels = output_face_labels.find(output_face); + labels != output_face_labels.end() && !labels->second.empty()) { + continue; + } + + add_adjacent_label_votes(output_face); + if (auto labels = output_face_labels.find(output_face); + labels != output_face_labels.end() && !labels->second.empty()) { + propagated_label = true; + } + } + } while (propagated_label); + for (auto output_face : arr_snap.face_handles()) { if (output_face->is_unbounded()) continue; const auto labels = output_face_labels.find(output_face);