diff --git a/include/roofer/reconstruction/LineRegulariserBase.hpp b/include/roofer/reconstruction/LineRegulariserBase.hpp index f5553ce6..97d65d95 100644 --- a/include/roofer/reconstruction/LineRegulariserBase.hpp +++ b/include/roofer/reconstruction/LineRegulariserBase.hpp @@ -45,6 +45,9 @@ namespace roofer::linereg { template struct Cluster { + // creation order: sets of clusters are ordered by it, not by address, so + // that the clustering does not depend on where the clusters land in memory + size_t id = 0; T value; bool has_intersection_line; std::vector lines; @@ -60,6 +63,17 @@ namespace roofer::linereg { void calc_mean_value() override; }; + // orders cluster handles by their creation order + struct ClusterHLess { + template + bool operator()(const ClusterH& lhs, const ClusterH& rhs) const { + return lhs->id < rhs->id; + } + }; + + template + using ClusterSet = std::set; + template struct DistanceTable { using ClusterPair = std::pair; @@ -103,9 +117,10 @@ namespace roofer::linereg { Cluster2DistPairMap cluster_to_dist_pairs; DistanceHeap distances; - std::set& clusters; + ClusterSet& clusters; - DistanceTable(std::set& clusters); // computes initial distances + DistanceTable( + ClusterSet& clusters); // computes initial distances void merge( ClusterH lhs, ClusterH rhs); // merges two clusters, then removes one from the @@ -166,8 +181,8 @@ namespace roofer::linereg { double angle_threshold, dist_threshold; std::unordered_map segments; - std::set angle_clusters; - std::set dist_clusters; + ClusterSet angle_clusters; + ClusterSet dist_clusters; LineRegulariser() = default; diff --git a/include/roofer/reconstruction/RegionGrower_DS_CGAL.hpp b/include/roofer/reconstruction/RegionGrower_DS_CGAL.hpp index 9430a8fd..bc69f06d 100644 --- a/include/roofer/reconstruction/RegionGrower_DS_CGAL.hpp +++ b/include/roofer/reconstruction/RegionGrower_DS_CGAL.hpp @@ -90,8 +90,8 @@ namespace roofer { for (size_t i = 0; i < size; ++i) { seeds.push_back(i); } - std::random_device rd; - std::mt19937 g(rd()); + // a fixed seed: the same input gives the same planes, run after run + std::mt19937 g(31415); std::shuffle(seeds.begin(), seeds.end(), g); return seeds; } diff --git a/src/core/reconstruction/ArrangementBase.cpp b/src/core/reconstruction/ArrangementBase.cpp index 0cf32804..09ddd01a 100644 --- a/src/core/reconstruction/ArrangementBase.cpp +++ b/src/core/reconstruction/ArrangementBase.cpp @@ -175,19 +175,27 @@ namespace roofer::reconstruction { KeyHash, KeyEqual> step_boundaries; + // face pairs in the order the arrangement's edges first meet them: the hash + // map's own order follows memory addresses, and decides which of two equal + // steps is removed first + std::vector step_order; + while (true) { double d_min = step_height_threshold; step_boundaries.clear(); + step_order.clear(); for (auto& edge : arr.edge_handles()) { auto f1 = edge->face(); auto f2 = edge->twin()->face(); if ((f1->data().in_footprint && f2->data().in_footprint) && (f1->data().segid != 0 && f2->data().segid != 0)) { - step_boundaries[FacePair(f1, f2)].push_back(edge); + auto [it, inserted] = step_boundaries.try_emplace(FacePair(f1, f2)); + if (inserted) step_order.push_back(it->first); + it->second.push_back(edge); } } FacePair facepair_min; - for (auto& [faces, edges] : step_boundaries) { + for (auto& faces : step_order) { double d = std::abs(faces.f_hi->data().elevation_50p - faces.f_lo->data().elevation_50p); if (d < d_min) { diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index 14620409..1a04fdc9 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -471,9 +471,15 @@ namespace roofer::reconstruction { adjacency[second].insert(first); } + // seed the leaves in vertex order, not in the hash map's, which follows + // memory addresses: the order constraints are removed in can change the + // triangulation std::queue leaves; - for (const auto& [vertex, neighbours] : adjacency) { - if (neighbours.size() == 1) leaves.push(vertex); + for (auto vertex = tri.finite_vertices_begin(); + vertex != tri.finite_vertices_end(); ++vertex) { + if (auto it = adjacency.find(vertex); + it != adjacency.end() && it->second.size() == 1) + leaves.push(vertex); } while (!leaves.empty()) { @@ -1265,17 +1271,44 @@ namespace roofer::reconstruction { bool found_short_edge; do { found_short_edge = false; + // collapse the shortest short edge, ties broken by its endpoints: + // CGAL's edge iterator reports each edge from whichever of its two + // faces sits lower in memory, so "the first short edge" used to + // depend on the heap layout + auto lex_less = [](const T::Point_2& a, const T::Point_2& b) { + return a.x() < b.x() || (a.x() == b.x() && a.y() < b.y()); + }; + bool have_shortest = false; + Vertex_handle shortest1, shortest2; + double shortest_d = 0; for (Finite_edges_iterator ceit = tri.finite_edges_begin(); ceit != tri.finite_edges_end(); ++ceit) { - auto v1 = ceit->first->vertex(tri.cw(ceit->second)); - auto v2 = ceit->first->vertex(tri.ccw(ceit->second)); + auto a = ceit->first->vertex(tri.cw(ceit->second)); + auto b = ceit->first->vertex(tri.ccw(ceit->second)); + if (lex_less(b->point(), a->point())) std::swap(a, b); + double d = CGAL::squared_distance(a->point(), b->point()); + if (d >= sq_dist_thres) continue; + if (!have_shortest || d < shortest_d || + (d == shortest_d && + (lex_less(a->point(), shortest1->point()) || + (a->point() == shortest1->point() && + lex_less(b->point(), shortest2->point()))))) { + have_shortest = true; + shortest_d = d; + shortest1 = a; + shortest2 = b; + } + } + if (have_shortest) { + auto v1 = shortest1; + auto v2 = shortest2; auto& p1 = v1->point(); auto& p2 = v2->point(); // do not collapse if this edge is on the footprint boundary // if (v1->info() && v2->info()) continue; - if (CGAL::squared_distance(p1, p2) < sq_dist_thres) { + { // std::cout << "short edge between " << p1 << " and " << p2 // << std::endl; @@ -1309,7 +1342,6 @@ namespace roofer::reconstruction { unbounded_label); found_short_edge = true; - break; } } } while (found_short_edge); @@ -1407,8 +1439,13 @@ namespace roofer::reconstruction { std::unordered_map vertex2arr_map; + // vertices in iteration order, which is deterministic; the edges are + // inserted in the order of their endpoints' indices, not in the order + // CGAL's edge iterator reports them, which follows memory addresses + std::unordered_map vertex_index; for (auto vh = tri.finite_vertices_begin(); vh != tri.finite_vertices_end(); ++vh) { + vertex_index[vh] = vertex_index.size(); // make sure not to add isolated vertices if (tri.are_there_incident_constraints(vh)) { vertex2arr_map[vh] = insert_point( @@ -1417,9 +1454,21 @@ namespace roofer::reconstruction { } } + std::vector> cedges; for (auto ce : tri.constrained_edges()) { - auto v1 = ce.first->vertex(tri.cw(ce.second)); - auto v2 = ce.first->vertex(tri.ccw(ce.second)); + auto a = ce.first->vertex(tri.cw(ce.second)); + auto b = ce.first->vertex(tri.ccw(ce.second)); + if (vertex_index[b] < vertex_index[a]) std::swap(a, b); + cedges.emplace_back(a, b); + } + std::sort(cedges.begin(), cedges.end(), + [&vertex_index](const auto& l, const auto& r) { + return std::make_pair(vertex_index[l.first], + vertex_index[l.second]) < + std::make_pair(vertex_index[r.first], + vertex_index[r.second]); + }); + for (auto& [v1, v2] : cedges) { auto& p1_ = v1->point(); auto& p2_ = v2->point(); auto p1 = Arrangement_2::Point_2(p1_.x(), p1_.y()); diff --git a/src/core/reconstruction/LineRegulariser.cpp b/src/core/reconstruction/LineRegulariser.cpp index f6a18ec5..1455ace6 100644 --- a/src/core/reconstruction/LineRegulariser.cpp +++ b/src/core/reconstruction/LineRegulariser.cpp @@ -93,8 +93,13 @@ namespace roofer::reconstruction { for (auto line : prio_lines) { double mean_angle = line->angle; auto centroid = line->midpoint; - segment = linereg::calc_segment(centroid, mean_angle, other_lines, - cfg.extension); + // a cluster of intersection lines only has no other line to bound + // them: calc_segment would then read its uninitialised extremes, so + // the intersection line bounds itself + std::vector own_line{line}; + segment = linereg::calc_segment( + centroid, mean_angle, + other_lines.empty() ? own_line : other_lines, cfg.extension); auto new_seg = Segment(); new_seg[0] = {float(CGAL::to_double(segment.source().x())), float(CGAL::to_double(segment.source().y())), 0}; diff --git a/src/core/reconstruction/LineRegulariserBase.cpp b/src/core/reconstruction/LineRegulariserBase.cpp index 30209e07..e84fbc2b 100644 --- a/src/core/reconstruction/LineRegulariserBase.cpp +++ b/src/core/reconstruction/LineRegulariserBase.cpp @@ -59,8 +59,8 @@ namespace roofer::linereg { lv = lv / CGAL::sqrt(CGAL::to_double(lv.squared_length())); bool setminmax = false; - Point_2 pmin, pmax; - double dmin, dmax; + Point_2 pmin = centroid, pmax = centroid; + double dmin = 0, dmax = 0; for (auto& line : lines) { auto p = line->segment.source(); auto d = CGAL::to_double(Vector_2(p.x(), p.y()) * lv); @@ -112,7 +112,7 @@ namespace roofer::linereg { } template - DistanceTable::DistanceTable(std::set& clusters) + DistanceTable::DistanceTable(ClusterSet& clusters) : clusters(clusters) { // compute only half of the distance table, since the other half will be // exactly the same @@ -192,8 +192,10 @@ namespace roofer::linereg { void LineRegulariser::perform_angle_clustering() { // make clusters angle_clusters.clear(); + size_t next_id = 0; for (auto& line : lines) { auto aclusterh = std::make_shared(); + aclusterh->id = next_id++; aclusterh->value = line.angle; // aclusterh->has_intersection_line = line.priority==2; aclusterh->lines.push_back(&line); @@ -225,10 +227,13 @@ namespace roofer::linereg { dist_clusters.clear(); // perform distance clustering for each angle cluster + // ids are unique across all angle clusters, as dist_clusters gathers them + size_t next_id = 0; for (auto& aclusterh : angle_clusters) { - std::set dclusters; + ClusterSet dclusters; for (auto& line : aclusterh->lines) { auto dclusterh = std::make_shared(); + dclusterh->id = next_id++; dclusterh->value = line->segment; // dclusterh->has_intersection_line = line->priority==2; dclusterh->lines.push_back(line); diff --git a/src/extra/io/StreamCropper.cpp b/src/extra/io/StreamCropper.cpp index 4f45c80d..1d02a0b3 100644 --- a/src/extra/io/StreamCropper.cpp +++ b/src/extra/io/StreamCropper.cpp @@ -501,7 +501,7 @@ namespace roofer::io { // merge buffer ground points into regular point_clouds now that the // proper counts have been established - for (size_t poly_i; poly_i < polygons.size(); poly_i++) { + for (size_t poly_i = 0; poly_i < polygons.size(); poly_i++) { auto& point_cloud = point_clouds.at(poly_i); auto classification = point_cloud.attributes.get_if("classification"); diff --git a/src/extra/misc/PointcloudRasteriser.cpp b/src/extra/misc/PointcloudRasteriser.cpp index 2845e533..90a3cb51 100644 --- a/src/extra/misc/PointcloudRasteriser.cpp +++ b/src/extra/misc/PointcloudRasteriser.cpp @@ -185,7 +185,7 @@ namespace roofer::misc { auto& p = pointcloud[pi]; auto& c = (*classification)[pi]; auto i = getLinearCoord(cnt_image, p[0], p[1]); - if (i < 0 || i > cnt_image.array.size()) continue; + if (i < 0 || i >= cnt_image.array.size()) continue; // check if we have too many points in this cell if (cnt_image.array[i] > max_cnt_per_cell) {