Calculate box-width as orthogonal distance; aabb code for segment intersection (disabled)

2026-08-10 09:48:32 +00:00 · 2026-05-14 14:17:10 +02:00
parent 10f93545da
commit 3a14786a5b
1 changed files with 175 additions and 38 deletions
@@ -267,11 +267,12 @@ void clean_polygon(Polygon_2& poly) {

 void smooth_polygon(double factor, Polygon_2& poly) {
    auto ps = create_and_convert_offset_polygon(-factor, poly);
-    if (ps.size() == 1) {
-        auto r2 = ps.front();
-        ps = create_and_convert_offset_polygon(+factor, r2);
-        if (ps.size() == 1) {
-            poly = ps.front();
+    auto it = std::max_element(ps.begin(), ps.end(), [&](const auto& p, const auto& q) { return p.area() < q.area(); });
+    if (it != ps.end()) {
+        auto qs = create_and_convert_offset_polygon(+factor, *it);
+        auto jt = std::max_element(qs.begin(), qs.end(), [&](const auto& p, const auto& q) { return p.area() < q.area(); });
+        if (jt != qs.end()) {
+            poly = *jt;
        }
    }
 }
@@ -884,8 +885,7 @@ Polygon_with_holes_2 subdivide_polygon_on_same_input(SegmentLookup& segment_look
 std::tuple<
    std::map<Point_2, std::vector<Point_2>>, 
    std::map<Point_2, std::pair<Point_2, Point_2>>,
-    std::map<std::pair<Point_2, Point_2>, std::vector<const CGAL::Polygon_2<K>*>>,
-    std::map<Point_2, double>
+    std::map<std::pair<Point_2, Point_2>, std::vector<const CGAL::Polygon_2<K>*>>
 >
 build_line_graph(const std::vector<Polygon_2>& input_polygons, const std::map<Point_2, SegmentLookup::PolygonIt>& point_lookup, const std::vector<Polygon_2>& triangular_polygons)
 {
@@ -896,7 +896,9 @@ build_line_graph(const std::vector<Polygon_2>& input_polygons, const std::map<Po
    std::map<std::pair<Point_2, Point_2>, Point_2> segment_to_midpoint;
    std::map<Point_2, std::pair<Point_2, Point_2>> midpoint_to_segment;
    std::map<const CGAL::Polygon_2<K>*, std::vector<std::pair<Point_2, Point_2>>> facet_to_segment;
-    std::map<Point_2, double> midpoint_to_edge_length;
+
+
+    // std::map<Point_2, double> midpoint_to_edge_length;

    for (auto& tri : triangular_polygons) {
        for (size_t i = 0; i < 3; ++i) {
@@ -929,7 +931,7 @@ build_line_graph(const std::vector<Polygon_2>& input_polygons, const std::map<Po
        if (p1index->second != input_polygons.end() && p2index->second != input_polygons.end() && p1index->second != p2index->second) {
            segment_to_midpoint[p.first] = center;
            midpoint_to_segment[center] = p.first;
-            midpoint_to_edge_length[center] = std::sqrt(CGAL::to_double(CGAL::squared_distance(p.first.first, p.first.second)));
+            // midpoint_to_edge_length[center] = std::sqrt(CGAL::to_double(CGAL::squared_distance(p.first.first, p.first.second)));
        }
    }

@@ -949,7 +951,7 @@ build_line_graph(const std::vector<Polygon_2>& input_polygons, const std::map<Po
        }
    }

-    return {line_graph, midpoint_to_segment, segment_to_input_facet, midpoint_to_edge_length};
+    return {line_graph, midpoint_to_segment, segment_to_input_facet}; // } , midpoint_to_edge_length};
 }

 using DPoint = CGAL::Simple_cartesian<double>::Point_2;
@@ -958,8 +960,8 @@ using DBox = std::array<DPoint, 2>;

 struct CenterLineGraphData {
    std::vector<Point_2> points;
+    std::vector<std::optional<std::pair<Point_2, Point_2>>> orig_segments;
    std::vector<DPoint> points_double;
-    std::vector<double> widths;
    std::vector<std::pair<size_t, size_t>> edges;
    std::vector<std::vector<size_t>> incident_edges;
 };
@@ -1126,7 +1128,7 @@ bool obb_overlap(const T& a, const U& b, double eps = 1.e-9) {

 CenterLineGraphData make_center_line_graph_data(
    const std::map<Point_2, std::vector<Point_2>>& line_graph,
-    const std::map<Point_2, double>& midpoint_to_edge_length)
+    const std::map<Point_2, std::pair<Point_2, Point_2>>& midpoint_to_segment)
 {
    CenterLineGraphData graph;
    std::map<Point_2, size_t> point_to_index;
@@ -1139,9 +1141,13 @@ CenterLineGraphData make_center_line_graph_data(
        auto i = graph.points.size();
        point_to_index[p] = i;
        graph.points.push_back(p);
+        auto mit = midpoint_to_segment.find(p);
+        if (mit == midpoint_to_segment.end()) {
+            graph.orig_segments.emplace_back();
+        } else {
+            graph.orig_segments.emplace_back(mit->second);
+        }
        graph.points_double.push_back(to_double_point(p));
-        auto wt = midpoint_to_edge_length.find(p);
-        graph.widths.push_back(wt == midpoint_to_edge_length.end() ? 0. : wt->second);
        graph.incident_edges.emplace_back();
        return i;
    };
@@ -1175,7 +1181,41 @@ CenterLineGraphData make_center_line_graph_data(
 }

 double segment_width(const CenterLineGraphData& graph, const std::pair<size_t, size_t>& edge) {
-    return 0.5 * (graph.widths[edge.first] + graph.widths[edge.second]);
+    auto s1 = graph.orig_segments[edge.first];
+    auto s2 = graph.orig_segments[edge.second];
+    if (!s1 || !s2) {
+        throw std::runtime_error("!!!");
+    }
+
+    // A line segment between two points is expected to span a triangle, which means that one of the
+    // segment points ought to be shared.
+    Point_2 refpoint;
+    if (s1->first == s2->first) {
+        refpoint = s1->first;
+    } else if (s1->second == s2->first) {
+        refpoint = s1->second;
+    } else if (s1->first == s2->second) {
+        refpoint = s1->first;
+    } else if (s1->second == s2->second) {
+        refpoint = s1->second;
+    } else {
+        throw std::runtime_error("!!!!!");
+    }
+
+    auto p1 = graph.points_double[edge.first];
+    auto p2 = graph.points_double[edge.second];
+    auto v = p2 - p1;
+
+    if (v.squared_length() < 1.e-9) {
+        throw std::runtime_error("!!!!!!!");
+    }
+
+    v /= std::sqrt(v.squared_length());
+    auto n = perpendicular(v);
+    auto P = to_double_point(refpoint);
+    auto l = CGAL::abs((P - p1) * n);
+
+    return 2 * l;
 }

 bool edge_supports_same_line(
@@ -1266,6 +1306,7 @@ std::vector<LineRun> runs_from_graph(const CenterLineGraphData& graph, double an
            auto len = std::sqrt(d.squared_length());
            total_length += len;
            weighted_width_sum += len * segment_width(graph, edge);
+            // std::cout << " l: " << len << " w: " << segment_width(graph, edge) << " p1: " << graph.points_double[edge.first] << " p2: " << graph.points_double[edge.second] << std::endl;
        }

        auto run_direction = direction_sum.squared_length() < 1.e-18 ? ref : unit(direction_sum);
@@ -1288,6 +1329,8 @@ std::vector<LineRun> runs_from_graph(const CenterLineGraphData& graph, double an

        auto avg_width = total_length < 1.e-9 ? segment_width(graph, seed_edge) : weighted_width_sum / total_length;

+        // std::cout << "avg_width: " << avg_width << std::endl;
+
        runs.push_back({
            graph.points[start_index],
            graph.points[end_index],
@@ -1469,7 +1512,7 @@ std::vector<MergedBoxRecord> merge_intersecting_parallel_boxes_iterative(const s
        std::vector<size_t> members = clusters[i].members;
        members.insert(members.end(), clusters[j].members.begin(), clusters[j].members.end());
        auto merged = BoxCluster{members, merge_cluster_to_box(members, records)};
-        std::cout << "Result width: " << merged.box.avg_width << " fromt " << clusters[i].box.avg_width << " & " << clusters[j].box.avg_width << std::endl;
+        // std::cout << "Result width: " << merged.box.avg_width << "; from " << clusters[i].box.avg_width << " & " << clusters[j].box.avg_width << std::endl;

        std::vector<BoxCluster> next_clusters;
        next_clusters.reserve(clusters.size() - 1);
@@ -1549,6 +1592,7 @@ double point_to_oriented_box_distance(const DPoint& p, const MergedBoxRecord& bo
 }

 std::map<Point_2, std::vector<Point_2>> snap_points_to_box_axes(
+    DebugWriter& debug,
    const CenterLineGraphData& graph,
    const std::vector<MergedBoxRecord>& boxes,
    const K::FT& max_projection_distance) {
@@ -1592,15 +1636,18 @@ std::map<Point_2, std::vector<Point_2>> snap_points_to_box_axes(
            if (angle_between_dirs_deg(boxes[c1.box_index].direction, boxes[c2.box_index].direction) > 8.) {
                if (auto x = intersect_infinite_lines_exact(boxes[c1.box_index], boxes[c2.box_index])) {
                    snapped_points[i] = *x;
+                    debug.write_segment(graph.points[i], *x, "snap_candidate_1");
                    continue;
                }
            }
            snapped_points[i] = c1.projection;
+            debug.write_segment(graph.points[i], c1.projection, "snap_candidate_2");
            continue;
        }

        if (containing.size() == 1) {
            snapped_points[i] = containing[0].projection;
+            debug.write_segment(graph.points[i], containing[0].projection, "snap_candidate_3");
            continue;
        }

@@ -1613,6 +1660,7 @@ std::map<Point_2, std::vector<Point_2>> snap_points_to_box_axes(

        if ((graph.points[i] - best.projection).squared_length() < (max_projection_distance * max_projection_distance)) {
            snapped_points[i] = best.projection;
+            debug.write_segment(graph.points[i], best.projection, "snap_candidate_4");
        } else {
            snapped_points[i] = graph.points[i];
            std::cout << "Warning: snapping distance exceeding distance: " << std::sqrt(CGAL::to_double((snapped_points[i] - best.projection).squared_length())) << " > " << max_projection_distance << std::endl;
@@ -1640,9 +1688,9 @@ std::map<Point_2, std::vector<Point_2>> snap_points_to_box_axes(
 Graph2D<K> join_segment_runs(
    DebugWriter& debug,
    const std::map<Point_2, std::vector<Point_2>>& line_graph,
-    const std::map<Point_2, double>& midpoint_to_edge_length,
+    const std::map<Point_2, std::pair<Point_2, Point_2>>& midpoint_to_segment,
    const K::FT& max_projection_distance) {
-    auto graph = make_center_line_graph_data(line_graph, midpoint_to_edge_length);
+    auto graph = make_center_line_graph_data(line_graph, midpoint_to_segment);
    auto runs = runs_from_graph(graph);
    runs.erase(std::remove_if(runs.begin(), runs.end(), [](const LineRun& run) {
        return run.vertex_count <= 5;
@@ -1672,7 +1720,7 @@ Graph2D<K> join_segment_runs(
    }
    debug.write_polygons(run_polygons, "merged_boxes");

-    auto snapped_graph = snap_points_to_box_axes(graph, boxes, max_projection_distance);
+    auto snapped_graph = snap_points_to_box_axes(debug, graph, boxes, max_projection_distance);
    return Graph2D<K>(snapped_graph);
 }

@@ -2166,17 +2214,69 @@ std::list<std::pair<Point_2, Point_2>> extend_end_vertices_based_on_input(

 std::list<std::pair<Point_2, Point_2>>
 extend_end_vertices_based_on_input_simple(
+    DebugWriter& debug_output,
    const Graph2D<K>& G,
    const Polygon_list& outer_perimiter,
-    const K::FT& max_projection_distance)
+    const K::FT& max_projection_distance, int pass)
 {
    auto max_intersection_distance = max_projection_distance / 4;

+    using ValidationSegmentList = std::list<CGAL::Segment_3<K>>;
+    using ValidationSegmentIt = ValidationSegmentList::iterator;
+    using ValidationTreeTraits = CGAL::AABB_traits<K, CGAL::AABB_segment_primitive<K, ValidationSegmentIt>>;
+    using ValidationTree = CGAL::AABB_tree<ValidationTreeTraits>;
+
+    const auto& to_3d = [](const Point_2& p) {
+        return CGAL::Point_3<K>(p.x(), p.y(), 0);
+    };
+
+    const auto& to_2d = [](const CGAL::Point_3<K>& p) {
+        return CGAL::Point_2<K>(p.x(), p.y());
+    };
+
+    ValidationSegmentList validation_segments;
+    for (auto it = G.edges_begin(); it != G.edges_end(); ++it) {
+        if (it->first != it->second) {
+            validation_segments.emplace_back(to_3d(it->first), to_3d(it->second));
+        }
+    }
+
+    ValidationTree validation_tree(validation_segments.begin(), validation_segments.end());
+
+    const auto has_intersection = [&](const Segment_2& candidate) {
+        // @nb still disabled.
+        return false;
+        std::vector<ValidationSegmentIt> intersected_segments;
+        validation_tree.all_intersected_primitives(CGAL::Segment_3<K>(to_3d(candidate.source()), to_3d(candidate.target())), std::back_inserter(intersected_segments));
+
+        for (auto it : intersected_segments) {
+            auto existing = CGAL::Segment_2<K>(to_2d(it->source()), to_2d(it->target()));
+            auto intersection = CGAL::intersection(candidate, existing);
+            if (!intersection) {
+                continue;
+            }
+
+            if (auto* point = variant_get<Point_2>(&*intersection)) {
+                const bool candidate_endpoint = *point == candidate.source() || *point == candidate.target();
+                const bool existing_endpoint = *point == existing.source() || *point == existing.target();
+                if (candidate_endpoint && existing_endpoint) {
+                    continue;
+                }
+            }
+
+            return true;
+        }
+
+        return false;
+    };
+
    const auto& process_point = [&](const Point_2& M, const Point_2& incoming) {
+        bool within_any_perimeter = false;
        for (auto& bnd : outer_perimiter) {
            // if point M is contained in bnd interior:
            // if (!bnd.has_on_unbounded_side(M)) {
            if (bnd.has_on_bounded_side(M)) {
+                within_any_perimeter = true;
                // create ray incoming -> M
                CGAL::Ray_2<K> ray(incoming, M - incoming);

@@ -2192,9 +2292,13 @@ extend_end_vertices_based_on_input_simple(
                            auto dist = ((*xp) - M).squared_length();
                            if (dist < sq_distance_along_ray) {
                                if (dist < (max_intersection_distance * max_intersection_distance)) {
-                                    closest_segment = seg;
-                                    closest_intersection_point = *xp;
-                                    sq_distance_along_ray = dist;
+                                    if (has_intersection(CGAL::Segment_2<K>(M, *xp))) {
+                                        debug_output.write_segment(M, *xp, "exterior_extension_intersection");
+                                    } else {
+                                        closest_segment = seg;
+                                        closest_intersection_point = *xp;
+                                        sq_distance_along_ray = dist;
+                                    }
                                } else {
                                }
                            }
@@ -2218,8 +2322,12 @@ extend_end_vertices_based_on_input_simple(
                                auto d = CGAL::squared_distance(Pp, M);
                                if (d < (max_projection_distance * max_projection_distance)) {
                                    if (d < closest_distance) {
-                                        closest_distance = d;
-                                        closest_point = Pp;
+                                        if (has_intersection(CGAL::Segment_2<K>(M, Pp))) {
+                                            debug_output.write_segment(M, Pp, "exterior_projection_intersection");
+                                        } else {
+                                            closest_distance = d;
+                                            closest_point = Pp;
+                                        }
                                    }
                                }
                            }
@@ -2236,9 +2344,13 @@ extend_end_vertices_based_on_input_simple(
                                auto Pp = *it;
                                auto d = CGAL::squared_distance(Pp, M);
                                if (d < (max_projection_distance * max_projection_distance)) {
-                                    if (d < closest_distance) {
-                                        closest_distance = d;
-                                        closest_point = Pp;
+                                    if (has_intersection(CGAL::Segment_2<K>(M, Pp))) {
+                                        debug_output.write_segment(M, Pp, "exterior_nearby_intersection");
+                                    } else {
+                                        if (d < closest_distance) {
+                                            closest_distance = d;
+                                            closest_point = Pp;
+                                        }
                                    }
                                }
                            }
@@ -2246,13 +2358,19 @@ extend_end_vertices_based_on_input_simple(

                        if (closest_point) {
                            return closest_point;
-                            // constructed_segments.push_front({M, *closest_point});
                        } else {
                        }
                    }
                }
+            } else if (bnd.has_on_boundary(M)) {
+                return boost::optional<Point_2>{M};
            }
        }
+        if (within_any_perimeter) {
+            std::cout << "Within boundary but still no solution given" << std::endl;
+        } else {
+            std::cout << "Outside of all boundaries" << std::endl;
+        }
        return boost::optional<Point_2>{};
    };

@@ -2263,10 +2381,14 @@ extend_end_vertices_based_on_input_simple(
        if (it->second.size() == 1) {
            auto& M = it->first;
            if (auto result = process_point(M, *it->second.begin())) {
+                if (*result == M) {
+                    std::cout << "Point already on perimeter (" << M.x() << " " << M.y() << ")" << std::endl;
+                    continue;
+                }
                auto d = (M - *result).squared_length();
                solutions.emplace_back(d, *result, *it->second.begin());
            } else {
-                std::cout << "Unable to find projection or intersection point for interior boundary (" << M.x() << " " << M.y() << ")" << std::endl;
+                std::cout << "Unable to find projection or intersection point for interior boundary pass " << pass << " [round 1] (" << M.x() << " " << M.y() << ")" << std::endl;
            }
        }
    }
@@ -2277,8 +2399,15 @@ extend_end_vertices_based_on_input_simple(
    for (auto& [d, point, incoming] : solutions) {
        if (auto result = process_point(point, incoming)) {
            constructed_segments.push_front({point, *result});
+            debug_output.write_segment(point, *result, "exterior_constructed_segment");
+
+            auto d = CGAL::squared_distance(point, *result);
+            std::cout << "Distance: " << std::sqrt(CGAL::to_double(d)) << std::endl;
+            validation_segments.emplace_back(to_3d(point), to_3d(*result));
+            auto inserted_it = std::prev(validation_segments.end());
+            validation_tree.insert(inserted_it, validation_segments.end());
        } else {
-            std::cout << "Unable to find projection or intersection point for interior boundary (" << M.x() << " " << M.y() << ")" << std::endl;
+            std::cout << "Unable to find projection or intersection point for interior boundary pass " << pass << " [round 2] (" << point.x() << " " << point.y() << ")" << std::endl;
        }
    }

@@ -3245,6 +3374,14 @@ void arrange_cgal_polygons(svgfill::arrange_polygon_settings settings, const std
        std::swap(input_polygons, split_polygons);
    }

+    // before overlap elimition we can (and should) still smooth
+    /*
+    * @todo
+    for (auto& r : input_polygons) {
+        smooth_polygon(polygon_offset_distance / 100., r);
+    }
+    */
+
    t0.stop();
    t0 = timer.start("overlap elimination");

@@ -3398,7 +3535,7 @@ void arrange_cgal_polygons(svgfill::arrange_polygon_settings settings, const std

    debug_output.write_polygons(triangular_polygons, "triangulated_corridor");

-    auto [line_graph, midpoint_to_segment, segment_to_input_facet, midpoint_to_edge_length] = build_line_graph(input_polygons, point_lookup, triangular_polygons);
+    auto [line_graph, midpoint_to_segment, segment_to_input_facet] = build_line_graph(input_polygons, point_lookup, triangular_polygons);
    for (auto& p : line_graph) {
        for (auto& q : p.second) {
            debug_output.write_segment(p.first, q, "network_1");
@@ -3438,17 +3575,17 @@ void arrange_cgal_polygons(svgfill::arrange_polygon_settings settings, const std
        Graph2D<K> G2(line_graph);
        G = G2.weld_vertices();
        for (auto it = G.edges_begin(); it != G.edges_end(); ++it) {
-            debug_output.write_segment(it->first, it->second, "network_2");
+            debug_output.write_segment(it->first, it->second, "network_b_2");
        }
        eliminate_colinear_vertices(G);
        edge_slide(G);
        for (auto it = G.edges_begin(); it != G.edges_end(); ++it) {
-            debug_output.write_segment(it->first, it->second, "network_3");
+            debug_output.write_segment(it->first, it->second, "network_b_3");
        }
    };

    if (settings.line_cleaning_algo == 0) {
-        G = join_segment_runs(debug_output, line_graph, midpoint_to_edge_length, subdivision_length * 4);
+        G = join_segment_runs(debug_output, line_graph, midpoint_to_segment, subdivision_length * 4);
        Arrangement_2 arr;
        G.to_arrangement(arr);
        Graph2D<K> G2;
@@ -3456,7 +3593,7 @@ void arrange_cgal_polygons(svgfill::arrange_polygon_settings settings, const std
        eliminate_colinear_vertices(G2);
        G = G2;
        for (auto it = G.edges_begin(); it != G.edges_end(); ++it) {
-            debug_output.write_segment(it->first, it->second, "network_2");
+            debug_output.write_segment(it->first, it->second, "network_a_2");
        }
    } else {
        apply_line_cleaning_algo_1();
@@ -3470,8 +3607,8 @@ void arrange_cgal_polygons(svgfill::arrange_polygon_settings settings, const std
    bool fallback_to_line_cleaning_algo_1 = false;
    
    if (settings.line_cleaning_algo == 0) {
-        segments1 = extend_end_vertices_based_on_input_simple(G, outer_perimiter, subdivision_length * 16);
-        segments2 = extend_end_vertices_based_on_input_simple(G_orig, outer_perimiter, subdivision_length * 16);
+        segments1 = extend_end_vertices_based_on_input_simple(debug_output, G, outer_perimiter, subdivision_length * 16, 0);
+        segments2 = extend_end_vertices_based_on_input_simple(debug_output, G_orig, outer_perimiter, subdivision_length * 16, 1);
        
        Arrangement_2 arr_clean;
        G.to_arrangement(arr_clean);