From 60a015ae5bbca608723280e64069e1a1dc219eae Mon Sep 17 00:00:00 2001 From: Pierre-Louis Guhur Date: Sat, 3 Oct 2026 18:03:41 +0200 Subject: [PATCH 1/6] fix(reconstruction): seed the region grower's shuffle, so runs are reproducible --- include/roofer/reconstruction/RegionGrower_DS_CGAL.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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; } From f785cecf4a35ae3eaf2e55c77335d6fdf5c867dc Mon Sep 17 00:00:00 2001 From: Pierre-Louis Guhur Date: Sun, 4 Oct 2026 11:03:12 +0200 Subject: [PATCH 2/6] fix(reconstruction): order line clusters and step merges by creation, not by address Several containers keyed by pointers or CGAL handles were iterated in hash or address order, so a building's model depended on where its objects landed in memory: the output changed with the length of the output path, and from run to run when the I/O threads' allocations shifted the heap. - line regularisation: clusters get a creation id; their sets are ordered by it, which fixes the distance heap's tie order, the merge direction and the cluster ids handed to the lines; - step merging: face pairs are visited in the order the arrangement's edges first meet them; Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LZBFfd55Us1zg7chCUDs5C --- .../reconstruction/LineRegulariserBase.hpp | 23 +++++++++++++++---- src/core/reconstruction/ArrangementBase.cpp | 12 ++++++++-- .../reconstruction/LineRegulariserBase.cpp | 9 ++++++-- 3 files changed, 36 insertions(+), 8 deletions(-) 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/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/LineRegulariserBase.cpp b/src/core/reconstruction/LineRegulariserBase.cpp index 30209e07..b2051291 100644 --- a/src/core/reconstruction/LineRegulariserBase.cpp +++ b/src/core/reconstruction/LineRegulariserBase.cpp @@ -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); From 0f73222fb984a97d8078329d7452747847a9cd37 Mon Sep 17 00:00:00 2001 From: Pierre-Louis Guhur Date: Sun, 4 Oct 2026 11:28:01 +0200 Subject: [PATCH 3/6] fix(reconstruction): bound a lone intersection line by itself, not by uninitialised extremes When a distance cluster holds only intersection lines, there is no other line to bound them, and calc_segment read its uninitialised pmin/pmax/dmin/dmax: the segment took whatever was left in memory, and the model changed from run to run. The intersection line now bounds itself, and calc_segment starts from the centroid so that an empty list can no longer read garbage. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LZBFfd55Us1zg7chCUDs5C --- src/core/reconstruction/LineRegulariser.cpp | 9 +++++++-- src/core/reconstruction/LineRegulariserBase.cpp | 4 ++-- 2 files changed, 9 insertions(+), 4 deletions(-) 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 b2051291..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); From efb5d135b66959cc8ccdadc60d7b93f082d25690 Mon Sep 17 00:00:00 2001 From: Pierre-Louis Guhur Date: Sun, 4 Oct 2026 16:28:51 +0200 Subject: [PATCH 4/6] fix(reconstruction): snap in an order that does not depend on memory addresses CGAL's edge iterator reports each edge from whichever of its two faces has the lower address, so the snapper's "first short edge" and the order in which the snapped arrangement was rebuilt depended on the heap layout. The snapper now collapses the shortest short edge (ties broken by its endpoints' coordinates), inserts the constrained edges in the order of their endpoints' vertex indices, and seeds the dangling-constraint peeling in vertex order rather than in the order of a hash map keyed by vertex handles. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LZBFfd55Us1zg7chCUDs5C --- .../reconstruction/ArrangementSnapper.cpp | 65 ++++++++++++++++--- 1 file changed, 57 insertions(+), 8 deletions(-) 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()); From 3921cb6c8903bf3c1eac00ba1c46209af9aa61ba Mon Sep 17 00:00:00 2001 From: Pierre-Louis Guhur Date: Sun, 4 Oct 2026 16:42:09 +0200 Subject: [PATCH 5/6] fix(cropper): initialise the loop index that merges buffer ground points The loop that merges each footprint's buffer ground points into its point cloud started from an uninitialised index, so whether those points were merged depended on whatever the stack held. It now starts at 0, so they are always merged, as the comment above the loop intends. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LZBFfd55Us1zg7chCUDs5C --- src/extra/io/StreamCropper.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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"); From 093908d73b8c41e702e0decf14a90fe08a335f3c Mon Sep 17 00:00:00 2001 From: Pierre-Louis Guhur Date: Sun, 4 Oct 2026 16:47:21 +0200 Subject: [PATCH 6/6] fix(thinning): skip points one past the count raster instead of reading past it gridthinPointcloud skipped points whose cell index was greater than the count raster's size, but not the one equal to it, so a point mapping to that cell read one element past the array. Whether it was kept, and whether a random draw was consumed, then depended on that memory, and so did every following point's draw. Footprint buffer ground points, which lie outside the footprint's raster, hit this case. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LZBFfd55Us1zg7chCUDs5C --- src/extra/misc/PointcloudRasteriser.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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) {