From 9775dfc2c1ff7d8b6eb8f7fa0cda7005a33d23fd Mon Sep 17 00:00:00 2001 From: Esteban DUGUEPEROUX Date: Tue, 25 Nov 2025 13:50:45 +0100 Subject: [PATCH] Add svgfill directly in IfcOpenShell git repo --- src/.gitattributes => .gitattributes | 0 .gitmodules | 3 + src/svgfill/3rdparty/svgpp | 1 + src/svgfill/CMakeLists.txt | 163 +++ src/svgfill/README.md | 65 + src/svgfill/examples/rects.svg | 80 ++ src/svgfill/examples/rects_output.svg | 1 + src/svgfill/src/arrange_polygons.cpp | 1774 +++++++++++++++++++++++++ src/svgfill/src/graph_2d.h | 509 +++++++ src/svgfill/src/main.cpp | 127 ++ src/svgfill/src/progress.h | 122 ++ src/svgfill/src/svgfill.cpp | 633 +++++++++ src/svgfill/src/svgfill.h | 119 ++ 13 files changed, 3597 insertions(+) rename src/.gitattributes => .gitattributes (100%) create mode 160000 src/svgfill/3rdparty/svgpp create mode 100644 src/svgfill/CMakeLists.txt create mode 100644 src/svgfill/README.md create mode 100644 src/svgfill/examples/rects.svg create mode 100644 src/svgfill/examples/rects_output.svg create mode 100644 src/svgfill/src/arrange_polygons.cpp create mode 100644 src/svgfill/src/graph_2d.h create mode 100644 src/svgfill/src/main.cpp create mode 100644 src/svgfill/src/progress.h create mode 100644 src/svgfill/src/svgfill.cpp create mode 100644 src/svgfill/src/svgfill.h diff --git a/src/.gitattributes b/.gitattributes similarity index 100% rename from src/.gitattributes rename to .gitattributes diff --git a/.gitmodules b/.gitmodules index 3c158b5748..ac12797d0b 100644 --- a/.gitmodules +++ b/.gitmodules @@ -17,3 +17,6 @@ [submodule "src/pyodide/demo-app/wheels"] path = src/pyodide/demo-app/wheels url = https://github.com/IfcOpenShell/wasm-wheels +[submodule "src/svgfill/3rdparty/svgpp"] + path = src/svgfill/3rdparty/svgpp + url = https://github.com/svgpp/svgpp diff --git a/src/svgfill/3rdparty/svgpp b/src/svgfill/3rdparty/svgpp new file mode 160000 index 0000000000..5f1870aa7b --- /dev/null +++ b/src/svgfill/3rdparty/svgpp @@ -0,0 +1 @@ +Subproject commit 5f1870aa7b757718ff5f86bdfb55966fa4f217f9 diff --git a/src/svgfill/CMakeLists.txt b/src/svgfill/CMakeLists.txt new file mode 100644 index 0000000000..155c016389 --- /dev/null +++ b/src/svgfill/CMakeLists.txt @@ -0,0 +1,163 @@ +cmake_minimum_required (VERSION 3.10) + +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) + +project (svgfill) + +cmake_policy(SET CMP0074 NEW) # find_package() uses _ROOT variables. +if (POLICY CMP0144) +cmake_policy(SET CMP0144 NEW) # find_package() uses upper-case _ROOT variables. +endif() +if(POLICY CMP0167) # 3.30 find_package(Boost) to use BoostConfig instead of FindBoost. + cmake_policy(SET CMP0167 OLD) +endif() + +include(GNUInstallDirs) + +# Specify paths to install files +if(NOT BINDIR) + set(BINDIR bin) +endif() +if(NOT IS_ABSOLUTE ${BINDIR}) + set(BINDIR ${CMAKE_INSTALL_BINDIR}) +endif() +message(STATUS "BINDIR: ${BINDIR}") + +if(NOT INCLUDEDIR) + set(INCLUDEDIR include) +endif() +if(NOT IS_ABSOLUTE ${INCLUDEDIR}) + set(INCLUDEDIR ${CMAKE_INSTALL_INCLUDEDIR}) +endif() +message(STATUS "INCLUDEDIR: ${INCLUDEDIR}") + +if(NOT LIBDIR) + set(LIBDIR lib) +endif() +if(NOT IS_ABSOLUTE ${LIBDIR}) + set(LIBDIR ${CMAKE_INSTALL_LIBDIR}) +endif() +message(STATUS "LIBDIR: ${LIBDIR}") + +set(CGAL_LIBRARY_NAMES libCGAL_Core libCGAL_ImageIO libCGAL) + +if(NOT CGAL_INCLUDE_DIR) + find_package(CGAL REQUIRED) + if(NOT CGAL_DIR) + message( + FATAL_ERROR + "CGAL_SUPPORT enabled, but CGAL_INCLUDE_DIR wasn't provided and CGAL package couldn't be found." + ) + endif() + message(STATUS "CGAL: found config at '${CGAL_DIR}'.") + set(CGAL_LIBRARIES CGAL::CGAL) +else() + set(CGAL_INCLUDE_DIR ${CGAL_INCLUDE_DIR} CACHE FILEPATH "CGAL header files") + message(STATUS "Looking for CGAL include files in: ${CGAL_INCLUDE_DIR}") + + if(NOT "${CGAL_LIBRARY_DIR}" STREQUAL "") + set(CGAL_LIBRARY_DIR ${CGAL_LIBRARY_DIR} CACHE FILEPATH "CGAL library files") + message(STATUS "Looking for CGAL library files in: ${CGAL_LIBRARY_DIR}") + endif() + + if(WASM_BUILD) + set(CMAKE_FIND_ROOT_PATH_BACKUP "${CMAKE_FIND_ROOT_PATH}") + set(CMAKE_FIND_ROOT_PATH "") + endif() + + find_library(libCGAL NAMES CGAL PATHS ${CGAL_LIBRARY_DIR} NO_DEFAULT_PATH) + + if(libCGAL) + message(STATUS "CGAL library files found") + foreach(lib ${CGAL_LIBRARY_NAMES}) + string(REPLACE libCGAL "${lib}" lib_path "${libCGAL}") + list(APPEND CGAL_LIBRARIES "${lib_path}") + endforeach() + else() + if(NOT "${CGAL_LIBRARY_DIR}" STREQUAL "") + file(GLOB CGAL_LIBRARIES ${CGAL_LIBRARY_DIR}/CGAL*.lib) + list(LENGTH CGAL_LIBRARY_NAMES num_cgal_library_names) + list(LENGTH CGAL_LIBRARIES num_cgal_libraries) + link_directories("${CGAL_LIBRARY_DIR}") + if(NOT "${num_cgal_library_names}" STREQUAL "${num_cgal_libraries}") + message(FATAL_ERROR "Unable to find CGAL library files, aborting") + endif() + message(STATUS "CGAL library files found") + endif() + endif() + find_library(libGMP NAMES gmp mpir PATHS ${GMP_LIBRARY_DIR} NO_DEFAULT_PATH) + find_library(libMPFR NAMES mpfr PATHS ${MPFR_LIBRARY_DIR} NO_DEFAULT_PATH) + if(NOT libGMP) + message(FATAL_ERROR "Unable to find GMP library files, aborting") + endif() + if(NOT libMPFR) + message(FATAL_ERROR "Unable to find MPFR library files, aborting") + endif() + + list(APPEND CGAL_LIBRARIES "${libMPFR}") + list(APPEND CGAL_LIBRARIES "${libGMP}") +endif(NOT CGAL_INCLUDE_DIR) + +if(WIN32 AND ("$ENV{CONDA_BUILD}" STREQUAL "")) + set(Boost_USE_STATIC_LIBS ON) + set(Boost_USE_MULTITHREADED ON) + if (USE_STATIC_MSVC_RUNTIME) + set(Boost_USE_STATIC_RUNTIME ON) + endif() +else() + # Disable Boost's autolinking as the libraries to be linked to are supplied + # already by CMake, and it's going to conflict if there are multiple, as is + # the case in conda-forge's libboost feedstock. + add_definitions(-DBOOST_ALL_NO_LIB) + if(WIN32) + # Necessary for boost version >= 1.67 + set(BCRYPT_LIBRARIES "bcrypt.lib") + endif() +endif() + +if (MSVC) + add_definitions(-bigobj) +endif() + +find_package(Boost) +message(STATUS "Boost include files found in ${Boost_INCLUDE_DIRS}") +find_package(LibXml2 REQUIRED) + +if(WASM_BUILD) + set(CMAKE_FIND_ROOT_PATH "${CMAKE_FIND_ROOT_PATH_BACKUP}") +endif() + + + +include_directories(${Boost_INCLUDE_DIRS} ${LIBXML2_INCLUDE_DIR} + ${CGAL_INCLUDE_DIR} ${GMP_INCLUDE_DIR} ${MPFR_INCLUDE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}/3rdparty/svgpp/include +) + +file(GLOB LIB_H_FILES src/*.h) +file(GLOB LIB_CPP_FILES src/svgfill.cpp src/arrange_polygons.cpp) +set(LIB_SRC_FILES ${LIB_H_FILES} ${LIB_CPP_FILES}) +add_library(svgfill ${LIB_SRC_FILES}) +if(LibXml2_DIR) + find_package(LibXml2 CONFIG REQUIRED) + target_compile_definitions(svgfill PRIVATE ${LIBXML2_DEFINITIONS}) +endif() +target_link_libraries(svgfill ${Boost_LIBRARIES} ${BCRYPT_LIBRARIES} ${LIBXML2_LIBRARIES} ${CGAL_LIBRARIES}) + +add_executable(svgfill_exe src/main.cpp) +target_link_libraries(svgfill_exe svgfill) +set_property(TARGET svgfill_exe PROPERTY OUTPUT_NAME svgfill) +if(WIN32) + # both the library and the executable now result in a file with basename svgfill, + # on linux the the library is prefixed with lib as libsvgfill.a. Windows does not + # have this mechanism, so on windows the linker would be created an import library + # for the executable, also named svgfill.lib. This naming conflict results in: + # LINK : fatal error LNK1149: output filename matches input filename + # This flag tells the linker not to generate an import library and therefore no + # conflict occurs. + target_link_options(svgfill_exe PRIVATE "/NOIMPLIB") +endif() + +install(TARGETS svgfill_exe DESTINATION ${BINDIR}) +install(TARGETS svgfill DESTINATION ${LIBDIR}) +install(FILES ${LIB_H_FILES} DESTINATION ${INCLUDEDIR}) diff --git a/src/svgfill/README.md b/src/svgfill/README.md new file mode 100644 index 0000000000..dc0348af7c --- /dev/null +++ b/src/svgfill/README.md @@ -0,0 +1,65 @@ +svgfill +======= + +An application to fill areas bounded by unconnected lines in SVG. + +Dependencies +------------ + +* [CGAL 2D Arrangements](https://doc.cgal.org/latest/Arrangement_on_surface_2/index.html) GPL +* [SVG++](http://svgpp.org/) Boost software license + +Compilation +----------- + +Installation is shown based on the IfcOpenShell build script directory output: + + git clone --recursive https://github.com/IfcOpenShell/svgfill + cd svgfill + mkdir build + cd build + + # windows + + set IFCOPENSHELL_ROOT=..\path\to\ifcopenshell\directory\ + cmake -DBOOST_ROOT=%IFCOPENSHELL_ROOT%\deps\boost_1_67_0 ^ + -DBOOST_LIBRARYDIR=%IFCOPENSHELL_ROOT%\deps\boost_1_67_0\stage\vs2017-Win32\lib ^ + -DLIBXML2_INCLUDE_DIR=%IFCOPENSHELL_ROOT%\deps\OpenCOLLADA\Externals\LibXML\include ^ + -DLIBXML2_LIBRARIES=%IFCOPENSHELL_ROOT%\deps-vs2017-x86-installed\OpenCOLLADA\lib\opencollada\xml.lib ^ + -DCGAL_INCLUDE_DIR=%IFCOPENSHELL_ROOT%\deps-vs2017-x86-installed\cgal\include ^ + -DCGAL_LIBRARY_DIR=%IFCOPENSHELL_ROOT%\deps-vs2017-x86-installed\cgal\lib ^ + -DGMP_INCLUDE_DIR=%IFCOPENSHELL_ROOT%\deps-vs2017-x86-installed\mpir ^ + -DGMP_LIBRARY_DIR=%IFCOPENSHELL_ROOT%\deps-vs2017-x86-installed\mpir ^ + -DMPFR_INCLUDE_DIR=%IFCOPENSHELL_ROOT%\deps-vs2017-x86-installed\mpfr ^ + -DMPFR_LIBRARY_DIR=%IFCOPENSHELL_ROOT%\deps-vs2017-x86-installed\mpfr ^ + .. + + # nix + + IFCOPENSHELL_INSTALL=~/IfcOpenShell/build/$(uname -s)/$(uname -m)/install + cmake -DBOOST_ROOT=${IFCOPENSHELL_INSTALL}/boost-1.69.0 \ + -DLIBXML2_INCLUDE_DIR=${IFCOPENSHELL_INSTALL}/libxml2-2.9.9/include/libxml2 \ + -DLIBXML2_LIBRARIES=${IFCOPENSHELL_INSTALL}/libxml2-2.9.9/lib/libxml2.a \ + -DCGAL_INCLUDE_DIR=${IFCOPENSHELL_INSTALL}/cgal-5.2/include \ + -DCGAL_LIBRARY_DIR=${IFCOPENSHELL_INSTALL}/cgal-5.2/lib \ + -DGMP_INCLUDE_DIR=${IFCOPENSHELL_INSTALL}/gmp-6.1.2/include \ + -DGMP_LIBRARY_DIR=${IFCOPENSHELL_INSTALL}/gmp-6.1.2/lib \ + -DMPFR_INCLUDE_DIR=${IFCOPENSHELL_INSTALL}/mpfr-3.1.5/include \ + -DMPFR_LIBRARY_DIR=${IFCOPENSHELL_INSTALL}/mpfr-3.1.5/lib \ + .. + +License +------- + +LGPL + +Example +------- + +in: + +![](examples/rects.svg) + +out: + +![](examples/rects_output.svg) diff --git a/src/svgfill/examples/rects.svg b/src/svgfill/examples/rects.svg new file mode 100644 index 0000000000..23f1bc841d --- /dev/null +++ b/src/svgfill/examples/rects.svg @@ -0,0 +1,80 @@ + + + + + + + + image/svg+xml + + + + + + + + + + + + diff --git a/src/svgfill/examples/rects_output.svg b/src/svgfill/examples/rects_output.svg new file mode 100644 index 0000000000..c16c588993 --- /dev/null +++ b/src/svgfill/examples/rects_output.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/src/svgfill/src/arrange_polygons.cpp b/src/svgfill/src/arrange_polygons.cpp new file mode 100644 index 0000000000..7bddc2e41a --- /dev/null +++ b/src/svgfill/src/arrange_polygons.cpp @@ -0,0 +1,1774 @@ +// #define SVGFILL_DEBUG +// #define SVGFILL_MAIN + +#ifndef SVGFILL_MAIN +#include "svgfill.h" +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "graph_2d.h" + +#if CGAL_VERSION_NR >= 1060000000 +#define variant_get std::get_if +#define my_shared_ptr std::shared_ptr +#else +#define variant_get boost::get +#define my_shared_ptr boost::shared_ptr +#endif + +typedef CGAL::Exact_predicates_exact_constructions_kernel K; +typedef CGAL::Polygon_2 Polygon_2; +typedef CGAL::Polygon_with_holes_2 Polygon_with_holes_2; +typedef K::Point_2 Point_2; +typedef K::Segment_2 Segment_2; +typedef std::vector Polygon_list; +typedef CGAL::Arr_segment_traits_2 Traits_2; +typedef typename Traits_2::Point_2 Arr_Point_2; +typedef typename Traits_2::X_monotone_curve_2 Arr_Segment_2; +typedef CGAL::Arrangement_2 Arrangement_2; + +template +using Triangle = std::array, 3>; + +template +Polygon_2 convert_polygon(const CGAL::Polygon_2& poly) { + Polygon_2 exact_poly; + typedef CGAL::Cartesian_converter Converter_Epick_to_Epeck; + Converter_Epick_to_Epeck converter; + for (auto vit = poly.vertices_begin(); vit != poly.vertices_end(); ++vit) { + exact_poly.push_back(converter(*vit)); // Convert each vertex + } + return exact_poly; +} + +template +void remove_close_points(P& p, double eps = 1.e-2) { + std::vector> ps; + 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; + } + // std::cout << "d " << std::sqrt(CGAL::to_double(CGAL::squared_distance(*I, *J))) << std::endl; + if (CGAL::squared_distance(*I, *J) > (eps * eps)) { + ps.push_back(*J); + I = J; + } + if (last) { + break; + } + } + if (ps.size() >= 2 && CGAL::squared_distance(ps.front(), ps.back()) <= eps * eps) { + // Remove the last point if it is too close to the first point + ps.pop_back(); + } + if (ps.size() != p.size()) { + // std::cerr << "Removed " << (p.size() - ps.size()) << " close points from polygon" << std::endl; + p = P(ps.begin(), ps.end()); + } +} + +std::vector create_and_convert_offset_polygon(double offset_distance, const Polygon_2& polygon_) { + auto polygon = polygon_; + if (!polygon.is_counterclockwise_oriented()) { + polygon.reverse_orientation(); + } + + remove_close_points(polygon); + + // Create the offset polygons using Epick kernel + // create_exterior_skeleton_and_offset_polygons_2() + std::vector>> offset_polygons; + + if (offset_distance >= 0.) { + offset_polygons = CGAL::create_exterior_skeleton_and_offset_polygons_2(offset_distance, polygon); + // erase the first outer frame + offset_polygons.erase(offset_polygons.begin()); + offset_polygons.front()->reverse_orientation(); + } else { + offset_polygons = CGAL::create_interior_skeleton_and_offset_polygons_2(-offset_distance, polygon); + } + + // Convert each offset polygon back to the Epeck kernel + std::vector exact_offset_polygons; + for (auto& inexact_poly_ptr : offset_polygons) { + remove_close_points(*inexact_poly_ptr); + Polygon_2 exact_poly = convert_polygon(*inexact_poly_ptr); + exact_offset_polygons.push_back(exact_poly); + } + + return exact_offset_polygons; +} + +template +T take_first_if_single_item(const std::vector& vec) { + if (vec.size() == 0) { + throw std::runtime_error("Expected at least one item"); + } + if (true || vec.size() == 1) { + return vec.front(); + } + throw std::runtime_error("Expected a single item"); +} + +template +boost::optional maybe_take_first_if_single_item(const std::vector& vec) { + if (vec.size() == 0) { + return boost::none; + } + if (true || vec.size() == 1) { + return vec.front(); + } +} + +template +boost::optional subtract_retain_largest(const T& lhs, const T& rhs) { + std::vector result; + boost::optional mp; + + CGAL::difference(lhs, rhs, std::back_inserter(result)); + + std::sort(result.begin(), result.end(), [](const Polygon_with_holes_2& a, const Polygon_with_holes_2& b) { + return a.outer_boundary().area() < b.outer_boundary().area(); + }); + + if (result.size() > 0) { + if (result.front().has_holes()) { + return boost::none; + } + return result.front().outer_boundary(); + } + + 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; + auto curr = circ; + do { + poly.push_back(curr->source()->point()); + } while (++curr != circ); + return poly; +} + +Polygon_with_holes_2 circ_to_poly(typename Arrangement_2::Ccb_halfedge_const_circulator circ, typename Arrangement_2::Inner_ccb_const_iterator a, typename Arrangement_2::Inner_ccb_const_iterator b) +{ + Polygon_with_holes_2 poly(circ_to_poly(circ)); + for (auto it = a; it != b; ++it) { + poly.add_hole(circ_to_poly(*it)); + } + 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; + for (auto& r : polygons) { + auto ps = create_and_convert_offset_polygon(polygon_offset_distance, r); + for (auto& p : ps) { + if (!p.is_simple()) { + /*{ + std::cerr << "["; + bool first = true; + for (auto& pp : r) { + if (!first) { + std::cerr << ","; + } + first = false; + std::cerr << "(" << pp.x() << "," << pp.y() << ")"; + } + std::cerr << "]" << std::endl; + } + { + 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()); + } + + // Perform Boolean union on the offset polygons + std::vector unioned_polygons; + CGAL::join(offset_polygons.begin(), offset_polygons.end(), std::back_inserter(unioned_polygons)); + Polygon_2 fused_removed_close_points = unioned_polygons.front().outer_boundary(); + remove_close_points(unioned_polygons.front().outer_boundary(), polygon_offset_distance); + + // Apply negative offset to get the outer perimeter polygon + auto inner_offset = create_and_convert_offset_polygon( + // Slightly smaller inset distance for non-manifold situs? + -polygon_offset_distance + 1.e-8, + fused_removed_close_points); + + if (inner_offset.size() != 1) { + throw std::runtime_error("Unexpected union outcome - num outer perimiters: " + std::to_string(inner_offset.size())); + } + + 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; + + 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; + } + } + 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 + + typedef CGAL::Box_intersection_d::Box_with_handle_d Box; + std::vector>> input_triangulated; + + std::vector boxes; + for (auto it = input_polygons.begin(); it != input_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"); + } + + CGAL::Polygon_triangulation_decomposition_2 decompositor; + std::vector temp; + decompositor(*it, std::back_inserter(temp)); + input_triangulated.emplace_back(); + for (auto& pol : temp) { + auto it = pol.vertices_circulator(); + const auto& p = *(it++); + const auto& q = *(it++); + const auto& r = *(it++); + input_triangulated.back().emplace_back(p, q, r); + } + } + + std::set> overlaps; + + CGAL::box_self_intersection_d(boxes.begin(), boxes.end(), [&input_triangulated, &overlaps](const Box& a, const Box& b) { + for (auto& t1 : input_triangulated[a.handle()]) { + 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() + }); + registered_overlap = true; + break; + } + } + if (registered_overlap) { + // no need to check other triangles + break; + } + } + }); + + if (true) { + // solve overlaps by means of subtraction + // loop over overlaps and subtract the smaller polygon from the larger one + + std::set eliminated_polies; + std::map overlap_counts; + for (auto& p : overlaps) { + overlap_counts[p.first]++; + overlap_counts[p.second]++; + } + + 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); + } + } + + /* + 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; + } + */ + + { + // [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(); + } + } + } + } + +#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); + } +#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; + + 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()) + { + // Find closest point & corresponding segment + 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()))); + if (d < (tol * tol)) { + return seg_to_poly.find(seg_ptr)->second; + } + return input_polygons.end(); + }; + + /* + 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) { + CGAL::Point_2 closest; + double closest_distance = std::numeric_limits::infinity(); + auto input_it = input_polygons.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& p : *it) { + auto d = std::sqrt(CGAL::to_double(CGAL::squared_distance(P, p))); + if (d < closest_distance) { + closest_distance = d; + closest = p; + input_it = it; + } + } + } + + return std::make_pair(input_it, closest); + }; + + auto project_input_point = [&input_polygons](const CGAL::Point_2& P) { + CGAL::Point_2 closest; + typename K::FT closest_sq_distance = std::numeric_limits::infinity(); + auto input_it = input_polygons.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 jt = it->edges_begin(); jt != it->edges_end(); ++jt) { + auto Pp = jt->supporting_line().projection(P); + auto d = CGAL::squared_distance(Pp, P); + if (d < closest_sq_distance) { + closest_sq_distance = d; + closest = Pp; + input_it = it; + } + } + } + + 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(); + } + + // 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()) { + /*{ + 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::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()); + } + } + +#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()); + } + + // 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 + + // 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, Point_2> segment_to_midpoint; + std::map> midpoint_to_segment; + std::map*, std::vector>> facet_to_segment; + + for (auto& tri : triangular_polygons) { + for (size_t i = 0; i < 3; ++i) { + size_t j = (i + 1) % 3; + auto& pi = tri.vertex(i); + auto& pj = tri.vertex(j); + const bool orientation = std::lexicographical_compare(pi.cartesian_begin(), pi.cartesian_end(), pj.cartesian_begin(), pj.cartesian_end()); + std::pair seg(orientation ? pi : pj, orientation ? pj : pi); + segment_to_facet[seg].push_back(&tri); + facet_to_segment[&tri].push_back(seg); + } + } + + // 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; + } + }; + + // 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); + + segment_to_input_facet[p.first].push_back(&*p1index); + segment_to_input_facet[p.first].push_back(&*p2index); + + if (p1index != input_polygons.end() && p2index != input_polygons.end() && p1index != p2index) { + 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) { + for (auto& q : segment_to_facet[p.first]) { + for (auto& r : facet_to_segment[q]) { + if (p.first == r) { + continue; + } + 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 + } + } + } + } + + // Find triangles in this network often occuring at junctions in the corridor mesh + std::set> triangles; + std::function&)> find_triangles_recursive; + find_triangles_recursive = [&](std::vector& path) -> void { + // If depth reaches 3, check for a triangle + if (path.size() == 3) { + // Check if we can complete the triangle by going from the current point back to the start + 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] }; + std::sort(triangle.begin(), triangle.end()); + triangles.insert(triangle); + } + return; + } + + // Otherwise, continue exploring neighbors + const std::vector& neighbors = line_graph.at(path.back()); + for (const Point_2& neighbor : neighbors) { + if (std::find(path.begin(), path.end(), neighbor) == path.end()) { + path.push_back(neighbor); + find_triangles_recursive(path); + path.pop_back(); // Backtrack + } + } + }; + + for (auto& p : line_graph) { + 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 + + // non exact time, we need sqrt + 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); + + double global_min_abs_dot = std::numeric_limits::infinity(); + size_t global_min_abs_dot_index; + + for (size_t i = 0; i < 3; ++i) { + auto j = (i + 2) % 3; + auto e0 = st[i] - st[j]; + e0 /= std::sqrt(e0.squared_length()); + + double max_abs_dot = 0.; + + { + auto& ni = line_graph[t[i]]; + for (auto& n : ni) { + if (std::find(t.begin(), t.end(), n) == t.end()) { + // not contained in triangle + auto sn = C(n); + auto en = sn - st[i]; + en /= std::sqrt(en.squared_length()); + auto dot = std::abs(en * e0); + + if (dot > max_abs_dot) { + max_abs_dot = dot; + } + } + } + } + + { + auto& nj = line_graph[t[j]]; + for (auto& n : nj) { + if (std::find(t.begin(), t.end(), n) == t.end()) { + // not contained in triangle + auto sn = C(n); + auto en = sn - st[j]; + en /= std::sqrt(en.squared_length()); + auto dot = std::abs(en * e0); + + if (dot > max_abs_dot) { + max_abs_dot = dot; + } + } + } + + } + + if (max_abs_dot < global_min_abs_dot) { + global_min_abs_dot = max_abs_dot; + global_min_abs_dot_index = i; + } + } + + { + 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 + } + + } + + Graph2D G2(line_graph); + for (auto& e : eliminated_segments) { + G2.remove_edge(e.first, e.second); + } + + 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; + } + } else { + ++vit; + } + } + // std::cout << "Eliminated " << n_vertices_removed << " vertices" << std::endl; + } + + // Ortho edge slide + { + std::list> edges_to_remove, edges_to_insert; + + 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 (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; + } + + 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; + } + } + } + } + } + + 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)); + } + + std::list> move_ops; + std::list> edge_ops; + + 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; + + if (midpoint_to_segment.find(M) == midpoint_to_segment.end()) { + typename K::FT min_sq_distance = std::numeric_limits::infinity(); + for (auto& pa : midpoint_to_segment) { + if (CGAL::squared_distance(pa.first, M) < min_sq_distance) { + q = &pa.second; + min_sq_distance = CGAL::squared_distance(pa.first, M); + } + } + } else { + q = &midpoint_to_segment[M]; + } + + if (q == nullptr) { + continue; + } + + bool handled_as_graph_path = false; + + // distance from unioned - shoot ray? + if (segment_to_input_facet[*q].size() == 2) { + for (auto& bnd : inner_offset) { + // if point M is contained in bnd interior: + if (bnd.has_on_bounded_side(M)) { + auto& incoming = *it->second.begin(); + // create ray incoming -> M + CGAL::Ray_2 ray(incoming, M - incoming); + // intersect ray with boundary + boost::optional> closest_segment; + boost::optional> closest_intersection_point; + K::FT sq_distance_along_ray = std::numeric_limits::infinity(); + for (auto jt = bnd.edges_begin(); jt != bnd.edges_end(); ++jt) { + const auto& seg = *jt; + auto x = CGAL::intersection(ray, seg); + if (x) { + if (auto* xp = variant_get>(&*x)) { + auto dist = ((*xp) - M).squared_length(); + if (dist < sq_distance_along_ray) { + closest_segment = seg; + closest_intersection_point = *xp; + sq_distance_along_ray = dist; + } + } + } + } + + if (closest_intersection_point) { + Graph2D GGG(bnd); + GGG.refine(*GGG.query(*closest_intersection_point, 0.01), *closest_intersection_point); + + std::array>, 2> input_points = { { {}, {} } }; + + size_t i = 0; + for (auto& fac : segment_to_input_facet[*q]) { + for (auto it = fac->vertices_begin(); it != fac->vertices_end(); ++it) { + auto seg = GGG.query(*it, 0.01); + if (seg) { + if (seg->source() != *it && seg->target() != *it) { + GGG.refine(*seg, *it); + } + input_points[i].insert(*it); + } + } + i++; + } + + auto a1 = GGG.shorted_path(*closest_intersection_point, input_points[0]); + auto a2 = GGG.shorted_path(*closest_intersection_point, input_points[1]); + + if (!a1.empty() && !a2.empty()) { + + if (M != *closest_intersection_point) { + edge_ops.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)) }); + } + for (auto it = a2.begin(); it != a2.end() && std::next(it) != a2.end(); ++it) { + edge_ops.push_front({ *it, *(std::next(it)) }); + } + + handled_as_graph_path = true; + break; + } + } + } + } + } + + 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 Q = pq.second; + auto R = pr.second; + + if (Q == R) { + // this can happen in situations like this: + // 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 + // + // ┌───────x───────────────o <---R' + // │ │ + // │ │ + // │ │ + // │ │ + // └───────────────────────┘ + + // @todo is this projection actually necessary or is it already 'exact enough'? + R = 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 }); + + } + } + } + +#ifdef SVGFILL_DEBUG + obj << "o network_4\n"; +#endif + + // 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; + std::set edges_to_remove; + + for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it) { + if (it->is_unbounded()) { + continue; + } + bool is_corridor = false; + { + auto curr = it->outer_ccb(); + do { + auto& p = curr++->source()->point(); + if (G.find(p) != G.vertices_end()) { + is_corridor = true; + break; + } + } while (curr != it->outer_ccb()); + + for (auto jt = it->inner_ccbs_begin(); jt != it->inner_ccbs_end(); ++jt) { + curr = *jt; + do { + auto& p = curr++->source()->point(); + if (G.find(p) != G.vertices_end()) { + is_corridor = true; + break; + } + } while (curr != *jt); + if (is_corridor) { + break; + } + } + } + + if (is_corridor) { + auto curr = it->outer_ccb(); + do { + 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); + const bool on_orig_bound = p1index != input_polygons.end(); + if (on_orig_bound) { + if (edges_to_remove.find(curr->twin()) != edges_to_remove.end()) { + // std::cerr << "Warning trying to delete edge twice" << std::endl; + } else { + edges_to_remove.insert(curr); + } + } + curr++; + } while (curr != it->outer_ccb()); + + for (auto jt = it->inner_ccbs_begin(); jt != it->inner_ccbs_end(); ++jt) { + curr = *jt; + do { + 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); + const bool on_orig_bound = p1index != input_polygons.end(); + if (on_orig_bound) { + if (edges_to_remove.find(curr->twin()) != edges_to_remove.end()) { + // std::cerr << "Warning trying to delete edge twice" << std::endl; + } else { + edges_to_remove.insert(curr); + } + } + curr++; + } while (curr != *jt); + if (is_corridor) { + break; + } + } + } + +#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 + CGAL::remove_edge(arr, e); + } + + 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 +} + +#ifndef SVGFILL_MAIN + +bool svgfill::arrange_polygons(const std::vector& polygons, std::vector& arranged) +{ + std::vector cgal_polygons, cgal_polygons_out; + std::transform(polygons.begin(), polygons.end(), std::back_inserter(cgal_polygons), [](auto& poly) { + Polygon_2 result; + std::transform(poly.boundary.begin(), poly.boundary.end(), std::back_inserter(result), [](auto& p) { + return Point_2(p[0], p[1]); + }); + return result; + }); + arrange_cgal_polygons(cgal_polygons, cgal_polygons_out); + std::transform(cgal_polygons_out.begin(), cgal_polygons_out.end(), std::back_inserter(arranged), [](auto& poly) { + svgfill::polygon_2 result; + std::transform(poly.begin(), poly.end(), std::back_inserter(result.boundary), [](auto& pt) { + return svgfill::point_2{ + CGAL::to_double(pt.cartesian(0)), + CGAL::to_double(pt.cartesian(1)), + }; + }); + return result; + }); + return true; +} + +#else + +template +Polygon_2 create_rectangle(T x_min, T y_min, T x_max, T y_max) { + Polygon_2 rectangle; + rectangle.push_back(Point_2(x_min, y_min)); + rectangle.push_back(Point_2(x_max, y_min)); + rectangle.push_back(Point_2(x_max, y_max)); + rectangle.push_back(Point_2(x_min, y_max)); + return rectangle; +} + +#include + +int main(int argc, char** argv) { + std::vector input_polygons, output; + + if (argc == 2) { + using json = nlohmann::json; + std::ifstream file(argv[1]); + json jsonData; + file >> jsonData; + size_t i = 0; + for (const auto& item : jsonData.items()) { + std::cout << "i " << i << std::endl; + i++; + input_polygons.clear(); + const auto& polygonsData = item.value(); + for (const auto& polygonData : polygonsData) { + input_polygons.emplace_back(); + for (const auto& pointData : polygonData) { + double x = pointData[0]; + double y = pointData[1]; + input_polygons.back().push_back(CGAL::Point_2(x, y)); + } + } + arrange_cgal_polygons(input_polygons, output); + break; + } + return 0; + } else { + Polygon_2 rect1 = create_rectangle(0, 0, 2, 1); + Polygon_2 rect2 = create_rectangle(2.2, 0, 4, 1.1); + Polygon_2 rect3 = create_rectangle(0, 1.2, 2, 4); + Polygon_2 rect4 = create_rectangle(2.2, 1.2, 6, 4); + Polygon_2 rect5 = create_rectangle(4.2, 0, 6, 1.1); + + input_polygons = { rect1, rect2, rect3, rect4, rect5 }; + } + arrange_cgal_polygons(input_polygons, output); + + return 0; +} + +#endif diff --git a/src/svgfill/src/graph_2d.h b/src/svgfill/src/graph_2d.h new file mode 100644 index 0000000000..aaaa059b4e --- /dev/null +++ b/src/svgfill/src/graph_2d.h @@ -0,0 +1,509 @@ +#ifndef GRAPH_2D_H +#define GRAPH_2D_H + +#ifdef SVGFILL_DEBUG +#include +#endif + +template +class Graph2D { +public: + typedef typename Kernel::Point_2 Point_2; + typedef std::pair Edge; + + Graph2D() {} + + Graph2D(const std::map>& input_adjacency_list) { + for (const auto& kv : input_adjacency_list) { + const Point_2& u = kv.first; + const std::vector& neighbors = kv.second; + for (const Point_2& v : neighbors) { + if (u != v) { // no self-edges + adjacency_list[u].insert(v); + adjacency_list[v].insert(u); // Since it's undirected + } + } + } + assert_symmetric(); + } + + Graph2D(const CGAL::Polygon_2& loop) { + for (auto it = loop.vertices_begin(); it != loop.vertices_end(); ++it) { + auto next_it = std::next(it); + if (next_it == loop.vertices_end()) { + next_it = loop.vertices_begin(); + } + adjacency_list[*it].insert(*next_it); + adjacency_list[*next_it].insert(*it); + } + } + + bool is_loop() const { + if (adjacency_list.size() < 3) { + return false; + } + for (const auto& p : adjacency_list) { + if (p.second.size() != 2) { + return false; + } + } + return true; + } + + auto find(const Point_2& p) { + return adjacency_list.find(p); + } + + boost::optional< CGAL::Segment_2 > query(const Point_2& p, typename Kernel::FT eps) { + boost::optional< CGAL::Segment_2 > closest_segment; + typename Kernel::FT closest_distance = std::numeric_limits::infinity(); + for (auto& p1 : adjacency_list) { + for (auto& p2 : p1.second) { + CGAL::Segment_2seg(p1.first, p2); + auto dist = CGAL::squared_distance(p, seg); + if (dist < eps * eps && dist < closest_distance) { + closest_distance = dist; + closest_segment = seg; + } + } + } + return closest_segment; + } + + void refine(const CGAL::Segment_2& seg, const Point_2& p) { + remove_edge(seg.source(), seg.target()); + insert(seg.source(), p); + insert(p, seg.target()); + } + + std::vector shorted_path(const Point_2& start, const Point_2& goal) const { + auto& adj = adjacency_list; + if (adj.count(start) == 0 || adj.count(goal) == 0) return {}; + + std::map predecessor; + std::queue q; + + // seed BFS + predecessor[start] = start; // mark start as "seen" + q.push(start); + + // BFS + bool found = false; + while (!q.empty() && !found) { + Point_2 u = q.front(); q.pop(); + for (auto& v : adj.at(u)) { + // if v has no predecessor yet, it's unseen + if (!predecessor.count(v)) { + predecessor[v] = u; + q.push(v); + if (v == goal) { found = true; break; } + } + } + } + + if (!found) return {}; + + // reconstruct path + std::vector path; + for (Point_2 cur = goal; cur != start; cur = predecessor[cur]) + path.push_back(cur); + path.push_back(start); + std::reverse(path.begin(), path.end()); + return path; + } + + + std::vector shorted_path(const Point_2& start, const std::set& goal) const { + auto& adj = adjacency_list; + if (adj.count(start) == 0) return {}; + + std::map predecessor; + std::queue q; + + // seed BFS + predecessor[start] = start; // mark start as "seen" + q.push(start); + + Point_2 used_goal; + + // BFS + bool found = false; + while (!q.empty() && !found) { + Point_2 u = q.front(); q.pop(); + for (auto& v : adj.at(u)) { + // if v has no predecessor yet, it's unseen + if (!predecessor.count(v)) { + predecessor[v] = u; + q.push(v); + if (goal.find(v) != goal.end()) { + found = true; + used_goal = v; + break; + } + } + } + } + + if (!found) return {}; + + // reconstruct path + std::vector path; + for (Point_2 cur = used_goal; cur != start; cur = predecessor[cur]) + path.push_back(cur); + path.push_back(start); + std::reverse(path.begin(), path.end()); + return path; + } + + + void move(const Point_2& from, const Point_2& to) { + // @todo should check for intersections? + auto it = adjacency_list.find(from); + if (it != adjacency_list.end()) { + auto neighbours = it->second; + for (auto& n : neighbours) { + adjacency_list[n].erase(from); + adjacency_list[n].insert(to); + } + adjacency_list.erase(it); + adjacency_list.insert({ to, neighbours }); + } + } + + bool is_valid() const { + typedef CGAL::Box_intersection_d::Box_with_handle_d Box; + std::vector boxes; + std::vector> segments; + for (const auto& p : adjacency_list) { + for (const auto& q : p.second) { + if (p.first < q) { + segments.emplace_back(p.first, q); + } + } + } + for (auto it = segments.begin(); it != segments.end(); ++it) { + boxes.emplace_back(it->bbox(), std::distance(segments.begin(), it)); + } + bool any = false; + CGAL::box_self_intersection_d(boxes.begin(), boxes.end(), [this, &segments, &any](const Box& a, const Box& b) { + auto& seg1 = segments[a.handle()]; + auto& seg2 = segments[b.handle()]; + // Skip topologically connected segments + if (seg1.source() == seg2.source() || seg1.source() == seg2.target() || seg1.target() == seg2.source() || seg1.target() == seg2.target()) { + return; + } + if (CGAL::do_intersect(seg1, seg2)) { + any = true; + } + }); + return any; + } + + // Eliminates a vertex with exactly two neighbors by connecting its neighbors + typename std::map>::iterator eliminate_vertex(typename std::map>::iterator it) { + if (it == adjacency_list.end()) { + // Vertex not found + return adjacency_list.end(); + } + const std::set& neighbors = it->second; + if (neighbors.size() != 2) { + // Not exactly two neighbors + return adjacency_list.end(); + } + // Get the two neighbors + auto neighbor_it = neighbors.begin(); + auto u = *neighbor_it++; + auto w = *neighbor_it; + auto v = it->first; + + // Remove all edge between v and u + adjacency_list[u].erase(v); + adjacency_list[v].erase(u); + + // Remove all edge between v and w + adjacency_list[w].erase(v); + adjacency_list[v].erase(w); + + // Add edge between u and w + adjacency_list[u].insert(w); + adjacency_list[w].insert(u); + + // Remove v from adjacency_list + auto jt = adjacency_list.erase(it); + + assert_symmetric(); + + return jt; + } + + Graph2D weld_vertices() const { + std::set points; + for (auto& p : adjacency_list) { + points.insert(p.first); + for (auto& q : p.second) { + points.insert(q); + } + } + + using It = typename std::set::iterator; + using Box = CGAL::Box_intersection_d::Box_with_handle_d; + + std::vector boxes; + for (auto it = points.begin(); it != points.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), + &*it + ); + } + + std::vector> overlaps; + + CGAL::box_self_intersection_d(boxes.begin(), boxes.end(), [&overlaps](const Box& a, const Box& b) { + overlaps.emplace_back(a.handle(), b.handle()); + }); + + std::map> input_adjacency_list; + + { + std::map> adj; + for (const auto& edge : overlaps) { + adj[edge.first].push_back(edge.second); + adj[edge.second].push_back(edge.first); + } + for (auto& p : points) { + adj[&p]; + } + + std::map visited; + std::vector> connected_components; + + for (auto& p : adj) { + if (!visited[p.first]) { + connected_components.emplace_back(); + + std::stack stack; + stack.push(p.first); + visited[p.first] = true; + + while (!stack.empty()) { + auto u = stack.top(); + stack.pop(); + connected_components.back().push_back(u); + + for (auto& neighbor : adj[u]) { + if (!visited[neighbor]) { + visited[neighbor] = true; + stack.push(neighbor); + } + } + } + } + } + + std::map mapping; + + for (auto& comp : connected_components) { + Point_2 avg(0, 0); + for (auto& c : comp) { + avg += (*c - CGAL::ORIGIN); + } + avg = CGAL::ORIGIN + ((avg - CGAL::ORIGIN) / comp.size()); + + for (auto& c : comp) { + mapping[*c] = avg; + } + } + + for (auto& comp : connected_components) { + for (auto& c : comp) { + const auto& C = mapping[*c]; + for (auto& n : adjacency_list.find(*c)->second) { + const auto& N = mapping[n]; + if (C != N) { + if (std::find(input_adjacency_list[C].begin(), input_adjacency_list[C].end(), N) == input_adjacency_list[C].end()) { + input_adjacency_list[C].push_back(N); + } + } + } + } + } + } + + return Graph2D(input_adjacency_list); + } + + void assert_symmetric() { +#ifdef SVGFILL_DEBUG +#if 0 + for (auto& p : adjacency_list) { + if (p.second.find(p.first) != p.second.end()) { + std::cout << "!! " << p.first << " self-edge" << std::endl; + throw std::runtime_error("self-edge"); + } + } + nlohmann::json json_obj = nlohmann::json::array(); + for (auto& p : adjacency_list) { + if (p.second.size() == 0) { + continue; + } + + nlohmann::json pair = nlohmann::json::array(); + nlohmann::json coord = nlohmann::json::array(); + nlohmann::json values = nlohmann::json::array(); + coord.push_back( + CGAL::to_double(p.first.x()) + ); + coord.push_back( + CGAL::to_double(p.first.y()) + ); + pair.push_back(coord); + for (auto& q : p.second) { + auto& r = adjacency_list[q]; + for (auto& s : r) { + nlohmann::json coord = nlohmann::json::array(); + coord.push_back( + CGAL::to_double(s.x()) + ); + coord.push_back( + CGAL::to_double(s.y()) + ); + values.push_back(coord); + } + if (r.find(p.first) == r.end()) { + for (auto& s : r) { + std::cout << " " << s << std::endl; + } + std::cout << "!! " << p.first << " not found in neighbours of " << q << std::endl; + throw std::runtime_error("internal error"); + } + } + pair.push_back(values); + json_obj.push_back( + pair + ); + } + + static int fff = 0; + std::ofstream file("graph_" + std::to_string(fff++) + ".json"); + file << json_obj.dump(2); +#endif +#endif + } + + void remove_edge(const Point_2& u, const Point_2& v) { + adjacency_list[u].erase(v); + adjacency_list[v].erase(u); + assert_symmetric(); + } + + void insert(const Point_2& u, const Point_2& v) { + adjacency_list[u].insert(v); + adjacency_list[v].insert(u); + assert_symmetric(); + } + + // Iterators over vertices + typedef typename std::map>::const_iterator vertex_const_iterator; + typedef typename std::map>::iterator vertex_iterator; + + vertex_const_iterator vertices_begin() const { + return adjacency_list.cbegin(); + } + + vertex_const_iterator vertices_end() const { + return adjacency_list.cend(); + } + + vertex_iterator vertices_begin() { + return adjacency_list.begin(); + } + + vertex_iterator vertices_end() { + return adjacency_list.end(); + } + + // Edge iterator class + class EdgeIterator { + public: + typedef std::forward_iterator_tag iterator_category; + typedef Edge value_type; + typedef ptrdiff_t difference_type; + typedef const Edge* pointer; + typedef const Edge& reference; + + EdgeIterator() : outer_it_(), inner_it_(), graph_(nullptr) {} + EdgeIterator(const Graph2D* graph, typename std::map>::const_iterator outer_it) + : outer_it_(outer_it), graph_(graph) { + if (outer_it_ != graph_->adjacency_list.end()) { + inner_it_ = outer_it_->second.begin(); + advance_to_valid(); + } + } + + reference operator*() const { + current_edge_ = Edge(outer_it_->first, *inner_it_); + return current_edge_; + } + + pointer operator->() const { + current_edge_ = Edge(outer_it_->first, *inner_it_); + return ¤t_edge_; + } + + EdgeIterator& operator++() { + ++inner_it_; + advance_to_valid(); + return *this; + } + + EdgeIterator operator++(int) { + EdgeIterator tmp = *this; + ++(*this); + return tmp; + } + + bool operator==(const EdgeIterator& other) const { + return outer_it_ == other.outer_it_ && (outer_it_ == graph_->adjacency_list.end() || inner_it_ == other.inner_it_); + } + + bool operator!=(const EdgeIterator& other) const { + return !(*this == other); + } + + private: + void advance_to_valid() { + while (outer_it_ != graph_->adjacency_list.end()) { + while (inner_it_ != outer_it_->second.end() && *inner_it_ < outer_it_->first) { + ++inner_it_; + } + if (inner_it_ != outer_it_->second.end()) { + break; + } + ++outer_it_; + if (outer_it_ != graph_->adjacency_list.end()) { + inner_it_ = outer_it_->second.begin(); + } + } + } + + mutable Edge current_edge_; + typename std::map>::const_iterator outer_it_; + typename std::set::const_iterator inner_it_; + const Graph2D* graph_; + }; + + EdgeIterator edges_begin() const { + return EdgeIterator(this, adjacency_list.begin()); + } + + EdgeIterator edges_end() const { + return EdgeIterator(this, adjacency_list.end()); + } + +private: + std::map> adjacency_list; +}; + +#endif diff --git a/src/svgfill/src/main.cpp b/src/svgfill/src/main.cpp new file mode 100644 index 0000000000..93ecfb6462 --- /dev/null +++ b/src/svgfill/src/main.cpp @@ -0,0 +1,127 @@ +/**************************************************************************** + * SVG fill * + * * + * Copyright(C) 2020 AECgeeks and Bimforce * + * * + * This program is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Lesser General Public * + * License as published by the Free Software Foundation; either * + * version 3 of the License, or (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * Lesser General Public License for more details. * + * * + * You should have received a copy of the GNU Lesser General Public License * + * along with this program; if not, write to the Free Software Foundation, * + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * + ****************************************************************************/ + +#include "svgfill.h" +#include "progress.h" + +#include +#include +#include + +#include +#include +#include +#include +#include + +int main(int argc, char** argv) { + bool valid_command_line = false; + bool random_color = false; + double eps = 1.e-5; + boost::optional class_name; + + std::vector flags; + std::vector args; + svgfill::solver s = svgfill::FILTERED_CARTESIAN_QUOTIENT; + std::map solver_mapping { + {"cartesian_double", svgfill::CARTESIAN_DOUBLE}, + {"cartesian_quotient", svgfill::CARTESIAN_QUOTIENT}, + {"filtered_cartesian_quotient", svgfill::FILTERED_CARTESIAN_QUOTIENT}, + {"exact_predicates", svgfill::EXACT_PREDICATES}, + {"exact_constructions", svgfill::EXACT_CONSTRUCTIONS}, + }; + progress_bar::style progress_style = progress_bar::BAR; + + for (int i = 1; i < argc; ++i) { + std::string a = argv[i]; + if (boost::starts_with(a, "-")) { + flags.push_back(a); + } + else { + args.push_back(a); + } + } + std::string fn, ofn; + + if (args.size() == 2) { + fn = args[0]; + ofn = args[1]; + valid_command_line = true; + } + + for (auto& f : flags) { + if (f == "--random-color") { + random_color = true; + } + else if (f == "-q") { + progress_style = progress_bar::DOTS; + } + else if (boost::starts_with(f, "--class=")) { + class_name = f.substr(strlen("--class=")); + } + else if (boost::starts_with(f, "--solver=")) { + std::string solver_str = f.substr(strlen("--solver=")); + auto it = solver_mapping.find(solver_str); + if (it == solver_mapping.end()) { + valid_command_line = false; + } + else { + s = it->second; + } + } + else if (boost::starts_with(f, "--eps=")) { + std::string eps_str = f.substr(strlen("--eps=")); + eps = boost::lexical_cast(eps_str); + } + else { + valid_command_line = false; + } + } + + if (!valid_command_line) { + std::cerr << "Usage: " << argv[0] << " [--random-color] [--class=...] " << std::endl; + return 1; + } + + std::vector> segments; + std::vector> polygons; + + progress_bar p(std::cout, progress_style); + application_progress ap({1., 10., 1.}, p); + std::function pfn = [&ap](float f) { ap(f); }; + + std::ifstream fs(fn.c_str()); + std::string data(std::istreambuf_iterator{fs}, {}); + fs.close(); + + if (!svgfill::svg_to_line_segments(data, class_name, segments)) { + return 1; + } + ap.finished(); + + if (!svgfill::line_segments_to_polygons(s, eps, segments, polygons, pfn)) { + return 1; + } + + ap.finished(); + + std::ofstream ofs(ofn.c_str()); + ofs << svgfill::polygons_to_svg(polygons, random_color); +} diff --git a/src/svgfill/src/progress.h b/src/svgfill/src/progress.h new file mode 100644 index 0000000000..daf54a3e59 --- /dev/null +++ b/src/svgfill/src/progress.h @@ -0,0 +1,122 @@ +/*************************************************************************** +/* * +/* Copyright 2021 AECgeeks * +/* * +/* Permission is hereby granted, free of charge, to any person obtaining a * +/* copy of this software and associated documentation files (the * +/* "Software"), to deal in the Software without restriction, including * +/* without limitation the rights to use, copy, modify, merge, publish, * +/* distribute, sublicense, and/or sell copies of the Software, and to * +/* permit persons to whom the Software is furnished to do so, subject to * +/* the following conditions: * +/* * +/* The above copyright notice and this permission notice shall be included * +/* in all copies or substantial portions of the Software. * +/* * +/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS * +/* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF * +/* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. * +/* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY * +/* CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, * +/* TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE * +/* SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. * +/* * +/***************************************************************************/ + +// A trivial C++ progress bar + +#ifndef PROGRESS_H +#define PROGRESS_H + +#include +#include +#include +#include +#include +#include + +class progress_bar { + std::ostream& s_; + float max_; + size_t width_; + size_t* last_emitted_p_ = nullptr; + +public: + enum style { + BAR, DOTS + }; + +private: + style style_; + +public: + progress_bar(std::ostream& s = std::cerr, style st = BAR, float max = 1., size_t width = 50) + : s_(s) + , max_(max) + , width_(st == BAR ? width : 100U) + , style_(st) + {} + + void operator()(size_t p) { + if (last_emitted_p_ && p <= *last_emitted_p_) { + return; + } + + p = p > width_ ? width_ : p; + if (style_ == BAR) { + s_ << "\r[" + std::string(p, '#') + std::string(width_ - p, ' ') + "]" << std::flush; + } else { + s_ << std::string(p - (last_emitted_p_ ? *last_emitted_p_ : 0U), '.') << std::flush; + } + + if (last_emitted_p_) { + *last_emitted_p_ = p; + } + else { + last_emitted_p_ = new size_t(p); + } + } + + void operator()(float p) { + (*this)((size_t) (p / max_ * width_)); + } + + ~progress_bar() { + delete last_emitted_p_; + } +}; + +class application_progress { + std::vector estimates_; + size_t phase_ = 0; + std::function callback_; + float total_; + +public: + void operator()(float p) { + auto progress = std::accumulate(estimates_.begin(), estimates_.begin() + phase_, 0.f); + progress += p * estimates_[phase_]; + callback_(progress / total_); + } + + application_progress(const std::vector& estimates, const std::function& callback) + : estimates_(estimates) + , callback_(callback) + { + total_ = std::accumulate(estimates_.begin(), estimates_.end(), 0.f); + (*this)(0.); + } + + void finished() { + ++phase_; + (*this)(0.); + } + + ~application_progress() { + phase_ = estimates_.size() - 2; + finished(); + } + +}; + +#endif \ No newline at end of file diff --git a/src/svgfill/src/svgfill.cpp b/src/svgfill/src/svgfill.cpp new file mode 100644 index 0000000000..8a2a1bb008 --- /dev/null +++ b/src/svgfill/src/svgfill.cpp @@ -0,0 +1,633 @@ +/**************************************************************************** + * SVG fill * + * * + * Copyright(C) 2020 AECgeeks and Bimforce * + * * + * This program is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Lesser General Public * + * License as published by the Free Software Foundation; either * + * version 3 of the License, or (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * Lesser General Public License for more details. * + * * + * You should have received a copy of the GNU Lesser General Public License * + * along with this program; if not, write to the Free Software Foundation, * + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * + ****************************************************************************/ + +#include "svgfill.h" + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace svgpp; + +class Context +{ +private: + size_t depth_ = 0; + int enabled_at_ = -1; + svgfill::point_2 start_, xy_; + +public: + boost::optional class_name; + std::vector> segments; + + void on_enter_element(tag::element::any) + { + ++depth_; + } + + void on_enter_element(tag::element::g) + { + ++depth_; + if (enabled_at_ == -1 && !class_name.is_initialized()) { + enabled_at_ = depth_; + segments.emplace_back(); + } + } + + void on_exit_element() + { + if (depth_-- == enabled_at_) { + enabled_at_ = -1; + } + } + + template + void set(tag::attribute::id, Str const & s) { + } + + template + void set(tag::attribute::class_, Str const & s) { + if (enabled_at_ == -1 && class_name.is_initialized() && std::string(s.begin(), s.size()).find(*class_name) != std::string::npos) { + enabled_at_ = depth_; + segments.emplace_back(); + } + } + + void transform_matrix(const boost::array & matrix) + {} + + void path_move_to(double x, double y, tag::coordinate::absolute) + { + start_ = xy_ = { x, y }; + } + + void path_line_to(double x, double y, tag::coordinate::absolute) + { + if (enabled_at_ != -1) { + svgfill::point_2 next{ x, y }; + segments.back().push_back({ xy_, next }); + xy_ = next; + } + } + + void path_cubic_bezier_to( + double x1, double y1, + double x2, double y2, + double x, double y, + tag::coordinate::absolute) {} + + void path_quadratic_bezier_to( + double x1, double y1, + double x, double y, + tag::coordinate::absolute) {} + + void path_elliptical_arc_to( + double rx, double ry, double x_axis_rotation, + bool large_arc_flag, bool sweep_flag, + double x, double y, + tag::coordinate::absolute) {} + + void path_close_subpath() { + if (enabled_at_ != -1) { + segments.back().push_back({ xy_, start_ }); + } + } + + void path_exit() {} +}; + +typedef +boost::mpl::set< + // SVG Structural Elements + tag::element::svg, + tag::element::g, + // SVG Shape Elements + tag::element::circle, + tag::element::ellipse, + tag::element::line, + tag::element::path, + tag::element::polygon, + tag::element::polyline, + tag::element::rect +>::type processed_elements_t; + +// This cryptic code just merges predefined sequences traits::shapes_attributes_by_element +// and traits::viewport_attributes with tag::attribute::transform and tag::attribute::xlink::href +// attributes into single MPL sequence +typedef +boost::mpl::fold< + boost::mpl::protect< + traits::shapes_attributes_by_element + >, + boost::mpl::set< + tag::attribute::id, + tag::attribute::class_ + >::type, + boost::mpl::insert +>::type processed_attributes_t; + +bool svgfill::svg_to_line_segments(const std::string& data, const boost::optional& class_name, std::vector>& segments) +{ + Context context; + context.class_name = class_name; + + xmlDoc* doc = xmlReadMemory(data.c_str(), data.size(), nullptr, nullptr, 0); + xmlNode* elem = xmlDocGetRootElement(doc); + + try { + document_traversal< + processed_elements, + processed_attributes + >::load_document(elem, context); + } + catch (std::exception& e) { + std::cerr << e.what() << std::endl; + return false; + } + + segments = context.segments; + + return true; +} + +bool svgfill::line_segments_to_polygons(solver s, double eps, const std::vector>& segments, std::vector>& polygons) +{ + std::function fn = [](float f) {}; + return line_segments_to_polygons(s, eps, segments, polygons, fn); +} + +bool svgfill::svg_to_polygons(const std::string& data, const boost::optional& class_name, std::vector& polygons) { + Context context; + context.class_name = class_name; + xmlDoc* doc = xmlReadMemory(data.c_str(), data.size(), nullptr, nullptr, 0); + xmlNode* elem = xmlDocGetRootElement(doc); + + try { + document_traversal< + processed_elements, + processed_attributes + >::load_document(elem, context); + } catch (std::exception& e) { + std::cerr << e.what() << std::endl; + return false; + } + std::function fn = [](float f) {}; + std::vector> ps; + if (!line_segments_to_polygons(svgfill::EXACT_PREDICATES, 0., context.segments, ps, fn)) { + return false; + } + if (ps.empty()) { + return false; + } + for (auto& p : ps) { + polygons.insert(polygons.end(), p.begin(), p.end()); + } + return true; +} + +template +class cgal_arrangement : public svgfill::abstract_arrangement { + typedef CGAL::Arr_segment_traits_2 Traits_2; + typedef typename Traits_2::Point_2 Point_2; + typedef typename Traits_2::X_monotone_curve_2 Segment_2; + typedef CGAL::Arrangement_2 Arrangement_2; + typedef CGAL::Polygon_2 Polygon_2; + typedef CGAL::Polygon_with_holes_2 Polygon_wh_2; + typedef typename Arrangement_2::Inner_ccb_const_iterator Inner_ccb_const_iterator; + typedef typename Arrangement_2::Ccb_halfedge_const_circulator Ccb_halfedge_const_circulator; + typedef typename Arrangement_2::Halfedge_handle Halfedge_handle; + typedef typename Arrangement_2::Face_handle Face_handle; + + Polygon_2 circ_to_poly(Ccb_halfedge_const_circulator circ) + { + Polygon_2 poly; + auto curr = circ; + do { + if (poly.size() == 0 || (*(poly.end() - 1)) != curr->source()->point()) { + poly.push_back(curr->source()->point()); + } + } while (++curr != circ); + return poly; + } + + CGAL::Triangle_2 poly_to_triangle(const Polygon_2& poly) + { + auto n = std::distance(poly.vertices_begin(), poly.vertices_end()); + if (n != 3) { + throw std::runtime_error("Unexpected number of points in polygon"); + } + auto p = *poly.vertices_begin(); + auto q = *next(poly.vertices_begin(), 1); + auto r = *next(poly.vertices_begin(), 2); + return CGAL::Triangle_2(p, q, r); + } + + Polygon_wh_2 circ_to_poly(Ccb_halfedge_const_circulator circ, Inner_ccb_const_iterator a, Inner_ccb_const_iterator b) + { + Polygon_wh_2 poly(circ_to_poly(circ)); + for (auto it = a; it != b; ++it) { + poly.add_hole(circ_to_poly(*it)); + } + return poly; + } + + svgfill::point_2 create_point(const Point_2& pt) + { + return svgfill::point_2{ + CGAL::to_double(pt.cartesian(0)), + CGAL::to_double(pt.cartesian(1)), + }; + } + + void set_point_inside(const Polygon_wh_2& inpoly, svgfill::polygon_2& outpoly) + { + /* + std::cout << std::endl; + for (auto& p : inpoly.outer_boundary()) { + std::cout << " " << p; + } + std::cout << std::endl; + */ + // create Delaunay triangulation and return the centroid of the largest triangle. + CGAL::Polygon_triangulation_decomposition_2 decompositor; + std::list decom_polies; + decompositor(inpoly, std::back_inserter(decom_polies)); + decom_polies.sort([](const Polygon_2& a, const Polygon_2& b) { + return a.area() > b.area(); + }); + if (!decom_polies.empty()) { + const Polygon_2& largest = decom_polies.front(); + auto triangle = poly_to_triangle(largest); + /* + for (auto& p : decom_polies) { + std::cout << "a " << CGAL::to_double(poly_to_triangle(largest).area()) << std::endl; + } + std::cout << "triangle area " << CGAL::to_double(triangle.area()) << std::endl; + */ + outpoly.point_inside = create_point(CGAL::centroid(triangle)); + } + } + + Arrangement_2 arr; + float total, i; + +public: + bool operator()(double eps, const std::vector& segments, std::function& progress) { + i = 0; + total = segments.size() + segments.size() / 2; + + for (auto& l : segments) { + Point_2 a(l[0][0], l[0][1]); + Point_2 b(l[1][0], l[1][1]); + if (a == b) { + continue; + } + if (eps != 0.) { + auto ab = b - a; + ab /= std::sqrt(CGAL::to_double(ab.squared_length())); + // This appears to work better generally, slightly nudge the + // end points to make sure segments intersect. + a -= ab * eps; + b += ab * eps; + } + Segment_2 seg(a, b); + CGAL::insert(arr, seg); + + if (progress) { + progress(i++ / total); + } + } + + return true; + } + + void remove_duplicates(svgfill::loop_2& l) { + auto norm2 = [](auto& a, auto& b) { + auto dx = a[0] - b[0]; + auto dy = a[1] - b[1]; + return std::sqrt(dx * dx + dy * dy); + }; + + while (l.size() > 1 && l.front() == l.back()) { + l.pop_back(); + } + + if (l.size() > 1) { + auto it = l.begin(); + auto next_it = std::next(it); + + while (next_it != l.end()) { + if (norm2(*it, *next_it) < 1.e-8) { + next_it = l.erase(next_it); + // 'it' remains the same; 'next_it' now points to the next element + } else { + // Move both iterators forward + ++it; + ++next_it; + } + } + } + } + + bool write(std::vector& polygons, std::function& progress) { + std::vector ps; + ps.reserve(std::distance(arr.faces_begin(), arr.faces_end())); + + for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it) { + const auto& f = *it; + if (!f.is_unbounded()) { + ps.push_back(circ_to_poly( + f.outer_ccb(), + f.inner_ccbs_begin(), + f.inner_ccbs_end() + )); + } + + if (progress) { + progress(i++ / total); + } + } + + // Sort polygons (only taking into account outer boundary) to have inner + // loops drawn over outer boundaries. In SVG draw order is defined by + // position in the tree. + // @nb we do now add the inner boundaries to the path as well. + /* + std::sort(ps.begin(), ps.end(), [](const Polygon_wh_2& a, const Polygon_wh_2& b) { + return a.outer_boundary().area() > b.outer_boundary().area(); + }); + */ + + polygons.reserve(ps.size()); + + + + std::transform(ps.begin(), ps.end(), std::back_inserter(polygons), [this, &progress](const Polygon_wh_2& p) { + svgfill::polygon_2 p2; + std::transform(p.outer_boundary().vertices_begin(), p.outer_boundary().vertices_end(), std::back_inserter(p2.boundary), [this](const Point_2& pt) { + return create_point(pt); + }); + + /* + static int NN = 0; + std::cout << NN++ << std::endl; + for (auto& p : p2.boundary) { + std::cout << std::setprecision(20) << p[0] << "," << p[1] << " "; + } + std::cout << std::endl; + */ + + // duplicates need to be removed after conversion from epeck to double + remove_duplicates(p2.boundary); + + /* + std::cout << "> " << std::endl; + for (auto& p : p2.boundary) { + std::cout << std::setprecision(20) << p[0] << "," << p[1] << " "; + } + std::cout << std::endl; + */ + + std::transform(p.holes_begin(), p.holes_end(), std::back_inserter(p2.inner_boundaries), [this](const Polygon_2& poly) { + svgfill::loop_2 lp; + std::transform(poly.vertices_begin(), poly.vertices_end(), std::back_inserter(lp), [this](const Point_2& pt) { + return create_point(pt); + }); + remove_duplicates(lp); + return lp; + }); + set_point_inside(p, p2); + + if (progress) { + progress(i++ / total); + } + + return p2; + }); + + return true; + } + + std::vector get_face_pairs() { + std::vector ps; + ps.reserve(arr.number_of_edges() * 2); + size_t n = 0; + + std::map face_to_bounded_index; + for (auto it = arr.faces_begin(); it != arr.faces_end(); ++it) { + if (!it->is_unbounded()) { + face_to_bounded_index[it] = face_to_bounded_index.size(); + } + } + + for (auto it = arr.edges_begin(); it != arr.edges_end(); ++it, ++n) { + auto v0 = it->source()->point(); + double v0x = CGAL::to_double(v0.cartesian(0)); + double v0y = CGAL::to_double(v0.cartesian(1)); + auto v1 = it->target()->point(); + double v1x = CGAL::to_double(v1.cartesian(0)); + double v1y = CGAL::to_double(v1.cartesian(1)); + double l = std::sqrt((v1x - v0x) * (v1x - v0x) + (v1y - v0y) * (v1y - v0y)); + // std::cout << n << " l " << l << std::endl; + bool emitted = false; + if (l > 1.) { + // std::cout << std::to_string(n) << " " << v0x << " " << v0y << std::endl; + // std::cout << std::string(std::to_string(n).size(), ' ') << " " << v1x << " " << v1y << std::endl; + + auto afit = face_to_bounded_index.find(it->face()); + auto bfit = face_to_bounded_index.find(it->twin()->face()); + + if (afit != face_to_bounded_index.end() && bfit != face_to_bounded_index.end()) { + ps.push_back(afit->second); + ps.push_back(bfit->second); + emitted = true; + } + } + + if (!emitted) { + ps.push_back(-1); + ps.push_back(-1); + } + } + return ps; + } + + void merge(const std::vector& edge_indices) { + if (edge_indices.empty()) { + return; + } + + std::list to_remove; + auto eit = edge_indices.begin(); + size_t n = 0; + for (auto it = arr.edges_begin(); it != arr.edges_end(); ++it, ++n) { + if (n == *eit) { + ++eit; + to_remove.push_back(it); + if (eit == edge_indices.end()) { + break; + } + } + } + + for (auto& h : to_remove) { + /* + auto v0 = h->source()->point(); + std::cout << CGAL::to_double(v0.cartesian(0)) << " " << CGAL::to_double(v0.cartesian(1)) << std::endl; + auto v1 = h->target()->point(); + std::cout << CGAL::to_double(v1.cartesian(0)) << " " << CGAL::to_double(v1.cartesian(1)) << std::endl << std::endl; + */ + arr.remove_edge(h); + } + } + + size_t num_edges() { + return arr.number_of_edges(); + } + + size_t num_faces() { + return arr.number_of_faces(); + } +}; + +bool svgfill::line_segments_to_polygons(solver s, double eps, const std::vector>& segment_groups, std::vector>& polygons, std::function& progress) +{ + bool b = false; + for (auto& segments : segment_groups) { + context ctx(s, eps, progress); + ctx.add(segments); + if (ctx.build()) { + ctx.write(polygons); + b = true; + } + } + return b; +} + +namespace { + std::string format_pt(const svgfill::point_2& p) { + std::ostringstream oss; + oss << std::setprecision(std::numeric_limits::max_digits10); + oss << p[0] << "," << p[1]; + return oss.str(); + } + + std::string format_poly(const svgfill::loop_2& p) { + std::ostringstream oss; + for (auto it = p.begin(); it != p.end(); ++it) { + oss << ((it == p.begin()) ? "M" : " L"); + oss << format_pt(*it); + } + oss << " Z"; + return oss.str(); + } +} + +std::string svgfill::polygons_to_svg(const std::vector>& polygons, bool random_color) { + std::random_device rd; + std::mt19937 mt(rd()); + std::uniform_int_distribution dist(0, 360); + + std::ostringstream oss; + + oss << ""; + oss << ""; + + for (auto& g : polygons) { + oss << ""; + for (auto& p : g) { + const int h = dist(mt); + const int s = 50; + const int l = 50; + std::string style; + if (random_color) { + std::ostringstream oss; + oss << "style = \"fill: hsl(" << h << "," << s << "%, " << l << "%)\""; + style = oss.str(); + } + oss << ""; + } + oss << ""; + } + + oss << ""; + + return oss.str(); +} + + +std::string svgfill::polygons_to_svg(const std::vector& polygons, bool random_color) { + std::vector> pps = { polygons }; + return polygons_to_svg(pps, random_color); +} + +void svgfill::context::add(const std::vector& segments) { + segments_.insert(segments_.end(), segments.begin(), segments.end()); +} + +bool svgfill::context::build() { + if (solver_ == CARTESIAN_DOUBLE) { + arr_ = new cgal_arrangement>; + } else if (solver_ == CARTESIAN_QUOTIENT) { + arr_ = new cgal_arrangement>>; + } else if (solver_ == FILTERED_CARTESIAN_QUOTIENT) { + arr_ = new cgal_arrangement>>>; + } else if (solver_ == EXACT_PREDICATES) { + arr_ = new cgal_arrangement; + } else if (solver_ == EXACT_CONSTRUCTIONS) { + arr_ = new cgal_arrangement; + } + return (*arr_)(eps_, segments_, progress_); +} + +void svgfill::context::merge(const std::vector& edge_indices) { + arr_->merge(edge_indices); +} + +void svgfill::context::write(std::vector>& p) { + std::vector polygons; + arr_->write(polygons, progress_); + p.push_back(polygons); +} diff --git a/src/svgfill/src/svgfill.h b/src/svgfill/src/svgfill.h new file mode 100644 index 0000000000..396fc9c924 --- /dev/null +++ b/src/svgfill/src/svgfill.h @@ -0,0 +1,119 @@ +/**************************************************************************** + * SVG fill * + * * + * Copyright(C) 2020 AECgeeks and Bimforce * + * * + * This program is free software; you can redistribute it and/or * + * modify it under the terms of the GNU Lesser General Public * + * License as published by the Free Software Foundation; either * + * version 3 of the License, or (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * Lesser General Public License for more details. * + * * + * You should have received a copy of the GNU Lesser General Public License * + * along with this program; if not, write to the Free Software Foundation, * + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. * + ****************************************************************************/ + +#ifndef SVGFILL_H +#define SVGFILL_H + +#ifdef IFC_SHARED_BUILD +#ifdef _WIN32 +#ifdef svgfill_EXPORTS +#define SVGFILL_API __declspec(dllexport) +#else +#define SVGFILL_API __declspec(dllimport) +#endif +#else // simply assume *nix + GCC-like compiler +#define SVGFILL_API __attribute__((visibility("default"))) +#endif +#else +#define SVGFILL_API +#endif + +#include + +#include +#include + +namespace svgfill { + typedef std::array point_2; + typedef std::array line_segment_2; + typedef std::vector loop_2; + struct SVGFILL_API polygon_2 { + loop_2 boundary; + std::vector inner_boundaries; + point_2 point_inside; + }; + + enum solver { + CARTESIAN_DOUBLE, + CARTESIAN_QUOTIENT, + FILTERED_CARTESIAN_QUOTIENT, + EXACT_PREDICATES, + EXACT_CONSTRUCTIONS + }; + + class SVGFILL_API abstract_arrangement { + public: + virtual ~abstract_arrangement() {} + virtual bool operator()(double eps, const std::vector& segments, std::function& progress) = 0; + virtual bool write(std::vector& polygons, std::function& progress) = 0; + virtual void merge(const std::vector& edge_indices) = 0; + virtual std::vector get_face_pairs() = 0; + virtual size_t num_edges() = 0; + virtual size_t num_faces() = 0; + }; + + class SVGFILL_API context { + private: + solver solver_; + double eps_; + std::vector segments_; + std::function progress_; + // std::vector polygons_; + abstract_arrangement* arr_; + + public: + context(solver s, double eps) + : solver_(s) + , eps_(eps) + , arr_(nullptr) + {} + + context(solver s, double eps, std::function& progress) + : solver_(s) + , eps_(eps) + , progress_(progress) + , arr_(nullptr) + {} + + void add(const std::vector& segments); + bool build(); + std::vector get_face_pairs() { + return arr_->get_face_pairs(); + } + void merge(const std::vector& edge_indices); + void write(std::vector>&); + size_t num_edges() { return arr_->num_edges(); } + size_t num_faces() { return arr_->num_faces(); } + + ~context() { + delete arr_; + } + }; + + SVGFILL_API bool svg_to_line_segments(const std::string& data, const boost::optional& class_name, std::vector>& segments); + SVGFILL_API bool line_segments_to_polygons(solver s, double eps, const std::vector>& segments, std::vector>& polygons); + SVGFILL_API bool line_segments_to_polygons(solver s, double eps, const std::vector>& segments, std::vector>& polygons, std::function& progress); + SVGFILL_API std::string polygons_to_svg(const std::vector>& polygons, bool random_color=false); + SVGFILL_API std::string polygons_to_svg(const std::vector& polygons, bool random_color = false); + SVGFILL_API bool svg_to_polygons(const std::string& data, const boost::optional& class_name, std::vector& polygons); + SVGFILL_API bool arrange_polygons(const std::vector& polygons, std::vector& arranged); +} + +#endif