Improvements in arrangement processing

CHANGELOG.md

@@ -7,6 +7,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+- Repair of non-manifold and self-intersecting roof junctions during arrangement snapping. This prevents certain geometric errors like non-watertight meshes, and highly non-planar faces that could prevriously happen in complex roof arrangements. See #168.
+
+### Fixed
+- Prevent roof-ground bow-ties by clipping roof faces against the terrain before extrusion.
+- Correct squared-distance tolerance comparisons during extrusion.
+
 ## [1.0.0] - 2026-04-20
 
 ## Added

apps/roofer-app/reconstruct_building.hpp

@@ -201,10 +201,13 @@ std::unordered_map<int, roofer::Mesh> extrude_lod22(
 
   auto& logger = roofer::logger::Logger::get_logger();
 
+  auto elevation_provider =
+      roofer::reconstruction::createElevationProvider(*building.h_ground);
+
   auto ArrangementDissolver =
       roofer::reconstruction::createArrangementDissolver();
   ArrangementDissolver->compute(
-      arrangement, SegmentRasteriser->heightfield,
+      arrangement, SegmentRasteriser->heightfield, *elevation_provider,
       {.dissolve_step_edges = dissolve_step_edges,
        .dissolve_all_interior = dissolve_all_interior,
        .step_height_threshold = cfg->lod13_step_height});
@@ -227,7 +230,7 @@ std::unordered_map<int, roofer::Mesh> extrude_lod22(
 
   auto ArrangementExtruder =
       roofer::reconstruction::createArrangementExtruder();
-  ArrangementExtruder->compute(arrangement, *building.h_ground,
+  ArrangementExtruder->compute(arrangement, *elevation_provider,
                                {.LoD2 = extrude_LoD2});
   // logger.debug("Completed ArrangementExtruder");
 #ifdef RF_USE_RERUN

flake.nix

@@ -191,6 +191,7 @@
               pkgsStatic.boost # need static for val3dity
               eigen
               fmt
+              catch2_3
 
               # val3dity
               pugixml

include/roofer/reconstruction/ArrangementBase.hpp

@@ -54,8 +54,6 @@ namespace roofer::reconstruction {
     FaceInfo operator()(const FaceInfo a, const FaceInfo b) const {
       auto r = FaceInfo();
       r.segid = 0;
-      // if (a.is_finite && b.is_finite)
-      r.is_finite = true;
 
       if (a.segid != 0 && b.segid == 0) {
         r = a;
@@ -105,18 +103,12 @@ namespace roofer::reconstruction {
           elevation(elevation),
           plane(plane) {
       CGAL_precondition(arr.is_empty());
-      for (auto uf = arr.unbounded_faces_begin();
-           uf != arr.unbounded_faces_end(); ++uf) {
-        uf->data().is_finite = false;
-      }
-      // arr.unbounded_face()->data().is_finite = false;
       n_faces++;
     };
     virtual void after_split_face(Face_handle old_face, Face_handle new_face,
                                   bool) {
       // Assign index to the new face.
       new_face->data().in_footprint = in_footprint;
-      new_face->data().is_finite = true;
       new_face->data().segid = plane_id;
       new_face->data().elevation_70p = elevation;
       new_face->data().plane = plane;

include/roofer/reconstruction/ArrangementDissolver.hpp

@@ -23,6 +23,7 @@
 #include <memory>
 #include <roofer/common/Raster.hpp>
 #include <roofer/common/datastructures.hpp>
+#include <roofer/reconstruction/ElevationProvider.hpp>
 #include <roofer/reconstruction/cgal_shared_definitions.hpp>
 
 namespace roofer::reconstruction {
@@ -51,6 +52,12 @@ namespace roofer::reconstruction {
     virtual ~ArrangementDissolverInterface() = default;
     virtual void compute(
         Arrangement_2& arrangement, const RasterTools::Raster& heightfield,
+        const ElevationProvider& elevation_provider,
+        ArrangementDissolverConfig config = ArrangementDissolverConfig()) = 0;
+
+    virtual void compute(
+        Arrangement_2& arrangement, const RasterTools::Raster& heightfield,
+        float base_elevation,
         ArrangementDissolverConfig config = ArrangementDissolverConfig()) = 0;
   };
 

include/roofer/reconstruction/ArrangementSnapper.hpp

@@ -28,6 +28,11 @@ namespace roofer::reconstruction {
 
   struct ArrangementSnapperConfig {
     float dist_thres = 0.05;
+    bool repair_non_manifold_vertices = true;
+    // manifold repairs with this radius, and down to half this radius in case
+    // clearance is limited. Below half this radius, repairs are skipped.
+    float manifold_repair_radius = 0.05;
+    float manifold_height_tolerance = 1e-4;
   };
 
   struct ArrangementSnapperInterface {

include/roofer/reconstruction/ElevationProvider.hpp

@@ -21,6 +21,8 @@
 
 #pragma once
 
+#include <CGAL/Polygon_with_holes_2.h>
+
 #include <roofer/reconstruction/cdt_util.hpp>
 #include <roofer/reconstruction/cgal_shared_definitions.hpp>
 
@@ -32,6 +34,10 @@ namespace roofer::reconstruction {
     virtual float get(const Point_2 pt) const = 0;
 
     virtual float get_percentile(float percentile) const = 0;
+
+    virtual std::vector<Segment_2> get_intersections(
+        const Plane& roof_plane,
+        const CGAL::Polygon_with_holes_2<EPECK>& bounds) const = 0;
   };
 
   std::unique_ptr<ElevationProvider> createElevationProvider(

include/roofer/reconstruction/cgal_shared_definitions.hpp

@@ -47,7 +47,6 @@ namespace roofer {
   typedef Traits_2::Point_2 Point_2;
 
   struct FaceInfo {
-    bool is_finite = false;
     bool is_ground = false;
     bool in_footprint = false;
     bool is_footprint_hole = false;

include/roofer/roofer.h

@@ -276,7 +276,7 @@ namespace roofer {
       auto ArrangementDissolver =
           roofer::reconstruction::createArrangementDissolver();
       ArrangementDissolver->compute(
-          arrangement, SegmentRasteriser->heightfield,
+          arrangement, SegmentRasteriser->heightfield, *elevation_provider,
           {.dissolve_step_edges = cfg.lod == 13,
            .dissolve_all_interior = cfg.lod == 12,
            .step_height_threshold = cfg.lod13_step_height});

src/core/reconstruction/ArrangementDissolver.cpp

@@ -24,15 +24,154 @@
 
 namespace roofer::reconstruction {
 
+  namespace {
+    using Polygon_2 = CGAL::Polygon_2<EPECK>;
+    using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2<EPECK>;
+
+    class Face_data_split_observer : public CGAL::Arr_observer<Arrangement_2> {
+     public:
+      explicit Face_data_split_observer(Arrangement_2& arrangement)
+          : CGAL::Arr_observer<Arrangement_2>(arrangement) {}
+
+      void after_split_face(Face_handle old_face, Face_handle new_face,
+                            bool) override {
+        new_face->set_data(old_face->data());
+      }
+    };
+
+    Polygon_with_holes_2 face_polygon(Face_handle face) {
+      Polygon_2 outer;
+      auto edge = face->outer_ccb();
+      auto first = edge;
+      do {
+        outer.push_back(edge->source()->point());
+        edge = edge->next();
+      } while (edge != first);
+
+      std::vector<Polygon_2> holes;
+      for (auto ccb = face->inner_ccbs_begin(); ccb != face->inner_ccbs_end();
+           ++ccb) {
+        Polygon_2 hole;
+        edge = *ccb;
+        first = edge;
+        do {
+          hole.push_back(edge->source()->point());
+          edge = edge->next();
+        } while (edge != first);
+        holes.push_back(std::move(hole));
+      }
+      return Polygon_with_holes_2(outer, holes.begin(), holes.end());
+    }
+
+    std::optional<Point_2> face_interior_point(Face_handle face) {
+      LinearRing polygon;
+      if (!arrangementface_to_polygon(face, polygon)) return std::nullopt;
+      auto triangulation = tri_util::create_from_polygon(polygon);
+      for (auto triangle = triangulation.finite_faces_begin();
+           triangle != triangulation.finite_faces_end(); ++triangle) {
+        if (!triangle->info().in_domain()) continue;
+        auto centroid = CGAL::centroid(triangulation.triangle(triangle));
+        return Point_2(centroid.x(), centroid.y());
+      }
+      return std::nullopt;
+    }
+
+    bool touches_footprint_boundary(Face_handle face) {
+      auto touches_boundary = [](Ccb_halfedge_circulator edge) {
+        auto first = edge;
+        do {
+          if (!edge->twin()->face()->data().in_footprint) return true;
+          edge = edge->next();
+        } while (edge != first);
+        return false;
+      };
+
+      if (touches_boundary(face->outer_ccb())) return true;
+      for (auto ccb = face->inner_ccbs_begin(); ccb != face->inner_ccbs_end();
+           ++ccb) {
+        if (touches_boundary(*ccb)) return true;
+      }
+      return false;
+    }
+
+    bool touches_unbounded_face(Face_handle face) {
+      auto touches_unbounded = [](Ccb_halfedge_circulator edge) {
+        auto first = edge;
+        do {
+          if (edge->twin()->face()->is_unbounded()) return true;
+          edge = edge->next();
+        } while (edge != first);
+        return false;
+      };
+
+      if (touches_unbounded(face->outer_ccb())) return true;
+      for (auto ccb = face->inner_ccbs_begin(); ccb != face->inner_ccbs_end();
+           ++ccb) {
+        if (touches_unbounded(*ccb)) return true;
+      }
+      return false;
+    }
+
+    double roof_height(const Plane& plane, const Point_2& point) {
+      return -(plane.a() * CGAL::to_double(point.x()) +
+               plane.b() * CGAL::to_double(point.y()) + plane.d()) /
+             plane.c();
+    }
+
+    void clip_roof_faces_to_terrain(
+        Arrangement_2& arrangement,
+        const ElevationProvider& elevation_provider) {
+      std::vector<Segment_2> intersections;
+      for (auto face : arrangement.face_handles()) {
+        if (face->is_unbounded() || !face->data().in_footprint ||
+            face->data().segid == 0) {
+          continue;
+        }
+        auto face_intersections = elevation_provider.get_intersections(
+            face->data().plane, face_polygon(face));
+        intersections.insert(intersections.end(), face_intersections.begin(),
+                             face_intersections.end());
+      }
+
+      Face_data_split_observer observer(arrangement);
+      for (const auto& intersection : intersections) {
+        CGAL::insert(arrangement, intersection);
+      }
+
+      constexpr double elevation_tolerance = 1e-6;
+      std::vector<std::pair<Face_handle, bool>> ground_faces;
+      for (auto face : arrangement.face_handles()) {
+        if (face->is_unbounded() || !face->data().in_footprint ||
+            face->data().segid == 0) {
+          continue;
+        }
+        auto sample = face_interior_point(face);
+        if (!sample) continue;
+        if (roof_height(face->data().plane, *sample) <
+            elevation_provider.get(*sample) - elevation_tolerance) {
+          ground_faces.emplace_back(face, touches_footprint_boundary(face));
+        }
+      }
+      for (auto [face, touches_boundary] : ground_faces) {
+        face->data().in_footprint = false;
+        face->data().is_ground = true;
+        face->data().is_footprint_hole = !touches_boundary;
+      }
+    }
+  }  // namespace
+
   class ArrangementDissolver : public ArrangementDissolverInterface {
    public:
     void compute(Arrangement_2& arr, const RasterTools::Raster& heightfield,
+                 const ElevationProvider& elevation_provider,
                  ArrangementDissolverConfig cfg) override {
       if (cfg.dissolve_seg_edges) {
         Face_merge_observer obs(arr);
         arr_dissolve_seg_edges(arr);
       }
 
+      clip_roof_faces_to_terrain(arr, elevation_provider);
+
       if (cfg.dissolve_all_interior) {
         arr_dissolve_fp(arr, true, false);
       }
@@ -113,7 +252,9 @@ namespace roofer::reconstruction {
       for (auto fh : arr.face_handles()) {
         if (fh != f_unb)
           if (!fh->data().in_footprint && !fh->data().is_footprint_hole) {
-            fh->data().is_footprint_hole = true;
+            if (!fh->data().is_ground || !touches_unbounded_face(fh)) {
+              fh->data().is_footprint_hole = true;
+            }
           }
       }
 
@@ -201,6 +342,13 @@ namespace roofer::reconstruction {
       // }
       // output("arrangement").set(arr);
     }
+
+    void compute(Arrangement_2& arr, const RasterTools::Raster& heightfield,
+                 float base_elevation,
+                 ArrangementDissolverConfig cfg) override {
+      auto elevation_provider = createElevationProvider(base_elevation);
+      compute(arr, heightfield, *elevation_provider, cfg);
+    }
   };
 
   std::unique_ptr<ArrangementDissolverInterface> createArrangementDissolver() {

src/core/reconstruction/ArrangementExtruder.cpp

@@ -76,12 +76,12 @@ namespace roofer::reconstruction {
       T& ring, Halfedge_handle hedge,
       std::unordered_map<Vertex_handle, std::vector<hf_pair>>& vertex_columns,
       std::unordered_map<Halfedge_handle, EPECK::Point_3>& extra_wall_points,
-      float& snap_tolerance) {
+      const float snap_tolerance, const float snap_tolerance_sq) {
     auto first = hedge;
     do {
       auto v = hedge->source();
       if (CGAL::squared_distance(v->point(), hedge->target()->point()) >
-          snap_tolerance) {
+          snap_tolerance_sq) {
         for (auto& [h, f_h] : vertex_columns[v]) {
           if (f_h == hedge->face()) {
             ring.push_back(v2p(v, h));
@@ -107,6 +107,7 @@ namespace roofer::reconstruction {
                  ArrangementExtruderConfig cfg) override {
       typedef Arrangement_2::Traits_2 AT;
       float snap_tolerance = std::pow(10, -cfg.snap_tolerance_exp);
+      float snap_tolerance_sq = snap_tolerance * snap_tolerance;
 
       // assume we have only one unbounded faces that just has the building
       // footprint as a hole
@@ -131,7 +132,7 @@ namespace roofer::reconstruction {
           do {
             if (CGAL::squared_distance(he->source()->point(),
                                        he->target()->point()) >
-                snap_tolerance) {
+                snap_tolerance_sq) {
               float pt_elevation =
                   elevation_provider.get(he->source()->point());
               floor.push_back(v2p(he->source(), pt_elevation));
@@ -157,7 +158,7 @@ namespace roofer::reconstruction {
             do {
               if (CGAL::squared_distance(he->source()->point(),
                                          he->target()->point()) >
-                  snap_tolerance) {
+                  snap_tolerance_sq) {
                 float pt_elevation =
                     elevation_provider.get(he->source()->point());
                 hole.push_back(v2p(he->source(), pt_elevation));
@@ -254,10 +255,10 @@ namespace roofer::reconstruction {
           auto& pl_b = f_b->data().plane;
 
           // edge has no length
-          if (CGAL::squared_distance(p1, p2) < snap_tolerance) continue;
+          if (CGAL::squared_distance(p1, p2) < snap_tolerance_sq) continue;
 
           // edge is not in footprint nor on its boundary
-          if (!fp_a & !fp_b) {
+          if (!fp_a && !fp_b) {
             continue;
           }
 
@@ -465,14 +466,14 @@ namespace roofer::reconstruction {
             auto he = face->outer_ccb();
 
             push_ccb(roofpart, he, vertex_columns, extra_wall_points,
-                     snap_tolerance);
+                     snap_tolerance, snap_tolerance_sq);
 
             for (Arrangement_2::Hole_iterator hole = face->holes_begin();
                  hole != face->holes_end(); ++hole) {
               vec3f roofpart_hole;
               auto he = *hole;
               push_ccb(roofpart_hole, he, vertex_columns, extra_wall_points,
-                       snap_tolerance);
+                       snap_tolerance, snap_tolerance_sq);
               if (roofpart_hole.size() > 2) {
                 roofpart.interior_rings().push_back(roofpart_hole);
               }