diff --git a/cartocrow/CMakeLists.txt b/cartocrow/CMakeLists.txt index 8d8b5b5e..686ae4a6 100644 --- a/cartocrow/CMakeLists.txt +++ b/cartocrow/CMakeLists.txt @@ -16,6 +16,7 @@ add_subdirectory(necklace_map) add_subdirectory(simplification) add_subdirectory(flow_map) add_subdirectory(isoline_simplification) +add_subdirectory(simplesets) add_library(cartocrow_lib INTERFACE) target_link_libraries( diff --git a/cartocrow/core/core.cpp b/cartocrow/core/core.cpp index 375e25bd..9361e5fe 100644 --- a/cartocrow/core/core.cpp +++ b/cartocrow/core/core.cpp @@ -24,6 +24,10 @@ Created by tvl (t.vanlankveld@esciencecenter.nl) on 05-12-2019 namespace cartocrow { +Color::Color() : r(0), g(0), b(0) {} +Color::Color(int r, int g, int b) : r(r), g(g), b(b) {} +Color::Color(int rgb) : r((rgb & 0xff0000) >> 16), g((rgb & 0x00ff00) >> 8), b(rgb & 0x0000ff) {} + Number wrapAngle(Number alpha, Number beta) { return wrap(alpha, beta, beta + M_2xPI); } diff --git a/cartocrow/core/core.h b/cartocrow/core/core.h index 377f3d93..f725c35f 100644 --- a/cartocrow/core/core.h +++ b/cartocrow/core/core.h @@ -28,6 +28,7 @@ Created by tvl (t.vanlankveld@esciencecenter.nl) on 07-11-2019 #include #include #include +#include #include #include #include @@ -60,6 +61,8 @@ template using Line = CGAL::Line_2; template using Segment = CGAL::Segment_2; /// A ray in the plane. See \ref CGAL::Ray_2. template using Ray = CGAL::Ray_2; +/// An axis-aligned rectangle in the plane. See \ref CGAL::Iso_rectangle_2. +template using Rectangle = CGAL::Iso_rectangle_2; /// A polygon in the plane. See \ref CGAL::Polygon_2. template using Polygon = CGAL::Polygon_2; @@ -170,6 +173,12 @@ struct Color { int g; /// Blue component (integer 0-255). int b; + /// Constructs the color black. + Color(); + /// Constructs a color. + Color(int r, int g, int b); + /// Constructs a color from a single integer (useful combined with hexadecimal literals, e.g. 0xFFFFFF). + Color(int rgb); }; /// Wraps the given number \f$n\f$ to the interval \f$[a, b)\f$. diff --git a/cartocrow/isoline_simplification/symmetric_difference.cpp b/cartocrow/isoline_simplification/symmetric_difference.cpp index 69fda924..fb5b096f 100644 --- a/cartocrow/isoline_simplification/symmetric_difference.cpp +++ b/cartocrow/isoline_simplification/symmetric_difference.cpp @@ -30,7 +30,7 @@ enum Side { TOP = 3, }; -Side closest_side(const Point& point, const Exact::Iso_rectangle_2& bb) { +Side closest_side(const Point& point, const Rectangle& bb) { std::vector dist({ to_double(point.x()) - to_double(bb.xmin()), to_double(point.y()) - to_double(bb.ymin()), to_double(bb.xmax()) - to_double(point.x()), to_double(bb.ymax()) - to_double(point.y()) }); auto it = std::min_element(dist.begin(), dist.end()); @@ -50,7 +50,7 @@ Vector side_direction(const Side& side) { } } -Point proj_on_side(Point p, Side side, const Exact::Iso_rectangle_2& rect) { +Point proj_on_side(Point p, Side side, const Rectangle& rect) { switch (side) { case LEFT: return { rect.xmin(), p.y() }; @@ -63,7 +63,7 @@ Point proj_on_side(Point p, Side side, const Exact::Iso_rectangle_ } } -Point corner(const Side& side1, const Side& side2, const Exact::Iso_rectangle_2& rect) { +Point corner(const Side& side1, const Side& side2, const Rectangle& rect) { if (side1 > side2) return corner(side2, side1, rect); auto dist = abs(side1 - side2); if (dist == 1) { @@ -77,7 +77,7 @@ Side next_side(const Side& side) { return static_cast((side + 1) % 4); } -Polygon close_isoline(const Isoline& isoline, Exact::Iso_rectangle_2& bb, Side source_side, Side target_side) { +Polygon close_isoline(const Isoline& isoline, Rectangle& bb, Side source_side, Side target_side) { std::vector> points; std::transform(isoline.m_points.begin(), isoline.m_points.end(), std::back_inserter(points), [](Point p) { return Point(p.x(), p.y()); }); @@ -133,7 +133,7 @@ double symmetric_difference(const Isoline& original, const Isoline& simpli std::transform(simplified.m_points.begin(), simplified.m_points.end(), std::back_inserter(all_points), [](Point p) { return Point(p.x(), p.y()); }); - Exact::Iso_rectangle_2 bb = CGAL::bounding_box(all_points.begin(), all_points.end()); + Rectangle bb = CGAL::bounding_box(all_points.begin(), all_points.end()); auto source_side = closest_side( Point(original.m_points.front().x(), original.m_points.front().y()), bb); diff --git a/cartocrow/renderer/CMakeLists.txt b/cartocrow/renderer/CMakeLists.txt index 5ea84fb3..abb324c1 100644 --- a/cartocrow/renderer/CMakeLists.txt +++ b/cartocrow/renderer/CMakeLists.txt @@ -3,6 +3,7 @@ set(SOURCES geometry_widget.cpp ipe_renderer.cpp painting_renderer.cpp + function_painting.cpp render_path.cpp svg_renderer.cpp ) @@ -12,6 +13,7 @@ set(HEADERS geometry_widget.h ipe_renderer.h painting_renderer.h + function_painting.h render_path.h svg_renderer.h ) diff --git a/cartocrow/renderer/function_painting.cpp b/cartocrow/renderer/function_painting.cpp new file mode 100644 index 00000000..e5259ab0 --- /dev/null +++ b/cartocrow/renderer/function_painting.cpp @@ -0,0 +1,10 @@ +#include "function_painting.h" + +namespace cartocrow::renderer { +FunctionPainting::FunctionPainting(const std::function& draw_function) + : m_draw_function(draw_function) {} + +void FunctionPainting::paint(GeometryRenderer& renderer) const { + m_draw_function(renderer); +} +} \ No newline at end of file diff --git a/cartocrow/renderer/function_painting.h b/cartocrow/renderer/function_painting.h new file mode 100644 index 00000000..681593cf --- /dev/null +++ b/cartocrow/renderer/function_painting.h @@ -0,0 +1,18 @@ +#ifndef CARTOCROW_FUNCTION_PAINTING_H +#define CARTOCROW_FUNCTION_PAINTING_H + +#include "geometry_renderer.h" +#include "geometry_painting.h" + +namespace cartocrow::renderer { +class FunctionPainting : public GeometryPainting { + public: + FunctionPainting(const std::function& draw_function); + void paint(renderer::GeometryRenderer& renderer) const override; + + private: + const std::function m_draw_function; +}; +} + +#endif //CARTOCROW_FUNCTION_PAINTING_H diff --git a/cartocrow/renderer/geometry_widget.cpp b/cartocrow/renderer/geometry_widget.cpp index da642ebe..2e7f70c5 100644 --- a/cartocrow/renderer/geometry_widget.cpp +++ b/cartocrow/renderer/geometry_widget.cpp @@ -23,6 +23,8 @@ along with this program. If not, see . #include "ipe_renderer.h" #include "svg_renderer.h" +#include "function_painting.h" + #include #include #include @@ -546,7 +548,7 @@ void GeometryWidget::draw(const BezierSpline& s) { void GeometryWidget::draw(const Ray& r) { Box bounds = inverseConvertBox(rect()); - auto result = intersection(r, CGAL::Iso_rectangle_2(Point(bounds.xmin(), bounds.ymin()), Point(bounds.xmax(), bounds.ymax()))); + auto result = intersection(r, Rectangle(Point(bounds.xmin(), bounds.ymin()), Point(bounds.xmax(), bounds.ymax()))); if (result) { if (const Segment* s = boost::get>(&*result)) { int oldMode = m_style.m_mode; @@ -685,6 +687,11 @@ void GeometryWidget::addPainting(std::shared_ptr painting, con updateLayerList(); } +void GeometryWidget::addPainting(const std::function& draw_function, const std::string& name) { + auto painting = std::make_shared(draw_function); + addPainting(painting, name); +} + void GeometryWidget::clear() { m_paintings.clear(); updateLayerList(); diff --git a/cartocrow/renderer/geometry_widget.h b/cartocrow/renderer/geometry_widget.h index beba95f6..67d81b16 100644 --- a/cartocrow/renderer/geometry_widget.h +++ b/cartocrow/renderer/geometry_widget.h @@ -180,6 +180,7 @@ class GeometryWidget : public QWidget, public GeometryRenderer { /// Adds a new painting to this widget. void addPainting(std::shared_ptr painting, const std::string& name); + void addPainting(const std::function& draw_function, const std::string& name); /// Removes all paintings from this widget. void clear(); diff --git a/cartocrow/renderer/ipe_renderer.cpp b/cartocrow/renderer/ipe_renderer.cpp index c4f3ee3e..e289e5eb 100644 --- a/cartocrow/renderer/ipe_renderer.cpp +++ b/cartocrow/renderer/ipe_renderer.cpp @@ -19,6 +19,7 @@ along with this program. If not, see . #include "ipe_renderer.h" +#include "function_painting.h" #include "geometry_renderer.h" #include @@ -74,9 +75,18 @@ void IpeRenderer::save(const std::filesystem::path& file) { m_alphaSheet->setName("alpha-values"); document.cascade()->insert(2, m_alphaSheet); setFillOpacity(255); // add default alpha to style sheet + setStrokeOpacity(255); // add default alpha to style sheet m_page = new ipe::Page(); - for (auto painting : m_paintings) { + document.push_back(m_page); + + int current_page = 0; + + for (auto painting : m_paintings) { // Assumes m_paintings are ordered in increasing page_index + while (painting.page_index > current_page) { + m_page = new ipe::Page(); + document.push_back(m_page); + } pushStyle(); if (auto name = painting.name) { m_page->addLayer(name->c_str()); @@ -88,7 +98,6 @@ void IpeRenderer::save(const std::filesystem::path& file) { popStyle(); } - document.push_back(m_page); document.save(file.string().c_str(), ipe::FileFormat::Xml, 0); } @@ -320,17 +329,37 @@ ipe::AllAttributes IpeRenderer::getAttributesForStyle() const { attributes.iStroke = ipe::Attribute(ipe::Color(m_style.m_strokeColor)); attributes.iFill = ipe::Attribute(ipe::Color(m_style.m_fillColor)); attributes.iOpacity = m_style.m_fillOpacity; + attributes.iStrokeOpacity = m_style.m_strokeOpacity; return attributes; } + +void IpeRenderer::addPainting(const std::function& draw_function) { + auto painting = std::make_shared(draw_function); + addPainting(painting); +} + +void IpeRenderer::addPainting(const std::function& draw_function, const std::string& name) { + auto painting = std::make_shared(draw_function); + addPainting(painting, name); +} + void IpeRenderer::addPainting(const std::shared_ptr& painting) { - m_paintings.push_back(DrawnPainting{painting}); + m_paintings.push_back(DrawnPainting{painting, std::nullopt, m_page_index}); } void IpeRenderer::addPainting(const std::shared_ptr& painting, const std::string& name) { std::string spaceless; std::replace_copy(name.begin(), name.end(), std::back_inserter(spaceless), ' ', '_'); - m_paintings.push_back(DrawnPainting{painting, spaceless}); + m_paintings.push_back(DrawnPainting{painting, spaceless, m_page_index}); +} + +void IpeRenderer::nextPage() { + ++m_page_index; +} + +int IpeRenderer::currentPage() { + return m_page_index; } std::string IpeRenderer::escapeForLaTeX(const std::string& text) const { diff --git a/cartocrow/renderer/ipe_renderer.h b/cartocrow/renderer/ipe_renderer.h index ecaf7bd6..916d1b58 100644 --- a/cartocrow/renderer/ipe_renderer.h +++ b/cartocrow/renderer/ipe_renderer.h @@ -98,9 +98,16 @@ class IpeRenderer : public GeometryRenderer { void setFill(Color color) override; void setFillOpacity(int alpha) override; + void addPainting(const std::function& draw_function); + void addPainting(const std::function& draw_function, const std::string& name); void addPainting(const std::shared_ptr& painting); void addPainting(const std::shared_ptr& painting, const std::string& name); + /// Paintings will be added to a new page. + void nextPage(); + /// Returns the current page index. + int currentPage(); + private: /// Converts a polygon to an Ipe curve. ipe::Curve* convertPolygonToCurve(const Polygon& p) const; @@ -119,6 +126,8 @@ class IpeRenderer : public GeometryRenderer { std::shared_ptr m_painting; /// The name of the painting displayed as a layer name in ipe. std::optional name; + /// The Ipe page the painting will be drawn onto. + int page_index; }; /// The paintings we're drawing. @@ -138,6 +147,8 @@ class IpeRenderer : public GeometryRenderer { ipe::StyleSheet* m_alphaSheet; /// The index of the Ipe layer we are currently drawing to. int m_layer; + /// The index of the Ipe page a painting will get drawn to. + int m_page_index = 0; }; } // namespace cartocrow::renderer diff --git a/cartocrow/simplesets/CMakeLists.txt b/cartocrow/simplesets/CMakeLists.txt new file mode 100644 index 00000000..19111279 --- /dev/null +++ b/cartocrow/simplesets/CMakeLists.txt @@ -0,0 +1,55 @@ +set(SOURCES + cat_point.cpp + types.cpp + parse_input.cpp + patterns/pattern.cpp + patterns/single_point.cpp + patterns/matching.cpp + patterns/island.cpp + patterns/bank.cpp + dilated/dilated_poly.cpp + helpers/approximate_convex_hull.cpp + helpers/cs_curve_helpers.cpp + helpers/cs_polygon_helpers.cpp + helpers/cs_polyline_helpers.cpp + helpers/poly_line_gon_intersection.cpp + partition_algorithm.cpp + partition_painting.cpp + drawing_algorithm.cpp + grow_circles.cpp +) +set(HEADERS + types.h + settings.h + cat_point.h + parse_input.h + patterns/pattern.h + patterns/poly_pattern.h + patterns/single_point.h + patterns/matching.h + patterns/island.h + patterns/bank.h + dilated/dilated_poly.h + helpers/approximate_convex_hull.h + helpers/arrangement_helpers.h + helpers/point_voronoi_helpers.h + helpers/poly_line_gon_intersection.h + helpers/cropped_voronoi.h + helpers/cs_curve_helpers.h + helpers/cs_polygon_helpers.h + helpers/cs_polyline_helpers.h + partition_algorithm.h + partition_painting.h + partition.h + drawing_algorithm.h + general_polyline.h + grow_circles.h +) + +add_library(simplesets ${SOURCES}) +target_link_libraries(simplesets + PUBLIC core +) + +cartocrow_install_module(simplesets) +install(FILES ${HEADERS} DESTINATION ${CARTOCROW_INSTALL_DIR}/simplesets) diff --git a/cartocrow/simplesets/cat_point.cpp b/cartocrow/simplesets/cat_point.cpp new file mode 100644 index 00000000..511cfa70 --- /dev/null +++ b/cartocrow/simplesets/cat_point.cpp @@ -0,0 +1,7 @@ +#include "cat_point.h" + +namespace cartocrow::simplesets { +std::ostream& operator<<(std::ostream& os, CatPoint const& catPoint) { + return os << catPoint.category << " " << catPoint.point; +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/cat_point.h b/cartocrow/simplesets/cat_point.h new file mode 100644 index 00000000..7776e5ec --- /dev/null +++ b/cartocrow/simplesets/cat_point.h @@ -0,0 +1,18 @@ +#ifndef CARTOCROW_CAT_POINT_H +#define CARTOCROW_CAT_POINT_H + +#include "types.h" + +namespace cartocrow::simplesets { +/// Categorical point +struct CatPoint { + unsigned int category; + Point point; + CatPoint(unsigned int category, Point point) : category(category), point(point) {}; + bool operator==(const CatPoint&) const = default; +}; + +std::ostream& operator<<(std::ostream& os, CatPoint const& catPoint); +} + +#endif //CARTOCROW_CAT_POINT_H diff --git a/cartocrow/simplesets/dilated/dilated_poly.cpp b/cartocrow/simplesets/dilated/dilated_poly.cpp new file mode 100644 index 00000000..28155ba1 --- /dev/null +++ b/cartocrow/simplesets/dilated/dilated_poly.cpp @@ -0,0 +1,73 @@ +#include "dilated_poly.h" +#include "../helpers/cs_polygon_helpers.h" +#include "../helpers/arrangement_helpers.h" +#include + +namespace cartocrow::simplesets { +CSPolygon dilateSegment(const Segment& segment, const Number& dilationRadius) { + std::vector> points({makeExact(segment.source()), makeExact(segment.target())}); + Polygon polygon(points.begin(), points.end()); + auto dilation = CGAL::approximated_offset_2(polygon, dilationRadius, M_EPSILON); + return dilation.outer_boundary(); +} + +Dilated::Dilated(const PolyPattern& polyPattern, const Number& dilationRadius) { + m_catPoints = polyPattern.catPoints(); + + auto cont = polyPattern.poly(); + + if (holds_alternative>(cont)) { + auto exactPolygon = makeExact(std::get>(cont)); + + if (exactPolygon.size() == 1) { + CSTraits::Rational_circle_2 circle(exactPolygon.vertex(0), dilationRadius * dilationRadius); + m_contour = circleToCSPolygon(circle); + if (m_contour.orientation() == CGAL::CLOCKWISE) { + m_contour.reverse_orientation(); + } + return; + } + + auto dilation = CGAL::approximated_offset_2(exactPolygon, dilationRadius, M_EPSILON); + if (dilation.has_holes()) { + throw std::runtime_error("Did not expect holes after dilating a polygonal pattern."); + } + m_contour = dilation.outer_boundary(); + if (m_contour.orientation() == CGAL::CLOCKWISE) { + m_contour.reverse_orientation(); + } + } else if (holds_alternative>(cont)) { + // 1. Dilate each segment + // 2. Make arrangement of dilated segments + // 3. Traverse and extract outer boundary + // (Assume that the dilation result has no holes. + // To ensure this we need to constrain the (relative) point size.) + + CSArrangement arr; + + auto polyline = std::get>(cont); + for (auto eit = polyline.edges_begin(); eit != polyline.edges_end(); ++eit) { + auto seg = *eit; + auto dilated = dilateSegment(seg, dilationRadius); + for (auto cit = dilated.curves_begin(); cit != dilated.curves_end(); ++cit) { + CGAL::insert(arr, *cit); + } + } + + m_contour = ccb_to_polygon(*(arr.unbounded_face()->inner_ccbs_begin())); + if (m_contour.orientation() == CGAL::CLOCKWISE) { + m_contour.reverse_orientation(); + } + } else { + throw std::runtime_error("Unknown pattern poly."); + } +} + +const std::vector& Dilated::catPoints() const { + return m_catPoints; +} + +std::variant, Polygon, CSPolygon> Dilated::contour() const { + return m_contour; +} +} diff --git a/cartocrow/simplesets/dilated/dilated_poly.h b/cartocrow/simplesets/dilated/dilated_poly.h new file mode 100644 index 00000000..66d14ca7 --- /dev/null +++ b/cartocrow/simplesets/dilated/dilated_poly.h @@ -0,0 +1,18 @@ +#ifndef CARTOCROW_DILATED_POLY_H +#define CARTOCROW_DILATED_POLY_H + +#include "../patterns/poly_pattern.h" + +namespace cartocrow::simplesets { +class Dilated : public Pattern { + public: + Dilated(const PolyPattern& polyPattern, const Number& dilationRadius); + std::variant, Polygon, CSPolygon> contour() const override; + const std::vector& catPoints() const override; + CSPolygon m_contour; + + private: + std::vector m_catPoints; +}; +} +#endif //CARTOCROW_DILATED_POLY_H diff --git a/cartocrow/simplesets/drawing_algorithm.cpp b/cartocrow/simplesets/drawing_algorithm.cpp new file mode 100644 index 00000000..3fd6cedc --- /dev/null +++ b/cartocrow/simplesets/drawing_algorithm.cpp @@ -0,0 +1,1153 @@ +#include "drawing_algorithm.h" + +#include "cartocrow/simplesets/helpers/poly_line_gon_intersection.h" +#include "grow_circles.h" +#include "helpers/approximate_convex_hull.h" +#include "helpers/arrangement_helpers.h" +#include "helpers/cs_curve_helpers.h" +#include "helpers/cs_polygon_helpers.h" +#include "helpers/cs_polyline_helpers.h" +#include +#include +#include +#include "cartocrow/renderer/ipe_renderer.h" + +using namespace cartocrow::renderer; + +namespace cartocrow::simplesets { +// todo: deal with holes +CSPolygon face_to_polygon(const Face& face) { + assert(face.number_of_holes() == 0); + return ccb_to_polygon(face.outer_ccb()); +} + +X_monotone_curve_2 make_x_monotone(const Segment& segment) { + CSTraits traits; + auto make_x_monotone = traits.make_x_monotone_2_object(); + std::vector> curves_and_points; + make_x_monotone(segment, std::back_inserter(curves_and_points)); + std::vector curves; + + // There should not be any isolated points + for (auto kinda_curve : curves_and_points) { + assert(kinda_curve.which() == 1); + auto curve = boost::get(kinda_curve); + curves.push_back(curve); + } + + assert(curves.size() == 1); + return curves[0]; +} + +Point get_approx_point_on_boundary(const Face& face) { + auto curve = face.outer_ccb()->curve(); + Rectangle bbox = curve.bbox(); + Rectangle rect({bbox.xmin() - 1, bbox.ymin() - 1}, {bbox.xmax() + 1, bbox.ymax() + 1}); + Point approx_source = makeExact(approximateAlgebraic(curve.source())); + Point approx_target = makeExact(approximateAlgebraic(curve.target())); + Point middle = CGAL::midpoint(approx_source, approx_target); + if (curve.is_linear()) { + return middle; + } else { + assert(curve.is_circular()); + Circle circle = curve.supporting_circle(); + Line l(approx_source, approx_target); + Line pl = l.perpendicular(middle); + auto inter = CGAL::intersection(pl, rect); + assert(inter.has_value()); + assert(inter->which() == 1); + Segment seg = boost::get>(*inter); +// CGAL::intersection(pl, circle); + CSTraits traits; + auto make_x_monotone = traits.make_x_monotone_2_object(); + std::vector> curves_and_points; +// make_x_monotone(circle, std::back_inserter(curves_and_points)); + make_x_monotone(seg, std::back_inserter(curves_and_points)); + std::vector curves; + + // There should not be any isolated points + for (auto kinda_curve : curves_and_points) { + auto curve = boost::get(kinda_curve); + curves.push_back(curve); + } + + typedef std::pair Intersection_point; + typedef boost::variant Intersection_result; + std::vector intersection_results; + assert(curves.size() == 1); + curve.intersect(curves[0], std::back_inserter(intersection_results)); + + std::vector intersection_pts; + for (const auto& result : intersection_results) { + assert(result.which() == 0); + intersection_pts.push_back(boost::get(result).first); + } + + assert(intersection_pts.size() == 1); + return makeExact(approximateAlgebraic(intersection_pts[0])); + } +} + +Point get_point_in(const Face& face) { + auto poly = face_to_polygon(face); + Rectangle bbox = poly.bbox(); + Rectangle rect({bbox.xmin() - 1, bbox.ymin() - 1}, {bbox.xmax() + 1, bbox.ymax() + 1}); + Point point_outside(rect.xmin(), rect.ymin()); + Point approx_point_on_boundary = get_approx_point_on_boundary(face); + Line line(point_outside, approx_point_on_boundary); + auto line_inter_box = CGAL::intersection(rect, line); + auto seg = boost::get>(*line_inter_box); + auto seg_curve = make_x_monotone(seg); + std::vector intersection_pts; + for (auto cit = poly.curves_begin(); cit != poly.curves_end(); cit++) { + auto curve = *cit; + typedef std::pair Intersection_point; + typedef boost::variant Intersection_result; + std::vector intersection_results; + curve.intersect(seg_curve, std::back_inserter(intersection_results)); + + for (const auto& result : intersection_results) { + assert(result.which() == 0); + intersection_pts.push_back(boost::get(result).first); + } + } + + std::sort(intersection_pts.begin(), intersection_pts.end(), [&seg](const auto& pt1, const auto& pt2) { + Vector v = seg.supporting_line().to_vector(); + return v * (makeExact(approximateAlgebraic(pt1)) - makeExact(approximateAlgebraic(pt2))) < 0; + }); + intersection_pts.erase(std::unique(intersection_pts.begin(), intersection_pts.end()), intersection_pts.end()); + + Point approx_source; + Point approx_target; + if (intersection_pts.size() >= 2) { + approx_source = makeExact(approximateAlgebraic(intersection_pts[0])); + approx_target = makeExact(approximateAlgebraic(intersection_pts[1])); + } else { + approx_source = makeExact(approximateAlgebraic(intersection_pts[0])); + approx_target = approx_point_on_boundary; + } + + return midpoint(Segment(approx_source, approx_target)); +} + +std::vector +connectedComponents(const DilatedPatternArrangement& arr, const std::function& in_component) { + std::vector remaining; + for (auto fit = arr.faces_begin(); fit != arr.faces_end(); ++fit) { + auto fh = fit.ptr(); + if (in_component(fh)) { + remaining.emplace_back(fh); + } + } + + std::vector components; + + while (!remaining.empty()) { + // We do a BFS + std::vector compFaces; + std::vector compBoundaryEdges; + auto first = remaining.front(); + std::deque q; + q.push_back(first); + + while (!q.empty()) { + auto f = q.front(); + q.pop_front(); + compFaces.push_back(f); + + // Go through boundaries of this face + std::vector ccbs; + std::copy(f->outer_ccbs_begin(), f->outer_ccbs_end(), std::back_inserter(ccbs)); + std::copy(f->inner_ccbs_begin(), f->inner_ccbs_end(), std::back_inserter(ccbs)); + for (auto ccb_start : ccbs) { + auto ccb_it = ccb_start; + + // Go through each neighbouring face + do { + auto candidate = ccb_it->twin()->face(); + if (!in_component(candidate)) { + compBoundaryEdges.emplace_back(ccb_it.ptr()); + } else { + // If this is one of the provided faces, and not yet added to queue or compFaces, add it to queue. + if (std::find(compFaces.begin(), compFaces.end(), candidate) == + compFaces.end() && + std::find(q.begin(), q.end(), candidate) == q.end()) { + q.push_back(candidate); + } + } + } while (++ccb_it != ccb_start); + } + } + + // Done with this connected component + remaining.erase(std::remove_if(remaining.begin(), remaining.end(), [&compFaces](const auto& f) { + return std::find(compFaces.begin(), compFaces.end(), f) != compFaces.end(); + }), remaining.end()); + components.emplace_back(std::move(compFaces), std::move(compBoundaryEdges), in_component); + } + + return components; +} + +Component::Component(std::vector faces, std::vector boundary_edges, std::function in_component) : + m_faces(std::move(faces)), m_in_component(std::move(in_component)) { + while (!boundary_edges.empty()) { + auto he = boundary_edges.front(); + auto circ_start = ComponentCcbCirculator(he, m_in_component); + auto circ = circ_start; + + std::vector xm_curves; + auto last_it = boundary_edges.end(); + do { + last_it = std::remove(boundary_edges.begin(), last_it, circ.handle()); + xm_curves.push_back(circ->curve()); + } while (++circ != circ_start); + boundary_edges.erase(last_it, boundary_edges.end()); + + CSPolygon polygon(xm_curves.begin(), xm_curves.end()); + auto orientation = polygon.orientation(); + if (orientation == CGAL::COUNTERCLOCKWISE) { + m_outer_ccbs.push_back(circ_start); + } else if (orientation == CGAL::CLOCKWISE) { + m_inner_ccbs.push_back(circ_start); + } else { + throw std::runtime_error("Face orientation is not clockwise nor counterclockwise."); + } + } +} + +std::ostream& operator<<(std::ostream& out, const Order& o) { + return out << to_string(o); +} + +std::ostream& operator<<(std::ostream& out, const Relation& r) { + return out << r.left << " " << r.ordering + << (r.preference != r.ordering ? "(" + to_string(r.preference) + ")" : "") + << " " << r.right; +} + +bool operator==(const Relation& lhs, const Relation& rhs) { + return lhs.left == rhs.left && lhs.right == rhs.right && lhs.preference == rhs.preference && lhs.ordering == rhs.ordering; +} + +bool operator==(const Hyperedge& lhs, const Hyperedge& rhs) { + return lhs.relations == rhs.relations && lhs.origins == rhs.origins; +} + +DilatedPatternDrawing::DilatedPatternDrawing(const Partition& partition, const GeneralSettings& gs, const ComputeDrawingSettings& cds) + : m_gs(gs), m_cds(cds) { + for (const auto& p : partition) { + m_dilated.emplace_back(*p, gs.dilationRadius()); + } + + std::vector curves; + std::vector> curves_data; + for (int i = 0; i < m_dilated.size(); i++) { + for (auto cit = m_dilated[i].m_contour.curves_begin(); cit != m_dilated[i].m_contour.curves_end(); ++cit) { + auto x_monotone = *cit; + CSTraits::Curve_2 curve = toCurve(x_monotone); + curves.emplace_back(curve); + curves_data.emplace_back(curve, i); + } + } + + CGAL::insert(m_arr, curves.begin(), curves.end()); + // Set for each face which pattern it is a subset of. + for (auto fit = m_arr.faces_begin(); fit != m_arr.faces_end(); ++fit) { + if (fit->is_unbounded()) continue; + auto pt = get_point_in(*fit); + std::vector origins; + for (int i = 0; i < m_dilated.size(); i++) { + if (on_or_inside(m_dilated[i].m_contour, pt)) { + origins.push_back(i); + m_iToFaces[i].push_back(fit); + } + } + fit->set_data(FaceData{origins}); + } + + // Store in each half-edges form which pattern it originates. + for (auto cit = m_arr.curves_begin(); cit != m_arr.curves_end(); ++cit) { + CSTraits::Curve_2 curve = *cit; + auto curve_data = std::find_if(curves_data.begin(), curves_data.end(), [&curve](const auto& pair) { + auto other = pair.first; + return curve.source() == other.source() && curve.target() == other.target() && + (curve.is_linear() && other.is_linear() || + curve.is_circular() && other.is_circular() && curve.supporting_circle() == other.supporting_circle()); + }); + for (auto eit = m_arr.induced_edges_begin(cit); eit != m_arr.induced_edges_end(cit); ++eit) { + DilatedPatternArrangement::Halfedge_handle eh = *eit; + eh->data().origins.push_back(curve_data->second); + eh->twin()->data().origins.push_back(curve_data->second); + } + } + + // Sort for std::set_intersection operation + for (int i = 0; i < m_dilated.size(); i++) { + auto& facesI = m_iToFaces[i]; + std::sort(facesI.begin(), facesI.end()); + } + + for (int i = 0; i < m_dilated.size(); i++) { + for (int j = i + 1; j < m_dilated.size(); j++) { + auto cs = intersectionComponents(i, j); + for (auto& c : cs) { + auto rel = computePreference(i, j, c); + for (auto fit = c.faces_begin(); fit != c.faces_end(); ++fit) { + fit->data().relations.push_back(rel); + } + } + } + } + + auto hEdges = hyperedges(); + for (auto edge : hEdges) { + auto order = getRelationOrder(*edge); + if (!order.has_value()) { + std::cerr << "Hyperedge has cycle; ignoring preferences in this hyperedge." << std::endl; + for (const auto& r : edge->relations) { + r->ordering = Order::EQUAL; + } + order = getRelationOrder(*edge); + assert(order.has_value()); + } + setRelationOrder(*edge, *order); + } + + for (auto fit = m_arr.faces_begin(); fit != m_arr.faces_end(); ++fit) { + auto& data = fit->data(); + if (data.origins.empty()) continue; + auto ordering = computeTotalOrder(data.origins, data.relations); + if (!ordering.has_value()) { + throw std::runtime_error("Impossible: no total order in a face"); + } + data.ordering = *ordering; + } + + for (int i = 0; i < m_dilated.size(); ++i) { + auto cs = intersectionComponents(i); + for (auto& c : cs) { + std::unordered_set avoidees; + for (auto fit = c.faces_begin(); fit != c.faces_end(); ++fit) { + for (int j : fit->data().ordering) { + if (j == i) break; + avoidees.insert(j); + } + } + if (avoidees.empty()) + continue; + + auto bpis = boundaryParts(c, i); + auto disks = includeExcludeDisks(i, avoidees, c); + auto inclDisks = disks.include; + auto exclDisks = disks.exclude; + + if (exclDisks.empty()) { + continue; + } + auto componentPolygon = ccb_to_polygon(c.outer_ccb()); + auto morphedComponentPolygon = morph(bpis, componentPolygon, inclDisks, exclDisks, m_gs, m_cds); + + // Compute the morphed version of the CSPolygon for this component. + // Set, for each face in component c, the morphed face to the intersection of this CSPolygon with the face. + // Set the morphed edges to the intersection of the boundary of the polygon with the face. + for (auto fit = c.faces_begin(); fit != c.faces_end(); ++fit) { + auto facePolygon = face_to_polygon(*fit); + + for (const auto& bp : bpis) { + CSPolyline mb = associatedBoundary(componentPolygon, morphedComponentPolygon, bp); + intersection(mb, facePolygon, std::back_inserter(fit->data().morphedEdges[i]), false, true); + } + + std::vector morphedFacePolygonsWithHoles; + + CGAL::intersection(morphedComponentPolygon, facePolygon, std::back_inserter(morphedFacePolygonsWithHoles)); + auto& mf = fit->data().morphedFace[i]; + for (const auto& morphedFacePolygonWithHoles : morphedFacePolygonsWithHoles) { + assert(!morphedFacePolygonWithHoles.has_holes()); // todo + mf.push_back(morphedFacePolygonWithHoles.outer_boundary()); + } + } + } + } +} + +bool overlap(const Circle& c1, const Circle& c2) { + return sqrt(CGAL::squared_distance(c1.center(), c2.center())) <= sqrt(c1.squared_radius()) + sqrt(c2.squared_radius()); +} + +std::vector>>> connectedDisks(const std::vector>& disks) { + typedef std::vector>> Component; + std::vector components; + + auto merge = [&components](std::vector& comps) { + Component newComponent; + for (auto& comp : comps) { + std::copy(comp->begin(), comp->end(), std::back_inserter(newComponent)); + comp->clear(); + } + components.erase(std::remove_if(components.begin(), components.end(), [](const Component& cs) { + return cs.empty(); + }), components.end()); + return newComponent; + }; + + for (int i = 0; i < disks.size(); ++i) { + const auto& d = disks[i]; + std::vector overlaps; + for (auto& comp : components) { + bool anyOverlap = false; + for (const auto& other : comp) { + if (overlap(d, other.second)) { + anyOverlap = true; + break; + } + } + if (anyOverlap) { + overlaps.push_back(&comp); + } + } + auto merged = merge(overlaps); + merged.emplace_back(i, d); + components.push_back(merged); + } + + return components; +} + +CSPolygon thinRectangle(const Point& p, const OneRootPoint& n, const Number& w) { + Point nApproxInexact = approximateAlgebraic(n); + Point nApprox(nApproxInexact.x(), nApproxInexact.y()); + Vector d = nApprox - p; + auto dl = sqrt(CGAL::to_double(d.squared_length())); + Vector normalized = d / dl; + auto perp = normalized.perpendicular(CGAL::COUNTERCLOCKWISE) * w / 2; + + Point p1 = p - perp; + Point p2 = nApprox + normalized * w / 10 - perp; + Point p3 = nApprox + normalized * w / 10 + perp; + Point p4 = p + perp; + + std::vector curves({ + Curve_2(p1, p2), + Curve_2(p2, p3), + Curve_2(p3, p4), + Curve_2(p4, p1) + }); + + std::vector xm_curves; + curvesToXMonotoneCurves(curves.begin(), curves.end(), std::back_inserter(xm_curves)); + + return {xm_curves.begin(), xm_curves.end()}; +} + +CSPolygon morph(const std::vector& boundaryParts, const CSPolygon& componentShape, const std::vector>& inclDisks, + const std::vector>& exclDisks, const GeneralSettings& gs, const ComputeDrawingSettings& cds) { + if (exclDisks.empty()) return componentShape; + + std::vector> lineCovering; + std::vector> arcCovering; + + for (const auto& d : exclDisks) { + std::vector inter; + for (const auto& boundaryPart : boundaryParts) { + intersection(boundaryPart, circleToCSPolygon(d), std::back_inserter(inter), false, true); + } + bool coversLine = true; + for (const auto& p : inter) { + if (!isStraight(p)) { + coversLine = false; + break; + } + } + if (coversLine) { + lineCovering.push_back(d); + } else { + arcCovering.push_back(d); + } + } + + auto dr = gs.dilationRadius(); + // Smoothing radius + auto sr = dr / 5; + + std::vector> expandedLineCoveringDisks; + for (const auto& d : lineCovering) { + expandedLineCoveringDisks.emplace_back(approximate(d.center()), squared(sqrt(CGAL::to_double(d.squared_radius())) + sr)); + } + auto diskComponents = connectedDisks(expandedLineCoveringDisks); + + auto nearestOnBoundary = [&boundaryParts](const Point& point) { + std::optional minSqrdDist; + std::optional closest; + for (const auto& bp : boundaryParts) { + auto n = nearest(bp, point); + auto sqrdDist = CGAL::square(n.x() - point.x()) + CGAL::square(n.y() - point.y()); + if (!minSqrdDist.has_value() || sqrdDist < *minSqrdDist) { + minSqrdDist = sqrdDist; + closest = n; + } + } + return *closest; + }; + + std::vector cuts; + std::vector> rectangleCutDisks; + for (const auto& comp : diskComponents) { + std::vector> disks; + for (const auto& [i, _] : comp) { + disks.push_back(lineCovering[i]); + } + auto hull = approximateConvexHull(disks); + cuts.push_back(hull); + + // Only cut out a rectangle to the nearest disk of the component + std::optional> minDist; + std::optional> closestDisk; + + for (const auto& d : disks) { + if (inside(componentShape, d.center())) { + auto n = nearestOnBoundary(d.center()); + auto dist = CGAL::squared_distance(approximateAlgebraic(n), approximate(d.center())); + if (!minDist.has_value() || dist < *minDist) { + minDist = dist; + closestDisk = d; + } + } + } + + if (closestDisk.has_value()) { + rectangleCutDisks.push_back(*closestDisk); + } + + CSTraitsBoolean traits; + if (!is_valid_unknown_polygon(hull, traits)) { + // export debug info to ipe. + IpeRenderer ipeRenderer; + ipeRenderer.addPainting([&hull, &disks](GeometryRenderer& renderer) { + renderer.setMode(GeometryRenderer::fill | GeometryRenderer::stroke); + renderer.setStroke(Color(0), 1.0); + renderer.setFill(Color(225, 225, 225)); + renderer.draw(renderPath(hull)); + for (const auto& d : disks) { + renderer.draw(d); + } + }); + ipeRenderer.save(std::filesystem::path("ch-debug.ipe")); + for (auto cit = hull.curves_begin(); cit != hull.curves_end(); ++cit) { + std::cerr << *cit << std::endl; + } + throw std::runtime_error("CH not simple; see ch-debug.ipe"); + } + } + + for (const auto& d : arcCovering) { + cuts.push_back(circleToCSPolygon(d)); + if (inside(componentShape, d.center())) { + rectangleCutDisks.push_back(d); + } + } + + CSPolygonSet polygonSet; + for (const auto& cut : cuts) { + polygonSet.join(cut); + } + + for (const auto& d : rectangleCutDisks) { + auto n = nearestOnBoundary(d.center()); + auto rect = thinRectangle(d.center(), n, gs.pointSize); + polygonSet.join(rect); + } + + for (const auto& d: inclDisks) { + polygonSet.difference(circleToCSPolygon(d)); + } + std::vector modifiedCuts; + polygonSet.polygons_with_holes(std::back_inserter(modifiedCuts)); + + CSPolygonSet polygonSet2(componentShape); + for (const auto& modifiedCut : modifiedCuts) { + polygonSet2.difference(modifiedCut); + } + + std::vector polys; + polygonSet2.polygons_with_holes(std::back_inserter(polys)); +; + CSPolygonWithHoles poly; + std::optional> max_area; + for (const auto& cp : polys) { + auto a = area(cp); + if (!max_area.has_value() || a > *max_area) { + max_area = a; + poly = cp; + } + } + + // If poly has holes ignore them (that is, cut them out as well). + auto outer = poly.outer_boundary(); + return outer; +} + +CSPolyline associatedBoundary(const CSPolygon& component, const CSPolygon& morphedComponent, const CSPolyline& boundaryPart) { + std::vector morphed_xm_curves; + for (auto cit = morphedComponent.curves_begin(); cit != morphedComponent.curves_end(); ++cit) { + morphed_xm_curves.push_back(*cit); + } + auto boundaryPartStart = boundaryPart.curves_begin()->source(); + auto endIt = boundaryPart.curves_end(); + --endIt; + auto boundaryPartEnd = endIt->target(); + int startIndex = -1; + int endIndex = -1; + + for (int i = 0; i < morphed_xm_curves.size(); i++) { + auto c = morphed_xm_curves[i]; + if (morphed_xm_curves[i].source() == boundaryPartStart) { + startIndex = i; + } + if (morphed_xm_curves[i].target() == boundaryPartEnd) { + endIndex = i; + } + } + + for (int i = 0; i < morphed_xm_curves.size() && (startIndex < 0 || endIndex < 0); i++) { + const auto& c = morphed_xm_curves[i]; + if (startIndex < 0 && liesOn(c.source(), component).has_value() && !liesOn(c, component)) { + startIndex = i; + } + if (endIndex < 0 && !liesOn(c, component) && liesOn(c.target(), component).has_value()) { + endIndex = i; + } + } + + assert(startIndex >= 0); + assert(endIndex >= 0); + + std::vector xm_curves; + + if (endIndex >= startIndex) { + for (int i = startIndex; i <= endIndex; ++i) { + xm_curves.push_back(morphed_xm_curves[i]); + } + } else { + for (int i = startIndex; i < morphed_xm_curves.size(); ++i) { + xm_curves.push_back(morphed_xm_curves[i]); + } + for (int i = 0; i <= endIndex; ++i) { + xm_curves.push_back(morphed_xm_curves[i]); + } + } + + return {xm_curves.begin(), xm_curves.end()}; +} + +/// Returns parts of the boundary of c that originate from i. +/// This function assumes that some part of the boundary, but not all of the boundary, originates from i. +std::vector boundaryParts(const Component& c, int i) { + std::vector ccbs; + std::copy(c.outer_ccbs_begin(), c.outer_ccbs_end(), std::back_inserter(ccbs)); + std::copy(c.inner_ccbs_begin(), c.inner_ccbs_end(), std::back_inserter(ccbs)); + + std::vector polylines; + + for (const auto& ccb : ccbs) { + boundaryParts(ccb, i, std::back_inserter(polylines)); + } + + return polylines; +} + +/// Returns parts of the boundary of c that originate from i. +/// This function assumes that some part of the boundary, but not all of the boundary, originates from i. +std::vector boundaryParts(FaceH fh, int i) { + std::vector ccbs; + std::copy(fh->outer_ccbs_begin(), fh->outer_ccbs_end(), std::back_inserter(ccbs)); + std::copy(fh->inner_ccbs_begin(), fh->inner_ccbs_end(), std::back_inserter(ccbs)); + + std::vector polylines; + + for (const auto& ccb : ccbs) { + boundaryParts(ccb, i, std::back_inserter(polylines)); + } + + return polylines; +} + +// todo: check if small circular arcs should be allowed +bool isStraight(const CSPolyline& polyline) { + for (auto cit = polyline.curves_begin(); cit != polyline.curves_end(); ++cit) { + if (cit->is_circular()) return false; + } + return true; +} + +/// The inclusion and exclusion disks for component \ref c when \ref i is stacked on top of \ref js. +IncludeExcludeDisks +DilatedPatternDrawing::includeExcludeDisks(int i, const std::unordered_set& js, const Component& c) const { + std::vector> ptsI; + const auto& catPointsI = m_dilated[i].catPoints(); + std::transform(catPointsI.begin(), catPointsI.end(), std::back_inserter(ptsI), [](const CatPoint& catPoint) { + return Point(catPoint.point.x(), catPoint.point.y()); + }); + std::vector> ptsJs; + for (int j : js) { + const auto& catPointsJ = m_dilated[j].catPoints(); + std::transform(catPointsJ.begin(), catPointsJ.end(), std::back_inserter(ptsJs), + [](const CatPoint& catPoint) { + return Point(catPoint.point.x(), catPoint.point.y()); + }); + } + + auto rSqrd = m_gs.dilationRadius() * m_gs.dilationRadius(); + auto result = approximateGrowCircles(ptsI, ptsJs, rSqrd, rSqrd * m_cds.cutoutRadiusFactor * m_cds.cutoutRadiusFactor); + + std::vector> relevantExclusionDisks; + for (const auto& d : result.second) { + if (CGAL::do_intersect(circleToCSPolygon(d), ccb_to_polygon(c.outer_ccb()))) { + relevantExclusionDisks.push_back(d); + } + } + return {result.first, relevantExclusionDisks}; +} + +/// The inclusion and exclusion disks for component \ref c when \ref i is stacked on top of \ref j. +IncludeExcludeDisks +DilatedPatternDrawing::includeExcludeDisks(int i, int j, const Component& c) const { + std::unordered_set js({j}); + return includeExcludeDisks(i, js, c); +} + +std::shared_ptr DilatedPatternDrawing::computePreference(int i, int j, const Component& c) { + // The preference indicates the relation R in iRj. + // If R is Order::GREATER then i > j and i is preferred to be on top of j. + auto pref = Order::EQUAL; + + // 3. Prefer to cover a line segment over covering a circular arc. + auto circArcIsCovered = [](const std::vector& bps) { + for (auto& polyline : bps) { + if (!isStraight(polyline)) { + return true; + } + } + return false; + }; + + // If j is stacked over i then it covers all edges of i. + // Check if any of them is a circular arc. + auto bpi = boundaryParts(c, i); + assert(!bpi.empty()); + bool iCircArcIsCovered = circArcIsCovered(bpi); + // Vice versa. + auto bpj = boundaryParts(c, j); + assert(!bpj.empty()); + bool jCircArcIsCovered = circArcIsCovered(bpj); + + if (iCircArcIsCovered && !jCircArcIsCovered) { + pref = Order::GREATER; + } + if (!iCircArcIsCovered && jCircArcIsCovered){ + pref = Order::SMALLER; + } + + // 2. Prefer to indent a line segment over indenting a circular arc. + // Disks that would be cut out of i to expose points in j. + auto jExclusion = includeExcludeDisks(i, j, c).exclude; + // Disks that would be cut out of j to expose points in i. + auto iExclusion = includeExcludeDisks(j, i, c).exclude; + + auto circularIndented = [](const std::vector>& exclusionDisks, const std::vector& bps) { + for (const auto& d : exclusionDisks) { + if (d.squared_radius() <= 0) continue; + for (auto& polyline : bps) { + auto inters = intersection(polyline, circleToCSPolygon(d), true); + for (auto& inter : inters) { + if (!isStraight(inter)) { + return true; + } + } + } + } + return false; + }; + bool iCircularIndented = circularIndented(jExclusion, bpi); + bool jCircularIndented = circularIndented(iExclusion, bpj); + + if (iCircularIndented && !jCircularIndented) { + pref = Order::SMALLER; + } + if (!iCircularIndented && jCircularIndented){ + pref = Order::GREATER; + } + + // 1. Prefer to avoid few points over many points. + // Fewer disks would be cut out of i than out of j, so prefer to stack i on top of j + if (jExclusion.size() < iExclusion.size()) { + pref = Order::GREATER; + } + if (iExclusion.size() < jExclusion.size()) { + pref = Order::SMALLER; + } + + return std::make_shared(i, j, pref, pref); +} + +Color whiten(const Color& color, double a) { + return {static_cast(255 * a + color.r * (1-a)), static_cast(255 * a + color.g * (1-a)), static_cast(255 * a + color.b * (1-a))}; +} + +void DilatedPatternDrawing::drawFaceFill(FaceH fh, renderer::GeometryRenderer& renderer, + const GeneralSettings& gs, const DrawSettings& ds) const { + auto& d = fh->data(); + for (int i : d.ordering) { + renderer.setMode(GeometryRenderer::fill | GeometryRenderer::stroke); + renderer.setFill(whiten(ds.colors[m_dilated[i].category()], ds.whiten)); + renderer.setStroke(whiten(ds.colors[m_dilated[i].category()], ds.whiten), ds.contourStrokeWeight(gs) / 1.5, true); + if (!d.morphedFace.contains(i)) { + auto poly = face_to_polygon(*fh); + renderer.draw(renderPath(poly)); + } else { + for (const auto& p : d.morphedFace[i]) { + renderer.draw(renderPath(p)); + } + } + } +} + +void DilatedPatternDrawing::drawFaceStroke(FaceH fh, renderer::GeometryRenderer& renderer, + const GeneralSettings& gs, const DrawSettings& ds) const { + auto& d = fh->data(); + for (int index = 0; index < d.ordering.size(); ++index) { + int i = d.ordering[index]; + + std::vector polylines; + if (d.morphedEdges[i].empty()) { + auto bps = boundaryParts(fh, i); + std::copy(bps.begin(), bps.end(), std::back_inserter(polylines)); + } + std::copy(d.morphedEdges[i].begin(), d.morphedEdges[i].end(), std::back_inserter(polylines)); + + + std::vector modifiedPolylines; + if (index == d.ordering.size() - 1) { + modifiedPolylines = polylines; + } else { + bool nothingVisible = false; + auto poly = face_to_polygon(*fh); + auto bb = poly.bbox(); + Rectangle bbX(bb.xmin() - 1, bb.ymin() - 1, bb.xmax() + 1, bb.ymax() + 1); + std::vector xm_cs; + for (int k = 0; k < 4; ++k) { + xm_cs.emplace_back(bbX.vertex(k), bbX.vertex((k + 1) % 4)); + } + CSPolygon bbXPoly(xm_cs.begin(), xm_cs.end()); + + CSPolygonSet polySet(bbXPoly); + for (int higherIndex = index + 1; higherIndex < d.ordering.size(); ++higherIndex) { + int j = d.ordering[higherIndex]; + // Pattern j is stacked on top of i and will cover the stroke of shape i. + if (!d.morphedFace.contains(j)) { + nothingVisible = true; + break; + } else { + for (const auto& p : d.morphedFace[j]) { + polySet.difference(p); + } + } + } + + if (nothingVisible) continue; + + std::vector polygonsWithHoles; + polySet.polygons_with_holes(std::back_inserter(polygonsWithHoles)); + + for (const auto& polyline : polylines) { + for (const auto& polygon : polygonsWithHoles) { + intersection(polyline, polygon, std::back_inserter(modifiedPolylines), false, false); + } + } + } + + renderer.setMode(GeometryRenderer::stroke); + renderer.setStroke(Color{0, 0, 0}, ds.contourStrokeWeight(gs), true); + for (const auto& polyline : modifiedPolylines) { + renderer.draw(renderPath(polyline)); + } + } +} + +std::optional> DilatedPatternDrawing::totalStackingOrder() const { + std::vector> relations; + std::vector origins; + for (int i = 0; i < m_dilated.size(); ++i) { + origins.push_back(i); + for (const auto& f : m_iToFaces.at(i)) { + for (const auto& r : f->data().relations) { + if (std::find_if(relations.begin(), relations.end(), [&r](const auto& it) { return it->left == r->left && it->right == r->right; }) == relations.end()) { + relations.push_back(r); + } + assert(r->ordering != Order::EQUAL); + } + } + } + return computeTotalOrder(origins, relations); +} + +std::vector +DilatedPatternDrawing::intersectionComponents(int i, int j) const { + return connectedComponents(m_arr, [i, j](FaceH fh) { + const auto& origins = fh->data().origins; + return std::find(origins.begin(), origins.end(), i) != origins.end() && + std::find(origins.begin(), origins.end(), j) != origins.end(); + }); +} + +std::vector +DilatedPatternDrawing::intersectionComponents(int i) const { + return connectedComponents(m_arr, [i](FaceH fh) { + const auto& origins = fh->data().origins; + return origins.size() > 1 && std::find(origins.begin(), origins.end(), i) != origins.end(); + }); +} + +std::vector> DilatedPatternDrawing::hyperedges() const { + std::vector interesting; + + for (auto fit = m_arr.faces_begin(); fit != m_arr.faces_end(); ++fit) { + if (fit->data().origins.size() >= 2) { + interesting.push_back(fit); + } + } + + std::sort(interesting.begin(), interesting.end(), [](const auto& fh1, const auto& fh2) { + return fh1->data().origins.size() < fh2->data().origins.size(); + }); + + std::vector>> hyperedgesGrouped; + + std::vector> currentGroup; + std::optional lastSize; + for (const auto& fh : interesting) { + if (lastSize.has_value() && fh->data().origins.size() != *lastSize && !currentGroup.empty()) { + hyperedgesGrouped.push_back(currentGroup); + currentGroup.clear(); + } + auto he = std::make_shared(fh->data().origins, fh->data().relations); + for (auto& r : he->relations) { + r->hyperedges.push_back(he); + } + currentGroup.push_back(he); + lastSize = fh->data().origins.size(); + } + if (!currentGroup.empty()) { + hyperedgesGrouped.push_back(currentGroup); + } + + std::vector>> trashCan; + for (int i = 0; i < hyperedgesGrouped.size(); i++) { + if (i + 1 >= hyperedgesGrouped.size()) break; + const auto& group = hyperedgesGrouped[i]; + for (const auto& hyperedge : group) { + for (const auto& larger : hyperedgesGrouped[i+1]) { + bool fullyContained = true; + for (const auto& r : hyperedge->relations) { + if (std::find(larger->relations.begin(), larger->relations.end(), r) == larger->relations.end()) { + fullyContained = false; + } + } + if (fullyContained) { + // store hyperedge for removal + trashCan.emplace_back(i, hyperedge); + break; + } + } + } + } + + for (const auto& [i, r] : trashCan) { + auto& group = hyperedgesGrouped[i]; + group.erase(std::remove(group.begin(), group.end(), r), group.end()); + } + + std::vector> hyperedges; + for (const auto& group : hyperedgesGrouped) { + std::copy(group.begin(), group.end(), std::back_inserter(hyperedges)); + } + + return hyperedges; +} + +std::optional> computeTotalOrder(const std::vector& origins, const std::vector>& relations) { + if (relations.empty()) { + return origins; + } + + struct Vertex { + int i; + std::vector neighbors; + bool hasIncoming = false; + int mark = 0; + }; + + std::unordered_map vertices; + + for (int i : origins) { + vertices[i] = Vertex{i}; + } + + for (const auto& r : relations) { + if (r->ordering == Order::EQUAL) continue; + auto& u = vertices.at(r->left); + auto& v = vertices.at(r->right); + + if (r->ordering == Order::SMALLER) { + v.neighbors.push_back(&u); + u.hasIncoming = true; + } else { + u.neighbors.push_back(&v); + v.hasIncoming = true; + } + } + + std::vector ordering; + + std::function visit; + + visit = [&ordering, &visit](Vertex& u) { + if (u.mark == 2) return true; + if (u.mark == 1) return false; + + u.mark = 1; + + for (auto& v : u.neighbors) { + bool success = visit(*v); + if (!success) return false; + } + + u.mark = 2; + ordering.push_back(u.i); + return true; + }; + + for (auto& [_, v] : vertices) { + if (!v.hasIncoming) { + bool success = visit(v); + if (!success) + return std::nullopt; + } + } + + if (ordering.size() != origins.size()) { + return std::nullopt; + } + + return ordering; +} + +std::optional> getRelationOrder(const Hyperedge& e) { + return computeTotalOrder(e.origins, e.relations); +} + +void setRelationOrder(Hyperedge& e, const std::vector& ordering) { + for (auto& r : e.relations) { + int i = std::find(ordering.begin(), ordering.end(), r->left) - ordering.begin(); + int j = std::find(ordering.begin(), ordering.end(), r->right) - ordering.begin(); + if (i < j) { + r->ordering = Order::SMALLER; + } else { + r->ordering = Order::GREATER; + } + } +} + +std::string to_string(Order ord) { + if (ord == Order::SMALLER) { + return "<"; + } else if (ord == Order::EQUAL) { + return "="; + } else { + return ">"; + } +} + +SimpleSetsPainting::SimpleSetsPainting(const DilatedPatternDrawing& dpd, const DrawSettings& ds) + : m_ds(ds), m_dpd(dpd) {} + +void SimpleSetsPainting::paint(renderer::GeometryRenderer& renderer) const { + auto stackingOrder = m_dpd.totalStackingOrder(); + // If there is a stacking order, draw the complete patterns stacked in that order + if (stackingOrder.has_value()) { + for (int i : *stackingOrder) { + auto comps = connectedComponents(m_dpd.m_arr, [i](const FaceH& fh) { + const auto& ors = fh->data().origins; + return std::find(ors.begin(), ors.end(), i) != ors.end(); + }); + assert(comps.size() == 1); + const auto& comp = comps[0]; + + std::vector boundaryPieces; + for (auto fit = comp.faces_begin(); fit != comp.faces_end(); ++fit) { + auto& data = fit->data(); + if (data.morphedFace.contains(i)) { + std::copy(data.morphedEdges[i].begin(), data.morphedEdges[i].end(), std::back_inserter(boundaryPieces)); + } else { + FaceH fh = fit.handle(); + auto bps = boundaryParts(fh, i); + std::copy(bps.begin(), bps.end(), std::back_inserter(boundaryPieces)); + } + } + + // no order or hash on OneRootPoint :( +// std::map sourceToI; +// for (int j = 0; j < boundaryPieces.size(); ++j) { +// auto& bp = boundaryPieces[j]; +// sourceToI[bp.curves_begin()->source()] = j; +// } + std::vector xm_curves; + + std::copy(boundaryPieces[0].curves_begin(), boundaryPieces[0].curves_end(), std::back_inserter(xm_curves)); + int count = 1; + while (count < boundaryPieces.size()) { + auto head = xm_curves.back().target(); + for (const auto& bp : boundaryPieces) { + if (bp.curves_begin()->source() == head) { + std::copy(bp.curves_begin(), bp.curves_end(), std::back_inserter(xm_curves)); + } + } + ++count; + } + + CSPolygon csPolygon(xm_curves.begin(), xm_curves.end()); + + renderer.setMode(GeometryRenderer::fill | GeometryRenderer::stroke); + renderer.setFill(whiten(m_ds.colors[m_dpd.m_dilated[i].category()], m_ds.whiten)); + renderer.setStroke(Color{0, 0, 0}, m_ds.contourStrokeWeight(m_dpd.m_gs), true); + renderer.draw(renderPath(csPolygon)); + } + } else { + // If there is no stacking order, draw each face of the arrangement separately + for (auto fit = m_dpd.m_arr.faces_begin(); fit != m_dpd.m_arr.faces_end(); ++fit) { + if (fit->is_unbounded()) + continue; + m_dpd.drawFaceFill(fit.ptr(), renderer, m_dpd.m_gs, m_ds); + } + for (auto fit = m_dpd.m_arr.faces_begin(); fit != m_dpd.m_arr.faces_end(); ++fit) { + if (fit->is_unbounded()) + continue; + m_dpd.drawFaceStroke(fit.ptr(), renderer, m_dpd.m_gs, m_ds); + } + } + + // Draw points + const auto& gs = m_dpd.m_gs; + renderer.setStroke(Color{0, 0, 0}, m_ds.pointStrokeWeight(gs), true); + renderer.setFillOpacity(255); + renderer.setMode(renderer::GeometryRenderer::fill | renderer::GeometryRenderer::stroke); + for (const auto& dp : m_dpd.m_dilated) { + for (const auto& cp : dp.catPoints()) { + renderer.setFill(m_ds.getColor(cp.category)); + renderer.draw(Circle{cp.point, gs.pointSize * gs.pointSize}); + } + } +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/drawing_algorithm.h b/cartocrow/simplesets/drawing_algorithm.h new file mode 100644 index 00000000..5c4ee69b --- /dev/null +++ b/cartocrow/simplesets/drawing_algorithm.h @@ -0,0 +1,383 @@ +#ifndef CARTOCROW_DRAWING_ALGORITHM_H +#define CARTOCROW_DRAWING_ALGORITHM_H + +#include "partition.h" +#include "dilated/dilated_poly.h" +#include "../renderer/geometry_painting.h" +#include + +#include + +namespace cartocrow::simplesets { +enum Order { + SMALLER = -1, + EQUAL = 0, + GREATER = 1, +}; + +std::string to_string(Order ord); + +std::ostream& operator<<(std::ostream& out, const Order& o); + +struct Hyperedge; + +struct Relation { + int left; + int right; + Order preference; + Order ordering; + std::vector> hyperedges; + Relation(int left, int right, Order preference, Order ordering) : + left(left), right(right), preference(preference), ordering(ordering) {} +}; + +bool operator==(const Relation& lhs, const Relation& rhs); + +std::ostream& operator<<(std::ostream& out, const Relation& r); + +struct Hyperedge { + std::vector origins; + std::vector> relations; + Hyperedge(std::vector origins, std::vector> relations) : + origins(origins), relations(relations) {}; +}; + +std::optional> getRelationOrder(const Hyperedge& e); +void setRelationOrder(Hyperedge& e, const std::vector& ordering); +std::optional> computeTotalOrder(const std::vector& origins, const std::vector>& relations); + +bool operator==(const Hyperedge& lhs, const Hyperedge& rhs); + +struct FaceData { + std::vector origins; + std::vector> relations; + std::vector ordering; + std::unordered_map> morphedEdges; + std::unordered_map> morphedFace; +}; + +// todo: edge case where edges of dilated patterns overlap, so a half-edge may have multiple origins. +struct HalfEdgeData { + std::vector origins; +}; + +struct VertexData { + +}; + +using DilatedPatternArrangement = + CGAL::Arrangement_with_history_2>; + +//DilatedPatternArrangement +//dilateAndArrange(const Partition& partition, const GeneralSettings& gs, const ComputeDrawingSettings& cds); + +using Face = DilatedPatternArrangement::Face; +using FaceH = DilatedPatternArrangement::Face_handle; +using FaceCH = DilatedPatternArrangement::Face_const_handle; +using VertexH = DilatedPatternArrangement::Vertex_handle; +using HalfEdgeH = DilatedPatternArrangement::Halfedge_handle; + +class Component { + public: + typedef DilatedPatternArrangement Arr; +// typedef std::vector::iterator Face_iterator; + typedef Arr::Size Size; + + class ComponentCcbCirculator { + private: + using Self = ComponentCcbCirculator; + Arr::Halfedge_handle m_halfedge; + std::function m_in_component; + + public: + using iterator_category = std::bidirectional_iterator_tag; + using value_type = Arr::Halfedge; + using difference_type = std::ptrdiff_t; + using pointer = Arr::Halfedge_handle; + using reference = value_type&; + + ComponentCcbCirculator(Arr::Halfedge_handle halfedge, std::function in_component) + : m_halfedge(halfedge), m_in_component(std::move(in_component)) {}; + + value_type operator*() const { + return *m_halfedge; + } + + pointer operator->() const { + return m_halfedge; + } + + pointer ptr() const { + return m_halfedge; + } + + pointer handle() const { + return m_halfedge; + } + + Self& operator++() { + m_halfedge = m_halfedge->next(); + while (m_in_component(m_halfedge->twin()->face())) { + m_halfedge = m_halfedge->twin()->next(); + } + return *this; + }; + + Self operator++(int) { + Self tmp = *this; + this->operator++(); + return tmp; + } + + Self& operator--() { + m_halfedge = m_halfedge->prev(); + while (m_in_component(m_halfedge->twin()->face())) { + m_halfedge = m_halfedge->twin()->prev(); + } + return *this; + }; + + Self operator--(int) { + Self tmp = *this; + this->operator--(); + return tmp; + } + + bool operator==(const Self& other) const { + return m_halfedge == other.m_halfedge; + } + + bool operator!=(const Self& other) const { + return m_halfedge != other.m_halfedge; + } + }; + + class Face_const_iterator { + private: + using Self = Face_const_iterator; + std::vector::const_iterator m_faceHandleIterator; + + public: + using iterator_category = std::input_iterator_tag; + using value_type = Arr::Face; + using difference_type = std::ptrdiff_t; + using pointer = Arr::Face_handle; + using reference = value_type&; + + Face_const_iterator(std::vector::const_iterator faceHandleIterator) : + m_faceHandleIterator(faceHandleIterator) {}; + + value_type operator*() const { + return *(*m_faceHandleIterator); + } + + pointer operator->() const { + return *m_faceHandleIterator; + } + + Self& operator++() { + ++m_faceHandleIterator; + return *this; + }; + + Self operator++(int) { + Self tmp = *this; + this->operator++(); + return tmp; + } + + pointer ptr() const { + return *m_faceHandleIterator; + } + + pointer handle() const { + return *m_faceHandleIterator; + } + + bool operator==(const Self& other) const { + return m_faceHandleIterator == other.m_faceHandleIterator; + } + + bool operator!=(const Self& other) const { + return m_faceHandleIterator != other.m_faceHandleIterator; + } + }; + + Component(std::vector faces, std::vector boundaryEdges, std::function inComponent); + + bool has_outer_ccb() const { + return !m_outer_ccbs.empty(); + } + Size number_of_outer_ccbs() const { + return m_outer_ccbs.size(); + } + std::vector::const_iterator outer_ccbs_begin() const { + return m_outer_ccbs.cbegin(); + } + std::vector::const_iterator outer_ccbs_end() const { + return m_outer_ccbs.cend(); + } + + ComponentCcbCirculator outer_ccb() const { + return m_outer_ccbs[0]; + } + std::vector::const_iterator inner_ccbs_begin() const { + return m_inner_ccbs.cbegin(); + } + std::vector::const_iterator inner_ccbs_end() const { + return m_inner_ccbs.cend(); + } + Size number_of_inner_ccbs() const { + return m_inner_ccbs.size(); + } + + // Hole is an alias for inner_ccb + std::vector::const_iterator holes_begin() const { + return inner_ccbs_begin(); + } + std::vector::const_iterator holes_end() const { + return inner_ccbs_end(); + } + Size number_of_holes() { + return number_of_inner_ccbs(); + } + + Face_const_iterator faces_begin() const { + return {m_faces.cbegin()}; + } + Face_const_iterator faces_end() const { + return {m_faces.cend()}; + } + + private: + std::vector m_faces; + std::function m_in_component; + std::vector m_outer_ccbs; + std::vector m_inner_ccbs; +}; + +template +void boundaryParts(Ccb ccb, int i, OutputIterator out) { + // First find a half-edge at the start of a 'part' of this CCB (connected component of the boundary). + auto circ = ccb; + + auto originates_from = [](Ccb circ, int i) { + auto& os = circ->data().origins; + return std::find(os.begin(), os.end(), i) != os.end(); + }; + + bool found = true; + while (!originates_from(circ, i)) { + ++circ; + if (circ == ccb) { + found = false; + break; + } + } + + if (!found) { + return; + } + + auto start = circ; + do { + auto prev = circ; + --prev; + if (originates_from(prev, i)) { + circ = prev; + } else { + break; + } + } while (circ != start); + assert(originates_from(circ, i)); + + Traits traits; + auto opposite = traits.construct_opposite_2_object(); + + // Next, make a polyline for every connected part of the boundary that originates from i. + std::vector xm_curves; + auto curr = circ; + do { + if (originates_from(curr, i)) { + auto& curve = curr->curve(); + if (curr->source()->point() == curve.source()) { + xm_curves.push_back(curve); + } else { + xm_curves.push_back(opposite(curve)); + } + } else { + if (!xm_curves.empty()) { + ++out = CSPolyline(xm_curves.begin(), xm_curves.end()); + xm_curves.clear(); + } + } + } while (++curr != circ); + if (!xm_curves.empty()) { + ++out = CSPolyline(xm_curves.begin(), xm_curves.end()); + xm_curves.clear(); + } +} + +std::vector boundaryParts(const Component& c, int i); +std::vector boundaryParts(FaceH f, int i); +std::vector>> originCcbs(const Component& c); + +// exposed for debugging purposes +std::vector>>> connectedDisks(const std::vector>& disks); +CSPolygon thinRectangle(const Point& p, const OneRootPoint& n, const Number& w); + +CSPolygon morph(const std::vector& boundaryParts, const CSPolygon& componentShape, const std::vector>& inclDisks, + const std::vector>& exclDisks, const GeneralSettings& gs, const ComputeDrawingSettings& cds); + +CSPolyline associatedBoundary(const CSPolygon& component, const CSPolygon& morphedComponent, const CSPolyline& boundaryPart); + +struct IncludeExcludeDisks { + std::vector> include; + std::vector> exclude; +}; + +bool isStraight(const CSPolyline& polyline); + +class DilatedPatternDrawing { + public: + DilatedPatternDrawing(const Partition& partition, const GeneralSettings& gs, const ComputeDrawingSettings& cds); + + std::vector intersectionComponents(int i) const; + std::vector intersectionComponents(int i, int j) const; + std::shared_ptr computePreference(int i, int j, const Component& c); + + IncludeExcludeDisks includeExcludeDisks(int i, int j, const Component& c) const; + IncludeExcludeDisks includeExcludeDisks(int i, const std::unordered_set& js, const Component& c) const; + + std::vector> hyperedges() const; + + void drawFaceFill(FaceH fh, renderer::GeometryRenderer& renderer, + const GeneralSettings& gs, const DrawSettings& ds) const; + void drawFaceStroke(FaceH fh, renderer::GeometryRenderer& renderer, + const GeneralSettings& gs, const DrawSettings& ds) const; + std::optional> totalStackingOrder() const; + + DilatedPatternArrangement m_arr; + std::unordered_map> m_iToFaces; + std::map m_curve_to_origin; + std::vector m_dilated; + const GeneralSettings& m_gs; + const ComputeDrawingSettings& m_cds; + + private: +}; + +class SimpleSetsPainting : public renderer::GeometryPainting { + public: + SimpleSetsPainting(const DilatedPatternDrawing& dpd, const DrawSettings& ds); + void paint(renderer::GeometryRenderer& renderer) const override; + + private: + const DrawSettings& m_ds; + const DilatedPatternDrawing& m_dpd; + +}; +} + +#endif //CARTOCROW_DRAWING_ALGORITHM_H diff --git a/cartocrow/simplesets/general_polyline.h b/cartocrow/simplesets/general_polyline.h new file mode 100644 index 00000000..33f06e90 --- /dev/null +++ b/cartocrow/simplesets/general_polyline.h @@ -0,0 +1,37 @@ +#ifndef CARTOCROW_GENERAL_POLYLINE_H +#define CARTOCROW_GENERAL_POLYLINE_H + +#include + +template +class General_polyline_2 { + public: + typedef typename std::vector::iterator Curve_iterator; + typedef typename std::vector::const_iterator Curve_const_iterator; + + template + General_polyline_2(InputIterator begin, InputIterator end) { + m_xm_curves = std::vector(begin, end); + } + + Curve_iterator curves_begin() { + return m_xm_curves.begin(); + } + Curve_iterator curves_end() { + return m_xm_curves.end(); + } + Curve_const_iterator curves_begin() const { + return m_xm_curves.cbegin(); + } + Curve_const_iterator curves_end() const { + return m_xm_curves.cend(); + } + [[nodiscard]] unsigned int size() const + { + return static_cast(m_xm_curves.size()); + } + private: + std::vector m_xm_curves; +}; + +#endif //CARTOCROW_GENERAL_POLYLINE_H diff --git a/cartocrow/simplesets/grow_circles.cpp b/cartocrow/simplesets/grow_circles.cpp new file mode 100644 index 00000000..d5945e17 --- /dev/null +++ b/cartocrow/simplesets/grow_circles.cpp @@ -0,0 +1,102 @@ +#include "grow_circles.h" + +namespace cartocrow::simplesets { +std::pair>, std::vector>> +approximateGrowCircles(const std::vector>& points1, + const std::vector>& points2, + const Number& maxSquaredRadius1, + const Number& maxSquaredRadius2) { + std::vector growingCircles1; + std::transform(points1.begin(), points1.end(), std::back_inserter(growingCircles1), [](const Point& pt) { + return GrowingCircle{pt, 0}; + }); + std::vector growingCircles2; + std::transform(points2.begin(), points2.end(), std::back_inserter(growingCircles2), [](const Point& pt) { + return GrowingCircle{pt, 0}; + }); + + // if statement and while(true) combined. + while (!growingCircles1.empty() && !growingCircles2.empty()) { + bool done = true; + for (const auto& c : growingCircles1) { + if (!c.frozen) { + done = false; + break; + } + } + for (const auto& c : growingCircles2) { + if (!c.frozen) { + done = false; + break; + } + } + if (done) break; + + // Note that all variables representing distances contain squared distances. + std::optional> minDist; + std::optional> minPair; + for (auto& c1 : growingCircles1) { + for (auto& c2 : growingCircles2) { + if (c1.frozen && c2.frozen) continue; + auto centerDist = squared_distance(c1.center, c2.center); + + Number growDist; + if (!c1.frozen && !c2.frozen) { + growDist = centerDist / 4; + } else if (c1.frozen) { + growDist = centerDist + c1.squaredRadius - 2 * sqrt(CGAL::to_double(centerDist)) * sqrt(CGAL::to_double(c1.squaredRadius)); + } else { + growDist = centerDist + c2.squaredRadius - 2 * sqrt(CGAL::to_double(centerDist)) * sqrt(CGAL::to_double(c2.squaredRadius)); + } + + if (!minDist.has_value() || growDist < *minDist) { + minDist = growDist; + minPair = {&c1, &c2}; + } + } + } + + auto& [c1, c2] = *minPair; + + if (!c1->frozen && !c2->frozen) { + auto d = std::min(*minDist, std::min(maxSquaredRadius1, maxSquaredRadius2)); + if (d == *minDist) { + c1->frozen = true; + c2->frozen = true; + } else if (d == maxSquaredRadius1) { + c1->frozen = true; + } else { + assert(d == maxSquaredRadius2); + c2->frozen = true; + } + c1->squaredRadius = d; + c2->squaredRadius = d; + } else if (c1->frozen) { + c2->squaredRadius = std::min(*minDist, maxSquaredRadius2); + c2->frozen = true; + } else { + assert(c2->frozen); + c1->squaredRadius = std::min(*minDist, maxSquaredRadius1); + c1->frozen = true; + } + } + // Trivial case + if(growingCircles1.empty() || growingCircles2.empty()) { + for (auto& c : growingCircles1) { + c.squaredRadius = maxSquaredRadius1; + } + for (auto& c : growingCircles2) { + c.squaredRadius = maxSquaredRadius2; + } + } + + std::vector> circles1; + std::vector> circles2; + std::transform(growingCircles1.begin(), growingCircles1.end(), std::back_inserter(circles1), + [](const GrowingCircle& gc) { return Circle(gc.center, gc.squaredRadius); }); + std::transform(growingCircles2.begin(), growingCircles2.end(), std::back_inserter(circles2), + [](const GrowingCircle& gc) { return Circle(gc.center, gc.squaredRadius); }); + + return {circles1, circles2}; +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/grow_circles.h b/cartocrow/simplesets/grow_circles.h new file mode 100644 index 00000000..fbb01861 --- /dev/null +++ b/cartocrow/simplesets/grow_circles.h @@ -0,0 +1,17 @@ +#ifndef CARTOCROW_GROW_CIRCLES_H +#define CARTOCROW_GROW_CIRCLES_H + +#include "../core/core.h" + +namespace cartocrow::simplesets { +struct GrowingCircle { + Point center; + Number squaredRadius; + bool frozen = false; +}; + +std::pair>, std::vector>> +approximateGrowCircles(const std::vector>& points1, const std::vector>& points2, + const Number& maxSquaredRadius1, const Number& maxSquaredRadius2); +} +#endif //CARTOCROW_GROW_CIRCLES_H diff --git a/cartocrow/simplesets/helpers/approximate_convex_hull.cpp b/cartocrow/simplesets/helpers/approximate_convex_hull.cpp new file mode 100644 index 00000000..f7aaff91 --- /dev/null +++ b/cartocrow/simplesets/helpers/approximate_convex_hull.cpp @@ -0,0 +1,302 @@ +#include "approximate_convex_hull.h" +#include "cs_curve_helpers.h" +#include "cs_polygon_helpers.h" + +namespace cartocrow::simplesets { +Segment tangent(const Circle& c1, const Circle& c2) { + auto distSq = CGAL::squared_distance(c1.center(), c2.center()); + auto hyp = c2.center() - c1.center(); + auto c1r = sqrt(c1.squared_radius()); + auto c2r = sqrt(c2.squared_radius()); + auto adj = c1r - c2r; + auto a = hyp * adj; + auto b = hyp.perpendicular(CGAL::POSITIVE) * sqrt(distSq - adj * adj); + auto v1 = (a - b) / distSq; + auto v2 = (a + b) / distSq; + + Segment t1(c1.center() + v1 * c1r, c2.center() + v1 * c2r); + Segment t2(c1.center() + v2 * c1r, c2.center() + v2 * c2r); + + if (CGAL::orientation(t1.source(), t1.target(), c2.center()) == CGAL::COUNTERCLOCKWISE) { + return t1; + } else { + return t2; + } +} + +std::pair, Point> +tangentPoints(const Circle& c1, const Circle& c2) { + auto distSq = CGAL::squared_distance(c1.center(), c2.center()); + auto hyp = c2.center() - c1.center(); + auto c1r = sqrt(c1.squared_radius()); + auto c2r = sqrt(c2.squared_radius()); + auto adj = c1r - c2r; + auto a = hyp * adj; + auto b = hyp.perpendicular(CGAL::POSITIVE) * sqrt(distSq - adj * adj); + auto v1 = (a - b) / distSq; + auto v2 = (a + b) / distSq; + + Segment t1(c1.center() + v1 * c1r, c2.center() + v1 * c2r); + Segment t2(c1.center() + v2 * c1r, c2.center() + v2 * c2r); + + if (CGAL::orientation(t1.source(), t1.target(), c2.center()) == CGAL::COUNTERCLOCKWISE) { + return {t1.source(), t1.target()}; + } else { + return {t2.source(), t2.target()}; + } +} + +std::pair +tangentPoints(const RationalRadiusCircle& c1, const RationalRadiusCircle& c2) { + Number distSq = CGAL::squared_distance(c1.center, c2.center); + Vector hyp = c2.center - c1.center; + Number c1r = c1.radius; + Number c2r = c2.radius; + Number adj = c1r - c2r; + Vector a = hyp * adj; + Vector bV = hyp.perpendicular(CGAL::POSITIVE); + Number bSSqr = distSq - adj * adj; + CSTraits::CoordNT bx(0, bV.x(), bSSqr); + CSTraits::CoordNT by(0, bV.y(), bSSqr); + + CSTraits::CoordNT v1x = (a.x() - bx) / distSq; + CSTraits::CoordNT v1y = (a.y() - by) / distSq; + CSTraits::CoordNT v2x = (a.x() + bx) / distSq; + CSTraits::CoordNT v2y = (a.y() + by) / distSq; + + CSTraits::CoordNT t1sx = c1.center.x() + v1x * c1r; + CSTraits::CoordNT t1sy = c1.center.y() + v1y * c1r; + CSTraits::Point_2 t1s(t1sx, t1sy); + CSTraits::CoordNT t1tx = c2.center.x() + v1x * c2r; + CSTraits::CoordNT t1ty = c2.center.y() + v1y * c2r; + CSTraits::Point_2 t1t(t1tx, t1ty); + + CSTraits::CoordNT t2sx = c1.center.x() + v2x * c1r; + CSTraits::CoordNT t2sy = c1.center.y() + v2y * c1r; + CSTraits::Point_2 t2s(t2sx, t2sy); + CSTraits::CoordNT t2tx = c2.center.x() + v2x * c2r; + CSTraits::CoordNT t2ty = c2.center.y() + v2y * c2r; + CSTraits::Point_2 t2t(t2tx, t2ty); + CSTraits::Point_2 c2c(c2.center.x(), c2.center.y()); + + return {t1s, t1t}; +} + +RationalRadiusCircle approximateRadiusCircle(const Circle& circle) { + Number r = sqrt(CGAL::to_double(circle.squared_radius())); + Number rExact = r; + return {circle.center(), rExact}; +} + +// Adapted from https://github.com/CGAL/cgal/blob/38871d9b125c5513ff8d14f9562795aa12681b38/Minkowski_sum_2/include/CGAL/Minkowski_sum_2/Approx_offset_base_2.h +// This function falls under the following license: +//// Copyright (c) 2006-2008 Tel-Aviv University (Israel). +//// All rights reserved. +//// +//// This function is part of CGAL (www.cgal.org). +//// +//// $URL$ +//// $Id$ +//// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +//// +//// Author(s) : Ron Wein +//// Andreas Fabri +//// Laurent Rineau +//// Efi Fogel +std::variant, std::pair, Segment>> +algebraicCircleTangentToRationalSegments(const CSTraits::Point_2& p1, const CSTraits::Point_2& p2, + const RationalRadiusCircle& c1, const RationalRadiusCircle& c2) { + const auto& x1 = p1.x(); + const auto& y1 = p1.y(); + const auto& x2 = p2.x(); + const auto& y2 = p2.y(); + + auto delta_x = x2 - x1; + auto delta_y = y2 - y1; + auto sqr_d = CGAL::square(delta_x) + CGAL::square(delta_y); + // This is hacky + Number app_delta_x = CGAL::to_double(x2-x1); + Number app_delta_y = CGAL::to_double(y2-y1); + Number app_d = sqrt(CGAL::to_double(sqr_d)); + + auto d_app_err = sqr_d - CGAL::square(app_d); + auto dx_app_err = app_delta_x - delta_x; + auto dy_app_err = app_delta_y - delta_y; + auto sign_d_app_err = CGAL::sign(d_app_err); + auto sign_dx_app_err = CGAL::sign(dx_app_err); + auto sign_dy_app_err = CGAL::sign(dy_app_err); + + if (sign_d_app_err == CGAL::ZERO && sign_dx_app_err == CGAL::ZERO && sign_dy_app_err == CGAL::ZERO) { + auto tp1 = Point (c1.center.x() + c1.radius * app_delta_y / app_d, c1.center.y() + c1.radius * (-app_delta_x) / app_d); + auto tp2 = Point (c2.center.x() + c2.radius * app_delta_y / app_d, c2.center.y() + c2.radius * (-app_delta_x) / app_d); + + auto seg1 = Segment (tp1, tp2); + return seg1; + } else { + // This is hacky + if (CGAL::sign(app_delta_x) == CGAL::ZERO) { + app_delta_x += M_EPSILON; + } + if (CGAL::sign(app_delta_y) == CGAL::ZERO) { + app_delta_y += M_EPSILON; + } + + bool rotate_pi2 = false; + if (CGAL::compare (CGAL::abs(delta_x), + CGAL::abs(delta_y)) == CGAL::SMALLER) + { + rotate_pi2 = true; + + // Swap the delta_x and delta_y values. + auto tmp_app = app_delta_x; + app_delta_x = -app_delta_y; + app_delta_y = tmp_app; + } + + auto lower_tan_half_phi = (app_d - app_delta_y) / (-app_delta_x); + auto upper_tan_half_phi = (-app_delta_x) / (app_d + app_delta_y); + if (upper_tan_half_phi < lower_tan_half_phi) { + auto temp = lower_tan_half_phi; + lower_tan_half_phi = upper_tan_half_phi; + upper_tan_half_phi = temp; + } + + // This is hacky + lower_tan_half_phi -= M_EPSILON; + upper_tan_half_phi += M_EPSILON; + + auto sqr_tan_half_phi = CGAL::square (lower_tan_half_phi); + auto sin_phi = 2 * lower_tan_half_phi / (1 + sqr_tan_half_phi); + auto cos_phi = (1 - sqr_tan_half_phi) / (1 + sqr_tan_half_phi); + + Point tp1; + if (! rotate_pi2) + { + tp1 = Point (c1.center.x() + c1.radius*cos_phi, c1.center.y() + c1.radius*sin_phi); + } + else + { + tp1 = Point (c1.center.x() + c1.radius*sin_phi, c1.center.y() - c1.radius*cos_phi); + } + + sqr_tan_half_phi = CGAL::square (upper_tan_half_phi); + sin_phi = 2 * upper_tan_half_phi / (1 + sqr_tan_half_phi); + cos_phi = (1 - sqr_tan_half_phi) / (1 + sqr_tan_half_phi); + + Point tp2; + if (! rotate_pi2) + { + tp2 = Point (c2.center.x() + c2.radius*cos_phi, c2.center.y() + c2.radius*sin_phi); + } + else + { + tp2 = Point (c2.center.x() + c2.radius*sin_phi, c2.center.y() - c2.radius*cos_phi); + } + + auto l1 = Line(c1.center, tp1).perpendicular(tp1); + auto l2 = Line(c2.center, tp2).perpendicular(tp2); + + // Intersect the two lines. The intersection point serves as a common + // end point for the two line segments we are about to introduce. + auto obj = CGAL::intersection(l1, l2); + + Point mid_p; + auto assign_success = CGAL::assign (mid_p, obj); + assert(assign_success); + + auto seg1 = Segment (tp1, mid_p); + auto seg2 = Segment (mid_p, tp2); + + return std::pair(seg1, seg2); + } +} + +std::variant, std::pair, Segment>> +approximateTangent(const RationalRadiusCircle& c1, const RationalRadiusCircle& c2) { + auto [source, target] = tangentPoints(c1, c2); + return algebraicCircleTangentToRationalSegments(source, target, c1, c2); +} + +std::vector circlesOnConvexHull(const std::vector& circles) { + if (circles.size() == 1) { + return {circles.front()}; + } + + Apollonius apo; + std::unordered_map vToCircle; + for (const auto& c : circles) { + auto vh = apo.insert(ASite(c.center, c.radius)); + vToCircle[vh] = c; + } + if (apo.number_of_vertices() == 1) { + auto site = apo.finite_vertices_begin()->site(); + return {{site.point(), site.weight()}}; + } + auto circ = apo.incident_vertices(apo.infinite_vertex()); + auto curr = circ; + std::vector hullCircles; + do { + hullCircles.push_back(vToCircle[curr]); + } while (++curr != circ); + + std::reverse(hullCircles.begin(), hullCircles.end()); + + return hullCircles; +} + +/// Precondition: the circle centers are distinct. +CSPolygon approximateConvexHull(const std::vector>& circles) { + // todo: approximating circle radii may cause problems when two circles overlap in a single point and one is contained in the other. + // solution? filter out any circle that is contained in another, before approximating the radii. + if (circles.size() == 1) { + return circleToCSPolygon(circles.front()); + } + std::vector rrCircles; + for (const auto& c : circles) { + rrCircles.push_back(approximateRadiusCircle(c)); + } + auto hullCircles = circlesOnConvexHull(rrCircles); + if (hullCircles.size() == 1) { + for (const auto& c : circles) { + if (c.center() == hullCircles[0].center) { + return circleToCSPolygon(c); + } + } + } + + std::vector>> tangents; + + for (int i = 0; i < hullCircles.size(); ++i) { + auto& c1 = hullCircles[i]; + auto& c2 = hullCircles[(i + 1) % hullCircles.size()]; + auto segOrPair = approximateTangent(c1, c2); + std::vector> segs; + if (segOrPair.index() == 0) { + segs.push_back(std::get>(segOrPair)); + } else { + auto pair = std::get, Segment>>(segOrPair); + segs.push_back(pair.first); + segs.push_back(pair.second); + } + tangents.push_back(segs); + } + + std::vector xm_curves; + for (int i = 0; i < hullCircles.size(); ++i) { + auto& c1 = hullCircles[i]; + auto& c2 = hullCircles[(i + 1) % hullCircles.size()]; + auto& t1 = tangents[i]; + auto& t2 = tangents[(i + 1) % tangents.size()]; + for (const auto& piece : t1) { + Curve_2 curve(piece); + curveToXMonotoneCurves(curve, std::back_inserter(xm_curves)); + } + OneRootPoint t1End(t1.back().target().x(), t1.back().target().y()); + OneRootPoint t2Start(t2.front().source().x(), t2.front().source().y()); + Curve_2 arc(Circle(c2.center, c2.radius * c2.radius), t1End, t2Start); + curveToXMonotoneCurves(arc, std::back_inserter(xm_curves)); + } + + return {xm_curves.begin(), xm_curves.end()}; +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/helpers/approximate_convex_hull.h b/cartocrow/simplesets/helpers/approximate_convex_hull.h new file mode 100644 index 00000000..518da640 --- /dev/null +++ b/cartocrow/simplesets/helpers/approximate_convex_hull.h @@ -0,0 +1,38 @@ +#ifndef CARTOCROW_APPROXIMATE_CONVEX_HULL_H +#define CARTOCROW_APPROXIMATE_CONVEX_HULL_H + +#include "../types.h" +#include +#include +#include +#include + +namespace cartocrow::simplesets { +typedef CGAL::Apollonius_graph_traits_2 AT; +typedef CGAL::Apollonius_graph_2 Apollonius; +typedef Apollonius::Site_2 ASite; + + +struct RationalRadiusCircle { + Point center; + Number radius; +}; + +Segment tangent(const Circle& c1, const Circle& c2); + +std::pair, Point> +tangentPoints(const Circle& c1, const Circle& c2); + +std::pair +tangentPoints(const RationalRadiusCircle& c1, const RationalRadiusCircle& c2); + +CSTraits::Point_2 closestOnCircle(const Circle& circle, const Point& point); + +std::variant, std::pair, Segment>> +algebraicCircleTangentToRationalSegments(const CSTraits::Point_2& p1, const CSTraits::Point_2& p2, + const Circle& c1, const Circle& c2); + +CSPolygon approximateConvexHull(const std::vector>& circles); +} + +#endif //CARTOCROW_APPROXIMATE_CONVEX_HULL_H diff --git a/cartocrow/simplesets/helpers/arrangement_helpers.h b/cartocrow/simplesets/helpers/arrangement_helpers.h new file mode 100644 index 00000000..9a17743c --- /dev/null +++ b/cartocrow/simplesets/helpers/arrangement_helpers.h @@ -0,0 +1,26 @@ +#ifndef CARTOCROW_ARRANGEMENT_HELPERS_H +#define CARTOCROW_ARRANGEMENT_HELPERS_H + +#include +#include + +template +CGAL::General_polygon_2 ccb_to_polygon(Ccb ccb) { + Traits traits; + auto opposite = traits.construct_opposite_2_object(); + auto curr = ccb; + + std::vector x_monotone_curves; + do { + auto curve = curr->curve(); + if (curr->source()->point() == curve.source()) { + x_monotone_curves.push_back(curve); + } else { + x_monotone_curves.push_back(opposite(curve)); + } + } while(++curr != ccb); + + return {x_monotone_curves.begin(), x_monotone_curves.end()}; +} + +#endif //CARTOCROW_ARRANGEMENT_HELPERS_H diff --git a/cartocrow/simplesets/helpers/cropped_voronoi.h b/cartocrow/simplesets/helpers/cropped_voronoi.h new file mode 100644 index 00000000..a2c130cb --- /dev/null +++ b/cartocrow/simplesets/helpers/cropped_voronoi.h @@ -0,0 +1,127 @@ +#ifndef CARTOCROW_CROPPED_VORONOI_H +#define CARTOCROW_CROPPED_VORONOI_H + +#include "cartocrow/core/core.h" +#include "cartocrow/simplesets/types.h" +#include + +using namespace cartocrow; +using namespace cartocrow::simplesets; + +std::optional> intersectionConvex(const Polygon& polygon, const Ray& ray) { + auto sourceInside = polygon.has_on_bounded_side(ray.source()); + + std::vector> inters; + for (auto eit = polygon.edges_begin(); eit != polygon.edges_end(); eit++) { + auto edge = *eit; + auto obj = CGAL::intersection(ray, edge); + if (obj.has_value()) { + Segment seg; + Point point; + if (CGAL::assign(seg, obj)) { + return seg; + } else if (CGAL::assign(point, obj)) { + inters.push_back(point); + } + } + } + + assert(!(!sourceInside && inters.size() == 1)); + if (inters.empty()) { + return std::nullopt; + } else if (inters.size() == 2) { + return Segment(inters[0], inters[1]); + } else { + return Segment(ray.source(), inters[0]); + } +} + +std::optional> intersectionConvex(const Polygon& polygon, const Line& line) { + std::vector> inters; + for (auto eit = polygon.edges_begin(); eit != polygon.edges_end(); eit++) { + auto edge = *eit; + auto obj = CGAL::intersection(line, edge); + if (obj.has_value()) { + Segment seg; + Point point; + if (CGAL::assign(seg, obj)) { + return seg; + } else if (CGAL::assign(point, obj)) { + inters.push_back(point); + } + } + } + + assert(inters.size() != 1); + if (inters.size() == 2) { + return Segment(inters[0], inters[1]); + } else { + return std::nullopt; + } +} + +std::optional> intersectionConvex(const Polygon& polygon, const Segment& segment) { + auto sourceInside = polygon.has_on_bounded_side(segment.source()); + auto targetInside = polygon.has_on_bounded_side(segment.target()); + if (sourceInside && targetInside) { + return segment; + } + + std::vector> inters; + for (auto eit = polygon.edges_begin(); eit != polygon.edges_end(); eit++) { + auto edge = *eit; + auto obj = CGAL::intersection(segment, edge); + if (obj.has_value()) { + Segment seg; + Point point; + if (CGAL::assign(seg, obj)) { + return seg; + } else if (CGAL::assign(point, obj)) { + inters.push_back(point); + } else { + throw std::runtime_error("Intersection between two line segments is not a point nor a line segment"); + } + } + } + + if (!sourceInside && !targetInside) { + if (inters.empty()) return std::nullopt; + assert(inters.size() == 2); + // note that here, the orientation of the result segment may not be the same as the original segment. + return Segment(inters[0], inters[1]); + } + + if (sourceInside) { + assert(!targetInside); + assert(inters.size() == 1); + return Segment(segment.source(), inters[0]); + } + + assert(targetInside); + assert(inters.size() == 1); + return Segment(inters[0], segment.target()); +} + +// This part is adapted from: +// https://github.com/CGAL/cgal/blob/master/Triangulation_2/examples/Triangulation_2/print_cropped_voronoi.cpp +// which falls under the CC0 license. + +//A class to recover Voronoi diagram from stream. +//Rays, lines and segments are cropped to a polygon +struct Cropped_voronoi_from_delaunay{ + std::list, Segment>> m_cropped_vd; + Polygon m_clipper; + Cropped_voronoi_from_delaunay(const Polygon& clipper):m_clipper(clipper){} + template + void crop_and_extract_segment(const Point& site, const RSL& rsl){ + auto s = intersectionConvex(m_clipper, rsl); + if (s.has_value()) { + m_cropped_vd.push_back({site, *s}); + } + } + void operator<<(std::pair&, const Ray&> site_ray) { crop_and_extract_segment(site_ray.first, site_ray.second); } + void operator<<(std::pair&, const Line&> site_line) { crop_and_extract_segment(site_line.first, site_line.second); } + void operator<<(std::pair&, const Segment&> site_seg){ crop_and_extract_segment(site_seg.first, site_seg.second); } +}; + +#endif //CARTOCROW_CROPPED_VORONOI_H diff --git a/cartocrow/simplesets/helpers/cs_curve_helpers.cpp b/cartocrow/simplesets/helpers/cs_curve_helpers.cpp new file mode 100644 index 00000000..b26a4b78 --- /dev/null +++ b/cartocrow/simplesets/helpers/cs_curve_helpers.cpp @@ -0,0 +1,194 @@ +#include "cs_curve_helpers.h" + +namespace cartocrow::simplesets { +OneRootPoint closestOnCircle(const Circle& circle, const Point& point) { + auto bb = circle.bbox(); + Rectangle bbX(bb.xmin() - 1, bb.ymin() - 1, bb.xmax() + 1, bb.ymax() + 1); + Ray ray(circle.center(), point); + auto inter = CGAL::intersection(bbX, ray); + auto seg = get>(*inter); + + using Arr = CGAL::Arrangement_2; + Arr arr; + CGAL::insert(arr, circle); + CGAL::insert(arr, seg); + + std::optional intersectionPoint; + for (auto vit = arr.vertices_begin(); vit != arr.vertices_end(); ++vit) { + if (vit->degree() == 4) { + intersectionPoint = vit->point(); + } + } + + if (!intersectionPoint.has_value()) { + std::stringstream ss; + ss << "Could not project point " << point << " on circle " << circle; + throw std::runtime_error(ss.str()); + } + + return *intersectionPoint; +} + +OneRootPoint nearest(const X_monotone_curve_2& xm_curve, const Point& point) { + if (xm_curve.is_linear()) { + auto l = xm_curve.supporting_line(); + auto k = l.perpendicular(point); + auto inter = get>((*CGAL::intersection(l, k))); + + if (xm_curve.is_vertical()) { + auto minY = std::min(xm_curve.left().y(), xm_curve.right().y()); + auto maxY = std::max(xm_curve.left().y(), xm_curve.right().y()); + auto x = xm_curve.left().x(); + if (inter.y() >= maxY) { + return {x, maxY}; + } else if (inter.y() <= minY) { + return {x, minY}; + } else { + return {inter.x(), inter.y()}; + } + } + + if (inter.x() <= xm_curve.left().x()) { + return xm_curve.left(); + } else if (inter.x() >= xm_curve.right().x()) { + return xm_curve.right(); + } else { + return {inter.x(), inter.y()}; + } + } else { + auto c = xm_curve.supporting_circle(); + xm_curve.point_position({point.x(), point.y()}); + auto inter = closestOnCircle(c, point); + if (inter.x() <= xm_curve.left().x()) { + return xm_curve.left(); + } else if (inter.x() >= xm_curve.right().x()) { + return xm_curve.right(); + } else { + if (liesOn(inter, xm_curve)) { + return inter; + } + OneRootPoint opposite(2 * c.center().x() - inter.x(), 2 * c.center().y() - inter.y()); + if (liesOn(opposite, xm_curve)) { + return opposite; + } + auto sdLeft = CGAL::square(point.x() - xm_curve.left().x()) + CGAL::square(point.y() - xm_curve.left().y()); + auto sdRight = CGAL::square(point.x() - xm_curve.right().x()) + CGAL::square(point.y() - xm_curve.right().y()); + if (sdLeft < sdRight) { + return xm_curve.left(); + } else { + return xm_curve.right(); + } + } + } +} + +bool liesOn(const Point& p, const X_monotone_curve_2& xm_curve) { + if (p.x() < xm_curve.left().x() || p.x() > xm_curve.right().x()) + return false; + + if (xm_curve.is_linear()) { + return xm_curve.supporting_line().has_on(p); + } else { + return xm_curve.point_position({p.x(), p.y()}) == CGAL::EQUAL; + } +} + +bool liesOn(const OneRootPoint& p, const X_monotone_curve_2& xm_curve) { + if (p.x() < xm_curve.left().x() || p.x() > xm_curve.right().x()) + return false; +// CSTraits traits; +// auto lies_on = traits.compare_y_at_x_2_object(); + return xm_curve.point_position(p) == CGAL::EQUAL; +} + +bool liesOn(const X_monotone_curve_2& c1, const X_monotone_curve_2& c2) { + auto sc = liesOn(c1.source(), c2); + auto tc = liesOn(c1.target(), c2); + if (!sc || !tc) { + return false; + } + if (c2.is_linear()) { + if (c1.is_circular()) return false; + if (c1.supporting_line() != c2.supporting_line()) return false; + } else { + if (c1.is_linear()) return false; + if (c1.supporting_circle() != c2.supporting_circle()) return false; + } + + return true; +} + +renderer::RenderPath renderPath(const X_monotone_curve_2& xm_curve) { + renderer::RenderPath path; + path.moveTo(approximateAlgebraic(xm_curve.source())); + + if (xm_curve.is_circular()) { + auto circle = xm_curve.supporting_circle(); + path.arcTo(approximate(circle.center()), xm_curve.orientation() == CGAL::CLOCKWISE, + approximateAlgebraic(xm_curve.target())); + } else { + auto line = xm_curve.supporting_line(); + path.lineTo(approximateAlgebraic(xm_curve.target())); + } + + return path; +} + +void addToRenderPath(const X_monotone_curve_2& xm_curve, renderer::RenderPath& path, bool& first) { + auto as = approximateAlgebraic(xm_curve.source()); + auto at = approximateAlgebraic(xm_curve.target()); + if (first) { + path.moveTo(as); + first = false; + } + if (xm_curve.is_linear()) { + path.lineTo(at); + } else if (xm_curve.is_circular()){ + if (CGAL::squared_distance(as, at) < M_EPSILON) { + return; + } + const auto& circle = xm_curve.supporting_circle(); + path.arcTo(approximate(circle.center()), xm_curve.orientation() == CGAL::CLOCKWISE, approximateAlgebraic(xm_curve.target())); + } +} + +void addToRenderPath(const Curve_2& curve, renderer::RenderPath& path, bool& first) { + if (curve.is_full()) { + const auto& circ = curve.supporting_circle(); + const auto& cent = approximate(circ.center()); + auto r = sqrt(CGAL::to_double(circ.squared_radius())); + Point s(cent.x() - r, cent.y()); + auto clockwise = circ.orientation() == CGAL::CLOCKWISE; + path.moveTo(s); + path.arcTo(cent, clockwise, {cent.x() + r, cent.y()}); + path.arcTo(cent, clockwise, s); + path.close(); + return; + } + auto as = approximateAlgebraic(curve.source()); + auto at = approximateAlgebraic(curve.target()); + if (first) { + path.moveTo(as); + first = false; + } + if (curve.is_linear()) { + path.lineTo(at); + } else if (curve.is_circular()){ + if (CGAL::squared_distance(as, at) < M_EPSILON) { + return; + } + const auto& circle = curve.supporting_circle(); + path.arcTo(approximate(circle.center()), curve.orientation() == CGAL::CLOCKWISE, approximateAlgebraic(curve.target())); + } +} + +Curve_2 toCurve(const X_monotone_curve_2& xmc) { + if (xmc.is_linear()) { + return {xmc.supporting_line(), xmc.source(), xmc.target()}; + } else if (xmc.is_circular()) { + return {xmc.supporting_circle(), xmc.source(), xmc.target()}; + } else { + throw std::runtime_error("Impossible: circle-segment x-monotone curve is neither linear nor circular."); + } +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/helpers/cs_curve_helpers.h b/cartocrow/simplesets/helpers/cs_curve_helpers.h new file mode 100644 index 00000000..a2d14c67 --- /dev/null +++ b/cartocrow/simplesets/helpers/cs_curve_helpers.h @@ -0,0 +1,87 @@ +#include "../types.h" +#include "cartocrow/renderer/render_path.h" + +#ifndef CARTOCROW_CS_CURVE_HELPERS_H +#define CARTOCROW_CS_CURVE_HELPERS_H + +namespace cartocrow::simplesets { +template +void curveToXMonotoneCurves(const Curve_2& curve, OutputIterator out) { + CSTraits traits; + auto make_x_monotone = traits.make_x_monotone_2_object(); + std::vector> curves_and_points; + make_x_monotone(curve, std::back_inserter(curves_and_points)); + + // There should not be any isolated points + for (auto kinda_curve : curves_and_points) { + if (kinda_curve.which() == 1) { + *out++ = boost::get(kinda_curve); + } else { + std::cout << "Converting curve into x-monotone curves results in isolated point." + << std::endl; + throw std::runtime_error("Cannot convert a degenerate curve into x-monotone curves."); + } + } +} + +template +void curvesToXMonotoneCurves(InputIterator begin, InputIterator end, OutputIterator out) { + for (auto it = begin; it != end; ++it) { + curveToXMonotoneCurves(*it, out); + } +} + +OneRootPoint nearest(const X_monotone_curve_2& xm_curve, const Point& point); + +bool liesOn(const Point& p, const X_monotone_curve_2& xm_curve); +bool liesOn(const OneRootPoint& p, const X_monotone_curve_2& xm_curve); +renderer::RenderPath renderPath(const X_monotone_curve_2& xm_curve); +void addToRenderPath(const X_monotone_curve_2& xm_curve, renderer::RenderPath& path, bool& first); +void addToRenderPath(const Curve_2& curve, renderer::RenderPath& path, bool& first); +Curve_2 toCurve(const X_monotone_curve_2& xmc); + +template +void toCurves(InputIterator begin, InputIterator end, OutputIterator out) { + std::optional lastCurve; + for (auto curr = begin; curr != end; ++curr) { + X_monotone_curve_2 xmc = *curr; + if (!lastCurve.has_value()) { + lastCurve = toCurve(xmc); + } else { + if (lastCurve->is_linear() && xmc.is_linear() && lastCurve->supporting_line() == xmc.supporting_line()) { + Curve_2 newCurve(lastCurve->supporting_line(), lastCurve->source(), xmc.target()); + lastCurve = newCurve; + } else if (lastCurve->is_circular() && xmc.is_circular() && lastCurve->supporting_circle() == xmc.supporting_circle()) { + Curve_2 newCurve; + if (xmc.target() == lastCurve->source()) { + newCurve = Curve_2(lastCurve->supporting_circle()); + } else { + newCurve = Curve_2(lastCurve->supporting_circle(), lastCurve->source(), xmc.target()); + } + lastCurve = newCurve; + } else { + ++out = *lastCurve; + lastCurve = toCurve(xmc); + } + } + } + if (lastCurve.has_value()) { + ++out = *lastCurve; + } +} + +template +CSPolycurve arrPolycurveFromXMCurves(InputIterator begin, InputIterator end) { + PolyCSTraits traits; + auto construct = traits.construct_curve_2_object(); + std::vector curves; + std::transform(begin, end, std::back_inserter(curves), [](const X_monotone_curve_2& xm_curve) { + return toCurve(xm_curve); + }); + return construct(curves.begin(), curves.end()); +} + +bool liesOn(const X_monotone_curve_2& c1, const X_monotone_curve_2& c2); +} + +#endif //CARTOCROW_CS_CURVE_HELPERS_H diff --git a/cartocrow/simplesets/helpers/cs_polygon_helpers.cpp b/cartocrow/simplesets/helpers/cs_polygon_helpers.cpp new file mode 100644 index 00000000..46dc4262 --- /dev/null +++ b/cartocrow/simplesets/helpers/cs_polygon_helpers.cpp @@ -0,0 +1,253 @@ +#include "cs_polygon_helpers.h" +#include "cs_curve_helpers.h" +#include + +// Code adapted from a Stack Overflow answer by HEKTO. +// Link: https://stackoverflow.com/questions/69399922/how-does-one-obtain-the-area-of-a-general-polygon-set-in-cgal +// License info: https://stackoverflow.com/help/licensing +// The only changes made were changing the auto return type to Number and using the +// typedefs for circle, point, polygon etc. + +namespace cartocrow::simplesets { +//For two circles of radii R and r and centered at (0,0) and (d,0) intersecting +//in a region shaped like an asymmetric lens. +constexpr double lens_area(const double r, const double R, const double d) { + return r * r * std::acos((d * d + r * r - R * R) / 2 / d / r) + + R * R * std::acos((d * d + R * R - r * r) / 2 / d / R) - + 0.5 * std::sqrt((-d + r + R) * (d + r - R) * (d - r + R) * (d + r + R)); +} + +// ------ return signed area under the linear segment (P1, P2) +Number area(const CSTraits::Point_2& P1, const CSTraits::Point_2& P2) { + auto const dx = CGAL::to_double(P1.x()) - CGAL::to_double(P2.x()); + auto const sy = CGAL::to_double(P1.y()) + CGAL::to_double(P2.y()); + return dx * sy / 2; +} + +// ------ return signed area under the circular segment (P1, P2, C) +Number area(const CSTraits::Point_2& P1, const CSTraits::Point_2& P2, const CSTraits::Rational_circle_2& C) { + auto const dx = CGAL::to_double(P1.x()) - CGAL::to_double(P2.x()); + auto const dy = CGAL::to_double(P1.y()) - CGAL::to_double(P2.y()); + auto const squaredChord = dx * dx + dy * dy; + auto const chord = std::sqrt(squaredChord); + auto const squaredRadius = CGAL::to_double(C.squared_radius()); + auto const areaSector = + squaredRadius * std::asin(std::min(1.0, chord / (std::sqrt(squaredRadius) * 2))); + auto const areaTriangle = chord * std::sqrt(std::max(0.0, squaredRadius * 4 - squaredChord)) / 4; + auto const areaCircularSegment = areaSector - areaTriangle; + int sign; + if (C.orientation() == CGAL::Sign::NEGATIVE) { + sign = -1; + } else if (C.orientation() == CGAL::Sign::POSITIVE) { + sign = 1; + } else { + assert(C.orientation() == CGAL::Sign::ZERO); + sign = 0; + } + return area(P1, P2) + sign * areaCircularSegment; +} + +// ------ return signed area under the X-monotone curve +Number area(const X_monotone_curve_2& XCV) { + if (XCV.is_linear()) { + return area(XCV.source(), XCV.target()); + } else if (XCV.is_circular()) { + return area(XCV.source(), XCV.target(), XCV.supporting_circle()); + } else { + return 0; + } +} + +// ------ return area of the simple polygon +Number area(const CSPolygon& P) { + Number res = 0; + for (auto it = P.curves_begin(); it != P.curves_end(); ++it) { + res += area(*it); + } + return res; +} + +// ------ return area of the polygon with (optional) holes +Number area(const CSPolygonWithHoles& P) { + auto res = area(P.outer_boundary()); + for (auto it = P.holes_begin(); it != P.holes_end(); ++it) { + res += area(*it); + } + return res; +} + +CSPolygon circleToCSPolygon(const Circle& circle) { + std::vector xm_curves; + curveToXMonotoneCurves(circle, std::back_inserter(xm_curves)); + return {xm_curves.begin(), xm_curves.end()}; +} + +std::optional liesOn(const Point& p, const CSPolygon& polygon) { + for (auto cit = polygon.curves_begin(); cit != polygon.curves_end(); ++cit) { + if (liesOn(p, *cit)) { + return cit; + } + } + return std::nullopt; +} + +std::optional liesOn(const OneRootPoint& p, const CSPolygon& polygon) { + for (auto cit = polygon.curves_begin(); cit != polygon.curves_end(); ++cit) { + if (liesOn(p, *cit)) { + return cit; + } + } + return std::nullopt; +} + +renderer::RenderPath operator<<(renderer::RenderPath& path, const CSPolygon& polygon) { + bool first = true; + std::vector mergedCurves; + toCurves(polygon.curves_begin(), polygon.curves_end(), std::back_inserter(mergedCurves)); + for (const auto& c : mergedCurves) { + addToRenderPath(c, path, first); + } + if (!holds_alternative(path.commands().back())) { + path.close(); + } + return path; +} + +renderer::RenderPath renderPath(const CSPolygon& polygon) { + renderer::RenderPath path; + path << polygon; + return path; +} + +renderer::RenderPath renderPath(const CSPolygonWithHoles& withHoles) { + renderer::RenderPath path; + path << withHoles.outer_boundary(); + for (auto hit = withHoles.holes_begin(); hit != withHoles.holes_end(); ++hit) { + path << *hit; + } + return path; +} + +bool on_or_inside(const CSPolygon& polygon, const Point& point) { + Ray ray(point, Vector(1, 0)); + + Rectangle bbox = polygon.bbox(); + Rectangle rect({bbox.xmin() - 1, bbox.ymin() - 1}, {bbox.xmax() + 1, bbox.ymax() + 1}); + + auto inter = CGAL::intersection(ray, rect); + if (!inter.has_value()) return false; + if (inter->type() == typeid(Point)) return true; + auto seg = boost::get>(*inter); + X_monotone_curve_2 seg_xm_curve(seg.source(), seg.target()); + + typedef std::pair Intersection_point; + typedef boost::variant Intersection_result; + std::vector intersection_results; + + for (auto cit = polygon.curves_begin(); cit != polygon.curves_end(); ++cit) { + const auto& curve = *cit; + curve.intersect(seg_xm_curve, std::back_inserter(intersection_results)); + } + + int count = 0; + for (const auto& ir : intersection_results) { + if (ir.which() == 0) { + auto ip = get(ir); + //(ir) Intersection points are double-counted, so increase count by half. + bool found = false; + for (auto cit = polygon.curves_begin(); cit != polygon.curves_end(); ++cit) { + if (cit->source() == ip.first) { + found = true; + break; + } + } + if (found) { + count += 1; + continue; + } + } + count += 2; + } + + return count % 4 != 0; +} + +bool liesOn(const X_monotone_curve_2& c, const CSPolygon& polygon) { + auto sc = liesOn(c.source(), polygon); + auto tc = liesOn(c.target(), polygon); + if (!sc.has_value() || !tc.has_value()) { + return false; + } + do { + auto next = *sc; + ++next; + if (next == polygon.curves_end()) { + next = polygon.curves_begin(); + } + if (liesOn(c.source(), *next)) { + sc = next; + } else { + break; + } + } while (true); + do { + auto next = *tc; + ++next; + if (next == polygon.curves_end()) { + next = polygon.curves_begin(); + } + if (liesOn(c.source(), *next)) { + tc = next; + } else { + break; + } + } while (true); + auto sit = *sc; + auto tit = *sc; + auto curr = sit; + do { + if (curr->is_linear()) { + if (c.is_circular()) return false; + if (curr->supporting_line() != c.supporting_line()) return false; + } else { + if (c.is_linear()) return false; + if (curr->supporting_circle() != c.supporting_circle()) return false; + } + } while (curr++ != tit); + + return true; +} + +bool inside(const CSPolygon& polygon, const Point& point) { + return on_or_inside(polygon, point) && !liesOn(point, polygon); +} + +CSPolycurve arrPolycurveFromCSPolygon(const CSPolygon& polygon) { + return arrPolycurveFromXMCurves(polygon.curves_begin(), polygon.curves_end()); +} + +Polygon linearSample(const CSPolygon& polygon, int n) { + std::vector> coords; + polygon.approximate(std::back_inserter(coords), n); + std::vector> points; + + // polygon.approximate returns duplicate points. + for (int i = 0; i < coords.size(); ++i) { + auto p = coords[i]; + if (i > 0) { + auto prev = coords[i - 1]; + if (p == prev) continue; + } + if (i == coords.size() - 1 && p == coords[0]) { + continue; + } + points.emplace_back(p.first, p.second); + } + return {points.begin(), points.end()}; +} + +CSPolygonWithHoles approximateDilate(const CSPolygon& polygon, double r, double eps, int n) { + auto poly = linearSample(polygon, n); + return CGAL::approximated_offset_2(poly, r, eps); +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/helpers/cs_polygon_helpers.h b/cartocrow/simplesets/helpers/cs_polygon_helpers.h new file mode 100644 index 00000000..0224ae6b --- /dev/null +++ b/cartocrow/simplesets/helpers/cs_polygon_helpers.h @@ -0,0 +1,41 @@ +#ifndef CARTOCROW_CS_POLYGON_HELPERS_H +#define CARTOCROW_CS_POLYGON_HELPERS_H + +#include "../types.h" +#include "../../renderer/render_path.h" + +namespace cartocrow::simplesets { +//For two circles of radii R and r and centered at (0,0) and (d,0) intersecting +//in a region shaped like an asymmetric lens. +constexpr double lens_area(const double r, const double R, const double d); + +// ------ return signed area under the linear segment (P1, P2) +Number area(const CSTraits::Point_2& P1, const CSTraits::Point_2& P2); + +// ------ return signed area under the circular segment (P1, P2, C) +Number area(const CSTraits::Point_2& P1, const CSTraits::Point_2& P2, const CSTraits::Rational_circle_2& C); + +// ------ return signed area under the X-monotone curve +Number area(const X_monotone_curve_2& XCV); + +// ------ return area of the simple polygon +Number area(const CSPolygon& P); + +// ------ return area of the polygon with (optional) holes +Number area(const CSPolygonWithHoles& P); + +CSPolygon circleToCSPolygon(const Circle& circle); + +std::optional liesOn(const Point& p, const CSPolygon& polygon); +std::optional liesOn(const OneRootPoint& p, const CSPolygon& polygon); +bool liesOn(const X_monotone_curve_2& c, const CSPolygon& polygon); +renderer::RenderPath renderPath(const CSPolygon& polygon); +renderer::RenderPath renderPath(const CSPolygonWithHoles& withHoles); +bool on_or_inside(const CSPolygon& polygon, const Point& point); +bool inside(const CSPolygon& polygon, const Point& point); +CSPolycurve arrPolycurveFromCSPolygon(const CSPolygon& polygon); +Polygon linearSample(const CSPolygon& polygon, int n); +CSPolygonWithHoles approximateDilate(const CSPolygon& polygon, double r, double eps, int n); +} + +#endif //CARTOCROW_CS_POLYGON_HELPERS_H diff --git a/cartocrow/simplesets/helpers/cs_polyline_helpers.cpp b/cartocrow/simplesets/helpers/cs_polyline_helpers.cpp new file mode 100644 index 00000000..444a9562 --- /dev/null +++ b/cartocrow/simplesets/helpers/cs_polyline_helpers.cpp @@ -0,0 +1,85 @@ +#include "cs_polyline_helpers.h" +#include "cs_curve_helpers.h" + +namespace cartocrow::simplesets { +OneRootPoint nearest(const CSPolyline& polyline, const Point& point) { + std::optional minSqrdDist; + std::optional closest; + for (auto cit = polyline.curves_begin(); cit != polyline.curves_end(); ++cit) { + const auto& curve = *cit; + auto n = nearest(curve, point); + auto sqrdDist = CGAL::square(n.x() - point.x()) + CGAL::square(n.y() - point.y()); + if (!minSqrdDist.has_value() || sqrdDist < *minSqrdDist) { + minSqrdDist = sqrdDist; + closest = n; + } + } + + if (!closest.has_value()) { + throw std::runtime_error("Cannot find closest point to empty polyline."); + } + + return *closest; +} + +std::optional liesOn(const Point& p, const CSPolyline& polyline) { + for (auto cit = polyline.curves_begin(); cit != polyline.curves_end(); ++cit) { + if (liesOn(p, *cit)) { + return cit; + } + } + return std::nullopt; +} + +std::optional liesOn(const OneRootPoint& p, const CSPolyline& polyline) { + for (auto cit = polyline.curves_begin(); cit != polyline.curves_end(); ++cit) { + if (liesOn(p, *cit)) { + return cit; + } + } + return std::nullopt; +} + +renderer::RenderPath renderPath(const CSPolyline& polyline) { + renderer::RenderPath path; + bool first = true; + for (auto cit = polyline.curves_begin(); cit != polyline.curves_end(); ++cit) { + addToRenderPath(*cit, path, first); + } + return path; +} + +bool liesOn(const X_monotone_curve_2& c, const CSPolyline& polyline) { + auto sc = liesOn(c.source(), polyline); + auto tc = liesOn(c.target(), polyline); + if (!sc.has_value() || !tc.has_value()) { + return false; + } + auto sit = *sc; + auto tit = *sc; + auto curr = sit; + do { + if (curr->is_linear()) { + if (c.is_circular()) return false; + if (curr->supporting_line() != c.supporting_line()) return false; + } else { + if (c.is_linear()) return false; + if (curr->supporting_circle() != c.supporting_circle()) return false; + } + } while (curr++ != tit); + + return true; +} + +CSPolycurve arrPolycurveFromCSPolyline(const CSPolyline& polyline) { + return arrPolycurveFromXMCurves(polyline.curves_begin(), polyline.curves_end()); +} + +CSPolyline polylineToCSPolyline(const Polyline& polyline) { + std::vector xm_curves; + for (auto eit = polyline.edges_begin(); eit != polyline.edges_end(); ++eit) { + xm_curves.emplace_back(eit->source(), eit->target()); + } + return {xm_curves.begin(), xm_curves.end()}; +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/helpers/cs_polyline_helpers.h b/cartocrow/simplesets/helpers/cs_polyline_helpers.h new file mode 100644 index 00000000..a410ea14 --- /dev/null +++ b/cartocrow/simplesets/helpers/cs_polyline_helpers.h @@ -0,0 +1,18 @@ +#ifndef CARTOCROW_CS_POLYLINE_HELPERS_H +#define CARTOCROW_CS_POLYLINE_HELPERS_H + +#include "../types.h" +#include "cartocrow/renderer/render_path.h" +#include + +namespace cartocrow::simplesets { +OneRootPoint nearest(const CSPolyline& polyline, const Point& point); +std::optional liesOn(const Point& p, const CSPolyline& polyline); +std::optional liesOn(const OneRootPoint& p, const CSPolyline& polyline); +bool liesOn(const X_monotone_curve_2& c, const CSPolyline& polyline); +renderer::RenderPath renderPath(const CSPolyline& polyline); +CSPolycurve arrPolycurveFromCSPolyline(const CSPolyline& polyline); +CSPolyline polylineToCSPolyline(const Polyline& polyline); +} + +#endif //CARTOCROW_CS_POLYLINE_HELPERS_H diff --git a/cartocrow/simplesets/helpers/point_voronoi_helpers.h b/cartocrow/simplesets/helpers/point_voronoi_helpers.h new file mode 100644 index 00000000..7f0fcf71 --- /dev/null +++ b/cartocrow/simplesets/helpers/point_voronoi_helpers.h @@ -0,0 +1,55 @@ +#ifndef CARTOCROW_POINT_VORONOI_HELPERS_H +#define CARTOCROW_POINT_VORONOI_HELPERS_H + +#include "cartocrow/core/core.h" +#include +#include +#include +#include +#include + +using namespace cartocrow; +//typedef CGAL::Delaunay_triangulation_2 DT; +//typedef CGAL::Delaunay_triangulation_adaptation_traits_2
AT; +//typedef CGAL::Delaunay_triangulation_caching_degeneracy_removal_policy_2
AP; +//typedef CGAL::Voronoi_diagram_2 VD; + +template, + class AT = CGAL::Delaunay_triangulation_adaptation_traits_2
, + class AP = CGAL::Delaunay_triangulation_caching_degeneracy_removal_policy_2
, + class VD = CGAL::Voronoi_diagram_2> +Polygon face_to_polygon(VD vd, typename VD::Face& face, Rectangle bbox) { + std::vector> pts; + auto circ_start = face.ccb(); + auto circ = circ_start; + + do { + auto he = *circ; + // Cannot easily access the geometry of the halfedge, so go via dual + const DT& dt = vd.dual(); + auto dte = he.dual(); + auto objE = dt.dual(dte); // line segment, ray, or line + Segment seg; + Ray ray; + Line line; + Segment segI; + if (CGAL::assign(seg, objE)) { + auto objI = CGAL::intersection(seg, bbox); + segI = CGAL::object_cast>(objI); + } else if (CGAL::assign(ray, objE)) { + auto objI = CGAL::intersection(ray, bbox); + segI = CGAL::object_cast>(objI); + } else if (CGAL::assign(line, objE)) { + auto objI = CGAL::intersection(line, bbox); + segI = CGAL::object_cast>(objI); + } else { + throw std::runtime_error("Voronoi edge is neither a segment, ray, nor a line."); + } + pts.push_back(segI.source()); + } while (++circ != circ_start); + + return {pts.begin(), pts.end()}; +} + +#endif //CARTOCROW_POINT_VORONOI_HELPERS_H diff --git a/cartocrow/simplesets/helpers/poly_line_gon_intersection.cpp b/cartocrow/simplesets/helpers/poly_line_gon_intersection.cpp new file mode 100644 index 00000000..d8827527 --- /dev/null +++ b/cartocrow/simplesets/helpers/poly_line_gon_intersection.cpp @@ -0,0 +1,25 @@ +#include "poly_line_gon_intersection.h" + +namespace cartocrow::simplesets { +std::vector intersection(const CSPolyline& line, const CSPolygon& gon, bool keepOverlap) { + CSPolygonWithHoles withHoles(gon); + return intersection(line, withHoles, keepOverlap); +} + +std::vector difference(const CSPolyline& line, const CSPolygon& gon, bool keepOverlap) { + CSPolygonWithHoles withHoles(gon); + return difference(line, withHoles, keepOverlap); +} + +std::vector intersection(const CSPolyline& line, const CSPolygonWithHoles& gon, bool keepOverlap) { + std::vector polylines; + intersection(line, gon, std::back_inserter(polylines), false, keepOverlap); + return polylines; +} + +std::vector difference(const CSPolyline& line, const CSPolygonWithHoles& gon, bool keepOverlap) { + std::vector polylines; + intersection(line, gon, std::back_inserter(polylines), true, keepOverlap); + return polylines; +} +} diff --git a/cartocrow/simplesets/helpers/poly_line_gon_intersection.h b/cartocrow/simplesets/helpers/poly_line_gon_intersection.h new file mode 100644 index 00000000..87221a1b --- /dev/null +++ b/cartocrow/simplesets/helpers/poly_line_gon_intersection.h @@ -0,0 +1,318 @@ +#ifndef CARTOCROW_POLY_LINE_GON_INTERSECTION_H +#define CARTOCROW_POLY_LINE_GON_INTERSECTION_H + +#include "cartocrow/simplesets/types.h" +#include "cs_curve_helpers.h" +#include "cs_polyline_helpers.h" +#include "cs_polygon_helpers.h" +#include +#include +#include + +namespace cartocrow::simplesets { +// todo: edge case where edges of dilated patterns overlap, so a half-edge may have multiple origins. +std::vector intersection(const CSPolyline& line, const CSPolygon& gon, bool keepOverlap); +std::vector difference(const CSPolyline& line, const CSPolygon& gon, bool keepOverlap); +std::vector intersection(const CSPolyline& line, const CSPolygonWithHoles& gon, bool keepOverlap); +std::vector difference(const CSPolyline& line, const CSPolygonWithHoles& gon, bool keepOverlap); + +template +void intersection(const CSPolyline& line, const CSPolygonWithHoles& gon, OutputIterator out, bool difference, bool keepOverlap) { + CSTraits traits; + auto equals = traits.equal_2_object(); + auto opposite = traits.construct_opposite_2_object(); + enum Origin { + Polyline, + PolygonOuter, + PolygonHole, + }; + struct HalfEdgeData { + std::vector origins; + }; + using Arr = CGAL::Arrangement_2>; + using FH = typename Arr::Face_handle; + using HeH = typename Arr::Halfedge_handle; + using VH = typename Arr::Vertex_handle; + + Arr arr; + class Observer : public CGAL::Arr_observer { + public: + Observer(Arr& arr) : CGAL::Arr_observer(arr) {} + + virtual void after_create_edge(HeH e) { + e->data().origins.push_back(currentOrigin); + e->twin()->data().origins.push_back(currentOrigin); + } + + virtual void before_split_edge(HeH e, VH v, const X_monotone_curve_2& c1, const X_monotone_curve_2& c2) { + beforeSplitData = e->data(); + } + + virtual void after_split_edge(HeH e1, HeH e2) { + e1->set_data(beforeSplitData); + e1->twin()->set_data(beforeSplitData); + e2->twin()->set_data(beforeSplitData); + e2->set_data(beforeSplitData); + CSTraits traits; + auto opposite = traits.construct_opposite_2_object(); + if (liesOn(e1->curve(), xmCurve) || liesOn(opposite(e1->curve()), xmCurve)) { + e1->data().origins.push_back(currentOrigin); + e1->twin()->data().origins.push_back(currentOrigin); + } + if (liesOn(e2->curve(), xmCurve) || liesOn(opposite(e2->curve()), xmCurve)) { + e2->data().origins.push_back(currentOrigin); + e2->twin()->data().origins.push_back(currentOrigin); + } + } + + virtual void after_modify_edge(HeH e) { + e->data().origins.push_back(currentOrigin); + e->twin()->data().origins.push_back(currentOrigin); + } + + Origin currentOrigin; + X_monotone_curve_2 xmCurve; + private: + HalfEdgeData beforeSplitData; + }; + + Observer observer(arr); + + observer.currentOrigin = Origin::Polyline; + for (auto cit = line.curves_begin(); cit != line.curves_end(); ++cit) { + observer.xmCurve = *cit; + CGAL::insert(arr, *cit); + } + observer.currentOrigin = Origin::PolygonOuter; + for (auto cit = gon.outer_boundary().curves_begin(); cit != gon.outer_boundary().curves_end(); ++cit) { + observer.xmCurve = *cit; + CGAL::insert(arr, *cit); + } + observer.currentOrigin = Origin::PolygonHole; + for (auto hit = gon.holes_begin(); hit != gon.holes_end(); ++hit) { + for (auto cit = hit->curves_begin(); cit != hit->curves_end(); ++cit) { + observer.xmCurve = *cit; + CGAL::insert(arr, *cit); + } + } + + std::vector line_edges_keep; + for (auto eit = arr.edges_begin(); eit != arr.edges_end(); ++eit) { + HeH edge = eit; + + bool onGonEdge = false; + bool onPolyline = false; + for (auto origin : edge->data().origins) { + if (origin == Origin::Polyline) { + onPolyline = true; + } + if (origin == Origin::PolygonOuter || origin == Origin::PolygonHole) { + onGonEdge = true; + } + } + if (!onPolyline) { + continue; + } + bool liesInGon = false; + if (onGonEdge) { + if (keepOverlap) { + liesInGon = true; + } + } else { + for (auto fh : {edge->face(), edge->twin()->face()}) { + if (!fh->has_outer_ccb()) { + continue; + } + auto ccb = fh->outer_ccb(); + auto ccbIt = ccb; + do { + if (!equals(ccbIt->source()->point(), ccbIt->curve().source())) + continue; + // if *ccbIt lies on outer face. + for (auto origin : ccbIt->data().origins) { + if (origin == Origin::PolygonOuter) { + liesInGon = true; + break; + } + } + } while (++ccbIt != ccb); + if (liesInGon) + break; + } + } + + if (!difference && liesInGon) { + if (equals(edge->source()->point(), edge->curve().source())) { + line_edges_keep.push_back(edge); + } else { + line_edges_keep.push_back(edge->twin()); + } + } + if (difference && !liesInGon) { + if (equals(edge->source()->point(), edge->curve().source())) { + line_edges_keep.push_back(edge); + } else { + line_edges_keep.push_back(edge->twin()); + } + } + } + + auto originatesFromPolyline = [&arr, &equals](HeH h) { + for (auto origin : h->data().origins) { + if (origin == Origin::Polyline && equals(h->source()->point(), h->curve().source())) { + return true; + } + } + return false; + }; + + while (!line_edges_keep.empty()) { + // Find first edge on connected component of polyline (in the intersection with polygon) + auto start = line_edges_keep.front(); + auto startStart = start; + while (originatesFromPolyline(startStart->prev())) { + startStart = startStart->prev(); + + if (startStart == start) { + break; + } + } + std::vector xmcs; + auto last_it = line_edges_keep.end(); + auto curr = startStart; + do { + last_it = std::remove(line_edges_keep.begin(), last_it, curr); + xmcs.push_back(curr->curve()); + curr = curr->next(); + } while (originatesFromPolyline(curr) && curr != startStart); + line_edges_keep.erase(last_it, line_edges_keep.end()); + + ++out = CSPolyline(xmcs.begin(), xmcs.end()); + } +} + +template +void intersection(const CSPolyline& line, const CSPolygon& gon, OutputIterator out, bool difference, bool keepOverlap) { + CSPolygonWithHoles withHoles(gon); + return intersection(line, withHoles, out, difference, keepOverlap); +} + +// May crash due to CGAL bug: https://github.com/CGAL/cgal/issues/8468 +#if 0 +template +void intersectionCrashes(const CSPolyline& line, const CSPolygonWithHoles& gon, OutputIterator out, bool difference, bool keepOverlap) { + PolyCSTraits traits; + auto equals = traits.equal_2_object(); + using Arr = CGAL::Arrangement_with_history_2; + Arr arr; + auto linePolycurve = arrPolycurveFromCSPolyline(line); + auto outerGonPolycurve = arrPolycurveFromCSPolygon(gon.outer_boundary()); + + auto ogch = CGAL::insert(arr, outerGonPolycurve); + auto lch = CGAL::insert(arr, linePolycurve); + std::vector hgchs; + + std::vector holesGonPolycurves; + for (auto hit = gon.holes_begin(); hit != gon.holes_end(); ++hit) { + hgchs.push_back(CGAL::insert(arr, arrPolycurveFromCSPolygon(*hit))); + } + + std::vector line_edges_keep; + for (auto eit = arr.induced_edges_begin(lch); eit != arr.induced_edges_end(lch); ++eit) { + Arr::Halfedge_handle edge = *eit; + + bool onGonEdge = false; + for (auto curveIt = arr.originating_curves_begin(edge); curveIt != arr.originating_curves_end(edge); ++curveIt) { + Arr::Curve_handle ch = curveIt; + if (ch == ogch) { + onGonEdge = true; + } + for (const auto& hgch : hgchs) { + if (ch == hgch) { + onGonEdge = true; + } + } + } + bool liesInGon = false; + if (onGonEdge) { + if (keepOverlap) { + liesInGon = true; + } + } else { + for (auto fh : {edge->face(), edge->twin()->face()}) { + if (!fh->has_outer_ccb()) { + continue; + } + auto ccb = fh->outer_ccb(); + auto ccbIt = ccb; + do { + if (!equals(ccbIt->source()->point(), ccbIt->curve().subcurves_begin()->source())) + continue; + // if *ccbIt lies on outer face. + for (auto curveIt = arr.originating_curves_begin(ccbIt); + curveIt != arr.originating_curves_end(ccbIt); ++curveIt) { + Arr::Curve_handle ch = curveIt; + if (ch == ogch) { + liesInGon = true; + break; + } + } + } while (++ccbIt != ccb); + if (liesInGon) + break; + } + } + + if (!difference && liesInGon) { + if (equals(edge->source()->point(), edge->curve().subcurves_begin()->source())) { + line_edges_keep.push_back(edge); + } else { + line_edges_keep.push_back(edge->twin()); + } + } + if (difference && !liesInGon) { + if (equals(edge->source()->point(), edge->curve().subcurves_begin()->source())) { + line_edges_keep.push_back(edge); + } else { + line_edges_keep.push_back(edge->twin()); + } + } + } + + auto originatesFromPolyline = [&arr, &lch, &equals](Arr::Halfedge_handle h) { + for (auto cit = arr.originating_curves_begin(h); cit != arr.originating_curves_end(h); ++cit) { + Arr::Curve_handle ch = cit; + if (ch == lch && equals(h->source()->point(), h->curve().subcurves_begin()->source())) { + return true; + } + } + return false; + }; + + while (!line_edges_keep.empty()) { + // Find first edge on connected component of polyline (in the intersection with polygon) + auto start = line_edges_keep.front(); + auto curr = start; + while (originatesFromPolyline(curr->prev())) { + curr = curr->prev(); + + if (curr == start) { + break; + } + } + std::vector xmcs; + auto last_it = line_edges_keep.end(); + do { + last_it = std::remove(line_edges_keep.begin(), last_it, curr); + std::copy(curr->curve().subcurves_begin(), curr->curve().subcurves_end(), std::back_inserter(xmcs)); + curr = curr->next(); + } while (originatesFromPolyline(curr)); + line_edges_keep.erase(last_it, line_edges_keep.end()); + + ++out = CSPolyline(xmcs.begin(), xmcs.end()); + } +} +#endif +} + +#endif //CARTOCROW_POLY_LINE_GON_INTERSECTION_H diff --git a/cartocrow/simplesets/parse_input.cpp b/cartocrow/simplesets/parse_input.cpp new file mode 100644 index 00000000..d3d6e767 --- /dev/null +++ b/cartocrow/simplesets/parse_input.cpp @@ -0,0 +1,44 @@ +#include "parse_input.h" + +namespace cartocrow::simplesets { +std::string getNextLine(std::istream& str) { + std::string line; + std::getline(str,line); + return line; +} + +std::vector splitIntoTokens(const std::string& line, char delimiter) { + std::vector result; + std::stringstream lineStream(line); + std::string cell; + + while(std::getline(lineStream, cell, delimiter)) + { + result.push_back(cell); + } + // This checks for a trailing comma with no data after it. + if (!lineStream && cell.empty()) + { + // If there was a trailing comma then add an empty element. + result.emplace_back(""); + } + return result; +} + +std::vector parseCatPoints(const std::string& s) { + std::stringstream ss(s); + + std::vector result; + + while (ss) { + auto parts = splitIntoTokens(getNextLine(ss), ' '); + if (parts.size() <= 1) break; + if (parts.size() != 3) { + throw std::runtime_error("Input has incorrect format."); + } + result.emplace_back(stoi(parts[0]), Point(stod(parts[1]), -stod(parts[2]))); + } + + return result; +} +} diff --git a/cartocrow/simplesets/parse_input.h b/cartocrow/simplesets/parse_input.h new file mode 100644 index 00000000..0831ee37 --- /dev/null +++ b/cartocrow/simplesets/parse_input.h @@ -0,0 +1,10 @@ +#ifndef CARTOCROW_PARSE_INPUT_H +#define CARTOCROW_PARSE_INPUT_H + +#include "cat_point.h" + +namespace cartocrow::simplesets { +std::vector parseCatPoints(const std::string& s); +} + +#endif //CARTOCROW_PARSE_INPUT_H diff --git a/cartocrow/simplesets/partition.h b/cartocrow/simplesets/partition.h new file mode 100644 index 00000000..6c563bdc --- /dev/null +++ b/cartocrow/simplesets/partition.h @@ -0,0 +1,10 @@ +#ifndef CARTOCROW_PARTITION_H +#define CARTOCROW_PARTITION_H + +#include "patterns/poly_pattern.h" + +namespace cartocrow::simplesets { +typedef std::vector> Partition; +} + +#endif //CARTOCROW_PARTITION_H diff --git a/cartocrow/simplesets/partition_algorithm.cpp b/cartocrow/simplesets/partition_algorithm.cpp new file mode 100644 index 00000000..245984bc --- /dev/null +++ b/cartocrow/simplesets/partition_algorithm.cpp @@ -0,0 +1,300 @@ +#include "partition_algorithm.h" +#include "dilated/dilated_poly.h" + +#include + +#include +#include "helpers/cs_polygon_helpers.h" + +namespace cartocrow::simplesets { + +std::variant to_bank_or_island(PolyPattern* polyPattern) { + if (auto bp = dynamic_cast(polyPattern)) { + return *bp; + } else if (auto ip = dynamic_cast(polyPattern)) { + return *ip; + } else if (auto spp = dynamic_cast(polyPattern)) { + return Bank(spp->catPoints()); + } else if (auto mp = dynamic_cast(polyPattern)) { + return Bank(mp->catPoints()); + } else { + throw std::runtime_error("Unknown PolyPattern."); + } +} + +Number squared_distance( + const std::variant, Polygon>& contour, + const Point& p) { + return std::visit([&p](auto& poly) { + Number min_sqrd_dist = std::numeric_limits::infinity(); + + for (auto eit = poly.edges_begin(); eit != poly.edges_end(); ++eit) { + auto seg = *eit; + auto dist = CGAL::squared_distance(seg, p); + if (dist < min_sqrd_dist) { + min_sqrd_dist = dist; + } + } + + return min_sqrd_dist; + }, contour); +} + +bool is_inside(const Point& point, const Polygon& polygon) { + return !polygon.has_on_unbounded_side(point); +} + +bool is_inside(const Point& point, const Polyline& polygon) { + return false; +} + +bool do_intersect( + const std::variant, Polygon>& cont1, + const std::variant, Polygon>& cont2) { + return std::visit([&cont2](auto& poly1) { + return std::visit([&poly1](auto& poly2) { + for (auto eit = poly1.edges_begin(); eit != poly1.edges_end(); ++eit) { + for (auto eitt = poly2.edges_begin(); eitt != poly2.edges_end(); ++eitt) { + if (CGAL::do_intersect(*eit, *eitt)) { + return true; + } + } + } + + return is_inside(poly1.vertex(0), poly2) || is_inside(poly2.vertex(0), poly1); + }, cont2); + }, cont1); +} + +Number intersectionDelay(const std::vector& points, const PolyPattern& p1, const PolyPattern& p2, + const PolyPattern& result, const GeneralSettings& gs, const PartitionSettings& ps) { + // todo: check consistency with paper + if (!ps.intersectionDelay) return 0; + Number intersectionArea = 0; + auto& resultPts = result.catPoints(); + auto resultPoly = result.poly(); + for (const auto& pt : points) { + if (std::find(resultPts.begin(), resultPts.end(), pt) == resultPts.end() && + squared_distance(resultPoly, pt.point) < squared(2 * gs.dilationRadius())) { + CSPolygon ptShape = Dilated(SinglePoint(pt), gs.dilationRadius()).m_contour; + CSPolygon rShape = Dilated(result, gs.dilationRadius()).m_contour; + std::vector inters; + CGAL::intersection(rShape, ptShape, std::back_inserter(inters)); + Number newArea = 0; + for (const auto& gp : inters) { + newArea += abs(area(gp)); + } + inters.clear(); + CSPolygon p1Shape = Dilated(p1, gs.dilationRadius()).m_contour; + CSPolygon p2Shape = Dilated(p2, gs.dilationRadius()).m_contour; + CGAL::intersection(p1Shape, ptShape, std::back_inserter(inters)); + CGAL::intersection(p2Shape, ptShape, std::back_inserter(inters)); + Number oldArea = 0; + for (const auto& gp : inters) { + oldArea += abs(area(gp)); + } + intersectionArea += newArea - oldArea; + } + } + + if (intersectionArea <= 0) { + return 0; + } + return sqrt(intersectionArea / M_PI); +} + +std::vector, Partition>> +partition(const std::vector& points, const GeneralSettings& gs, const PartitionSettings& ps, Number maxTime) { + // Create initial partition consisting of single points + std::vector> initialPatterns; + std::transform(points.begin(), points.end(), std::back_inserter(initialPatterns), [](const auto& pt) + { return std::make_shared(pt); }); + Partition partition(initialPatterns); + + std::vector, Partition>> history; + history.emplace_back(0, partition); + + // Priority queue storing events based on their time + auto comparison = [](const PossibleMergeEvent& e1, const PossibleMergeEvent& e2) { return e1.time > e2.time; }; + std::priority_queue, decltype(comparison)> events{comparison}; + + // Add SinglePoint--SinglePoint merges + for (int i = 0; i < partition.size(); ++i) { + auto& p = dynamic_cast(*partition[i]); + for (int j = i + 1; j < partition.size(); ++j) { + auto& q = dynamic_cast(*partition[j]); + if (p.category() != q.category()) continue; + if (squared_distance(p.catPoint().point, q.catPoint().point) > squared(2 * maxTime)) continue; + + auto newPattern = std::make_shared(p.catPoint(), q.catPoint()); + Segment seg(p.catPoint().point, q.catPoint().point); + + bool tooClose = false; + for (const CatPoint& pt : points) { + // Check if a point is too close to the segment + if (pt != p.catPoint() && pt != q.catPoint() && + CGAL::squared_distance(seg, pt.point) < squared(ps.admissibleRadiusFactor * gs.dilationRadius()) && + CGAL::squared_distance(seg, pt.point) < CGAL::min(CGAL::squared_distance(p.catPoint().point, pt.point), + CGAL::squared_distance(q.catPoint().point, pt.point)) - M_EPSILON) { + tooClose = true; + break; + } + } + if (tooClose) continue; + + PossibleMergeEvent event{newPattern->coverRadius(), partition[i], partition[j], newPattern, false}; + events.push(event); + } + } + + while (!events.empty()) { + auto ev = events.top(); + events.pop(); + + if (ev.time > maxTime) break; + + if (!ev.final) { + auto delay = intersectionDelay(points, *ev.p1, *ev.p2, *ev.result, gs, ps); + ev.time += delay; + ev.final = true; + events.push(ev); + continue; + } + + // Check if patterns that events wants to merge still exist + bool foundP1 = false; + bool foundP2 = false; + for (const auto& pattern : partition) { + if (pattern == ev.p1) { + foundP1 = true; + } + if (pattern == ev.p2) { + foundP2 = true; + } + } + if (!foundP1 || !foundP2) continue; + + auto& newPts = ev.result->catPoints(); + auto newPoly = ev.result->poly(); + + // Check for intersections with existing patterns + bool intersects = false; + for (const auto& pattern : partition) { + if (pattern != ev.p1 && pattern != ev.p2 && do_intersect(pattern->poly(), newPoly)) { + intersects = true; + break; + } + } + if (intersects) continue; + + bool tooClose = false; + // Check whether any point is too close to the result pattern + for (const auto& pt : points) { + if (std::find(newPts.begin(), newPts.end(), pt) == newPts.end()) { + auto polyPtDist = squared_distance(newPoly, pt.point); + std::optional> pointPtDist; + for (auto np : newPts) { + auto d = CGAL::squared_distance(np.point, pt.point); + if (!pointPtDist.has_value() || d < *pointPtDist) { + pointPtDist = d; + } + } + if (polyPtDist < squared(ps.admissibleRadiusFactor * gs.dilationRadius()) && + polyPtDist < pointPtDist) { + tooClose = true; + break; + } + } + } + if (tooClose) continue; + + // Merge + // Remove the two patterns that are merged + auto startTrash = std::remove_if(partition.begin(), partition.end(), + [&ev](const auto& part) { return part == ev.p1 || part == ev.p2; }); + partition.erase(startTrash, partition.end()); + // Add the result + partition.push_back(ev.result); + // Save this partition + history.emplace_back(ev.time, partition); + + // Create new merge events + for (const auto& pattern : partition) { + if (pattern == ev.result || pattern->category() != ev.result->category()) continue; + + if (ps.islands) { + // Do relatively cheap check first: are the points in the two patterns close enough to ever form a pattern? + Number min_sqrd_dist = std::numeric_limits::infinity(); + for (const auto& p : pattern->catPoints()) { + for (const auto& q : newPts) { + auto dist = squared_distance(p.point, q.point); + if (dist < min_sqrd_dist) { + min_sqrd_dist = dist; + } + } + } + + if (min_sqrd_dist <= squared(2 * maxTime)) { + std::vector mergedPoints; + std::copy(newPts.begin(), newPts.end(), std::back_inserter(mergedPoints)); + std::copy(pattern->catPoints().begin(), pattern->catPoints().end(), std::back_inserter(mergedPoints)); + auto newIsland = std::make_shared(mergedPoints); + + // todo? do intersection check here already? check how this affects performance + + Number regDelay = !ps.regularityDelay ? 0 : newIsland->coverRadius() - std::max(pattern->coverRadius(), ev.result->coverRadius()); + Number eventTime = newIsland->coverRadius() + regDelay; + PossibleMergeEvent newEvent{eventTime, ev.result, pattern, newIsland, false}; + if (eventTime <= maxTime) { + events.push(newEvent); + } + } + } + + if (ps.banks) { + auto pat1 = to_bank_or_island(&*pattern); + auto pat2 = to_bank_or_island(&*ev.result); + + if (std::holds_alternative(pat1) && std::holds_alternative(pat2)) { + Bank bank1 = std::get(pat1); + Bank bank2 = std::get(pat2); + auto pts1 = bank1.catPoints(); + auto pts2 = bank2.catPoints(); + + std::vector b1Pts; + std::copy(pts1.begin(), pts1.end(), std::back_inserter(b1Pts)); + std::copy(pts2.begin(), pts2.end(), std::back_inserter(b1Pts)); + auto b1 = std::make_shared(b1Pts); + + std::vector b2Pts; + std::copy(pts1.begin(), pts1.end(), std::back_inserter(b2Pts)); + std::copy(pts2.rbegin(), pts2.rend(), std::back_inserter(b2Pts)); + auto b2 = std::make_shared(b2Pts); + + std::vector b3Pts; + std::copy(pts1.rbegin(), pts1.rend(), std::back_inserter(b3Pts)); + std::copy(pts2.rbegin(), pts2.rend(), std::back_inserter(b3Pts)); + auto b3 = std::make_shared(b3Pts); + + std::vector b4Pts; + std::copy(pts1.rbegin(), pts1.rend(), std::back_inserter(b4Pts)); + std::copy(pts2.begin(), pts2.end(), std::back_inserter(b4Pts)); + auto b4 = std::make_shared(b4Pts); + + for (auto& b : {b1, b2, b3, b4}) { + if (!b->isValid(gs)) continue; + Number regDelay = !ps.regularityDelay ? 0 : b->coverRadius() - std::max(ev.result->coverRadius(), pattern->coverRadius()); + Number eventTime = b->coverRadius() + regDelay; + PossibleMergeEvent newEvent{eventTime * 1, ev.result, pattern, b, false}; + if (eventTime <= maxTime) { + events.push(newEvent); + } + } + } + } + } + } + + return history; +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/partition_algorithm.h b/cartocrow/simplesets/partition_algorithm.h new file mode 100644 index 00000000..c879ca3f --- /dev/null +++ b/cartocrow/simplesets/partition_algorithm.h @@ -0,0 +1,29 @@ +#ifndef CARTOCROW_PARTITION_ALGORITHM_H +#define CARTOCROW_PARTITION_ALGORITHM_H + +#include "types.h" +#include "cat_point.h" +#include "partition.h" +#include "patterns/bank.h" +#include "patterns/island.h" +#include "patterns/matching.h" +#include "patterns/single_point.h" +#include "settings.h" + +namespace cartocrow::simplesets { +struct PossibleMergeEvent { + Number time; + std::shared_ptr p1; + std::shared_ptr p2; + std::shared_ptr result; + /// We compute expensive delays lazily. This boolean indicates whether all delays have been computed. + bool final; +}; + +Number intersectionDelay(const std::vector& points, const PolyPattern& p1, const PolyPattern& p2, + const PolyPattern& result, const GeneralSettings& gs, const PartitionSettings& ps); + +std::vector, Partition>> +partition(const std::vector& points, const GeneralSettings& gs, const PartitionSettings& ps, Number maxTime); +} +#endif //CARTOCROW_PARTITION_ALGORITHM_H diff --git a/cartocrow/simplesets/partition_painting.cpp b/cartocrow/simplesets/partition_painting.cpp new file mode 100644 index 00000000..134120f8 --- /dev/null +++ b/cartocrow/simplesets/partition_painting.cpp @@ -0,0 +1,37 @@ +#include "partition_painting.h" + +namespace cartocrow::simplesets { + +void draw_poly_pattern(const PolyPattern& pattern, renderer::GeometryRenderer& renderer, const GeneralSettings& gs, const DrawSettings& ds) { + if (std::holds_alternative>(pattern.poly())) { + renderer.setMode(renderer::GeometryRenderer::stroke); + } else { + renderer.setMode(renderer::GeometryRenderer::fill | renderer::GeometryRenderer::stroke); + } + renderer.setFill(ds.getColor(pattern.category())); + renderer.setFillOpacity(100); + renderer.setStroke(Color{0, 0, 0}, ds.contourStrokeWeight(gs), true); + auto& pts = pattern.catPoints(); + if (pts.size() == 1) { +// renderer.draw(pts[0].point); + } else { + std::visit([&renderer](auto shape) { renderer.draw(shape); }, pattern.poly()); + } + for (const auto& pt : pattern.catPoints()) { + renderer.setStroke(Color{0, 0, 0}, ds.pointStrokeWeight(gs), true); + renderer.setFillOpacity(255); + renderer.setFill(ds.getColor(pt.category)); + renderer.setMode(renderer::GeometryRenderer::fill | renderer::GeometryRenderer::stroke); + renderer.draw(Circle{pt.point, gs.pointSize * gs.pointSize}); + } +} + +PartitionPainting::PartitionPainting(const Partition& partition, const GeneralSettings& gs, const DrawSettings& ds) + : m_partition(partition), m_gs(gs), m_ds(ds) {} + +void PartitionPainting::paint(renderer::GeometryRenderer& renderer) const { + for (const auto& pattern : m_partition) { + draw_poly_pattern(*pattern, renderer, m_gs, m_ds); + } +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/partition_painting.h b/cartocrow/simplesets/partition_painting.h new file mode 100644 index 00000000..141f1ced --- /dev/null +++ b/cartocrow/simplesets/partition_painting.h @@ -0,0 +1,24 @@ +#ifndef CARTOCROW_PARTITION_PAINTING_H +#define CARTOCROW_PARTITION_PAINTING_H + +#include "types.h" +#include "settings.h" +#include "partition.h" +#include "../renderer/geometry_painting.h" + +namespace cartocrow::simplesets { +class PartitionPainting : public renderer::GeometryPainting { + public: + PartitionPainting(const Partition& partition, const GeneralSettings& gs, const DrawSettings& ds); + void paint(renderer::GeometryRenderer& renderer) const override; + + private: + const Partition& m_partition; + const GeneralSettings& m_gs; + const DrawSettings& m_ds; +}; + +void draw_poly_pattern(const PolyPattern& pattern, renderer::GeometryRenderer& renderer, const GeneralSettings& gs, const DrawSettings& ds); +} + +#endif //CARTOCROW_PARTITION_PAINTING_H diff --git a/cartocrow/simplesets/patterns/bank.cpp b/cartocrow/simplesets/patterns/bank.cpp new file mode 100644 index 00000000..a1f9b65d --- /dev/null +++ b/cartocrow/simplesets/patterns/bank.cpp @@ -0,0 +1,90 @@ +#include "bank.h" + +namespace cartocrow::simplesets { +Bank::Bank(std::vector catPoints): m_catPoints(std::move(catPoints)) { + // Store the point positions separately, sometimes only the positions are needed. + std::transform(m_catPoints.begin(), m_catPoints.end(), std::back_inserter(m_points), [](const CatPoint& cp) { + return cp.point; + }); + + // Store polyline + m_polyline = Polyline(m_points); + + // Compute the cover radius + std::optional> maxSquaredDistance; + for (int i = 0; i < m_points.size() - 1; i++) { + auto p = m_points[i]; + auto q = m_points[i+1]; + auto dist = squared_distance(p, q); + if (!maxSquaredDistance.has_value() || dist > *maxSquaredDistance) { + maxSquaredDistance = dist; + } + } + m_coverRadius = sqrt(*maxSquaredDistance) / 2; + + // Compute bends + computeBends(); +} + +Number computeAngleBetween(const Vector& v, const Vector& w) { + return acos((v * w) / (sqrt(v.squared_length()) * sqrt(w.squared_length()))); +} + +void Bank::computeBends() { + m_bends.clear(); + + std::optional orientation; + Number bendTotalAngle = 0; + Number bendMaxAngle = 0; + int startIndex = 0; + + for (int i = 0; i < m_points.size(); i++) { + if (i + 2 > m_points.size() - 1) break; + auto orient = CGAL::orientation(m_points[i], m_points[i+1], m_points[i+2]); + auto angle = computeAngleBetween(m_points[i+1] - m_points[i], m_points[i+2] - m_points[i+1]); + if (orientation == -orient) { + // Switched orientation + m_bends.emplace_back(*orientation, bendMaxAngle, bendTotalAngle, startIndex, i+1); + orientation = orient; + bendTotalAngle = angle; + bendMaxAngle = angle; + startIndex = i; + } else { + orientation = orient; + bendTotalAngle += angle; + bendMaxAngle = std::max(bendMaxAngle, angle); + } + } + + if (orientation.has_value()) { + m_bends.emplace_back(*orientation, bendMaxAngle, bendTotalAngle, startIndex, static_cast(m_points.size()-1)); + } +} + +std::variant, Polygon> Bank::poly() const { + return m_polyline; +} + +const std::vector& Bank::catPoints() const { + return m_catPoints; +} + +bool Bank::isValid(const GeneralSettings& gs) const { + bool inflectionIsFine = m_bends.size() <= gs.inflectionLimit; + bool anglesAreFine = true; + for (const auto& bend : m_bends) { + anglesAreFine = anglesAreFine && bend.maxAngle <= gs.maxTurnAngle; + } + bool totalAngleIsFine = true; + for (const auto& bend : m_bends) { + if (bend.totalAngle > gs.maxBendAngle) { + totalAngleIsFine = false; + } + } + return inflectionIsFine && anglesAreFine && totalAngleIsFine; +} + +Number Bank::coverRadius() const { + return m_coverRadius; +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/patterns/bank.h b/cartocrow/simplesets/patterns/bank.h new file mode 100644 index 00000000..0acbf464 --- /dev/null +++ b/cartocrow/simplesets/patterns/bank.h @@ -0,0 +1,38 @@ +#ifndef CARTOCROW_BANK_H +#define CARTOCROW_BANK_H + +#include "pattern.h" +#include "poly_pattern.h" + +namespace cartocrow::simplesets { +struct Bend { + CGAL::Orientation orientation; + Number maxAngle; + Number totalAngle; + int startIndex; + int endIndex; + Bend(CGAL::Orientation orientation, Number maxAngle, Number totalAngle, int startIndex, int endIndex) : + orientation(orientation), maxAngle(maxAngle), totalAngle(totalAngle), startIndex(startIndex), endIndex(endIndex){} +}; + +class Bank : public PolyPattern { + public: + Bank(std::vector catPoints); + + std::variant, Polygon> poly() const override; + const std::vector& catPoints() const override; + Number coverRadius() const override; + bool isValid(const GeneralSettings& gs) const; + + private: + std::vector m_catPoints; + std::vector> m_points; + Number m_coverRadius; + Polyline m_polyline; + std::vector m_bends; + + void computeBends(); +}; +} + +#endif //CARTOCROW_BANK_H diff --git a/cartocrow/simplesets/patterns/island.cpp b/cartocrow/simplesets/patterns/island.cpp new file mode 100644 index 00000000..ba7fbbf4 --- /dev/null +++ b/cartocrow/simplesets/patterns/island.cpp @@ -0,0 +1,96 @@ +#include "island.h" +#include "bank.h" +#include "cartocrow/simplesets/helpers/cropped_voronoi.h" +#include "cartocrow/simplesets/helpers/point_voronoi_helpers.h" +#include +#include +#include + +namespace cartocrow::simplesets { +Polygon convexHull(const std::vector>& points) { + std::vector> chPoints; + CGAL::convex_hull_2(points.begin(), points.end(), std::back_inserter(chPoints)); + return {chPoints.begin(), chPoints.end()}; +} + +std::optional> convexIntersection(const Polygon& p, const Polygon& q) { + std::vector> intersection_result; + CGAL::intersection(p, q, std::back_inserter(intersection_result)); + if (intersection_result.empty()) { + return std::nullopt; + } else { + return intersection_result[0].outer_boundary(); + } +} + +Number coverRadiusOfPoints(const std::vector>& points) { + DT dt; + auto exact_points = makeExact(points); + dt.insert(exact_points.begin(), exact_points.end()); + + Rectangle bbox = CGAL::bounding_box(points.begin(), points.end()); + std::optional> squaredCoverRadius; + + auto hull = makeExact(convexHull(points)); + Cropped_voronoi_from_delaunay cropped_voronoi(hull); + + for (auto eit = dt.finite_edges_begin(); eit != dt.finite_edges_end(); ++eit) { + CGAL::Object o = dt.dual(eit); + Line l; + Ray r; + Segment s; + auto site = eit->first->vertex(dt.cw(eit->second))->point(); + if(CGAL::assign(s,o)) cropped_voronoi << std::pair(site, s); + if(CGAL::assign(r,o)) cropped_voronoi << std::pair(site, r); + if(CGAL::assign(l,o)) cropped_voronoi << std::pair(site, l); + } + + for (const auto& [site, seg] : cropped_voronoi.m_cropped_vd) { + for (const auto& v : {seg.source(), seg.target()}) { + auto d = squared_distance(approximate(v), approximate(site)); + if (!squaredCoverRadius.has_value() || d > *squaredCoverRadius) { + squaredCoverRadius = d; + } + } + } + + return sqrt(*squaredCoverRadius); +} + +Island::Island(std::vector catPoints): m_catPoints(std::move(catPoints)) { + // Store the point positions separately, sometimes only the positions are needed. + std::transform(m_catPoints.begin(), m_catPoints.end(), std::back_inserter(m_points), [](const CatPoint& cp) { + return cp.point; + }); + + if (m_catPoints.size() >= 3) { + bool collinear = true; + for (int i = 0; i < m_catPoints.size() - 2; ++i) { + if (!CGAL::collinear(m_catPoints[i].point, m_catPoints[i+1].point, m_catPoints[i+2].point)) { + collinear = false; + } + } + if (collinear) { + Bank bank(m_catPoints); + m_coverRadius = bank.coverRadius(); + m_poly = bank.poly(); + return; + } + } + + m_coverRadius = coverRadiusOfPoints(m_points); + m_poly = convexHull(m_points); +} + +std::variant, Polygon> Island::poly() const { + return m_poly; +} + +const std::vector& Island::catPoints() const { + return m_catPoints; +} + +Number Island::coverRadius() const { + return m_coverRadius; +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/patterns/island.h b/cartocrow/simplesets/patterns/island.h new file mode 100644 index 00000000..3668e4c3 --- /dev/null +++ b/cartocrow/simplesets/patterns/island.h @@ -0,0 +1,27 @@ +#ifndef CARTOCROW_ISLAND_H +#define CARTOCROW_ISLAND_H + +#include "pattern.h" +#include "poly_pattern.h" +#include + +namespace cartocrow::simplesets { +typedef CGAL::Delaunay_triangulation_2 DT; + +class Island : public PolyPattern { + public: + Island(std::vector catPoints); + + std::variant, Polygon> poly() const override; + const std::vector& catPoints() const override; + Number coverRadius() const override; + + private: + std::vector m_catPoints; + std::vector> m_points; + Number m_coverRadius; + std::variant, Polygon> m_poly; +}; +} + +#endif //CARTOCROW_ISLAND_H diff --git a/cartocrow/simplesets/patterns/matching.cpp b/cartocrow/simplesets/patterns/matching.cpp new file mode 100644 index 00000000..16651670 --- /dev/null +++ b/cartocrow/simplesets/patterns/matching.cpp @@ -0,0 +1,19 @@ +#include "matching.h" + +namespace cartocrow::simplesets { +Matching::Matching(const CatPoint& catPoint1, const CatPoint& catPoint2) : m_catPoints({catPoint1, catPoint2}) {} + +std::variant, Polygon> Matching::poly() const { +// std::vector> pts({m_catPoints.first.point, m_catPoints.second.point}); +// return Polygon(pts.begin(), pts.end()); + return Polyline(std::vector({m_catPoints[0].point, m_catPoints[1].point})); +} + +const std::vector& Matching::catPoints() const { + return m_catPoints; +} + +Number Matching::coverRadius() const { + return sqrt(squared_distance(m_catPoints[0].point, m_catPoints[1].point)) / 2; +} +} diff --git a/cartocrow/simplesets/patterns/matching.h b/cartocrow/simplesets/patterns/matching.h new file mode 100644 index 00000000..89cfe615 --- /dev/null +++ b/cartocrow/simplesets/patterns/matching.h @@ -0,0 +1,19 @@ +#ifndef CARTOCROW_MATCHING_H +#define CARTOCROW_MATCHING_H + +#include "pattern.h" +#include "poly_pattern.h" + +namespace cartocrow::simplesets { +class Matching : public PolyPattern { + public: + Matching(const CatPoint& catPoint1, const CatPoint& catPoint2); + std::variant, Polygon> poly() const override; + const std::vector& catPoints() const override; + Number coverRadius() const override; + + private: + std::vector m_catPoints; +}; +} +#endif //CARTOCROW_MATCHING_H diff --git a/cartocrow/simplesets/patterns/pattern.cpp b/cartocrow/simplesets/patterns/pattern.cpp new file mode 100644 index 00000000..f282a422 --- /dev/null +++ b/cartocrow/simplesets/patterns/pattern.cpp @@ -0,0 +1,2 @@ +#include "pattern.h" + diff --git a/cartocrow/simplesets/patterns/pattern.h b/cartocrow/simplesets/patterns/pattern.h new file mode 100644 index 00000000..ae419bad --- /dev/null +++ b/cartocrow/simplesets/patterns/pattern.h @@ -0,0 +1,23 @@ +#ifndef CARTOCROW_PATTERN_H +#define CARTOCROW_PATTERN_H + +#include "cartocrow/core/core.h" +#include "cartocrow/core/polyline.h" +#include "cartocrow/simplesets/cat_point.h" +#include "cartocrow/simplesets/settings.h" +#include "cartocrow/simplesets/types.h" +#include + +namespace cartocrow::simplesets { +class Pattern { + public: + virtual std::variant, Polygon, CSPolygon> contour() const = 0; + virtual const std::vector& catPoints() const = 0; + /// Return the category to which the points of this category belong.. + int category() const { + return catPoints()[0].category; + } +}; +} + +#endif //CARTOCROW_PATTERN_H diff --git a/cartocrow/simplesets/patterns/poly_pattern.h b/cartocrow/simplesets/patterns/poly_pattern.h new file mode 100644 index 00000000..15b0738a --- /dev/null +++ b/cartocrow/simplesets/patterns/poly_pattern.h @@ -0,0 +1,36 @@ +#ifndef CARTOCROW_POLYPATTERN_H +#define CARTOCROW_POLYPATTERN_H + +#include "pattern.h" + +namespace cartocrow::simplesets { +template +struct variant_cast_proxy +{ + std::variant v; + + template + operator std::variant() const + { + return std::visit([](auto&& arg) -> std::variant { return arg ; }, + v); + } +}; + +template +auto variant_cast(const std::variant& v) -> variant_cast_proxy +{ + return {v}; +} + +class PolyPattern : public Pattern { + public: + virtual std::variant, Polygon> poly() const = 0; + std::variant, Polygon, CSPolygon> contour() const override { + return variant_cast(poly()); + }; + virtual Number coverRadius() const = 0; +}; +} + +#endif //CARTOCROW_POLYPATTERN_H diff --git a/cartocrow/simplesets/patterns/single_point.cpp b/cartocrow/simplesets/patterns/single_point.cpp new file mode 100644 index 00000000..63a075f2 --- /dev/null +++ b/cartocrow/simplesets/patterns/single_point.cpp @@ -0,0 +1,22 @@ +#include "single_point.h" + +namespace cartocrow::simplesets { +SinglePoint::SinglePoint(CatPoint catPoint): m_catPoints({catPoint}) {} + +std::variant, Polygon> SinglePoint::poly() const { + std::vector> pts({m_catPoints[0].point}); + return Polygon(pts.begin(), pts.end()); +} + +const std::vector& SinglePoint::catPoints() const { + return m_catPoints; +} + +Number SinglePoint::coverRadius() const { + return 0; +} + +CatPoint SinglePoint::catPoint() const { + return m_catPoints[0]; +} +} \ No newline at end of file diff --git a/cartocrow/simplesets/patterns/single_point.h b/cartocrow/simplesets/patterns/single_point.h new file mode 100644 index 00000000..fa40d976 --- /dev/null +++ b/cartocrow/simplesets/patterns/single_point.h @@ -0,0 +1,21 @@ +#ifndef CARTOCROW_SINGLE_POINT_H +#define CARTOCROW_SINGLE_POINT_H + +#include "pattern.h" +#include "poly_pattern.h" + +namespace cartocrow::simplesets { +class SinglePoint : public PolyPattern { + public: + SinglePoint(CatPoint catPoint); + std::variant, Polygon> poly() const override; + const std::vector& catPoints() const override; + Number coverRadius() const override; + CatPoint catPoint() const; + + private: + std::vector m_catPoints; +}; +} + +#endif //CARTOCROW_SINGLE_POINT_H diff --git a/cartocrow/simplesets/settings.h b/cartocrow/simplesets/settings.h new file mode 100644 index 00000000..802271aa --- /dev/null +++ b/cartocrow/simplesets/settings.h @@ -0,0 +1,72 @@ +#ifndef CARTOCROW_SETTINGS_H +#define CARTOCROW_SETTINGS_H + +#include "types.h" + +namespace cartocrow::simplesets { +struct GeneralSettings { + /// Radius of circle that represents a point. + Number pointSize; + /// Maximum number of inflections a bank is allowed to have. + int inflectionLimit; + /// Maximum total angle of a bend (maximum monotone subsequence of a bank). + Number maxBendAngle; + /// Maximum turning angle in a bank. + Number maxTurnAngle; + + /// The distance each pattern is dilated. + Number dilationRadius() const { + return pointSize * 3; + } +}; + +struct PartitionSettings { + /// Create banks? + bool banks; + /// Create islands? + bool islands; + /// Delay merges that create patterns whose points are not distributed 'regularly'? + /// A pattern is not regular if it has clearly discernible sub-patterns. + bool regularityDelay; + /// Delay merges that create patterns that intersect points. + bool intersectionDelay; + /// Disallow merges that have a point within distance admissibleRadiusFactor * dilationRadius. + Number admissibleRadiusFactor; +}; + +struct ComputeDrawingSettings { + /// Aim to keep a disk around each point visible of radius cutoutRadiusFactor * dilationRadius. + Number cutoutRadiusFactor; +}; + +struct DrawSettings { + std::vector colors; + Number whiten; + Number pointStrokeWeight(GeneralSettings gs) const { + return gs.pointSize / 2.5; + } + Number contourStrokeWeight(GeneralSettings gs) const { + return gs.pointSize / 3.5; + } + Color getColor(int category) const { + Color fillColor; + if (category > colors.size() || category < 0) { + std::cerr << "Warning! No color specified for category " << category << std::endl; + fillColor = Color{240, 240, 240}; + } else { + fillColor = colors[category]; + } + return fillColor; + } +}; + +struct Settings { + GeneralSettings gs; + PartitionSettings ps; + ComputeDrawingSettings cds; + DrawSettings ds; +}; +} + + +#endif //CARTOCROW_SETTINGS_H diff --git a/cartocrow/simplesets/types.cpp b/cartocrow/simplesets/types.cpp new file mode 100644 index 00000000..1de371f3 --- /dev/null +++ b/cartocrow/simplesets/types.cpp @@ -0,0 +1,29 @@ +#include "types.h" + +namespace cartocrow::simplesets { +std::vector> makeExact(const std::vector>& points) { + std::vector> exact_points; + std::transform(points.begin(), points.end(), std::back_inserter(exact_points), + [](const Point& pt) { return makeExact(pt); }); + return exact_points; +} + +Point makeExact(const Point& point) { + return {point.x(), point.y()}; +} + +Circle makeExact(const Circle& circle) { + return {makeExact(circle.center()), circle.squared_radius()}; +} + +Polygon makeExact(const Polygon& polygon) { + std::vector> exact_points; + std::transform(polygon.vertices_begin(), polygon.vertices_end(), std::back_inserter(exact_points), + [](const Point& pt) { return makeExact(pt); }); + return {exact_points.begin(), exact_points.end()}; +} + +Point approximateAlgebraic(const CSTraits::Point_2& algebraic_point) { + return {CGAL::to_double(algebraic_point.x()), CGAL::to_double(algebraic_point.y())}; +} +} diff --git a/cartocrow/simplesets/types.h b/cartocrow/simplesets/types.h new file mode 100644 index 00000000..146df7c7 --- /dev/null +++ b/cartocrow/simplesets/types.h @@ -0,0 +1,43 @@ +#ifndef CARTOCROW_TYPES_H +#define CARTOCROW_TYPES_H + +#include "../core/core.h" +#include +#include +#include +#include +#include +#include +#include + +#include "general_polyline.h" + +namespace cartocrow::simplesets { +typedef CGAL::Arr_circle_segment_traits_2 CSTraits; +typedef CGAL::Gps_circle_segment_traits_2 CSTraitsBoolean; +typedef CGAL::Arr_polycurve_traits_2 PolyCSTraits; +typedef CSTraitsBoolean::Polygon_2 CSPolygon; +typedef CSTraitsBoolean::Polygon_with_holes_2 CSPolygonWithHoles; +typedef CGAL::General_polygon_set_2 CSPolygonSet; +typedef General_polyline_2 CSPolyline; +typedef PolyCSTraits::Curve_2 CSPolycurve; +typedef PolyCSTraits::X_monotone_curve_2 CSPolycurveXM; +typedef CGAL::Arrangement_2 CSArrangement; +typedef CSTraits::X_monotone_curve_2 X_monotone_curve_2; +typedef CSTraits::Curve_2 Curve_2; +typedef CSTraits::CoordNT OneRootNumber; +typedef CSTraits::Point_2 OneRootPoint; + +Point makeExact(const Point& point); +Circle makeExact(const Circle& circle); +std::vector> makeExact(const std::vector>& points); +Polygon makeExact(const Polygon& polygon); +Point approximateAlgebraic(const CSTraits::Point_2& algebraic_point); + +template +T squared(T x) { + return x * x; +} +} + +#endif //CARTOCROW_TYPES_H diff --git a/data/SOURCES.md b/data/SOURCES.md index 542c2478..7e845268 100644 --- a/data/SOURCES.md +++ b/data/SOURCES.md @@ -4,14 +4,31 @@ The europe.ipe file is based on data from [naturalearthdata.com](https://www.nat The file example_isolines.ipe contains isolines generated from [The Reference Elevation Model of Antarctica (REMA)](https://www.pgc.umn.edu/data/rema/). These DEMs were provided by the Byrd Polar and Climate Research Center and the Polar Geospatial Center under NSF-OPP awards 1043681, 1542736, 1543501, 1559691, 1810976, and 2129685. -See the end of the file for the relevant bibliography. +See the end of the file for the relevant bibliography [3, 4]. + +The nyc.txt file contains the locations of hotels, subway entrances, and medical clinics (green) in lower Manhattan. +It is a common benchmark dataset that originates from the paper that introduced Bubble Sets [1]. +The point locations were manually traced from the image in the paper. + +The diseasome.txt file contains vertices of an embedded graph of disorders, from the human +disease network constructed by Goh et al. [2]. The dataset consists of 516 vertices (disorders) of +twenty-one disorder classes. We use the graph layout of a poster by Bastian and Heymann archived at +https://web.archive.org/web/20121116145141/http://diseasome.eu/data/diseasome_poster.pdf. ## References -Howat, Ian; Porter, Claire; Noh, Myoung-Jon; Husby, Erik; Khuvis, Samuel; Danish, Evan; +[1] C. Collins, G. Penn, and S. Carpendale. Bubble Sets: Revealing set relations with isocontours over existing visualizations. +IEEE Transactions on Visualization and Computer Graphics, 15(6):1009–1016, 2009. +doi: https://doi.org/10.1109/TVCG.2009.122 + +[2] K.-I. Goh, M. E. Cusick, D. Valle, B. Childs, M. Vidal, and A.-L. Barabási. +The human disease network. Proceedings of the National Academy of +Sciences, 104(21):8685–8690, 2007. doi: https://doi.org/10.1073/pnas.0701361104 + +[3] Howat, Ian; Porter, Claire; Noh, Myoung-Jon; Husby, Erik; Khuvis, Samuel; Danish, Evan; Tomko, Karen; Gardiner, Judith; Negrete, Adelaide; Yadav, Bidhyananda; Klassen, James; Kelleher, Cole; Cloutier, Michael; Bakker, Jesse; Enos, Jeremy; Arnold, Galen; Bauer, Greg; Morin, Paul, 2022, "The Reference Elevation Model of Antarctica - Mosaics, Version 2", https://doi.org/10.7910/DVN/EBW8UC, Harvard Dataverse, V1 -Howat, I. M., Porter, C., Smith, B. E., Noh, M. J., & Morin, P. (2019). The reference elevation model of Antarctica. +[4] Howat, I. M., Porter, C., Smith, B. E., Noh, M. J., & Morin, P. (2019). The reference elevation model of Antarctica. The Cryosphere, 13(2), 665-674. diff --git a/data/diseasome.txt b/data/diseasome.txt new file mode 100644 index 00000000..24d00932 --- /dev/null +++ b/data/diseasome.txt @@ -0,0 +1,516 @@ +0 2502.7 -327.4 +1 2501.7 789.6 +1 2132.6 603.4 +2 2643.3 -760.5 +3 2396 -65.3 +4 1260.8 448.1 +5 2060.7 662.1 +1 2156.3 391.1 +6 1532.5 -743.2 +7 2472 425.2 +8 2740.6 -593.3 +3 2490.9 670.3 +9 2782.4 640.8 +10 1579.8 279.9 +0 2353.9 139.5 +5 2016.5 845.8 +11 1376.1 192.2 +12 2654.4 53.3 +2 2570.7 -782.1 +10 2451.3 -560.5 +2 2638.4 -476.4 +11 3194.1 -457 +13 2549.8 -803.2 +14 2572.2 -548 +0 2710.6 806.5 +15 2448.9 -253.7 +2 2469 -716.8 +16 1386.3 -315.6 +12 2801.6 -713 +15 2325.2 652.5 +1 2133 263.9 +0 2733.1 450.4 +15 3250.3 378.5 +14 1280.9 -529.7 +14 1901.4 794.4 +12 2145.4 334.8 +8 2561.7 -107 +2 2927.2 -415 +2 2660.3 -713.6 +6 1843.9 774.1 +5 1243.6 -314 +8 3079.3 -222.8 +17 3173.1 397.8 +4 1703.6 596 +8 2945.1 618.1 +11 1841 36.3 +11 2287.2 346.2 +5 1123.8 -423 +16 1867.5 289.2 +4 2399.5 -509.8 +2 2520.7 -500.8 +0 2354.6 782.5 +11 1825.3 134.6 +18 1701.9 497.1 +1 2320.8 -346.1 +2 3068.3 617.2 +19 2457.6 -529.2 +11 1631.1 -236.2 +10 2597.9 152.3 +2 3182 -546.5 +1 2621.2 864.4 +16 2798.2 -827 +8 2991.8 640.1 +5 1550.6 -314.3 +20 1444.5 -294.4 +7 2522.2 722.4 +2 2244.9 -493.2 +4 2869.9 -661 +0 2043.2 8.2 +15 3271.3 437.5 +4 2702.9 -146.3 +14 1325.8 -508.1 +13 1784.5 -117.7 +1 2556.8 886.9 +0 2347.1 759.1 +0 2456.7 576.8 +1 2344.8 424.9 +1 2834.9 -146.4 +2 2630.2 -571.1 +15 2287 -314.8 +16 1986.8 310.7 +0 1588.3 -637 +11 1902.1 -12.4 +12 2330.5 33.4 +8 3155.7 -76.4 +16 2799.4 -802.4 +16 1208.3 -340.4 +5 2072.6 736.4 +4 2029.3 566.7 +4 2040.5 -40.6 +12 2653.1 -225 +12 2815.7 -417.1 +4 2476.4 -958.1 +6 2493.1 353 +10 1375.8 278.7 +1 2625.4 728.5 +7 1502.1 -821.3 +11 1269.2 -573.3 +10 1523.2 236.8 +1 2564.5 842 +5 1344.1 -190.9 +21 2413.2 282.2 +16 2670.8 -902.8 +15 2591.8 319.3 +2 2879.7 -458.4 +0 1741.2 -100.5 +16 1766 735.2 +2 2985.3 -502.8 +4 2382.6 71.3 +5 2233.9 952.9 +2 2662.4 -372.8 +7 2311.1 582.5 +2 2441.8 -692.2 +7 2609.2 130.6 +11 2477.8 304.9 +5 2112.3 993 +10 1490.8 74.9 +2 2964 -737.8 +2 2688.7 -352.1 +2 3147.4 -653.4 +2 2652.6 -644.2 +1 1593.9 471.7 +15 2854.6 226.1 +2 3002.1 -683.2 +6 1448.7 -633.2 +6 1774 449 +5 2045.3 957.3 +4 2529.7 -909 +15 2695.8 318 +15 2690 567.6 +0 2950.8 186.5 +17 2565.5 596.1 +2 2693.3 -547.4 +4 1607 617.1 +11 2094.7 510.8 +0 2391.1 13.7 +7 2632 668.8 +19 2787.5 728 +1 2669.7 107.4 +5 1663.7 -722.5 +15 2904.6 122.7 +4 2368.9 -490.4 +10 1413 340.4 +20 3292.6 398.2 +18 1817.3 467.7 +2 2964 -357.5 +4 2730.9 585.1 +15 2342.5 -212.4 +7 1563 -541.9 +15 3082.6 258.5 +11 2763.5 470.5 +16 2256.8 -254.9 +2 3170.5 -432.6 +15 1781 -465.1 +10 1447 171.9 +2 2644.3 -427 +5 1165.4 -395.7 +8 3198.4 -97.8 +15 3242.3 345.2 +2 2963.5 533 +15 2515.2 334.4 +2 2472.9 -579.3 +21 3031.8 -896.4 +20 2328.6 115.3 +1 2284.1 -366.6 +11 2440.3 240.4 +14 1493.8 -607.7 +15 3222.5 238.5 +13 1830 -443.8 +2 2713.7 -399.4 +11 2432.1 445.2 +16 1941 754.8 +20 2594.9 -1060.8 +4 2692.3 -982.4 +2 2992.3 -797.9 +2 2951.1 702.2 +2 3177.1 -677.7 +2 2734.9 -453.2 +11 1724.3 195.6 +2 2648.9 -183.3 +1 2381.8 602 +8 3176.5 -131.7 +2 2871.4 664.9 +1 2201.7 412 +12 2619.8 -164.3 +1 2634.3 614.8 +4 1774.1 9.2 +2 1565.5 -254.4 +13 1831.8 -35.5 +8 1203.6 470.2 +13 2107.6 -78.1 +0 2741.6 407.7 +15 2851.3 617.4 +1 1674.3 406.5 +5 2161.2 869.1 +8 2759 -164.3 +21 3013.4 -832.7 +17 2802 71.4 +2 3018.7 565.3 +4 1359.8 -87.4 +10 1407.2 215 +2 3160.3 -738.1 +11 2553.3 194.5 +4 1601.7 528.9 +2 2746.7 667.6 +21 2132.8 641.3 +5 2280.1 622 +11 2529.8 402.8 +9 2477 259.8 +12 2582.7 -691.2 +15 2858.9 53.8 +12 2567.5 228.6 +6 1178.6 384.1 +2 2604.5 -593.1 +2 2487.2 -384.2 +2 3156.3 -570.8 +5 2074.8 800.2 +10 1885 370.8 +12 2411.5 324.1 +15 2951 163.3 +12 2359.8 259.7 +15 1821.8 717.4 +3 2465.6 -67.2 +2 3035.7 728.5 +5 1432.5 -64.2 +15 2047.9 529.8 +5 2275.4 715.2 +13 1725 -10.9 +8 2510.9 -88.3 +4 3143.7 321.9 +5 1377.6 -150.6 +15 2689.2 705.6 +2 2252.9 -423.7 +4 2211.6 -273.8 +1 2605.5 635.3 +4 2115.1 547.6 +4 2178.8 -294.2 +5 2077.2 894.3 +16 1852.7 818.7 +15 1951.8 -33.6 +2 2796.1 -566.4 +15 2868.2 165.4 +5 2152.3 914.3 +4 2365.9 158.5 +10 1442.9 119 +8 3099.3 -368.5 +0 2353.4 -134.1 +1 2358.3 -387.9 +8 3056.9 674.5 +15 2589.7 298.9 +15 2364 -273.5 +8 2967.6 598.4 +11 1855.3 230 +0 2655.1 788.9 +8 3133.8 -272.4 +4 2459.6 -598.8 +12 2654.5 340.3 +6 1689.3 657.8 +15 2560.7 378 +0 2717.7 34.6 +6 1754.6 428.1 +21 2604.6 805.9 +0 2233.8 780.7 +16 2594.5 -933.3 +14 2649 -1035.5 +12 2214.3 13.9 +16 1894.9 250.9 +8 1804.8 -298.8 +4 1999.3 -58.4 +2 2776.7 -498.5 +2 2471 -477.8 +2 3027.6 -428 +4 1315.1 -551.4 +11 2301.2 303 +8 3140.8 -181 +13 1677 -78.1 +15 2468.1 30.7 +15 3120.8 369.9 +2 2603 -304.9 +2 2819.8 -693.1 +16 2527.3 -1034.7 +2 2622.9 -837.6 +20 1236.9 -423.7 +2 3046.7 -623.6 +15 3258.6 259.6 +2 2799 -476.3 +15 2773.8 429.5 +16 1926.2 732.4 +5 2103.2 712 +13 1769.4 -401.8 +15 2739.5 14 +9 2452.3 -25.1 +15 2894.9 245.8 +8 3097.8 772.5 +13 1717.5 -422.5 +15 2204.9 -334.1 +17 2548.3 692.2 +4 1605.6 427 +5 2103 756.3 +2 3157.8 -594.6 +1 2468.6 -782.7 +13 1734.6 54.9 +0 2624.5 769.8 +10 1477.5 301.2 +4 1667.1 712.3 +1 2180.2 576.5 +4 2776.3 684.5 +9 2843.4 750.2 +15 2811 11.2 +5 1224.8 -475.3 +4 2913.1 -709.6 +4 2717 549.9 +21 1774.2 -355.1 +2 2588.6 -452.2 +2 2719.2 -202 +15 1794.4 -253.9 +15 3161.1 280.8 +4 2925.6 -778.4 +15 1452.2 -170.2 +11 1324.1 299.6 +15 2407 401.5 +13 1806.6 -75.1 +12 2592.2 88.9 +4 1526.1 -274.8 +5 1662.3 839.7 +8 1151.3 427.1 +8 2388.1 218.2 +4 2271.7 676.8 +5 1996.6 914.6 +12 2638.3 170.5 +11 1803.3 203 +11 1664.3 153.5 +10 1533.6 93.9 +15 1947.1 595 +11 1816.4 271.2 +4 1600.4 132.3 +1 2766.3 86.9 +2 3018.6 -533.6 +5 1175.4 -447.5 +13 1661.5 30 +0 1620.4 -293.5 +8 2474.3 100.8 +2 2801.3 -349.9 +15 2368.7 559.4 +0 2391.3 734.5 +15 3153.8 214.8 +1 2527.8 818.3 +0 1592.2 257.7 +6 1490.6 -714.8 +15 1809 529 +5 2065.3 776.2 +0 2389.8 697 +8 3177.7 -56.2 +1 2849.1 -127.2 +5 1615 -765.9 +10 1748.1 572.9 +11 1905 66.2 +11 1711.5 112.2 +1 2917.3 204.8 +15 1649.3 -705.9 +4 2449.9 -429.7 +5 2199.2 694.2 +2 2712.5 -289.7 +2 3122.6 649.7 +15 2886 303.3 +17 2504 653.2 +7 3124.4 -763.6 +4 1914.7 644 +1 2277 541.2 +2 3010.4 -391 +4 1769.4 387.1 +13 2541.6 -350.4 +2 2821.8 -329 +10 1427.2 257.4 +1 2253.3 521.1 +12 2160.4 53 +5 1987.1 874.5 +5 2039.4 933.5 +12 2610.5 71.4 +17 2479.9 631.7 +1 1660 451.2 +1 2241.2 492.1 +16 1612.7 806.5 +8 2813.8 -537.4 +2 2830.3 -518.4 +0 1957.5 -78.6 +2 2912.2 -627.5 +13 1730.3 -380.1 +15 2922.1 283.2 +10 1578.5 320.1 +8 1220 405.7 +2 2894.4 -478.5 +18 1821.8 493 +4 2764.4 604.3 +5 1414.8 -109.1 +6 1561.8 -844.2 +6 1793.5 796.1 +16 2467.9 -937.7 +2 3117.9 -625.4 +2 2764.9 -732.6 +6 1624.9 761.9 +19 2619.4 204.6 +4 3174.3 -392.4 +6 1547.1 -801.3 +14 1734.3 636.8 +20 1982.7 331.3 +1 2238.2 452.3 +2 3009.4 -658.7 +4 1396.4 151.5 +0 2426.5 180.3 +14 2628.5 -5.5 +0 1739.3 -55.1 +4 1630.3 551.2 +4 1913.9 10.3 +15 1806.6 548.9 +12 2281.2 92.2 +2 2797 -256.2 +0 1997.6 32.7 +12 2579.8 174.3 +4 3249.3 308.5 +0 2457.4 -109.8 +2 2985.8 514.1 +16 1868.9 697.8 +10 2703.1 -247.3 +5 1160.4 -314.3 +15 2858.3 9.1 +4 2618.3 -520.5 +2 2863.8 -309.1 +2 3031.8 -454.3 +4 3130.5 -345.6 +16 2017.1 286.8 +19 2200.6 369.9 +21 2492.7 467.4 +17 3219.4 417.1 +4 2742.5 -127.5 +2 2205.4 -445.3 +8 2542.4 445.6 +8 3100.7 300.8 +7 2783.8 -182.7 +15 2854.8 143.6 +0 2237.2 431.3 +16 2799.9 -894.4 +15 2607.2 276.7 +10 1555.3 573.7 +20 2781 -107 +14 1358.1 -597.2 +20 1950.1 351.8 +2 2874.3 -288 +4 1831.5 153.2 +2 3213.9 -413.3 +16 2738.7 -870.7 +6 1460 -778.7 +14 2403.9 -5.7 +5 2165.2 971.5 +17 2854 265.1 +4 3103.7 343.5 +5 1270.7 -369 +0 2588.9 749.8 +11 1776.7 341.5 +2 2804.1 -755.6 +2 2415 -410 +14 1713.1 818.9 +2 2839.2 -380.8 +2 2778.4 -435.3 +4 1728.2 678.7 +2 2704.2 -781.5 +6 1820.3 678.2 +2 2935.8 -566.4 +5 1334.7 -129.6 +0 2241.6 737.5 +12 2165 281 +5 2163.8 819.2 +4 2988.9 -874.9 +5 1120.1 -368.6 +5 1424.6 -348.7 +1 2659.8 829.5 +15 1786.5 409.1 +15 2286 -190.9 +10 1545 58.1 +16 1664.5 784.6 +10 1884.1 408.8 +0 2539.8 -25.7 +10 2337.1 -446.6 +12 2581 14.6 +12 2659.2 279.5 +12 2512.9 52.5 +4 1928.9 389.1 +2 2533.9 -740.6 +4 3084.6 -711.5 +14 1584.6 -560.9 +15 2306.4 -236.5 +21 2506.2 -45.8 +4 3055.5 -853.7 +2 2462.2 -826.8 +2 3200.2 -493.7 +2 2488.4 -761.9 +13 1689.3 -34.4 +14 2067.6 585.1 +8 2719.3 388.3 +15 2726.8 358.5 +2 1655.3 -217.7 +14 1518.6 -656.4 +1 2631.8 -125.9 +11 1758.6 175.7 +0 2349.5 372.5 +2 1235.3 361.2 +2 2613.4 -405.7 +2 2864.7 -439.7 +4 1610.4 -201 +14 1486.2 -676 +16 2575 -1008.9 +1 2142.6 308.8 +15 2262.7 -172.9 +2 2301.4 -468.5 +11 1676.4 216.9 +12 2282.8 -26.6 \ No newline at end of file diff --git a/data/nyc.txt b/data/nyc.txt new file mode 100644 index 00000000..c44d07af --- /dev/null +++ b/data/nyc.txt @@ -0,0 +1,96 @@ +0 39.211 556.19 +0 46.179 526.507 +0 61.915 524.604 +0 65.72 531.695 +0 80.419 528.928 +0 12.264 762.294 +0 41.919 776.673 +0 43.419 783.344 +0 65.187 738.367 +0 64.72 723.974 +0 73.372 708.833 +0 64.32 693.092 +0 77.498 693.492 +0 56.467 656.354 +0 115.645 682.323 +0 113.988 729.952 +0 137.154 742.503 +0 111.103 608.229 +0 143.582 593.054 +0 155.529 654.019 +0 174.164 638.711 +0 166.311 667.463 +0 178.723 694.75 +0 184.846 697.679 +0 195.096 702.338 +0 193.765 730.757 +0 174.597 745.798 +0 157.011 802.159 +0 191.335 580.202 +0 205.844 587.923 +0 221.684 576.475 +0 203.714 568.755 +0 190.803 560.236 +0 202.911 540.915 +0 231.618 535.208 +0 217.61 518.78 +2 230.821 566.587 +2 221.37 584.59 +2 228.558 582.327 +2 236.412 584.723 +2 238.142 602.693 +2 229.224 600.43 +2 223.101 610.813 +2 190.974 645.466 +2 253.43 715.841 +2 187.604 754.033 +2 199.062 797.392 +2 228.043 820.308 +2 67.299 827.721 +2 80.104 852.209 +2 180.258 578.719 +2 160.488 565.688 +2 136.224 539.403 +2 95.336 591.973 +2 97.134 561.42 +2 98.032 549.962 +1 107.019 543.391 +1 120.273 516.431 +1 187.447 511.489 +1 209.014 526.541 +1 87.698 532.382 +1 56.021 598.432 +1 101.851 602.251 +1 173.743 578.887 +1 203.398 595.062 +1 162.285 622.695 +1 163.858 628.761 +1 83.879 645.555 +1 76.015 676.558 +1 142.515 669.368 +1 136 683.185 +1 169.923 691.273 +1 124.317 685.881 +1 124.542 706.269 +1 150.378 723.792 +1 151.277 748.729 +1 109.265 754.346 +1 104.772 748.28 +1 88.372 734.127 +1 194.861 741.596 +1 177.786 784.057 +1 63.435 809.893 +1 52.271 835.573 +1 46.391 830.731 +1 41.204 826.753 +1 35.324 822.257 +1 12.324 835.919 +1 22.527 818.28 +1 43.798 804.1 +1 49.677 798.22 +1 36.016 779.5 +1 31.52 774.312 +1 51.752 738.343 +1 64.722 711.064 +1 55.557 702.59 +1 233.768 618.887 \ No newline at end of file diff --git a/demos/CMakeLists.txt b/demos/CMakeLists.txt index 9b0f4c90..44c2c34b 100644 --- a/demos/CMakeLists.txt +++ b/demos/CMakeLists.txt @@ -1,3 +1,4 @@ add_subdirectory(flow_map) add_subdirectory(isoline_simplification) +add_subdirectory(simplesets) add_subdirectory(renderer) diff --git a/demos/isoline_simplification/CMakeLists.txt b/demos/isoline_simplification/CMakeLists.txt index bb8d5ab4..5d953b36 100644 --- a/demos/isoline_simplification/CMakeLists.txt +++ b/demos/isoline_simplification/CMakeLists.txt @@ -15,4 +15,4 @@ target_link_libraries( Qt5::Widgets ) -install(TARGETS isoline_simplification DESTINATION ${INSTALL_BINARY_DIR}) +install(TARGETS isoline_simplification_demo DESTINATION ${INSTALL_BINARY_DIR}) diff --git a/demos/simplesets/CMakeLists.txt b/demos/simplesets/CMakeLists.txt new file mode 100644 index 00000000..b3662be5 --- /dev/null +++ b/demos/simplesets/CMakeLists.txt @@ -0,0 +1,21 @@ +set(SOURCES + simplesets_demo.cpp + colors.cpp +) + +add_executable(simplesets_demo ${SOURCES}) + +target_link_libraries( + simplesets_demo + PRIVATE + ${COMMON_CLA_TARGET} + core + renderer + simplesets + CGAL::CGAL + Qt5::Widgets +) + +install(TARGETS simplesets_demo DESTINATION ${INSTALL_BINARY_DIR}) + +add_subdirectory(helpers) \ No newline at end of file diff --git a/demos/simplesets/colors.cpp b/demos/simplesets/colors.cpp new file mode 100644 index 00000000..28857f6f --- /dev/null +++ b/demos/simplesets/colors.cpp @@ -0,0 +1,41 @@ +#include "colors.h" + +namespace cartocrow::CB { +Color light_blue(166,206,227); +Color blue(31,120,180); +Color light_green(178,223,138); +Color green(51,160,44); +Color light_red(251,154,153); +Color red(227,26,28); +Color light_orange(253,191,111); +Color orange(255,127,0); +Color light_purple(202,178,214); +Color purple(106,61,154); +} + +namespace cartocrow::diseasome { +std::vector colors({ + {0xB9E1EE}, + {0x9AC019}, + {0xCD6814}, + {0xE53389}, + {0xC1BC56}, + {0x923B8B}, + {0xFBD2AA}, + {0x999999}, + {0xFECD0F}, + {0xCB9A03}, + {0xF3983B}, + {0x4B8EC7}, + {0x2E9A67}, + {0xE95937}, + {0xF8EE82}, + {0xE74646}, + {0xCBBC9D}, + {0x6699CD}, + {0x6FC4C6}, + {0xF1979A}, + {0x8F5A9C}, + {0xBB3087}, +}); +} \ No newline at end of file diff --git a/demos/simplesets/colors.h b/demos/simplesets/colors.h new file mode 100644 index 00000000..bca1c07d --- /dev/null +++ b/demos/simplesets/colors.h @@ -0,0 +1,24 @@ +#ifndef CARTOCROW_COLORS_H +#define CARTOCROW_COLORS_H + +#include "cartocrow/core/core.h" + +namespace cartocrow::CB { + extern Color light_blue; + extern Color blue; + extern Color light_green; + extern Color green; + extern Color light_red; + extern Color red; + extern Color light_orange; + extern Color orange; + extern Color light_purple; + extern Color purple; +} + +namespace cartocrow::diseasome { + extern std::vector colors; +} + + +#endif //CARTOCROW_COLORS_H diff --git a/demos/simplesets/helpers/CMakeLists.txt b/demos/simplesets/helpers/CMakeLists.txt new file mode 100644 index 00000000..88100abf --- /dev/null +++ b/demos/simplesets/helpers/CMakeLists.txt @@ -0,0 +1,18 @@ +set(SOURCES + circle_convex_hull.cpp +) + +add_executable(circle_convex_hull_demo ${SOURCES}) + +target_link_libraries( + circle_convex_hull_demo + PRIVATE + ${COMMON_CLA_TARGET} + core + renderer + simplesets + CGAL::CGAL + Qt5::Widgets +) + +install(TARGETS circle_convex_hull_demo DESTINATION ${INSTALL_BINARY_DIR}) diff --git a/demos/simplesets/helpers/circle_convex_hull.cpp b/demos/simplesets/helpers/circle_convex_hull.cpp new file mode 100644 index 00000000..7c1936c4 --- /dev/null +++ b/demos/simplesets/helpers/circle_convex_hull.cpp @@ -0,0 +1,46 @@ +#include "circle_convex_hull.h" +#include +#include "cartocrow/simplesets/helpers/approximate_convex_hull.h" +#include "cartocrow/simplesets/helpers/cs_polygon_helpers.h" + +CircleConvexHullDemo::CircleConvexHullDemo() { + setWindowTitle("Convex hull of circles"); + auto renderer = new GeometryWidget(); + renderer->setDrawAxes(false); + setCentralWidget(renderer); + + std::vector> cs({ + {{0, 0}, 2}, + {{10, 4}, 12}, + {{7, -6}, 8}, + {{5, -8}, 1}, + {{3, 3}, 3}, + {{15, -4}, 9}, + {{5, -4}, 8}, + {{0, -1}, 5}, + {{5, -3}, 12}, + {{8, -9}, 16}, + }); + + auto hull = approximateConvexHull(cs); + + std::function drawFunc = [cs, hull](GeometryRenderer& renderer) { + RenderPath path = renderPath(hull); + renderer.setMode(GeometryRenderer::fill); + renderer.setFill(Color{150, 150, 150}); + for (const auto& c : cs) { + renderer.draw(c); + } + renderer.setMode(GeometryRenderer::stroke); + renderer.setStroke(Color{0, 0, 0}, 3.0); + renderer.draw(path); + }; + renderer->addPainting(drawFunc, "Disks"); +} + +int main(int argc, char* argv[]) { + QApplication app(argc, argv); + CircleConvexHullDemo demo; + demo.show(); + app.exec(); +} diff --git a/demos/simplesets/helpers/circle_convex_hull.h b/demos/simplesets/helpers/circle_convex_hull.h new file mode 100644 index 00000000..716e35dc --- /dev/null +++ b/demos/simplesets/helpers/circle_convex_hull.h @@ -0,0 +1,25 @@ +#ifndef CARTOCROW_CIRCLE_CONVEX_HULL_H +#define CARTOCROW_CIRCLE_CONVEX_HULL_H + +#include "cartocrow/core/core.h" +#include "cartocrow/core/ipe_reader.h" +#include "cartocrow/renderer/geometry_painting.h" +#include "cartocrow/renderer/geometry_widget.h" +#include "cartocrow/simplesets/types.h" +#include "cartocrow/simplesets/partition.h" +#include "cartocrow/simplesets/drawing_algorithm.h" +#include + +using namespace cartocrow; +using namespace cartocrow::renderer; +using namespace cartocrow::simplesets; + +class CircleConvexHullDemo: public QMainWindow { + Q_OBJECT + + public: + CircleConvexHullDemo(); + +}; + +#endif //CARTOCROW_CIRCLE_CONVEX_HULL_H diff --git a/demos/simplesets/simplesets_demo.cpp b/demos/simplesets/simplesets_demo.cpp new file mode 100644 index 00000000..08422ebf --- /dev/null +++ b/demos/simplesets/simplesets_demo.cpp @@ -0,0 +1,218 @@ +/* +The CartoCrow library implements algorithmic geo-visualization methods, +developed at TU Eindhoven. +Copyright (C) 2024 TU Eindhoven + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU 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 General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . +*/ + +#include "simplesets_demo.h" +#include "cartocrow/core/ipe_reader.h" +#include "cartocrow/renderer/ipe_renderer.h" +#include "cartocrow/simplesets/dilated/dilated_poly.h" +#include "cartocrow/simplesets/drawing_algorithm.h" +#include "cartocrow/simplesets/helpers/approximate_convex_hull.h" +#include "cartocrow/simplesets/helpers/arrangement_helpers.h" +#include "cartocrow/simplesets/helpers/cs_polygon_helpers.h" +#include "cartocrow/simplesets/helpers/cs_polyline_helpers.h" +#include "cartocrow/simplesets/helpers/poly_line_gon_intersection.h" +#include "cartocrow/simplesets/parse_input.h" +#include "cartocrow/simplesets/partition_algorithm.h" +#include "cartocrow/simplesets/partition_painting.h" +#include "cartocrow/simplesets/patterns/bank.h" +#include "cartocrow/simplesets/patterns/island.h" +#include "cartocrow/simplesets/patterns/matching.h" +#include "cartocrow/simplesets/patterns/single_point.h" +#include "colors.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace fs = std::filesystem; +using namespace cartocrow; +using namespace cartocrow::renderer; +using namespace cartocrow::simplesets; + +SimpleSetsDemo::SimpleSetsDemo() { + setWindowTitle("SimpleSets"); + + m_cds = ComputeDrawingSettings{0.675}; + + // Initial file + std::filesystem::path filePath("data/nyc.txt"); + + // nyc + m_gs = GeneralSettings{2.1, 2, M_PI, 70.0 / 180 * M_PI}; + m_ds = DrawSettings{{CB::light_blue, CB::light_red, CB::light_green, CB::light_purple, CB::light_orange}, 0.7}; + m_ps = PartitionSettings{true, true, true, true, 0.1}; + + // diseasome +// m_gs = GeneralSettings{5.204, 2, M_PI, 70.0 / 180 * M_PI}; +// m_ds = DrawSettings{diseasome::colors, 0.7}; +// m_ps = PartitionSettings{true, true, true, false, 0.5}; + + auto* dockWidget = new QDockWidget(); + addDockWidget(Qt::RightDockWidgetArea, dockWidget); + auto* vWidget = new QWidget(); + auto* vLayout = new QVBoxLayout(vWidget); + vLayout->setAlignment(Qt::AlignTop); + dockWidget->setWidget(vWidget); + + auto* basicOptions = new QLabel("

Input

"); + vLayout->addWidget(basicOptions); + auto* fileSelector = new QPushButton("Select file"); + vLayout->addWidget(fileSelector); + + auto* fitToScreenButton = new QPushButton("Fit to screen"); + vLayout->addWidget(fitToScreenButton); + + auto* settingsLabel = new QLabel("

Settings

"); + vLayout->addWidget(settingsLabel); + auto* coverLabel = new QLabel("Cover"); + vLayout->addWidget(coverLabel); + auto* coverSlider = new QSlider(Qt::Orientation::Horizontal); + vLayout->addWidget(coverSlider); + coverSlider->setMinimum(0); + coverSlider->setMaximum(80); + coverSlider->setValue(47); + + auto* ptSizeLabel = new QLabel("Point size"); + vLayout->addWidget(ptSizeLabel); + auto* ptSizeSlider = new QSlider(Qt::Orientation::Horizontal); + vLayout->addWidget(ptSizeSlider); + ptSizeSlider->setMinimum(1); + ptSizeSlider->setMaximum(80); + ptSizeSlider->setValue(static_cast(m_gs.pointSize * 10)); + + m_renderer = new GeometryWidget(); + m_renderer->setDrawAxes(false); + setCentralWidget(m_renderer); + + m_renderer->setMinZoom(0.01); + m_renderer->setMaxZoom(1000.0); + + loadFile(filePath); + computePartitions(); + computeDrawing(coverSlider->value() / 10.0); + fitToScreen(); + + connect(fitToScreenButton, &QPushButton::clicked, [this]() { + fitToScreen(); + }); + connect(fileSelector, &QPushButton::clicked, [this, fileSelector, coverSlider]() { + QString startDir = "data/"; + std::filesystem::path filePath = QFileDialog::getOpenFileName(this, tr("Select SimpleSets input"), startDir).toStdString(); + if (filePath == "") return; + loadFile(filePath); + computePartitions(); + computeDrawing(coverSlider->value() / 10.0); + fitToScreen(); + fileSelector->setText(QString::fromStdString(filePath.filename())); + }); + connect(coverSlider, &QSlider::valueChanged, [this, coverSlider] { + double value = coverSlider->value() / 10.0; + computeDrawing(value); + }); + connect(ptSizeSlider, &QSlider::valueChanged, [this, ptSizeSlider, coverSlider] { + double value = ptSizeSlider->value() / 10.0; + m_gs.pointSize = value; + computePartitions(); + computeDrawing(coverSlider->value() / 10.0); + fitToScreen(); + }); + + fitToScreenButton->click(); +} + +void SimpleSetsDemo::loadFile(const std::filesystem::path& filePath) { + std::ifstream inputStream(filePath, std::ios_base::in); + if (!inputStream.good()) { + throw std::runtime_error("Failed to read input"); + } + std::stringstream buffer; + buffer << inputStream.rdbuf(); + m_points = parseCatPoints(buffer.str()); +} + +void SimpleSetsDemo::fitToScreen() { + std::vector> pts; + std::transform(m_points.begin(), m_points.end(), std::back_inserter(pts), + [](const CatPoint& pt) { return pt.point; }); + Box box = CGAL::bbox_2(pts.begin(), pts.end()); + auto delta = 2 * m_gs.dilationRadius(); + Box expanded(box.xmin() - delta, box.ymin() - delta, box.xmax() + delta, box.ymax() + delta); + m_renderer->fitInView(expanded); +} + +void SimpleSetsDemo::resizeEvent(QResizeEvent *event) { + fitToScreen(); +} + +void SimpleSetsDemo::computePartitions(){ + m_partitions = partition(m_points, m_gs, m_ps, 8 * m_gs.dilationRadius()); +} + +void SimpleSetsDemo::computeDrawing(double cover) { + Partition* thePartition; + bool found = false; + for (auto& [time, partition] : m_partitions) { + if (time < cover * m_gs.dilationRadius()) { + thePartition = &partition; + found = true; + } + } + if (found) { + m_partition = *thePartition; + } else { + m_partition = m_partitions.front().second; + } + + m_renderer->clear(); + + auto pp = std::make_shared(m_partition, m_gs, m_ds); + m_renderer->addPainting(pp, "Partition"); + + bool wellSeparated = true; + for (const auto& p : m_points) { + for (const auto& q : m_points) { + if (p.category == q.category) continue; + if (CGAL::squared_distance(p.point, q.point) < 4 * m_gs.pointSize * m_gs.pointSize) { + wellSeparated = false; + } + } + } + if (wellSeparated) { + m_dpd = std::make_shared(m_partition, m_gs, m_cds); + auto ap = std::make_shared(*m_dpd, m_ds); + m_renderer->addPainting(ap, "Arrangement"); + } else { + std::cerr << "Points of different category are too close together; not computing a drawing." << std::endl; + } +} + +int main(int argc, char* argv[]) { + QApplication app(argc, argv); + SimpleSetsDemo demo; + demo.show(); + app.exec(); +} diff --git a/demos/simplesets/simplesets_demo.h b/demos/simplesets/simplesets_demo.h new file mode 100644 index 00000000..cd23b91b --- /dev/null +++ b/demos/simplesets/simplesets_demo.h @@ -0,0 +1,62 @@ +/* +The CartoCrow library implements algorithmic geo-visualization methods, +developed at TU Eindhoven. +Copyright (C) 2024 TU Eindhoven + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU 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 General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see . +*/ + +#ifndef CARTOCROW_SIMPLESETS_DEMO_H +#define CARTOCROW_SIMPLESETS_DEMO_H + +#include "cartocrow/core/core.h" +#include "cartocrow/core/ipe_reader.h" +#include "cartocrow/renderer/geometry_painting.h" +#include "cartocrow/renderer/geometry_widget.h" +#include "cartocrow/simplesets/types.h" +#include "cartocrow/simplesets/partition.h" +#include "cartocrow/simplesets/drawing_algorithm.h" +#include +#include + +using namespace cartocrow; +using namespace cartocrow::renderer; +using namespace cartocrow::simplesets; + +class SimpleSetsDemo : public QMainWindow { + Q_OBJECT + + public: + SimpleSetsDemo(); + void resizeEvent(QResizeEvent *event) override; + + private: + std::vector m_points; + Partition m_partition; + std::shared_ptr m_dpd; + GeneralSettings m_gs; + DrawSettings m_ds; + PartitionSettings m_ps; + ComputeDrawingSettings m_cds; + GeometryWidget* m_renderer; + std::vector, Partition>> m_partitions; + + std::shared_ptr> m_cc; + void fitToScreen(); + void loadFile(const std::filesystem::path& filePath); + void computePartitions(); + void computeDrawing(double cover); +}; + +#endif //CARTOCROW_SIMPLESETS_DEMO_H diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d93ee33d..7f5c1bfb 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -19,6 +19,9 @@ set(TEST_SOURCES "cartocrow_test.cpp" "necklace_map/range.cpp" "renderer/ipe_renderer.cpp" "simplification/vw_simplification.cpp" + "simplesets/poly_line_gon_intersection.cpp" + "simplesets/partition_algorithm.cpp" + "simplesets/collinear_island.cpp" ) add_executable(cartocrow_test cartocrow_test.cpp ${TEST_SOURCES}) @@ -29,4 +32,5 @@ target_link_libraries(cartocrow_test necklace_map renderer simplification + simplesets ) diff --git a/test/simplesets/collinear_island.cpp b/test/simplesets/collinear_island.cpp new file mode 100644 index 00000000..81b2ae8d --- /dev/null +++ b/test/simplesets/collinear_island.cpp @@ -0,0 +1,11 @@ +#include "../catch.hpp" +#include "cartocrow/simplesets/patterns/island.h" + +using namespace cartocrow; +using namespace cartocrow::simplesets; + +TEST_CASE("Collinear island") { + Island island({{0, {0, 0}}, {0, {1, 0}}, {0, {2, 0}}, {0, {3, 0}}}); + CHECK(island.coverRadius() == 0.5); + CHECK(std::holds_alternative>(island.poly())); +} \ No newline at end of file diff --git a/test/simplesets/partition_algorithm.cpp b/test/simplesets/partition_algorithm.cpp new file mode 100644 index 00000000..5e5ffbaa --- /dev/null +++ b/test/simplesets/partition_algorithm.cpp @@ -0,0 +1,24 @@ +#include "../catch.hpp" +#include "cartocrow/simplesets/partition_algorithm.h" + +using namespace cartocrow; +using namespace cartocrow::simplesets; + +TEST_CASE("Intersection delay") { + GeneralSettings gs{1, 0, 0, 0}; + PartitionSettings ps{false, true, true, true, 0}; + std::vector points({ + {0, {0, 0}}, + {0, {0, 15}}, + {0, {15, 0}}, + {0, {15, 15}}, + {1, {7.5, 18}}, + {2, {-3, 15}}, + }); + + Island p1({points[0], points[1]}); + Island p2({points[2], points[3]}); + Island p3({points[0], points[1], points[2], points[3]}); + CHECK(intersectionDelay(p3.catPoints(), p1, p2, p3, gs, ps) == 0); + CHECK(abs(intersectionDelay(points, p1, p2, p3, gs, ps) - 3 / sqrt(2)) < M_EPSILON); +} \ No newline at end of file diff --git a/test/simplesets/poly_line_gon_intersection.cpp b/test/simplesets/poly_line_gon_intersection.cpp new file mode 100644 index 00000000..0abec824 --- /dev/null +++ b/test/simplesets/poly_line_gon_intersection.cpp @@ -0,0 +1,117 @@ +#include "../catch.hpp" +#include "cartocrow/simplesets/helpers/poly_line_gon_intersection.h" +#include "cartocrow/simplesets/helpers/cs_curve_helpers.h" +#include "cartocrow/simplesets/helpers/arrangement_helpers.h" +#include "cartocrow/renderer/ipe_renderer.h" + +using namespace cartocrow; +using namespace cartocrow::renderer; +using namespace cartocrow::simplesets; + +TEST_CASE("Intersection lies in polygon and has correct orientation") { + X_monotone_curve_2 xm_curve({-2, 0}, {2, 0}); + std::vector xm_curves{xm_curve}; + CSPolyline polyline(xm_curves.begin(), xm_curves.end()); + CSPolygon disk = circleToCSPolygon({{0, 0}, 1}); + auto result = intersection(polyline, disk, false); + CHECK(result.size() == 1); + CHECK(result[0].size() == 1); + CHECK(result[0].curves_begin()->source() == OneRootPoint(-1, 0)); + CHECK(result[0].curves_begin()->target() == OneRootPoint(1, 0)); +} + +TEST_CASE("Boundary overlap") { + X_monotone_curve_2 xm_curve({-2, 0}, {2, 0}); + std::vector xm_curves{xm_curve}; + CSPolyline polyline(xm_curves.begin(), xm_curves.end()); + std::vector xm_curves_pgn{{{-4, 0}, {4, 0}}, {{4, 0}, {4, 2}}, {{4, 2}, {-4, 2}}, {{-4, 2}, {-4, 0}}}; + CSPolygon polygon(xm_curves_pgn.begin(), xm_curves_pgn.end()); + auto result1 = intersection(polyline, polygon, false); + auto result2 = intersection(polyline, polygon, true); + CHECK(result1.empty()); + CHECK(!result2.empty()); + CHECK(result2[0].curves_begin()->source() == OneRootPoint(-2, 0)); + CHECK(result2[0].curves_begin()->target() == OneRootPoint(2, 0)); +} + +TEST_CASE("Multiple and connected parts of intersection") { + Number half = 1; + half /= 2; + std::vector> points({{-1, -1}, {0, 0}, {0, -1 - half}, {1 + half, -1 - half}, {1 + half, 0}, {half, 0}}); + Polyline pl(points.begin(), points.end()); + CSPolyline polyline = polylineToCSPolyline(pl); + CSPolygon disk = circleToCSPolygon({{0, 0}, 1}); + auto result = intersection(polyline, disk, false); + CHECK(result.size() == 2); + auto r1 = std::find_if(result.begin(), result.end(), [half](const CSPolyline& pl) { + OneRootNumber negHalfSqrt2(0, -half, 2); + return pl.curves_begin()->source() == OneRootPoint(negHalfSqrt2, negHalfSqrt2); + }); + CHECK(r1 != result.end()); + CHECK(r1->size() == 2); + CHECK((++(r1->curves_begin()))->target() == OneRootPoint(0, -1)); + auto r2 = std::find_if(result.begin(), result.end(), [half](const CSPolyline& pl) { + return pl.curves_begin()->source() == OneRootPoint(1, 0); + }); + CHECK(r2 != result.end()); + CHECK(r2->size() == 1); + CHECK(r2->curves_begin()->target() == OneRootPoint(half, 0)); +} + +// Test case in https://github.com/CGAL/cgal/issues/8468 +TEST_CASE("Poly-circular-arcs that partially overlap") { + auto r = 5.204 * 3; + auto r_ = 5.204 * 3 * 0.675; + auto r2 = r * r; + auto r2_ = r_ * r_; + Circle c1({2597.9, -364.3}, r2, CGAL::CLOCKWISE); + Circle c2({2609.2, -342.6}, r2, CGAL::CLOCKWISE); + Circle c2_({2609.2, -342.6}, r2_, CGAL::CLOCKWISE); + + auto getIntersections = [](const Circle& one, const Circle& two) { + CSArrangement arr; + CGAL::insert(arr, one); + CGAL::insert(arr, two); + std::vector intersectionPoints; + for (auto vit = arr.vertices_begin(); vit != arr.vertices_end(); ++vit) { + if (vit->degree() == 4) { + intersectionPoints.push_back(vit->point()); + } + } + std::sort(intersectionPoints.begin(), intersectionPoints.end(), [](const OneRootPoint& p1, const OneRootPoint& p2) { return p1.x() < p2.x(); }); + assert(intersectionPoints.size() == 2); + return intersectionPoints; + }; + + auto isp12 = getIntersections(c1, c2); + + Curve_2 arc1(c1, isp12[0], isp12[1]); + Curve_2 arc2(c2, isp12[1], isp12[0]); + std::vector pgnArcs({arc1, arc2}); + + PolyCSTraits traits; + auto construct = traits.construct_curve_2_object(); + CSPolycurve pgnCurve(pgnArcs.begin(), pgnArcs.end()); + + CSArrangement arr; + CGAL::insert(arr, c1); + CGAL::insert(arr, c2); + CGAL::insert(arr, c2_); + + CSArrangement::Face_handle fh; + for (auto eit = arr.halfedges_begin(); eit != arr.halfedges_end(); ++eit) { + if (eit->source()->point() == isp12[0]) { + fh = eit->face(); + } + } + + std::vector xm_curves; + curvesToXMonotoneCurves(pgnArcs.begin(), pgnArcs.end(), std::back_inserter(xm_curves)); + CSPolygon pgn(xm_curves.begin(), xm_curves.end()); + + auto plnPoly = ccb_to_polygon(fh->outer_ccb()); + CSPolyline pln(plnPoly.curves_begin(), plnPoly.curves_end()); + intersection(pln, pgn, true); +} + +// todo: test difference and holes