Make roofer deterministic: seed, address-independent orders, uninitialised reads

include/roofer/reconstruction/LineRegulariserBase.hpp

@@ -45,6 +45,9 @@ namespace roofer::linereg {
 
   template <typename T>
   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<linetype*> lines;
@@ -60,6 +63,17 @@ namespace roofer::linereg {
     void calc_mean_value() override;
   };
 
+  // orders cluster handles by their creation order
+  struct ClusterHLess {
+    template <typename ClusterH>
+    bool operator()(const ClusterH& lhs, const ClusterH& rhs) const {
+      return lhs->id < rhs->id;
+    }
+  };
+
+  template <typename ClusterH>
+  using ClusterSet = std::set<ClusterH, ClusterHLess>;
+
   template <typename ClusterH>
   struct DistanceTable {
     using ClusterPair = std::pair<ClusterH, ClusterH>;
@@ -103,9 +117,10 @@ namespace roofer::linereg {
 
     Cluster2DistPairMap cluster_to_dist_pairs;
     DistanceHeap distances;
-    std::set<ClusterH>& clusters;
+    ClusterSet<ClusterH>& clusters;
 
-    DistanceTable(std::set<ClusterH>& clusters);  // computes initial distances
+    DistanceTable(
+        ClusterSet<ClusterH>& 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<size_t, SegmentVec> segments;
-    std::set<AngleClusterH> angle_clusters;
-    std::set<DistClusterH> dist_clusters;
+    ClusterSet<AngleClusterH> angle_clusters;
+    ClusterSet<DistClusterH> dist_clusters;
 
     LineRegulariser() = default;
 

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;
       }

src/core/reconstruction/ArrangementBase.cpp

@@ -1,3 +1,3 @@
 // Copyright (c) 2018-2026 TU Delft 3D geoinformation group, Ravi Peters (3DGI),
 // and Balazs Dukai (3DGI)
 
@@ -19,7 +19,7 @@
 // Author(s):
 // Ravi Peters
 
 #include <roofer/reconstruction/ArrangementBase.hpp>
 #include <stack>
 
 namespace roofer::reconstruction {
@@ -175,19 +175,27 @@
                        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<FacePair> 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) {

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<Vertex_handle> 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<T::Vertex_handle, Arrangement_2::Vertex_handle>
             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<T::Vertex_handle, std::size_t> 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<std::pair<T::Vertex_handle, T::Vertex_handle>> 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());

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<linereg::linetype*> 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};

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 <typename ClusterH>
-  DistanceTable<ClusterH>::DistanceTable(std::set<ClusterH>& clusters)
+  DistanceTable<ClusterH>::DistanceTable(ClusterSet<ClusterH>& 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<AngleCluster>();
+      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<DistClusterH> dclusters;
+      ClusterSet<DistClusterH> dclusters;
       for (auto& line : aclusterh->lines) {
         auto dclusterh = std::make_shared<DistCluster>();
+        dclusterh->id = next_id++;
         dclusterh->value = line->segment;
         // dclusterh->has_intersection_line = line->priority==2;
         dclusterh->lines.push_back(line);

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<int>("classification");

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) {