From dffa3515c0a87d717b66ee0588d3087206de52bf Mon Sep 17 00:00:00 2001 From: Dion Moult Date: Sun, 8 Feb 2026 19:17:12 +1100 Subject: [PATCH 1/8] Reimplement sort / reverse / join function to format language, simplify text annotation variables, add tests Previously, sort, reverse list, and join functionality was implemented as special cases in Bonsai itself. Given that it has usecases (especially in material lists, but any sort of list applies) I've moved this function into the IOS formatting language. The IOS formatting language previously wasn't capable of this, but the awesome addition by @falken10vdl made the formatting language accept queries inline, so that means it can handle lists. I also added tests for all the new functions and expression syntax (+-*/ operators). I simplified the code that gets the evaluated text literal - previously it seems to call format() multiple times. --- src/bonsai/bonsai/bim/module/drawing/data.py | 22 +----------- src/bonsai/bonsai/tool/drawing.py | 13 +++---- .../ifcopenshell-python/selector_syntax.rst | 4 +++ .../ifcopenshell/util/selector.py | 33 +++++++++++------- .../test/util/test_selector.py | 34 ++++++++++++++++++- 5 files changed, 62 insertions(+), 44 deletions(-) diff --git a/src/bonsai/bonsai/bim/module/drawing/data.py b/src/bonsai/bonsai/bim/module/drawing/data.py index a4396d84c5..1668d386d3 100644 --- a/src/bonsai/bonsai/bim/module/drawing/data.py +++ b/src/bonsai/bonsai/bim/module/drawing/data.py @@ -365,15 +365,10 @@ class DecoratorData: for literal in literals: literal_value = literal.Literal - try: - eval_value = cls.evaluate_formatting_expressions(literal_value, product) - current_value = tool.Drawing.replace_text_literal_variables(eval_value, product) - except Exception: - current_value = literal_value literal_data = { "Literal": literal_value, "BoxAlignment": literal.BoxAlignment, - "CurrentValue": current_value, + "CurrentValue": tool.Drawing.replace_text_literal_variables(literal_value, product), } literals_data.append(literal_data) @@ -399,21 +394,6 @@ class DecoratorData: return element - @classmethod - def evaluate_formatting_expressions(cls, text: str, element=None) -> str: - """Evaluate formatting expressions wrapped in backticks using ifcopenshell.util.selector.format, always passing element context""" - import re - - def evaluate_expression(match): - try: - expression = match.group(1) - result = ifcopenshell.util.selector.format(expression, element) - return str(result) - except Exception as e: - return match.group(0) - - return re.sub(r"``([^`]+)``", evaluate_expression, text) - @classmethod def get_element_value_by_key(cls, element: ifcopenshell.entity_instance, key: str): """Get element value by its key using IfcOpenShell selector syntax""" diff --git a/src/bonsai/bonsai/tool/drawing.py b/src/bonsai/bonsai/tool/drawing.py index 3b98b92611..a7238c75dc 100644 --- a/src/bonsai/bonsai/tool/drawing.py +++ b/src/bonsai/bonsai/tool/drawing.py @@ -1306,9 +1306,7 @@ class Drawing(bonsai.core.tool.Drawing): element = tool.Ifc.get_entity(obj) assert element # updating text font size in EPset_Annotation.Classes - print("we got", font_size, repr(font_size)) font_size_str = next((key for key in FONT_SIZES if FONT_SIZES[key] == font_size), None) - print("so", font_size_str) classes = ifcopenshell.util.element.get_pset(element, "EPset_Annotation", "Classes") assert isinstance(classes, Union[str, None]) classes_split = classes.split() if classes else [] @@ -2101,16 +2099,13 @@ class Drawing(bonsai.core.tool.Drawing): if not product: return text - for command in re.findall("``.*?``", text): + for command in re.findall("``.+?``", text): original_command = command command_content = command[2:-2] - if command_content is None or str(command_content).strip().lower() == "none": + try: + text = text.replace(original_command, ifcopenshell.util.selector.format(command_content, product)) + except Exception: text = text.replace(original_command, "") - else: - try: - text = text.replace(original_command, ifcopenshell.util.selector.format(command_content, product)) - except Exception: - text = text.replace(original_command, "") for variable in re.findall("{{.*?}}", text): value = ifcopenshell.util.selector.get_element_value(product, variable[2:-2]) if isinstance(value, (list, tuple)): diff --git a/src/ifcopenshell-python/docs/ifcopenshell-python/selector_syntax.rst b/src/ifcopenshell-python/docs/ifcopenshell-python/selector_syntax.rst index 67f7b792b3..2d0cc9e071 100644 --- a/src/ifcopenshell-python/docs/ifcopenshell-python/selector_syntax.rst +++ b/src/ifcopenshell-python/docs/ifcopenshell-python/selector_syntax.rst @@ -252,6 +252,10 @@ nest formulas, for example ``concat(title("foo"), lower("Bar"))`` will produce "``number({{value}}[, {{decimal_separator}}[, {{thousands_separator}}]])``", "``number(1234.56, "","", ""."")``", "``1.234,56``", "Formats {{value}} with an optional custom {{decimal_separator}} and {{thousands_separator}}. The default separators are ``.`` and ``,``." "``metric_length({{value}}, {{precision}}, {{decimals}})``", "``metric_length(3.123, 0.1, 2)``", "``3.10``", "Rounds ``{{value}}`` to the nearest ``{{precision}}`` then displays using a certain amount of decimal places." "``imperial_length({{value}}, {{precision}}, {{input_unit}}, {{output_unit}}, {{suppress_zero_inches}})``", "``imperial_length(3.0, 4, ""foot"", ""foot"", true)`` OR ``imperial_length(3.0, 4, ""foot"", ""foot"", false)``", "``3'`` OR ``3' - 0""``", "The ``{{value}}`` may be specified either as ``foot`` or ``inch`` depending on ``{{input_unit}}``. The ``{{value}}`` is then rounded to the nearest ``1/{{precision}}`` inch, then formatted using fractional feet and inches if ``{{output_unit}}`` is set to ``foot``, or just inches if ``{{output_unit}}`` is set to ``inch``. When ``{{suppress_zero_inches}}`` is ``true`` (default), measurements with zero inches will omit the inch portion (e.g., ``3'`` instead of ``3' - 0""``)." + "``sort({{values}})``", "``sort({{mats.Name}})``", "``Name1, Name2``", "Sorts a list of items." + "``reverse({{values}})``", "``reverse({{mats.Name}})``", "``Name2, Name1``", "Reverses a list of items." + "``join({{separator}}, {{values}})``", "``join("-", {{mats.Name}})``", "``Name1-Name2``", "Joins a list of items with a custom separator. By default, all lists a rendered as comma separated." + "``{{value1}}[+-*/]{{value2}}``", "``{{z}}+3``", "``5``", "Does arithmetic. Typical operators such as +, -, \*, and / are allowed and can be mixed with other variables and formatting functions." When using queries in an IfcAnnotation tag surround with backticks. Examples: diff --git a/src/ifcopenshell-python/ifcopenshell/util/selector.py b/src/ifcopenshell-python/ifcopenshell/util/selector.py index 17051d6917..d3bef9758f 100644 --- a/src/ifcopenshell-python/ifcopenshell/util/selector.py +++ b/src/ifcopenshell-python/ifcopenshell/util/selector.py @@ -144,7 +144,7 @@ format_grammar = lark.Lark( | mul_div "*" function -> multiply | mul_div "/" function -> divide - function: round | number | int | format_length | lower | upper | title | concat | substr | variable | ESCAPED_STRING | SIGNED_NUMBER | "(" expression ")" + function: round | number | int | format_length | lower | upper | title | concat | substr | sort | reverse | join | variable | ESCAPED_STRING | SIGNED_NUMBER | "(" expression ")" variable: "{{" query_path "}}" query_path: /[^}]+/ @@ -160,6 +160,9 @@ format_grammar = lark.Lark( title: "title(" expression ")" concat: "concat(" expression ("," expression)* ")" substr: "substr(" expression "," SIGNED_INT ["," SIGNED_INT] ")" + sort: "sort(" expression ")" + reverse: "reverse(" expression ")" + join: "join(" ESCAPED_STRING "," expression ")" boolean: TRUE | FALSE TRUE: "true" | "True" | "TRUE" @@ -201,6 +204,8 @@ class FormatTransformer(lark.Transformer): self.element = element def start(self, args): + if isinstance(args[0], (list, tuple)): + return ", ".join(args[0]) return args[0] def expression(self, args): @@ -208,18 +213,11 @@ class FormatTransformer(lark.Transformer): def variable(self, args): """Handle variable substitution like {{z}} or {{Pset_Wall.FireRating}}""" - if self.element is None: - return "0" # Default value if no element context - - query_path = args[0] - try: - value = get_element_value(self.element, query_path) - if value is None: - return "0" - # Convert to string for further processing - return str(value) - except: - return "0" # Return default on error + if self.element: + try: + return get_element_value(self.element, args[0]) + except: + pass def query_path(self, args): """Extract the query path from variable""" @@ -301,6 +299,15 @@ class FormatTransformer(lark.Transformer): elif len(args) == 2: return str(args[0])[int(args[1]) :] + def sort(self, args): + return sorted(args[0]) + + def reverse(self, args): + return list(reversed(args[0])) + + def join(self, args): + return args[0].join(args[1]) + def boolean(self, args): if not args: return True diff --git a/src/ifcopenshell-python/test/util/test_selector.py b/src/ifcopenshell-python/test/util/test_selector.py index 959cacd02d..3c94876a66 100644 --- a/src/ifcopenshell-python/test/util/test_selector.py +++ b/src/ifcopenshell-python/test/util/test_selector.py @@ -35,7 +35,7 @@ import ifcopenshell.util.selector as subject import test.bootstrap -class TestFormat: +class TestFormat(test.bootstrap.IFC4): def test_no_formatting(self): assert subject.format("123") == "123" assert subject.format('"123"') == "123" @@ -78,6 +78,38 @@ class TestFormat: assert subject.format('imperial_length(3.0, 4, "foot", "foot", false)') == "3' - 0\"" assert subject.format('imperial_length(3.0, 4, "foot", "foot", False)') == "3' - 0\"" + def test_variable_formatting(self): + assert subject.format('{{undefined}}') is None + assert subject.format('upper({{undefined}})') == "NONE" + assert subject.format('int({{undefined}})') == "0" + element = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall") + assert subject.format('{{undefined}}', element) is None + assert subject.format('{{class}}', element) == "IfcWall" + assert subject.format('{{ class }}', element) == "IfcWall" + assert subject.format('upper({{ class }})', element) == "IFCWALL" + + def test_list_formatting(self): + element = ifcopenshell.api.root.create_entity(self.file, ifc_class="IfcWall") + material = ifcopenshell.api.material.add_material(self.file, name="CON01") + material2 = ifcopenshell.api.material.add_material(self.file, name="CON03") + material3 = ifcopenshell.api.material.add_material(self.file, name="CON02") + material_set = ifcopenshell.api.material.add_material_set(self.file, set_type="IfcMaterialLayerSet") + layer = ifcopenshell.api.material.add_layer(self.file, layer_set=material_set, material=material) + layer = ifcopenshell.api.material.add_layer(self.file, layer_set=material_set, material=material2) + layer = ifcopenshell.api.material.add_layer(self.file, layer_set=material_set, material=material3) + ifcopenshell.api.material.assign_material(self.file, products=[element], material=material_set) + assert subject.format('{{materials.Name}}', element) == "CON01, CON03, CON02" + assert subject.format('sort({{materials.Name}})', element) == "CON01, CON02, CON03" + assert subject.format('reverse({{materials.Name}})', element) == "CON02, CON03, CON01" + assert subject.format('join("-", {{materials.Name}})', element) == "CON01-CON03-CON02" + + def test_expressions(self): + assert subject.format('2+3') == "5" + assert subject.format('-2+3') == "1" + assert subject.format('2-3') == "-1" + assert subject.format('3*2') == "6" + assert subject.format('3/2') == "1.5" + class TestGetElementValue(test.bootstrap.IFC4): def test_selecting_an_elements_class_or_id_using_a_query(self): From 3fe9980aa41cfa90dc9e7d27b4b9bab4e30b4880 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Fri, 6 Feb 2026 22:45:56 +0000 Subject: [PATCH 2/8] Bump ruff from 0.14.14 to 0.15.0 Bumps [ruff](https://github.com/astral-sh/ruff) from 0.14.14 to 0.15.0. - [Release notes](https://github.com/astral-sh/ruff/releases) - [Changelog](https://github.com/astral-sh/ruff/blob/main/CHANGELOG.md) - [Commits](https://github.com/astral-sh/ruff/compare/0.14.14...0.15.0) --- updated-dependencies: - dependency-name: ruff dependency-version: 0.15.0 dependency-type: direct:production update-type: version-update:semver-minor ... Signed-off-by: dependabot[bot] --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 4e9d9687c5..74c3b0ee18 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -3,7 +3,7 @@ name = "IfcOpenShell" version = "0.0.0" dependencies = [ "black==26.1.0", - "ruff==0.14.14", + "ruff==0.15.0", "poethepoet", "gersemi==0.25.4", ] From 17b88fff229ebf14d93af3957488753fce3c50bb Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Fri, 6 Feb 2026 22:45:40 +0000 Subject: [PATCH 3/8] Bump aws-actions/configure-aws-credentials from 5 to 6 Bumps [aws-actions/configure-aws-credentials](https://github.com/aws-actions/configure-aws-credentials) from 5 to 6. - [Release notes](https://github.com/aws-actions/configure-aws-credentials/releases) - [Changelog](https://github.com/aws-actions/configure-aws-credentials/blob/main/CHANGELOG.md) - [Commits](https://github.com/aws-actions/configure-aws-credentials/compare/v5...v6) --- updated-dependencies: - dependency-name: aws-actions/configure-aws-credentials dependency-version: '6' dependency-type: direct:production update-type: version-update:semver-major ... Signed-off-by: dependabot[bot] --- .github/workflows/build_osx.yml | 2 +- .github/workflows/build_pyodide.yml | 2 +- .github/workflows/build_rocky.yml | 2 +- .github/workflows/build_rocky_arm.yml | 2 +- .github/workflows/build_win.yml | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/build_osx.yml b/.github/workflows/build_osx.yml index f53d91b98d..056cf13763 100644 --- a/.github/workflows/build_osx.yml +++ b/.github/workflows/build_osx.yml @@ -139,7 +139,7 @@ jobs: cd .. - name: Configure AWS credentials - uses: aws-actions/configure-aws-credentials@v5 + uses: aws-actions/configure-aws-credentials@v6 with: aws-access-key-id: ${{ secrets.AWS_UPLOAD_ACCESS_KEY_ID }} aws-secret-access-key: ${{ secrets.AWS_UPLOAD_SECRET_ACCESS_KEY }} diff --git a/.github/workflows/build_pyodide.yml b/.github/workflows/build_pyodide.yml index b5692b18f4..744bd9e772 100644 --- a/.github/workflows/build_pyodide.yml +++ b/.github/workflows/build_pyodide.yml @@ -77,7 +77,7 @@ jobs: git push || echo "Push failed" - name: Configure AWS credentials - uses: aws-actions/configure-aws-credentials@v5 + uses: aws-actions/configure-aws-credentials@v6 with: aws-access-key-id: ${{ secrets.AWS_UPLOAD_ACCESS_KEY_ID }} aws-secret-access-key: ${{ secrets.AWS_UPLOAD_SECRET_ACCESS_KEY }} diff --git a/.github/workflows/build_rocky.yml b/.github/workflows/build_rocky.yml index 6a21a9dfc3..662aa04f26 100644 --- a/.github/workflows/build_rocky.yml +++ b/.github/workflows/build_rocky.yml @@ -118,7 +118,7 @@ jobs: cd .. - name: Configure AWS credentials - uses: aws-actions/configure-aws-credentials@v5 + uses: aws-actions/configure-aws-credentials@v6 with: aws-access-key-id: ${{ secrets.AWS_UPLOAD_ACCESS_KEY_ID }} aws-secret-access-key: ${{ secrets.AWS_UPLOAD_SECRET_ACCESS_KEY }} diff --git a/.github/workflows/build_rocky_arm.yml b/.github/workflows/build_rocky_arm.yml index c595c63676..445f49658a 100644 --- a/.github/workflows/build_rocky_arm.yml +++ b/.github/workflows/build_rocky_arm.yml @@ -118,7 +118,7 @@ jobs: cd .. - name: Configure AWS credentials - uses: aws-actions/configure-aws-credentials@v5 + uses: aws-actions/configure-aws-credentials@v6 with: aws-access-key-id: ${{ secrets.AWS_UPLOAD_ACCESS_KEY_ID }} aws-secret-access-key: ${{ secrets.AWS_UPLOAD_SECRET_ACCESS_KEY }} diff --git a/.github/workflows/build_win.yml b/.github/workflows/build_win.yml index 6ef112da61..ab9d7b6f80 100644 --- a/.github/workflows/build_win.yml +++ b/.github/workflows/build_win.yml @@ -73,7 +73,7 @@ jobs: git push || echo "Push failed" - name: Configure AWS Credentials - uses: aws-actions/configure-aws-credentials@v5 + uses: aws-actions/configure-aws-credentials@v6 with: aws-access-key-id: ${{ secrets.AWS_UPLOAD_ACCESS_KEY_ID }} aws-secret-access-key: ${{ secrets.AWS_UPLOAD_SECRET_ACCESS_KEY }} From a4d5d4e19ad576ce237362802479401078ae7763 Mon Sep 17 00:00:00 2001 From: Thomas Krijnen Date: Mon, 9 Feb 2026 21:21:10 +0100 Subject: [PATCH 4/8] Ignore site placement also for site geometry #7654 --- src/ifcgeom/mapping/IfcObjectPlacement.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/ifcgeom/mapping/IfcObjectPlacement.cpp b/src/ifcgeom/mapping/IfcObjectPlacement.cpp index be0c222e08..098236d7c1 100644 --- a/src/ifcgeom/mapping/IfcObjectPlacement.cpp +++ b/src/ifcgeom/mapping/IfcObjectPlacement.cpp @@ -22,6 +22,18 @@ using namespace ifcopenshell::geometry; taxonomy::ptr mapping::map_impl(const IfcSchema::IfcObjectPlacement* inst) { + if (placement_rel_to_type_ || placement_rel_to_instance_) { + // @nb this is not a full solution because we only look for the direct PlacesObject relationships of the current placement, + // a more complete solution should track whether this element sits above the element of which the placement is being ignored. + auto self_places = inst->PlacesObject(); + for (auto iter = self_places->begin(); iter != self_places->end(); ++iter) { + if ((placement_rel_to_type_ && (*iter)->declaration().is(*placement_rel_to_type_)) || + (placement_rel_to_instance_ && (*iter)->as() == placement_rel_to_instance_)){ + return taxonomy::make(); + } + } + } + const IfcSchema::IfcObjectPlacement* relative_to = nullptr; const IfcUtil::IfcBaseInterface* transform; From a46cdbb907974c03be91bad002ce4f88aae682dd Mon Sep 17 00:00:00 2001 From: Thomas Krijnen Date: Mon, 9 Feb 2026 21:37:21 +0100 Subject: [PATCH 5/8] Ignore site placement also for site geometry - add queue #7654 --- src/ifcgeom/mapping/IfcObjectPlacement.cpp | 36 +++++++++++++++++----- 1 file changed, 29 insertions(+), 7 deletions(-) diff --git a/src/ifcgeom/mapping/IfcObjectPlacement.cpp b/src/ifcgeom/mapping/IfcObjectPlacement.cpp index 098236d7c1..a0a37dd30f 100644 --- a/src/ifcgeom/mapping/IfcObjectPlacement.cpp +++ b/src/ifcgeom/mapping/IfcObjectPlacement.cpp @@ -21,15 +21,37 @@ #define mapping POSTFIX_SCHEMA(mapping) using namespace ifcopenshell::geometry; +#include + taxonomy::ptr mapping::map_impl(const IfcSchema::IfcObjectPlacement* inst) { if (placement_rel_to_type_ || placement_rel_to_instance_) { - // @nb this is not a full solution because we only look for the direct PlacesObject relationships of the current placement, - // a more complete solution should track whether this element sits above the element of which the placement is being ignored. - auto self_places = inst->PlacesObject(); - for (auto iter = self_places->begin(); iter != self_places->end(); ++iter) { - if ((placement_rel_to_type_ && (*iter)->declaration().is(*placement_rel_to_type_)) || - (placement_rel_to_instance_ && (*iter)->as() == placement_rel_to_instance_)){ - return taxonomy::make(); + using QueueItem = std::pair; + std::deque q = {{inst, 0}}; + while (!q.empty()) { + auto [placement_entity, depth] = q.front(); + q.pop_front(); + + auto placement = placement_entity->as(); + if (!placement) { + continue; + } + + auto self_places = placement->PlacesObject(); + inst->ReferencedByPlacements(); + for (auto iter = self_places->begin(); iter != self_places->end(); ++iter) { + if ((placement_rel_to_type_ && (*iter)->declaration().is(*placement_rel_to_type_)) || + (placement_rel_to_instance_ && (*iter)->as() == placement_rel_to_instance_)) { + return taxonomy::make(); + } + } + + // Look for two levels deep, we want to know if we're at or *above* the + // element we're ignoring, but we don't want to traverse the entire model. + if (depth < 2) { + auto refs = placement->ReferencedByPlacements(); + for (auto& ref : *refs) { + q.emplace_back(ref, depth + 1); + } } } } From f5686fff2619578b5a3fb486cdb599f6e74e6566 Mon Sep 17 00:00:00 2001 From: Thomas Krijnen Date: Tue, 10 Feb 2026 09:25:02 +0100 Subject: [PATCH 6/8] Simplify destructor by removing null check #7650 --- src/ifcgeom/Converter.cpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/ifcgeom/Converter.cpp b/src/ifcgeom/Converter.cpp index 6b70ac7f51..a50a43cc66 100644 --- a/src/ifcgeom/Converter.cpp +++ b/src/ifcgeom/Converter.cpp @@ -12,11 +12,8 @@ ifcopenshell::geometry::Converter::Converter(std::unique_ptrsettings(); } -ifcopenshell::geometry::Converter::~Converter() -{ - if (mapping_ != nullptr) { - delete mapping_; - } +ifcopenshell::geometry::Converter::~Converter() { + delete mapping_; } namespace { From c6072e416cf3cc05ade4fa669e2afece292ee607 Mon Sep 17 00:00:00 2001 From: Thomas Krijnen Date: Tue, 10 Feb 2026 11:04:43 +0100 Subject: [PATCH 7/8] N-Section Lofting for Non-Polygonal (Curved) Shapes #7658 --- src/ifcgeom/kernels/opencascade/loft.cpp | 100 +++++++++++++++-------- 1 file changed, 65 insertions(+), 35 deletions(-) diff --git a/src/ifcgeom/kernels/opencascade/loft.cpp b/src/ifcgeom/kernels/opencascade/loft.cpp index 60997d24ad..74c9c187df 100644 --- a/src/ifcgeom/kernels/opencascade/loft.cpp +++ b/src/ifcgeom/kernels/opencascade/loft.cpp @@ -82,49 +82,79 @@ bool OpenCascadeKernel::convert(const taxonomy::loft::ptr loft, TopoDS_Shape& re } if (non_polygonal) { - if (loft->children.size() == 2) { - BRep_Builder BB; - TopoDS_Shell comp; - BB.MakeShell(comp); + if (loft->children.size() < 2) { + Logger::Error("Not enough sections to loft"); + return false; + } + std::vector> sections; + sections.reserve(loft->children.size()); - TopoDS_Shape f0, f1; - if (!convert(std::static_pointer_cast(loft->children.front()), f0) || - !convert(std::static_pointer_cast(loft->children.back()), f1)) - { + TopoDS_Shape f0, f1; + + // Convert all children to vectors of wires + for (const auto& child : loft->children) { + TopoDS_Shape shape; + if (!convert(std::static_pointer_cast(child), shape)) { + return false; + } + if (shape.ShapeType() != TopAbs_FACE) { + return false; + } + // At least make sure to have outer wire consistent, but in reality + // this is probably not a concern given how to build up these faces + auto f = TopoDS::Face(shape); + + if (child == loft->children.front()) { + f0 = f; + } else if (child == loft->children.back()) { + f1 = f; + } + + auto outer = BRepTools::OuterWire(f); + sections.emplace_back(); + sections.back().push_back(outer); + for (TopoDS_Iterator it(f); it.More(); it.Next()) { + if (outer != it.Value()) { + sections.back().push_back(TopoDS::Wire(it.Value())); + } + } + } + + auto first_wire_count = sections.front().size(); + for (auto& section : sections) { + if (section.size() != first_wire_count) { + Logger::Error("Inconsistent number of wires in sections"); return false; } - if (f0.ShapeType() != TopAbs_FACE || f1.ShapeType() != TopAbs_FACE) { + } + + BRep_Builder BB; + TopoDS_Shell comp; + BB.MakeShell(comp); + + for (size_t i = 0; i < first_wire_count; ++i) { + // Rule=True uses linear interpolation. + // This is critical for preventing twists in roads/railings. + BRepOffsetAPI_ThruSections builder(false, true); + for (auto& ws : sections) { + builder.AddWire(ws[i]); + } + builder.Build(); + if (!builder.IsDone()) { return false; } - - TopExp_Explorer exp1(f0, TopAbs_WIRE); - TopExp_Explorer exp2(f1, TopAbs_WIRE); - for (; exp1.More() && exp2.More(); exp1.Next(), exp2.Next()) { - const auto& w1 = TopoDS::Wire(exp1.Current()); - const auto& w2 = TopoDS::Wire(exp2.Current()); - BRepOffsetAPI_ThruSections builder; - builder.AddWire(w1); - builder.AddWire(w2); - builder.Build(); - if (!builder.IsDone()) { - return false; - } - for (TopExp_Explorer exp(builder.Shape(), TopAbs_FACE); exp.More(); exp.Next()) { - BB.Add(comp, exp.Current()); - } + for (TopExp_Explorer exp(builder.Shape(), TopAbs_FACE); exp.More(); exp.Next()) { + BB.Add(comp, exp.Current()); } - - BB.Add(comp, f0.Reversed()); - BB.Add(comp, f1); - - result = BRepBuilderAPI_MakeSolid(comp).Solid(); - - return true; - } else { - Logger::Error("Lofting more than two sections is not supported"); - return false; } + + BB.Add(comp, f0.Reversed()); + BB.Add(comp, f1); + + result = BRepBuilderAPI_MakeSolid(comp).Solid(); + + return true; } TopTools_ListOfShape faces; From e6780973da0d08771734aab5019998148706ef21 Mon Sep 17 00:00:00 2001 From: Thomas Krijnen Date: Tue, 10 Feb 2026 14:01:10 +0100 Subject: [PATCH 8/8] arrange_poly: Refactor into logical blocks; add timing --- src/svgfill/src/arrange_polygons.cpp | 1778 ++++++++++++-------------- src/svgfill/src/graph_2d.h | 12 + 2 files changed, 825 insertions(+), 965 deletions(-) diff --git a/src/svgfill/src/arrange_polygons.cpp b/src/svgfill/src/arrange_polygons.cpp index 7bddc2e41a..0271d76146 100644 --- a/src/svgfill/src/arrange_polygons.cpp +++ b/src/svgfill/src/arrange_polygons.cpp @@ -164,31 +164,6 @@ boost::optional subtract_retain_largest(const T& lhs, const T& rhs) { return boost::none; } -// Function to write polygons as line segments in OBJ format -void write_polygon_to_obj(std::ofstream& ofs, size_t& vertex_index, bool as_line, const Polygon_2& polygon, const std::string& name) { - ofs << "o " << name << "\n"; // Object name - - // Write vertices - for (auto vit = polygon.vertices_begin(); vit != polygon.vertices_end(); ++vit) { - ofs << "v " << CGAL::to_double(vit->x()) << " " << CGAL::to_double(vit->y()) << " 0\n"; - } - - if (as_line) { - // Write line segments (edges) - for (size_t j = 0; j < polygon.size(); ++j) { - ofs << "l " << vertex_index + j << " " << vertex_index + (j + 1) % polygon.size() << "\n"; - } - } else { - ofs << "f"; - for (size_t j = 0; j < polygon.size(); ++j) { - ofs << " " << vertex_index + j; - } - ofs << "\n"; - } - - vertex_index += polygon.size(); -} - Polygon_2 circ_to_poly(typename Arrangement_2::Ccb_halfedge_const_circulator circ) { Polygon_2 poly; @@ -208,27 +183,6 @@ Polygon_with_holes_2 circ_to_poly(typename Arrangement_2::Ccb_halfedge_const_cir return poly; } -void write_polygon_to_svg(std::ostream& ofs, const Polygon_2& polygon) { - ofs << "x()) << "," << CGAL::to_double(vit->y()) << " "; - } - ofs << "\" style=\"fill:none;stroke-width:1\" />\n"; -} - -// Function to write a Polygon_with_holes_2 to an SVG file -void write_polygon_with_holes_to_svg(std::ostream& ofs, const Polygon_with_holes_2& polygon_with_holes) { - // Write the outer boundary (main polygon) - if (!polygon_with_holes.is_unbounded()) { - write_polygon_to_svg(ofs, polygon_with_holes.outer_boundary()); - } - - // Write the holes (if any) with a different color (e.g., red) - for (auto hit = polygon_with_holes.holes_begin(); hit != polygon_with_holes.holes_end(); ++hit) { - write_polygon_to_svg(ofs, *hit); - } -} - Polygon_2 fuse_with_offset(const std::vector& polygons, double polygon_offset_distance) { // Find the outer perimeter using offset - union - negative offset std::vector offset_polygons; @@ -285,78 +239,77 @@ Polygon_2 fuse_with_offset(const std::vector& polygons, double polygo return inner_offset.front(); } -void arrange_cgal_polygons(const std::vector& input_polygons_, std::vector& output_polygons, double polygon_offset_distance = -1.) { - static const double OVERLAP_RESOLUTION_DISTANCE = 1.e-2; - // even larger amount of inset so that outer perimeter is safely within all input polygons even when overlap resolution is applied - // no, `1.e-2 + 1.e-5` creates issues with the outer perimeter, are there other tolerances in play? - static const double OUTER_PERIMITER_ADDITIONAL_INSET_AMOUNT = 1.e-5; +double estimate_polygon_offset_distance(const std::vector& polygons) { + double total_edge_length = 0.; + size_t num_edges = 0; + for (auto& p : polygons) { + for (auto it = p.edges_begin(); it != p.edges_end(); ++it) { + total_edge_length += std::sqrt(CGAL::to_double(CGAL::squared_distance(it->start(), it->end()))); + num_edges += 1; + } + } + return total_edge_length / num_edges / 2; +} - if (polygon_offset_distance < 0.) { - double total_edge_length = 0.; - size_t num_edges = 0; - for (auto& p : input_polygons_) { - for (auto it = p.edges_begin(); it != p.edges_end(); ++it) { - total_edge_length += std::sqrt(CGAL::to_double(CGAL::squared_distance(it->start(), it->end()))); - num_edges += 1; +void clean_polygon(Polygon_2& poly) { + // Ensure counterclockwise orientation and remove duplicate last point if present also remove close points + if (!poly.is_counterclockwise_oriented()) { + poly.reverse_orientation(); + } + std::vector> ps(poly.begin(), poly.end()); + if (ps.front() == ps.back()) { + ps.pop_back(); + } + poly = Polygon_2(ps.begin(), ps.end()); + remove_close_points(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(); + } + } +} + +template +void split_self_intersecting_polygon(const CGAL::Polygon_2& poly, OutIt output_it) { + if (poly.is_simple()) { + *output_it++ = poly; + return; + } + Arrangement_2 arr; + for (auto it = poly.edges_begin(); it != poly.edges_end(); ++it) { + CGAL::insert(arr, Segment_2(it->start(), it->end())); + } + for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it) { + if (it->is_unbounded()) { + for (auto jt = it->inner_ccbs_begin(); jt != it->inner_ccbs_end(); ++jt) { + auto inner = circ_to_poly(*jt); + // reverse because it's an inner bound to the infinite outer facet + inner.reverse_orientation(); + *output_it++ = inner; } } - polygon_offset_distance = total_edge_length / num_edges / 2; } +} - auto input_polygons__ = input_polygons_; - decltype(input_polygons__) input_polygons; - - for (auto& i : input_polygons__) { - std::vector> ps(i.begin(), i.end()); - if (ps.front() == ps.back()) { - ps.pop_back(); - } - input_polygons.emplace_back(ps.begin(), ps.end()); - } - - for (auto& polygon : input_polygons) { - if (!polygon.is_counterclockwise_oriented()) { - polygon.reverse_orientation(); - } - } - - for (auto& polygon : input_polygons) { - remove_close_points(polygon); - } - -#ifdef SVGFILL_DEBUG - std::ofstream obj("obj.obj"); - size_t vi = 1; - - std::ofstream svg("svg.svg"); - svg << "\n"; - - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { - write_polygon_to_obj(obj, vi, true, *it, "input_poly_" + std::to_string(std::distance(input_polygons.begin(), it))); - write_polygon_to_svg(svg, *it); - } - - obj << std::flush; -#endif - +std::set> +find_overlaps(const std::vector& polygons) { typedef CGAL::Box_intersection_d::Box_with_handle_d Box; + std::vector boxes; + std::vector>> input_triangulated; - std::vector boxes; - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { + for (auto it = polygons.begin(); it != polygons.end(); ++it) { constexpr double offset = 1.e-3; auto b = it->bbox(); boxes.emplace_back( CGAL::Bbox_2(b.xmin() - offset, b.ymin() - offset, b.xmax() + offset, b.ymax() + offset), - std::distance(input_polygons.begin(), it) - ); - - if (!it->is_simple()) { -#ifdef SVGFILL_DEBUG - write_polygon_to_obj(obj, vi, true, *it, "self-intersecting"); -#endif - throw std::runtime_error("Self-intersecting input"); - } + std::distance(polygons.begin(), it)); CGAL::Polygon_triangulation_decomposition_2 decompositor; std::vector temp; @@ -378,10 +331,8 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v bool registered_overlap = false; for (auto& t2 : input_triangulated[b.handle()]) { if (CGAL::squared_distance(t1, t2) < (1.e-3 * 1.e-3)) { - overlaps.insert({ - (a.handle() < b.handle()) ? a.handle() : b.handle(), - (a.handle() < b.handle()) ? b.handle() : a.handle() - }); + overlaps.insert({(a.handle() < b.handle()) ? a.handle() : b.handle(), + (a.handle() < b.handle()) ? b.handle() : a.handle()}); registered_overlap = true; break; } @@ -393,297 +344,319 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } }); - if (true) { - // solve overlaps by means of subtraction - // loop over overlaps and subtract the smaller polygon from the larger one + return overlaps; +} - std::set eliminated_polies; - std::map overlap_counts; - for (auto& p : overlaps) { - overlap_counts[p.first]++; - overlap_counts[p.second]++; +class DebugWriter { + public: + DebugWriter(bool enabled, const std::string& filename_prefix) + : enabled_(enabled) { + if (enabled_) { + obj.open(filename_prefix + ".obj"); + vi = 1; + svg.open(filename_prefix + ".svg"); + svg << "\n"; } - - for (const auto& edge : overlaps) { - // Skip eliminated - if (eliminated_polies.find(edge.first) != eliminated_polies.end() || - eliminated_polies.find(edge.second) != eliminated_polies.end()) { - continue; - } - - // Many overlaps indicate an aggregated polygon, skip them - /* - if (overlap_counts[edge.first] > 10 || overlap_counts[edge.second] > 10) { - if (overlap_counts[edge.first] > 10) { - eliminated_polies.insert(edge.first); - } - if (overlap_counts[edge.second] > 10) { - eliminated_polies.insert(edge.second); - } - continue; - } - */ - - // these are pointers now, because otherwise swap would not work? - auto* poly1 = &input_polygons[edge.first]; - auto* poly2 = &input_polygons[edge.second]; - - // Populate eliminated_polies with small polygons - // This can happen over time when modifications are made to the polygons to solve overlaps - bool skip = false; - if (poly1->area() < 1.e-2) { - eliminated_polies.insert(edge.first); - skip = true; - } - if (poly2->area() < 1.e-2) { - eliminated_polies.insert(edge.second); - skip = true; - } - // Small slivers are also just eliminated - if (!maybe_take_first_if_single_item(create_and_convert_offset_polygon(-1.e-1, *poly1))) { - eliminated_polies.insert(edge.first); - skip = true; - } - if (!maybe_take_first_if_single_item(create_and_convert_offset_polygon(-1.e-1, *poly2))) { - eliminated_polies.insert(edge.second); - skip = true; - } - if (skip) { - continue; - } - - // Skip polygons that have a very high intersection over union - // ratio, which indicates that they are very likely duplicates - if (CGAL::do_intersect(*poly1, *poly2)) { - std::vector result; - CGAL::intersection(*poly1, *poly2, std::back_inserter(result)); - typename K::FT intersection_area = 0; - for (auto& r : result) { - auto poly_area = r.outer_boundary().area(); - for (auto& h : r.holes()) { - poly_area -= h.area(); - } - intersection_area += poly_area; - } - CGAL::Polygon_with_holes_2 poly12; - CGAL::join(*poly1, *poly2, poly12); - typename K::FT union_area = poly12.outer_boundary().area(); - for (auto& h : poly12.holes()) { - union_area -= h.area(); - } - if (union_area > 0 && intersection_area / union_area > 0.99) { - // std::cerr << intersection_area / union_area << std::endl; - eliminated_polies.insert(edge.first); - continue; - } - } - - if (!(poly1->is_simple() && poly2->is_simple())) { - continue; - } - - { - std::vector result; - // std::cerr << poly1.area() << " " << poly2.area() << std::endl; - // std::cerr.flush(); - - boost::optional mp1, mp2, mp3, mp4; - bool swap = false; - - swap = poly1->area() <= poly2->area(); - if (swap) { - std::swap(poly1, poly2); - } - - bool success = false; - if ((mp1 = maybe_take_first_if_single_item(create_and_convert_offset_polygon(OVERLAP_RESOLUTION_DISTANCE, *poly2)))) { - if ((mp2 = subtract_retain_largest(*poly1, *mp1))) { - if ((mp3 = maybe_take_first_if_single_item(create_and_convert_offset_polygon(OVERLAP_RESOLUTION_DISTANCE * 2, *mp2)))) { - if ((mp4 = subtract_retain_largest(*poly2, *mp3))) { - *poly1 = *mp2; - *poly2 = *mp4; - success = true; - } - } - } - } - - /* - if (swap) { - // swap back to retain original ordering - // what's the point in swapping back here? - std::swap(poly1, poly2); - } - */ - - if (!success) { - eliminated_polies.insert(swap ? edge.first : edge.second); - continue; - } - } - } - - // iterate over the eliminated polygons and remove them from the input polygons - for (auto it = eliminated_polies.rbegin(); it != eliminated_polies.rend(); ++it) { - input_polygons.erase(input_polygons.begin() + *it); + } + ~DebugWriter() { + if (enabled_) { + svg << "\n"; + obj << std::flush; + obj.close(); + svg.close(); } } + void write_polygon(const Polygon_2& polygon, const std::string& name) { + if (enabled_) { + write_polygon_to_obj_(obj, vi, true, polygon, name); + write_polygon_to_svg_(svg, polygon, name); + obj << std::flush; + } + } + + void write_segment(const Point_2& p, const Point_2& q, const std::string& name) { + if (enabled_) { + if (last_segment_name_ != name) { + last_segment_name_ = name; + obj << "o " << name << "\n"; + } + obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; + obj << "v " << CGAL::to_double(q.x()) << " " << CGAL::to_double(q.y()) << " 0\n"; + obj << "l " << vi++; + obj << " " << vi++ << "\n"; + + svg << ""; + + obj << std::flush; + } + } + + void write_polygon(const Polygon_with_holes_2& polygon, const std::string& name) { + if (enabled_) { + write_polygon(polygon.outer_boundary(), name); + for (auto hit = polygon.holes_begin(); hit != polygon.holes_end(); ++hit) { + write_polygon(*hit, name); + } + } + } + + void write_polygons(const std::vector& polygons, const std::string& name) { + if (enabled_) { + size_t i = 0; + for (auto& polygon : polygons) { + write_polygon_to_obj_(obj, vi, true, polygon, name + "_" + std::to_string(i++)); + write_polygon_to_svg_(svg, polygon, name); + } + obj << std::flush; + } + } + + void write_polygons(const std::vector& polygons, const std::string& name) { + if (enabled_) { + size_t i = 0; + for (auto& polygon : polygons) { + write_polygon_to_obj_(obj, vi, true, polygon.outer_boundary(), name + "_" + std::to_string(i)); + write_polygon_to_svg_(svg, polygon.outer_boundary(), name); + for (auto hit = polygon.holes_begin(); hit != polygon.holes_end(); ++hit) { + write_polygon_to_obj_(obj, vi, true, *hit, name + "_" + std::to_string(i)); + write_polygon_to_svg_(svg, *hit, name); + } + } + obj << std::flush; + } + } + + private: + std::ofstream obj; + size_t vi; + std::ofstream svg; + bool enabled_; + std::string last_segment_name_; + + void write_polygon_to_svg_(std::ostream& ofs, const Polygon_2& polygon, const std::string& class_name = "") { + ofs << "x()) << "," << -CGAL::to_double(vit->y()) << " "; + } + ofs << "\"/>\n"; + } + + void write_polygon_to_obj_(std::ofstream& ofs, size_t& vertex_index, bool as_line, const Polygon_2& polygon, const std::string& name) { + ofs << "o " << name << "\n"; // Object name + + // Write vertices + for (auto vit = polygon.vertices_begin(); vit != polygon.vertices_end(); ++vit) { + ofs << "v " << CGAL::to_double(vit->x()) << " " << CGAL::to_double(vit->y()) << " 0\n"; + } + + if (as_line) { + // Write line segments (edges) + for (size_t j = 0; j < polygon.size(); ++j) { + ofs << "l " << vertex_index + j << " " << vertex_index + (j + 1) % polygon.size() << "\n"; + } + } else { + ofs << "f"; + for (size_t j = 0; j < polygon.size(); ++j) { + ofs << " " << vertex_index + j; + } + ofs << "\n"; + } + + vertex_index += polygon.size(); + } +}; + +void eliminate_overlaps(double OVERLAP_RESOLUTION_DISTANCE, std::vector& polygons) { + // solve overlaps by means of subtraction + // loop over overlaps and subtract the smaller polygon from the larger one + + std::set eliminated_polies; + /* - if constexpr (false) { - // solve overlap by means of union into components - std::vector> adj(input_polygons.size()); - for (const auto& edge : overlaps) { - adj[edge.first].push_back(edge.second); - adj[edge.second].push_back(edge.first); - } - - std::vector visited(input_polygons.size(), false); - std::vector> connected_components; - - for (size_t v = 0; v < input_polygons.size(); ++v) { - if (!visited[v]) { - connected_components.emplace_back(); - - std::stack stack; - stack.push(v); - visited[v] = true; - - while (!stack.empty()) { - size_t u = stack.top(); - stack.pop(); - connected_components.back().push_back(u); - - for (size_t neighbor : adj[u]) { - if (!visited[neighbor]) { - visited[neighbor] = true; - stack.push(neighbor); - } - } - } - } - } - - std::vector fused_polies; - - for (auto& comp : connected_components) { - std::vector comp_polies; - if (comp.size() == 1) { - fused_polies.push_back(input_polygons[comp.front()]); - } else { - for (auto& c : comp) { - comp_polies.push_back(input_polygons[c]); - } - fused_polies.push_back(fuse_with_offset(comp_polies, 1.e-2)); - } - } - -#ifdef SVGFILL_DEBUG - for (auto it = fused_polies.begin(); it != fused_polies.end(); ++it) { - write_polygon_to_obj(obj, vi, true, *it, "fused_poly_" + std::to_string(std::distance(fused_polies.begin(), it))); - write_polygon_to_svg(svg, *it); - } -#endif - - input_polygons = fused_polies; + std::map overlap_counts; + for (auto& p : overlaps) { + overlap_counts[p.first]++; + overlap_counts[p.second]++; } */ - { - // [NB Nov 6] we cannot do this anymore because it could revert the spacing between input polygons - // that touch in the corner. - // Now that overlaps/touches at corners are handled more locally only a small indent is produced - // which would be undone by means of an inset+offset. - // - // [NB Nov 10] this is actually still necessary though, but we apply a much smaller distance now - // to keep the overlap eliminations in tact - // - // Inset-offset to remove tiny details that may cause enourmous spikes in offsets - for (auto& r : input_polygons) { - auto ps = create_and_convert_offset_polygon(-polygon_offset_distance / 10000., r); - if (ps.size() == 1) { - auto r2 = ps.front(); - ps = create_and_convert_offset_polygon(+polygon_offset_distance / 10000., r2); - if (ps.size() == 1) { - r = ps.front(); + auto overlaps = find_overlaps(polygons); + + for (const auto& edge : overlaps) { + // Skip eliminated + if (eliminated_polies.find(edge.first) != eliminated_polies.end() || + eliminated_polies.find(edge.second) != eliminated_polies.end()) { + continue; + } + + // Many overlaps indicate an aggregated polygon, skip them + /* + if (overlap_counts[edge.first] > 10 || overlap_counts[edge.second] > 10) { + if (overlap_counts[edge.first] > 10) { + eliminated_polies.insert(edge.first); + } + if (overlap_counts[edge.second] > 10) { + eliminated_polies.insert(edge.second); + } + continue; + } + */ + + // these are pointers now, because otherwise swap would not work? + auto* poly1 = &polygons[edge.first]; + auto* poly2 = &polygons[edge.second]; + + // @todo this is applied during overlap processing, maybe better after the boolean operation, + // because they can be come small or narrow when overlaps are resolved + + // Populate eliminated_polies with small polygons + // This can happen over time when modifications are made to the polygons to solve overlaps + bool skip = false; + if (poly1->area() < 1.e-2) { + eliminated_polies.insert(edge.first); + skip = true; + } + if (poly2->area() < 1.e-2) { + eliminated_polies.insert(edge.second); + skip = true; + } + // Small slivers are also just eliminated + if (!maybe_take_first_if_single_item(create_and_convert_offset_polygon(-1.e-1, *poly1))) { + eliminated_polies.insert(edge.first); + skip = true; + } + if (!maybe_take_first_if_single_item(create_and_convert_offset_polygon(-1.e-1, *poly2))) { + eliminated_polies.insert(edge.second); + skip = true; + } + if (skip) { + continue; + } + + // Skip polygons that have a very high intersection over union + // ratio, which indicates that they are very likely duplicates + if (CGAL::do_intersect(*poly1, *poly2)) { + std::vector result; + CGAL::intersection(*poly1, *poly2, std::back_inserter(result)); + typename K::FT intersection_area = 0; + for (auto& r : result) { + auto poly_area = r.outer_boundary().area(); + for (auto& h : r.holes()) { + poly_area -= h.area(); } + intersection_area += poly_area; + } + CGAL::Polygon_with_holes_2 poly12; + CGAL::join(*poly1, *poly2, poly12); + typename K::FT union_area = poly12.outer_boundary().area(); + for (auto& h : poly12.holes()) { + union_area -= h.area(); + } + if (union_area > 0 && intersection_area / union_area > 0.99) { + // std::cerr << intersection_area / union_area << std::endl; + eliminated_polies.insert(edge.first); + continue; + } + } + + if (!(poly1->is_simple() && poly2->is_simple())) { + continue; + } + + { + std::vector result; + + boost::optional mp1, mp2, mp3, mp4; + bool swap = false; + + swap = poly1->area() <= poly2->area(); + if (swap) { + std::swap(poly1, poly2); + } + + bool success = false; + if ((mp1 = maybe_take_first_if_single_item(create_and_convert_offset_polygon(OVERLAP_RESOLUTION_DISTANCE, *poly2)))) { + if ((mp2 = subtract_retain_largest(*poly1, *mp1))) { + if ((mp3 = maybe_take_first_if_single_item(create_and_convert_offset_polygon(OVERLAP_RESOLUTION_DISTANCE * 2, *mp2)))) { + if ((mp4 = subtract_retain_largest(*poly2, *mp3))) { + *poly1 = *mp2; + *poly2 = *mp4; + success = true; + } + } + } + } + + if (!success) { + eliminated_polies.insert(swap ? edge.first : edge.second); + continue; } } } -#ifdef SVGFILL_DEBUG - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { - write_polygon_to_obj(obj, vi, true, *it, "processed_input_poly_" + std::to_string(std::distance(input_polygons.begin(), it))); - write_polygon_to_svg(svg, *it); + // iterate over the eliminated polygons and remove them from the input polygons + for (auto it = eliminated_polies.rbegin(); it != eliminated_polies.rend(); ++it) { + polygons.erase(polygons.begin() + *it); } -#endif +} - // Unfortunately CGAL does not seem to have a ready to use aabb primitive for segments in 2D, - // so we have to use 3D segments and aabb tree for 2D polygons. - std::list> all_segs; - std::unordered_map*, decltype(input_polygons.begin())> seg_to_poly; +class SegmentLookup { + public: + typedef std::vector::const_iterator PolygonIt; - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { - for (auto eit = it->edges_begin(); eit != it->edges_end(); ++eit) { - CGAL::Segment_3 seg3d( - CGAL::Point_3(eit->source().x(), eit->source().y(), 0), - CGAL::Point_3(eit->target().x(), eit->target().y(), 0) - ); - all_segs.push_back(seg3d); - seg_to_poly[&all_segs.back()] = it; - } - } - - using TreeTraits = CGAL::AABB_traits>::iterator>>; - using Tree = CGAL::AABB_tree; - - Tree tree(all_segs.begin(), all_segs.end()); - tree.accelerate_distance_queries(); - - auto input_polygon_boundary = - [&](const Point_2& p, double tol = 1e-5) -> decltype(input_polygons.begin()) + SegmentLookup(const std::vector& polygons) + : polygons_ref_(polygons) { + // Unfortunately CGAL does not seem to have a ready to use aabb primitive for segments in 2D, + // so we have to use 3D segments and aabb tree for 2D polygons. + for (auto it = polygons.begin(); it != polygons.end(); ++it) { + for (auto eit = it->edges_begin(); eit != it->edges_end(); ++eit) { + CGAL::Segment_3 seg3d( + CGAL::Point_3(eit->source().x(), eit->source().y(), 0), + CGAL::Point_3(eit->target().x(), eit->target().y(), 0)); + all_segs.push_back(seg3d); + seg_to_poly[&all_segs.back()] = it; + } + } + tree_ = Tree(all_segs.begin(), all_segs.end()); + tree_.accelerate_distance_queries(); + } + + // This part is the most computationally expensive. Caching effectively halves the lookup time here, since every vertex on the subdivided corridor mesh has on average two outgoing edges. + PolygonIt input_polygon_boundary(const Point_2& p, double tol = 1e-5) { + auto it = input_polygon_boundary_cache_.find(p); + if (it != input_polygon_boundary_cache_.end()) { + return it->second; + } + // Find closest point & corresponding segment - auto closest = tree.closest_point_and_primitive(CGAL::Point_3(p.x(), p.y(), 0)); + auto closest = tree_.closest_point_and_primitive(CGAL::Point_3(p.x(), p.y(), 0)); const auto& closest_pt = closest.first; auto seg_ptr = &*closest.second; double d = CGAL::to_double(CGAL::squared_distance(p, Point_2(closest_pt.x(), closest_pt.y()))); + + PolygonIt res; if (d < (tol * tol)) { - return seg_to_poly.find(seg_ptr)->second; + res = seg_to_poly.find(seg_ptr)->second; + } else { + res = polygons_ref_.end(); } - return input_polygons.end(); + + input_polygon_boundary_cache_[p] = res; + return res; }; - /* - auto input_polygon_boundary = [&input_polygons](const CGAL::Point_2& p, double tol = 1.e-5) { - // unfortunately some imprecision slept into the code so we can't - // so we can't just use has_on_boundary() anymore - double D = std::numeric_limits::infinity(); - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { - for (auto jt = it->edges_begin(); jt != it->edges_end(); ++jt) { - const auto& seg = *jt; - auto d = std::sqrt(CGAL::to_double(CGAL::squared_distance(seg, p))); - if (d < D) { - D = d; - } - if (d < tol) { - return it; - } - } - } - return input_polygons.end(); - }; - */ - - auto close_input_point = [&input_polygons](const CGAL::Point_2& P) { + std::pair> close_input_point(const CGAL::Point_2& P) const { + // @todo use tree CGAL::Point_2 closest; double closest_distance = std::numeric_limits::infinity(); - auto input_it = input_polygons.end(); + auto input_it = polygons_ref_.end(); // unfortunately some imprecision slept into the code so we can't // so we can't just use has_on_boundary() anymore - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { + for (auto it = polygons_ref_.begin(); it != polygons_ref_.end(); ++it) { for (auto& p : *it) { auto d = std::sqrt(CGAL::to_double(CGAL::squared_distance(P, p))); if (d < closest_distance) { @@ -697,14 +670,16 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v return std::make_pair(input_it, closest); }; - auto project_input_point = [&input_polygons](const CGAL::Point_2& P) { + std::pair> project_input_point(const CGAL::Point_2& P) const { + // @todo use tree + CGAL::Point_2 closest; typename K::FT closest_sq_distance = std::numeric_limits::infinity(); - auto input_it = input_polygons.end(); + auto input_it = polygons_ref_.end(); // unfortunately some imprecision slept into the code so we can't // so we can't just use has_on_boundary() anymore - for (auto it = input_polygons.begin(); it != input_polygons.end(); ++it) { + for (auto it = polygons_ref_.begin(); it != polygons_ref_.end(); ++it) { for (auto jt = it->edges_begin(); jt != it->edges_end(); ++jt) { auto Pp = jt->supporting_line().projection(P); auto d = CGAL::squared_distance(Pp, P); @@ -719,224 +694,52 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v return std::make_pair(input_it, closest); }; - // Find the outer perimeter using offset - union - negative offset - std::vector offset_polygons; - for (auto& r : input_polygons) { - auto R = r; - if (!R.is_counterclockwise_oriented()) { - R.reverse_orientation(); - } +private: + using TreeTraits = CGAL::AABB_traits>::iterator>>; + using Tree = CGAL::AABB_tree; - // Overlap removal can also result in close points causing problems when converted into non-exact nt - remove_close_points(R); + const std::vector& polygons_ref_; + std::list> all_segs; + std::unordered_map*, PolygonIt> seg_to_poly; + Tree tree_; - auto ps = create_and_convert_offset_polygon(polygon_offset_distance, R); - for (auto& p : ps) { - if (!p.is_simple()) { - /*{ - std::cerr << "input ["; - bool first = true; - for (auto& pp : r) { - if (!first) { - std::cerr << ","; - } - first = false; - std::cerr << "(" << pp.x() << "," << pp.y() << ")"; - } - std::cerr << "]" << std::endl; - } + std::map::const_iterator> input_polygon_boundary_cache_; +}; - { - std::cerr << "["; - bool first = true; - for (auto& pp : p) { - if (!first) { - std::cerr << ","; - } - first = false; - std::cerr << "(" << pp.x() << "," << pp.y() << ")"; - } - std::cerr << "]" << std::endl; - }*/ - - throw std::runtime_error("Complex polygon originated from offset"); - } - } - offset_polygons.insert(offset_polygons.end(), ps.begin(), ps.end()); - } - -#ifdef SVGFILL_DEBUG - for (auto it = offset_polygons.begin(); it != offset_polygons.end(); ++it) { - write_polygon_to_obj(obj, vi, true, *it, "offset_poly_" + std::to_string(std::distance(offset_polygons.begin(), it))); - write_polygon_to_svg(svg, *it); - } -#endif - - // Perform Boolean union on the offset polygons - std::vector unioned_polygons; - CGAL::join(offset_polygons.begin(), offset_polygons.end(), std::back_inserter(unioned_polygons)); - - if (unioned_polygons.size() > 1) { - // @todo this is currently one of the major limitations in the code that still can be eliminated - // by grouping the input polygons by their perimiter polygon in unioned_polygons - std::sort(unioned_polygons.begin(), unioned_polygons.end(), [](auto& p, auto& q) { return p.outer_boundary().area() > q.outer_boundary().area(); }); - } - -#ifdef SVGFILL_DEBUG - write_polygon_to_obj(obj, vi, true, unioned_polygons.front().outer_boundary(), "offset_poly_joined"); - write_polygon_to_svg(svg, unioned_polygons.front().outer_boundary()); - -#endif - - Polygon_2 fused_removed_close_points; - { - std::vector> ps; - auto& p = unioned_polygons.front().outer_boundary(); - ps.reserve(p.size()); - auto I = p.begin(); - auto J = I + 1; - for (;; ++J) { - bool last = false; - if (J == p.end()) { - J = p.begin(); - last = true; - } - // if (CGAL::squared_distance(*I, *J) > (polygon_offset_distance * polygon_offset_distance)) { - if (CGAL::squared_distance(*I, *J) > (1.e-4 * 1.e-4)) { - ps.push_back(*J); - I = J; - } - if (last) { - break; - } - } - fused_removed_close_points = Polygon_2(ps.begin(), ps.end()); - } - - // Apply negative offset to get the outer perimeter polygon - auto inner_offset = create_and_convert_offset_polygon( - // Because polygon_offset is inexact, make sure our inset distance is slightly larger - // std::nexttoward(-polygon_offset_distance, -std::numeric_limits::infinity()), - - // 1.e-8 even was too little and still resulted in slivers of triangle around the perimeter - -polygon_offset_distance - OUTER_PERIMITER_ADDITIONAL_INSET_AMOUNT, - fused_removed_close_points); - -#ifdef SVGFILL_DEBUG - write_polygon_to_obj(obj, vi, true, inner_offset.front(), "joined_inset"); - write_polygon_to_svg(svg, inner_offset.front()); -#endif - - /* - // there is non-insignificant chance that around the outer boundary, vertices are located in - // between of the input polyhedra, but intermediate vertices result in triangles that will no longer - // span between the two spaces with two edges and therefore cause the topological centre line - // to no run up to the center. Eliminate all vertices that are not on the polyhedral boundary of polygon. - - // this theory proved to be false. once we have topological end points in our graph that are - // connected to input polyhedra to form closed cells, we move those topological end points to - // the average of the input polyhedra corner points, thus effectively also moving them outwards. - { - for (auto& i : inner_offset) { - std::vector> ps; - for (auto& p : i) { - if (input_polygon_boundary(p, 1.e-3) != input_polygons.end()) { - ps.push_back(p); - } - } - i = Polygon_2(ps.begin(), ps.end()); +Polygon_2 subdivide_polygon(double max_distance, const Polygon_2 & p) { + std::vector points; + for (auto it = p.edges_begin(); it != p.edges_end(); ++it) { + const auto& seg = *it; + auto num_splits = (int)std::ceil(std::sqrt(CGAL::to_double(seg.squared_length())) / max_distance) - 1; + points.push_back(seg.source()); + for (auto i = 0; i < num_splits; ++i) { + auto d = (seg.target() - seg.source()) / (num_splits + 1) * (i + 1); + points.push_back(seg.source() + d); } } + return Polygon_2(points.begin(), points.end()); +}; -#ifdef SVGFILL_DEBUG - write_polygon_to_obj(obj, vi, true, inner_offset.front(), "joined_inset_cleaned"); - write_polygon_to_svg(svg, inner_offset.front()); -#endif - */ - - // Subtract original polygons from outer perimeter - std::vector difference_result, difference_result_subdivided; - for (auto& i : inner_offset) { - std::vector working_copy; - working_copy.emplace_back(i); - - for (auto& r : input_polygons) { - std::vector temp_working_copy; - for (auto& wc : working_copy) { - CGAL::difference(wc, r, std::back_inserter(temp_working_copy)); - } - working_copy = temp_working_copy; - } - difference_result.insert(difference_result.end(), working_copy.begin(), working_copy.end()); +Polygon_with_holes_2 subdivide_polygon(double max_distance, const Polygon_with_holes_2& pwh) { + Polygon_2 outer = subdivide_polygon(max_distance, pwh.outer_boundary()); + std::vector holes; + for (auto hit = pwh.holes_begin(); hit != pwh.holes_end(); ++hit) { + holes.push_back(subdivide_polygon(max_distance, *hit)); } + return Polygon_with_holes_2(outer, holes.begin(), holes.end()); +}; - // subdivide difference_result to have better behave triangulation - - { - const double max_distance = polygon_offset_distance / 8.; - auto subdivide_polygon = [max_distance](const Polygon_2& p) { - std::vector points; - for (auto it = p.edges_begin(); it != p.edges_end(); ++it) { - const auto& seg = *it; - auto num_splits = (int)std::ceil(std::sqrt(CGAL::to_double(seg.squared_length())) / max_distance) - 1; - points.push_back(seg.source()); - for (auto i = 0; i < num_splits; ++i) { - auto d = (seg.target() - seg.source()) / (num_splits + 1) * (i + 1); - points.push_back(seg.source() + d); - } - } - return Polygon_2(points.begin(), points.end()); - }; - - for (auto& pwh : difference_result) { - // Subdivide outer boundary - Polygon_2 outer = subdivide_polygon(pwh.outer_boundary()); - // Subdivide holes - std::vector holes; - for (auto hit = pwh.holes_begin(); hit != pwh.holes_end(); ++hit) { - holes.push_back(subdivide_polygon(*hit)); - } - // Construct new Polygon_with_holes_2 - difference_result_subdivided.push_back(Polygon_with_holes_2(outer, holes.begin(), holes.end())); - } - } - -#ifdef SVGFILL_DEBUG - for (auto it = difference_result_subdivided.begin(); it != difference_result_subdivided.end(); ++it) { - auto i = std::distance(difference_result_subdivided.begin(), it); - write_polygon_to_obj(obj, vi, true, it->outer_boundary(), "difference_result_subdivided_" + std::to_string(i)); - write_polygon_to_svg(svg, it->outer_boundary()); - for (auto& p : it->holes()) { - write_polygon_to_obj(obj, vi, true, p, "difference_result_subdivided_" + std::to_string(i)); - write_polygon_to_svg(svg, p); - } - } -#endif - - std::list> triangular_polygons; - - for (auto& pwh : difference_result_subdivided) { - CGAL::Polygon_triangulation_decomposition_2 decompositor; - decompositor(pwh, std::back_inserter(triangular_polygons)); - } - - triangular_polygons.erase(std::remove_if(triangular_polygons.begin(), triangular_polygons.end(), [](const CGAL::Polygon_2& p) { - return CGAL::to_double(p.area()) < 1.e-8; - }), triangular_polygons.end()); - -#ifdef SVGFILL_DEBUG - for (auto it = triangular_polygons.begin(); it != triangular_polygons.end(); ++it) { - write_polygon_to_obj(obj, vi, false, *it, "tri_" + std::to_string(std::distance(triangular_polygons.begin(), it))); - write_polygon_to_svg(svg, *it); - } -#endif - +std::tuple< + std::map>, + std::map>, + std::map, std::vector*>>> +build_line_graph(const std::vector& input_polygons, SegmentLookup& segment_lookup, const std::vector& triangular_polygons) { // Build maps of triangle -> edge and edge -> triangle in order to do traversal on the 'corridor mesh' - std::map, std::vector*>> segment_to_facet; - std::map, std::vector*>> segment_to_input_facet; + std::map, std::vector*>> segment_to_facet; + std::map, std::vector*>> segment_to_input_facet; std::map, Point_2> segment_to_midpoint; std::map> midpoint_to_segment; - std::map*, std::vector>> facet_to_segment; + std::map*, std::vector>> facet_to_segment; for (auto& tri : triangular_polygons) { for (size_t i = 0; i < 3; ++i) { @@ -950,26 +753,14 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } } - // This part is the most computationally expensive. Caching effectively halves the lookup time here, since every vertex has two outgoing edges. - std::map input_polygon_boundary_cache; - auto cached_input_polygon_boundary = [&](const Point_2& p, double tol = 1e-5) -> decltype(input_polygons.begin()) - { - auto it = input_polygon_boundary_cache.find(p); - if (it == input_polygon_boundary_cache.end()) { - auto index = input_polygon_boundary(p, tol); - input_polygon_boundary_cache[p] = index; - return index; - } else { - return it->second; - } - }; + // @todo The smarter thing to do probably after creating the corridor mesh, register segments wrt to originating input polygon(s) and maintain that mapping when subdividing // Register midpoints on the edges within the 'corridor mesh' that span multiple input polygons for (auto& p : segment_to_facet) { auto center = CGAL::ORIGIN + (((p.first.first - CGAL::ORIGIN) + (p.first.second - CGAL::ORIGIN)) / 2); - auto p1index = cached_input_polygon_boundary(p.first.first); - auto p2index = cached_input_polygon_boundary(p.first.second); + auto p1index = segment_lookup.input_polygon_boundary(p.first.first); + auto p2index = segment_lookup.input_polygon_boundary(p.first.second); segment_to_input_facet[p.first].push_back(&*p1index); segment_to_input_facet[p.first].push_back(&*p2index); @@ -978,17 +769,8 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v segment_to_midpoint[p.first] = center; midpoint_to_segment[center] = p.first; } - - if (p1index != input_polygons.end() && p2index != input_polygons.end() && p1index != p2index) { - segment_to_midpoint[p.first] = center; - midpoint_to_segment[center] = p.first; - } } -#ifdef SVGFILL_DEBUG - obj << "o network_1\n"; -#endif - // Observe corridor mesh topology to join edge midpoints into a network std::map> line_graph; for (auto& p : segment_to_midpoint) { @@ -1000,20 +782,15 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v decltype(segment_to_midpoint)::const_iterator it; if ((it = segment_to_midpoint.find(r)) != segment_to_midpoint.end()) { line_graph[p.second].push_back(it->second); - -#ifdef SVGFILL_DEBUG - obj << "v " << CGAL::to_double(p.second.x()) << " " << CGAL::to_double(p.second.y()) << " 0\n"; - obj << "v " << CGAL::to_double(it->second.x()) << " " << CGAL::to_double(it->second.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; - - svg << "second.x()) << "\" y2=\"" << CGAL::to_double(it->second.y()) << "\" />"; -#endif } } } } + return {line_graph, midpoint_to_segment, segment_to_input_facet}; +} + +std::set> find_triangles(const std::map>& line_graph) { // Find triangles in this network often occuring at junctions in the corridor mesh std::set> triangles; std::function&)> find_triangles_recursive; @@ -1024,7 +801,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v const std::vector& neighbors_current = line_graph.at(path.back()); if (std::find(neighbors_current.begin(), neighbors_current.end(), path.front()) != neighbors_current.end()) { // We found a triangle, add it to the set - Triangle triangle = { path[0], path[1], path[2] }; + Triangle triangle = {path[0], path[1], path[2]}; std::sort(triangle.begin(), triangle.end()); triangles.insert(triangle); } @@ -1037,29 +814,28 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v if (std::find(path.begin(), path.end(), neighbor) == path.end()) { path.push_back(neighbor); find_triangles_recursive(path); - path.pop_back(); // Backtrack + path.pop_back(); // Backtrack } } }; for (auto& p : line_graph) { - std::vector ps = { p.first }; + std::vector ps = {p.first}; find_triangles_recursive(ps); } - // For every triangle found in the network we eliminate one edge to break the cycle - // The edge we eliminate is the edge with the greatest angle with any of it's neighbours + return triangles; +} - // non exact time, we need sqrt +std::set> eliminate_triangles(const std::map>& line_graph) { + auto triangles = find_triangles(line_graph); + + // @todo this currently uses a simple cartesian kernel for performance for support of sqrt, but + // this should be possible to rewrite as ratios/slopes in the exact kernel as well using SK = CGAL::Simple_cartesian; CGAL::Cartesian_converter C{}; -#ifdef SVGFILL_DEBUG - obj << "o eliminated\n"; -#endif - std::set> eliminated_segments; - for (auto& t : triangles) { Triangle st; std::transform(t.begin(), t.end(), st.begin(), C); @@ -1075,7 +851,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v double max_abs_dot = 0.; { - auto& ni = line_graph[t[i]]; + auto& ni = line_graph.find(t[i])->second; for (auto& n : ni) { if (std::find(t.begin(), t.end(), n) == t.end()) { // not contained in triangle @@ -1092,7 +868,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } { - auto& nj = line_graph[t[j]]; + auto& nj = line_graph.find(t[j])->second; for (auto& n : nj) { if (std::find(t.begin(), t.end(), n) == t.end()) { // not contained in triangle @@ -1106,7 +882,6 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } } } - } if (max_abs_dot < global_min_abs_dot) { @@ -1119,176 +894,132 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v auto i = global_min_abs_dot_index; auto j = (i + 2) % 3; - eliminated_segments.insert({ t[i], t[j] }); - eliminated_segments.insert({ t[j], t[i] }); - -#ifdef SVGFILL_DEBUG - obj << "v " << st[j].x() << " " << st[j].y() << " 0\n"; - obj << "v " << st[i].x() << " " << st[i].y() << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; - - svg << ""; -#endif + eliminated_segments.insert({t[i], t[j]}); + eliminated_segments.insert({t[j], t[i]}); } - } - Graph2D G2(line_graph); - for (auto& e : eliminated_segments) { - G2.remove_edge(e.first, e.second); + return eliminated_segments; +} + +bool is_parallel_2degree_node(Graph2D::vertex_const_iterator vit) { + auto it = vit->second.begin(); + auto& P = *it++; + auto& Q = *it++; + auto e1 = P - vit->first; + auto e2 = vit->first - Q; + if (e1.squared_length() == 0 || e2.squared_length() == 0) { + // @todo why does this happen? + return false; } + e1 /= std::sqrt(CGAL::to_double(e1.squared_length())); + e2 /= std::sqrt(CGAL::to_double(e2.squared_length())); + return std::abs(CGAL::to_double(e1 * e2)) > (1. - 1.e-5); +}; - auto G = G2.weld_vertices(); -#ifdef SVGFILL_DEBUG - obj << "o network_2\n"; - for (auto it = G.edges_begin(); it != G.edges_end(); ++it) { - obj << "v " << CGAL::to_double(it->first.x()) << " " << CGAL::to_double(it->first.y()) << " 0\n"; - obj << "v " << CGAL::to_double(it->second.x()) << " " << CGAL::to_double(it->second.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; - } - obj << std::flush; -#endif - - auto is_parallel_2degree_node = [](decltype(G)::vertex_const_iterator vit) { - auto it = vit->second.begin(); - auto& P = *it++; - auto& Q = *it++; - auto e1 = P - vit->first; - auto e2 = vit->first - Q; - if (e1.squared_length() == 0 || e2.squared_length() == 0) { - // @todo why does this happen? - return false; - } - e1 /= std::sqrt(CGAL::to_double(e1.squared_length())); - e2 /= std::sqrt(CGAL::to_double(e2.squared_length())); - return std::abs(CGAL::to_double(e1 * e2)) > (1. - 1.e-5); - }; - - { - // Remove colinear vertices - size_t n_vertices_removed = 0; - for (auto vit = G.vertices_begin(); vit != G.vertices_end();) { - if (vit->second.size() == 2) { - if (is_parallel_2degree_node(vit)) { - vit = G.eliminate_vertex(vit); - ++n_vertices_removed; - } else { - ++vit; - } +void eliminate_colinear_vertices(Graph2D& G) { + size_t n_vertices_removed = 0; + for (auto vit = G.vertices_begin(); vit != G.vertices_end();) { + if (vit->second.size() == 2) { + if (is_parallel_2degree_node(vit)) { + vit = G.eliminate_vertex(vit); + ++n_vertices_removed; } else { ++vit; } + } else { + ++vit; } - // std::cout << "Eliminated " << n_vertices_removed << " vertices" << std::endl; } +} - // Ortho edge slide - { - std::list> edges_to_remove, edges_to_insert; +void edge_slide(Graph2D& G) { + std::list> edges_to_remove, edges_to_insert; - for (auto vit = G.vertices_begin(); vit != G.vertices_end(); ++vit) { - auto& selected = vit->first; + for (auto vit = G.vertices_begin(); vit != G.vertices_end(); ++vit) { + auto& selected = vit->first; - if (vit->second.size() >= 3) { - for (auto vjt = vit->second.begin(); vjt != vit->second.end(); ++vjt) { - auto& neighbour = *vjt; - bool processed_neighbour = false; + if (vit->second.size() >= 3) { + for (auto vjt = vit->second.begin(); vjt != vit->second.end(); ++vjt) { + auto& neighbour = *vjt; + bool processed_neighbour = false; - if (G.find(neighbour)->second.size() == 2 && !is_parallel_2degree_node(G.find(neighbour))) { - auto vkt = G.find(neighbour)->second.begin(); - if (selected == *vkt) { - vkt++; - } - auto& other = *vkt; + if (G.find(neighbour)->second.size() == 2 && !is_parallel_2degree_node(G.find(neighbour))) { + auto vkt = G.find(neighbour)->second.begin(); + if (selected == *vkt) { + vkt++; + } + auto& other = *vkt; - if ((other - neighbour).squared_length() < (neighbour - selected).squared_length()) { - continue; - } + if ((other - neighbour).squared_length() < (neighbour - selected).squared_length()) { + continue; + } - auto incoming = CGAL::Ray_2(other, neighbour - other); - boost::optional> closest_neighbouring_segment; - boost::optional> closest_intersection_point; - K::FT sq_distance_along_ray = std::numeric_limits::infinity(); + auto incoming = CGAL::Ray_2(other, neighbour - other); + boost::optional> closest_neighbouring_segment; + boost::optional> closest_intersection_point; + K::FT sq_distance_along_ray = std::numeric_limits::infinity(); - for (auto vlt = vit->second.begin(); vlt != vit->second.end(); ++vlt) { - auto& other_neighbour = *vlt; - if (vlt != vjt) { - CGAL::Segment_2 neighbouring_segment(selected, other_neighbour); - auto x = CGAL::intersection(incoming, neighbouring_segment); - if (x) { - if (auto* xp = variant_get>(&*x)) { - auto dist = ((*xp) - other).squared_length(); - if (dist < sq_distance_along_ray) { - closest_neighbouring_segment = neighbouring_segment; - closest_intersection_point = *xp; - sq_distance_along_ray = dist; - } + for (auto vlt = vit->second.begin(); vlt != vit->second.end(); ++vlt) { + auto& other_neighbour = *vlt; + if (vlt != vjt) { + CGAL::Segment_2 neighbouring_segment(selected, other_neighbour); + auto x = CGAL::intersection(incoming, neighbouring_segment); + if (x) { + if (auto* xp = variant_get>(&*x)) { + auto dist = ((*xp) - other).squared_length(); + if (dist < sq_distance_along_ray) { + closest_neighbouring_segment = neighbouring_segment; + closest_intersection_point = *xp; + sq_distance_along_ray = dist; } } } } - - if (closest_intersection_point && closest_neighbouring_segment) { - edges_to_remove.push_back(*closest_neighbouring_segment); - edges_to_remove.push_back({ neighbour, selected }); - edges_to_insert.push_back({ closest_neighbouring_segment->source(), *closest_intersection_point }); - edges_to_insert.push_back({ closest_neighbouring_segment->target(), *closest_intersection_point }); - edges_to_insert.push_back({ neighbour, *closest_intersection_point }); - - processed_neighbour = true; - } } - if (processed_neighbour) { - // Only one neigbour is processed because otherwise we obtain intersections - break; + + if (closest_intersection_point && closest_neighbouring_segment) { + edges_to_remove.push_back(*closest_neighbouring_segment); + edges_to_remove.push_back({neighbour, selected}); + edges_to_insert.push_back({closest_neighbouring_segment->source(), *closest_intersection_point}); + edges_to_insert.push_back({closest_neighbouring_segment->target(), *closest_intersection_point}); + edges_to_insert.push_back({neighbour, *closest_intersection_point}); + + processed_neighbour = true; } } + if (processed_neighbour) { + // Only one neigbour is processed because otherwise we obtain intersections + break; + } } } - - for (auto& s : edges_to_remove) { - G.remove_edge(s.source(), s.target()); - } - - - for (auto& s : edges_to_insert) { - G.insert(s.source(), s.target()); - } - -#ifdef SVGFILL_DEBUG - obj << "o network_3\n"; - for (auto it = G.edges_begin(); it != G.edges_end(); ++it) { - obj << "v " << CGAL::to_double(it->first.x()) << " " << CGAL::to_double(it->first.y()) << " 0\n"; - obj << "v " << CGAL::to_double(it->second.x()) << " " << CGAL::to_double(it->second.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; - } -#endif } - // Now plot the edges on an arrangement in order to find planar cycles - // and merge the corridor-halves with their neighbouring input polygon - - Arrangement_2 arr; - - for (auto it = G.edges_begin(); it != G.edges_end(); ++it) { - if (it->first == it->second) { - continue; - } - CGAL::insert(arr, Segment_2(it->first, it->second)); + for (auto& s : edges_to_remove) { + G.remove_edge(s.source(), s.target()); } - std::list> move_ops; - std::list> edge_ops; + for (auto& s : edges_to_insert) { + G.insert(s.source(), s.target()); + } +} + +std::list> extend_end_vertices_based_on_input( + const Graph2D& G, + const std::map>& midpoint_to_segment, + const std::map, std::vector*>>& segment_to_input_facet, + const Polygon_list& inner_offset, + const SegmentLookup& segment_lookup +){ + std::list> constructed_segments; for (auto it = G.vertices_begin(); it != G.vertices_end(); ++it) { if (it->second.size() == 1) { auto& M = it->first; - decltype(midpoint_to_segment)::mapped_type* q = nullptr; + const std::pair* q = nullptr; if (midpoint_to_segment.find(M) == midpoint_to_segment.end()) { typename K::FT min_sq_distance = std::numeric_limits::infinity(); @@ -1299,7 +1030,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } } } else { - q = &midpoint_to_segment[M]; + q = &midpoint_to_segment.find(M)->second; } if (q == nullptr) { @@ -1309,7 +1040,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v bool handled_as_graph_path = false; // distance from unioned - shoot ray? - if (segment_to_input_facet[*q].size() == 2) { + if (segment_to_input_facet.find(*q)->second.size() == 2) { for (auto& bnd : inner_offset) { // if point M is contained in bnd interior: if (bnd.has_on_bounded_side(M)) { @@ -1339,10 +1070,10 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v Graph2D GGG(bnd); GGG.refine(*GGG.query(*closest_intersection_point, 0.01), *closest_intersection_point); - std::array>, 2> input_points = { { {}, {} } }; + std::array>, 2> input_points = {{{}, {}}}; size_t i = 0; - for (auto& fac : segment_to_input_facet[*q]) { + for (auto& fac : segment_to_input_facet.find(*q)->second) { for (auto it = fac->vertices_begin(); it != fac->vertices_end(); ++it) { auto seg = GGG.query(*it, 0.01); if (seg) { @@ -1361,13 +1092,13 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v if (!a1.empty() && !a2.empty()) { if (M != *closest_intersection_point) { - edge_ops.push_front({ M, *closest_intersection_point }); + constructed_segments.push_front({M, *closest_intersection_point}); } for (auto it = a1.begin(); it != a1.end() && std::next(it) != a1.end(); ++it) { - edge_ops.push_front({ *it, *(std::next(it)) }); + constructed_segments.push_front({*it, *(std::next(it))}); } for (auto it = a2.begin(); it != a2.end() && std::next(it) != a2.end(); ++it) { - edge_ops.push_front({ *it, *(std::next(it)) }); + constructed_segments.push_front({*it, *(std::next(it))}); } handled_as_graph_path = true; @@ -1381,8 +1112,8 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v if (!handled_as_graph_path) { // else we choose to map point to the midpoint of the found two close points. - auto pq = close_input_point(q->first); - auto pr = close_input_point(q->second); + auto pq = segment_lookup.close_input_point(q->first); + auto pr = segment_lookup.close_input_point(q->second); auto Q = pq.second; auto R = pr.second; @@ -1392,16 +1123,16 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v // where Q and R are co-located, because the point R' is further away // in that case M + M-Q should gives is x that we then project onto the // input boundary - // - // - // ┌───────┐ - // │ │ - // │ │ - // │ │ - // └───────o <--Q,R - // - // ────────o <--M - // + // + // + // ┌───────┐ + // │ │ + // │ │ + // │ │ + // └───────o <--Q,R + // + // ────────o <--M + // // ┌───────x───────────────o <---R' // │ │ // │ │ @@ -1410,93 +1141,22 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v // └───────────────────────┘ // @todo is this projection actually necessary or is it already 'exact enough'? - R = project_input_point(M + (M - Q)).second; + R = segment_lookup.project_input_point(M + (M - Q)).second; } auto avg = CGAL::ORIGIN + ((Q - CGAL::ORIGIN) + (R - CGAL::ORIGIN)) / 2; - move_ops.push_front({ M, avg }); - edge_ops.push_front({ avg, Q }); - edge_ops.push_front({ avg, R }); - + constructed_segments.push_front({M, avg}); + constructed_segments.push_front({avg, Q}); + constructed_segments.push_front({avg, R}); } } } -#ifdef SVGFILL_DEBUG - obj << "o network_4\n"; -#endif + return constructed_segments; +} - // note that we actually don't move but draw an edge - for (auto& pq : move_ops) { - if (pq.first == pq.second) { - continue; - } - CGAL::insert(arr, Segment_2(pq.first, pq.second)); - -#ifdef SVGFILL_DEBUG - obj << "v " << CGAL::to_double(pq.first.x()) << " " << CGAL::to_double(pq.first.y()) << " 0\n"; - obj << "v " << CGAL::to_double(pq.second.x()) << " " << CGAL::to_double(pq.second.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; -#endif - } - - - for (auto& pq : edge_ops) { - if (pq.first == pq.second) { - continue; - } - CGAL::insert(arr, Segment_2(pq.first, pq.second)); - -#ifdef SVGFILL_DEBUG - obj << "v " << CGAL::to_double(pq.first.x()) << " " << CGAL::to_double(pq.first.y()) << " 0\n"; - obj << "v " << CGAL::to_double(pq.second.x()) << " " << CGAL::to_double(pq.second.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; -#endif - } - - // Plot input polygons - for (auto& poly : input_polygons) { - for (size_t i = 0; i != poly.size(); ++i) { - auto j = (i + 1) % poly.size(); - if (poly.vertex(i) == poly.vertex(j)) { - continue; - } - CGAL::insert(arr, Segment_2(poly.vertex(i), poly.vertex(j))); - } - } - - -#ifdef SVGFILL_DEBUG - { - obj << "o arrangement_1\n"; - for (auto it = arr.edges_begin(); it != arr.edges_end(); ++it) { - auto& p = it->source()->point(); - auto& q = it->target()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - obj << "v " << CGAL::to_double(q.x()) << " " << CGAL::to_double(q.y()) << " 0\n"; - obj << "l " << vi++; - obj << " " << vi++ << "\n"; - } - } -#endif - - /* { - // debug, add outer bounds so that we can plot the face for any remaining edges - auto poly = unioned_polygons.front().outer_boundary(); - for (size_t i = 0; i != poly.size(); ++i) { - auto j = (i + 1) % poly.size(); - CGAL::insert(arr, Segment_2(poly.vertex(i), poly.vertex(j))); - } - } */ - - // Now loop over the arrangement faces, when a face coincides with a point on the - // corridor network we know it needs to be joined with an input polygon. In that - // case the edges need to be eliminated that correspond to original geometry. - - size_t face_id = 0; +void fuse_corridor_halves_with_input(Arrangement_2& arr, Graph2D& G, SegmentLookup& segment_lookup, const Polygon_list& input_polygons, DebugWriter& debug_output) { std::set edges_to_remove; for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it) { @@ -1535,7 +1195,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v auto& p = curr->source()->point(); auto& q = curr->target()->point(); auto center = CGAL::ORIGIN + (((p - CGAL::ORIGIN) + (q - CGAL::ORIGIN)) / 2); - auto p1index = input_polygon_boundary(center); + auto p1index = segment_lookup.input_polygon_boundary(center); const bool on_orig_bound = p1index != input_polygons.end(); if (on_orig_bound) { if (edges_to_remove.find(curr->twin()) != edges_to_remove.end()) { @@ -1553,7 +1213,7 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v auto& p = curr->source()->point(); auto& q = curr->target()->point(); auto center = CGAL::ORIGIN + (((p - CGAL::ORIGIN) + (q - CGAL::ORIGIN)) / 2); - auto p1index = input_polygon_boundary(center); + auto p1index = segment_lookup.input_polygon_boundary(center); const bool on_orig_bound = p1index != input_polygons.end(); if (on_orig_bound) { if (edges_to_remove.find(curr->twin()) != edges_to_remove.end()) { @@ -1569,126 +1229,314 @@ void arrange_cgal_polygons(const std::vector& input_polygons_, std::v } } } - -#ifdef SVGFILL_DEBUG - write_polygon_to_svg(svg, circ_to_poly(it->outer_ccb())); - - obj << "o " << "face_"; - if (is_corridor) { - obj << "corri_"; - } - obj << face_id++ << "\n"; - - std::ostringstream oss; - - { - auto vv = vi; - auto curr = it->outer_ccb(); - do { - auto& p = curr->source()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - oss << "l " << vi++; - ++curr; - if (curr == it->outer_ccb()) { - oss << " " << vv << "\n"; - } else { - oss << " " << vi << "\n"; - } - } while (curr != it->outer_ccb()); - } - - for (auto jt = it->inner_ccbs_begin(); jt != it->inner_ccbs_end(); ++jt) { - auto vv = vi; - auto curr = *jt; - do { - auto& p = curr->source()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - oss << "l " << vi++; - ++curr; - if (curr == *jt) { - oss << " " << vv << "\n"; - } else { - oss << " " << vi << "\n"; - } - } while (curr != *jt); - } - - obj << oss.str(); -#endif } size_t remove_id = 0; for (auto& e : edges_to_remove) { -#ifdef SVGFILL_DEBUG - obj << "o " << "remove_" << remove_id++ << "\n"; - { - auto& p = e->source()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - } - { - auto& p = e->target()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - } - obj << "l " << vi++; - obj << " " << vi++ << std::endl; -#endif + debug_output.write_segment(e->source()->point(), e->target()->point(), "arr_remove_edge_" + std::to_string(remove_id++)); CGAL::remove_edge(arr, e); } +} + +class timer { + class entry { + public: + entry(std::map::const_iterator start_it) + : start_it(start_it) {} + void stop() { + auto end = std::chrono::high_resolution_clock::now(); + auto duration = std::chrono::duration(end - start_it->second).count(); + std::cerr << "Timing for " << start_it->first << ": " << duration << " ms" << std::endl; + } + + private: + std::map::const_iterator start_it; + }; + + public: + entry start(const std::string& name) { + return timings_.insert({name, std::chrono::high_resolution_clock::now()}).first; + } + + private: + std::map< + std::string, + std::chrono::high_resolution_clock::time_point> + timings_; +}; + +void arrange_cgal_polygons(const std::vector& input_polygons_, std::vector& output_polygons, double polygon_offset_distance = -1.) { + static const double OVERLAP_RESOLUTION_DISTANCE = 1.e-2; + // even larger amount of inset so that outer perimeter is safely within all input polygons even when overlap resolution is applied + // no, `1.e-2 + 1.e-5` creates issues with the outer perimeter, are there other tolerances in play? + static const double OUTER_PERIMITER_ADDITIONAL_INSET_AMOUNT = 1.e-5; + +#ifdef SVGFILL_DEBUG + auto t = std::time(nullptr); + auto tm = *std::localtime(&t); + + std::ostringstream oss; + oss << std::put_time(&tm, "arrangement_%Y%m%d%H%M%S"); + auto now = oss.str(); + DebugWriter debug_output(true, now); +#else + DebugWriter debug_output(false, ""); +#endif + + timer timer; + + auto t0 = timer.start("input"); + + debug_output.write_polygons(input_polygons_, "input"); + + if (polygon_offset_distance < 0.) { + polygon_offset_distance = estimate_polygon_offset_distance(input_polygons_); + } + + // Create copy to make mutable for cleaning + auto input_polygons = input_polygons_; + + for (auto& polygon : input_polygons) { + clean_polygon(polygon); + } + + { + decltype(input_polygons) split_polygons; + for (auto& poly : input_polygons) { + split_self_intersecting_polygon(poly, std::back_inserter(split_polygons)); + } + std::swap(input_polygons, split_polygons); + } + + t0.stop(); + t0 = timer.start("overlap elimination"); + + eliminate_overlaps(OVERLAP_RESOLUTION_DISTANCE, input_polygons); + + t0.stop(); + + // [NB Nov 6] we cannot do this anymore because it could revert the spacing between input polygons + // that touch in the corner. + // Now that overlaps/touches at corners are handled more locally only a small indent is produced + // which would be undone by means of an inset+offset. + // + // [NB Nov 10] this is actually still necessary though, but we apply a much smaller distance now + // to keep the overlap eliminations in tact + // + // Inset-offset to remove tiny details that may cause enourmous spikes in offsets + for (auto& r : input_polygons) { + smooth_polygon(-polygon_offset_distance / 10000., r); + } + + debug_output.write_polygons(input_polygons, "processed_input"); + + SegmentLookup segment_lookup(input_polygons); + + t0 = timer.start("outer perimeter"); + + // Find the outer perimeter using offset - union - negative offset + std::vector offset_polygons; + for (auto& r : input_polygons) { + auto R = r; + if (!R.is_counterclockwise_oriented()) { + R.reverse_orientation(); + } + + // Overlap removal can also result in close points causing problems when converted into non-exact nt + remove_close_points(R); + + auto ps = create_and_convert_offset_polygon(polygon_offset_distance, R); + for (auto& p : ps) { + if (!p.is_simple()) { + throw std::runtime_error("Complex polygon originated from offset"); + } + } + offset_polygons.insert(offset_polygons.end(), ps.begin(), ps.end()); + } + + debug_output.write_polygons(offset_polygons, "offset_input"); + + // Perform Boolean union on the offset polygons + std::vector unioned_polygons; + CGAL::join(offset_polygons.begin(), offset_polygons.end(), std::back_inserter(unioned_polygons)); + + if (unioned_polygons.size() > 1) { + // @todo this is currently one of the major limitations in the code that still can be eliminated + // by grouping the input polygons by their perimiter polygon in unioned_polygons + std::sort(unioned_polygons.begin(), unioned_polygons.end(), [](auto& p, auto& q) { return p.outer_boundary().area() > q.outer_boundary().area(); }); + } + + debug_output.write_polygon(unioned_polygons.front().outer_boundary(), "offset_joined"); + + Polygon_2 fused_removed_close_points = unioned_polygons.front().outer_boundary(); + remove_close_points(fused_removed_close_points, 1.e-4); + + // Apply negative offset to get the outer perimeter polygon + auto inner_offset = create_and_convert_offset_polygon( + // Because polygon_offset is inexact, make sure our inset distance is slightly larger + // std::nexttoward(-polygon_offset_distance, -std::numeric_limits::infinity()), + + // 1.e-8 even was too little and still resulted in slivers of triangle around the perimeter + -polygon_offset_distance - OUTER_PERIMITER_ADDITIONAL_INSET_AMOUNT, + fused_removed_close_points); + + debug_output.write_polygons(inner_offset, "outer_perimiter"); + + t0.stop(); + t0 = timer.start("corridor creation"); + + // Subtract original polygons from outer perimeter + std::vector difference_result, difference_result_subdivided; + for (auto& i : inner_offset) { + std::vector working_copy; + working_copy.emplace_back(i); + + for (auto& r : input_polygons) { + std::vector temp_working_copy; + for (auto& wc : working_copy) { + CGAL::difference(wc, r, std::back_inserter(temp_working_copy)); + } + working_copy = temp_working_copy; + } + difference_result.insert(difference_result.end(), working_copy.begin(), working_copy.end()); + } + + t0.stop(); + t0 = timer.start("corridor triangulation"); + + // subdivide difference_result to have better more detailed triangulation and therefore less-pronounced artefacts in midpoint network + + for (auto& pwh : difference_result) { + difference_result_subdivided.push_back(subdivide_polygon(polygon_offset_distance / 8., pwh)); + // difference_result_subdivided.push_back(subdivide_polygon(polygon_offset_distance / 64., pwh)); + } + + debug_output.write_polygons(difference_result_subdivided, "corridor_subdivided"); + + std::vector> triangular_polygons; + for (auto& pwh : difference_result_subdivided) { + CGAL::Polygon_triangulation_decomposition_2 decompositor; + decompositor(pwh, std::back_inserter(triangular_polygons)); + } + + t0.stop(); + + /* + * // @todo decide whether this is smart or not + * // Would this not hurt topology too much? + triangular_polygons.erase(std::remove_if(triangular_polygons.begin(), triangular_polygons.end(), [](const CGAL::Polygon_2& p) { + return CGAL::to_double(p.area()) < 1.e-8; + }), triangular_polygons.end()); + */ + + t0 = timer.start("center line"); + + debug_output.write_polygons(triangular_polygons, "triangulated_corridor"); + + auto [line_graph, midpoint_to_segment, segment_to_input_facet] = build_line_graph(input_polygons, segment_lookup, triangular_polygons); + for (auto& p : line_graph) { + for (auto& q : p.second) { + debug_output.write_segment(p.first, q, "network_1"); + } + } + + t0.stop(); + + t0 = timer.start("center line cleaning"); + + auto triangles = find_triangles(line_graph); + + // For every triangle found in the network we eliminate one edge to break the cycle + // The edge we eliminate is the edge with the greatest angle with any of it's neighbours + auto eliminated_segments = eliminate_triangles(line_graph); + + Graph2D G2(line_graph); + for (auto& e : eliminated_segments) { + debug_output.write_segment(e.first, e.second, "eliminated"); + G2.remove_edge(e.first, e.second); + } + + auto 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"); + } + + 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"); + } + + t0.stop(); + + t0 = timer.start("topology"); + + auto segments = extend_end_vertices_based_on_input(G, midpoint_to_segment, segment_to_input_facet, inner_offset, segment_lookup); + + // Now plot the edges on an arrangement in order to find planar cycles + // and merge the corridor-halves with their neighbouring input polygon + Arrangement_2 arr; + G.to_arrangement(arr); + + for (auto& pq : segments) { + if (pq.first == pq.second) { + continue; + } + CGAL::insert(arr, Segment_2(pq.first, pq.second)); + + debug_output.write_segment(pq.first, pq.second, "extended_segments"); + } + + // Write input polygons to arrangement_2 + for (auto& poly : input_polygons) { + for (size_t i = 0; i != poly.size(); ++i) { + auto j = (i + 1) % poly.size(); + if (poly.vertex(i) == poly.vertex(j)) { + continue; + } + CGAL::insert(arr, Segment_2(poly.vertex(i), poly.vertex(j))); + } + } + + // Just for the automatic numbering, create a full vector + std::vector temp; + for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it) { + if (it->is_unbounded()) { + continue; + } + temp.push_back(circ_to_poly(it->outer_ccb())); + } + debug_output.write_polygons(temp, "arr_faces"); + + + /* { + // debug, add outer bounds so that we can plot the face for any remaining edges + auto poly = unioned_polygons.front().outer_boundary(); + for (size_t i = 0; i != poly.size(); ++i) { + auto j = (i + 1) % poly.size(); + CGAL::insert(arr, Segment_2(poly.vertex(i), poly.vertex(j))); + } + } */ + + // Now loop over the arrangement faces, when a face coincides with a point on the + // corridor network we know it needs to be joined with an input polygon. In that + // case the edges need to be eliminated that correspond to original geometry. + + fuse_corridor_halves_with_input(arr, G, segment_lookup, input_polygons, debug_output); + + t0.stop(); for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it) { if (it->is_unbounded()) { continue; } - output_polygons.push_back(circ_to_poly(it->outer_ccb())); - -#ifdef SVGFILL_DEBUG - write_polygon_to_svg(svg, circ_to_poly(it->outer_ccb())); - - obj << "o " << "merged_face_"; - obj << face_id++ << "\n"; - - std::ostringstream oss; - - { - auto vv = vi; - auto curr = it->outer_ccb(); - do { - auto& p = curr->source()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - oss << "l " << vi++; - ++curr; - if (curr == it->outer_ccb()) { - oss << " " << vv << "\n"; - } else { - oss << " " << vi << "\n"; - } - } while (curr != it->outer_ccb()); - } - - for (auto jt = it->inner_ccbs_begin(); jt != it->inner_ccbs_end(); ++jt) { - auto vv = vi; - auto curr = *jt; - do { - auto& p = curr->source()->point(); - obj << "v " << CGAL::to_double(p.x()) << " " << CGAL::to_double(p.y()) << " 0\n"; - oss << "l " << vi++; - ++curr; - if (curr == *jt) { - oss << " " << vv << "\n"; - } else { - oss << " " << vi << "\n"; - } - } while (curr != *jt); - } - - obj << oss.str(); -#endif } -#ifdef SVGFILL_DEBUG - svg << "\n"; -#endif + debug_output.write_polygons(output_polygons, "arr_faces_merged"); } #ifndef SVGFILL_MAIN diff --git a/src/svgfill/src/graph_2d.h b/src/svgfill/src/graph_2d.h index aaaa059b4e..da2b5014ec 100644 --- a/src/svgfill/src/graph_2d.h +++ b/src/svgfill/src/graph_2d.h @@ -2,8 +2,10 @@ #define GRAPH_2D_H #ifdef SVGFILL_DEBUG +#if 0 #include #endif +#endif template class Graph2D { @@ -334,6 +336,16 @@ public: return Graph2D(input_adjacency_list); } + template + void to_arrangement(T& arr) { + for (auto it = edges_begin(); it != edges_end(); ++it) { + if (it->first == it->second) { + continue; + } + CGAL::insert(arr, CGAL::Segment_2(it->first, it->second)); + } + } + void assert_symmetric() { #ifdef SVGFILL_DEBUG #if 0