Merge branch 'v0.8.0' into ifcmax/initial-refresh

.github/workflows/ci-lint.yaml

@@ -30,7 +30,7 @@ jobs:
           uv tool install ruff
           uv tool install black
           uv tool install poethepoet
-          uv tool install ty
+          uv tool install ty==0.0.34
 
       # black doesn't catch all syntax errors, so we check them explicitly.
       - name: Check syntax errors

.github/workflows/ci.yml

@@ -51,7 +51,7 @@ jobs:
       - name: Install dependencies
         run: |
           python -m pip install --upgrade pip
-          pip install xmlschema xsdata numpy lxml pytest isodate lark networkx tabulate python-dateutil shapely
+          pip install xmlschema xsdata numpy lxml pytest isodate lark networkx tabulate python-dateutil shapely pyparsing
           pip install src/bcf --no-deps
           pip install pytest-xdist==3.8.0
 

nix/build-all.py

@@ -126,9 +126,8 @@ ssl._create_default_https_context = ssl._create_unverified_context
 import time
 from collections.abc import Generator, Sequence
 from pathlib import Path
-from urllib.request import urlretrieve
-
 from typing import Literal, Union
+from urllib.request import urlretrieve
 
 logger = logging.getLogger(__name__)
 logger.setLevel(logging.INFO)

src/bonsai/bonsai/bim/module/drawing/helper.py

@@ -313,7 +313,7 @@ def format_distance(
             if not feet and not add_inches:
                 tx_dist += str(feet) + "'"
 
-            if not feet and add_inches:
+            if not feet and add_inches and unit_length != "INCHES":
                 if value < 0:
                     tx_dist += "-0' - "
                 else:

src/bonsai/bonsai/bim/module/model/wall.py

@@ -468,7 +468,11 @@ class ChangeExtrusionXAngle(bpy.types.Operator, tool.Ifc.Operator):
                     existing_x_angle = 0 if tool.Cad.is_x(existing_x_angle, 0, tolerance=0.001) else existing_x_angle
                     existing_x_angle = 0 if tool.Cad.is_x(existing_x_angle, pi, tolerance=0.001) else existing_x_angle
 
-                    profiles = extrusion.SweptArea.Profiles if extrusion.SweptArea.is_a("IfcCompositeProfileDef") else [extrusion.SweptArea]
+                    profiles = (
+                        extrusion.SweptArea.Profiles
+                        if extrusion.SweptArea.is_a("IfcCompositeProfileDef")
+                        else [extrusion.SweptArea]
+                    )
                     for profile in profiles:
                         coord_list = builder.get_polyline_coords(profile.OuterCurve)
                         coord_list = [

src/bonsai/test/tool/test_cost.py

@@ -17,20 +17,20 @@
 # along with Bonsai.  If not, see <http://www.gnu.org/licenses/>.
 
 
-import test.bim.bootstrap
 import ifcopenshell.api.cost
 
 import bonsai.core.tool
 import bonsai.tool as tool
 import test.bim.bootstrap
+from bonsai.tool.cost import Cost as subject
 from test.bim.bootstrap import NewFile
 
-from bonsai.tool.cost import Cost as subject
 
 class TestImplementsTool(NewFile):
     def test_run(self):
         assert isinstance(subject(), bonsai.core.tool.Cost)
 
+
 class TestDisableEditingCostItemParent(NewFile):
     def test_avoid_recursion_error(newfile, monkeypatch):
         class DummyProps:
@@ -39,11 +39,7 @@ class TestDisableEditingCostItemParent(NewFile):
                 self.active_cost_item_id = 5
 
         props = DummyProps()
-        monkeypatch.setattr(
+        monkeypatch.setattr("bonsai.tool.Cost.get_cost_props", lambda: props)
-            "bonsai.tool.Cost.get_cost_props",
-            lambda: props
-        )
         subject.disable_editing_cost_item_parent()
         assert props.active_cost_item_id == 0
         assert props.change_cost_item_parent is not False
-

src/ifcgeom/mapping/mapping.cpp

@@ -562,6 +562,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcMaterial* material) {
             }
             // Check if it's failed or just some unsupported case.
             if (failed_on_purpose_.find(styled_item) == failed_on_purpose_.end()) {
+                failed_on_purpose_.insert(material);
                 return nullptr;
             }
             Logger::Warning("Skipping unsupported material style for material: ", material);
@@ -569,6 +570,7 @@ taxonomy::ptr mapping::map_impl(const IfcSchema::IfcMaterial* material) {
     }
 
     // When material does not have a representation we don't create a style from it
+    failed_on_purpose_.insert(material);
     return nullptr;
 
     /*

src/ifcopenshell-python/ifcopenshell/draw.py

@@ -538,7 +538,10 @@ def main(
                     *(tup for i, tup in enumerate(zip(path_objects, section_polies, polies)) if has_relevant_zone(i))
                 )
 
-            arranged = W.arrange_polygons(*filter(None, (ARRANGE_POLYGON_SETTINGS,)), polies)
+            arranged = W.arrange_polygons(
+                *filter(None, (ARRANGE_POLYGON_SETTINGS,)),
+                polies,  # ty: ignore[too-many-positional-arguments]
+            )
             svg_data_3 = W.polygons_to_svg(arranged, False)
             dom3 = parseString(svg_data_3)
             svg3 = dom3.childNodes[0]

src/ifcopenshell-python/ifcopenshell/express/schema_class.py

@@ -420,7 +420,15 @@ class SchemaClass(codegen.Base):
 
             if isinstance(type, nodes.AggregationType):
                 aggr_type = type.aggregate_type
-                make_bound = lambda b: -1 if b == "?" else int(b)
+
+                def make_bound(b):
+                    # `?` and non-literal bounds (attribute references, arithmetic expressions) collapse to -1.
+                    #
+                    try:
+                        return int(b)
+                    except (TypeError, ValueError):
+                        return -1
+
                 bound1, bound2 = map(make_bound, (type.bounds.lower, type.bounds.upper))
                 decl_type = get_declared_type(type.type, emitted_names)
                 return x.aggregation_type(aggr_type, bound1, bound2, decl_type)
@@ -547,7 +555,16 @@ class SchemaClass(codegen.Base):
                 inv_attrs = []
                 for attr in type.inverse:
                     if attr.bounds:
-                        make_bound = lambda b: -1 if b == "?" else int(b)
+
+                        def make_bound(b):
+                            # `?` and non-literal bounds (attribute references, arithmetic
+                            # expressions) collapse to -1 (unbounded) — the C++ runtime has
+                            # no third state for "dynamic cardinality".
+                            try:
+                                return int(b)
+                            except (TypeError, ValueError):
+                                return -1
+
                         bound1, bound2 = map(make_bound, (attr.bounds.lower, attr.bounds.upper))
                     else:
                         bound1, bound2 = -1, -1

src/ifcopenshell-python/ifcopenshell/ifcopenshell_wrapper.pyi

@@ -1695,7 +1695,7 @@ class type_declaration(declaration):
 
 class uninitialized_tag: ...
 
-def arrange_polygons(polygons): ...
+def arrange_polygons(settings, polygons): ...
 def clear_schemas(): ...
 def construct_iterator(geometry_library, settings, file, num_threads): ...
 def construct_iterator_with_include_exclude(geometry_library, settings, file, elems, include, num_threads): ...

src/ifcopenshell-python/pyproject.toml

@@ -21,6 +21,7 @@ dependencies = [
     "isodate",
     "python-dateutil",
     "lark",
+    "pyparsing",
     "typing-extensions",
 ]
 

src/ifcopenshell-python/test/test_express_aggregate_bounds.py

@@ -0,0 +1,74 @@
+import os
+import sys
+import tempfile
+import unittest
+
+import ifcopenshell.express
+
+sys.path.insert(0, os.path.dirname(ifcopenshell.express.__file__))
+
+
+def _parse(schema_text):
+    with tempfile.NamedTemporaryFile(mode="w", suffix=".exp", delete=False) as f:
+        f.write(schema_text)
+        path = f.name
+    try:
+        return ifcopenshell.express.parse(path)
+    finally:
+        os.unlink(path)
+        cache = path + ".cache.dat"
+        if os.path.exists(cache):
+            os.unlink(cache)
+
+
+class TestAggregateBounds(unittest.TestCase):
+    def test_literal_bounds_preserved(self):
+        """After loading [1;3] -> (1, 3)?"""
+        s = _parse("SCHEMA t; ENTITY E; v : ARRAY [1:3] OF REAL; END_ENTITY; END_SCHEMA;")
+        agg = (
+            next(d for d in s.schema.declarations() if d.name() == "E")
+            .attributes()[0]
+            .type_of_attribute()
+            .as_aggregation_type()
+        )
+        self.assertEqual((agg.bound1(), agg.bound2()), (1, 3))
+        s.disown()
+
+    def test_unbounded_marker(self):
+        """[0:?] -> (0, -1)?"""
+        s = _parse("SCHEMA t; ENTITY E; v : LIST [0:?] OF REAL; END_ENTITY; END_SCHEMA;")
+        agg = (
+            next(d for d in s.schema.declarations() if d.name() == "E")
+            .attributes()[0]
+            .type_of_attribute()
+            .as_aggregation_type()
+        )
+        # import pdb; pdb.set_trace()
+        self.assertEqual((agg.bound1(), agg.bound2()), (0, -1))
+        s.disown()
+
+    def test_voxel_grid_with_dynamic_bound_loads(self):
+        """
+        Array that is an expression :  [1:dim_x*dim_y*dim_z]
+        Parsing must not crash, Bbund must be (1, -1)
+        """
+        s = _parse("""
+            SCHEMA t;
+            TYPE IfcBoolean = BOOLEAN; END_TYPE;
+
+            ENTITY IfcVoxelHolder;
+              NumberOfVoxelsX : INTEGER;
+              NumberOfVoxelsY : INTEGER;
+              NumberOfVoxelsZ : INTEGER;
+              Voxels : ARRAY [1:NumberOfVoxelsX*NumberOfVoxelsY*NumberOfVoxelsZ] OF IfcBoolean;
+            END_ENTITY;
+            END_SCHEMA;
+            """)
+        holder = next(d for d in s.schema.declarations() if d.name() == "IfcVoxelHolder")
+        voxels = holder.attributes()[-1].type_of_attribute().as_aggregation_type()
+        self.assertEqual((voxels.bound1(), voxels.bound2()), (1, -1))
+        s.disown()
+
+
+if __name__ == "__main__":
+    unittest.main()

src/ifcopenshell-python/test/util/test_cost.py

@@ -16,7 +16,6 @@
 # You should have received a copy of the GNU Lesser General Public License
 # along with IfcOpenShell.  If not, see <http://www.gnu.org/licenses/>.
 
-import pytest
 
 import ifcopenshell.api.control
 import ifcopenshell.api.cost
@@ -25,6 +24,7 @@ import ifcopenshell.api.root
 
 import ifcopenshell.util.cost as subject
 
+
 class TestGetCostItemForProduct(test.bootstrap.IFC4):
     def test_run(self):
         model = self.file
@@ -40,7 +40,7 @@ class TestGetCostItemForProduct(test.bootstrap.IFC4):
         cost_schedule = ifcopenshell.api.cost.add_cost_schedule(model)
         item1 = ifcopenshell.api.cost.add_cost_item(model, cost_schedule=cost_schedule)
         ifcopenshell.api.control.assign_control(model, related_objects=[element], relating_control=item1)
-        ifcopenshell.api.cost.remove_cost_item(model, cost_item = item1)
+        ifcopenshell.api.cost.remove_cost_item(model, cost_item=item1)
         assert list(subject.get_cost_items_for_product(element)) == []
 
     def test_no_assigned_cost_items(self):
@@ -49,4 +49,3 @@ class TestGetCostItemForProduct(test.bootstrap.IFC4):
         cost_schedule = ifcopenshell.api.cost.add_cost_schedule(model)
         item1 = ifcopenshell.api.cost.add_cost_item(model, cost_schedule=cost_schedule)
         assert list(subject.get_cost_items_for_product(element)) == []
-

src/ifcpatch/ifcpatch/recipes/ExtractElements.py

@@ -33,7 +33,7 @@ class Patcher(ifcpatch.BasePatcher):
         file: ifcopenshell.file,
         logger: Union[Logger, None] = None,
         query: str = "IfcWall",
-        assume_asset_uniqueness_by_name: bool = True,
+        assume_asset_uniqueness_by_name: bool = False,
     ):
         """Extract certain elements into a new model
 

src/svgfill/src/arrange_polygons.cpp

@@ -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 it = std::max_element(ps.begin(), ps.end(), [&](const auto& p, const auto& q) { return p.area() < q.area(); });
-        auto r2 = ps.front();
+    if (it != ps.end()) {
-        ps = create_and_convert_offset_polygon(+factor, r2);
+        auto qs = create_and_convert_offset_polygon(+factor, *it);
-        if (ps.size() == 1) {
+        auto jt = std::max_element(qs.begin(), qs.end(), [&](const auto& p, const auto& q) { return p.area() < q.area(); });
-            poly = ps.front();
+        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<std::pair<Point_2, Point_2>, std::vector<const CGAL::Polygon_2<K>*>>
-    std::map<Point_2, double>
 >
 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;
 };
@@ -1085,9 +1087,48 @@ bool aabb_overlap(const DBox& a, const DBox& b, double eps = 1.e-9) {
            a[1].y() + eps >= b[0].y();
 }
 
+std::pair<double, double> projected_interval_on_axis(const std::array<DPoint, 4>& points, const DDir& axis_u) {
+    auto u = unit(axis_u);
+    auto t0 = (points.front() - CGAL::ORIGIN) * u;
+    auto interval = std::make_pair(t0, t0);
+    for (auto& p : points) {
+        auto t = (p - CGAL::ORIGIN) * u;
+        interval.first = std::min(interval.first, t);
+        interval.second = std::max(interval.second, t);
+    }
+    return interval;
+}
+
+bool intervals_overlap(const std::pair<double, double>& a, const std::pair<double, double>& b, double eps = 1.e-9) {
+    return a.first <= b.second + eps && b.first <= a.second + eps;
+}
+
+bool obb_overlap(const std::array<DPoint, 4>& a, const std::array<DPoint, 4>& b, double eps = 1.e-9) {
+    auto has_separating_axis = [&](const std::array<DPoint, 4>& points) {
+        for (size_t i = 0; i < points.size(); ++i) {
+            auto edge = points[(i + 1) % points.size()] - points[i];
+            auto axis = unit(perpendicular(edge));
+            if (axis.squared_length() < 1.e-18) {
+                continue;
+            }
+            if (!intervals_overlap(projected_interval_on_axis(a, axis), projected_interval_on_axis(b, axis), eps)) {
+                return true;
+            }
+        }
+        return false;
+    };
+
+    return !has_separating_axis(a) && !has_separating_axis(b);
+}
+
+template <typename T, typename U>
+bool obb_overlap(const T& a, const U& b, double eps = 1.e-9) {
+    return obb_overlap(a.corners, b.corners, eps);
+}
+
 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;
@@ -1100,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;
     };
@@ -1136,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(
@@ -1227,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);
@@ -1249,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],
@@ -1375,9 +1457,19 @@ std::pair<double, double> merge_score(const MergedBoxRecord& a, const MergedBoxR
 }
 
 bool clusters_can_merge(const BoxCluster& a, const BoxCluster& b, double angle_tol_deg = 5., double axis_overlap_ratio_limit = 0.5) {
+    auto min_width = a.box.avg_width < b.box.avg_width ? a.box.avg_width : b.box.avg_width;
+    auto max_width = a.box.avg_width > b.box.avg_width ? a.box.avg_width : b.box.avg_width;
+    if (min_width > 1.e-9) {
+        if (max_width / min_width > 5) {
+            return false;
+        }
+    }
     if (!aabb_overlap(a.box.bbox, b.box.bbox)) {
         return false;
     }
+    if (!obb_overlap(a.box, b.box)) {
+        return false;
+    }
     if (angle_between_dirs_deg(a.box.direction, b.box.direction) > angle_tol_deg) {
         return false;
     }
@@ -1427,6 +1519,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 << "; from " << clusters[i].box.avg_width << " & " << clusters[j].box.avg_width << std::endl;
 
         std::vector<BoxCluster> next_clusters;
         next_clusters.reserve(clusters.size() - 1);
@@ -1506,6 +1599,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) {
@@ -1548,16 +1642,34 @@ std::map<Point_2, std::vector<Point_2>> snap_points_to_box_axes(
             auto& c2 = containing[1];
             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;
+                    auto seg = CGAL::Segment_2<K>(graph.points[i], *x);
-                    continue;
+                    bool intersects_with_other_box_axis = false;
+                    for (size_t j = 0; j < boxes.size(); ++j) {
+                        if (j == c1.box_index || j == c2.box_index) {
+                            continue;
+                        }
+                        auto& box = boxes[j];
+                        auto box_seg = CGAL::Segment_2<K>(box.exact_start, box.exact_end);
+                        if (CGAL::do_intersect(seg, box_seg)) {
+                            intersects_with_other_box_axis = true;
+                            break;
+                        }
+                    }
+                    if (!intersects_with_other_box_axis) {
+                        snapped_points[i] = *x;
+                        debug.write_segment(graph.points[i], *x, "snap_candidate_1");
+                        continue;
+                    }
                 }
             }
-            snapped_points[i] = c1.projection;
+            snapped_points[i] = (c1.projection - graph.points[i]).squared_length() < (c2.projection - graph.points[i]).squared_length() ? c1.projection : c2.projection;
+            debug.write_segment(graph.points[i], snapped_points[i], "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;
         }
 
@@ -1570,6 +1682,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;
@@ -1597,9 +1710,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;
@@ -1629,7 +1742,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);
 }
 
@@ -2123,61 +2236,139 @@ 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;
-    std::list<std::pair<Point_2, Point_2>> constructed_segments;
 
-        for (auto it = G.vertices_begin(); it != G.vertices_end(); ++it) {
+    using ValidationSegmentList = std::list<CGAL::Segment_3<K>>;
-        if (it->second.size() == 1) {
+    using ValidationSegmentIt = ValidationSegmentList::iterator;
-            auto& M = it->first;
+    using ValidationTreeTraits = CGAL::AABB_traits<K, CGAL::AABB_segment_primitive<K, ValidationSegmentIt>>;
+    using ValidationTree = CGAL::AABB_tree<ValidationTreeTraits>;
 
-            for (auto& bnd : outer_perimiter) {
+    const auto& to_3d = [](const Point_2& p) {
-                // if point M is contained in bnd interior:
+        return CGAL::Point_3<K>(p.x(), p.y(), 0);
-                // if (!bnd.has_on_unbounded_side(M)) {
+    };
-                if (bnd.has_on_bounded_side(M)) {
-                    auto& incoming = *it->second.begin();
-                    // create ray incoming -> M
-                    CGAL::Ray_2<K> ray(incoming, M - incoming);
 
-                    // intersect ray with boundary
+    const auto& to_2d = [](const CGAL::Point_3<K>& p) {
-                    boost::optional<CGAL::Segment_2<K>> closest_segment;
+        return CGAL::Point_2<K>(p.x(), p.y());
-                    boost::optional<CGAL::Point_2<K>> closest_intersection_point;
+    };
-                    K::FT sq_distance_along_ray = std::numeric_limits<double>::infinity();
+
-                    for (auto jt = bnd.edges_begin(); jt != bnd.edges_end(); ++jt) {
+    ValidationSegmentList validation_segments;
-                        const auto& seg = *jt;
+    for (auto it = G.edges_begin(); it != G.edges_end(); ++it) {
-                        auto x = CGAL::intersection(ray, seg);
+        if (it->first != it->second) {
-                        if (x) {
+            validation_segments.emplace_back(to_3d(it->first), to_3d(it->second));
-                            if (auto* xp = variant_get<CGAL::Point_2<K>>(&*x)) {
+        }
-                                auto dist = ((*xp) - M).squared_length();
+    }
-                                if (dist < sq_distance_along_ray) {
+
-                                    if (dist < (max_intersection_distance * max_intersection_distance)) {
+    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);
+
+                // intersect ray with boundary
+                boost::optional<CGAL::Segment_2<K>> closest_segment;
+                boost::optional<CGAL::Point_2<K>> closest_intersection_point;
+                K::FT sq_distance_along_ray = std::numeric_limits<double>::infinity();
+                for (auto jt = bnd.edges_begin(); jt != bnd.edges_end(); ++jt) {
+                    const auto& seg = *jt;
+                    auto x = CGAL::intersection(ray, seg);
+                    if (x) {
+                        if (auto* xp = variant_get<CGAL::Point_2<K>>(&*x)) {
+                            auto dist = ((*xp) - M).squared_length();
+                            if (dist < sq_distance_along_ray) {
+                                if (dist < (max_intersection_distance * max_intersection_distance)) {
+                                    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 {
+                                    }
+                                } else {
+                                }
+                            }
+                        }
+                    }
+                }
+
+                if (closest_intersection_point) {
+                    return closest_intersection_point;
+                    // constructed_segments.push_front({M, *closest_intersection_point});
+                } else {
+
+                    // Loop over boundary segments, and project point onto it, take the closest
+                    K::FT closest_distance = std::numeric_limits<double>::infinity();
+                    boost::optional<CGAL::Point_2<K>> closest_point;
+                    for (auto& poly : outer_perimiter) {
+                        for (auto jt = poly.edges_begin(); jt != poly.edges_end(); ++jt) {
+                            auto seg = *jt;
+                            auto Pp = seg.supporting_line().projection(M);
+                            if (seg.has_on(Pp)) {
+                                auto d = CGAL::squared_distance(Pp, M);
+                                if (d < (max_projection_distance * max_projection_distance)) {
+                                    if (d < closest_distance) {
+                                        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;
+                                        }
                                     }
                                 }
                             }
                         }
                     }
 
-                    if (closest_intersection_point) {
+                    if (closest_point) {
-                        constructed_segments.push_front({M, *closest_intersection_point});
+                        return closest_point;
+                        // constructed_segments.push_front({M, *closest_point});
                     } else {
 
-                        // Loop over boundary segments, and project point onto it, take the closest
-                        K::FT closest_distance = std::numeric_limits<double>::infinity();
-                        boost::optional<CGAL::Point_2<K>> closest_point;
                         for (auto& poly : outer_perimiter) {
-                            for (auto jt = poly.edges_begin(); jt != poly.edges_end(); ++jt) {
+                            for (auto it = poly.begin(); it != poly.end(); ++it) {
-                                auto seg = *jt;
+                                auto Pp = *it;
-                                auto Pp = seg.supporting_line().projection(M);
+                                auto d = CGAL::squared_distance(Pp, M);
-                                if (seg.has_on(Pp)) {
+                                if (d < (max_projection_distance * max_projection_distance)) {
-                                    auto d = CGAL::squared_distance(Pp, M);
+                                    if (has_intersection(CGAL::Segment_2<K>(M, Pp))) {
-                                    if (d < (max_projection_distance * max_projection_distance)) {
+                                        debug_output.write_segment(M, Pp, "exterior_nearby_intersection");
+                                    } else {
                                         if (d < closest_distance) {
                                             closest_distance = d;
                                             closest_point = Pp;
@@ -2188,32 +2379,58 @@ extend_end_vertices_based_on_input_simple(
                         }
 
                         if (closest_point) {
-                            constructed_segments.push_front({M, *closest_point});
+                            return closest_point;
                         } else {
-
-                            for (auto& poly : outer_perimiter) {
-                                for (auto it = poly.begin(); it != poly.end(); ++it) {
-                                    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 (closest_point) {
-                                constructed_segments.push_front({M, *closest_point});
-                            } else {
-                                std::cout << "Unable to find projection or intersection point for interior boundary (" << M.x() << " " << M.y() << ")" << std::endl;
-                            }
                         }
                     }
                 }
+            } 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>{};
+    };
+
+    using solution_length_point_incoming = std::tuple<K::FT, Point_2, Point_2>;
+    std::vector<solution_length_point_incoming> solutions;
+
+    for (auto it = G.vertices_begin(); it != G.vertices_end(); ++it) {
+        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, M, *it->second.begin());
+            } else {
+                std::cout << "Unable to find projection or intersection point for interior boundary pass " << pass << " [round 1] (" << M.x() << " " << M.y() << ")" << std::endl;
+            }
+        }
+    }
+
+    std::sort(solutions.begin(), solutions.end());
+    std::list<std::pair<Point_2, Point_2>> constructed_segments;
+
+    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 pass " << pass << " [round 2] (" << point.x() << " " << point.y() << ")" << std::endl;
+        }
     }
 
     return constructed_segments;
@@ -3179,6 +3396,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");
 
@@ -3332,7 +3557,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");
@@ -3372,17 +3597,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;
@@ -3390,7 +3615,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();
@@ -3404,8 +3629,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);
+        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(G_orig, outer_perimiter, subdivision_length * 16);
+        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);