From 9e6cc805c7b9d2b180d7541bf668cc5b49afc15c Mon Sep 17 00:00:00 2001 From: ylannl Date: Wed, 27 May 2026 13:21:13 +0200 Subject: [PATCH 01/21] properly interpret sq tolerances --- src/core/reconstruction/ArrangementExtruder.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/core/reconstruction/ArrangementExtruder.cpp b/src/core/reconstruction/ArrangementExtruder.cpp index ec6099e4..b1f8e8b2 100644 --- a/src/core/reconstruction/ArrangementExtruder.cpp +++ b/src/core/reconstruction/ArrangementExtruder.cpp @@ -76,12 +76,12 @@ namespace roofer::reconstruction { T& ring, Halfedge_handle hedge, std::unordered_map>& vertex_columns, std::unordered_map& 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); } From d2fa7b1ab6cf00b91fd4ca7609ec6345cad672c1 Mon Sep 17 00:00:00 2001 From: ylannl Date: Wed, 10 Jun 2026 18:25:13 +0200 Subject: [PATCH 02/21] add catch2 to dev shell --- flake.nix | 1 + 1 file changed, 1 insertion(+) diff --git a/flake.nix b/flake.nix index f0bdd7c5..99101373 100644 --- a/flake.nix +++ b/flake.nix @@ -191,6 +191,7 @@ pkgsStatic.boost # need static for val3dity eigen fmt + catch2_3 # val3dity pugixml From 6bd0b825bf769eda3064603c95aeb123569cae83 Mon Sep 17 00:00:00 2001 From: ylannl Date: Wed, 10 Jun 2026 22:26:54 +0200 Subject: [PATCH 03/21] Do not create roof-ground bowties. Fixes #17 --- apps/roofer-app/reconstruct_building.hpp | 7 +- .../reconstruction/ArrangementDissolver.hpp | 7 + .../reconstruction/ElevationProvider.hpp | 6 + include/roofer/roofer.h | 2 +- .../reconstruction/ArrangementDissolver.cpp | 151 +++++++++++++++++- src/core/reconstruction/ElevationProvider.cpp | 142 +++++++++++++++- tests/CMakeLists.txt | 6 + tests/test_arrangement_dissolver.cpp | 93 +++++++++++ 8 files changed, 407 insertions(+), 7 deletions(-) create mode 100644 tests/test_arrangement_dissolver.cpp diff --git a/apps/roofer-app/reconstruct_building.hpp b/apps/roofer-app/reconstruct_building.hpp index b0e7a822..2794b9e4 100644 --- a/apps/roofer-app/reconstruct_building.hpp +++ b/apps/roofer-app/reconstruct_building.hpp @@ -201,10 +201,13 @@ std::unordered_map 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 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 diff --git a/include/roofer/reconstruction/ArrangementDissolver.hpp b/include/roofer/reconstruction/ArrangementDissolver.hpp index f4e4a328..ef3357a0 100644 --- a/include/roofer/reconstruction/ArrangementDissolver.hpp +++ b/include/roofer/reconstruction/ArrangementDissolver.hpp @@ -23,6 +23,7 @@ #include #include #include +#include #include 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; }; diff --git a/include/roofer/reconstruction/ElevationProvider.hpp b/include/roofer/reconstruction/ElevationProvider.hpp index ab30ef52..9e5f5479 100644 --- a/include/roofer/reconstruction/ElevationProvider.hpp +++ b/include/roofer/reconstruction/ElevationProvider.hpp @@ -21,6 +21,8 @@ #pragma once +#include + #include #include @@ -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 get_intersections( + const Plane& roof_plane, + const CGAL::Polygon_with_holes_2& bounds) const = 0; }; std::unique_ptr createElevationProvider( diff --git a/include/roofer/roofer.h b/include/roofer/roofer.h index 5ad5e584..43eae0b0 100644 --- a/include/roofer/roofer.h +++ b/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}); diff --git a/src/core/reconstruction/ArrangementDissolver.cpp b/src/core/reconstruction/ArrangementDissolver.cpp index d7fc4d18..edb32566 100644 --- a/src/core/reconstruction/ArrangementDissolver.cpp +++ b/src/core/reconstruction/ArrangementDissolver.cpp @@ -24,15 +24,155 @@ namespace roofer::reconstruction { + namespace { + using Polygon_2 = CGAL::Polygon_2; + using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; + + class Face_data_split_observer : public CGAL::Arr_observer { + public: + explicit Face_data_split_observer(Arrangement_2& arrangement) + : CGAL::Arr_observer(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 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 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 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> 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 +253,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 +343,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 createArrangementDissolver() { diff --git a/src/core/reconstruction/ElevationProvider.cpp b/src/core/reconstruction/ElevationProvider.cpp index 3b391d93..59f0577f 100644 --- a/src/core/reconstruction/ElevationProvider.cpp +++ b/src/core/reconstruction/ElevationProvider.cpp @@ -17,18 +17,124 @@ // . // Author(s): -// Ivan Paden +// Ivan Paden, Ravi Peters #include #include +#include + namespace roofer::reconstruction { + namespace { + using Polygon_2 = CGAL::Polygon_2; + using Polygon_with_holes_2 = CGAL::Polygon_with_holes_2; + + bool is_inside(const Point_2& point, const Polygon_with_holes_2& polygon) { + if (polygon.outer_boundary().bounded_side(point) != + CGAL::ON_BOUNDED_SIDE) { + return false; + } + for (auto hole = polygon.holes_begin(); hole != polygon.holes_end(); + ++hole) { + if (hole->bounded_side(point) != CGAL::ON_UNBOUNDED_SIDE) { + return false; + } + } + return true; + } + + template + void add_boundary_intersections(const Curve& curve, + const Polygon_2& boundary, + std::vector& points) { + for (auto edge = boundary.edges_begin(); edge != boundary.edges_end(); + ++edge) { + auto result = CGAL::intersection(curve, *edge); + if (!result) continue; + if (const auto* point = std::get_if(&*result)) { + points.push_back(*point); + } + } + } + + template + void add_boundary_intersections(const Curve& curve, + const Polygon_with_holes_2& polygon, + std::vector& points) { + add_boundary_intersections(curve, polygon.outer_boundary(), points); + for (auto hole = polygon.holes_begin(); hole != polygon.holes_end(); + ++hole) { + add_boundary_intersections(curve, *hole, points); + } + } + + std::vector clip_line(const EPECK::Line_2& line, + const Polygon_with_holes_2& polygon) { + std::vector points; + add_boundary_intersections(line, polygon, points); + + auto direction = line.to_vector(); + std::sort( + points.begin(), points.end(), + [&](const Point_2& lhs, const Point_2& rhs) { + auto lhs_t = lhs.x() * direction.x() + lhs.y() * direction.y(); + auto rhs_t = rhs.x() * direction.x() + rhs.y() * direction.y(); + return lhs_t < rhs_t; + }); + + std::vector result; + for (size_t i = 1; i < points.size(); ++i) { + if (points[i - 1] == points[i]) continue; + Segment_2 segment(points[i - 1], points[i]); + if (is_inside(CGAL::midpoint(segment.source(), segment.target()), + polygon)) { + result.push_back(segment); + } + } + return result; + } + + std::vector clip_segment(const Segment_2& segment, + const Polygon_with_holes_2& polygon) { + std::vector points = {segment.source(), segment.target()}; + add_boundary_intersections(segment, polygon, points); + std::sort(points.begin(), points.end(), + [&](const Point_2& lhs, const Point_2& rhs) { + return CGAL::squared_distance(segment.source(), lhs) < + CGAL::squared_distance(segment.source(), rhs); + }); + + std::vector result; + for (size_t i = 1; i < points.size(); ++i) { + if (points[i - 1] == points[i]) continue; + Segment_2 clipped(points[i - 1], points[i]); + if (is_inside(CGAL::midpoint(clipped.source(), clipped.target()), + polygon)) { + result.push_back(clipped); + } + } + return result; + } + + std::optional triangle_intersection( + const Plane& roof_plane, const proj_tri_util::DT::Triangle& triangle) { + auto intersection = CGAL::intersection(roof_plane, triangle); + if (!intersection) return std::nullopt; + + const auto* segment = std::get_if(&*intersection); + if (!segment) return std::nullopt; + + return Segment_2(Point_2(segment->source().x(), segment->source().y()), + Point_2(segment->target().x(), segment->target().y())); + } + } // namespace + struct ConstantElevationProvider : public ElevationProvider { const float floor_elevation_; ConstantElevationProvider(const float floor_elevation) - : floor_elevation_(floor_elevation){}; + : floor_elevation_(floor_elevation) {}; virtual float get(const Point_2 /* pt */) const override { return floor_elevation_; @@ -37,13 +143,28 @@ namespace roofer::reconstruction { virtual float get_percentile(float /* percentile */) const override { return floor_elevation_; } + + std::vector get_intersections( + const Plane& roof_plane, + const CGAL::Polygon_with_holes_2& bounds) const override { + constexpr double epsilon = 1e-8; + // check if plane is approximately horizontal + if (std::abs(roof_plane.a()) < epsilon && + std::abs(roof_plane.b()) < epsilon) { + return {}; + } + EPECK::Line_2 intersection( + roof_plane.a(), roof_plane.b(), + roof_plane.c() * floor_elevation_ + roof_plane.d()); + return clip_line(intersection, bounds); + } }; struct InterpolatedElevationProvider : public ElevationProvider { std::shared_ptr base_cdt_ptr_; InterpolatedElevationProvider(const proj_tri_util::DT& base_cdt) - : base_cdt_ptr_(std::make_shared(base_cdt)){}; + : base_cdt_ptr_(std::make_shared(base_cdt)) {}; virtual float get(const Point_2 pt) const override { return proj_tri_util::interpolate_from_cdt(pt, *base_cdt_ptr_); @@ -63,6 +184,21 @@ namespace roofer::reconstruction { return compute_percentile(elevations, percentile); } + std::vector get_intersections( + const Plane& roof_plane, + const CGAL::Polygon_with_holes_2& bounds) const override { + std::vector result; + for (auto face = base_cdt_ptr_->finite_faces_begin(); + face != base_cdt_ptr_->finite_faces_end(); ++face) { + auto intersection = triangle_intersection( + roof_plane, base_cdt_ptr_->triangle(face)); + if (!intersection) continue; + auto clipped = clip_segment(*intersection, bounds); + result.insert(result.end(), clipped.begin(), clipped.end()); + } + return result; + } + private: float compute_percentile(std::vector& z_vec, float percentile) const { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1728d60e..337ab145 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -52,6 +52,12 @@ add_executable("test_logger" "${CMAKE_CURRENT_SOURCE_DIR}/test_logger.cpp") target_link_libraries("test_logger" PUBLIC logger) target_link_libraries("test_logger" PRIVATE Catch2::Catch2WithMain) +add_executable("test_arrangement_dissolver" + "${CMAKE_CURRENT_SOURCE_DIR}/test_arrangement_dissolver.cpp") +target_link_libraries("test_arrangement_dissolver" + PRIVATE roofer-core Catch2::Catch2WithMain) +catch_discover_tests("test_arrangement_dissolver") + add_executable("test_vector_reader" "${CMAKE_CURRENT_SOURCE_DIR}/test_vector_reader.cpp") target_link_libraries("test_vector_reader" PRIVATE Catch2::Catch2WithMain) diff --git a/tests/test_arrangement_dissolver.cpp b/tests/test_arrangement_dissolver.cpp new file mode 100644 index 00000000..bd6bd5f6 --- /dev/null +++ b/tests/test_arrangement_dissolver.cpp @@ -0,0 +1,93 @@ +#include + +#include + +#include +#include + +namespace { + roofer::Arrangement_2 square_arrangement() { + using roofer::Point_2; + using roofer::Segment_2; + + roofer::Arrangement_2 arrangement; + std::array points = {Point_2(0, 0), Point_2(10, 0), + Point_2(10, 10), Point_2(0, 10)}; + for (size_t i = 0; i < points.size(); ++i) { + CGAL::insert(arrangement, + Segment_2(points[i], points[(i + 1) % points.size()])); + } + + for (auto face : arrangement.face_handles()) { + if (face->is_unbounded()) continue; + face->data().is_finite = true; + face->data().in_footprint = true; + face->data().segid = 1; + } + return arrangement; + } + + roofer::RasterTools::Raster empty_heightfield() { + roofer::RasterTools::Raster raster(1, -1, 11, -1, 11); + raster.prefill_arrays(roofer::RasterTools::MIN); + return raster; + } +} // namespace + +TEST_CASE("terrain clipping marks a boundary subface as ground") { + auto arrangement = square_arrangement(); + for (auto face : arrangement.face_handles()) { + if (!face->is_unbounded()) face->data().plane = roofer::Plane(1, 0, -1, 0); + } + + auto elevation_provider = + roofer::reconstruction::createElevationProvider(5.0F); + auto dissolver = roofer::reconstruction::createArrangementDissolver(); + auto heightfield = empty_heightfield(); + dissolver->compute( + arrangement, heightfield, *elevation_provider, + {.dissolve_seg_edges = false, .dissolve_outside_fp = false}); + + size_t roof_faces = 0; + size_t ground_faces = 0; + for (auto face : arrangement.face_handles()) { + if (face->is_unbounded()) continue; + if (face->data().in_footprint) ++roof_faces; + if (face->data().is_ground) { + ++ground_faces; + CHECK_FALSE(face->data().in_footprint); + CHECK_FALSE(face->data().is_footprint_hole); + } + } + CHECK(roof_faces == 1); + CHECK(ground_faces == 1); +} + +TEST_CASE("terrain clipping marks an enclosed subface as a footprint hole") { + auto arrangement = square_arrangement(); + for (auto face : arrangement.face_handles()) { + if (!face->is_unbounded()) face->data().plane = roofer::Plane(0, 0, 1, -5); + } + + roofer::LinearRing terrain; + terrain.insert(terrain.end(), + {{0, 0, 0}, {10, 0, 0}, {10, 10, 0}, {0, 10, 0}, {5, 5, 10}}); + auto terrain_cdt = roofer::proj_tri_util::cdt_from_linearing(terrain); + auto elevation_provider = + roofer::reconstruction::createElevationProvider(terrain_cdt); + auto dissolver = roofer::reconstruction::createArrangementDissolver(); + auto heightfield = empty_heightfield(); + dissolver->compute( + arrangement, heightfield, *elevation_provider, + {.dissolve_seg_edges = false, .dissolve_outside_fp = false}); + + size_t enclosed_ground_faces = 0; + for (auto face : arrangement.face_handles()) { + if (face->is_unbounded()) continue; + if (face->data().is_ground && face->data().is_footprint_hole) { + ++enclosed_ground_faces; + CHECK_FALSE(face->data().in_footprint); + } + } + CHECK(enclosed_ground_faces == 1); +} From 464829a18f345efb539086ab275432cebd9ad191 Mon Sep 17 00:00:00 2001 From: ylannl Date: Thu, 11 Jun 2026 09:11:55 +0200 Subject: [PATCH 04/21] remove FaceInfo is_finite --- include/roofer/reconstruction/ArrangementBase.hpp | 8 -------- include/roofer/reconstruction/cgal_shared_definitions.hpp | 1 - src/core/reconstruction/ArrangementSnapper.cpp | 1 - tests/test_arrangement_dissolver.cpp | 1 - 4 files changed, 11 deletions(-) diff --git a/include/roofer/reconstruction/ArrangementBase.hpp b/include/roofer/reconstruction/ArrangementBase.hpp index 16d32939..f430047e 100644 --- a/include/roofer/reconstruction/ArrangementBase.hpp +++ b/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; diff --git a/include/roofer/reconstruction/cgal_shared_definitions.hpp b/include/roofer/reconstruction/cgal_shared_definitions.hpp index b9fe5d3c..0c2ef51b 100644 --- a/include/roofer/reconstruction/cgal_shared_definitions.hpp +++ b/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; diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index 43e6313b..b6b417ed 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -540,7 +540,6 @@ namespace roofer::reconstruction { arrFace->data() = best_face->first->data(); } else { std::cout << "Unable to locate overlapping triangle\n"; - arrFace->data().is_finite = true; arrFace->data().is_ground = true; arrFace->data().in_footprint = false; arrFace->data().is_footprint_hole = false; diff --git a/tests/test_arrangement_dissolver.cpp b/tests/test_arrangement_dissolver.cpp index bd6bd5f6..590dc04e 100644 --- a/tests/test_arrangement_dissolver.cpp +++ b/tests/test_arrangement_dissolver.cpp @@ -20,7 +20,6 @@ namespace { for (auto face : arrangement.face_handles()) { if (face->is_unbounded()) continue; - face->data().is_finite = true; face->data().in_footprint = true; face->data().segid = 1; } From e4d158bd6208641c197d010a92a69c32a012c2be Mon Sep 17 00:00:00 2001 From: ylannl Date: Fri, 12 Jun 2026 12:59:08 +0200 Subject: [PATCH 05/21] formatting --- .../reconstruction/ArrangementDissolver.cpp | 383 +++++++++--------- src/core/reconstruction/ElevationProvider.cpp | 8 +- 2 files changed, 195 insertions(+), 196 deletions(-) diff --git a/src/core/reconstruction/ArrangementDissolver.cpp b/src/core/reconstruction/ArrangementDissolver.cpp index edb32566..eb462b9d 100644 --- a/src/core/reconstruction/ArrangementDissolver.cpp +++ b/src/core/reconstruction/ArrangementDissolver.cpp @@ -135,225 +135,224 @@ namespace roofer::reconstruction { Face_data_split_observer observer(arrangement); for (const auto& intersection : intersections) { - CGAL::insert(arrangement, intersection); - } + CGAL::insert(arrangement, intersection); } + } - constexpr double elevation_tolerance = 1e-6; - std::vector> 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)); - } + constexpr double elevation_tolerance = 1e-6; + std::vector> ground_faces; + for (auto face : arrangement.face_handles()) { + if (face->is_unbounded() || !face->data().in_footprint || + face->data().segid == 0) { + continue; } - for (auto [face, touches_boundary] : ground_faces) { - face->data().in_footprint = false; - face->data().is_ground = true; - face->data().is_footprint_hole = !touches_boundary; + 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)); } } - } // 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); - } + 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); + clip_roof_faces_to_terrain(arr, elevation_provider); - if (cfg.dissolve_all_interior) { - arr_dissolve_fp(arr, true, false); - } + if (cfg.dissolve_all_interior) { + arr_dissolve_fp(arr, true, false); + } - if (cfg.dissolve_step_edges) { - // compute elevations so we can do perform the generalisation lod2.2 -> - // lod1.3 - for (auto face : arr.face_handles()) { - if (face->data().in_footprint) { - LinearRing polygon; - arrangementface_to_polygon(face, polygon); - - auto height_points = heightfield.rasterise_polygon(polygon, true); - size_t datasize = height_points.size(), data_cnt = 0; - if (datasize == - 0) { // polygon was too small to yield any pixels/height_points - auto& pz = polygon[0][2]; - face->data().elevation_50p = pz; - face->data().elevation_70p = pz; - face->data().elevation_97p = pz; - face->data().elevation_min = pz; - face->data().elevation_max = pz; - face->data().data_coverage = 0; - } else { - auto& plane = face->data().plane; - // fall back to elevation based on the plane of this face. This is - // more reliable in case the data_coverage is limited - for (auto& p : height_points) { - if (p[2] != heightfield.noDataVal_) - data_cnt++; - else - p[2] = -plane.a() / plane.c() * p[0] - - plane.b() / plane.c() * p[1] - plane.d() / plane.c(); - } - std::sort(height_points.begin(), height_points.end(), - [](auto& p1, auto& p2) { return p1[2] < p2[2]; }); - int elevation_id = std::floor(0.5 * float(datasize - 1)); - face->data().elevation_50p = height_points[elevation_id][2]; - elevation_id = std::floor(0.7 * float(datasize - 1)); - face->data().elevation_70p = height_points[elevation_id][2]; - elevation_id = std::floor(0.97 * float(datasize - 1)); - face->data().elevation_97p = height_points[elevation_id][2]; - face->data().elevation_min = height_points[0][2]; - face->data().elevation_max = height_points[datasize - 1][2]; - face->data().data_coverage = float(data_cnt) / float(datasize); + if (cfg.dissolve_step_edges) { + // compute elevations so we can do perform the generalisation lod2.2 -> + // lod1.3 + for (auto face : arr.face_handles()) { + if (face->data().in_footprint) { + LinearRing polygon; + arrangementface_to_polygon(face, polygon); + + auto height_points = heightfield.rasterise_polygon(polygon, true); + size_t datasize = height_points.size(), data_cnt = 0; + if (datasize == + 0) { // polygon was too small to yield any pixels/height_points + auto& pz = polygon[0][2]; + face->data().elevation_50p = pz; + face->data().elevation_70p = pz; + face->data().elevation_97p = pz; + face->data().elevation_min = pz; + face->data().elevation_max = pz; + face->data().data_coverage = 0; + } else { + auto& plane = face->data().plane; + // fall back to elevation based on the plane of this face. This is + // more reliable in case the data_coverage is limited + for (auto& p : height_points) { + if (p[2] != heightfield.noDataVal_) + data_cnt++; + else + p[2] = -plane.a() / plane.c() * p[0] - + plane.b() / plane.c() * p[1] - plane.d() / plane.c(); } - face->data().pixel_count = datasize; + std::sort(height_points.begin(), height_points.end(), + [](auto& p1, auto& p2) { return p1[2] < p2[2]; }); + int elevation_id = std::floor(0.5 * float(datasize - 1)); + face->data().elevation_50p = height_points[elevation_id][2]; + elevation_id = std::floor(0.7 * float(datasize - 1)); + face->data().elevation_70p = height_points[elevation_id][2]; + elevation_id = std::floor(0.97 * float(datasize - 1)); + face->data().elevation_97p = height_points[elevation_id][2]; + face->data().elevation_min = height_points[0][2]; + face->data().elevation_max = height_points[datasize - 1][2]; + face->data().data_coverage = float(data_cnt) / float(datasize); } + face->data().pixel_count = datasize; } - - Face_merge_observer obs(arr); - arr_dissolve_step_edges(arr, cfg.step_height_threshold); } - // remove any dangling edges - { - std::vector to_remove; - for (auto he : arr.edge_handles()) { - if (he->face() == he->twin()->face()) to_remove.push_back(he); - } - for (auto he : to_remove) { - arr.remove_edge(he); - } - } + Face_merge_observer obs(arr); + arr_dissolve_step_edges(arr, cfg.step_height_threshold); + } - // remove all lines not inside footprint - if (cfg.dissolve_outside_fp) { - arr_dissolve_fp(arr, false, true); + // remove any dangling edges + { + std::vector to_remove; + for (auto he : arr.edge_handles()) { + if (he->face() == he->twin()->face()) to_remove.push_back(he); + } + for (auto he : to_remove) { + arr.remove_edge(he); } + } - arr_remove_redundant_vertices(arr); + // remove all lines not inside footprint + if (cfg.dissolve_outside_fp) { + arr_dissolve_fp(arr, false, true); + } - // label different building parts - arr_label_buildingparts(arr); + arr_remove_redundant_vertices(arr); - // ensure all holes are marked as such - auto f_unb = arr.unbounded_face(); - for (auto fh : arr.face_handles()) { - if (fh != f_unb) - if (!fh->data().in_footprint && !fh->data().is_footprint_hole) { - if (!fh->data().is_ground || !touches_unbounded_face(fh)) { - fh->data().is_footprint_hole = true; - } - } - } + // label different building parts + arr_label_buildingparts(arr); - // if (remove_duplicates) { - // arr_snap_duplicates(arr, (double) std::pow(10,-dupe_threshold_exp)); - // } - // snap edge shorter than snap_tolerance - // for (auto he : arr.edge_handles()) { - // if(CGAL::squared_distance(he->source()->point(), - // he->target()->point()) < snap_tolerance) { - - // arr.remove_edge(he->next()); - // arr.remove_edge(he->twin()->previous()); - // } - // } - // compute final data_area and elevation stats for each face - std::vector all_elevations; - for (auto face : arr.face_handles()) { - if (face->data().in_footprint) { - LinearRing polygon; - arrangementface_to_polygon(face, polygon); + // ensure all holes are marked as such + auto f_unb = arr.unbounded_face(); + for (auto fh : arr.face_handles()) { + if (fh != f_unb) + if (!fh->data().in_footprint && !fh->data().is_footprint_hole) { + if (!fh->data().is_ground || !touches_unbounded_face(fh)) { + fh->data().is_footprint_hole = true; + } + } + } - auto height_points = heightfield.rasterise_polygon(polygon, true); - size_t data_cnt = 0; - bool skip_face = height_points.size() == 0; - { - // auto& plane = face->data().plane; - // compute elevations based on the assigned plane on this face - std::vector pts; - if (!skip_face) { - for (auto& p : height_points) { - if (p[2] != heightfield.noDataVal_) { - data_cnt++; - pts.push_back(&p); - all_elevations.push_back(p[2]); - }; - } - } - if (data_cnt == 0) { - auto& pz = polygon[0][2]; - face->data().elevation_50p = heightfield.noDataVal_; - face->data().elevation_70p = heightfield.noDataVal_; - face->data().elevation_97p = heightfield.noDataVal_; - face->data().elevation_min = heightfield.noDataVal_; - face->data().elevation_max = heightfield.noDataVal_; - face->data().data_coverage = 0; - } else { - std::sort(pts.begin(), pts.end(), - [](auto& p1, auto& p2) { return (*p1)[2] < (*p2)[2]; }); - size_t datasize = pts.size(); - int elevation_id = std::floor(0.5 * float(datasize - 1)); - face->data().elevation_50p = (*(pts[elevation_id]))[2]; - elevation_id = std::floor(0.7 * float(datasize - 1)); - face->data().elevation_70p = (*(pts[elevation_id]))[2]; - elevation_id = std::floor(0.97 * float(datasize - 1)); - face->data().elevation_97p = (*(pts[elevation_id]))[2]; - face->data().elevation_min = (*(pts[0]))[2]; - face->data().elevation_max = (*(pts[datasize - 1]))[2]; - face->data().data_coverage = - float(data_cnt) / float(height_points.size()); + // if (remove_duplicates) { + // arr_snap_duplicates(arr, (double) std::pow(10,-dupe_threshold_exp)); + // } + // snap edge shorter than snap_tolerance + // for (auto he : arr.edge_handles()) { + // if(CGAL::squared_distance(he->source()->point(), + // he->target()->point()) < snap_tolerance) { + + // arr.remove_edge(he->next()); + // arr.remove_edge(he->twin()->previous()); + // } + // } + // compute final data_area and elevation stats for each face + std::vector all_elevations; + for (auto face : arr.face_handles()) { + if (face->data().in_footprint) { + LinearRing polygon; + arrangementface_to_polygon(face, polygon); + + auto height_points = heightfield.rasterise_polygon(polygon, true); + size_t data_cnt = 0; + bool skip_face = height_points.size() == 0; + { + // auto& plane = face->data().plane; + // compute elevations based on the assigned plane on this face + std::vector pts; + if (!skip_face) { + for (auto& p : height_points) { + if (p[2] != heightfield.noDataVal_) { + data_cnt++; + pts.push_back(&p); + all_elevations.push_back(p[2]); + }; } - face->data().pixel_count = data_cnt; } + if (data_cnt == 0) { + auto& pz = polygon[0][2]; + face->data().elevation_50p = heightfield.noDataVal_; + face->data().elevation_70p = heightfield.noDataVal_; + face->data().elevation_97p = heightfield.noDataVal_; + face->data().elevation_min = heightfield.noDataVal_; + face->data().elevation_max = heightfield.noDataVal_; + face->data().data_coverage = 0; + } else { + std::sort(pts.begin(), pts.end(), + [](auto& p1, auto& p2) { return (*p1)[2] < (*p2)[2]; }); + size_t datasize = pts.size(); + int elevation_id = std::floor(0.5 * float(datasize - 1)); + face->data().elevation_50p = (*(pts[elevation_id]))[2]; + elevation_id = std::floor(0.7 * float(datasize - 1)); + face->data().elevation_70p = (*(pts[elevation_id]))[2]; + elevation_id = std::floor(0.97 * float(datasize - 1)); + face->data().elevation_97p = (*(pts[elevation_id]))[2]; + face->data().elevation_min = (*(pts[0]))[2]; + face->data().elevation_max = (*(pts[datasize - 1]))[2]; + face->data().data_coverage = + float(data_cnt) / float(height_points.size()); + } + face->data().pixel_count = data_cnt; } } - - // compute golbal height statistics - // if (all_elevations.size()==0) { - // output("global_elevation_50p").set_from_any(std::any()); - // output("global_elevation_70p").set_from_any(std::any()); - // output("global_elevation_97p").set_from_any(std::any()); - // output("global_elevation_min").set_from_any(std::any()); - // output("global_elevation_max").set_from_any(std::any()); - // } else { - // std::sort(all_elevations.begin(), all_elevations.end()); - // size_t last_id = all_elevations.size()-1; - // int elevation_id = std::floor(0.5*float(last_id)); - // output("global_elevation_50p").set(all_elevations[elevation_id]); - // elevation_id = std::floor(0.7*float(last_id)); - // output("global_elevation_70p").set(all_elevations[elevation_id]); - // elevation_id = std::floor(0.99*float(last_id)); - // output("global_elevation_97p").set(all_elevations[elevation_id]); - // output("global_elevation_min").set(all_elevations[0]); - // output("global_elevation_max").set(all_elevations[last_id]); - // } - // 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); - } - }; + // compute golbal height statistics + // if (all_elevations.size()==0) { + // output("global_elevation_50p").set_from_any(std::any()); + // output("global_elevation_70p").set_from_any(std::any()); + // output("global_elevation_97p").set_from_any(std::any()); + // output("global_elevation_min").set_from_any(std::any()); + // output("global_elevation_max").set_from_any(std::any()); + // } else { + // std::sort(all_elevations.begin(), all_elevations.end()); + // size_t last_id = all_elevations.size()-1; + // int elevation_id = std::floor(0.5*float(last_id)); + // output("global_elevation_50p").set(all_elevations[elevation_id]); + // elevation_id = std::floor(0.7*float(last_id)); + // output("global_elevation_70p").set(all_elevations[elevation_id]); + // elevation_id = std::floor(0.99*float(last_id)); + // output("global_elevation_97p").set(all_elevations[elevation_id]); + // output("global_elevation_min").set(all_elevations[0]); + // output("global_elevation_max").set(all_elevations[last_id]); + // } + // output("arrangement").set(arr); + } - std::unique_ptr createArrangementDissolver() { - return std::make_unique(); + 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 createArrangementDissolver() { + return std::make_unique(); +} } // namespace roofer::reconstruction diff --git a/src/core/reconstruction/ElevationProvider.cpp b/src/core/reconstruction/ElevationProvider.cpp index 59f0577f..ac4934db 100644 --- a/src/core/reconstruction/ElevationProvider.cpp +++ b/src/core/reconstruction/ElevationProvider.cpp @@ -134,7 +134,7 @@ namespace roofer::reconstruction { const float floor_elevation_; ConstantElevationProvider(const float floor_elevation) - : floor_elevation_(floor_elevation) {}; + : floor_elevation_(floor_elevation){}; virtual float get(const Point_2 /* pt */) const override { return floor_elevation_; @@ -164,7 +164,7 @@ namespace roofer::reconstruction { std::shared_ptr base_cdt_ptr_; InterpolatedElevationProvider(const proj_tri_util::DT& base_cdt) - : base_cdt_ptr_(std::make_shared(base_cdt)) {}; + : base_cdt_ptr_(std::make_shared(base_cdt)){}; virtual float get(const Point_2 pt) const override { return proj_tri_util::interpolate_from_cdt(pt, *base_cdt_ptr_); @@ -190,8 +190,8 @@ namespace roofer::reconstruction { std::vector result; for (auto face = base_cdt_ptr_->finite_faces_begin(); face != base_cdt_ptr_->finite_faces_end(); ++face) { - auto intersection = triangle_intersection( - roof_plane, base_cdt_ptr_->triangle(face)); + auto intersection = + triangle_intersection(roof_plane, base_cdt_ptr_->triangle(face)); if (!intersection) continue; auto clipped = clip_segment(*intersection, bounds); result.insert(result.end(), clipped.begin(), clipped.end()); From 8853a656960c0c552bf3733e08e0eb906c604d9e Mon Sep 17 00:00:00 2001 From: ylannl Date: Fri, 12 Jun 2026 13:14:00 +0200 Subject: [PATCH 06/21] fix syntax --- .../reconstruction/ArrangementDissolver.cpp | 380 +++++++++--------- 1 file changed, 190 insertions(+), 190 deletions(-) diff --git a/src/core/reconstruction/ArrangementDissolver.cpp b/src/core/reconstruction/ArrangementDissolver.cpp index eb462b9d..b73cae1f 100644 --- a/src/core/reconstruction/ArrangementDissolver.cpp +++ b/src/core/reconstruction/ArrangementDissolver.cpp @@ -137,222 +137,222 @@ namespace roofer::reconstruction { for (const auto& intersection : intersections) { CGAL::insert(arrangement, intersection); } - } - constexpr double elevation_tolerance = 1e-6; - std::vector> ground_faces; - for (auto face : arrangement.face_handles()) { - if (face->is_unbounded() || !face->data().in_footprint || - face->data().segid == 0) { - continue; + constexpr double elevation_tolerance = 1e-6; + std::vector> 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)); + } } - 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; } } - 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); + } // 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); + } - if (cfg.dissolve_all_interior) { - arr_dissolve_fp(arr, true, false); - } + clip_roof_faces_to_terrain(arr, elevation_provider); - if (cfg.dissolve_step_edges) { - // compute elevations so we can do perform the generalisation lod2.2 -> - // lod1.3 - for (auto face : arr.face_handles()) { - if (face->data().in_footprint) { - LinearRing polygon; - arrangementface_to_polygon(face, polygon); + if (cfg.dissolve_all_interior) { + arr_dissolve_fp(arr, true, false); + } - auto height_points = heightfield.rasterise_polygon(polygon, true); - size_t datasize = height_points.size(), data_cnt = 0; - if (datasize == - 0) { // polygon was too small to yield any pixels/height_points - auto& pz = polygon[0][2]; - face->data().elevation_50p = pz; - face->data().elevation_70p = pz; - face->data().elevation_97p = pz; - face->data().elevation_min = pz; - face->data().elevation_max = pz; - face->data().data_coverage = 0; - } else { - auto& plane = face->data().plane; - // fall back to elevation based on the plane of this face. This is - // more reliable in case the data_coverage is limited - for (auto& p : height_points) { - if (p[2] != heightfield.noDataVal_) - data_cnt++; - else - p[2] = -plane.a() / plane.c() * p[0] - - plane.b() / plane.c() * p[1] - plane.d() / plane.c(); + if (cfg.dissolve_step_edges) { + // compute elevations so we can do perform the generalisation lod2.2 -> + // lod1.3 + for (auto face : arr.face_handles()) { + if (face->data().in_footprint) { + LinearRing polygon; + arrangementface_to_polygon(face, polygon); + + auto height_points = heightfield.rasterise_polygon(polygon, true); + size_t datasize = height_points.size(), data_cnt = 0; + if (datasize == + 0) { // polygon was too small to yield any pixels/height_points + auto& pz = polygon[0][2]; + face->data().elevation_50p = pz; + face->data().elevation_70p = pz; + face->data().elevation_97p = pz; + face->data().elevation_min = pz; + face->data().elevation_max = pz; + face->data().data_coverage = 0; + } else { + auto& plane = face->data().plane; + // fall back to elevation based on the plane of this face. This is + // more reliable in case the data_coverage is limited + for (auto& p : height_points) { + if (p[2] != heightfield.noDataVal_) + data_cnt++; + else + p[2] = -plane.a() / plane.c() * p[0] - + plane.b() / plane.c() * p[1] - plane.d() / plane.c(); + } + std::sort(height_points.begin(), height_points.end(), + [](auto& p1, auto& p2) { return p1[2] < p2[2]; }); + int elevation_id = std::floor(0.5 * float(datasize - 1)); + face->data().elevation_50p = height_points[elevation_id][2]; + elevation_id = std::floor(0.7 * float(datasize - 1)); + face->data().elevation_70p = height_points[elevation_id][2]; + elevation_id = std::floor(0.97 * float(datasize - 1)); + face->data().elevation_97p = height_points[elevation_id][2]; + face->data().elevation_min = height_points[0][2]; + face->data().elevation_max = height_points[datasize - 1][2]; + face->data().data_coverage = float(data_cnt) / float(datasize); } - std::sort(height_points.begin(), height_points.end(), - [](auto& p1, auto& p2) { return p1[2] < p2[2]; }); - int elevation_id = std::floor(0.5 * float(datasize - 1)); - face->data().elevation_50p = height_points[elevation_id][2]; - elevation_id = std::floor(0.7 * float(datasize - 1)); - face->data().elevation_70p = height_points[elevation_id][2]; - elevation_id = std::floor(0.97 * float(datasize - 1)); - face->data().elevation_97p = height_points[elevation_id][2]; - face->data().elevation_min = height_points[0][2]; - face->data().elevation_max = height_points[datasize - 1][2]; - face->data().data_coverage = float(data_cnt) / float(datasize); + face->data().pixel_count = datasize; } - face->data().pixel_count = datasize; } - } - - Face_merge_observer obs(arr); - arr_dissolve_step_edges(arr, cfg.step_height_threshold); - } - // remove any dangling edges - { - std::vector to_remove; - for (auto he : arr.edge_handles()) { - if (he->face() == he->twin()->face()) to_remove.push_back(he); + Face_merge_observer obs(arr); + arr_dissolve_step_edges(arr, cfg.step_height_threshold); } - for (auto he : to_remove) { - arr.remove_edge(he); + + // remove any dangling edges + { + std::vector to_remove; + for (auto he : arr.edge_handles()) { + if (he->face() == he->twin()->face()) to_remove.push_back(he); + } + for (auto he : to_remove) { + arr.remove_edge(he); + } } - } - // remove all lines not inside footprint - if (cfg.dissolve_outside_fp) { - arr_dissolve_fp(arr, false, true); - } + // remove all lines not inside footprint + if (cfg.dissolve_outside_fp) { + arr_dissolve_fp(arr, false, true); + } - arr_remove_redundant_vertices(arr); + arr_remove_redundant_vertices(arr); - // label different building parts - arr_label_buildingparts(arr); + // label different building parts + arr_label_buildingparts(arr); - // ensure all holes are marked as such - auto f_unb = arr.unbounded_face(); - for (auto fh : arr.face_handles()) { - if (fh != f_unb) - if (!fh->data().in_footprint && !fh->data().is_footprint_hole) { - if (!fh->data().is_ground || !touches_unbounded_face(fh)) { - fh->data().is_footprint_hole = true; + // ensure all holes are marked as such + auto f_unb = arr.unbounded_face(); + for (auto fh : arr.face_handles()) { + if (fh != f_unb) + if (!fh->data().in_footprint && !fh->data().is_footprint_hole) { + if (!fh->data().is_ground || !touches_unbounded_face(fh)) { + fh->data().is_footprint_hole = true; + } } - } - } + } + + // if (remove_duplicates) { + // arr_snap_duplicates(arr, (double) std::pow(10,-dupe_threshold_exp)); + // } + // snap edge shorter than snap_tolerance + // for (auto he : arr.edge_handles()) { + // if(CGAL::squared_distance(he->source()->point(), + // he->target()->point()) < snap_tolerance) { + + // arr.remove_edge(he->next()); + // arr.remove_edge(he->twin()->previous()); + // } + // } + // compute final data_area and elevation stats for each face + std::vector all_elevations; + for (auto face : arr.face_handles()) { + if (face->data().in_footprint) { + LinearRing polygon; + arrangementface_to_polygon(face, polygon); - // if (remove_duplicates) { - // arr_snap_duplicates(arr, (double) std::pow(10,-dupe_threshold_exp)); - // } - // snap edge shorter than snap_tolerance - // for (auto he : arr.edge_handles()) { - // if(CGAL::squared_distance(he->source()->point(), - // he->target()->point()) < snap_tolerance) { - - // arr.remove_edge(he->next()); - // arr.remove_edge(he->twin()->previous()); - // } - // } - // compute final data_area and elevation stats for each face - std::vector all_elevations; - for (auto face : arr.face_handles()) { - if (face->data().in_footprint) { - LinearRing polygon; - arrangementface_to_polygon(face, polygon); - - auto height_points = heightfield.rasterise_polygon(polygon, true); - size_t data_cnt = 0; - bool skip_face = height_points.size() == 0; - { - // auto& plane = face->data().plane; - // compute elevations based on the assigned plane on this face - std::vector pts; - if (!skip_face) { - for (auto& p : height_points) { - if (p[2] != heightfield.noDataVal_) { - data_cnt++; - pts.push_back(&p); - all_elevations.push_back(p[2]); - }; + auto height_points = heightfield.rasterise_polygon(polygon, true); + size_t data_cnt = 0; + bool skip_face = height_points.size() == 0; + { + // auto& plane = face->data().plane; + // compute elevations based on the assigned plane on this face + std::vector pts; + if (!skip_face) { + for (auto& p : height_points) { + if (p[2] != heightfield.noDataVal_) { + data_cnt++; + pts.push_back(&p); + all_elevations.push_back(p[2]); + }; + } } + if (data_cnt == 0) { + auto& pz = polygon[0][2]; + face->data().elevation_50p = heightfield.noDataVal_; + face->data().elevation_70p = heightfield.noDataVal_; + face->data().elevation_97p = heightfield.noDataVal_; + face->data().elevation_min = heightfield.noDataVal_; + face->data().elevation_max = heightfield.noDataVal_; + face->data().data_coverage = 0; + } else { + std::sort(pts.begin(), pts.end(), + [](auto& p1, auto& p2) { return (*p1)[2] < (*p2)[2]; }); + size_t datasize = pts.size(); + int elevation_id = std::floor(0.5 * float(datasize - 1)); + face->data().elevation_50p = (*(pts[elevation_id]))[2]; + elevation_id = std::floor(0.7 * float(datasize - 1)); + face->data().elevation_70p = (*(pts[elevation_id]))[2]; + elevation_id = std::floor(0.97 * float(datasize - 1)); + face->data().elevation_97p = (*(pts[elevation_id]))[2]; + face->data().elevation_min = (*(pts[0]))[2]; + face->data().elevation_max = (*(pts[datasize - 1]))[2]; + face->data().data_coverage = + float(data_cnt) / float(height_points.size()); + } + face->data().pixel_count = data_cnt; } - if (data_cnt == 0) { - auto& pz = polygon[0][2]; - face->data().elevation_50p = heightfield.noDataVal_; - face->data().elevation_70p = heightfield.noDataVal_; - face->data().elevation_97p = heightfield.noDataVal_; - face->data().elevation_min = heightfield.noDataVal_; - face->data().elevation_max = heightfield.noDataVal_; - face->data().data_coverage = 0; - } else { - std::sort(pts.begin(), pts.end(), - [](auto& p1, auto& p2) { return (*p1)[2] < (*p2)[2]; }); - size_t datasize = pts.size(); - int elevation_id = std::floor(0.5 * float(datasize - 1)); - face->data().elevation_50p = (*(pts[elevation_id]))[2]; - elevation_id = std::floor(0.7 * float(datasize - 1)); - face->data().elevation_70p = (*(pts[elevation_id]))[2]; - elevation_id = std::floor(0.97 * float(datasize - 1)); - face->data().elevation_97p = (*(pts[elevation_id]))[2]; - face->data().elevation_min = (*(pts[0]))[2]; - face->data().elevation_max = (*(pts[datasize - 1]))[2]; - face->data().data_coverage = - float(data_cnt) / float(height_points.size()); - } - face->data().pixel_count = data_cnt; } } + + // compute golbal height statistics + // if (all_elevations.size()==0) { + // output("global_elevation_50p").set_from_any(std::any()); + // output("global_elevation_70p").set_from_any(std::any()); + // output("global_elevation_97p").set_from_any(std::any()); + // output("global_elevation_min").set_from_any(std::any()); + // output("global_elevation_max").set_from_any(std::any()); + // } else { + // std::sort(all_elevations.begin(), all_elevations.end()); + // size_t last_id = all_elevations.size()-1; + // int elevation_id = std::floor(0.5*float(last_id)); + // output("global_elevation_50p").set(all_elevations[elevation_id]); + // elevation_id = std::floor(0.7*float(last_id)); + // output("global_elevation_70p").set(all_elevations[elevation_id]); + // elevation_id = std::floor(0.99*float(last_id)); + // output("global_elevation_97p").set(all_elevations[elevation_id]); + // output("global_elevation_min").set(all_elevations[0]); + // output("global_elevation_max").set(all_elevations[last_id]); + // } + // output("arrangement").set(arr); } - // compute golbal height statistics - // if (all_elevations.size()==0) { - // output("global_elevation_50p").set_from_any(std::any()); - // output("global_elevation_70p").set_from_any(std::any()); - // output("global_elevation_97p").set_from_any(std::any()); - // output("global_elevation_min").set_from_any(std::any()); - // output("global_elevation_max").set_from_any(std::any()); - // } else { - // std::sort(all_elevations.begin(), all_elevations.end()); - // size_t last_id = all_elevations.size()-1; - // int elevation_id = std::floor(0.5*float(last_id)); - // output("global_elevation_50p").set(all_elevations[elevation_id]); - // elevation_id = std::floor(0.7*float(last_id)); - // output("global_elevation_70p").set(all_elevations[elevation_id]); - // elevation_id = std::floor(0.99*float(last_id)); - // output("global_elevation_97p").set(all_elevations[elevation_id]); - // output("global_elevation_min").set(all_elevations[0]); - // output("global_elevation_max").set(all_elevations[last_id]); - // } - // 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); + } + }; - 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 createArrangementDissolver() { + return std::make_unique(); } -}; - -std::unique_ptr createArrangementDissolver() { - return std::make_unique(); -} } // namespace roofer::reconstruction From ea069e70db2a0945b4cb345b0f69c056b4d7c037 Mon Sep 17 00:00:00 2001 From: ylannl Date: Fri, 12 Jun 2026 13:49:43 +0200 Subject: [PATCH 07/21] fix linking issue with arrangement bowtie test --- tests/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 337ab145..47d139aa 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -55,7 +55,7 @@ target_link_libraries("test_logger" PRIVATE Catch2::Catch2WithMain) add_executable("test_arrangement_dissolver" "${CMAKE_CURRENT_SOURCE_DIR}/test_arrangement_dissolver.cpp") target_link_libraries("test_arrangement_dissolver" - PRIVATE roofer-core Catch2::Catch2WithMain) + PRIVATE Catch2::Catch2WithMain roofer-core) catch_discover_tests("test_arrangement_dissolver") add_executable("test_vector_reader" From 371549b4561dd5c9bfbbfcceafa4e1b0d62d6f63 Mon Sep 17 00:00:00 2001 From: ylannl Date: Mon, 15 Jun 2026 22:20:50 +0200 Subject: [PATCH 08/21] Prevent non manifold edges --- .../reconstruction/ArrangementSnapper.hpp | 3 + .../reconstruction/ArrangementSnapper.cpp | 396 +++++++++++++++--- tests/CMakeLists.txt | 6 + tests/test_arrangement_snapper.cpp | 129 ++++++ 4 files changed, 487 insertions(+), 47 deletions(-) create mode 100644 tests/test_arrangement_snapper.cpp diff --git a/include/roofer/reconstruction/ArrangementSnapper.hpp b/include/roofer/reconstruction/ArrangementSnapper.hpp index ee050b84..a4d9aeda 100644 --- a/include/roofer/reconstruction/ArrangementSnapper.hpp +++ b/include/roofer/reconstruction/ArrangementSnapper.hpp @@ -28,6 +28,9 @@ namespace roofer::reconstruction { struct ArrangementSnapperConfig { float dist_thres = 0.05; + bool repair_non_manifold_vertices = true; + float manifold_repair_radius = 0.05; + float manifold_height_tolerance = 1e-4; }; struct ArrangementSnapperInterface { diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index b6b417ed..851d0ba6 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -30,6 +30,12 @@ #include #include +#include +#include +#include +#include +#include + namespace roofer::reconstruction { namespace arragementsnapper { @@ -60,6 +66,296 @@ namespace roofer::reconstruction { typedef std::unordered_map ConstraintMap; + struct ForcedRegionLabel { + T::Point_2 seed; + FaceInfo* face_info; + }; + + using SourceFaceIds = std::unordered_map; + + struct IncidentSector { + Vertex_handle neighbour; + FaceInfo* face_info; + double angle; + double height; + }; + + double plane_height(const FaceInfo& face_info, const T::Point_2& point) { + const auto& plane = face_info.plane; + if (plane.c() == 0) return 0; + return -(plane.a() * point.x() + plane.b() * point.y() + plane.d()) / + plane.c(); + } + + FaceInfo* source_face_at( + const T::Point_2& point, + CGAL::Arr_walk_along_line_point_location& walk_pl, + Arrangement_2& source_arrangement) { + auto object = + walk_pl.locate(Arrangement_2::Point_2(point.x(), point.y())); + if (auto face = std::get_if(&object)) { + return &source_arrangement.non_const_handle(*face)->data(); + } + return nullptr; + } + + std::vector> constrained_regions(T& tri) { + std::vector> regions; + std::unordered_set visited; + + for (auto face : tri.finite_face_handles()) { + if (visited.contains(face)) continue; + + std::vector region; + std::queue queue; + queue.push(face); + visited.insert(face); + while (!queue.empty()) { + auto current = queue.front(); + queue.pop(); + region.push_back(current); + + for (int i = 0; i < 3; ++i) { + Edge edge(current, i); + auto neighbour = current->neighbor(i); + if (!tri.is_constrained(edge) && !tri.is_infinite(neighbour) && + visited.insert(neighbour).second) { + queue.push(neighbour); + } + } + } + regions.push_back(std::move(region)); + } + return regions; + } + + void label_final_regions( + T& tri, + CGAL::Arr_walk_along_line_point_location& walk_pl, + Arrangement_2& source_arrangement, const SourceFaceIds& source_ids, + const std::vector& forced_labels) { + for (auto& region : constrained_regions(tri)) { + std::unordered_map overlap_area; + for (auto face : region) { + auto centroid = CGAL::centroid(tri.triangle(face)); + if (auto* source = + source_face_at(centroid, walk_pl, source_arrangement)) { + overlap_area[source] += std::abs(tri.triangle(face).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; + } + } + + for (const auto& forced : forced_labels) { + auto located = tri.locate(forced.seed); + if (!tri.is_infinite(located) && + std::find(region.begin(), region.end(), located) != + region.end()) { + selected = forced.face_info; + break; + } + } + + for (auto face : region) face->info() = selected; + } + } + + std::optional vertex_clearance(T& tri, Vertex_handle vertex) { + double clearance_sq = std::numeric_limits::max(); + Face_circulator face = tri.incident_faces(vertex), done(face); + if (face == nullptr) return std::nullopt; + do { + if (tri.is_infinite(face)) return std::nullopt; + auto opposite_edge = Edge(face, face->index(vertex)); + clearance_sq = std::min( + clearance_sq, CGAL::to_double(CGAL::squared_distance( + vertex->point(), tri.segment(opposite_edge)))); + } while (++face != done); + if (!(clearance_sq > 0)) return std::nullopt; + return std::sqrt(clearance_sq); + } + + std::vector incident_sectors(T& tri, Vertex_handle vertex) { + std::vector sectors; + Edge_circulator edge = tri.incident_edges(vertex), done(edge); + if (edge == nullptr) return sectors; + do { + if (!tri.is_constrained(*edge)) continue; + auto first = edge->first->vertex(tri.cw(edge->second)); + auto second = edge->first->vertex(tri.ccw(edge->second)); + auto neighbour = first == vertex ? second : first; + const double dx = neighbour->point().x() - vertex->point().x(); + const double dy = neighbour->point().y() - vertex->point().y(); + sectors.push_back({neighbour, nullptr, std::atan2(dy, dx), 0}); + } while (++edge != done); + + std::sort(sectors.begin(), sectors.end(), + [](const auto& lhs, const auto& rhs) { + return lhs.angle < rhs.angle; + }); + if (sectors.size() < 2) return sectors; + + auto clearance = vertex_clearance(tri, vertex); + if (!clearance) return {}; + const double sample_radius = *clearance * 0.25; + constexpr double two_pi = 2 * CGAL_PI; + for (std::size_t i = 0; i < sectors.size(); ++i) { + double next_angle = sectors[(i + 1) % sectors.size()].angle; + if (next_angle <= sectors[i].angle) next_angle += two_pi; + const double sample_angle = (sectors[i].angle + next_angle) * 0.5; + T::Point_2 sample( + vertex->point().x() + sample_radius * std::cos(sample_angle), + vertex->point().y() + sample_radius * std::sin(sample_angle)); + auto face = tri.locate(sample); + if (tri.is_infinite(face) || face->info() == nullptr) return {}; + sectors[i].face_info = face->info(); + sectors[i].height = plane_height(*face->info(), vertex->point()); + } + return sectors; + } + + bool has_non_manifold_height_order( + const std::vector& sectors, double tolerance) { + std::vector heights; + heights.reserve(sectors.size()); + for (const auto& sector : sectors) heights.push_back(sector.height); + std::sort(heights.begin(), heights.end()); + + std::vector distinct_heights; + for (double height : heights) { + if (distinct_heights.empty() || + height - distinct_heights.back() >= tolerance) { + distinct_heights.push_back(height); + } + } + + for (std::size_t level = 1; level < distinct_heights.size(); ++level) { + const double sample_height = + (distinct_heights[level - 1] + distinct_heights[level]) * 0.5; + std::size_t crossings = 0; + for (std::size_t i = 0; i < sectors.size(); ++i) { + const double lhs = + sectors[(i + sectors.size() - 1) % sectors.size()].height; + const double rhs = sectors[i].height; + if ((lhs < sample_height && rhs > sample_height) || + (rhs < sample_height && lhs > sample_height)) { + ++crossings; + } + } + if (crossings > 2) return true; + } + return false; + } + + std::optional find_problematic_vertex( + T& tri, double height_tolerance) { + for (auto vertex = tri.finite_vertices_begin(); + vertex != tri.finite_vertices_end(); ++vertex) { + if (vertex->info()) continue; + auto sectors = incident_sectors(tri, vertex); + if (sectors.size() < 4) continue; + if (std::any_of(sectors.begin(), sectors.end(), [](const auto& sector) { + return !sector.face_info->in_footprint || + sector.face_info->segid == 0; + })) { + continue; + } + if (has_non_manifold_height_order(sectors, height_tolerance)) { + return vertex; + } + } + return std::nullopt; + } + + ForcedRegionLabel repair_vertex(T& tri, Vertex_handle vertex, + const SourceFaceIds& source_ids, + double maximum_radius, + double height_tolerance) { + auto sectors = incident_sectors(tri, vertex); + auto clearance = vertex_clearance(tri, vertex); + if (sectors.size() < 4 || !clearance) { + throw roofer::rooferException( + "Unable to compute a safe non-manifold junction repair"); + } + + auto selected = std::min_element( + sectors.begin(), sectors.end(), + [&](const auto& lhs, const auto& rhs) { + if (std::abs(lhs.height - rhs.height) >= height_tolerance) + return lhs.height < rhs.height; + return source_ids.at(lhs.face_info) < source_ids.at(rhs.face_info); + }); + FaceInfo* selected_face = selected->face_info; + + double shortest_edge = std::numeric_limits::max(); + for (const auto& sector : sectors) { + shortest_edge = std::min( + shortest_edge, std::sqrt(CGAL::to_double(CGAL::squared_distance( + vertex->point(), sector.neighbour->point())))); + } + const double radius = + std::min({maximum_radius, *clearance * 0.8, shortest_edge * 0.45}); + if (!(radius > std::numeric_limits::epsilon())) { + throw roofer::rooferException( + "Non-manifold junction repair radius is degenerate"); + } + + const auto original_point = vertex->point(); + std::vector split_vertices; + split_vertices.reserve(sectors.size()); + for (const auto& sector : sectors) { + const double dx = sector.neighbour->point().x() - original_point.x(); + const double dy = sector.neighbour->point().y() - original_point.y(); + const double length = std::sqrt(dx * dx + dy * dy); + T::Point_2 split_point(original_point.x() + radius * dx / length, + original_point.y() + radius * dy / length); + auto split_vertex = tri.insert(split_point); + split_vertex->info() = false; + split_vertices.push_back(split_vertex); + } + + tri.remove_incident_constraints(vertex); + tri.remove(vertex); + for (std::size_t i = 0; i < sectors.size(); ++i) { + tri.insert_constraint(split_vertices[i], sectors[i].neighbour); + } + + std::optional forced_seed; + for (std::size_t i = 0; i < sectors.size(); ++i) { + auto next = (i + 1) % sectors.size(); + if (sectors[i].face_info == selected_face) { + if (!forced_seed) { + forced_seed = T::Point_2( + (2 * original_point.x() + split_vertices[i]->point().x() + + split_vertices[next]->point().x()) / + 4, + (2 * original_point.y() + split_vertices[i]->point().y() + + split_vertices[next]->point().y()) / + 4); + } + continue; + } + tri.insert_constraint(split_vertices[i], split_vertices[next]); + } + + if (!forced_seed) { + throw roofer::rooferException( + "Non-manifold junction repair has no merge sector"); + } + return {*forced_seed, selected_face}; + } + // tri_util::CDT triangulate_polygon(LinearRing& poly, float // dupe_threshold_exp=3) { // tri_util::CDT triangulation; @@ -163,6 +459,12 @@ namespace roofer::reconstruction { T tri; float sq_dist_thres = cfg.dist_thres * cfg.dist_thres; + SourceFaceIds source_face_ids; + std::size_t next_source_face_id = 0; + for (auto face : arr.face_handles()) { + source_face_ids[&face->data()] = next_source_face_id++; + } + // map from arr vertices to tri vertices std::unordered_map vertex_map; @@ -443,6 +745,20 @@ namespace roofer::reconstruction { } // } while (found_short_edge); + std::vector forced_region_labels; + label_final_regions(tri, walk_pl, arr, source_face_ids, + forced_region_labels); + if (cfg.repair_non_manifold_vertices) { + while (auto vertex = find_problematic_vertex( + tri, cfg.manifold_height_tolerance)) { + forced_region_labels.push_back(repair_vertex( + tri, *vertex, source_face_ids, cfg.manifold_repair_radius, + cfg.manifold_height_tolerance)); + label_final_regions(tri, walk_pl, arr, source_face_ids, + forced_region_labels); + } + } + // TriangleCollection triangles_snapped; // vec1i segment_ids_snapped; // for (auto fh = tri.finite_faces_begin(); fh != @@ -496,55 +812,41 @@ namespace roofer::reconstruction { // } } - for (auto& arrFace : arr_snap.face_handles()) { - LinearRing poly; - if (arrFace->is_fictitious()) continue; - - if (arrangementface_to_polygon(arrFace, poly)) { - // std::cout << "poly size: " << poly.size() << std::endl; - tri_util::CDT cdt = tri_util::create_from_polygon(poly); - - std::unordered_map - canidate_faces; - for (tri_util::CDT::Finite_faces_iterator fit = - cdt.finite_faces_begin(); - fit != cdt.finite_faces_end(); ++fit) { - if (!fit->info().in_domain()) continue; - auto p = CGAL::centroid(cdt.triangle(fit)); - - // !! this turns out to be messy - // auto snap_triangle = tri.locate(p); - // arrFace->data() = *snap_triangle->info(); - - auto obj = - walk_pl.locate(Walk_pl::Arrangement_2::Point_2(p.x(), p.y())); - - if (auto f = std::get_if( - &obj)) { // located inside a face - // arrFace->data() = (*f)->data(); - canidate_faces[arr.non_const_handle(*f)] += - cdt.triangle(fit).area(); - } - // break; - } + typedef CGAL::Arr_walk_along_line_point_location + Snap_walk_pl; + Snap_walk_pl snap_walk_pl(arr_snap); + std::unordered_map> + output_face_labels; + for (auto face : tri.finite_face_handles()) { + if (face->info() == nullptr) continue; + auto centroid = CGAL::centroid(tri.triangle(face)); + auto object = snap_walk_pl.locate( + Arrangement_2::Point_2(centroid.x(), centroid.y())); + if (auto output_face = std::get_if(&object)) { + output_face_labels[arr_snap.non_const_handle(*output_face)] + [face->info()] += + std::abs(tri.triangle(face).area()); + } + } - // pick the candidate with the largest overlapping area - // std::cerr << "Size=" << canidate_faces.size() << std::endl; - if (canidate_faces.size()) { - auto best_face = std::max_element( - canidate_faces.begin(), canidate_faces.end(), - [](const std::pair& p1, - const std::pair& p2) { - return p1.second < p2.second; - }); - arrFace->data() = best_face->first->data(); - } else { - std::cout << "Unable to locate overlapping triangle\n"; - arrFace->data().is_ground = true; - arrFace->data().in_footprint = false; - arrFace->data().is_footprint_hole = false; - } + arr_snap.unbounded_face()->data() = arr.unbounded_face()->data(); + for (auto output_face : arr_snap.face_handles()) { + if (output_face->is_unbounded()) continue; + const auto labels = output_face_labels.find(output_face); + if (labels == output_face_labels.end() || labels->second.empty()) { + throw roofer::rooferException( + "Unable to transfer snapped arrangement face label"); } + auto best = + std::max_element(labels->second.begin(), labels->second.end(), + [&](const auto& lhs, const auto& rhs) { + if (lhs.second != rhs.second) + return lhs.second < rhs.second; + return source_face_ids.at(lhs.first) > + source_face_ids.at(rhs.first); + }); + output_face->data() = *best->first; } // remove dangling edges if any, eg holes that collapse to a single edge diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 47d139aa..00ac5f8e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -58,6 +58,12 @@ target_link_libraries("test_arrangement_dissolver" PRIVATE Catch2::Catch2WithMain roofer-core) catch_discover_tests("test_arrangement_dissolver") +add_executable("test_arrangement_snapper" + "${CMAKE_CURRENT_SOURCE_DIR}/test_arrangement_snapper.cpp") +target_link_libraries("test_arrangement_snapper" + PRIVATE Catch2::Catch2WithMain roofer-core) +catch_discover_tests("test_arrangement_snapper") + add_executable("test_vector_reader" "${CMAKE_CURRENT_SOURCE_DIR}/test_vector_reader.cpp") target_link_libraries("test_vector_reader" PRIVATE Catch2::Catch2WithMain) diff --git a/tests/test_arrangement_snapper.cpp b/tests/test_arrangement_snapper.cpp new file mode 100644 index 00000000..28de4ef4 --- /dev/null +++ b/tests/test_arrangement_snapper.cpp @@ -0,0 +1,129 @@ +#include +#include +#include + +#include + +#include + +#include +#include + +namespace { + using roofer::Arrangement_2; + using roofer::Point_2; + using roofer::Segment_2; + using Face_const_handle = Arrangement_2::Face_const_handle; + + Arrangement_2 cross_arrangement(const std::array& heights) { + Arrangement_2 arrangement; + const std::array corners = {Point_2(0, 0), Point_2(10, 0), + Point_2(10, 10), Point_2(0, 10)}; + for (std::size_t i = 0; i < corners.size(); ++i) { + CGAL::insert(arrangement, + Segment_2(corners[i], corners[(i + 1) % corners.size()])); + } + CGAL::insert(arrangement, Segment_2(Point_2(0, 5), Point_2(10, 5))); + CGAL::insert(arrangement, Segment_2(Point_2(5, 0), Point_2(5, 10))); + + using PointLocation = + CGAL::Arr_walk_along_line_point_location; + PointLocation point_location(arrangement); + const std::array samples = { + Point_2(7.5, 7.5), Point_2(2.5, 7.5), Point_2(2.5, 2.5), + Point_2(7.5, 2.5)}; + for (std::size_t i = 0; i < samples.size(); ++i) { + auto object = point_location.locate(samples[i]); + auto face = std::get_if(&object); + REQUIRE(face != nullptr); + auto mutable_face = arrangement.non_const_handle(*face); + mutable_face->data().in_footprint = true; + mutable_face->data().segid = static_cast(i + 1); + mutable_face->data().plane = roofer::Plane(0, 0, 1, -heights[i]); + } + return arrangement; + } + + Arrangement_2::Face_handle face_at(Arrangement_2& arrangement, + const Point_2& point) { + using PointLocation = + CGAL::Arr_walk_along_line_point_location; + PointLocation point_location(arrangement); + auto object = point_location.locate(point); + auto face = std::get_if(&object); + REQUIRE(face != nullptr); + return arrangement.non_const_handle(*face); + } + + Arrangement_2::Vertex_handle vertex_at(Arrangement_2& arrangement, + const Point_2& point) { + for (auto vertex : arrangement.vertex_handles()) { + if (vertex->point() == point) return vertex; + } + return {}; + } + + void check_two_manifold_edges(const roofer::Mesh& mesh) { + using Edge = std::pair; + std::map edge_incidence; + + auto add_ring = [&](const roofer::vec3f& ring) { + for (std::size_t i = 0; i < ring.size(); ++i) { + auto first = ring[i]; + auto second = ring[(i + 1) % ring.size()]; + if (second < first) std::swap(first, second); + if (first != second) ++edge_incidence[{first, second}]; + } + }; + for (const auto& polygon : mesh.get_polygons()) { + add_ring(polygon); + for (const auto& hole : polygon.interior_rings()) add_ring(hole); + } + for (const auto& [edge, incidence] : edge_incidence) { + CAPTURE(edge.first, edge.second); + CHECK(incidence == 2); + } + } +} // namespace + +TEST_CASE("snapper repairs an alternating four-way roof junction") { + // Heights in NE, NW, SW, SE order. The uniquely lowest NW face must absorb + // the repair area. + auto arrangement = cross_arrangement({15, 5, 15, 6}); + auto snapper = roofer::reconstruction::createArrangementSnapper(); + snapper->compute(arrangement, {.dist_thres = 0.001F, + .repair_non_manifold_vertices = true, + .manifold_repair_radius = 0.5F, + .manifold_height_tolerance = 1e-4F}); + + CHECK(vertex_at(arrangement, Point_2(5, 5)) == + Arrangement_2::Vertex_handle()); + CHECK(face_at(arrangement, Point_2(5, 5))->data().segid == 2); + + std::size_t roof_faces = 0; + for (auto face : arrangement.face_handles()) { + if (face->data().in_footprint) ++roof_faces; + } + CHECK(roof_faces == 4); + for (auto vertex : arrangement.vertex_handles()) { + CHECK(vertex->degree() <= 3); + } + + auto extruder = roofer::reconstruction::createArrangementExtruder(); + extruder->compute(arrangement, 0.0F); + REQUIRE(extruder->meshes.size() == 1); + check_two_manifold_edges(extruder->meshes.front()); +} + +TEST_CASE("snapper preserves a manifold four-way roof junction") { + auto arrangement = cross_arrangement({0, 1, 2, 3}); + auto snapper = roofer::reconstruction::createArrangementSnapper(); + snapper->compute(arrangement, {.dist_thres = 0.001F, + .repair_non_manifold_vertices = true, + .manifold_repair_radius = 0.5F, + .manifold_height_tolerance = 1e-4F}); + + auto center = vertex_at(arrangement, Point_2(5, 5)); + REQUIRE(center != Arrangement_2::Vertex_handle()); + CHECK(center->degree() == 4); +} From 3c652bd3c5fe90bd5d1c8af915da701dccc054f9 Mon Sep 17 00:00:00 2001 From: ylannl Date: Tue, 16 Jun 2026 10:25:08 +0200 Subject: [PATCH 09/21] simplify manifold edge fixing code --- .../reconstruction/ArrangementSnapper.cpp | 186 ++++++++---------- 1 file changed, 80 insertions(+), 106 deletions(-) diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index 851d0ba6..309c477b 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -64,7 +64,7 @@ namespace roofer::reconstruction { typedef T::Face_handle Face_handle; typedef std::pair Edge; - typedef std::unordered_map ConstraintMap; + typedef std::unordered_set ConstraintSet; struct ForcedRegionLabel { T::Point_2 seed; @@ -80,6 +80,12 @@ namespace roofer::reconstruction { double height; }; + struct RepairCandidate { + Vertex_handle vertex; + std::vector sectors; + double clearance; + }; + double plane_height(const FaceInfo& face_info, const T::Point_2& point) { const auto& plane = face_info.plane; if (plane.c() == 0) return 0; @@ -134,6 +140,14 @@ namespace roofer::reconstruction { CGAL::Arr_walk_along_line_point_location& walk_pl, Arrangement_2& source_arrangement, const SourceFaceIds& source_ids, const std::vector& forced_labels) { + std::vector> forced_face_labels; + for (const auto& forced : forced_labels) { + auto located = tri.locate(forced.seed); + if (!tri.is_infinite(located)) { + forced_face_labels.emplace_back(located, forced.face_info); + } + } + for (auto& region : constrained_regions(tri)) { std::unordered_map overlap_area; for (auto face : region) { @@ -157,12 +171,10 @@ namespace roofer::reconstruction { } } - for (const auto& forced : forced_labels) { - auto located = tri.locate(forced.seed); - if (!tri.is_infinite(located) && - std::find(region.begin(), region.end(), located) != - region.end()) { - selected = forced.face_info; + 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; } } @@ -186,7 +198,18 @@ namespace roofer::reconstruction { return std::sqrt(clearance_sq); } - std::vector incident_sectors(T& tri, Vertex_handle vertex) { + std::size_t constrained_incident_edge_count(T& tri, Vertex_handle vertex) { + std::size_t count = 0; + Edge_circulator edge = tri.incident_edges(vertex), done(edge); + if (edge == nullptr) return count; + do { + if (tri.is_constrained(*edge)) ++count; + } while (++edge != done); + return count; + } + + std::vector incident_sectors(T& tri, Vertex_handle vertex, + double clearance) { std::vector sectors; Edge_circulator edge = tri.incident_edges(vertex), done(edge); if (edge == nullptr) return sectors; @@ -206,9 +229,7 @@ namespace roofer::reconstruction { }); if (sectors.size() < 2) return sectors; - auto clearance = vertex_clearance(tri, vertex); - if (!clearance) return {}; - const double sample_radius = *clearance * 0.25; + const double sample_radius = clearance * 0.25; constexpr double two_pi = 2 * CGAL_PI; for (std::size_t i = 0; i < sectors.size(); ++i) { double next_angle = sectors[(i + 1) % sectors.size()].angle; @@ -258,12 +279,15 @@ namespace roofer::reconstruction { return false; } - std::optional find_problematic_vertex( + std::optional find_problematic_vertex( T& tri, double height_tolerance) { for (auto vertex = tri.finite_vertices_begin(); vertex != tri.finite_vertices_end(); ++vertex) { if (vertex->info()) continue; - auto sectors = incident_sectors(tri, vertex); + if (constrained_incident_edge_count(tri, vertex) < 4) continue; + auto clearance = vertex_clearance(tri, vertex); + if (!clearance) continue; + auto sectors = incident_sectors(tri, vertex, *clearance); if (sectors.size() < 4) continue; if (std::any_of(sectors.begin(), sectors.end(), [](const auto& sector) { return !sector.face_info->in_footprint || @@ -272,31 +296,40 @@ namespace roofer::reconstruction { continue; } if (has_non_manifold_height_order(sectors, height_tolerance)) { - return vertex; + return RepairCandidate{vertex, std::move(sectors), *clearance}; } } return std::nullopt; } - ForcedRegionLabel repair_vertex(T& tri, Vertex_handle vertex, + ForcedRegionLabel repair_vertex(T& tri, const RepairCandidate& candidate, const SourceFaceIds& source_ids, double maximum_radius, double height_tolerance) { - auto sectors = incident_sectors(tri, vertex); - auto clearance = vertex_clearance(tri, vertex); - if (sectors.size() < 4 || !clearance) { + const auto& sectors = candidate.sectors; + auto vertex = candidate.vertex; + if (sectors.size() < 4) { throw roofer::rooferException( "Unable to compute a safe non-manifold junction repair"); } - auto selected = std::min_element( - sectors.begin(), sectors.end(), - [&](const auto& lhs, const auto& rhs) { - if (std::abs(lhs.height - rhs.height) >= height_tolerance) - return lhs.height < rhs.height; - return source_ids.at(lhs.face_info) < source_ids.at(rhs.face_info); - }); - FaceInfo* selected_face = selected->face_info; + const double min_height = + std::min_element(sectors.begin(), sectors.end(), + [](const auto& lhs, const auto& rhs) { + return lhs.height < rhs.height; + }) + ->height; + FaceInfo* selected_face = nullptr; + std::size_t selected_id = std::numeric_limits::max(); + for (const auto& sector : sectors) { + if (sector.height <= min_height + height_tolerance) { + auto id = source_ids.at(sector.face_info); + if (id < selected_id) { + selected_face = sector.face_info; + selected_id = id; + } + } + } double shortest_edge = std::numeric_limits::max(); for (const auto& sector : sectors) { @@ -304,8 +337,8 @@ namespace roofer::reconstruction { shortest_edge, std::sqrt(CGAL::to_double(CGAL::squared_distance( vertex->point(), sector.neighbour->point())))); } - const double radius = - std::min({maximum_radius, *clearance * 0.8, shortest_edge * 0.45}); + const double radius = std::min( + {maximum_radius, candidate.clearance * 0.8, shortest_edge * 0.45}); if (!(radius > std::numeric_limits::epsilon())) { throw roofer::rooferException( "Non-manifold junction repair radius is degenerate"); @@ -378,77 +411,40 @@ namespace roofer::reconstruction { void get_incident_constraints(T& tri, Vertex_handle& vthis, Vertex_handle& vexcept1, Vertex_handle vexcept2, - ConstraintMap& constraints_to_restore) { + ConstraintSet& constraints_to_restore) { // std::cout << "vthis degree=" << tri.degree(vthis) << std::endl; Edge_circulator ec = tri.incident_edges(vthis), done(ec); if (ec != nullptr) { do { if (tri.is_constrained(*ec)) { - // figure out which is the other vertex on this edge - // and what is the face counter clockwise to the edge vthis -> - // vother auto va = ec->first->vertex(tri.cw(ec->second)); auto vb = ec->first->vertex(tri.ccw(ec->second)); Vertex_handle vother; - Face_handle face_ccw; if (va == vthis) { vother = vb; - face_ccw = ec->first->neighbor(ec->second); } else { vother = va; - face_ccw = ec->first; }; - // check if face_ccw will be collapsed. If so we should not restore - // it. vmirror is the vertex opposing ec in face_ccw face_ccw will - // be collapsed if vmirror is one of the vexcept vertices auto i = - // face_ccw->index(vthis); auto vmirror = face_ccw->vertex( - // tri.cw(i) ); if ( - // vmirror == vexcept1 || - // vmirror == vexcept2 - // ) { - // continue; - // } - - // check if edge is not one that we will collapse if (vother != vexcept1 && vother != vexcept2) { - constraints_to_restore[vother] = face_ccw->info(); + constraints_to_restore.insert(vother); } - // to_remove.push_back(*ec); } } while (++ec != done); } } void restore_constraints(T& tri, T::Point_2& pnew, - ConstraintMap& constraints_to_restore) { + ConstraintSet& constraints_to_restore) { auto vnew = tri.insert(pnew); // restore constraints // std::cout << "restoring " << constraints_to_restore.size() << " // constraints\n"; - for (auto& [vh, finfo] : constraints_to_restore) { + for (auto& vh : constraints_to_restore) { // std::cout << "reinsert constrained " << *vnew << " - " << *vh << // std::endl; tri.insert_constraint(vnew, vh); - - // find the ccw face to the edge vnew->vh - // Face_circulator fc = tri.incident_faces(vnew), - // done(fc); - // if (fc != nullptr) { - // do { - // if (fc->has_vertex(vh)) { - // // check if fc is ccw to edge vnew->vh - // if (vh == fc->vertex( tri.ccw( fc->index(vnew) ) ) ) { - // fc->info() = finfo; - // } else { - // // fc is cw to edge, so find the neighbor on the other side - // of edge fc->neighbor( tri.cw( fc->index(vh) ) )->info() = - // finfo; - // } - // } - // } while (++fc != done); - // } } } @@ -517,34 +513,6 @@ namespace roofer::reconstruction { // } // } - // copy semantics - for (auto fit : tri.finite_face_handles()) { - auto p = CGAL::centroid(tri.triangle(fit)); - auto obj = - walk_pl.locate(Walk_pl::Arrangement_2::Point_2(p.x(), p.y())); - - // std::cout << "The point (" << p << ") is located "; - if (auto f = std::get_if( - &obj)) { // located inside a face - // std::cout << "inside " - // << (((*f)->is_unbounded()) ? "the unbounded" : "a - // bounded") - // << " face." << std::endl; - fit->info() = &(arr.non_const_handle(*f))->data(); - } else { - // std::cout << "Failed point location in arrangement.\n"; - } - // else if (auto e = boost::get(&obj)) // - // located on an edge - // std::cout << "on an edge: " << (*e)->curve() << std::endl; - // else if (auto v = boost::get(&obj)) // located - // on a vertex - // std::cout << "on " << (((*v)->is_isolated()) ? "an isolated" : - // "a") - // << " vertex: " << (*v)->point() << std::endl; - // else CGAL_error_msg("Invalid object."); - } - // TriangleCollection triangles_og; // vec1i segment_ids_og; // for (auto fh = tri.finite_faces_begin(); fh != @@ -595,7 +563,7 @@ namespace roofer::reconstruction { ) { // std::cout << "small triangle between " << p0 << " and " << p1 // << " and " << p2 << std::endl; - ConstraintMap constraints_to_restore; + ConstraintSet constraints_to_restore; // collect incident constraint edges // std::vector to_remove; @@ -660,7 +628,7 @@ namespace roofer::reconstruction { // << std::endl; // auto vi = ceit->second; - ConstraintMap constraints_to_restore; + ConstraintSet constraints_to_restore; get_incident_constraints(tri, v1, v1, v2, constraints_to_restore); get_incident_constraints(tri, v2, v1, v2, @@ -749,10 +717,10 @@ namespace roofer::reconstruction { label_final_regions(tri, walk_pl, arr, source_face_ids, forced_region_labels); if (cfg.repair_non_manifold_vertices) { - while (auto vertex = find_problematic_vertex( + while (auto candidate = find_problematic_vertex( tri, cfg.manifold_height_tolerance)) { forced_region_labels.push_back(repair_vertex( - tri, *vertex, source_face_ids, cfg.manifold_repair_radius, + tri, *candidate, source_face_ids, cfg.manifold_repair_radius, cfg.manifold_height_tolerance)); label_final_regions(tri, walk_pl, arr, source_face_ids, forced_region_labels); @@ -818,15 +786,21 @@ namespace roofer::reconstruction { std::unordered_map> output_face_labels; - for (auto face : tri.finite_face_handles()) { - if (face->info() == nullptr) continue; - auto centroid = CGAL::centroid(tri.triangle(face)); + for (const auto& region : constrained_regions(tri)) { + FaceInfo* region_label = nullptr; + double region_area = 0; + for (auto face : region) { + if (region_label == nullptr) region_label = face->info(); + region_area += std::abs(tri.triangle(face).area()); + } + if (region_label == nullptr || region_area == 0) continue; + + auto centroid = CGAL::centroid(tri.triangle(region.front())); auto object = snap_walk_pl.locate( Arrangement_2::Point_2(centroid.x(), centroid.y())); if (auto output_face = std::get_if(&object)) { output_face_labels[arr_snap.non_const_handle(*output_face)] - [face->info()] += - std::abs(tri.triangle(face).area()); + [region_label] += region_area; } } From e89746bd65a5906515000bcef6c5e5c44e314353 Mon Sep 17 00:00:00 2001 From: ylannl Date: Tue, 16 Jun 2026 10:40:54 +0200 Subject: [PATCH 10/21] apply non-manifold edge fix approach also to prevent face self-intersect --- .../reconstruction/ArrangementSnapper.cpp | 63 +++++++++++++------ tests/test_arrangement_snapper.cpp | 24 +++++++ 2 files changed, 69 insertions(+), 18 deletions(-) diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index 309c477b..8dc26bcb 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -84,6 +84,7 @@ namespace roofer::reconstruction { Vertex_handle vertex; std::vector sectors; double clearance; + FaceInfo* preferred_merge_face = nullptr; }; double plane_height(const FaceInfo& face_info, const T::Point_2& point) { @@ -279,8 +280,29 @@ namespace roofer::reconstruction { return false; } + FaceInfo* repeated_incident_face(const std::vector& sectors, + const SourceFaceIds& source_ids) { + std::unordered_map counts; + for (const auto& sector : sectors) ++counts[sector.face_info]; + + FaceInfo* selected = nullptr; + std::size_t selected_count = 0; + std::size_t selected_id = std::numeric_limits::max(); + for (const auto& [face_info, count] : counts) { + if (count < 2) continue; + auto id = source_ids.at(face_info); + if (count > selected_count || + (count == selected_count && id < selected_id)) { + selected = face_info; + selected_count = count; + selected_id = id; + } + } + return selected; + } + std::optional find_problematic_vertex( - T& tri, double height_tolerance) { + T& tri, const SourceFaceIds& source_ids, double height_tolerance) { for (auto vertex = tri.finite_vertices_begin(); vertex != tri.finite_vertices_end(); ++vertex) { if (vertex->info()) continue; @@ -295,8 +317,11 @@ namespace roofer::reconstruction { })) { continue; } - if (has_non_manifold_height_order(sectors, height_tolerance)) { - return RepairCandidate{vertex, std::move(sectors), *clearance}; + auto repeated_face = repeated_incident_face(sectors, source_ids); + if (repeated_face || + has_non_manifold_height_order(sectors, height_tolerance)) { + return RepairCandidate{vertex, std::move(sectors), *clearance, + repeated_face}; } } return std::nullopt; @@ -313,20 +338,22 @@ namespace roofer::reconstruction { "Unable to compute a safe non-manifold junction repair"); } - const double min_height = - std::min_element(sectors.begin(), sectors.end(), - [](const auto& lhs, const auto& rhs) { - return lhs.height < rhs.height; - }) - ->height; - FaceInfo* selected_face = nullptr; - std::size_t selected_id = std::numeric_limits::max(); - for (const auto& sector : sectors) { - if (sector.height <= min_height + height_tolerance) { - auto id = source_ids.at(sector.face_info); - if (id < selected_id) { - selected_face = sector.face_info; - selected_id = id; + FaceInfo* selected_face = candidate.preferred_merge_face; + if (selected_face == nullptr) { + const double min_height = + std::min_element(sectors.begin(), sectors.end(), + [](const auto& lhs, const auto& rhs) { + return lhs.height < rhs.height; + }) + ->height; + std::size_t selected_id = std::numeric_limits::max(); + for (const auto& sector : sectors) { + if (sector.height <= min_height + height_tolerance) { + auto id = source_ids.at(sector.face_info); + if (id < selected_id) { + selected_face = sector.face_info; + selected_id = id; + } } } } @@ -718,7 +745,7 @@ namespace roofer::reconstruction { forced_region_labels); if (cfg.repair_non_manifold_vertices) { while (auto candidate = find_problematic_vertex( - tri, cfg.manifold_height_tolerance)) { + tri, source_face_ids, cfg.manifold_height_tolerance)) { forced_region_labels.push_back(repair_vertex( tri, *candidate, source_face_ids, cfg.manifold_repair_radius, cfg.manifold_height_tolerance)); diff --git a/tests/test_arrangement_snapper.cpp b/tests/test_arrangement_snapper.cpp index 28de4ef4..08f69bf7 100644 --- a/tests/test_arrangement_snapper.cpp +++ b/tests/test_arrangement_snapper.cpp @@ -55,6 +55,15 @@ namespace { return arrangement.non_const_handle(*face); } + Arrangement_2 repeated_face_arrangement() { + auto arrangement = cross_arrangement({1, 2, 1, 3}); + auto repeated_face = face_at(arrangement, Point_2(7.5, 7.5)); + auto other_face = face_at(arrangement, Point_2(2.5, 2.5)); + other_face->data().segid = repeated_face->data().segid; + other_face->data().plane = repeated_face->data().plane; + return arrangement; + } + Arrangement_2::Vertex_handle vertex_at(Arrangement_2& arrangement, const Point_2& point) { for (auto vertex : arrangement.vertex_handles()) { @@ -127,3 +136,18 @@ TEST_CASE("snapper preserves a manifold four-way roof junction") { REQUIRE(center != Arrangement_2::Vertex_handle()); CHECK(center->degree() == 4); } + +TEST_CASE("snapper repairs a repeated-face junction") { + auto arrangement = repeated_face_arrangement(); + auto repeated_segid = face_at(arrangement, Point_2(7.5, 7.5))->data().segid; + + auto snapper = roofer::reconstruction::createArrangementSnapper(); + snapper->compute(arrangement, {.dist_thres = 0.001F, + .repair_non_manifold_vertices = true, + .manifold_repair_radius = 0.5F, + .manifold_height_tolerance = 1e-4F}); + + CHECK(vertex_at(arrangement, Point_2(5, 5)) == + Arrangement_2::Vertex_handle()); + CHECK(face_at(arrangement, Point_2(5, 5))->data().segid == repeated_segid); +} From 957600d8ef6e095f8994c40947d7e8465ecdf9be Mon Sep 17 00:00:00 2001 From: ylannl Date: Tue, 16 Jun 2026 12:16:49 +0200 Subject: [PATCH 11/21] attempt to handle issues on roof boundary --- .../reconstruction/ArrangementSnapper.cpp | 72 ++++++++++++------- tests/test_arrangement_snapper.cpp | 72 +++++++++++++++++++ 2 files changed, 119 insertions(+), 25 deletions(-) diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index 8dc26bcb..6a58fdd0 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -94,6 +94,11 @@ namespace roofer::reconstruction { plane.c(); } + bool is_roof_face(const FaceInfo* face_info) { + return face_info != nullptr && face_info->in_footprint && + face_info->segid != 0; + } + FaceInfo* source_face_at( const T::Point_2& point, CGAL::Arr_walk_along_line_point_location& walk_pl, @@ -189,13 +194,15 @@ namespace roofer::reconstruction { Face_circulator face = tri.incident_faces(vertex), done(face); if (face == nullptr) return std::nullopt; do { - if (tri.is_infinite(face)) return std::nullopt; + if (tri.is_infinite(face)) continue; auto opposite_edge = Edge(face, face->index(vertex)); clearance_sq = std::min( clearance_sq, CGAL::to_double(CGAL::squared_distance( vertex->point(), tri.segment(opposite_edge)))); } while (++face != done); - if (!(clearance_sq > 0)) return std::nullopt; + if (!(clearance_sq > 0) || + clearance_sq == std::numeric_limits::max()) + return std::nullopt; return std::sqrt(clearance_sq); } @@ -209,8 +216,10 @@ namespace roofer::reconstruction { return count; } - std::vector incident_sectors(T& tri, Vertex_handle vertex, - double clearance) { + std::vector incident_sectors( + T& tri, Vertex_handle vertex, double clearance, + Arrangement_2& source_arrangement, + CGAL::Arr_walk_along_line_point_location& walk_pl) { std::vector sectors; Edge_circulator edge = tri.incident_edges(vertex), done(edge); if (edge == nullptr) return sectors; @@ -239,10 +248,10 @@ namespace roofer::reconstruction { T::Point_2 sample( vertex->point().x() + sample_radius * std::cos(sample_angle), vertex->point().y() + sample_radius * std::sin(sample_angle)); - auto face = tri.locate(sample); - if (tri.is_infinite(face) || face->info() == nullptr) return {}; - sectors[i].face_info = face->info(); - sectors[i].height = plane_height(*face->info(), vertex->point()); + auto face_info = source_face_at(sample, walk_pl, source_arrangement); + if (face_info == nullptr) return {}; + sectors[i].face_info = face_info; + sectors[i].height = plane_height(*face_info, vertex->point()); } return sectors; } @@ -283,7 +292,9 @@ namespace roofer::reconstruction { FaceInfo* repeated_incident_face(const std::vector& sectors, const SourceFaceIds& source_ids) { std::unordered_map counts; - for (const auto& sector : sectors) ++counts[sector.face_info]; + for (const auto& sector : sectors) { + if (is_roof_face(sector.face_info)) ++counts[sector.face_info]; + } FaceInfo* selected = nullptr; std::size_t selected_count = 0; @@ -302,19 +313,21 @@ namespace roofer::reconstruction { } std::optional find_problematic_vertex( - T& tri, const SourceFaceIds& source_ids, double height_tolerance) { + T& tri, Arrangement_2& source_arrangement, + CGAL::Arr_walk_along_line_point_location& walk_pl, + const SourceFaceIds& source_ids, double height_tolerance) { for (auto vertex = tri.finite_vertices_begin(); vertex != tri.finite_vertices_end(); ++vertex) { - if (vertex->info()) continue; if (constrained_incident_edge_count(tri, vertex) < 4) continue; auto clearance = vertex_clearance(tri, vertex); if (!clearance) continue; - auto sectors = incident_sectors(tri, vertex, *clearance); + auto sectors = incident_sectors(tri, vertex, *clearance, + source_arrangement, walk_pl); if (sectors.size() < 4) continue; - if (std::any_of(sectors.begin(), sectors.end(), [](const auto& sector) { - return !sector.face_info->in_footprint || - sector.face_info->segid == 0; - })) { + if (std::none_of(sectors.begin(), sectors.end(), + [](const auto& sector) { + return is_roof_face(sector.face_info); + })) { continue; } auto repeated_face = repeated_incident_face(sectors, source_ids); @@ -340,15 +353,23 @@ namespace roofer::reconstruction { FaceInfo* selected_face = candidate.preferred_merge_face; if (selected_face == nullptr) { - const double min_height = - std::min_element(sectors.begin(), sectors.end(), - [](const auto& lhs, const auto& rhs) { - return lhs.height < rhs.height; - }) - ->height; + auto min_sector = std::min_element( + sectors.begin(), sectors.end(), + [](const auto& lhs, const auto& rhs) { + if (is_roof_face(lhs.face_info) != is_roof_face(rhs.face_info)) + return is_roof_face(lhs.face_info); + return lhs.height < rhs.height; + }); + if (min_sector == sectors.end() || + !is_roof_face(min_sector->face_info)) { + throw roofer::rooferException( + "Non-manifold junction repair has no roof face to merge into"); + } + const double min_height = min_sector->height; std::size_t selected_id = std::numeric_limits::max(); for (const auto& sector : sectors) { - if (sector.height <= min_height + height_tolerance) { + if (is_roof_face(sector.face_info) && + sector.height <= min_height + height_tolerance) { auto id = source_ids.at(sector.face_info); if (id < selected_id) { selected_face = sector.face_info; @@ -744,8 +765,9 @@ namespace roofer::reconstruction { label_final_regions(tri, walk_pl, arr, source_face_ids, forced_region_labels); if (cfg.repair_non_manifold_vertices) { - while (auto candidate = find_problematic_vertex( - tri, source_face_ids, cfg.manifold_height_tolerance)) { + while (auto candidate = + find_problematic_vertex(tri, arr, walk_pl, source_face_ids, + cfg.manifold_height_tolerance)) { forced_region_labels.push_back(repair_vertex( tri, *candidate, source_face_ids, cfg.manifold_repair_radius, cfg.manifold_height_tolerance)); diff --git a/tests/test_arrangement_snapper.cpp b/tests/test_arrangement_snapper.cpp index 08f69bf7..2af69366 100644 --- a/tests/test_arrangement_snapper.cpp +++ b/tests/test_arrangement_snapper.cpp @@ -64,6 +64,50 @@ namespace { return arrangement; } + Arrangement_2 boundary_junction_arrangement() { + Arrangement_2 arrangement; + const std::array corners = {Point_2(0, 0), Point_2(10, 0), + Point_2(10, 10), Point_2(0, 10)}; + for (std::size_t i = 0; i < corners.size(); ++i) { + CGAL::insert(arrangement, + Segment_2(corners[i], corners[(i + 1) % corners.size()])); + } + CGAL::insert(arrangement, Segment_2(Point_2(5, 0), Point_2(2, 10))); + CGAL::insert(arrangement, Segment_2(Point_2(5, 0), Point_2(8, 10))); + + auto left = face_at(arrangement, Point_2(1, 5)); + left->data().in_footprint = true; + left->data().segid = 1; + left->data().plane = roofer::Plane(0, 0, 1, -15); + + auto middle = face_at(arrangement, Point_2(5, 5)); + middle->data().in_footprint = true; + middle->data().segid = 2; + middle->data().plane = roofer::Plane(0, 0, 1, -5); + + auto right = face_at(arrangement, Point_2(9, 5)); + right->data().in_footprint = true; + right->data().segid = 3; + right->data().plane = roofer::Plane(0, 0, 1, -15); + + arrangement.unbounded_face()->data().in_footprint = false; + arrangement.unbounded_face()->data().segid = 0; + arrangement.unbounded_face()->data().plane = roofer::Plane(0, 0, 1, 10); + + return arrangement; + } + + Arrangement_2 finite_exterior_junction_arrangement() { + auto arrangement = cross_arrangement({15, 5, 15, 0}); + + auto hole = face_at(arrangement, Point_2(7.5, 2.5)); + hole->data().in_footprint = false; + hole->data().is_footprint_hole = true; + hole->data().segid = 0; + hole->data().plane = roofer::Plane(0, 0, 1, 10); + return arrangement; + } + Arrangement_2::Vertex_handle vertex_at(Arrangement_2& arrangement, const Point_2& point) { for (auto vertex : arrangement.vertex_handles()) { @@ -151,3 +195,31 @@ TEST_CASE("snapper repairs a repeated-face junction") { Arrangement_2::Vertex_handle()); CHECK(face_at(arrangement, Point_2(5, 5))->data().segid == repeated_segid); } + +TEST_CASE("snapper keeps a boundary repair cell inside the footprint") { + auto arrangement = boundary_junction_arrangement(); + + auto snapper = roofer::reconstruction::createArrangementSnapper(); + snapper->compute(arrangement, {.dist_thres = 0.001F, + .repair_non_manifold_vertices = true, + .manifold_repair_radius = 0.5F, + .manifold_height_tolerance = 1e-4F}); + + CHECK(vertex_at(arrangement, Point_2(5, 0)) == + Arrangement_2::Vertex_handle()); + CHECK(face_at(arrangement, Point_2(5, 0.05))->data().in_footprint); +} + +TEST_CASE("snapper keeps a finite exterior repair cell inside the footprint") { + auto arrangement = finite_exterior_junction_arrangement(); + + auto snapper = roofer::reconstruction::createArrangementSnapper(); + snapper->compute(arrangement, {.dist_thres = 0.001F, + .repair_non_manifold_vertices = true, + .manifold_repair_radius = 0.5F, + .manifold_height_tolerance = 1e-4F}); + + CHECK(vertex_at(arrangement, Point_2(5, 5)) == + Arrangement_2::Vertex_handle()); + CHECK(face_at(arrangement, Point_2(5, 5))->data().segid == 2); +} From ccc5614c663dcfd0e86190e7e83910abd5d261bd Mon Sep 17 00:00:00 2001 From: ylannl Date: Tue, 16 Jun 2026 23:26:12 +0200 Subject: [PATCH 12/21] fix for edge case of self intersection detection and repair --- .../reconstruction/ArrangementSnapper.cpp | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index 6a58fdd0..dff29bcb 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -289,6 +289,27 @@ namespace roofer::reconstruction { return false; } + std::size_t cyclic_face_run_count( + const std::vector& sectors, const FaceInfo* face_info) { + if (sectors.empty()) return 0; + if (std::all_of(sectors.begin(), sectors.end(), + [face_info](const auto& sector) { + return sector.face_info == face_info; + })) { + return 1; + } + + std::size_t runs = 0; + for (std::size_t i = 0; i < sectors.size(); ++i) { + const auto previous = (i + sectors.size() - 1) % sectors.size(); + if (sectors[i].face_info == face_info && + sectors[previous].face_info != face_info) { + ++runs; + } + } + return runs; + } + FaceInfo* repeated_incident_face(const std::vector& sectors, const SourceFaceIds& source_ids) { std::unordered_map counts; @@ -301,6 +322,7 @@ namespace roofer::reconstruction { std::size_t selected_id = std::numeric_limits::max(); for (const auto& [face_info, count] : counts) { if (count < 2) continue; + if (cyclic_face_run_count(sectors, face_info) < 2) continue; auto id = source_ids.at(face_info); if (count > selected_count || (count == selected_count && id < selected_id)) { @@ -321,6 +343,7 @@ namespace roofer::reconstruction { if (constrained_incident_edge_count(tri, vertex) < 4) continue; auto clearance = vertex_clearance(tri, vertex); if (!clearance) continue; + if (*clearance <= 1e-6) continue; auto sectors = incident_sectors(tri, vertex, *clearance, source_arrangement, walk_pl); if (sectors.size() < 4) continue; From 3e27ebfbbbc4cbafa8a492eb3670dcf93ce2632e Mon Sep 17 00:00:00 2001 From: ylannl Date: Wed, 17 Jun 2026 21:05:21 +0200 Subject: [PATCH 13/21] do not use vertex references after removal --- .../reconstruction/ArrangementSnapper.cpp | 36 +++++++++---------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index dff29bcb..1d9f1d1c 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -479,8 +479,8 @@ namespace roofer::reconstruction { // return triangulation; // } - void get_incident_constraints(T& tri, Vertex_handle& vthis, - Vertex_handle& vexcept1, + void get_incident_constraints(T& tri, Vertex_handle vthis, + Vertex_handle vexcept1, Vertex_handle vexcept2, ConstraintSet& constraints_to_restore) { // std::cout << "vthis degree=" << tri.degree(vthis) << std::endl; @@ -637,7 +637,6 @@ namespace roofer::reconstruction { ConstraintSet constraints_to_restore; // collect incident constraint edges - // std::vector to_remove; for (size_t i = 0; i < 3; ++i) { auto vthis = fit->vertex(i); auto vexcept1 = fit->vertex(tri.cw(i)); @@ -646,16 +645,6 @@ namespace roofer::reconstruction { constraints_to_restore); } - // if (constraints_to_restore.size() > 50) continue; - - // collapse the triangle to centroid - tri.remove_incident_constraints(v0); - tri.remove(v0); - tri.remove_incident_constraints(v1); - tri.remove(v1); - tri.remove_incident_constraints(v2); - tri.remove(v2); - // but first check points for being on the footprint boundary T::Point_2 pnew; if (v0->info()) { @@ -667,6 +656,14 @@ namespace roofer::reconstruction { } else { pnew = CGAL::centroid(p0, p1, p2); } + // collapse the triangle to centroid + tri.remove_incident_constraints(v0); + tri.remove(v0); + tri.remove_incident_constraints(v1); + tri.remove(v1); + tri.remove_incident_constraints(v2); + tri.remove(v2); + restore_constraints(tri, pnew, constraints_to_restore); found_small_face = true; @@ -705,12 +702,6 @@ namespace roofer::reconstruction { get_incident_constraints(tri, v2, v1, v2, constraints_to_restore); - // remove edge - tri.remove_incident_constraints(v1); - tri.remove(v1); - tri.remove_incident_constraints(v2); - tri.remove(v2); - // insert midpoint and restore constraints // first check points for being on the footprint boundary T::Point_2 pnew; @@ -721,6 +712,13 @@ namespace roofer::reconstruction { } else { pnew = CGAL::midpoint(p1, p2); } + + // remove edge + tri.remove_incident_constraints(v1); + tri.remove(v1); + tri.remove_incident_constraints(v2); + tri.remove(v2); + restore_constraints(tri, pnew, constraints_to_restore); found_short_edge = true; From 9587be044401898730b67391cf7a5174ddbaf07a Mon Sep 17 00:00:00 2001 From: ylannl Date: Fri, 19 Jun 2026 17:37:11 +0200 Subject: [PATCH 14/21] add some comments --- .../reconstruction/ArrangementSnapper.cpp | 38 +++++++------------ 1 file changed, 14 insertions(+), 24 deletions(-) diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index 1d9f1d1c..d6ea73b6 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -111,6 +111,7 @@ namespace roofer::reconstruction { return nullptr; } + // region grows constrained regions of triangles std::vector> constrained_regions(T& tri) { std::vector> regions; std::unordered_set visited; @@ -141,6 +142,7 @@ namespace roofer::reconstruction { return regions; } + // labels the triangulation with region IDs based on the arrangement void label_final_regions( T& tri, CGAL::Arr_walk_along_line_point_location& walk_pl, @@ -155,6 +157,8 @@ namespace roofer::reconstruction { } for (auto& region : constrained_regions(tri)) { + // for this triangle region collect total overlap area for each + // arrangement face std::unordered_map overlap_area; for (auto face : region) { auto centroid = CGAL::centroid(tri.triangle(face)); @@ -164,6 +168,7 @@ 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(); @@ -177,6 +182,7 @@ namespace roofer::reconstruction { } } + // apply the forced regions (these are determined upstream) for (const auto& [forced_face, forced_label] : forced_face_labels) { if (std::find(region.begin(), region.end(), forced_face) != region.end()) { @@ -189,6 +195,9 @@ namespace roofer::reconstruction { } } + // Heuristically compute clearance: the minimum distance to and opposing + // edge from current vertex This is a conservative estimate, and a lower + // bound on the actual clearance std::optional vertex_clearance(T& tri, Vertex_handle vertex) { double clearance_sq = std::numeric_limits::max(); Face_circulator face = tri.incident_faces(vertex), done(face); @@ -216,6 +225,8 @@ namespace roofer::reconstruction { return count; } + // calculate azimuth for each incident triangulation edge, and sample + // corresponding face from arrangement std::vector incident_sectors( T& tri, Vertex_handle vertex, double clearance, Arrangement_2& source_arrangement, @@ -353,6 +364,7 @@ namespace roofer::reconstruction { })) { continue; } + // we check repeated incident faces for self intersection at the vertex auto repeated_face = repeated_incident_face(sectors, source_ids); if (repeated_face || has_non_manifold_height_order(sectors, height_tolerance)) { @@ -734,9 +746,6 @@ namespace roofer::reconstruction { // TODO: move the vertex opposed to long edge to be on the long edge ?? // TODO: detect if an edge of the triangle is on the footprint boundary - // bool found_small_face; - // do { - // found_small_face = false; for (Finite_faces_iterator fit = tri.finite_faces_begin(); fit != tri.finite_faces_end(); ++fit) { auto v0 = fit->vertex(0); @@ -780,7 +789,6 @@ namespace roofer::reconstruction { } } } - // } while (found_short_edge); std::vector forced_region_labels; label_final_regions(tri, walk_pl, arr, source_face_ids, @@ -797,27 +805,8 @@ namespace roofer::reconstruction { } } - // TriangleCollection triangles_snapped; - // vec1i segment_ids_snapped; - // for (auto fh = tri.finite_faces_begin(); fh != - // tri.finite_faces_end(); ++fh) { - // // only export triangles in the interior of a shape (thus excluding - // holes and exterior) - - // arr3f p0 = {float (fh->vertex(0)->point().x()), float - // (fh->vertex(0)->point().y()), 0}; arr3f p1 = {float - // (fh->vertex(1)->point().x()), float (fh->vertex(1)->point().y()), - // 0}; arr3f p2 = {float (fh->vertex(2)->point().x()), float - // (fh->vertex(2)->point().y()), 0}; triangles_snapped.push_back({ - // p0,p1,p2 }); - // // segment_ids_snapped.push_back(fh->info()->segid); - // // segment_ids_snapped.push_back(fh->info()->segid); - // // segment_ids_snapped.push_back(fh->info()->segid); - // } - // output("triangles_snapped").set(triangles_snapped); - // output("segment_ids_snapped").set(segment_ids_snapped); - // convert back from triangulation to arrangement + // 1 recreate vertices and faces Arrangement_2 arr_snap; std::unordered_map vertex2arr_map; @@ -850,6 +839,7 @@ namespace roofer::reconstruction { // } } + // 2 transfer labels from triangulation to new arrangement typedef CGAL::Arr_walk_along_line_point_location Snap_walk_pl; Snap_walk_pl snap_walk_pl(arr_snap); From b4c87fe5b97cb8909fcae30a0c4ae8eb2e207248 Mon Sep 17 00:00:00 2001 From: ylannl Date: Fri, 19 Jun 2026 18:43:40 +0200 Subject: [PATCH 15/21] simplify min face height selection --- .../reconstruction/ArrangementSnapper.cpp | 37 +++++++++---------- 1 file changed, 17 insertions(+), 20 deletions(-) diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index d6ea73b6..40d28f47 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -388,30 +388,24 @@ namespace roofer::reconstruction { FaceInfo* selected_face = candidate.preferred_merge_face; if (selected_face == nullptr) { - auto min_sector = std::min_element( - sectors.begin(), sectors.end(), - [](const auto& lhs, const auto& rhs) { - if (is_roof_face(lhs.face_info) != is_roof_face(rhs.face_info)) - return is_roof_face(lhs.face_info); - return lhs.height < rhs.height; - }); - if (min_sector == sectors.end() || - !is_roof_face(min_sector->face_info)) { - throw roofer::rooferException( - "Non-manifold junction repair has no roof face to merge into"); - } - const double min_height = min_sector->height; + double selected_height = std::numeric_limits::infinity(); std::size_t selected_id = std::numeric_limits::max(); + for (const auto& sector : sectors) { - if (is_roof_face(sector.face_info) && - sector.height <= min_height + height_tolerance) { - auto id = source_ids.at(sector.face_info); - if (id < selected_id) { - selected_face = sector.face_info; - selected_id = id; - } + if (!is_roof_face(sector.face_info)) continue; + + auto id = source_ids.at(sector.face_info); + if (sector.height < selected_height || + (sector.height == selected_height && id < selected_id)) { + selected_face = sector.face_info; + selected_height = sector.height; + selected_id = id; } } + if (selected_face == nullptr) { + throw roofer::rooferException( + "Non-manifold junction repair has no roof face to merge into"); + } } double shortest_edge = std::numeric_limits::max(); @@ -427,6 +421,7 @@ namespace roofer::reconstruction { "Non-manifold junction repair radius is degenerate"); } + // insert split vertices along each incident edge const auto original_point = vertex->point(); std::vector split_vertices; split_vertices.reserve(sectors.size()); @@ -441,12 +436,14 @@ namespace roofer::reconstruction { split_vertices.push_back(split_vertex); } + // 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) { tri.insert_constraint(split_vertices[i], sectors[i].neighbour); } + // std::optional forced_seed; for (std::size_t i = 0; i < sectors.size(); ++i) { auto next = (i + 1) % sectors.size(); From 2a668e817e6c4feae8fc4356e9ff581d1ad3fb82 Mon Sep 17 00:00:00 2001 From: ylannl Date: Fri, 19 Jun 2026 20:04:44 +0200 Subject: [PATCH 16/21] remove excessive relabeling --- src/core/reconstruction/ArrangementSnapper.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index 40d28f47..7dd3dd2b 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -788,8 +788,6 @@ namespace roofer::reconstruction { } std::vector forced_region_labels; - label_final_regions(tri, walk_pl, arr, source_face_ids, - forced_region_labels); if (cfg.repair_non_manifold_vertices) { while (auto candidate = find_problematic_vertex(tri, arr, walk_pl, source_face_ids, @@ -797,10 +795,10 @@ namespace roofer::reconstruction { forced_region_labels.push_back(repair_vertex( tri, *candidate, source_face_ids, cfg.manifold_repair_radius, cfg.manifold_height_tolerance)); - label_final_regions(tri, walk_pl, arr, source_face_ids, - forced_region_labels); } } + label_final_regions(tri, walk_pl, arr, source_face_ids, + forced_region_labels); // convert back from triangulation to arrangement // 1 recreate vertices and faces From 6005da69b30d158f79fa515aeb706ebe8a42727c Mon Sep 17 00:00:00 2001 From: ylannl Date: Fri, 19 Jun 2026 20:40:00 +0200 Subject: [PATCH 17/21] simplify final arrangement labeling in snapper --- .../reconstruction/ArrangementSnapper.cpp | 126 ++++++++---------- 1 file changed, 57 insertions(+), 69 deletions(-) diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index 7dd3dd2b..85217613 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -27,7 +27,6 @@ // #include #include #include -#include #include #include @@ -48,10 +47,7 @@ namespace roofer::reconstruction { typedef CGAL::Triangulation_vertex_base_with_info_2 VertexBaseWithInfo; typedef CGAL::Constrained_triangulation_face_base_2 FaceBase; - typedef CGAL::Triangulation_face_base_with_info_2 - FaceBaseWithInfo; - typedef CGAL::Triangulation_data_structure_2 + typedef CGAL::Triangulation_data_structure_2 TriangulationDataStructure; typedef CGAL::Constrained_Delaunay_triangulation_2< K, TriangulationDataStructure, Tag> @@ -142,57 +138,59 @@ namespace roofer::reconstruction { return regions; } - // labels the triangulation with region IDs based on the arrangement - void label_final_regions( - T& tri, - CGAL::Arr_walk_along_line_point_location& walk_pl, - Arrangement_2& source_arrangement, const SourceFaceIds& source_ids, - const std::vector& forced_labels) { - std::vector> forced_face_labels; + using ForcedFaceLabels = std::vector>; + + ForcedFaceLabels locate_forced_labels( + T& tri, const std::vector& forced_labels) { + ForcedFaceLabels forced_face_labels; for (const auto& forced : forced_labels) { auto located = tri.locate(forced.seed); if (!tri.is_infinite(located)) { forced_face_labels.emplace_back(located, forced.face_info); } } + return forced_face_labels; + } - for (auto& region : constrained_regions(tri)) { - // for this triangle region collect total overlap area for each - // arrangement face - std::unordered_map overlap_area; - for (auto face : region) { - auto centroid = CGAL::centroid(tri.triangle(face)); - if (auto* source = - source_face_at(centroid, walk_pl, source_arrangement)) { - overlap_area[source] += std::abs(tri.triangle(face).area()); - } + FaceInfo* select_region_label( + 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) { + // For this triangle region collect total overlap area for each source + // arrangement face. + std::unordered_map overlap_area; + for (auto face : region) { + auto centroid = CGAL::centroid(tri.triangle(face)); + if (auto* source = + source_face_at(centroid, walk_pl, source_arrangement)) { + overlap_area[source] += std::abs(tri.triangle(face).area()); } + } - // 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; - } + // 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; } + } - // apply the forced regions (these are determined upstream) - 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; - } + // Apply the forced regions determined by non-manifold repairs. + 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; } - - for (auto face : region) face->info() = selected; } + return selected; } // Heuristically compute clearance: the minimum distance to and opposing @@ -797,8 +795,6 @@ namespace roofer::reconstruction { cfg.manifold_height_tolerance)); } } - label_final_regions(tri, walk_pl, arr, source_face_ids, - forced_region_labels); // convert back from triangulation to arrangement // 1 recreate vertices and faces @@ -838,44 +834,36 @@ namespace roofer::reconstruction { typedef CGAL::Arr_walk_along_line_point_location Snap_walk_pl; Snap_walk_pl snap_walk_pl(arr_snap); - std::unordered_map> - output_face_labels; + std::unordered_set labelled_output_faces; + auto forced_face_labels = + locate_forced_labels(tri, forced_region_labels); for (const auto& region : constrained_regions(tri)) { - FaceInfo* region_label = nullptr; - double region_area = 0; - for (auto face : region) { - if (region_label == nullptr) region_label = face->info(); - region_area += std::abs(tri.triangle(face).area()); - } - if (region_label == nullptr || region_area == 0) continue; + FaceInfo* region_label = select_region_label( + tri, region, walk_pl, arr, source_face_ids, forced_face_labels); + if (region_label == nullptr) continue; auto centroid = CGAL::centroid(tri.triangle(region.front())); auto object = snap_walk_pl.locate( Arrangement_2::Point_2(centroid.x(), centroid.y())); - if (auto output_face = std::get_if(&object)) { - output_face_labels[arr_snap.non_const_handle(*output_face)] - [region_label] += region_area; + if (auto located_face = std::get_if(&object)) { + auto output_face = arr_snap.non_const_handle(*located_face); + if (output_face->is_unbounded()) continue; + if (!labelled_output_faces.insert(output_face).second) { + // throw roofer::rooferException( + // "Multiple snapped triangulation regions map to one " + // "arrangement face"); + } + output_face->data() = *region_label; } } arr_snap.unbounded_face()->data() = arr.unbounded_face()->data(); for (auto output_face : arr_snap.face_handles()) { if (output_face->is_unbounded()) continue; - const auto labels = output_face_labels.find(output_face); - if (labels == output_face_labels.end() || labels->second.empty()) { + if (!labelled_output_faces.contains(output_face)) { throw roofer::rooferException( "Unable to transfer snapped arrangement face label"); } - auto best = - std::max_element(labels->second.begin(), labels->second.end(), - [&](const auto& lhs, const auto& rhs) { - if (lhs.second != rhs.second) - return lhs.second < rhs.second; - return source_face_ids.at(lhs.first) > - source_face_ids.at(rhs.first); - }); - output_face->data() = *best->first; } // remove dangling edges if any, eg holes that collapse to a single edge From de5c65b31751ca502784088dbb0858ee383ed1ea Mon Sep 17 00:00:00 2001 From: ylannl Date: Mon, 22 Jun 2026 11:00:30 +0200 Subject: [PATCH 18/21] Do more rigorous repair clearance calculation + do not repair when lower than half repair radius --- .../reconstruction/ArrangementSnapper.hpp | 2 + .../reconstruction/ArrangementSnapper.cpp | 116 +++++++++++++++--- tests/test_arrangement_snapper.cpp | 33 +++++ 3 files changed, 133 insertions(+), 18 deletions(-) diff --git a/include/roofer/reconstruction/ArrangementSnapper.hpp b/include/roofer/reconstruction/ArrangementSnapper.hpp index a4d9aeda..e6b14d19 100644 --- a/include/roofer/reconstruction/ArrangementSnapper.hpp +++ b/include/roofer/reconstruction/ArrangementSnapper.hpp @@ -29,6 +29,8 @@ 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; }; diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index 85217613..9499da6c 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -193,23 +193,99 @@ namespace roofer::reconstruction { return selected; } - // Heuristically compute clearance: the minimum distance to and opposing - // edge from current vertex This is a conservative estimate, and a lower - // bound on the actual clearance - std::optional vertex_clearance(T& tri, Vertex_handle vertex) { - double clearance_sq = std::numeric_limits::max(); + double squared_distance_to_face(T& tri, Face_handle face, + const T::Point_2& point) { + auto triangle = tri.triangle(face); + if (triangle.has_on_boundary(point) || + triangle.has_on_bounded_side(point)) { + return 0; + } + + double distance_sq = std::numeric_limits::max(); + for (int i = 0; i < 3; ++i) { + distance_sq = + std::min(distance_sq, CGAL::to_double(CGAL::squared_distance( + point, tri.segment({face, i})))); + } + return distance_sq; + } + + std::vector get_repair_region(T& tri, Vertex_handle vertex, + double radius) { + std::vector region; + if (!(radius > 0)) return region; + + const double radius_sq = radius * radius; + const double tolerance = std::max(1.0, radius_sq) * 1e-12; + std::unordered_set visited; + std::queue queue; + + auto push_face = [&](Face_handle face) { + if (tri.is_infinite(face)) return; + if (visited.contains(face)) return; + if (squared_distance_to_face(tri, face, vertex->point()) > + radius_sq + tolerance) { + return; + } + visited.insert(face); + queue.push(face); + }; + Face_circulator face = tri.incident_faces(vertex), done(face); - if (face == nullptr) return std::nullopt; + if (face == nullptr) return region; do { - if (tri.is_infinite(face)) continue; - auto opposite_edge = Edge(face, face->index(vertex)); - clearance_sq = std::min( - clearance_sq, CGAL::to_double(CGAL::squared_distance( - vertex->point(), tri.segment(opposite_edge)))); + push_face(face); } while (++face != done); - if (!(clearance_sq > 0) || - clearance_sq == std::numeric_limits::max()) - return std::nullopt; + + while (!queue.empty()) { + auto current = queue.front(); + queue.pop(); + region.push_back(current); + + for (int i = 0; i < 3; ++i) { + push_face(current->neighbor(i)); + } + } + return region; + } + + bool edge_is_incident_to_vertex(T& tri, const Edge& edge, + Vertex_handle vertex) { + return edge.first->vertex(tri.cw(edge.second)) == vertex || + edge.first->vertex(tri.ccw(edge.second)) == vertex; + } + + // Compute the real local clearance up to the configured repair radius: + // grow over all triangulation faces touched by that radius, then measure + // only constrained edges from the grown patch. + std::optional calculate_vertex_clearance(T& tri, + Vertex_handle vertex, + double repair_radius) { + if (!(repair_radius > 0)) return std::nullopt; + + const double repair_radius_sq = repair_radius * repair_radius; + const double tolerance = std::max(1.0, repair_radius_sq) * 1e-12; + double clearance_sq = std::numeric_limits::max(); + for (auto face : get_repair_region(tri, vertex, repair_radius)) { + for (int i = 0; i < 3; ++i) { + Edge edge(face, i); + if (!tri.is_constrained(edge)) continue; + if (edge_is_incident_to_vertex(tri, edge, vertex)) continue; + + const double distance_sq = CGAL::to_double( + CGAL::squared_distance(vertex->point(), tri.segment(edge))); + if (distance_sq <= repair_radius_sq + tolerance) { + clearance_sq = std::min(clearance_sq, distance_sq); + } + } + } + + if (clearance_sq == std::numeric_limits::max()) { + // No constrained edge was found inside the search disk, so clearance + // should not reduce the configured maximum repair radius. + return std::numeric_limits::infinity(); + } + if (!(clearance_sq > 0)) return std::nullopt; return std::sqrt(clearance_sq); } @@ -346,14 +422,17 @@ namespace roofer::reconstruction { std::optional find_problematic_vertex( T& tri, Arrangement_2& source_arrangement, CGAL::Arr_walk_along_line_point_location& walk_pl, - const SourceFaceIds& source_ids, double height_tolerance) { + const SourceFaceIds& source_ids, double repair_radius, + double height_tolerance) { for (auto vertex = tri.finite_vertices_begin(); vertex != tri.finite_vertices_end(); ++vertex) { if (constrained_incident_edge_count(tri, vertex) < 4) continue; - auto clearance = vertex_clearance(tri, vertex); + auto clearance = calculate_vertex_clearance(tri, vertex, repair_radius); if (!clearance) continue; - if (*clearance <= 1e-6) continue; - auto sectors = incident_sectors(tri, vertex, *clearance, + if (*clearance < repair_radius / 2) continue; + const double sector_clearance = + std::isfinite(*clearance) ? *clearance : repair_radius; + auto sectors = incident_sectors(tri, vertex, sector_clearance, source_arrangement, walk_pl); if (sectors.size() < 4) continue; if (std::none_of(sectors.begin(), sectors.end(), @@ -789,6 +868,7 @@ namespace roofer::reconstruction { if (cfg.repair_non_manifold_vertices) { while (auto candidate = find_problematic_vertex(tri, arr, walk_pl, source_face_ids, + cfg.manifold_repair_radius, cfg.manifold_height_tolerance)) { forced_region_labels.push_back(repair_vertex( tri, *candidate, source_face_ids, cfg.manifold_repair_radius, diff --git a/tests/test_arrangement_snapper.cpp b/tests/test_arrangement_snapper.cpp index 2af69366..1dcb73c1 100644 --- a/tests/test_arrangement_snapper.cpp +++ b/tests/test_arrangement_snapper.cpp @@ -108,6 +108,13 @@ namespace { return arrangement; } + Arrangement_2 clearance_regression_arrangement() { + auto arrangement = cross_arrangement({15, 5, 15, 6}); + CGAL::insert(arrangement, Segment_2(Point_2(5.05, 5.55), Point_2(5.05, 9))); + CGAL::insert(arrangement, Segment_2(Point_2(5.55, 5.05), Point_2(9, 5.05))); + return arrangement; + } + Arrangement_2::Vertex_handle vertex_at(Arrangement_2& arrangement, const Point_2& point) { for (auto vertex : arrangement.vertex_handles()) { @@ -116,6 +123,18 @@ namespace { return {}; } + bool has_vertex_near(Arrangement_2& arrangement, const Point_2& point, + double tolerance) { + const double tolerance_sq = tolerance * tolerance; + for (auto vertex : arrangement.vertex_handles()) { + if (CGAL::to_double(CGAL::squared_distance(vertex->point(), point)) <= + tolerance_sq) { + return true; + } + } + return false; + } + void check_two_manifold_edges(const roofer::Mesh& mesh) { using Edge = std::pair; std::map edge_incidence; @@ -168,6 +187,20 @@ TEST_CASE("snapper repairs an alternating four-way roof junction") { check_two_manifold_edges(extruder->meshes.front()); } +TEST_CASE("snapper ignores non-constrained triangulation edges for clearance") { + auto arrangement = clearance_regression_arrangement(); + auto snapper = roofer::reconstruction::createArrangementSnapper(); + snapper->compute(arrangement, {.dist_thres = 0.001F, + .repair_non_manifold_vertices = true, + .manifold_repair_radius = 0.5F, + .manifold_height_tolerance = 1e-4F}); + + CHECK(vertex_at(arrangement, Point_2(5, 5)) == + Arrangement_2::Vertex_handle()); + CHECK(has_vertex_near(arrangement, Point_2(5.5, 5), 1e-9)); + CHECK(has_vertex_near(arrangement, Point_2(5, 5.5), 1e-9)); +} + TEST_CASE("snapper preserves a manifold four-way roof junction") { auto arrangement = cross_arrangement({0, 1, 2, 3}); auto snapper = roofer::reconstruction::createArrangementSnapper(); From 109bdb247ca07a9358f2b0bce18ecad8a2f74fa0 Mon Sep 17 00:00:00 2001 From: ylannl Date: Mon, 22 Jun 2026 12:16:16 +0200 Subject: [PATCH 19/21] perform dangling edge cleanup prior in cdt prior to manifold vertex repair --- .../reconstruction/ArrangementSnapper.cpp | 82 ++++++++++++++++--- tests/test_arrangement_snapper.cpp | 29 +++++++ 2 files changed, 100 insertions(+), 11 deletions(-) diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index 9499da6c..dd5184c7 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -299,6 +299,57 @@ namespace roofer::reconstruction { return count; } + 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; + 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)); + adjacency[first].insert(second); + adjacency[second].insert(first); + } + + std::queue leaves; + for (const auto& [vertex, neighbours] : adjacency) { + if (neighbours.size() == 1) leaves.push(vertex); + } + + while (!leaves.empty()) { + auto vertex = leaves.front(); + leaves.pop(); + + auto& neighbours = adjacency[vertex]; + if (neighbours.size() != 1) continue; + + auto neighbour = *neighbours.begin(); + neighbours.clear(); + + auto& neighbour_adjacency = adjacency[neighbour]; + neighbour_adjacency.erase(vertex); + if (neighbour_adjacency.size() == 1) leaves.push(neighbour); + + // Locate the edge from its stable endpoint handles instead of retaining + // a face handle across triangulation updates. + Face_handle face; + int index; + if (tri.is_edge(vertex, neighbour, face, index) && + tri.is_constrained({face, index})) { + tri.remove_constrained_edge(face, index); + } + } + + std::vector vertices_to_remove; + for (auto vertex = tri.finite_vertices_begin(); + vertex != tri.finite_vertices_end(); ++vertex) { + if (!tri.are_there_incident_constraints(vertex)) { + vertices_to_remove.push_back(vertex); + } + } + for (auto vertex : vertices_to_remove) tri.remove(vertex); + } + // calculate azimuth for each incident triangulation edge, and sample // corresponding face from arrangement std::vector incident_sectors( @@ -864,6 +915,12 @@ namespace roofer::reconstruction { } } + // Remove dangling constraint trees and the unconstrained vertices left + // behind by snapping. + remove_dangling_constraints_and_vertices(tri); + + // Detect and repair non-manifold vertices (ie. leading to a + // non-manifold edge during extrusion) and self-intersecting faces std::vector forced_region_labels; if (cfg.repair_non_manifold_vertices) { while (auto candidate = @@ -946,17 +1003,20 @@ namespace roofer::reconstruction { } } - // remove dangling edges if any, eg holes that collapse to a single edge - // after snapping - { - std::vector to_remove; - for (auto he : arr_snap.edge_handles()) { - if (he->face() == he->twin()->face()) to_remove.push_back(he); - } - for (auto he : to_remove) { - arr_snap.remove_edge(he); - } - } + // This should not be necessary, unless there is the above code still + // produces dangling edges (which it shoudldnt) + // // remove dangling edges if any, eg holes that collapse to a single + // edge + // // after snapping + // { + // std::vector to_remove; + // for (auto he : arr_snap.edge_handles()) { + // if (he->face() == he->twin()->face()) to_remove.push_back(he); + // } + // for (auto he : to_remove) { + // arr_snap.remove_edge(he); + // } + // } arr = arr_snap; } diff --git a/tests/test_arrangement_snapper.cpp b/tests/test_arrangement_snapper.cpp index 1dcb73c1..195d222c 100644 --- a/tests/test_arrangement_snapper.cpp +++ b/tests/test_arrangement_snapper.cpp @@ -187,6 +187,35 @@ TEST_CASE("snapper repairs an alternating four-way roof junction") { check_two_manifold_edges(extruder->meshes.front()); } +TEST_CASE("snapper removes dangling constraint chains before repair") { + auto expected = cross_arrangement({15, 5, 15, 6}); + auto arrangement = cross_arrangement({15, 5, 15, 6}); + CGAL::insert(arrangement, Segment_2(Point_2(5, 5), Point_2(6, 5.5))); + CGAL::insert(arrangement, Segment_2(Point_2(6, 5.5), Point_2(7, 6))); + + const roofer::reconstruction::ArrangementSnapperConfig config{ + .dist_thres = 0.001F, + .repair_non_manifold_vertices = true, + .manifold_repair_radius = 0.5F, + .manifold_height_tolerance = 1e-4F}; + auto snapper = roofer::reconstruction::createArrangementSnapper(); + snapper->compute(expected, config); + snapper->compute(arrangement, config); + + CHECK(arrangement.number_of_vertices() == expected.number_of_vertices()); + CHECK(arrangement.number_of_edges() == expected.number_of_edges()); + CHECK(arrangement.number_of_faces() == expected.number_of_faces()); + for (auto vertex : expected.vertex_handles()) { + CHECK(vertex_at(arrangement, vertex->point()) != + Arrangement_2::Vertex_handle()); + } + for (auto vertex : arrangement.vertex_handles()) { + CAPTURE(vertex->point()); + CHECK(vertex_at(expected, vertex->point()) != + Arrangement_2::Vertex_handle()); + } +} + TEST_CASE("snapper ignores non-constrained triangulation edges for clearance") { auto arrangement = clearance_regression_arrangement(); auto snapper = roofer::reconstruction::createArrangementSnapper(); From 7f7a86df761d9a93a89523ea0358a5fb95cd3f7b Mon Sep 17 00:00:00 2001 From: ylannl Date: Mon, 22 Jun 2026 12:27:01 +0200 Subject: [PATCH 20/21] cleanup comments --- .../reconstruction/ArrangementSnapper.cpp | 54 ------------------- 1 file changed, 54 deletions(-) diff --git a/src/core/reconstruction/ArrangementSnapper.cpp b/src/core/reconstruction/ArrangementSnapper.cpp index dd5184c7..05b1bb84 100644 --- a/src/core/reconstruction/ArrangementSnapper.cpp +++ b/src/core/reconstruction/ArrangementSnapper.cpp @@ -597,25 +597,6 @@ namespace roofer::reconstruction { return {*forced_seed, selected_face}; } - // tri_util::CDT triangulate_polygon(LinearRing& poly, float - // dupe_threshold_exp=3) { - // tri_util::CDT triangulation; - - // float dupe_threshold = (float) std::pow(10,-dupe_threshold_exp); - - // tri_util::insert_ring(poly, triangulation); - // for (auto& ring : poly.interior_rings()) { - // tri_util::insert_ring(ring, triangulation); - // } - - // if (triangulation.number_of_faces()==0) - // return triangulation; - - // mark_domains(triangulation); - - // return triangulation; - // } - void get_incident_constraints(T& tri, Vertex_handle vthis, Vertex_handle vexcept1, Vertex_handle vexcept2, @@ -705,41 +686,6 @@ namespace roofer::reconstruction { } } - // remove isolated vertices (sometimes these result from input - // arrangements with edge between 2 vertices with the same coordinates) - // { - // std::vector to_remove; - // for (auto vh = tri.finite_vertices_begin(); vh != - // tri.finite_vertices_end(); ++vh) { - // if (!tri.are_there_incident_constraints(vh)) { - // to_remove.push_back(vh); - // } - // } - // for (auto& v: to_remove) { - // std:: cout << "removing vertex without incident constraints...\n" - // tri.remove(v); - // } - // } - - // TriangleCollection triangles_og; - // vec1i segment_ids_og; - // for (auto fh = tri.finite_faces_begin(); fh != - // tri.finite_faces_end(); ++fh) { - // // only export triangles in the interior of a shape (thus excluding - // holes and exterior) - - // arr3f p0 = {float (fh->vertex(0)->point().x()), float - // (fh->vertex(0)->point().y()), 0}; arr3f p1 = {float - // (fh->vertex(1)->point().x()), float (fh->vertex(1)->point().y()), - // 0}; arr3f p2 = {float (fh->vertex(2)->point().x()), float - // (fh->vertex(2)->point().y()), 0}; triangles_og.push_back({ - // p0,p1,p2 }); segment_ids_og.push_back(fh->info()->segid); - // segment_ids_og.push_back(fh->info()->segid); - // segment_ids_og.push_back(fh->info()->segid); - // } - // output("triangles_og").set(triangles_og); - // output("segment_ids_og").set(segment_ids_og); - // Detect triangles with 3 short edges => collapse triangle to point // (remove 2 vertices) bool found_small_face; From bc0ac02dd79a12478aa76aad237fa14e1fd3fdbd Mon Sep 17 00:00:00 2001 From: ylannl Date: Mon, 22 Jun 2026 12:32:52 +0200 Subject: [PATCH 21/21] update changelog --- CHANGELOG.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 67615f5a..57c6d724 100644 --- a/CHANGELOG.md +++ b/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