From c9038d7c7386555f1fd0ad406be26f3fe04fc5ed Mon Sep 17 00:00:00 2001 From: Shepard Liu Date: Wed, 20 Aug 2025 10:11:59 +0800 Subject: [PATCH] fixed drawing out of bounds triangles near polar & cleanup --- .../Arrangement_on_surface_2/CMakeLists.txt | 6 +- .../Arrangement_on_surface_2/draw_agas.cpp | 105 ++++ .../Arrangement_on_surface_2/draw_arr.cpp | 158 +++++- .../draw_arr_degen_holes.cpp | 44 -- .../draw_arr_spherical_lakes.cpp | 35 -- .../draw_arr_spherical_polar_faces.cpp | 28 -- .../spherical_insert.cpp | 2 - .../spherical_overlay.cpp | 5 +- .../unbounded_linear.cpp | 51 -- .../Arr_geodesic_arc_on_sphere_traits_2.h | 21 +- .../include/CGAL/Arr_linear_traits_2.h | 42 +- .../CGAL/Draw_aos/Arr_approximation_cache.h | 95 +--- .../Draw_aos/Arr_bounded_approximate_face.h | 358 ++++++-------- .../Arr_bounded_approximate_halfedge.h | 418 ++++++++++------ .../Draw_aos/Arr_bounded_approximate_vertex.h | 22 +- .../Draw_aos/Arr_bounded_face_triangulator.h | 367 +++++++------- .../CGAL/Draw_aos/Arr_bounded_renderer.h | 26 +- .../CGAL/Draw_aos/Arr_coordinate_converter.h | 117 +++++ .../CGAL/Draw_aos/Arr_face_point_generator.h | 101 ++++ .../include/CGAL/Draw_aos/Arr_graph_conn.h | 85 ---- .../include/CGAL/Draw_aos/Arr_portals.h | 96 ---- .../include/CGAL/Draw_aos/Arr_projector.h | 88 ---- .../CGAL/Draw_aos/Arr_render_context.h | 55 ++- .../include/CGAL/Draw_aos/Arr_viewer.h | 378 +++++++------- .../include/CGAL/Draw_aos/type_utils.h | 225 ++++++--- .../include/CGAL/draw_arrangement_2.h | 466 +++++++++++++----- 26 files changed, 1915 insertions(+), 1479 deletions(-) create mode 100644 Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_agas.cpp delete mode 100644 Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_degen_holes.cpp delete mode 100644 Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_spherical_lakes.cpp delete mode 100644 Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_spherical_polar_faces.cpp delete mode 100644 Arrangement_on_surface_2/examples/Arrangement_on_surface_2/unbounded_linear.cpp create mode 100644 Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_coordinate_converter.h create mode 100644 Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_face_point_generator.h delete mode 100644 Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_graph_conn.h delete mode 100644 Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_portals.h delete mode 100644 Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_projector.h diff --git a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/CMakeLists.txt b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/CMakeLists.txt index 4690e4467b7..c8af32eca1a 100644 --- a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/CMakeLists.txt +++ b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/CMakeLists.txt @@ -18,10 +18,7 @@ endforeach() if(CGAL_Qt6_FOUND) target_link_libraries(draw_arr PRIVATE CGAL::CGAL_Basic_viewer) - target_link_libraries(draw_arr_degen_holes PRIVATE CGAL::CGAL_Basic_viewer) - target_link_libraries(draw_arr_spherical_lakes PRIVATE CGAL::CGAL_Basic_viewer) - target_link_libraries(draw_arr_spherical_polar_faces PRIVATE CGAL::CGAL_Basic_viewer) - target_link_libraries(unbounded_linear PRIVATE CGAL::CGAL_Basic_viewer) + target_link_libraries(draw_agas PRIVATE CGAL::CGAL_Basic_viewer) target_link_libraries(unbounded_non_intersecting PRIVATE CGAL::CGAL_Basic_viewer) target_link_libraries(linear_conics PRIVATE CGAL::CGAL_Basic_viewer) target_link_libraries(parabolas PRIVATE CGAL::CGAL_Basic_viewer) @@ -32,6 +29,7 @@ if(CGAL_Qt6_FOUND) target_link_libraries(circular_arcs PRIVATE CGAL::CGAL_Basic_viewer) target_link_libraries(rational_functions PRIVATE CGAL::CGAL_Basic_viewer) target_link_libraries(spherical_insert PRIVATE CGAL::CGAL_Basic_viewer) + target_link_libraries(spherical_overlay PRIVATE CGAL::CGAL_Basic_viewer) else() message( STATUS diff --git a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_agas.cpp b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_agas.cpp new file mode 100644 index 00000000000..ec665ebc894 --- /dev/null +++ b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_agas.cpp @@ -0,0 +1,105 @@ +#include + +#include + +#include "arr_geodesic.h" +#include "arr_print.h" + +void draw_face_crossing_boundary() { + Arrangement arr; + const auto& traits = *arr.geometry_traits(); + auto ctr_p = traits.construct_point_2_object(); + auto ctr_cv = traits.construct_curve_2_object(); + + auto p1 = ctr_p(-0.95, 0.32, 0), p2 = ctr_p(-0.87, 0.02, 0.49), p3 = ctr_p(-0.93, -0.36, 0), + p4 = ctr_p(-0.81, -0.03, -0.59); + auto arcs = {ctr_cv(p1, p2), ctr_cv(p2, p3), ctr_cv(p3, p4), ctr_cv(p4, p1)}; + CGAL::insert(arr, arcs.begin(), arcs.end()); + + print_arrangement_size(arr); + CGAL::draw(arr, "face crossing boundary"); +} + +void draw_lakes() { + Arrangement arr; + const auto& traits = *arr.geometry_traits(); + auto ctr_p = traits.construct_point_2_object(); + auto ctr_cv = traits.construct_curve_2_object(); + + auto poly1 = { + ctr_p(-0.27, -0.053, -0.96), ctr_p(-0.76, -0.15, -0.63), ctr_p(-0.98, -0.19, -0.063), ctr_p(-0.98, -0.098, 0.2), + ctr_p(-0.44, -0.18, 0.88), ctr_p(0.39, -0.0049, 0.92), ctr_p(-0.01, 0.39, 0.92), ctr_p(-0.54, 0.66, 0.53), + ctr_p(-0.83, 0.56, 0.025), ctr_p(-0.57, 0.32, -0.75), ctr_p(-0.087, 0.048, -1), ctr_p(-0.048, 0.088, -1), + ctr_p(0.12, -0.14, -0.98), ctr_p(-0.12, -0.14, -0.98), + }; + auto poly2 = {ctr_p(-0.24, -0.53, -0.81), ctr_p(-0.47, -0.54, -0.69), ctr_p(-0.68, -0.65, -0.32), + ctr_p(-0.71, -0.68, 0.2), ctr_p(-0.54, -0.52, 0.67), ctr_p(-0.18, -0.72, 0.67), + ctr_p(0.31, -0.68, 0.67), ctr_p(0.71, -0.69, 0.11), ctr_p(0.6, -0.58, -0.56), + ctr_p(0.21, -0.62, -0.75)}; + auto poly3 = {ctr_p(0.44, 0.27, -0.86), ctr_p(0.58, -0.063, -0.81), ctr_p(0.87, -0.094, -0.48), + ctr_p(0.97, -0.1, 0.2), ctr_p(0.46, 0.77, 0.45), ctr_p(-0.023, 0.89, 0.45), + ctr_p(-0.3, 0.95, 0.11), ctr_p(-0.22, 0.69, -0.69), ctr_p(-0.076, 0.35, -0.93)}; + auto poly4 = { + ctr_p(0.4, 0.67, -0.63), ctr_p(0.78, 0.39, -0.48), ctr_p(0.92, 0.35, -0.15), + ctr_p(0.52, 0.86, 0.025), ctr_p(0.068, 0.99, -0.15), ctr_p(0.22, 0.85, -0.48), + }; + std::vector arcs; + std::vector> polygons{poly1, poly2, poly3, poly4}; + for(const auto& poly : polygons) { + for(size_t i = 0; i < poly.size(); ++i) { + size_t next = (i + 1) % poly.size(); + arcs.push_back(ctr_cv(poly[i], poly[next])); + } + } + + CGAL::insert(arr, arcs.begin(), arcs.end()); + print_arrangement_size(arr); + CGAL::draw(arr, "lakes"); +} + +void draw_guassian_map() { + Arrangement arr; + const auto& traits = *arr.geometry_traits(); + auto ctr_p = traits.construct_point_2_object(); + auto ctr_cv = traits.construct_curve_2_object(); + + auto p1 = ctr_p(1, 1, 1), p2 = ctr_p(-1, -1, 1), p3 = ctr_p(-1, 1, -1), p4 = ctr_p(1, -1, -1); + auto arcs = {ctr_cv(p1, p2), ctr_cv(p2, p3), ctr_cv(p3, p1), ctr_cv(p1, p4), ctr_cv(p2, p4), ctr_cv(p3, p4)}; + CGAL::insert(arr, arcs.begin(), arcs.end()); + print_arrangement_size(arr); + CGAL::draw(arr, "guassian map of a tetrahedron"); +} + +void draw_random_arcs(int n) { + Arrangement arr; + const auto& traits = *arr.geometry_traits(); + auto ctr_p = traits.construct_point_2_object(); + auto ctr_cv = traits.construct_curve_2_object(); + + CGAL::Random random; + std::vector points; + for(int i = 0; i < n; ++i) { + double x = random.get_double(-1.0, 1.0); + double y = random.get_double(-1.0, 1.0); + double z = random.get_double(-1.0, 1.0); + points.push_back(ctr_p(x, y, z)); + } + std::vector curves; + for(int i = 0; i < n; ++i) { + int j = random.get_int(0, n - 1); + if(i == j) j = (j + 1) % n; + curves.push_back(ctr_cv(points[i], points[j])); + } + + CGAL::insert(arr, curves.begin(), curves.end()); + print_arrangement_size(arr); + CGAL::draw(arr, (std::to_string(n) + " random arcs").c_str()); +} + +int main() { + draw_face_crossing_boundary(); + draw_lakes(); + draw_guassian_map(); + draw_random_arcs(100); + return 0; +} diff --git a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr.cpp b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr.cpp index 0476648efa9..a04cff18de3 100644 --- a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr.cpp +++ b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr.cpp @@ -1,15 +1,13 @@ -#include "CGAL/Random.h" +#include + +#include +#include #include #include #include #include #include -using Kernel = CGAL::Exact_predicates_exact_constructions_kernel; -using Traits = CGAL::Arr_non_caching_segment_traits_2; -using Point = Traits::Point_2; -using Arrangement_2 = CGAL::Arrangement_2; - /*! Convert HSV to RGB color space * Converts a given set of HSV values `h', `s', `v' into RGB coordinates. * The output RGB values are in the range [0, 255], and the input HSV values @@ -70,9 +68,14 @@ std::tuple hsv_to_rgb(float hue, fl return std::make_tuple(redc, greenc, bluec); } -int main() { +void draw_rect() { + using Kernel = CGAL::Exact_predicates_exact_constructions_kernel; + using Traits = CGAL::Arr_non_caching_segment_traits_2; + using Point = Traits::Point_2; + using Arrangement = CGAL::Arrangement_2; + Traits traits; - Arrangement_2 arr(&traits); + Arrangement arr(&traits); auto ctr_xcv = traits.construct_x_monotone_curve_2_object(); CGAL::insert(arr, ctr_xcv(Point(-2, -2), Point(2, -2))); @@ -92,12 +95,12 @@ int main() { std::cout << arr.number_of_vertices() << ", " << arr.number_of_edges() << ", " << arr.number_of_faces() << std::endl; - CGAL::Graphics_scene_options + CGAL::Graphics_scene_options gso; - gso.colored_face = [](const Arrangement_2&, Arrangement_2::Face_const_handle) -> bool { return true; }; + gso.colored_face = [](const Arrangement&, Arrangement::Face_const_handle) -> bool { return true; }; - gso.face_color = [](const Arrangement_2& arr, Arrangement_2::Face_const_handle fh) -> CGAL::IO::Color { + gso.face_color = [](const Arrangement& arr, Arrangement::Face_const_handle fh) -> CGAL::IO::Color { CGAL::Random random((size_t(fh.ptr()))); float h = 360.0f * random.get_double(0, 1); float s = 0.5; @@ -106,7 +109,136 @@ int main() { return CGAL::IO::Color(r, g, b); }; - CGAL::draw(arr, gso, "hsv colors"); + CGAL::draw(arr, gso, "rect with hsv colors"); +} +void draw_nested() { + using Kernel = CGAL::Exact_predicates_exact_constructions_kernel; + using Traits = CGAL::Arr_segment_traits_2; + using Point = Traits::Point_2; + using Arrangement = CGAL::Arrangement_2; + using X_monotone_curve = Traits::X_monotone_curve_2; + + Arrangement arr; + auto traits = arr.traits(); + auto ctr_xcv = traits->construct_x_monotone_curve_2_object(); + + std::vector curves; + { + // a hexagon centered at the origin + const double r = 10.0; + for(int i = 0; i < 6; ++i) { + int next = (i + 1) % 6; + Point source(r * cos(i * CGAL_PI / 3), r * sin(i * CGAL_PI / 3)); + Point target(r * cos(next * CGAL_PI / 3), r * sin(next * CGAL_PI / 3)); + curves.push_back(ctr_xcv(source, target)); + } + } + { + // a square inside the hexagon + const double size = 4.0; + Point p1(-size, size), p2(size, size), p3(size, -size), p4(-size, -size); + auto rect = {ctr_xcv(p1, p2), ctr_xcv(p2, p3), ctr_xcv(p3, p4), ctr_xcv(p4, p1)}; + curves.insert(curves.end(), rect.begin(), rect.end()); + } + { + // two adjacent triangle inside the square + Point p1(-1, 0), p2(1, 0), p3(0, sqrt(3)), p4(0, -sqrt(3)); + auto tri1 = {ctr_xcv(p1, p2), ctr_xcv(p2, p3), ctr_xcv(p3, p1)}; + auto tri2 = {ctr_xcv(p1, p2), ctr_xcv(p2, p4), ctr_xcv(p4, p1)}; + curves.insert(curves.end(), tri1.begin(), tri1.end()); + curves.insert(curves.end(), tri2.begin(), tri2.end()); + } + // a degenerate hole inside the square + auto degen_seg = ctr_xcv({-1, -3}, {1, -3}); + curves.push_back(degen_seg); + // an isolated vertex inside the square + auto iso_point = Point{1, -1}; + + CGAL::insert(arr, curves.begin(), curves.end()); + CGAL::insert_point(arr, iso_point); + CGAL::draw(arr, "nested polygons"); +} + +void draw_unbounded_linear_grid() { + using Kernel = CGAL::Exact_predicates_exact_constructions_kernel; + using Traits = CGAL::Arr_linear_traits_2; + using Point = Traits::Point_2; + using Segment = Traits::Segment_2; + using Ray = Traits::Ray_2; + using Line = Traits::Line_2; + using X_monotone_curve = Traits::X_monotone_curve_2; + using Arrangement = CGAL::Arrangement_2; + + Arrangement arr; + auto& traits = *arr.traits(); + + // Insert a n*n grid + int n = 5; + for(int i = 0; i < n; ++i) { + Point p1(i * 5, 0); + Point p2(i * 5, 1); + CGAL::insert(arr, X_monotone_curve(Line(p1, p2))); + } + for(int i = 0; i < n; ++i) { + Point p1(0, i * 5); + Point p2(1, i * 5); + CGAL::insert(arr, X_monotone_curve(Line(p1, p2))); + } + // Generate a inner square(2*2) for all cells + // And an inner triangle for each square + for(int i = 0; i < n; ++i) { + for(int j = 0; j < n; ++j) { + Point p1(i * 5 + 1, j * 5 + 1); + Point p2(i * 5 + 4, j * 5 + 4); + CGAL::insert(arr, X_monotone_curve(Segment(p1, Point(p2.x(), p1.y())))); + CGAL::insert(arr, X_monotone_curve(Segment(Point(p1.x(), p2.y()), p2))); + CGAL::insert(arr, X_monotone_curve(Segment(p1, Point(p1.x(), p2.y())))); + CGAL::insert(arr, X_monotone_curve(Segment(Point(p2.x(), p1.y()), p2))); + // Insert a triangle inside the square + Point tri_p1(i * 5 + 2, j * 5 + 2); + Point tri_p2(i * 5 + 3, j * 5 + 2); + Point tri_p3(i * 5 + 2.5, j * 5 + 3); + CGAL::insert(arr, X_monotone_curve(Segment(tri_p1, tri_p2))); + CGAL::insert(arr, X_monotone_curve(Segment(tri_p2, tri_p3))); + CGAL::insert(arr, X_monotone_curve(Segment(tri_p3, tri_p1))); + // Connect the triangle to the square + Point top(i * 5 + 2.5, j * 5 + 4); + CGAL::insert(arr, X_monotone_curve(Segment(tri_p1, top))); + } + } + + CGAL::draw(arr, "unbounded linear grid"); +} + +void draw_random_segments(int n) { + using Kernel = CGAL::Exact_predicates_exact_constructions_kernel; + using Traits = CGAL::Arr_segment_traits_2; + using Point = Traits::Point_2; + using Arrangement = CGAL::Arrangement_2; + using X_monotone_curve = Traits::X_monotone_curve_2; + + Arrangement arr; + auto traits = arr.traits(); + auto ctr_xcv = traits->construct_x_monotone_curve_2_object(); + CGAL::Random random; + + std::vector curves; + for(int i = 0; i < n; ++i) { + double x1 = random.get_double(-100, 100); + double y1 = random.get_double(-100, 100); + double x2 = random.get_double(-100, 100); + double y2 = random.get_double(-100, 100); + curves.push_back(ctr_xcv(Point(x1, y1), Point(x2, y2))); + } + CGAL::insert(arr, curves.begin(), curves.end()); + CGAL::draw(arr, (std::to_string(n) + " random segments").c_str()); +} + +int main() { + draw_rect(); + draw_nested(); + draw_unbounded_linear_grid(); + draw_random_segments(100); return EXIT_SUCCESS; } diff --git a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_degen_holes.cpp b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_degen_holes.cpp deleted file mode 100644 index cb12cc07947..00000000000 --- a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_degen_holes.cpp +++ /dev/null @@ -1,44 +0,0 @@ -#include -#include -#include - -using Kernel = CGAL::Exact_predicates_exact_constructions_kernel; -using Traits = CGAL::Arr_segment_traits_2; -using Point = Traits::Point_2; -using Arrangement = CGAL::Arrangement_2; -using X_monotone_curve = Traits::X_monotone_curve_2; - -int main() { - Arrangement arr; - auto traits = arr.traits(); - auto cst_x_curve = traits->construct_x_monotone_curve_2_object(); - - // make a hexagon centered at the origin with radius 10 - double radius = 10.0; - std::array hexagon_points = {{ - {radius * cos(0 * CGAL_PI / 3), radius * sin(0 * CGAL_PI / 3)}, // 0° - {radius * cos(1 * CGAL_PI / 3), radius * sin(1 * CGAL_PI / 3)}, // 60° - {radius * cos(2 * CGAL_PI / 3), radius * sin(2 * CGAL_PI / 3)}, // 120° - {radius * cos(3 * CGAL_PI / 3), radius * sin(3 * CGAL_PI / 3)}, // 180° - {radius * cos(4 * CGAL_PI / 3), radius * sin(4 * CGAL_PI / 3)}, // 240° - {radius * cos(5 * CGAL_PI / 3), radius * sin(5 * CGAL_PI / 3)}, // 300° - }}; - std::array hexagon; - for(size_t i = 0; i < hexagon_points.size(); ++i) { - size_t next_i = (i + 1) % hexagon_points.size(); - hexagon[i] = cst_x_curve(hexagon_points[i], hexagon_points[next_i]); - } - // rect hole - auto hole_rectangle = {cst_x_curve({-2, -2}, {2, -2}), cst_x_curve({2, -2}, {2, 2}), cst_x_curve({2, 2}, {-2, 2}), - cst_x_curve({-2, 2}, {-2, -2})}; - // iso vertex inside rect hole - auto iso_vertex_inside_hole = Point{0.5, 0.5}; - // degenerate segment below the rect hole - auto degenerate_segment = cst_x_curve({0, -3}, {1, -3}); - - CGAL::insert_point(arr, iso_vertex_inside_hole); - CGAL::insert(arr, hexagon.begin(), hexagon.end()); - CGAL::insert(arr, hole_rectangle.begin(), hole_rectangle.end()); - CGAL::insert(arr, degenerate_segment); - CGAL::draw(arr); -} \ No newline at end of file diff --git a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_spherical_lakes.cpp b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_spherical_lakes.cpp deleted file mode 100644 index 32a67ad63a9..00000000000 --- a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_spherical_lakes.cpp +++ /dev/null @@ -1,35 +0,0 @@ -#include -#include -#include -#include -#include - -#include "arr_geodesic.h" -#include "arr_print.h" - -int main() { - Geom_traits traits; - auto ctr_p = traits.construct_point_2_object(); - auto ctr_cv = traits.construct_curve_2_object(); - Arrangement arr(&traits); - // identification curve in two parts - CGAL::insert(arr, ctr_cv(ctr_p(0, 0, -1), ctr_p(-1, 0, 0))); - CGAL::insert(arr, ctr_cv(ctr_p(-1, 0, 0), ctr_p(0, 0, 1))); - - { - // outer face - Point p1 = ctr_p(0, 1, 0), p2 = ctr_p(0, 0, 1), p3 = ctr_p(-1, 0, 0); - auto arcs = {ctr_cv(p1, p2), ctr_cv(p1, p3), ctr_cv(p2, p3)}; - CGAL::insert(arr, arcs.begin(), arcs.end()); - } - { - // lake - Point p1 = ctr_p(0.45, -0.45, -0.75), p2 = ctr_p(0.45, -0.85, -0.30), p3 = ctr_p(0.85, -0.35, -0.35); - auto arcs = {ctr_cv(p1, p2), ctr_cv(p2, p3), ctr_cv(p3, p1)}; - CGAL::insert(arr, arcs.begin(), arcs.end()); - } - - print_arrangement_size(arr); - CGAL::draw(arr, "spherical lakes"); - return 0; -} diff --git a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_spherical_polar_faces.cpp b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_spherical_polar_faces.cpp deleted file mode 100644 index 97049570215..00000000000 --- a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/draw_arr_spherical_polar_faces.cpp +++ /dev/null @@ -1,28 +0,0 @@ -#include -#include -#include -#include -#include - -#include "arr_geodesic.h" -#include "arr_print.h" - -int main() { - Geom_traits traits; - auto ctr_p = traits.construct_point_2_object(); - auto ctr_cv = traits.construct_curve_2_object(); - Arrangement arr(&traits); - // identification curve in two parts - CGAL::insert(arr, ctr_cv(ctr_p(0, 0, -1), ctr_p(-1, 0, 0))); - CGAL::insert(arr, ctr_cv(ctr_p(-1, 0, 0), ctr_p(0, 0, 1))); - - CGAL::insert(arr, ctr_cv(ctr_p(0, 0, 1), ctr_p(-1, 0, 0))); - Point p1 = ctr_p(1, 0, CGAL_PI / 4.0), p2 = ctr_p(0, 1, CGAL_PI / 4.0), p3 = ctr_p(-1, 0, CGAL_PI / 4.0), - p4 = ctr_p(0, -1, CGAL_PI / 4.0); - Curve arcs[] = {ctr_cv(p1, p2), ctr_cv(p2, p3), ctr_cv(p3, p4), ctr_cv(p4, p1)}; - CGAL::insert(arr, arcs, arcs + sizeof(arcs) / sizeof(Curve)); - - print_arrangement_size(arr); - CGAL::draw(arr, "spherical faces containing one of the pole"); - return 0; -} diff --git a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/spherical_insert.cpp b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/spherical_insert.cpp index cdce7e358c9..ad54d5890b6 100644 --- a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/spherical_insert.cpp +++ b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/spherical_insert.cpp @@ -1,7 +1,5 @@ //! \file examples/Arrangement_on_surface_2/spherical_insert.cpp // Constructing an arrangement of arcs of great circles. -#define CGAL_DRAW_AOS_DEBUG -#define CGAL_DRAW_AOS_TRIANGULATOR_DEBUG_FILE_DIR "/Users/shep/Codes/aos_2_js_helper" #include #include #include diff --git a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/spherical_overlay.cpp b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/spherical_overlay.cpp index e43b1df9ab8..3ff03015e93 100644 --- a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/spherical_overlay.cpp +++ b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/spherical_overlay.cpp @@ -3,6 +3,7 @@ #include #include +#include #include "arr_geodesic_on_sphere.h" @@ -34,8 +35,8 @@ int main() { Arrangement overlay_arr; Overlay_traits overlay_traits; overlay(arr1, arr2, overlay_arr, overlay_traits); - std::cout << "No. of vertices: " << overlay_arr.number_of_vertices() - << std::endl; + std::cout << "No. of vertices: " << overlay_arr.number_of_vertices() << std::endl; + CGAL::draw(overlay_arr, "spherical overlay"); return 0; } diff --git a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/unbounded_linear.cpp b/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/unbounded_linear.cpp deleted file mode 100644 index 55c281e5132..00000000000 --- a/Arrangement_on_surface_2/examples/Arrangement_on_surface_2/unbounded_linear.cpp +++ /dev/null @@ -1,51 +0,0 @@ -#include "arr_linear.h" -#include "arr_print.h" - -#include "CGAL/draw_arrangement_2.h" - -int main() { - Arrangement arr; - auto& traits = *arr.traits(); - - // Insert a n*n grid, each cell is a square of size 5 - int n = 5; - for(int i = 0; i < n; ++i) { - Point p1(i * 5, 0); - Point p2(i * 5, 1); - CGAL::insert(arr, X_monotone_curve(Line(p1, p2))); - } - - for(int i = 0; i < n; ++i) { - Point p1(0, i * 5); - Point p2(1, i * 5); - CGAL::insert(arr, X_monotone_curve(Line(p1, p2))); - } - - Vertex_const_handle vh; - // Generate a inner square(2*2) for all cells - // And an inner triangle for each square - for(int i = 0; i < n; ++i) { - for(int j = 0; j < n; ++j) { - Point p1(i * 5 + 1, j * 5 + 1); - Point p2(i * 5 + 4, j * 5 + 4); - CGAL::insert(arr, X_monotone_curve(Segment(p1, Point(p2.x(), p1.y())))); - CGAL::insert(arr, X_monotone_curve(Segment(Point(p1.x(), p2.y()), p2))); - CGAL::insert(arr, X_monotone_curve(Segment(p1, Point(p1.x(), p2.y())))); - CGAL::insert(arr, X_monotone_curve(Segment(Point(p2.x(), p1.y()), p2))); - // Insert a triangle inside the square - Point tri_p1(i * 5 + 2, j * 5 + 2); - Point tri_p2(i * 5 + 3, j * 5 + 2); - Point tri_p3(i * 5 + 2.5, j * 5 + 3); - CGAL::insert(arr, X_monotone_curve(Segment(tri_p1, tri_p2))); - CGAL::insert(arr, X_monotone_curve(Segment(tri_p2, tri_p3))); - CGAL::insert(arr, X_monotone_curve(Segment(tri_p3, tri_p1))); - // Connect the triangle to the square - Point top(i * 5 + 2.5, j * 5 + 4); - CGAL::insert(arr, X_monotone_curve(Segment(tri_p1, top))); - } - } - - print_arrangement_size(arr); - - CGAL::draw(arr); -} \ No newline at end of file diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_geodesic_arc_on_sphere_traits_2.h b/Arrangement_on_surface_2/include/CGAL/Arr_geodesic_arc_on_sphere_traits_2.h index b893940543d..e16e0823660 100644 --- a/Arrangement_on_surface_2/include/CGAL/Arr_geodesic_arc_on_sphere_traits_2.h +++ b/Arrangement_on_surface_2/include/CGAL/Arr_geodesic_arc_on_sphere_traits_2.h @@ -35,6 +35,8 @@ #include #include #include +#include "CGAL/Kernel/global_functions_3.h" +#include "CGAL/enum.h" namespace CGAL { @@ -362,6 +364,11 @@ public: */ inline Comparison_result compare_x(const Direction_3& d1, const Direction_3& d2) const { + bool is_d1_on_boundary = Direction_3(CGAL::cross_product(d1.vector(), Direction_3(0, 0, 1).vector())) == identification_normal(); + bool is_d2_on_boundary = Direction_3(CGAL::cross_product(d2.vector(), Direction_3(0, 0, 1).vector())) == identification_normal(); + if(is_d1_on_boundary && is_d2_on_boundary) return EQUAL; + if(is_d1_on_boundary) return SMALLER; + if(is_d2_on_boundary) return LARGER; // Compare the projections onto the xy plane: Direction_2 d1_2 = project_xy(d1); Direction_2 d2_2 = project_xy(d2); @@ -397,8 +404,7 @@ public: const Point_2& point) const { using Traits = Arr_geodesic_arc_on_sphere_traits_2; - CGAL_precondition(!point.is_min_boundary()); - CGAL_precondition(!point.is_max_boundary()); + if(point.is_max_boundary() || point.is_min_boundary()) return true; Direction_2 p = project_xy(point); if (xcv.is_vertical()) { @@ -503,7 +509,7 @@ public: Direction_3& d(p); d = Direction_3(x, y, z); init(p, std::integral_constant()); - return p; + return p; } /*! constructs a point on the sphere from a (not necessarily normalized) @@ -1035,9 +1041,6 @@ public: * \pre p2 does not lie on the boundary. */ Comparison_result operator()(const Point_2& p1, const Point_2& p2) const { - if(!p1.is_no_boundary() && !p2.is_no_boundary()) return EQUAL; - if(!p1.is_no_boundary()) return SMALLER; - if(p2.is_no_boundary()) return LARGER; return m_traits.compare_x(p1, p2); } }; @@ -1093,9 +1096,6 @@ public: * \pre p2 does not lie on the boundary. */ Comparison_result operator()(const Point_2& p1, const Point_2& p2) const { - CGAL_precondition(p1.is_no_boundary()); - CGAL_precondition(p2.is_no_boundary()); - return m_traits.compare_xy(p1, p2); } }; @@ -1179,7 +1179,6 @@ public: */ Comparison_result operator()(const Point_2& p, const X_monotone_curve_2& xc) const { - CGAL_precondition(!p.is_min_boundary() && !p.is_max_boundary()); CGAL_precondition(m_traits.is_in_x_range(xc, p)); if (xc.is_vertical()) { @@ -1653,10 +1652,8 @@ public: const X_monotone_curve_2& xcv, Arr_curve_end CGAL_precondition_code(ce)) const { - CGAL_precondition(point.is_no_boundary()); CGAL_precondition_code (const Point_2& p2 = (ce == ARR_MIN_END) ? xcv.left() : xcv.right();); - CGAL_precondition(!p2.is_no_boundary()); CGAL_precondition(xcv.is_vertical()); CGAL_precondition(!m_traits.is_on_y_identification_2_object()(xcv)); diff --git a/Arrangement_on_surface_2/include/CGAL/Arr_linear_traits_2.h b/Arrangement_on_surface_2/include/CGAL/Arr_linear_traits_2.h index 36ca18b980b..609abf0ec2b 100644 --- a/Arrangement_on_surface_2/include/CGAL/Arr_linear_traits_2.h +++ b/Arrangement_on_surface_2/include/CGAL/Arr_linear_traits_2.h @@ -34,6 +34,7 @@ #include #include #include +#include "CGAL/number_utils.h" namespace CGAL { @@ -1560,6 +1561,7 @@ public: template OutputIterator operator()(const X_monotone_curve_2& xcv, double /* error */, OutputIterator oi, bool l2r = true) const { + if(xcv.is_ray() || xcv.is_line()) return oi; auto min_vertex = m_traits.construct_min_vertex_2_object(); auto max_vertex = m_traits.construct_max_vertex_2_object(); const auto& src = (l2r) ? min_vertex(xcv) : max_vertex(xcv); @@ -1582,24 +1584,25 @@ public: { using Approx_pnt = Approximate_point_2; using Approx_seg = Approximate_kernel::Segment_2; + using Approx_ray = Approximate_kernel::Ray_2; + using Approx_lin = Approximate_kernel::Line_2; auto xmin = bbox.xmin(); auto ymin = bbox.ymin(); auto xmax = bbox.xmax(); auto ymax = bbox.ymax(); - Approximate_kernel::Iso_rectangle_2 rect(xmin, ymin, xmax, ymax); - auto xs = CGAL::to_double(xcv.source().x()); - auto ys = CGAL::to_double(xcv.source().y()); - auto xt = CGAL::to_double(xcv.target().x()); - auto yt = CGAL::to_double(xcv.target().y()); + Approximate_kernel::Iso_rectangle_2 rect(xmin, ymin, xmax, ymax); if (xcv.is_ray()) { - using Approx_ray = Approximate_kernel::Ray_2; - Approx_ray ray(Approx_pnt(xs, ys), Approx_pnt(xt, yt)); - const auto result = CGAL::intersection(rect, ray); + auto ray = xcv.ray(); + Kernel kernel; + auto construct_vertex = kernel.construct_point_on_2_object(); + Approx_pnt s = this->operator()(construct_vertex(ray, 0)); + Approx_pnt t = this->operator()(construct_vertex(ray, 1)); + const auto result = CGAL::intersection(rect, Approx_ray(s, t)); if (! result) return oi; if (const auto* res_seg = std::get_if(&*result)) { - *oi++ = res_seg->source(); - *oi++ = res_seg->target(); + *oi++ = l2r ? res_seg->min() : res_seg->max(); + *oi++ = l2r ? res_seg->max() : res_seg->min(); return oi; } const auto* res_pnt = std::get_if(&*result); @@ -1608,14 +1611,17 @@ public: return oi; } if (xcv.is_line()) { - using Approx_lin = Approximate_kernel::Line_2; - Approx_lin line(Approx_pnt(xs, ys), Approx_pnt(xt, yt)); - const auto result = CGAL::intersection(rect, line); + const Line_2 & supp_line = xcv.supp_line(); + Approx_lin approx_supp_line( + CGAL::to_double(supp_line.a()), + CGAL::to_double(supp_line.b()), + CGAL::to_double(supp_line.c())); + const auto result = CGAL::intersection(rect, approx_supp_line); if (! result) return oi; if (const auto* res_seg = std::get_if(&*result)) { - *oi++ = res_seg->source(); - *oi++ = res_seg->target(); + *oi++ = l2r ? res_seg->min() : res_seg->max(); + *oi++ = l2r ? res_seg->max() : res_seg->min(); return oi; } const auto* res_pnt = std::get_if(&*result); @@ -1623,13 +1629,13 @@ public: *oi++ = *res_pnt; return oi; } - Approx_seg seg(Approx_pnt(xs, ys), Approx_pnt(xt, yt)); + Approx_seg seg(this->operator()(xcv.source()), this->operator()(xcv.target())); const auto result = CGAL::intersection(rect, seg); if (! result) return oi; if (const auto* res_seg = std::get_if(&*result)) { - *oi++ = res_seg->source(); - *oi++ = res_seg->target(); + *oi++ = l2r ? res_seg->min() : res_seg->max(); + *oi++ = l2r ? res_seg->max() : res_seg->min(); return oi; } const auto* res_pnt = std::get_if(&*result); diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_approximation_cache.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_approximation_cache.h index 28bf64e005a..3f6a43ab402 100644 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_approximation_cache.h +++ b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_approximation_cache.h @@ -1,7 +1,20 @@ +// Copyright (c) 2025 +// Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), +// INRIA Sophia-Antipolis (France), +// Max-Planck-Institute Saarbruecken (Germany), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s): Shepard Liu + #ifndef CGAL_DRAW_AOS_ARR_APPROXIMATION_CACHE_H #define CGAL_DRAW_AOS_ARR_APPROXIMATION_CACHE_H -#include - #include #include @@ -24,87 +37,27 @@ class Arr_approximation_cache using Geom_traits = typename Arrangement::Geometry_traits_2; using Approx_traits = Arr_approximate_traits; -public: - using Vertex_cache_obj = typename Approx_traits::Point_geom; - using Halfedge_cache_obj = typename Approx_traits::Polyline_geom; - using Face_cache_obj = typename Approx_traits::Triangulated_face; + using Vertex_cache_obj = typename Approx_traits::Point; + using Halfedge_cache_obj = typename Approx_traits::Polyline; + using Face_cache_obj = typename Approx_traits::Triangle_soup; -private: using Vertex_const_handle = typename Arrangement::Vertex_const_handle; - using Edge_const_handle = typename Arrangement::Edge_const_iterator; using Halfedge_const_handle = typename Arrangement::Halfedge_const_iterator; using Face_const_handle = typename Arrangement::Face_const_handle; using Vertex_cache = unordered_flat_map; using Halfedge_cache = unordered_flat_map; using Face_cache = unordered_flat_map; - using Vertex_cache_const_iterator = typename Vertex_cache::const_iterator; - using Halfedge_cache_const_iterator = typename Halfedge_cache::const_iterator; - using Face_cache_const_iterator = typename Face_cache::const_iterator; - using Vertex_cache_range = boost::iterator_range; - using Halfedge_cache_range = boost::iterator_range; - using Face_cache_range = boost::iterator_range; - public: Arr_approximation_cache() = default; - void reserve_vertices(std::size_t size) { m_vertices.reserve(size); } - void reserve_halfedges(std::size_t size) { m_halfedges.reserve(size); } - void reserve_faces(std::size_t size) { m_faces.reserve(size); } + const Vertex_cache& vertices() const { return m_vertices; } + const Halfedge_cache& halfedges() const { return m_halfedges; } + const Face_cache& faces() const { return m_faces; } - std::pair try_emplace(const Vertex_const_handle& vh) { - const auto& [it, inserted] = m_vertices.try_emplace(vh, Vertex_cache_obj()); - return {it->second, inserted}; - } - std::pair try_emplace(const Halfedge_const_handle& he) { - const auto& [it, inserted] = m_halfedges.try_emplace(he, Halfedge_cache_obj()); - return {it->second, inserted}; - } - std::pair try_emplace(const Face_const_handle& fh) { - const auto& [it, inserted] = m_faces.try_emplace(fh, Face_cache_obj()); - return {it->second, inserted}; - } - - std::pair get(const Vertex_const_handle& vh) const { - auto it = m_vertices.find(vh); - if(it != m_vertices.end()) { - return {it->second, true}; - } - return {Vertex_cache_obj(), false}; - } - std::pair get(const Halfedge_const_handle& he) const { - auto it = m_halfedges.find(he); - if(it != m_halfedges.end()) { - return {it->second, true}; - } - return {Halfedge_cache_obj(), false}; - } - std::pair get(const Face_const_handle& fh) const { - auto it = m_faces.find(fh); - if(it != m_faces.end()) { - return {it->second, true}; - } - return {Face_cache_obj(), false}; - } - - bool has(const Vertex_const_handle& vh) const { return m_vertices.find(vh) != m_vertices.end(); } - bool has(const Halfedge_const_handle& he) const { return m_halfedges.find(he) != m_halfedges.end(); } - bool has(const Face_const_handle& fh) const { return m_faces.find(fh) != m_faces.end(); } - - Vertex_cache_const_iterator vertex_begin() const { return m_vertices.begin(); } - Vertex_cache_const_iterator vertex_end() const { return m_vertices.end(); } - Vertex_cache_range vertices() const { return boost::make_iterator_range(vertex_begin(), vertex_end()); } - std::size_t number_of_vertices() const { return m_vertices.size(); } - - Halfedge_cache_const_iterator halfedge_begin() const { return m_halfedges.begin(); } - Halfedge_cache_const_iterator halfedge_end() const { return m_halfedges.end(); } - Halfedge_cache_range halfedges() const { return boost::make_iterator_range(halfedge_begin(), halfedge_end()); } - std::size_t number_of_halfedges() const { return m_halfedges.size(); } - - Face_cache_const_iterator face_begin() const { return m_faces.begin(); } - Face_cache_const_iterator face_end() const { return m_faces.end(); } - Face_cache_range faces() const { return boost::make_iterator_range(face_begin(), face_end()); } - std::size_t number_of_faces() const { return m_faces.size(); } + Vertex_cache& vertices() { return m_vertices; } + Halfedge_cache& halfedges() { return m_halfedges; } + Face_cache& faces() { return m_faces; } private: Vertex_cache m_vertices; diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_face.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_face.h index 4052a6c17cb..af14acbdfa8 100644 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_face.h +++ b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_face.h @@ -1,12 +1,23 @@ +// Copyright (c) 2025 +// Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), +// INRIA Sophia-Antipolis (France), +// Max-Planck-Institute Saarbruecken (Germany), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s): Shepard Liu + #ifndef CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_FACE_H #define CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_FACE_H - -#include -#include +#include #include -#include #include -#include #include #include @@ -15,264 +26,209 @@ #include #include #include -#include #include #include #include -#include "CGAL/basic.h" -#include "CGAL/enum.h" -#include "CGAL/number_utils.h" namespace CGAL { namespace draw_aos { -/** - * @brief Bounded face approximation for arrangements. - * @note Member functions are not thread-safe. +/*! + * \brief Functor to approximate arrangement face with triangles within a bounding box. + * + * \tparam Arrangement */ template class Arr_bounded_approximate_face { - using Geom_traits = typename Arrangement::Geometry_traits_2; - using Approx_traits = Arr_approximate_traits; using Face_const_handle = typename Arrangement::Face_const_handle; using Halfedge_const_handle = typename Arrangement::Halfedge_const_handle; using Vertex_const_handle = typename Arrangement::Vertex_const_handle; - using Point_geom = typename Approx_traits::Point_geom; - using Polyline_geom = typename Approx_traits::Polyline_geom; - using Point_2 = typename Geom_traits::Point_2; using Ccb_halfedge_const_circulator = typename Arrangement::Ccb_halfedge_const_circulator; - using Approx_nt = typename Approx_traits::Approx_nt; - using Approx_kernel = typename Approx_traits::Approx_kernel; - using Approx_line_2 = typename Approx_kernel::Line_2; - using Triangulated_face = typename Approx_traits::Triangulated_face; - using Portal_exit = typename Arr_portals::Portal_exit; - using Portal_exit_vector = typename Arr_portals::Portal_exit_vector; + + using Geom_traits = typename Arrangement::Geometry_traits_2; + using Approx_traits = Arr_approximate_traits; + using Point = typename Approx_traits::Point; + using Polyline = typename Approx_traits::Polyline; + using Triangle_soup = typename Approx_traits::Triangle_soup; using Bounded_approximate_vertex = Arr_bounded_approximate_vertex; using Bounded_approximate_halfedge = Arr_bounded_approximate_halfedge; using Bounded_render_context = Arr_bounded_render_context; - using Triangulator = Arr_bounded_face_triangulator; + static constexpr bool Is_on_curved_surface = is_or_derived_from_curved_surf_traits_v; + struct Left_to_right_tag {}; struct Right_to_left_tag {}; - struct Inner_ccb_tag - {}; - struct Outer_ccb_tag - {}; - private: + /*! + * \brief A stateful geometry simplifier that simplifies horizontal and vertical segments + * + * \tparam OutputIterator + */ + template + class Colinear_simplifier + { + public: + Colinear_simplifier(OutputIterator out_it) + : m_out_it(out_it) {} + + void dump() { + if(m_start.has_value()) { + *m_out_it++ = m_start.value(); + m_start.reset(); + } + if(m_mid.has_value()) { + *m_out_it++ = m_mid.value(); + m_mid.reset(); + } + } + + void push_back(Point p) { + if(m_mid.has_value()) { + if(p.y() == m_mid->y() && p.y() == m_start->y() || p.x() == m_mid->x() && p.x() == m_start->x()) + // Three points are collinear horizontally or vertically. + m_mid = p; + else { + *m_out_it++ = m_start.value(); + m_start = m_mid; + m_mid = p; + } + return; + } + + if(m_start.has_value()) + m_mid = p; + else + m_start = p; + } + + ~Colinear_simplifier() { dump(); } + + private: + OutputIterator m_out_it; + std::optional m_start, m_mid; + }; + class Context : public Bounded_render_context { - private: - using Output_iterator = typename Triangulator::Insert_iterator; + using Simplifier = Colinear_simplifier>; public: - using Insert_iterator = boost::function_output_iterator>; - - Context(const Bounded_render_context& ctx, - Output_iterator out_it, - const Bounded_approximate_vertex& bounded_approx_vertex, - const Bounded_approximate_halfedge& bounded_approx_halfedge) + Context(const Bounded_render_context& ctx, Triangulator& triangulator) : Bounded_render_context(ctx) - , m_base_out_it(out_it) - , m_approx(ctx.m_traits.approximate_2_object()) - , m_proj(ctx.m_traits) - , m_bounded_approx_vertex(bounded_approx_vertex) - , m_bounded_approx_halfedge(bounded_approx_halfedge) { - m_out_it = boost::make_function_output_iterator(std::function([this](Point_geom pt) { - if(!this->contains_x(pt.x())) return; - pt = Point_geom(pt.x(), std::clamp(pt.y(), this->ymin(), this->ymax())); - m_last_pt = pt; - *m_base_out_it++ = pt; - })); + , m_triangulator(triangulator) { + if constexpr(!Is_on_curved_surface) m_simplifier.emplace(std::back_inserter(m_triangulator)); } // Let's not accidentally copy this object. Context(const Context&) = delete; Context& operator=(const Context&) = delete; - Point_geom approx(const Point_2& pt) const { return m_proj.project(m_approx(pt)); } + void insert(Point pt) { + if(Approx_traits::is_null(pt) || pt == m_last_pt) return; + pt = Point(pt.x(), std::clamp(pt.y(), this->ymin(), this->ymax())); + if constexpr(!Is_on_curved_surface) { + m_simplifier->push_back(pt); + return; + } + m_triangulator.push_back(pt); + m_last_pt = pt; + } + + void start_ccb() { m_triangulator.start_constraint(); } + + void end_ccb() { + if constexpr(!Is_on_curved_surface) m_simplifier->dump(); + m_triangulator.end_constraint(); + } + + const std::optional& last_pt() const { return m_last_pt; } private: - Output_iterator m_base_out_it; - const typename Geom_traits::Approximate_2 m_approx; - const Arr_projector m_proj; - - public: - const Bounded_approximate_vertex& m_bounded_approx_vertex; - const Bounded_approximate_halfedge& m_bounded_approx_halfedge; - std::optional m_last_pt; - Insert_iterator m_out_it; + Triangulator& m_triangulator; + // Colinear simplifier is only used for optimizing planar arrangements. + std::optional m_simplifier; + std::optional m_last_pt; }; private: - /** - * @brief Get the first halfedge on the inner ccb that should be traversed. - * - * @precondition: The vertex must not be isolated. - * @param vh - */ - static Halfedge_const_handle get_inner_ccb_first_halfedge(const Vertex_const_handle& vh) { - auto circ = vh->incident_halfedges(); - if(vh->degree() == 1) return circ->twin(); - - auto curr = circ; - auto next = curr; - ++next; - // Traverse the halfedges incident to the vertex in clockwise direction. - // Note that circulator around a vertex points to a halfedge whose target is the vertex. - do { - // If "curr" crosses 12 o'clock and reaches "next", we found the wanted halfedge. - if(curr->direction() == ARR_LEFT_TO_RIGHT && next->direction() == ARR_RIGHT_TO_LEFT) break; - ++curr; - ++next; - } while(curr != circ); - // The opposite halfedge is the one we want. - return next->twin(); - } - - static Arr_parameter_space side_of_fictitious_edge(const Halfedge_const_handle& he) { + static Arr_parameter_space side_of_fict_edge(const Halfedge_const_handle& he) { const auto& source = he->source(); const auto& target = he->target(); - Arr_parameter_space src_x_space = source->parameter_space_in_x(); - Arr_parameter_space src_y_space = source->parameter_space_in_y(); - Arr_parameter_space tgt_x_space = target->parameter_space_in_x(); - Arr_parameter_space tgt_y_space = target->parameter_space_in_y(); - if(src_x_space == tgt_x_space && src_x_space != ARR_INTERIOR) return src_x_space; - if(src_y_space == tgt_y_space && src_y_space != ARR_INTERIOR) return src_y_space; - CGAL_assertion(false && "Unexpected parameter space for fictitious edge vertices."); + auto sx = source->parameter_space_in_x(); + auto sy = source->parameter_space_in_y(); + auto tx = target->parameter_space_in_x(); + auto ty = target->parameter_space_in_y(); + if(sx == tx && sx != ARR_INTERIOR) return sx; + if(sy == ty && sy != ARR_INTERIOR) return sy; + CGAL_assertion(false && "Unexpected parameter space for fictitious edge ends."); return ARR_INTERIOR; } - // Generate dummy curve(directed left to right) for the fictitious edge. - static Polyline_geom approximate_fictitious_edge(const Context& ctx, const Halfedge_const_handle& he) { - Arr_parameter_space side = side_of_fictitious_edge(he); + // Generate dummy segment for fictitious edge he at its corresponding boundary. + static Polyline approximate_fict_edge(const Context& ctx, const Halfedge_const_handle& he) { + auto side = side_of_fict_edge(he); // There's no need to handle fictitious edges on left or right boundaries. - if(side == ARR_LEFT_BOUNDARY || side == ARR_RIGHT_BOUNDARY) return Polyline_geom{}; - if(side == ARR_BOTTOM_BOUNDARY) { - Point_geom from_pt(ctx.m_last_pt.has_value() ? ctx.m_last_pt->x() : ctx.xmin(), ctx.ymin()); - Point_geom to_pt(ctx.xmax(), ctx.ymin()); - return Polyline_geom{from_pt, to_pt}; - } - if(side == ARR_TOP_BOUNDARY) { - Point_geom from_pt(ctx.xmin(), ctx.ymax()); - Point_geom to_pt(ctx.m_last_pt.has_value() ? ctx.m_last_pt->x() : ctx.xmax(), ctx.ymax()); - return Polyline_geom{from_pt, to_pt}; - } + if(side == ARR_LEFT_BOUNDARY || side == ARR_RIGHT_BOUNDARY) return Polyline{}; + if(side == ARR_BOTTOM_BOUNDARY) return Polyline{ctx.bottom_left(), ctx.bottom_right()}; + if(side == ARR_TOP_BOUNDARY) return Polyline{ctx.top_right(), ctx.top_left()}; CGAL_assertion(false && "Unexpected side for a fictitious edge."); - return Polyline_geom{}; + return Polyline{}; } - static void approximate_vertex(Context& ctx, const Vertex_const_handle& vh) { + void approximate_vertex(Context& ctx, const Vertex_const_handle& vh) const { if(vh->is_at_open_boundary()) return; - ctx.m_bounded_approx_vertex(vh); - // Handle portal from this vertex if there is one. - auto portals_it = ctx.m_portals_map.find(vh); - if(portals_it == ctx.m_portals_map.end()) return; - const auto& portals = portals_it->second; - if(portals.empty()) return; - CGAL_assertion_msg(portals.size() == 1, "Vertex should have at most one portal."); - traverse_portal(ctx, ctx.approx(vh->point()), portals.front()); + m_bounded_approx_vertex(vh); } - template - static void approximate_halfedge(Context& ctx, const Halfedge_const_handle& he, DirectionTag dir_tag) { - const Polyline_geom& polyline = - he->is_fictitious() ? approximate_fictitious_edge(ctx, he) : ctx.m_bounded_approx_halfedge(he); - auto portal_map_it = ctx.m_portals_map.find(he); - const Portal_exit_vector& portals = - portal_map_it != ctx.m_portals_map.end() ? portal_map_it->second : Portal_exit_vector{}; - - // A variant of two pointers algorithm. std::merge() can't fit in our purpose so we implement it ourselves. - auto compare_x = [](const Point_geom& a, const Point_geom& b) -> Comparison_result { - return std::is_same_v ? sign(a.x() - b.x()) : sign(b.x() - a.x()); - }; - std::optional prev_pt; - auto portals_iter = portals.begin(), portals_end = portals.end(); - for(const auto& curr_pt : polyline) { - for(; portals_iter != portals_end; ++portals_iter) { - const Portal_exit& exit = *portals_iter; - Point_geom exit_pt = ctx.approx(exit->point()); - if(compare_x(exit_pt, curr_pt) >= 0 || !prev_pt.has_value()) break; - if(compare_x(exit_pt, *prev_pt) < 0) { - traverse_portal(ctx, exit_pt, exit); - continue; - } - // Compute the portal entry point. - std::optional> res = CGAL::intersection( - Approx_line_2(*prev_pt, curr_pt), Approx_line_2(exit_pt, Point_geom(exit_pt.x(), exit_pt.y() + 1))); - CGAL_assertion(res.has_value() && std::holds_alternative(*res)); - traverse_portal(ctx, std::get(*res), exit); - } - - *ctx.m_out_it++ = curr_pt; - prev_pt = curr_pt; - } - - for(; portals_iter != portals_end; ++portals_iter) { - const Portal_exit& exit = *portals_iter; - traverse_portal(ctx, ctx.approx(exit->point()), exit); - } + void approximate_halfedge(Context& ctx, const Halfedge_const_handle& he) const { + const Polyline& polyline = he->is_fictitious() ? approximate_fict_edge(ctx, he) : m_bounded_approx_halfedge(he); + for(const auto& curr_pt : polyline) ctx.insert(curr_pt); } - static void approximate_inner_ccb(Context& ctx, const Vertex_const_handle& vh) { - if(vh->is_isolated()) { - approximate_vertex(ctx, vh); - return; - } - approximate_ccb(ctx, get_inner_ccb_first_halfedge(vh)); - } - - static void traverse_portal(Context& ctx, Point_geom entry_pt, const Vertex_const_handle& exit_vh) { - Point_geom exit_pt = ctx.approx(exit_vh->point()); - Approx_nt portal_x = exit_pt.x(); - // Try enforce a vertical segment. - entry_pt = Point_geom(portal_x, entry_pt.y()); - *ctx.m_out_it++ = entry_pt; - - if constexpr(is_or_derived_from_curved_surf_traits) - // Approximate the vertical segment, only meaningful when dealing with arranagements on curved surfaces. - for(Approx_nt y = entry_pt.y() - ctx.m_approx_error; y > exit_pt.y(); y -= ctx.m_approx_error) - *ctx.m_out_it++ = Point_geom(portal_x, y); - - approximate_inner_ccb(ctx, exit_vh); - // Come back here after inner ccb is approximated. - *ctx.m_out_it++ = entry_pt; - } - - static void approximate_ccb(Context& ctx, Ccb_halfedge_const_circulator start_circ) { + void approximate_ccb(Context& ctx, Ccb_halfedge_const_circulator start) const { // Try to start on a concrete halfedge. // For any unbounded face, there can't be more than 4 adjacent fictitious edges. - for(int i = 0; i < 4 && start_circ->is_fictitious(); ++i) ++start_circ; + for(int i = 0; i < 4 && start->is_fictitious(); ++i) ++start; - auto circ = start_circ; + ctx.start_ccb(); + auto circ = start; do { - circ->direction() == ARR_LEFT_TO_RIGHT ? approximate_halfedge(ctx, circ, Left_to_right_tag{}) - : approximate_halfedge(ctx, circ, Right_to_left_tag{}); + approximate_halfedge(ctx, circ); approximate_vertex(ctx, circ->target()); - } while(++circ != start_circ); + } while(++circ != start); + ctx.end_ccb(); } public: Arr_bounded_approximate_face(const Bounded_render_context& ctx) : m_ctx(ctx) - , m_bounded_approx_vertex(ctx) - , m_bounded_approx_halfedge(ctx) {} + , m_bounded_approx_halfedge(ctx) + , m_bounded_approx_vertex(ctx) {} - const Triangulated_face& operator()(const Face_const_handle& fh) const { + /*! + * \brief Approximate an arrangement face with a bunch of triangles. + * + * \param fh + * \return const Triangulated_face& + */ + const Triangle_soup& operator()(const Face_const_handle& fh) const { CGAL_precondition_msg(!fh->is_fictitious(), "Cannot approximate a fictitious face."); - auto [triangulated_face, inserted] = m_ctx.m_cache.try_emplace(fh); - if(!inserted) return triangulated_face; - if(m_ctx.is_cancelled()) return triangulated_face; + auto [iter, inserted] = m_ctx.m_cache.faces().try_emplace(fh); + Triangle_soup& ts = iter->second; + if(!inserted || m_ctx.is_cancelled()) return ts; + auto triangulator = Triangulator(m_ctx, fh); + auto ctx = Context(m_ctx, triangulator); - if(!fh->has_outer_ccb()) { - // The face is the unbounded face of bounded arrangements, we skip approximating any non degenerate features. + if(!Is_on_curved_surface && !fh->has_outer_ccb()) { + // Skip approximation of the unbounded face in planar arrangements. + // However, degenerate holes still need to be approximated. for(auto inner_ccb = fh->inner_ccbs_begin(); inner_ccb != fh->inner_ccbs_end(); ++inner_ccb) { auto circ = *inner_ccb; do { @@ -281,19 +237,23 @@ public: } while(++circ != *inner_ccb); } for(auto vh = fh->isolated_vertices_begin(); vh != fh->isolated_vertices_end(); ++vh) m_bounded_approx_vertex(vh); - return triangulated_face; + return ts; } - auto triangulator = Triangulator(m_ctx); - auto ctx = Context(m_ctx, triangulator.insert_iterator(), m_bounded_approx_vertex, m_bounded_approx_halfedge); - approximate_ccb(ctx, fh->outer_ccb()); - return triangulated_face = std::move(triangulator); + for(auto outer_ccb = fh->outer_ccbs_begin(); outer_ccb != fh->outer_ccbs_end(); ++outer_ccb) + approximate_ccb(ctx, *outer_ccb); + for(auto inner_ccb = fh->inner_ccbs_begin(); inner_ccb != fh->inner_ccbs_end(); ++inner_ccb) + approximate_ccb(ctx, *inner_ccb); + for(auto iso_vertex = fh->isolated_vertices_begin(); iso_vertex != fh->isolated_vertices_end(); ++iso_vertex) + approximate_vertex(ctx, iso_vertex); + + return ts = std::move(triangulator); } private: const Bounded_render_context& m_ctx; - const Bounded_approximate_vertex m_bounded_approx_vertex; - const Bounded_approximate_halfedge m_bounded_approx_halfedge; + Bounded_approximate_halfedge m_bounded_approx_halfedge; + Bounded_approximate_vertex m_bounded_approx_vertex; }; } // namespace draw_aos diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_halfedge.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_halfedge.h index 91a440dcadd..9654d4fdc5c 100644 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_halfedge.h +++ b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_halfedge.h @@ -1,22 +1,33 @@ +// Copyright (c) 2025 +// Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), +// INRIA Sophia-Antipolis (France), +// Max-Planck-Institute Saarbruecken (Germany), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s): Shepard Liu + #ifndef CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_HALFEDGE_H #define CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_HALFEDGE_H - #include #include -#include -#include #include -#include -#include #include #include -#include +#include +#include #include +#include #include #include -#include "CGAL/Draw_aos/Arr_projector.h" namespace CGAL { namespace draw_aos { @@ -33,196 +44,297 @@ class Arr_bounded_approximate_halfedge { using Geom_traits = typename Arrangement::Geometry_traits_2; using Halfedge_const_handle = typename Arrangement::Halfedge_const_handle; + using Gt_point = typename Geom_traits::Point_2; using Approx_traits = Arr_approximate_traits; using Approx_nt = typename Approx_traits::Approx_nt; using Approx_point = typename Approx_traits::Approx_point; - using Point_geom = typename Approx_traits::Point_geom; - using Polyline_geom = typename Approx_traits::Polyline_geom; + using Point = typename Approx_traits::Point; + using Polyline = typename Approx_traits::Polyline; using Approx_kernel = typename Approx_traits::Approx_kernel; using Approx_line_2 = typename Approx_kernel::Line_2; - using X_monotone_curve_2 = typename Geom_traits::X_monotone_curve_2; - using Point_2 = typename Geom_traits::Point_2; - using Bounded_render_context = Arr_bounded_render_context; using Boundary_lines = std::array; + constexpr static bool Has_approximate_xcv_with_bounds = + has_approximate_xcv_with_bounds_v; + private: struct Context : public Bounded_render_context { - Context(const Bounded_render_context& ctx, - const X_monotone_curve_2& curve, - Polyline_geom& polyline, - const Boundary_lines& boundary_lines) + Context(const Bounded_render_context& ctx, const X_monotone_curve_2& curve, Polyline& polyline) : Bounded_render_context(ctx) - , m_curve(curve) - , m_boundary_lines(boundary_lines) - , m_proj(ctx.m_traits) - , m_base_out_it(std::back_inserter(polyline)) - , m_out_it(boost::make_function_output_iterator(std::function([this](Point_geom pt) { - if(pt.x() < this->xmin()) { - // We need the last point if not yet x-inbound. - m_last_pt = pt; - return; - } else if(pt.x() > this->xmax()) - return; - - *m_base_out_it++ = pt; - m_last_pt = pt; - }))) {} - - private: - std::back_insert_iterator m_base_out_it; + , m_polyline(polyline) + , m_curve(curve) {} + // Prevent accidental copying. + Context(const Context&) = delete; + Context& operator=(const Context&) = delete; public: - Arr_projector m_proj; + /*! + * \brief Insert a point to the polyline if it is within the x-range of the curve + * \note Will be replaced after AosApproximateUnboundedTraits_2 is fully available. + * \param pt + */ + void insert(Point pt) { + if(pt.x() < this->xmin()) { + // We need the last point if not yet x-inbound. + m_last_pt = pt; + return; + } else if(pt.x() > this->xmax()) + return; + + m_polyline.push_back(pt); + m_last_pt = pt; + } + + const std::optional& last_pt() const { return m_last_pt; } + + private: + std::optional m_last_pt; + + public: + Polyline& m_polyline; const X_monotone_curve_2& m_curve; - const Boundary_lines& m_boundary_lines; - std::optional m_last_pt; - boost::function_output_iterator> m_out_it; }; - static Point_geom trace_boundary_inter(const Context& ctx, Point_geom pt, Side_of_boundary side) { - Point_geom inter = std::get( - *CGAL::intersection(Approx_line_2(*ctx.m_last_pt, pt), ctx.m_boundary_lines[static_cast(side)])); - // Prevent floating point errors. + /*! + * \brief Computes the intersection point between the given boundary side and the line segment from last_pt to pt. + */ + Point boundary_intersection(const Context& ctx, Point pt, Boundary_side side) const { + std::optional x, y; + const Approx_line_2* line; switch(side) { - case Side_of_boundary::Left: - return Point_geom(ctx.xmin(), inter.y()); - case Side_of_boundary::Right: - return Point_geom(ctx.xmax(), inter.y()); - case Side_of_boundary::Top: - return Point_geom(inter.x(), ctx.ymax()); - case Side_of_boundary::Bottom: - return Point_geom(inter.x(), ctx.ymin()); + case Boundary_side::Left: + x = ctx.xmin(); + line = &m_left; + break; + case Boundary_side::Right: + x = ctx.xmax(); + line = &m_right; + break; + case Boundary_side::Top: + y = ctx.ymax(); + line = &m_top; + break; + case Boundary_side::Bottom: + y = ctx.ymin(); + line = &m_bottom; + break; default: CGAL_assertion(false && "Unexpected side of boundary."); - return Point_geom(); + } + Point inter = std::get(*CGAL::intersection(Approx_line_2(*ctx.last_pt(), pt), *line)); + if(x.has_value()) return Point(*x, inter.y()); + return Point(inter.x(), *y); + } + + /*! + * \brief Trace approximated curve point in ltr ordering, adding boundary intersections if necessary. + * + * \note This method will eventually be replaced by AosApproximateUnboundedTraits_2. + */ + void trace_add(Context& ctx, Point pt) const { + if(!ctx.last_pt().has_value()) { + ctx.insert(pt); + return; + } + if(ctx.last_pt()->x() < ctx.xmin() && pt.x() >= ctx.xmin()) + ctx.insert(boundary_intersection(ctx, pt, Boundary_side::Left)); + if(ctx.last_pt()->y() < ctx.ymin()) { + if(pt.y() > ctx.ymin()) ctx.insert(boundary_intersection(ctx, pt, Boundary_side::Bottom)); + if(pt.y() > ctx.ymax()) ctx.insert(boundary_intersection(ctx, pt, Boundary_side::Top)); + } else if(ctx.last_pt()->y() > ctx.ymax()) { + if(pt.y() < ctx.ymax()) ctx.insert(boundary_intersection(ctx, pt, Boundary_side::Top)); + if(pt.y() < ctx.ymin()) ctx.insert(boundary_intersection(ctx, pt, Boundary_side::Bottom)); + } else { + if(pt.y() < ctx.ymin()) + ctx.insert(boundary_intersection(ctx, pt, Boundary_side::Bottom)); + else if(pt.y() > ctx.ymax()) + ctx.insert(boundary_intersection(ctx, pt, Boundary_side::Top)); + } + if(ctx.last_pt()->x() <= ctx.xmax() && pt.x() > ctx.xmax()) + ctx.insert(boundary_intersection(ctx, pt, Boundary_side::Right)); + ctx.insert(pt); + } + + /*! + * \brief Check if the point is within the x-range of the curve. + */ + static bool is_in_x_range(const Context& ctx, const Gt_point& pt) { + const Geom_traits& traits = ctx.m_traits; + const X_monotone_curve_2& curve = ctx.m_curve; + + if constexpr(has_is_in_x_range_v) return curve.is_in_x_range(pt); + if constexpr(!has_parameter_space_in_x_2::value) { + const auto& min_pt = traits.construct_min_vertex_2_object()(curve); + const auto& max_pt = traits.construct_max_vertex_2_object()(curve); + return traits.compare_x_2_object()(pt, min_pt) != CGAL::SMALLER && + traits.compare_x_2_object()(pt, max_pt) != CGAL::LARGER; + } + + Comparison_result left_cmp; + if(auto left_loc = traits.parameter_space_in_x_2_object()(curve, ARR_MIN_END); left_loc == ARR_INTERIOR) + left_cmp = traits.compare_x_2_object()(pt, traits.construct_min_vertex_2_object()(curve)); + else if(left_loc == ARR_LEFT_BOUNDARY) + left_cmp = CGAL::LARGER; + else + left_cmp = traits.compare_x_on_boundary_2_object()(pt, curve, ARR_MIN_END); + if(left_cmp == CGAL::SMALLER) return false; + if(left_cmp == CGAL::EQUAL) return true; + + Comparison_result right_cmp; + if(auto right_loc = traits.parameter_space_in_x_2_object()(curve, ARR_MAX_END); right_loc == ARR_INTERIOR) + right_cmp = traits.compare_x_2_object()(pt, traits.construct_max_vertex_2_object()(curve)); + else if(right_loc == ARR_RIGHT_BOUNDARY) + right_cmp = CGAL::SMALLER; + else + right_cmp = traits.compare_x_on_boundary_2_object()(pt, curve, ARR_MAX_END); + return right_cmp != CGAL::LARGER; + } + + /*! + * \brief transform approximated curve points(ltr ordering) in place based on the halfedge, giving correct + * ordering, continuity, etc. + */ + static void transform_polyline(Context& ctx, Polyline& polyline, const Halfedge_const_handle& he) { + transform_polyline_impl(ctx, polyline, he); + } + + // For planar arrangements, we only need to reverse the polyline if the halfedge is rtl. + template , int> = 0> + static void transform_polyline_impl(Context&, Polyline& polyline, const Halfedge_const_handle& he) { + if(he->direction() == CGAL::ARR_LEFT_TO_RIGHT) return; + std::reverse(polyline.begin(), polyline.end()); + } + + template , int> = 0> + static void transform_polyline_impl(Context& ctx, Polyline& polyline, const Halfedge_const_handle& he) { + using Direction_3 = typename Geom_traits::Direction_3; + using Vector_3 = typename Geom_traits::Vector_3; + + if(polyline.size() < 2) return; + const X_monotone_curve_2& curve = he->curve(); + const auto& traits = ctx.m_traits; + if(curve.is_vertical()) { + Direction_3 normal_dir = curve.is_directed_right() ? curve.normal() : -curve.normal(); + Direction_3 azimuth_dir(CGAL::cross_product(Vector_3(0, 0, 1), normal_dir.vector())); + Approx_nt azimuth = ctx.to_uv(traits.approximate_2_object()(traits.construct_point_2_object()(azimuth_dir))).x(); + if(azimuth == 0 && he->direction() == ARR_LEFT_TO_RIGHT) azimuth = 2 * CGAL_PI; + std::transform(polyline.begin(), polyline.end(), polyline.begin(), + [azimuth](Point pt) { return Point(azimuth, pt.y()); }); + } else if(polyline.back().x() == 0) { + // For strictly x-monotone arcs whose target point sits on the boundary, the x should be set to 2 * CGAL_PI + polyline.back() = Point(2 * CGAL_PI, polyline.back().y()); + } + if(he->direction() == CGAL::ARR_LEFT_TO_RIGHT) return; + std::reverse(polyline.begin(), polyline.end()); + } + + void approximate_curve(Context& ctx) const { approximate_curve_impl(ctx); } + + // If Approximate_2 supports curve approximation with bounding box + template , int> = 0> + void approximate_curve_impl(Context& ctx) const { + const Geom_traits& traits = ctx.m_traits; + const X_monotone_curve_2& curve = ctx.m_curve; + Polyline& polyline = ctx.m_polyline; + auto compare_y_at_x_2 = traits.compare_y_at_x_2_object(); + + if(is_in_x_range(ctx, m_top_left)) { + if(compare_y_at_x_2(m_top_left, curve) == CGAL::SMALLER) { + polyline.insert(polyline.end(), {Approx_traits::Null_point, Point(ctx.xmin(), ctx.ymax())}); + } else if(compare_y_at_x_2(m_bottom_left, curve) == CGAL::LARGER) { + polyline.insert(polyline.end(), {Approx_traits::Null_point, Point(ctx.xmin(), ctx.ymin())}); + } + } + traits.approximate_2_object()(curve, ctx.m_approx_error, + boost::make_function_output_iterator([&ctx, this](Approx_point approx_pt) { + ctx.m_polyline.push_back(snap_to_boundary(ctx, ctx.to_uv(approx_pt))); + }), + ctx.bbox(), true); + if(is_in_x_range(ctx, m_top_right)) { + if(compare_y_at_x_2(m_top_right, curve) == CGAL::SMALLER) { + polyline.insert(polyline.end(), {Point(ctx.xmax(), ctx.ymax()), Approx_traits::Null_point}); + } else if(compare_y_at_x_2(m_bottom_right, curve) == CGAL::LARGER) { + polyline.insert(polyline.end(), {Point(ctx.xmax(), ctx.ymin()), Approx_traits::Null_point}); + } } } - static void update(Context& ctx, Point_geom pt) { - if(!ctx.m_last_pt.has_value()) { - *ctx.m_out_it++ = pt; - return; - } + // If Approximate_2 does not support curve approximation with bounding box + template , int> = 0> + void approximate_curve_impl(Context& ctx) const { + m_ctx.m_traits.approximate_2_object()( + ctx.m_curve, ctx.m_approx_error, + boost::make_function_output_iterator([&ctx, this](Approx_point pt) { trace_add(ctx, ctx.to_uv(pt)); }), true); + } - if(ctx.m_last_pt->x() < ctx.xmin() && pt.x() >= ctx.xmin()) { - *ctx.m_out_it++ = trace_boundary_inter(ctx, pt, Side_of_boundary::Left); - } - - if(ctx.m_last_pt->y() < ctx.ymin()) { - if(pt.y() > ctx.ymin()) { - *ctx.m_out_it++ = trace_boundary_inter(ctx, pt, Side_of_boundary::Bottom); - } - if(pt.y() > ctx.ymax()) { - *ctx.m_out_it++ = trace_boundary_inter(ctx, pt, Side_of_boundary::Top); - } - } else if(ctx.m_last_pt->y() > ctx.ymax()) { - if(pt.y() < ctx.ymax()) { - *ctx.m_out_it++ = trace_boundary_inter(ctx, pt, Side_of_boundary::Top); - } - if(pt.y() < ctx.ymin()) { - *ctx.m_out_it++ = trace_boundary_inter(ctx, pt, Side_of_boundary::Bottom); - } - } else { - if(pt.y() < ctx.ymin()) { - *ctx.m_out_it++ = trace_boundary_inter(ctx, pt, Side_of_boundary::Bottom); - } else if(pt.y() > ctx.ymax()) { - *ctx.m_out_it++ = trace_boundary_inter(ctx, pt, Side_of_boundary::Top); - } - } - - if(ctx.m_last_pt->x() <= ctx.xmax() && pt.x() > ctx.xmax()) { - *ctx.m_out_it++ = trace_boundary_inter(ctx, pt, Side_of_boundary::Right); - } - - *ctx.m_out_it++ = pt; + /*! + * \brief Adjusts a point by snapping it to the nearest boundary to reduce floating-point error. + * + * \return The adjusted (snapped) point if it lies within snapping tolerance, or the original point otherwise. + */ + Point snap_to_boundary(const Context& ctx, Point pt) const { + Approx_nt x = pt.x(), y = pt.y(); + if(std::abs(x - ctx.xmin()) < m_ep_left) + x = ctx.xmin(); + else if(std::abs(x - ctx.xmax()) < m_ep_right) + x = ctx.xmax(); + if(std::abs(y - ctx.ymin()) < m_ep_bottom) + y = ctx.ymin(); + else if(std::abs(y - ctx.ymax()) < m_ep_top) + y = ctx.ymax(); + return Point(x, y); } public: Arr_bounded_approximate_halfedge(const Bounded_render_context& ctx) : m_ctx(ctx) - , m_boundary_lines({ - Approx_line_2(Point_geom(ctx.xmin(), ctx.ymax()), Point_geom(ctx.xmax(), ctx.ymax())), // Top = 0 - Approx_line_2(Point_geom(ctx.xmin(), ctx.ymin()), Point_geom(ctx.xmin(), ctx.ymax())), // Left = 1 - Approx_line_2(Point_geom(ctx.xmin(), ctx.ymin()), Point_geom(ctx.xmax(), ctx.ymin())), // Bottom = 2 - Approx_line_2(Point_geom(ctx.xmax(), ctx.ymin()), Point_geom(ctx.xmax(), ctx.ymax())), // Right = 3 - }) {} + , m_top(ctx.top_left(), ctx.top_right()) + , m_bottom(ctx.bottom_left(), ctx.bottom_right()) + , m_left(ctx.bottom_left(), ctx.top_left()) + , m_right(ctx.bottom_right(), ctx.top_right()) { + Construct_gt_point_2 ctr_p; + m_top_left = ctr_p(ctx.to_cartesian(ctx.top_left())); + m_top_right = ctr_p(ctx.to_cartesian(ctx.top_right())); + m_bottom_left = ctr_p(ctx.to_cartesian(ctx.bottom_left())); + m_bottom_right = ctr_p(ctx.to_cartesian(ctx.bottom_right())); + Approx_nt ep_base = std::numeric_limits::epsilon(); + m_ep_left = std::max(std::abs(ep_base * ctx.xmin()), ep_base); + m_ep_right = std::max(std::abs(ep_base * ctx.xmax()), ep_base); + m_ep_bottom = std::max(std::abs(ep_base * ctx.ymin()), ep_base); + m_ep_top = std::max(std::abs(ep_base * ctx.ymax()), ep_base); + } - const Polyline_geom& operator()(const Halfedge_const_handle& he) const { + const Polyline& operator()(const Halfedge_const_handle& he) const { CGAL_assertion(!he->is_fictitious()); - auto [polyline, inserted] = m_ctx.m_cache.try_emplace(he); + auto& cache = m_ctx.m_cache.halfedges(); + auto [iter, inserted] = cache.try_emplace(he, Polyline()); + Polyline& polyline = iter->second; if(!inserted) return polyline; if(m_ctx.is_cancelled()) return polyline; const X_monotone_curve_2& curve = he->curve(); - Context ctx(m_ctx, curve, polyline, m_boundary_lines); - m_ctx.m_traits.approximate_2_object()( - curve, m_ctx.m_approx_error, - boost::make_function_output_iterator([&ctx](Approx_point pt) { update(ctx, ctx.m_proj.project(pt)); }), - true); // ltr ordering + Context ctx(m_ctx, curve, polyline); + approximate_curve(ctx); + Polyline poly_copy(polyline); + transform_polyline(ctx, polyline, he); // also approximate the twin halfedge - auto [twin_poly, twin_inserted] = m_ctx.m_cache.try_emplace(he->twin()); - twin_poly = polyline; - if(twin_inserted) adapt_polyline(ctx, twin_poly, he->twin()); - - adapt_polyline(ctx, polyline, he); - return polyline; - } - - /** - * @brief Functor to adapt approximated curve points based on the actual halfedge, giving correct ordering, - * continuity, etc. - */ - static void adapt_polyline(Context& ctx, Polyline_geom& polyline, const Halfedge_const_handle& he) { - adapt_polyline_impl(ctx, polyline, he); - } - - template , int> = 0> - static void adapt_polyline_impl(Context& ctx, Polyline_geom& polyline, const Halfedge_const_handle& he) { - if(he->direction() == CGAL::ARR_LEFT_TO_RIGHT) return; - std::reverse(polyline.begin(), polyline.end()); - } - - template , int> = 0> - static void adapt_polyline_impl(Context& ctx, Polyline_geom& polyline, const Halfedge_const_handle& he) { - if(polyline.size() < 2) return; - - auto azimuth_of_vertical_curve = [&ctx](const X_monotone_curve_2& curve) -> Approx_nt { - using Point_2 = typename Geom_traits::Point_2; - using Direction_3 = typename Geom_traits::Direction_3; - - const auto& traits = ctx.m_traits; - const Direction_3& normal_dir = curve.normal(); - if(normal_dir.dx() == 0 && normal_dir.dy() == 1) return 0; // overlaps with the identification curve - Point_2 normal_pt(normal_dir, Point_2::NO_BOUNDARY_LOC); - Approx_point normal_approx_pt = traits.approximate_2_object()(normal_pt); - Point_geom normal_proj_pt = ctx.m_proj.project(normal_approx_pt); - return std::fmod(normal_proj_pt.x() + 0.5 * CGAL_PI, CGAL_PI * 2.0); - }; - - const X_monotone_curve_2& curve = he->curve(); - - if(curve.is_vertical()) { - Approx_nt azimuth = azimuth_of_vertical_curve(curve); - if(azimuth == 0 && he->direction() == ARR_LEFT_TO_RIGHT) azimuth = 2 * CGAL_PI; - std::transform(polyline.begin(), polyline.end(), polyline.begin(), - [azimuth](Point_geom pt) { return Point_geom(azimuth, pt.y()); }); - } else if(polyline.back().x() == 0) { - // For strictly x-monotone arcs, if the target point sits on the boundary, the x should be set to 2 * CGAL_PI - polyline.back() = Point_geom(2 * CGAL_PI, polyline.back().y()); - } - - if(he->direction() == ARR_RIGHT_TO_LEFT) std::reverse(polyline.begin(), polyline.end()); + auto [twin_iter, twin_inserted] = cache.try_emplace(he->twin(), std::move(poly_copy)); + if(twin_inserted) transform_polyline(ctx, twin_iter->second, he->twin()); + // The previous iterator might have been invalidated by the second try_emplace call, so we do an extra lookup. + return cache.at(he); } private: const Bounded_render_context& m_ctx; - const Boundary_lines m_boundary_lines; + Approx_line_2 m_left, m_right, m_top, m_bottom; + Gt_point m_top_left, m_top_right, m_bottom_left, m_bottom_right; + Approx_nt m_ep_left, m_ep_right, m_ep_bottom, m_ep_top; }; } // namespace draw_aos diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_vertex.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_vertex.h index 709fd1f747d..80729fa18e1 100644 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_vertex.h +++ b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_vertex.h @@ -1,10 +1,23 @@ +// Copyright (c) 2025 +// Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), +// INRIA Sophia-Antipolis (France), +// Max-Planck-Institute Saarbruecken (Germany), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s): Shepard Liu #ifndef CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_VERTEX_H #define CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_VERTEX_H #include #include -#include "CGAL/Draw_aos/Arr_projector.h" namespace CGAL { namespace draw_aos { @@ -15,7 +28,7 @@ class Arr_bounded_approximate_vertex using Geom_traits = typename Arrangement::Geometry_traits_2; using Point_2 = typename Geom_traits::Point_2; using Vertex_const_handle = typename Arrangement::Vertex_const_handle; - using Point_geom = typename Arr_approximate_traits::Point_geom; + using Point_geom = typename Arr_approximate_traits::Point; using Bounded_render_context = Arr_bounded_render_context; public: @@ -32,9 +45,10 @@ public: * @return const Point_geom& */ const Point_geom& operator()(const Vertex_const_handle& vh) const { - auto [point, inserted] = m_ctx.m_cache.try_emplace(vh); + auto [iter, inserted] = m_ctx.m_cache.vertices().try_emplace(vh); + Point_geom& point = iter->second; if(!inserted) return point; - return point = Arr_projector(m_ctx.m_traits).project(m_ctx.m_traits.approximate_2_object()(vh->point())); + return point = m_ctx.to_uv(m_ctx.m_traits.approximate_2_object()(vh->point())); } private: diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_face_triangulator.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_face_triangulator.h index a8f88ae6f3e..4384190df86 100644 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_face_triangulator.h +++ b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_face_triangulator.h @@ -1,10 +1,23 @@ +// Copyright (c) 2025 +// Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), +// INRIA Sophia-Antipolis (France), +// Max-Planck-Institute Saarbruecken (Germany), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s): Shepard Liu + #ifndef CGAL_DRAW_AOS_ARR_FACE_TRIANGULATOR_H #define CGAL_DRAW_AOS_ARR_FACE_TRIANGULATOR_H #include -#include -#include -#include +#include #include #include #include @@ -22,7 +35,6 @@ #include #include #include -#include #if defined(CGAL_DRAW_AOS_DEBUG) && defined(CGAL_DRAW_AOS_TRIANGULATOR_DEBUG_FILE_DIR) #include @@ -40,186 +52,197 @@ namespace draw_aos { /** * @brief Triangulator for a face of an arrangement within a bounding box. - * - * The original topology of a face is preserved, but the geometry will be bounded by the bbox. */ template class Arr_bounded_face_triangulator { using Geom_traits = typename Arrangement::Geometry_traits_2; + constexpr static bool Is_on_curved_surface = is_or_derived_from_curved_surf_traits_v; + using Approx_traits = Arr_approximate_traits; - using Point_geom = typename Approx_traits::Point_geom; - using Triangle = typename Approx_traits::Triangle; - using Triangulated_face = typename Approx_traits::Triangulated_face; + using Point = typename Approx_traits::Point; + using Approx_point = typename Approx_traits::Approx_point; + using Approx_kernel = typename Approx_traits::Approx_kernel; + using Triangle_soup = typename Approx_traits::Triangle_soup; + using Triangle = typename Triangle_soup::Triangle; + using Face_const_handle = typename Arrangement::Face_const_handle; #if defined(CGAL_DRAW_AOS_DEBUG) template friend void debug_print(const Arr_bounded_face_triangulator& triangulator); #endif - struct Point_index + enum Point_type { Vertex_only, Constraint_only, Vertex_and_constraint }; + + /*! + * \brief A index wrapper defaulted to invalid. + */ + class Index { - constexpr static std::size_t Invalid_index = std::numeric_limits::max(); - std::size_t index{Invalid_index}; - Point_index() = default; - Point_index(std::size_t idx) - : index(idx) {} - bool is_valid() const { return index != Invalid_index; } - operator std::size_t() const { return index; } + public: + Index() = default; + Index(int idx) + : m_index(idx) {} + + bool is_valid() const { return m_index != Invalid_index; } + operator int() const { return m_index; } + + private: + constexpr static int Invalid_index = -1; + int m_index{Invalid_index}; }; + using Epick = Exact_predicates_inexact_constructions_kernel; - using Vb = Triangulation_vertex_base_with_info_2; + using Vb = Triangulation_vertex_base_with_info_2; using Fb = Constrained_triangulation_face_base_2; using Tds = Triangulation_data_structure_2; // For planar arrangements, Constrained_triangulation_2 is enough. - using Ct = std::conditional_t, + using Ct = std::conditional_t, Constrained_triangulation_2>; - using KPoint = Epick::Point_2; - using KPoint_with_info = std::pair; + using KPoint = Epick::Point_2; + using KPoint_with_index = std::pair; using Bounded_render_context = Arr_bounded_render_context; public: - using Insert_iterator = boost::function_output_iterator>; + using value_type = Point; private: - static KPoint transform_point(Point_geom pt) { return KPoint(pt.x(), pt.y()); } + static KPoint to_kpoint(Point pt) { return KPoint(pt.x(), pt.y()); } + + /*! + * \brief Offset a point on a specific boundary outward by a given offset. + * + * \pre side != Boundary_side::None + */ + static Point offset_boundary_point(Point pt, Boundary_side side, double offset) { + CGAL_precondition(side != Boundary_side::None); - static Point_geom offset_boundary_point(Point_geom pt, Side_of_boundary side, double offset) { - CGAL_precondition(side != Side_of_boundary::None); switch(side) { - case Side_of_boundary::Left: - return Point_geom(pt.x() - offset, pt.y()); - case Side_of_boundary::Right: - return Point_geom(pt.x() + offset, pt.y()); - case Side_of_boundary::Top: - return Point_geom(pt.x(), pt.y() + offset); - case Side_of_boundary::Bottom: - return Point_geom(pt.x(), pt.y() - offset); + case Boundary_side::Left: + return Point(pt.x() - offset, pt.y()); + case Boundary_side::Right: + return Point(pt.x() + offset, pt.y()); + case Boundary_side::Top: + return Point(pt.x(), pt.y() + offset); + case Boundary_side::Bottom: + return Point(pt.x(), pt.y() - offset); default: return pt; // Should not reach here } } - /** - * Sample the interior of a face for geometry traits on curved surfaces + /*! + * \brief Find the shared boundary side of two points, or None if they are not on the same boundary. */ - void sample_interior() { sample_interior_impl(); } - - template , int> = 0> - void sample_interior_impl() {} - - template , int> = 0> - void sample_interior_impl() { - const double scale_factor = 20.0; - double cell_size = m_ctx.m_approx_error * scale_factor; - std::size_t grid_x_min = m_face_bbox.xmin() / cell_size; - std::size_t grid_x_max = std::ceil(m_face_bbox.xmax() / cell_size); - std::size_t grid_y_min = m_face_bbox.ymin() / cell_size; - std::size_t grid_y_max = std::ceil(m_face_bbox.ymax() / cell_size); - for(std::size_t x = grid_x_min + 1; x < grid_x_max; ++x) - for(std::size_t y = grid_y_min + 1; y < grid_y_max; ++y) - m_points.emplace_back(Point_geom(x * cell_size, y * cell_size)); - // sample the bbox boundary - for(double x = m_face_bbox.xmin(); x <= m_face_bbox.xmax(); x += cell_size) { - m_points.emplace_back(Point_geom(x, m_face_bbox.ymin())); - m_points.emplace_back(Point_geom(x, m_face_bbox.ymax())); - } - for(double y = m_face_bbox.ymin(); y <= m_face_bbox.ymax(); y += cell_size) { - m_points.emplace_back(Point_geom(m_face_bbox.xmin(), y)); - m_points.emplace_back(Point_geom(m_face_bbox.xmax(), y)); - } + Boundary_side shared_boundary(const Point& pt1, const Point& pt2) const { + if(m_ctx.is_on_left(pt1) && m_ctx.is_on_left(pt2)) return Boundary_side::Left; + if(m_ctx.is_on_right(pt1) && m_ctx.is_on_right(pt2)) return Boundary_side::Right; + if(m_ctx.is_on_bottom(pt1) && m_ctx.is_on_bottom(pt2)) return Boundary_side::Bottom; + if(m_ctx.is_on_top(pt1) && m_ctx.is_on_top(pt2)) return Boundary_side::Top; + return Boundary_side::None; } - void insert_ccb() { - auto begin = m_points.begin(); - auto end = m_points.end(); - auto helpers_iter = m_helper_indices.begin(); - auto helpers_end = m_helper_indices.end(); - // A two pointers algorithm implemented with boost filters. - auto point_filter = [&helpers_iter, helpers_end](std::size_t idx) { - if(helpers_iter == helpers_end) { - return true; - } - if(idx == *helpers_iter) { - ++helpers_iter; - return false; - } - return true; - }; - auto index_to_point_with_info = [this](std::size_t idx) -> KPoint_with_info { - return std::make_pair(transform_point(m_points[idx]), idx); - }; - auto indexes_begin = boost::make_counting_iterator(0); - auto indexes_end = boost::make_counting_iterator(m_points.size()); - auto filtered_begin = boost::make_filter_iterator(point_filter, indexes_begin, indexes_end); - auto filtered_end = boost::make_filter_iterator(point_filter, indexes_end, indexes_end); - auto transformed_begin = boost::make_transform_iterator(filtered_begin, index_to_point_with_info); - auto transformed_end = boost::make_transform_iterator(filtered_end, index_to_point_with_info); - // Constrained_triangulation_2 and Constrained_Delaunay_triangulation_2 have slightly different interfaces. - if constexpr(is_or_derived_from_curved_surf_traits) - m_ct.insert(transformed_begin, transformed_end); - else - m_ct.template insert_with_info(transformed_begin, transformed_end); - } - - Side_of_boundary shared_boundary(const Point_geom& pt1, const Point_geom& pt2) const { - if(pt1.x() == m_ctx.xmin() && pt2.x() == m_ctx.xmin() && m_ctx.contains_y(pt1.y()) && m_ctx.contains_y(pt2.y())) - return Side_of_boundary::Left; - if(pt1.x() == m_ctx.xmax() && pt2.x() == m_ctx.xmax() && m_ctx.contains_y(pt1.y()) && m_ctx.contains_y(pt2.y())) - return Side_of_boundary::Right; - if(pt1.y() == m_ctx.ymin() && pt2.y() == m_ctx.ymin() && m_ctx.contains_x(pt1.x()) && m_ctx.contains_x(pt2.x())) - return Side_of_boundary::Bottom; - if(pt1.y() == m_ctx.ymax() && pt2.y() == m_ctx.ymax() && m_ctx.contains_x(pt1.x()) && m_ctx.contains_x(pt2.x())) - return Side_of_boundary::Top; - return Side_of_boundary::None; - } - - void add_boundary_helper_point(Point_geom from, Point_geom to) { + /*! + * \brief Add a helper point on the shared boundary of two points if they are on the same boundary side. + * + * When triangulating a arrangement face within a bounding box, curves outside the bounding box are projected on the + * four sides of the bbox. Topological errors could be introduced if several polylines are lying on the same side. + * Thus we add the midpoint in between the two points on boundary and move it outward with an increasing offset. + */ + void add_boundary_helper_point(Point from, Point to) { + // Arrangements on curved surfaces currently draws the entire parameter space, so there's no need to add + // helper points. + if constexpr(Is_on_curved_surface) return; if(from == to) return; auto shared_side = shared_boundary(from, to); - if(shared_side == Side_of_boundary::None) return; - m_helper_indices.push_back(m_points.size()); - m_points.emplace_back( - offset_boundary_point(Point_geom((from.x() + to.x()) / 2, (from.y() + to.y()) / 2), shared_side, m_offset)); - m_offset += 0.5; + if(shared_side == Boundary_side::None) return; + Point mid = CGAL::midpoint(from, to); + m_points.push_back(offset_boundary_point(mid, shared_side, m_offset += 0.1)); + m_point_types.push_back(Constraint_only); + } + + void insert_all_vertices() { + auto vertex_filter = [this](int idx) { return m_point_types[idx] != Constraint_only; }; + auto index_to_point_with_info = [this](int idx) -> KPoint_with_index { + return std::make_pair(to_kpoint(m_points[idx]), idx); + }; + auto indexes_begin = boost::make_counting_iterator(0); + auto indexes_end = boost::make_counting_iterator(m_points.size()); + auto filtered_begin = boost::make_filter_iterator(vertex_filter, indexes_begin, indexes_end); + auto filtered_end = boost::make_filter_iterator(vertex_filter, indexes_end, indexes_end); + auto transformed_begin = boost::make_transform_iterator(filtered_begin, index_to_point_with_info); + auto transformed_end = boost::make_transform_iterator(filtered_end, index_to_point_with_info); + + // Constrained_triangulation_2 and Constrained_Delaunay_triangulation_2 have slightly different interfaces. + if constexpr(Is_on_curved_surface) + m_ct.insert(transformed_begin, transformed_end); + else + m_ct.template insert_with_info(transformed_begin, transformed_end); + } + + void insert_all_constraints() { + auto constraint_filter = [this](int idx) { return m_point_types[idx] != Vertex_only; }; + auto index_to_point = [this](int idx) -> KPoint { return to_kpoint(m_points[idx]); }; + for(auto [start_idx, end_idx] : m_cst_ranges) { + auto indexes_begin = boost::make_counting_iterator(start_idx); + auto indexes_end = boost::make_counting_iterator(end_idx); + auto filtered_begin = boost::make_filter_iterator(constraint_filter, indexes_begin, indexes_end); + auto filtered_end = boost::make_filter_iterator(constraint_filter, indexes_end, indexes_end); + auto transformed_begin = boost::make_transform_iterator(filtered_begin, index_to_point); + auto transformed_end = boost::make_transform_iterator(filtered_end, index_to_point); + m_ct.insert_constraint(transformed_begin, transformed_end, true); + } } public: - Arr_bounded_face_triangulator(const Bounded_render_context& ctx) - : m_ctx(ctx) {} + Arr_bounded_face_triangulator(const Bounded_render_context& ctx, Face_const_handle fh) + : m_ctx(ctx) + , m_fh(fh) {} - Insert_iterator insert_iterator() { - return boost::make_function_output_iterator(std::function([this](Point_geom pt) { - CGAL_assertion_msg(m_ctx.contains(pt), "Point should be within the bounding box."); + void push_back(Point pt) { + CGAL_assertion_msg(m_curr_cst_begin.has_value(), "Call start_constraint() before push_back()."); - if(m_points.size() != 0) add_boundary_helper_point(m_points.back(), pt); - m_points.emplace_back(pt); - this->m_face_bbox += pt.bbox(); - })); + if(m_points.size() - *m_curr_cst_begin >= 2) add_boundary_helper_point(m_points.back(), pt); + m_points.push_back(pt); + m_point_types.push_back(Vertex_and_constraint); } - /** - * @brief Converts the triangulator to a triangulated face, moving internal data to the result. - * - * @return Triangulated_face - */ - operator Triangulated_face() && { - if(m_points.empty()) { - return Triangulated_face(); - } + void start_constraint() { m_curr_cst_begin = m_points.size(); } - m_constraint_end = m_points.size(); - add_boundary_helper_point(m_points.back(), m_points.front()); - sample_interior(); - insert_ccb(); - // insert_constraint() should be called after insert_with_info(), or info will not be set correctly. - m_ct.insert_constraint(boost::make_transform_iterator(m_points.begin(), transform_point), - boost::make_transform_iterator(m_points.begin() + m_constraint_end, transform_point), true); - if(m_ct.number_of_faces() == 0) { - return Triangulated_face(); + void end_constraint() { + CGAL_assertion_msg(m_curr_cst_begin.has_value(), "Call start_constraint() before end_constraint()."); + + int cst_begin = *m_curr_cst_begin; + m_curr_cst_begin.reset(); + if(m_points.size() - cst_begin <= 2) { + m_points.erase(m_points.begin() + cst_begin, m_points.end()); + m_point_types.erase(m_point_types.begin() + cst_begin, m_point_types.end()); + return; } + add_boundary_helper_point(m_points.back(), m_points[cst_begin]); + m_cst_ranges.emplace_back(cst_begin, m_points.size()); + } + + /*! + * \brief Converts the triangulator to a triangulated face, moving internal data to the result. + * + * \return Triangulated_face + */ + operator Triangle_soup() && { + CGAL_assertion_msg(!m_curr_cst_begin.has_value(), "Call end_constraint() before conversion"); + + if(m_points.empty()) return Triangle_soup(); + if constexpr(Is_on_curved_surface) { + if(auto it = m_ctx.m_face_points.find(m_fh); it != m_ctx.m_face_points.end()) { + m_points.insert(m_points.end(), it->second.begin(), it->second.end()); + m_point_types.insert(m_point_types.end(), it->second.size(), Vertex_only); + } + } + insert_all_vertices(); + insert_all_constraints(); + if(m_ct.number_of_faces() == 0) return Triangle_soup(); #if defined(CGAL_DRAW_AOS_DEBUG) debug_print(*this); @@ -230,29 +253,29 @@ public: boost::associative_property_map in_domain(in_domain_map); CGAL::mark_domain_in_triangulation(m_ct, in_domain); // Collect triangles within the constrained domain. - Triangulated_face tf; - tf.triangles.reserve(m_ct.number_of_faces()); + Triangle_soup ts; + ts.triangles.reserve(m_ct.number_of_faces()); for(auto fit = m_ct.finite_faces_begin(); fit != m_ct.finite_faces_end(); ++fit) { - if(!get(in_domain, fit)) continue; - Point_index v1 = fit->vertex(0)->info(); - Point_index v2 = fit->vertex(1)->info(); - Point_index v3 = fit->vertex(2)->info(); + Index v1 = fit->vertex(0)->info(); + Index v2 = fit->vertex(1)->info(); + Index v3 = fit->vertex(2)->info(); if(!v1.is_valid() || !v2.is_valid() || !v3.is_valid()) continue; - Triangle tri{v1, v2, v3}; - tf.triangles.emplace_back(tri); + if(!get(in_domain, fit)) continue; + ts.triangles.push_back(Triangle{v1, v2, v3}); } - tf.points = std::move(m_points); - return tf; + ts.points = std::move(m_points); + return ts; } private: const Bounded_render_context& m_ctx; + Face_const_handle m_fh; Ct m_ct; - std::vector m_points; - Bbox_2 m_face_bbox; // The bounding box of the face. - std::vector m_helper_indices; // The extra helper point indices when inserting outer ccb constraint. - std::size_t m_constraint_end; // The index past the last point in the outer ccb constraint. - double m_offset{0.5}; + std::vector m_points; + std::vector m_point_types; + std::vector> m_cst_ranges; + std::optional m_curr_cst_begin; + double m_offset{0}; }; #if defined(CGAL_DRAW_AOS_DEBUG) && defined(CGAL_DRAW_AOS_TRIANGULATOR_DEBUG_FILE_DIR) @@ -260,36 +283,40 @@ template void debug_print(const Arr_bounded_face_triangulator& triangulator) { const auto& ctx = triangulator.m_ctx; const auto& m_points = triangulator.m_points; - const auto& m_helper_indices = triangulator.m_helper_indices; + const auto& m_point_types = triangulator.m_point_types; + using Point_type = typename Arr_bounded_face_triangulator::Point_type; using Path = std::filesystem::path; Path debug_dir(CGAL_DRAW_AOS_TRIANGULATOR_DEBUG_FILE_DIR); std::string index_file_name = "index.txt"; Path index_file_path = debug_dir / index_file_name; - std::string ccb_file_name = "ccb_" + std::to_string(*ctx.debug_counter) + ".txt"; - std::string ccb_constraint_file_name = "ccb_constraint_" + std::to_string(*ctx.debug_counter) + ".txt"; - Path ccb_file_path = debug_dir / ccb_file_name; - Path ccb_constraint_file_path = debug_dir / ccb_constraint_file_name; + std::string points_file_name_prefix = "face_" + std::to_string(*ctx.debug_counter) + "_points"; + std::string ccb_constraint_file_name_prefix = "face_" + std::to_string(*ctx.debug_counter) + "_constraint"; const_cast(*ctx.debug_counter)++; std::ofstream ofs_index(index_file_path, std::ios::app); - ofs_index << ccb_file_name << "\n" << ccb_constraint_file_name << std::endl; - std::ofstream ofs_ccb(ccb_file_path); - std::size_t helper_indices_index = 0; - for(std::size_t i = 0; i < triangulator.m_constraint_end; ++i) { + auto points_filename = points_file_name_prefix + ".txt"; + auto points_path = debug_dir / points_filename; + std::ofstream ofs_points(points_path); + ofs_index << points_filename << std::endl; + for(int i = 0; i < triangulator.m_points.size(); ++i) { + if(m_point_types[i] == Point_type::Constraint_only) continue; const auto& pt = m_points[i]; - if(helper_indices_index < m_helper_indices.size() && i == m_helper_indices[helper_indices_index]) { - helper_indices_index++; - continue; - } - ofs_ccb << pt.x() << " " << pt.y() << "\n"; + ofs_points << pt.x() << " " << pt.y() << "\n"; } - std::ofstream ofs_ccb_constraint(ccb_constraint_file_path); - for(std::size_t i = 0; i < triangulator.m_constraint_end; ++i) { - const auto& pt = m_points[i]; - ofs_ccb_constraint << pt.x() << " " << pt.y() << "\n"; + int counter = 0; + for(auto [start_idx, end_idx] : triangulator.m_cst_ranges) { + auto filename = ccb_constraint_file_name_prefix + "_" + std::to_string(counter++) + ".txt"; + auto filepath = debug_dir / filename; + ofs_index << filename << std::endl; + std::ofstream ofs_ccb_constraint(filepath); + for(int i = start_idx; i < end_idx; ++i) { + if(m_point_types[i] == Point_type::Vertex_only) continue; + const auto& pt = m_points[i]; + ofs_ccb_constraint << pt.x() << " " << pt.y() << "\n"; + } } } #endif diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_renderer.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_renderer.h index 637cfad708b..4dedeae7c6b 100644 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_renderer.h +++ b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_renderer.h @@ -1,3 +1,18 @@ +// Copyright (c) 2025 +// Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), +// INRIA Sophia-Antipolis (France), +// Max-Planck-Institute Saarbruecken (Germany), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s): Shepard Liu + #ifndef CGAL_DRAW_AOS_ARR_BOUNDED_RENDERER_H #define CGAL_DRAW_AOS_ARR_BOUNDED_RENDERER_H @@ -5,7 +20,6 @@ #include #include #include -#include #include namespace CGAL { @@ -30,12 +44,12 @@ public: Approx_cache render() const { Approx_cache cache; if(m_ctx.is_cancelled()) return cache; - cache.reserve_faces(m_ctx.m_arr.number_of_faces()); - cache.reserve_halfedges(m_ctx.m_arr.number_of_halfedges()); - cache.reserve_vertices(m_ctx.m_arr.number_of_vertices()); + cache.vertices().reserve(m_ctx.m_arr.number_of_vertices()); + cache.halfedges().reserve(m_ctx.m_arr.number_of_halfedges()); + cache.faces().reserve(m_ctx.m_arr.number_of_faces()); - Arr_bounded_render_context ctx(m_ctx, m_bbox, cache); - Arr_bounded_approximate_face bounded_approx_face(ctx); + Arr_bounded_render_context derived_ctx(m_ctx, m_bbox, cache); + Arr_bounded_approximate_face bounded_approx_face(derived_ctx); for(Face_const_handle fh = m_ctx.m_arr.faces_begin(); fh != m_ctx.m_arr.faces_end(); ++fh) bounded_approx_face(fh); return cache; diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_coordinate_converter.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_coordinate_converter.h new file mode 100644 index 00000000000..98bbeba4405 --- /dev/null +++ b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_coordinate_converter.h @@ -0,0 +1,117 @@ +// Copyright (c) 2025 +// Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), +// INRIA Sophia-Antipolis (France), +// Max-Planck-Institute Saarbruecken (Germany), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s): Shepard Liu + +#ifndef CGAL_DRAW_AOS_ARR_COORDINATE_CONVERTER_H +#define CGAL_DRAW_AOS_ARR_COORDINATE_CONVERTER_H +#include + +#include +#include +#include + +namespace CGAL { +namespace draw_aos { + +/*! + * \brief class handling coordinate conversion between 2D parameterized surface coordinates and cartesian coordinates. + * + * \tparam GeomTraits + */ +template +class Arr_coordinate_converter +{ + using Geom_traits = GeomTraits; + using Approx_traits = Arr_approximate_traits; + using Approx_point = typename Approx_traits::Approx_point; + using Point = typename Approx_traits::Point; + +public: + Arr_coordinate_converter(const GeomTraits& traits) + : m_traits(traits) {} + + /*! + * \brief Converts a point in cartesian coordinates to parameterized surface coordinates. + * + * \param pt + * \return Point + */ + Point to_uv(Approx_point pt) const { return pt; } + + /*! + * \brief Converts a point in parameterized surface coordinates to cartesian coordinates. + * + * \param pt + * \return Approx_point + */ + Approx_point to_cartesian(Point pt) const { return pt; } + +private: + const GeomTraits& m_traits; +}; + +/*! + * \brief Converter specialization for geodesic arc on sphere traits. + * + * Provides conversions between spherical coordinates and right-handed Cartesian coordinates. Sphercial coordinates are + * represented as azimuth ( [0, 2 Pi) ) and polar ( [0, Pi] ) angle in radians. Points on the identification curve have + * azimuth == 0. The south pole has polar == 0. + * + * \tparam Kernel + * \tparam atanX + * \tparam atanY + */ +template +class Arr_coordinate_converter> +{ + using Geom_traits = Arr_geodesic_arc_on_sphere_traits_2; + using Approx_traits = Arr_approximate_traits; + using Approx_point = typename Approx_traits::Approx_point; + using Approx_nt = typename Approx_traits::Approx_nt; + using Point = typename Approx_traits::Point; + +public: + Arr_coordinate_converter(const Geom_traits& traits) + : m_traits(traits) {} + + Point to_uv(Approx_point point) const { + if(point.location() == Approx_point::MAX_BOUNDARY_LOC) return Point(0, CGAL_PI); + if(point.location() == Approx_point::MIN_BOUNDARY_LOC) return Point(0, 0); + Approx_nt azimuth_from_id = + std::fmod(std::atan2(point.dy(), point.dx()) - std::atan2(atanY, atanX) + 2 * CGAL_PI, 2 * CGAL_PI); + return Point(azimuth_from_id, std::acos(-point.dz())); + } + + Approx_point to_cartesian(Point point) const { + using Direction_3 = typename Geom_traits::Approximate_kernel::Direction_3; + + Approx_nt polar = point.y(); + if(point.y() == CGAL_PI) return Approx_point(Direction_3(0, 0, 1), Approx_point::MAX_BOUNDARY_LOC); + if(point.y() == 0) return Approx_point(Direction_3(0, 0, -1), Approx_point::MIN_BOUNDARY_LOC); + Approx_nt azimuth = point.x() + std::atan2(atanY, atanX); + Approx_nt x = std::sin(polar) * std::cos(azimuth); + Approx_nt y = std::sin(polar) * std::sin(azimuth); + Approx_nt z = -std::cos(polar); + Direction_3 dir(x, y, z); + return Approx_point(dir, azimuth == 0 ? Approx_point::MID_BOUNDARY_LOC : Approx_point::NO_BOUNDARY_LOC); + } + +private: + const Geom_traits& m_traits; +}; + +} // namespace draw_aos +} // namespace CGAL + +#endif \ No newline at end of file diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_face_point_generator.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_face_point_generator.h new file mode 100644 index 00000000000..dea43538151 --- /dev/null +++ b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_face_point_generator.h @@ -0,0 +1,101 @@ +// Copyright (c) 2025 +// Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), +// INRIA Sophia-Antipolis (France), +// Max-Planck-Institute Saarbruecken (Germany), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s): Shepard Liu + +#ifndef CGAL_DRAW_AOS_ARR_FACE_POINT_GENERATOR_H +#define CGAL_DRAW_AOS_ARR_FACE_POINT_GENERATOR_H +#include +#include +#include + +#include + +#include "CGAL/unordered_flat_map.h" +#include "CGAL/Arr_batched_point_location.h" +#include "CGAL/Arr_point_location_result.h" +#include "CGAL/Draw_aos/Arr_coordinate_converter.h" +#include "CGAL/Draw_aos/type_utils.h" + +namespace CGAL { +namespace draw_aos { + +/*! + * \brief Generate face interior points. + * + * \tparam Arrangement + */ +template +class Arr_face_point_generator; + +template +class Arr_face_point_generator< + Arrangement, + std::enable_if_t>> +{ + using Point_geom = typename Arr_approximate_traits::Point; + using Face_const_handle = typename Arrangement::Face_const_handle; + +public: + using Face_points_map = unordered_flat_map>; + + // No-op implementation for non-curved surface arrangements. + Face_points_map operator()(const Arrangement&, double) { return {}; } +}; + +template +class Arr_face_point_generator>> +{ + using Geom_traits = typename Arrangement::Geometry_traits_2; + using Approx_traits = Arr_approximate_traits; + using Approx_nt = typename Approx_traits::Approx_nt; + using Point = typename Approx_traits::Point; + using Face_const_handle = typename Arrangement::Face_const_handle; + using Gt_point = typename Geom_traits::Point_2; + using Query_result = std::pair::Type>; + +public: + using Face_points_map = unordered_flat_map>; + + Face_points_map operator()(const Arrangement& arr, double error) { + const Geom_traits& traits = *arr.geometry_traits(); + + // Grid sampling in parameter space. + Approx_nt cell_size = 2.0 * std::acos(1 - error); + std::vector points; + Arr_coordinate_converter coords(traits); + points.reserve(2 * CGAL_PI / cell_size * CGAL_PI / cell_size); + for(Approx_nt x = 0; x < 2 * CGAL_PI; x += cell_size) { + for(Approx_nt y = 0; y < CGAL_PI; y += cell_size) { + auto pt = coords.to_cartesian(Point(x, y)); + points.push_back(traits.construct_point_2_object()(pt.dx(), pt.dy(), pt.dz())); + } + } + + unordered_flat_map> face_points; + CGAL::locate(arr, points.begin(), points.end(), + boost::make_function_output_iterator([&face_points, &traits, &coords](const Query_result& res) { + if(!std::holds_alternative(res.second)) return; + Face_const_handle fh = std::get(res.second); + auto [it, _] = face_points.try_emplace(fh, std::vector()); + it->second.push_back(coords.to_uv(traits.approximate_2_object()(res.first))); + })); + return face_points; + } +}; + +} // namespace draw_aos +} // namespace CGAL + +#endif \ No newline at end of file diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_graph_conn.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_graph_conn.h deleted file mode 100644 index 39dc75ccc46..00000000000 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_graph_conn.h +++ /dev/null @@ -1,85 +0,0 @@ - -#ifndef CGAL_DRAW_AOS_ARR_GRAPH_CONN_H -#define CGAL_DRAW_AOS_ARR_GRAPH_CONN_H -#include - -#include -#include -#include - -namespace CGAL { -/** - * @brief Arr_graph_conn provides fast connectivity queries for arrangement vertices - * based on union-find data structure. - */ -template -class Arr_graph_conn -{ - using Vertex_const_handle = typename Arr::Vertex_const_handle; - using Edge_const_handle = typename Arr::Edge_const_iterator; - using Halfedge_const_handle = typename Arr::Halfedge_const_iterator; - using Union_find_handle = typename Union_find::handle; - -private: - void insert_halfedge(Halfedge_const_handle he) { - const auto& source = he->source(); - const auto& target = he->target(); - auto [source_it, source_inserted] = m_lookup.try_emplace(source, Union_find_handle()); - if(source_inserted) source_it->second = m_uf.make_set(source); - auto [target_it, target_inserted] = m_lookup.try_emplace(target, Union_find_handle()); - if(target_inserted) target_it->second = m_uf.make_set(target); - m_uf.unify_sets(source_it->second, target_it->second); - } - - void insert_isolated_vertex(const Vertex_const_handle& vh) { m_lookup[vh] = m_uf.make_set(vh); } - -public: - Arr_graph_conn(const Arr& arr) { - m_lookup.reserve(arr.number_of_vertices()); - for(const auto& he : arr.halfedge_handles()) { - if(he->direction() != ARR_LEFT_TO_RIGHT) continue; - insert_halfedge(he); - } - for(const auto& vh : arr.vertex_handles()) { - if(!vh->is_isolated()) continue; - insert_isolated_vertex(vh); - } - // add fictitious edges and open vertices - for(auto fh = arr.unbounded_faces_begin(); fh != arr.unbounded_faces_end(); ++fh) { - if(!fh->has_outer_ccb()) continue; - auto outer_ccb = fh->outer_ccb(); - auto curr = outer_ccb; - do { - if(!curr->is_fictitious()) continue; - insert_halfedge(curr); - } while(++curr != outer_ccb); - } - } - - bool is_connected(const Vertex_const_handle& v1, const Vertex_const_handle& v2) const { - auto it1 = m_lookup.find(v1); - auto it2 = m_lookup.find(v2); - // This can't happen - CGAL_assertion(it1 != m_lookup.end() && it2 != m_lookup.end()); - return m_uf.same_set(it1->second, it2->second); - } - - /** - * @brief Returns the representative vertex of the connected component containing the given vertex. - * For each connected component boundary (CCB), the same representative vertex is consistently returned. - * @param vh a vertex handle in the ccb - * @return Vertex_const_handle - */ - Vertex_const_handle ccb_representative_vertex(const Vertex_const_handle& vh) const { - // path compression typically does not mutate the representative member of a set. - auto it = m_lookup.find(vh); - CGAL_assertion(it != m_lookup.end()); - return *m_uf.find(it->second); - } - -private: - Union_find m_uf; - unordered_flat_map m_lookup; -}; -} // namespace CGAL -#endif // CGAL_DRAW_AOS_ARR_GRAPH_CONN_H \ No newline at end of file diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_portals.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_portals.h deleted file mode 100644 index 955b2405baf..00000000000 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_portals.h +++ /dev/null @@ -1,96 +0,0 @@ -#ifndef CGAL_DRAW_AOS_ARR_CREATE_PORTALS_H -#define CGAL_DRAW_AOS_ARR_CREATE_PORTALS_H -#include -#include - -#include - -#include -#include -#include -#include -#include - -namespace CGAL { -namespace draw_aos { - -/** - * @brief Portals are virtual vertical segments that connect the outer - * connected component boundary (ccb) of a face with its inner ccbs. - * - * We use vertical decomposition to cast upward rays in a batched manner. - * For each inner ccb, only one "portal" is created at a vertex whose upward ray - * hits another ccb. Eventually, faces of the arrangement become hole-free and can - * be drawn with Graphics_scene. - */ -template -class Arr_portals -{ - using Geom_traits = typename Arrangement::Geometry_traits_2; - using Vertex_const_handle = typename Arrangement::Vertex_const_handle; - using Halfedge_const_handle = typename Arrangement::Halfedge_const_handle; - using Face_const_handle = typename Arrangement::Face_const_handle; - using Point_2 = typename Geom_traits::Point_2; - using X_monotone_curve_2 = typename Geom_traits::X_monotone_curve_2; - using Approx_point = typename Arr_approximate_traits::Approx_point; - using Feature_const = std::variant; - -public: - using Portal_exit = Vertex_const_handle; - using Portal_exit_vector = std::vector; - // Map from a feature to its portals sorted by the x coordinate of the virtual vertical segments. - using Feature_portals_map = unordered_flat_map; - -public: - Arr_portals(const Arrangement& arr) { init(arr); } - -private: - void init(const Arrangement& arr) { - using Object_pair = std::pair; - using Vert_decomp_entry = std::pair; - - Arr_graph_conn conn(arr); - auto visited_ccbs = std::unordered_set(); - - auto func_out_iter = boost::make_function_output_iterator([&](const Vert_decomp_entry& entry) { - const auto& [vh, obj_pair] = entry; - const auto& above_feat = obj_pair.second; - Vertex_const_handle ccb_main_vertex = conn.ccb_representative_vertex(vh); - // Skip processed ccb. - if(visited_ccbs.find(ccb_main_vertex) != visited_ccbs.end()) return; - - if(Vertex_const_handle above_vh; CGAL::assign(above_vh, above_feat)) { - // Skip vertex connected to vh - if(conn.is_connected(above_vh, vh)) return; - const auto& [it, _] = m_map.try_emplace(above_vh, Portal_exit_vector{}); - it->second.emplace_back(vh); - } else if(Halfedge_const_handle above_he; CGAL::assign(above_he, above_feat)) { - // The given halfedge is always directed from right to left (exactly what we need). - if(conn.is_connected((above_he)->source(), vh)) return; - const auto& [it, _] = m_map.try_emplace(above_he, Portal_exit_vector{}); - it->second.emplace_back(vh); - } else // Don't handle faces. - return; - - visited_ccbs.insert(ccb_main_vertex); - }); - - decompose(arr, func_out_iter); - - // reverse portals on each halfedge, as they are stored in left-to-right order - for(auto& [feat, portals] : m_map) { - if(!std::holds_alternative(feat)) continue; - std::reverse(portals.begin(), portals.end()); - } - } - -public: - const Feature_portals_map& get() const { return m_map; } - -private: - Feature_portals_map m_map; -}; - -} // namespace draw_aos -} // namespace CGAL -#endif \ No newline at end of file diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_projector.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_projector.h deleted file mode 100644 index 9e2519bceee..00000000000 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_projector.h +++ /dev/null @@ -1,88 +0,0 @@ -#ifndef CGAL_DRAW_AOS_ARR_PROJECTION_H -#define CGAL_DRAW_AOS_ARR_PROJECTION_H -#include - -#include -#include -#include - -namespace CGAL { -namespace draw_aos { - -/** - * @brief class handling projection between 2D parameter space and coordinate space - * - * @tparam GeomTraits - */ -template -class Arr_projector -{ - using Geom_traits = GeomTraits; - using Approx_traits = Arr_approximate_traits; - using Approx_point = typename Approx_traits::Approx_point; - using Approx_proj_point = typename Approx_traits::Approx_proj_point; - -public: - Arr_projector(const GeomTraits& traits) - : m_traits(traits) {} - - Approx_proj_point project(Approx_point pt) const { return pt; } - - Approx_point unproject(Approx_proj_point pt) const { return pt; } - -private: - const GeomTraits& m_traits; -}; - -/** - * @brief Projector specialization for geodesic arc on sphere traits. - * - * The projection process is essentially a conversion between spherical coordinates and right-handed Cartesian - * coordinates. Sphercial coordinates are represented as azimuth and polar angle. Note that the polar angle starts from - * north pole. - * - * @tparam Kernel - * @tparam atanX - * @tparam atanY - */ -template -class Arr_projector> -{ - using Geom_traits = Arr_geodesic_arc_on_sphere_traits_2; - using Approx_traits = Arr_approximate_traits; - using Approx_point = typename Approx_traits::Approx_point; - using Approx_nt = typename Approx_traits::Approx_nt; - using Approx_proj_point = typename Approx_traits::Approx_proj_point; - -public: - Arr_projector(const Geom_traits& traits) - : m_traits(traits) {} - - Approx_proj_point project(Approx_point point) const { - if(point.location() == Approx_point::MAX_BOUNDARY_LOC) return Approx_proj_point(0, CGAL_PI); - if(point.location() == Approx_point::MIN_BOUNDARY_LOC) return Approx_proj_point(0, 0); - Approx_nt azimuth_from_id = - std::fmod(std::atan2(point.dy(), point.dx()) - std::atan2(atanY, atanX) + 2 * CGAL_PI, 2 * CGAL_PI); - return Approx_proj_point(azimuth_from_id, std::acos(-point.dz())); - } - - Approx_point unproject(Approx_proj_point point) const { - using Direction_3 = typename Geom_traits::Approximate_kernel::Direction_3; - - Approx_nt polar = point.y(); - if(point.y() == CGAL_PI) return Approx_point(Direction_3(0, 0, 1), Approx_point::MAX_BOUNDARY_LOC); - if(point.y() == 0) return Approx_point(Direction_3(0, 0, -1), Approx_point::MIN_BOUNDARY_LOC); - Approx_nt azimuth = point.x() + std::atan2(atanY, atanX); - return Approx_point( - Direction_3(std::sin(polar) * std::cos(azimuth), std::sin(polar) * std::sin(azimuth), -std::cos(polar)), - azimuth == 0 ? Approx_point::MID_BOUNDARY_LOC : Approx_point::NO_BOUNDARY_LOC); - } - -private: - const Geom_traits& m_traits; -}; - -} // namespace draw_aos -} // namespace CGAL - -#endif \ No newline at end of file diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_render_context.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_render_context.h index c2522ae9573..144b97fb3bc 100644 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_render_context.h +++ b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_render_context.h @@ -1,3 +1,18 @@ +// Copyright (c) 2025 +// Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), +// INRIA Sophia-Antipolis (France), +// Max-Planck-Institute Saarbruecken (Germany), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s): Shepard Liu + #ifndef CGAL_DRAW_AOS_ARR_RENDER_CONTEXT_H #define CGAL_DRAW_AOS_ARR_RENDER_CONTEXT_H #include @@ -10,8 +25,9 @@ #include #include #include -#include #include +#include +#include #if defined(CGAL_DRAW_AOS_DEBUG) #include @@ -65,7 +81,7 @@ class Arr_bounds_context_mixin { using Geom_traits = GeomTraits; using Approx_traits = Arr_approximate_traits; - using Point_geom = typename Approx_traits::Point_geom; + using Point = typename Approx_traits::Point; using Approx_nt = typename Approx_traits::Approx_nt; protected: @@ -81,32 +97,43 @@ public: bool contains_x(Approx_nt x) const { return xmin() <= x && x <= xmax(); } bool contains_y(Approx_nt y) const { return ymin() <= y && y <= ymax(); } - bool contains(Point_geom pt) const { return contains_x(pt.x()) && contains_y(pt.y()); } + bool contains(Point pt) const { return contains_x(pt.x()) && contains_y(pt.y()); } - bool is_on_boundary(Point_geom pt) const { - return (pt.x() == xmin() || pt.x() == xmax()) && contains_y(pt.y()) || - (pt.y() == ymin() || pt.y() == ymax()) && contains_x(pt.x()); - } + Point top_left() const { return Point(xmin(), ymax()); } + Point top_right() const { return Point(xmax(), ymax()); } + Point bottom_left() const { return Point(xmin(), ymin()); } + Point bottom_right() const { return Point(xmax(), ymin()); } + + bool is_on_left(Point pt) const { return pt.x() == xmin() && contains_y(pt.y()); } + bool is_on_right(Point pt) const { return pt.x() == xmax() && contains_y(pt.y()); } + bool is_on_bottom(Point pt) const { return pt.y() == ymin() && contains_x(pt.x()); } + bool is_on_top(Point pt) const { return pt.y() == ymax() && contains_x(pt.x()); } + bool is_on_boundary(Point pt) const { return is_on_left(pt) || is_on_right(pt) || is_on_bottom(pt) || is_on_top(pt); } private: const Bbox_2 m_bbox; }; +template +using Arr_parameterization_context_mixin = Arr_coordinate_converter; + template -class Arr_render_context : public Arr_cancellable_context_mixin +class Arr_render_context : public Arr_cancellable_context_mixin, + public Arr_parameterization_context_mixin { using Cancellable_context_mixin = Arr_cancellable_context_mixin; + using Param_context_mixin = Arr_parameterization_context_mixin; using Geom_traits = typename Arrangement::Geometry_traits_2; - using Portals = Arr_portals; - using Feature_portals_map = typename Portals::Feature_portals_map; + using Face_points_map = typename Arr_face_point_generator::Face_points_map; public: - Arr_render_context(const Arrangement& arr, const Feature_portals_map& portals_map, double approx_error) + Arr_render_context(const Arrangement& arr, double approx_error, Face_points_map& face_points) : Cancellable_context_mixin() + , Param_context_mixin(*arr.geometry_traits()) , m_arr(arr) , m_traits(*arr.geometry_traits()) - , m_portals_map(portals_map) - , m_approx_error(approx_error) { + , m_approx_error(approx_error) + , m_face_points(face_points) { #if defined(CGAL_DRAW_AOS_DEBUG) && defined(CGAL_DRAW_AOS_TRIANGULATOR_DEBUG_FILE_DIR) std::filesystem::path debug_file_dir(CGAL_DRAW_AOS_TRIANGULATOR_DEBUG_FILE_DIR); // clear the index file. @@ -118,7 +145,7 @@ public: const double m_approx_error; const Arrangement& m_arr; const Geom_traits& m_traits; - const Feature_portals_map& m_portals_map; + const Face_points_map& m_face_points; #if defined(CGAL_DRAW_AOS_DEBUG) std::shared_ptr debug_counter = std::make_shared(0); diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_viewer.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_viewer.h index 1527933dd87..5be7a151bd3 100644 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_viewer.h +++ b/Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_viewer.h @@ -16,15 +16,11 @@ #ifndef ARR_VIEWER_H #define ARR_VIEWER_H +#include #include -#include #include -#include #include -#include -#include - #include #include #include @@ -38,98 +34,75 @@ #include #include #include -#include #include #include #include #include -#include +#include #include -#include -#include #include +#include +#include +#include namespace CGAL { namespace draw_aos { -/** - * @brief Viewer for visualizing arrangements on surface. +/*! + * \brief Viewport helper functions * - * @tparam Arrangement should be parameterized by a geometry traits that models Approximate_2, or it silently fails. - * @tparam GSOptions Graphics scene options + * \tparam Arrangement */ -template -class Arr_viewer; +template +class Arr_viewport_helpers; -template -class Arr_viewer>> - : public Qt::Basic_viewer +// Specialization for planar arrangements +template +class Arr_viewport_helpers< + Arrangement, + std::enable_if_t>> { -public: - Arr_viewer(QWidget* parent, const Arrangement& arr, GSOptions options, const char* title = "Arrangement Viewer") - : Qt::Basic_viewer(parent, Graphics_scene(), title) { - std::cerr << "Arr_viewer requires geometry traits that supports approximation of Point_2 and " - "X_monotone_curve_2." - << std::endl; - exit(1); - } - - virtual ~Arr_viewer() = default; -}; - -template -class Arr_viewer>> - : public Qt::Basic_viewer -{ - using Basic_viewer = Qt::Basic_viewer; - using Vertex_const_handle = typename Arrangement::Vertex_const_handle; - using Halfedge_const_handle = typename Arrangement::Halfedge_const_handle; - using Face_const_handle = typename Arrangement::Face_const_handle; - using Feature_portal_map = typename Arr_portals::Feature_portals_map; using Geom_traits = typename Arrangement::Geometry_traits_2; using Approx_traits = Arr_approximate_traits; using Approx_point = typename Approx_traits::Approx_point; - using Point_geom = typename Approx_traits::Point_geom; + using Camera = qglviewer::Camera; + using Point = typename Approx_traits::Point; + using Local_point = Buffer_for_vao::Local_point; - struct Render_params - { - bool operator==(const Render_params& other) const { - return bbox == other.bbox && approx_error == other.approx_error; - } - Bbox_2 bbox; - double approx_error{0}; - }; +protected: + Arr_viewport_helpers(const Arrangement& arr) + : m_arr(arr) {} -private: - static bool contains(const Bbox_2& bbox, const Point_geom& pt) { - return bbox.xmin() <= pt.x() && pt.x() <= bbox.xmax() && bbox.ymin() <= pt.y() && pt.y() <= bbox.ymax(); + /*! + * \brief Computes a subpixel-level approximation error based on the bounding box and viewport width. + * + * \param bbox + * \param viewport_width width of the viewport in pixels + * \return double + */ + double approximation_error(const Bbox_2& bbox, int viewport_width) const { + return bbox.x_span() / viewport_width; } - /** - * @brief Computes the initial bounding box of the arrangement. + /*! + * \brief Computes a parameter space bounding box that contains everything in the arrangement with some margin. * - * @return Bbox_2 + * \note For arrangement induced by unbounded curves, the bounding box only fits all vertices. + * \return Bbox_2 */ - Bbox_2 initial_bbox() const { return initial_bbox_impl(); } - - template , int> = 0> - Bbox_2 initial_bbox_impl() const { + Bbox_2 arr_bbox() const { const auto& traits = *m_arr.geometry_traits(); Bbox_2 bbox; // Computes a rough bounding box from the vertices. for(const auto& vh : m_arr.vertex_handles()) { - bbox += m_proj.project(traits.approximate_2_object()(vh->point())).bbox(); + bbox += traits.approximate_2_object()(vh->point()).bbox(); } - double approx_error = get_approx_error(bbox); + double approx_error = approximation_error(bbox, 100); // Computes a more precise bounding box from the halfedges. for(const auto& he : m_arr.halfedge_handles()) { - traits.approximate_2_object()(he->curve(), approx_error, - boost::make_function_output_iterator( - [this, &bbox](Approx_point pt) { bbox += this->m_proj.project(pt).bbox(); })); + traits.approximate_2_object()( + he->curve(), approx_error, + boost::make_function_output_iterator([this, &bbox](Approx_point pt) { bbox += pt.bbox(); })); } // Place margin around the bbox. double dx = bbox.x_span() * 0.1; @@ -141,22 +114,26 @@ private: return bbox; } - // Specialization for geodesic arc on sphere traits - template , int> = 0> - Bbox_2 initial_bbox_impl() const { - return Bbox_2(0, 0, 2 * CGAL_PI, CGAL_PI); + /*! + * \brief Fits the camera to bbox. + * + * \param bbox + * \param camera + */ + void fit_camera(const Bbox_2& bbox, Camera& cam) const { + using Vec = qglviewer::Vec; + cam.fitBoundingBox(Vec(bbox.xmin(), bbox.ymin(), 0.0), Vec(bbox.xmax(), bbox.ymax(), 0.0)); } - /** - * @brief Computes the corresponding world bounding box of viewport. - * @return Bbox_2 + /*! + * \brief Computes parameter space axis aligned bounding box from camera parameters. + * + * \param cam + * \return Bbox_2 */ - Bbox_2 view_bbox_from_camera() const { return view_bbox_from_camera_impl(); } - - template , int> = 0> - Bbox_2 view_bbox_from_camera_impl() const { + Bbox_2 screen_to_world(const Camera& cam) const { QMatrix4x4 mvp; - camera_->getModelViewProjectionMatrix(mvp.data()); + cam.getModelViewProjectionMatrix(mvp.data()); QMatrix4x4 inverse_mvp = mvp.inverted(); // Define 4 corners of the near plane in NDC (-1 to 1 in x and y) std::array clip_space_corners{QVector4D(-1.0, -1.0, 0.0, 1.0), QVector4D(-1.0, 1.0, 0.0, 1.0), @@ -178,80 +155,45 @@ private: return Bbox_2(xmin, ymin, xmax, ymax); } - // Specialization for geodesic arc on sphere traits - template , int> = 0> - Bbox_2 view_bbox_from_camera_impl() const { - // We fix the bounding box for spherical arrangements to cover the whole sphere. - return Bbox_2(0, 0, 2 * CGAL_PI, CGAL_PI); - } - - /** - * @brief Converts a Point_geom point to euclidean point that can be handled by Buffer_for_vao. + /*! + * \brief Converts a parameter space point to a local point of the buffer object. * - * @param pt - * @return Buffer_for_vao::Local_point + * \param pt + * \return Local_point */ - Buffer_for_vao::Local_point transform_point(const Point_geom& proj_pt) const { - return transform_point_impl(proj_pt); - } + Local_point to_local_point(Point pt) const { return Local_point(pt.x(), pt.y(), 0.0); } - template , int> = 0> - Buffer_for_vao::Local_point transform_point_impl(Point_geom proj_pt) const { - return Buffer_for_vao::Local_point(proj_pt.x(), proj_pt.y(), 0.0); - } +private: + const Arrangement& m_arr; +}; - // Specialization for geodesic arc on sphere traits - template , int> = 0> - Buffer_for_vao::Local_point transform_point_impl(Point_geom proj_pt) const { - auto approx_pt = m_proj.unproject(proj_pt); - return Buffer_for_vao::Local_point(approx_pt.dx(), approx_pt.dy(), approx_pt.dz()); - } +// Spherical arrangement specialization +template +class Arr_viewport_helpers>> +{ + using Geom_traits = typename Arrangement::Geometry_traits_2; + using Approx_traits = Arr_approximate_traits; + using Approx_point = typename Approx_traits::Approx_point; + using Camera = qglviewer::Camera; + using Point = typename Approx_traits::Point; + using Local_point = Buffer_for_vao::Local_point; - /** - * @brief Fits the camera to the given bounding box. - * - * Works only for orthogonal camera directed onto a 2D plane by now. - * @param bbox - */ - void fit_camera(const Bbox_2& bbox) { fit_camera_impl(bbox); } +protected: + Arr_viewport_helpers(const Arrangement& arr) + : m_arr(arr) {} - template , int> = 0> - void fit_camera_impl(const Bbox_2& bbox) { - // CGAL_assertion(camera_->type() == qglviewer::Camera::ORTHOGRAPHIC); + Bbox_2 arr_bbox() const { return Bbox_2(0, 0, 2 * CGAL_PI, CGAL_PI); } + + Bbox_2 screen_to_world(const Camera& cam) const { return Bbox_2(0, 0, 2 * CGAL_PI, CGAL_PI); } + + void fit_camera(const Bbox_2&, Camera& cam) { using Vec = qglviewer::Vec; - camera_->fitBoundingBox(Vec(bbox.xmin(), bbox.ymin(), 0.0), Vec(bbox.xmax(), bbox.ymax(), 0.0)); + cam.setSceneCenter(Vec(0, 0, 0)); + cam.fitSphere(Vec(0, 0, 0), 1.1); // slightly larger than the unit sphere } - // Specialization for geodesic arc on sphere traits - template , int> = 0> - void fit_camera_impl(const Bbox_2&) { - using Vec = qglviewer::Vec; - camera_->setSceneCenter(Vec(0, 0, 0)); - camera_->fitSphere(Vec(0, 0, 0), 1.1); // slightly larger than the sphere - } - - /** - * @brief Computes a proper approximation error bound. - * - * @param bbox 2D parameter space bounding box - * @return double - */ - double get_approx_error(const Bbox_2& bbox) const { return get_approx_error_impl(bbox); } - - template , int> = 0> - double get_approx_error_impl(const Bbox_2& bbox) const { - std::array viewport; - camera_->getViewport(viewport.data()); - double viewport_width = static_cast(viewport[2]); - return bbox.x_span() / viewport_width; - } - - // Specialization for geodesic arc on sphere traits - template , int> = 0> - double get_approx_error_impl(const Bbox_2& bbox) const { - std::array viewport; - camera_->getViewport(viewport.data()); - double viewport_width = static_cast(viewport[2]); + double approximation_error(const Bbox_2& bbox, int viewport_width) const { // If crossing hemisphere if(bbox.x_span() >= CGAL_PI) return 1.0 / viewport_width; // Otherwise we evalute the error bound with respect to the longest longitude arc @@ -260,89 +202,143 @@ private: return bbox.x_span() * std::sin(theta) / viewport_width; } - /** - * @brief Computes the render parameters from camera states. - */ - Render_params get_arr_render_params() { + Buffer_for_vao::Local_point to_local_point(Point pt) const { + auto approx_pt = Arr_coordinate_converter(*m_arr.geometry_traits()).to_cartesian(pt); + return Buffer_for_vao::Local_point(approx_pt.dx(), approx_pt.dy(), approx_pt.dz()); + } + +private: + const Arrangement& m_arr; +}; + +/*! Viewer for visualizing arrangements on surface. + * + * \tparam Arrangement + * \tparam GSOptions + */ +template +class Arr_viewer : public Qt::Basic_viewer, Arr_viewport_helpers +{ + using Basic_viewer = Qt::Basic_viewer; + using Helpers = Arr_viewport_helpers; + using Vertex_const_handle = typename Arrangement::Vertex_const_handle; + using Halfedge_const_handle = typename Arrangement::Halfedge_const_handle; + using Face_const_handle = typename Arrangement::Face_const_handle; + using Geom_traits = typename Arrangement::Geometry_traits_2; + using Approx_traits = Arr_approximate_traits; + using Approx_point = typename Approx_traits::Approx_point; + using Point = typename Approx_traits::Point; + using Point_generator = Arr_face_point_generator; + using Faces_point_map = typename Point_generator::Face_points_map; + + struct Render_params + { + bool operator==(const Render_params& other) const { + return bbox == other.bbox && approx_error == other.approx_error; + } + Bbox_2 bbox; + double approx_error{0}; + }; + + constexpr static bool Is_on_curved_surface = is_or_derived_from_curved_surf_traits_v; + +private: + static bool contains(const Bbox_2& bbox, const Point& pt) { + return bbox.xmin() <= pt.x() && pt.x() <= bbox.xmax() && bbox.ymin() <= pt.y() && pt.y() <= bbox.ymax(); + } + + int viewport_width() const { + std::array viewport; + this->camera_->getViewport(viewport.data()); + return viewport[2]; + } + + Render_params compute_render_params() { Render_params params; - params.bbox = view_bbox_from_camera(); - params.approx_error = get_approx_error(params.bbox); + params.bbox = this->screen_to_world(*this->camera_); + params.approx_error = this->approximation_error(params.bbox, viewport_width()); return params; } - /** - * @brief Render arrangement within the given bounding box. - * - * @param bbox - */ void render_arr(const Render_params& params) { const Bbox_2& bbox = params.bbox; - Arr_render_context ctx(m_arr, m_portals.get(), params.approx_error); + auto face_points = Point_generator()(m_arr, params.approx_error); + Arr_render_context ctx(m_arr, params.approx_error, face_points); Arr_bounded_renderer renderer(ctx, bbox); auto cache = renderer.render(); + // add faces - for(const auto& [fh, face_tris] : cache.faces()) { - const auto& points = face_tris.points; - const auto& tris = face_tris.triangles; - bool draw_face = m_gso.colored_face(m_arr, fh); - for(const auto& t : tris) { - if(draw_face) - m_gs.face_begin(m_gso.face_color(m_arr, fh)); + for(const auto& [fh, tf] : cache.faces()) { + if(!m_gso.draw_face(m_arr, fh)) continue; + bool colored_face = m_gso.colored_face(m_arr, fh); + auto color = colored_face ? m_gso.face_color(m_arr, fh) : CGAL::IO::Color(); + for(const auto& tri : tf.triangles) { + if(colored_face) + m_gs.face_begin(color); else m_gs.face_begin(); - for(const auto idx : t) m_gs.add_point_in_face(transform_point(points[idx])); + for(const auto i : tri) m_gs.add_point_in_face(this->to_local_point(tf.points[i])); m_gs.face_end(); } } // add edges for(const auto& [he, polyline] : cache.halfedges()) { - if(polyline.size() < 2) continue; - bool draw_colored_edge = m_gso.colored_edge(m_arr, he); - auto color = draw_colored_edge ? m_gso.edge_color(m_arr, he) : CGAL::IO::Color(); - for(size_t i = 0; i < polyline.size() - 1; ++i) { - const auto& cur_pt = polyline[i]; - const auto& next_pt = polyline[i + 1]; - auto mid_pt = CGAL::midpoint(cur_pt, next_pt); - if(!contains(bbox, mid_pt)) continue; - if(draw_colored_edge) - m_gs.add_segment(transform_point(cur_pt), transform_point(next_pt), color); + if(he->direction() == ARR_RIGHT_TO_LEFT || !m_gso.draw_edge(m_arr, he) || polyline.size() < 2) continue; + bool colored_edge = m_gso.colored_edge(m_arr, he); + auto color = colored_edge ? m_gso.edge_color(m_arr, he) : CGAL::IO::Color(); + // skip first two if starts with a sep point. + int start_idx = Approx_traits::is_null(polyline.front()) ? 2 : 0; + // skip last two if ends with a sep point. + int end_idx = Approx_traits::is_null(polyline.back()) ? polyline.size() - 2 : polyline.size(); + for(int i = start_idx; i < end_idx - 1; ++i) { + const auto& src = polyline[i]; + const auto& tgt = polyline[i + 1]; + if(Approx_traits::is_null(src) || Approx_traits::is_null(tgt)) continue; + if(!contains(bbox, src) || !contains(bbox, tgt)) continue; + if(colored_edge) + m_gs.add_segment(this->to_local_point(src), this->to_local_point(tgt), color); else - m_gs.add_segment(transform_point(cur_pt), transform_point(next_pt)); + m_gs.add_segment(this->to_local_point(src), this->to_local_point(tgt)); } } // add vertices for(const auto& [vh, pt] : cache.vertices()) { - if(!contains(bbox, pt)) continue; + if(!m_gso.draw_vertex(m_arr, vh) || !contains(bbox, pt)) continue; if(m_gso.colored_vertex(m_arr, vh)) - m_gs.add_point(transform_point(pt), m_gso.vertex_color(m_arr, vh)); + m_gs.add_point(this->to_local_point(pt), m_gso.vertex_color(m_arr, vh)); else - m_gs.add_point(transform_point(pt)); + m_gs.add_point(this->to_local_point(pt)); } } - /** - * @brief Rerender scene within the given bounding box. + /*! + * \brief Rerender scene within the given bounding box. - * @param bbox + * \param bbox */ void rerender(const Render_params& params) { - if(params == m_last_render_params) return; - m_last_render_params = params; + if(params == m_last_params) return; + m_last_params = params; m_gs.clear(); render_arr(params); Basic_viewer::redraw(); } public: - Arr_viewer(QWidget* parent, const Arrangement& arr, GSOptions gso, const char* title = "Arrangement Viewer") + Arr_viewer(QWidget* parent, const Arrangement& arr, const GSOptions& gso, const char* title, Bbox_2 initial_bbox) : Basic_viewer(parent, m_gs, title) + , Helpers(arr) , m_gso(gso) , m_arr(arr) - , m_proj(*arr.geometry_traits()) - , m_portals(arr) {} + , m_coords(*arr.geometry_traits()) { + if(initial_bbox.x_span() == 0 || initial_bbox.y_span() == 0 || Is_on_curved_surface) + m_initial_bbox = this->arr_bbox(); + else + m_initial_bbox = initial_bbox; + } virtual void draw() override { - Render_params params = get_arr_render_params(); + Render_params params = compute_render_params(); #if defined(CGAL_DRAW_AOS_DEBUG) if constexpr(!is_or_derived_from_agas_v) { @@ -359,7 +355,7 @@ public: if(!m_initialized) { // The initial render must be done with original camera parameters or the width of edges gets exaggerated. // So we fit the camera after initial render. - fit_camera(initial_bbox()); + this->fit_camera(m_initial_bbox, *this->camera_); m_initialized = true; } @@ -372,10 +368,10 @@ private: Graphics_scene m_gs; GSOptions m_gso; const Arrangement& m_arr; - const Arr_portals m_portals; bool m_initialized{false}; - const Arr_projector m_proj; - Render_params m_last_render_params; + Bbox_2 m_initial_bbox; + const Arr_coordinate_converter m_coords; + Render_params m_last_params; }; } // namespace draw_aos diff --git a/Arrangement_on_surface_2/include/CGAL/Draw_aos/type_utils.h b/Arrangement_on_surface_2/include/CGAL/Draw_aos/type_utils.h index aa0c6b42366..c2d421c657f 100644 --- a/Arrangement_on_surface_2/include/CGAL/Draw_aos/type_utils.h +++ b/Arrangement_on_surface_2/include/CGAL/Draw_aos/type_utils.h @@ -1,14 +1,31 @@ +// Copyright (c) 2025 +// Utrecht University (The Netherlands), +// ETH Zurich (Switzerland), +// INRIA Sophia-Antipolis (France), +// Max-Planck-Institute Saarbruecken (Germany), +// and Tel-Aviv University (Israel). All rights reserved. +// +// This file is part of CGAL (www.cgal.org) +// +// $URL$ +// $Id$ +// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial +// +// Author(s): Shepard Liu + #ifndef CGAL_DRAW_AOS_TYPE_UTILS_H #define CGAL_DRAW_AOS_TYPE_UTILS_H +#include #include #include +#include #include namespace CGAL { namespace draw_aos { -enum class Side_of_boundary { +enum class Boundary_side { Top = 0, Left = 1, Bottom = 2, @@ -16,72 +33,95 @@ enum class Side_of_boundary { None = -1, }; -template > -struct has_construct_x_monotone_curve_2 : std::false_type -{}; - -template -struct has_construct_x_monotone_curve_2> : std::true_type -{}; - template > struct has_approximate_2_object : std::false_type {}; -// Specialization: detection succeeds if decltype(T::approximate_2_object()) is valid -template -struct has_approximate_2_object().approximate_2_object())>> : std::true_type +template +struct has_approximate_2_object().approximate_2_object())>> : std::true_type {}; -// Convenience variable -template -inline constexpr bool has_approximate_2_object_v = has_approximate_2_object::value; +// Detect whether Gt has defined a member function approximate_2_object() +template +inline constexpr bool has_approximate_2_object_v = has_approximate_2_object::value; -// Primary templates: detection fails by default -// Does a class have operator()(const Point&)? template > -struct has_operator_point : std::false_type +struct has_approximate_point : std::false_type {}; -// Specialization: detection succeeds if decltype works out -template -struct has_operator_point()(std::declval()))>> +template +struct has_approximate_point()(std::declval()))>> : std::true_type {}; -// Convenience variable -template -inline constexpr bool has_operator_point_v = has_operator_point::value; +// Detect whether A has operator()(const Gt::Point_2&) +template +inline constexpr bool has_approximate_point_v = has_approximate_point::value; -// Primary templates: detection fails by default -// Does a class have operator()(const X_monotone_curve&)? template > -struct has_operator_xcv : std::false_type +struct has_approximate_xcv : std::false_type {}; -/*! - */ -template -struct has_operator_xcv()(std::declval(), - std::declval(), - std::declval(), - std::declval()))>> : std::true_type +template +struct has_approximate_xcv< + Gt, + A, + O, + std::void_t()(std::declval(), + std::declval(), + std::declval(), + std::declval()))>> : std::true_type {}; -// Convenience variable -template -constexpr bool has_operator_xcv_v = has_operator_xcv::value; +// Detect whether A has operator()(const Gt::X_monotone_curve_2&, double, OutputIterator, bool)? +template +constexpr bool has_approximate_xcv_v = has_approximate_xcv::value; -template +template > +struct has_approximate_xcv_with_bounds : std::false_type +{}; + +template +struct has_approximate_xcv_with_bounds< + Gt, + A, + O, + std::void_t()(std::declval(), + std::declval(), + std::declval(), + std::declval(), + std::declval()))>> : std::true_type +{}; + +// Detect whether A has operator()(const X_monotone_curve&, double, OutputIterator, Bbox_2, bool) +template +inline constexpr bool has_approximate_xcv_with_bounds_v = has_approximate_xcv_with_bounds::value; + +// Detect whether a geometry traits has all the necessary types and functions for approximation +template constexpr bool has_approximate_traits_v = - has_approximate_2_object_v && has_operator_point_v && - has_operator_xcv_v; + has_approximate_2_object_v && has_approximate_point_v && + (has_approximate_xcv_v || + has_approximate_xcv_with_bounds_v); -// Detect whether T is or derives from Arr_geodesic_arc_on_sphere_traits_2<*, *, *> -template +template > +struct has_is_in_x_range : std::false_type +{}; + +template +struct has_is_in_x_range().is_in_x_range( + std::declval()))>> : std::true_type +{}; + +// Detect whether Gt::X_monotone_curve_2 has a member function bool is_in_x_range(const Gt::Point_2&) +template +inline constexpr bool has_is_in_x_range_v = has_is_in_x_range::value; + +// Detect whether Gt is or derives from Arr_geodesic_arc_on_sphere_traits_2<*, *, *> +template struct is_or_derived_from_agas { private: @@ -91,38 +131,97 @@ private: static std::false_type test(...); public: - static constexpr bool value = decltype(test(static_cast(nullptr)))::value; + static constexpr bool value = decltype(test(static_cast(nullptr)))::value; }; -template -inline constexpr bool is_or_derived_from_agas_v = is_or_derived_from_agas::value; +template +inline constexpr bool is_or_derived_from_agas_v = is_or_derived_from_agas::value; // Detect whether T is or derives from a geometry traits on curved surfaces -template -inline constexpr bool is_or_derived_from_curved_surf_traits = is_or_derived_from_agas_v; +template +inline constexpr bool is_or_derived_from_curved_surf_traits_v = is_or_derived_from_agas_v; + +// Static helpers to get template arguments from a geometry traits +template +struct tmpl_args +{}; + +template +struct tmpl_args> +{ + using Kernel = Kernel_; + static constexpr int atan_x = AtanX; + static constexpr int atan_y = AtanY; +}; /*! + * \brief Approximation data types + * + * \tparam Gt Geometry traits */ -template +template class Arr_approximate_traits { - using Geom_traits = GeomTraits; + using Geom_traits = Gt; + + template + struct Triangle_soup_ + { + using Index = I; + using Triangle = std::array; + using Point = P; + + std::vector points; + std::vector triangles; + }; public: using Approx_point = typename Geom_traits::Approximate_point_2; using Approx_nt = typename Geom_traits::Approximate_number_type; using Approx_kernel = typename Geom_traits::Approximate_kernel; - using Approx_proj_point = typename Approx_kernel::Point_2; - using Point_geom = Approx_proj_point; - using Point_geom_vec = std::vector; - using Polyline_geom = Point_geom_vec; - using Triangle = std::array; - using Triangle_vec = std::vector; - struct Triangulated_face - { - Point_geom_vec points; - Triangle_vec triangles; - }; + + // 2D parameter space point + using Point = typename Approx_kernel::Point_2; + using Polyline = std::vector; + + using Triangle_soup = Triangle_soup_; + + // A null point with NaN coordinates. Use ::is_null(pt) to check if a point is null. + inline static const Point Null_point = + Point(std::numeric_limits::signaling_NaN(), std::numeric_limits::signaling_NaN()); + static bool is_null(Point pt) { return std::isnan(pt.x()) || std::isnan(pt.y()); } +}; + +/*! + * \brief Functor to construct a Point_2 from an Approximate_point_2. + * + * \tparam Gt Geometry traits + */ +template +class Construct_gt_point_2 +{ + using Approx_traits = Arr_approximate_traits; + using Approx_point = typename Approx_traits::Approx_point; + using Gt_point = typename Gt::Point_2; + +public: + Gt_point operator()(const Approx_point& pt) const { return Gt_point(pt.x(), pt.y()); } +}; + +// Specialization for Arr_geodesic_arc_on_sphere_traits_2 +template +class Construct_gt_point_2> +{ + using Geom_traits = Arr_geodesic_arc_on_sphere_traits_2; + using Approx_traits = Arr_approximate_traits>; + using Approx_point = typename Approx_traits::Approx_point; + using Gt_point = typename Geom_traits::Point_2; + +public: + Gt_point operator()(const Approx_point& pt) const { + using Direction_3 = typename Kernel::Direction_3; + return Gt_point(Direction_3(pt.dx(), pt.dy(), pt.dz()), static_cast(pt.location())); + } }; } // namespace draw_aos diff --git a/Arrangement_on_surface_2/include/CGAL/draw_arrangement_2.h b/Arrangement_on_surface_2/include/CGAL/draw_arrangement_2.h index 9daa24b371b..98c914cbc55 100644 --- a/Arrangement_on_surface_2/include/CGAL/draw_arrangement_2.h +++ b/Arrangement_on_surface_2/include/CGAL/draw_arrangement_2.h @@ -11,16 +11,20 @@ // $Id$ // SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial // -// Author(s): Efi Fogel +// Author(s): Efi Fogel +// Shepard Liu #ifndef CGAL_DRAW_ARRANGEMENT_2_H #define CGAL_DRAW_ARRANGEMENT_2_H #include +#include #include +#include #include #include +#include #include #include @@ -30,99 +34,16 @@ #include #include #include +#include +#include +#include #include +#include "CGAL/Bbox_2.h" namespace CGAL { -namespace draw_function_for_arrangement_2 { +namespace draw_aos { -// ============================ -// Detection idiom using void_t -// ============================ - -// ======== - -// Primary templates: detection fails by default -// Does the traits have approximate_2_object()? -template > -struct has_approximate_2_object : std::false_type -{}; - -// Specialization: detection succeeds if decltype(T::approximate_2_object()) is valid -template -struct has_approximate_2_object().approximate_2_object())>> : std::true_type -{}; - -// Convenience variable -template -inline constexpr bool has_approximate_2_object_v = has_approximate_2_object::value; - -// ======== - -// Primary templates: detection fails by default -// Does a class have operator()(const Point&)? -template > -struct has_operator_point : std::false_type -{}; - -// Specialization: detection succeeds if decltype works out -template -struct has_operator_point()(std::declval()))>> - : std::true_type -{}; - -// Convenience variable -template -inline constexpr bool has_operator_point_v = has_operator_point::value; - -// ======== - -// Primary templates: detection fails by default -// Does a class have operator()(const X_monotone_curve&)? -template > -struct has_operator_xcv : std::false_type -{}; - -// Specialization: detection succeeds if decltype works out -struct Dummy_output -{}; -using Concrete_output_iterator = Dummy_output*; - -template -struct has_operator_xcv()(std::declval(), - std::declval, - std::declval(), - std::declval()))>> : std::true_type -{}; - -// Convenience variable -template -inline constexpr bool has_operator_xcv_v = has_operator_xcv::value; - -// ======== - -// Helper: detect whether T is or derives from Arr_geodesic_arc_on_sphere_traits_2<*, *, *> -template -struct is_or_derived_from_agas -{ -private: - template - static std::true_type test(const Arr_geodesic_arc_on_sphere_traits_2*); - - static std::false_type test(...); - -public: - static constexpr bool value = decltype(test(static_cast(nullptr)))::value; -}; - -template -inline constexpr bool is_or_derived_from_agas_v = is_or_derived_from_agas::value; - -// ======== - -/// template class Draw_arr_tool { @@ -149,7 +70,7 @@ public: // std::cout << "add_face()\n"; for(Inner_ccb_const_iterator it = face->inner_ccbs_begin(); it != face->inner_ccbs_end(); ++it) add_ccb(*it); - if (! face->is_unbounded()) { + if(!face->is_unbounded()) { for(Outer_ccb_const_iterator it = face->outer_ccbs_begin(); it != face->outer_ccbs_end(); ++it) { add_ccb(*it); draw_region(*it); @@ -210,13 +131,13 @@ public: /// Compile time dispatching /// - template , int> = 0> + template , int> = 0> void draw_region_impl2(const T& /* traits */, const A& /* approximate */, Halfedge_const_handle curr) { draw_exact_region(curr); } /// - template , int> = 0> + template , int> = 0> auto draw_region_impl2(const T& /* traits */, const A& approx, Halfedge_const_handle curr) { draw_approximate_region(curr, approx); } @@ -308,14 +229,14 @@ public: } /// - template , int> = 0> + template , int> = 0> void draw_point_impl2( const T& /* traits */, const A& /* approximate */, const Point& p, bool colored, const CGAL::IO::Color& c) { draw_exact_point(p, colored, c); } /// - template , int> = 0> + template , int> = 0> auto draw_point_impl2(const T& /* traits */, const A& approx, const Point& p, bool colored, const CGAL::IO::Color& c) { draw_approximate_point(p, approx, colored, c); @@ -475,7 +396,7 @@ public: } /// - template , int> = 0> + template , int> = 0> void draw_curve_impl2(const T& /* traits */, const A& /* approximate */, const X_monotone_curve& xcv, @@ -485,7 +406,7 @@ public: } /// - template , int> = 0> + template , int> = 0> auto draw_curve_impl2( const T& /* traits */, const A& approx, const X_monotone_curve& xcv, bool colored, const CGAL::IO::Color& c) { draw_approximate_curve(xcv, approx, colored, c); @@ -580,29 +501,273 @@ protected: std::unordered_map m_visited; }; -} // namespace draw_function_for_arrangement_2 +template +static auto map_from_pair_ranges(Range1 range1, Range2 range2) { + CGAL_assertion_msg(range1.size() == range2.size(), "The two ranges must have the same size."); + auto begin = boost::make_zip_iterator(boost::make_tuple(range1.begin(), range2.begin())); + auto end = boost::make_zip_iterator(boost::make_tuple(range1.end(), range2.end())); + auto tuple_to_pair = [](const auto& t) { return std::make_pair(boost::get<0>(t), boost::get<1>(t)); }; + return unordered_flat_map(boost::make_transform_iterator(begin, tuple_to_pair), + boost::make_transform_iterator(end, tuple_to_pair)); +} -#define CGAL_ARR_TYPE CGAL::Arrangement_on_surface_2 +/*! + * \brief tracking changes between an arrangement and its copy that will be later inserted to. + * + * \note tracks insertions only. If any other actions made(e.g. deletions, merging, etc), the state of the tracker + * instance may become invalid. + * + * \tparam Arrangement + */ +template +class Arr_insertion_tracker : Arr_observer +{ + using Base = Arr_observer; + using Halfedge_handle = typename Arrangement::Halfedge_handle; + using Halfedge_const_handle = typename Arrangement::Halfedge_const_handle; + using Face_handle = typename Arrangement::Face_handle; + using Face_const_handle = typename Arrangement::Face_const_handle; + using Vertex_handle = typename Arrangement::Vertex_handle; + using Vertex_const_handle = typename Arrangement::Vertex_const_handle; + using X_monotone_curve_2 = typename Arrangement::X_monotone_curve_2; -/// Draw an arrangement on surface. -template -void draw(const CGAL_ARR_TYPE& aos, - const GSOptions& gso, - const char* title = "2D Arrangement on Surface Basic Viewer") { - using Arrangement = CGAL_ARR_TYPE; +protected: + virtual void after_create_vertex(Vertex_handle v) override { m_vertex_map[v] = Vertex_const_handle(); } - Qt::init_ogl_context(4, 3); - int argc; - QApplication app(argc, nullptr); - auto viewer = draw_aos::Arr_viewer(app.activeWindow(), aos, gso, title); + virtual void after_create_edge(Halfedge_handle e) override { + m_halfedge_map[e] = Halfedge_const_handle(); + m_halfedge_map[e->twin()] = Halfedge_const_handle(); // twin is created as well + } + + virtual void before_split_edge(Halfedge_handle e, + Vertex_handle v, + const X_monotone_curve_2& c1, + const X_monotone_curve_2& c2) override { + if(m_vertex_map.find(v) == m_vertex_map.end()) m_vertex_map[v] = Vertex_const_handle(); // v is newly created + } + + virtual void after_split_edge(Halfedge_handle e1, Halfedge_handle e2) override { + if(auto it = m_halfedge_map.find(e1); it == m_halfedge_map.end()) + m_halfedge_map[e2] = e1; + else if(it->second == Halfedge_const_handle()) + m_halfedge_map[e2] = Halfedge_const_handle(); // e1 has no corresponding edge in the original arrangement + else + m_halfedge_map[e2] = it->second; // e1 is created by splitting an existing edge + } + + virtual void after_split_face(Face_handle f1, Face_handle f2, bool) override { + // Face cannot be created but by splitting an existing face. + if(auto it = m_face_map.find(f1); it == m_face_map.end()) + m_face_map[f2] = f1; + else + m_face_map[f2] = it->second; // f1 is created by splitting an existing face + } + +public: + Arr_insertion_tracker(Arrangement& arr) + : Base(arr) {} + + /*! + * \brief Query the original face of a given face. + * + * \param fh a valid face handle in the modified arrangement. + * \return Face_const_handle + */ + Face_const_handle original_face(Face_const_handle fh) const { + auto it = m_face_map.find(fh); + if(it == m_face_map.end()) return fh; + return it->second; // new face from splitting an existing face + } + + /*! + * \brief Query the original halfedge of a given halfedge. + * + * \param heh a valid halfedge handle in the modified arrangement. + * \return Halfedge_const_handle + */ + Halfedge_const_handle original_halfedge(Halfedge_const_handle he) const { + auto it = m_halfedge_map.find(he); + if(it == m_halfedge_map.end()) return he; + if(it->second == Halfedge_const_handle()) return Halfedge_const_handle(); // newly created halfedge + return it->second; + } + + /*! + * \brief Query the original vertex of a given vertex. + * + * \param vh a valid vertex handle in the modified arrangement. + * \return Vertex_const_handle + */ + Vertex_const_handle original_vertex(Vertex_const_handle vh) const { + auto it = m_vertex_map.find(vh); + if(it == m_vertex_map.end()) return vh; + if(it->second == Vertex_const_handle()) return Vertex_const_handle(); // newly created vertex + return it->second; // it will never reach here. + } + +private: + /*! + * Maps tracking the changes between the original arrangement and modified arrangement. + * The key is the current feature, and the value is the corresponding feature before modification. + * If there is no entry about a feature, the corresponding feature is itself. + * If the value is a invalid handle, it means that the feature is newly created and thus has no corresponding + * feature in the original arrangement. + */ + unordered_flat_map m_face_map; + unordered_flat_map m_halfedge_map; + unordered_flat_map m_vertex_map; +}; + +void draw_unimplemented() { + std::cerr << "Geometry traits type of arrangement is required to support approximation of Point_2 and " + "X_monotone_curve_2. Traits on curved surfaces needs additional support for parameterization." + << std::endl; + exit(1); +} + +template +void draw_impl_planar( + const Arrangement& arr, const GSOptions& gso, const char* title, Bbox_2 initial_bbox, QApplication& app) { + Arr_viewer viewer(app.activeWindow(), arr, gso, title, initial_bbox); viewer.show(); app.exec(); } -/// Draw an arrangement on surface. -template -void draw(const CGAL_ARR_TYPE& aos, const char* title = "2D Arrangement on Surface Basic Viewer") { - using Arrangement = CGAL_ARR_TYPE; +template +void draw_impl_agas( + const Arrangement& arr, const GSOptions& gso, const char* title, Bbox_2 initial_bbox, QApplication& app) { + using Halfedge_const_handle = typename Arrangement::Halfedge_const_handle; + using Face_const_handle = typename Arrangement::Face_const_handle; + using Vertex_const_handle = typename Arrangement::Vertex_const_handle; + using Geom_traits = typename Arrangement::Geometry_traits_2; + using X_monotone_curve_2 = typename Geom_traits::X_monotone_curve_2; + using Direction_3 = typename Geom_traits::Direction_3; + using Point_2 = typename Geom_traits::Point_2; + using Agas_template_args = tmpl_args; + + Arrangement derived_arr(arr); + auto vertex_map = map_from_pair_ranges(derived_arr.vertex_handles(), + arr.vertex_handles()); + auto halfedge_map = map_from_pair_ranges(derived_arr.halfedge_handles(), + arr.halfedge_handles()); + auto face_map = + map_from_pair_ranges(derived_arr.face_handles(), arr.face_handles()); + // setup tracker and insert the identification curve. + Arr_insertion_tracker tracker(derived_arr); + X_monotone_curve_2 id_curve = arr.geometry_traits()->construct_x_monotone_curve_2_object()( + Point_2(Direction_3(0, 0, -1), Point_2::MIN_BOUNDARY_LOC), + Point_2(Direction_3(0, 0, 1), Point_2::MAX_BOUNDARY_LOC), + Direction_3(Agas_template_args::atan_y, -Agas_template_args::atan_x, 0)); + insert(derived_arr, id_curve); + + // derived_gso proxies the call to the original gso + GSOptions derived_gso(gso); + derived_gso.draw_vertex = [&](const Arrangement&, const Vertex_const_handle& vh) { + Vertex_const_handle original_vh = tracker.original_vertex(vh); + if(original_vh == Vertex_const_handle() || vertex_map.find(original_vh) == vertex_map.end()) return false; + return gso.draw_vertex(arr, vertex_map.at(original_vh)); + }; + derived_gso.colored_vertex = [&](const Arrangement&, const Vertex_const_handle& vh) { + Vertex_const_handle original_vh = tracker.original_vertex(vh); + if(original_vh == Vertex_const_handle() || vertex_map.find(original_vh) == vertex_map.end()) return false; + return gso.colored_vertex(arr, vertex_map.at(original_vh)); + }; + derived_gso.vertex_color = [&](const Arrangement&, const Vertex_const_handle& vh) -> CGAL::IO::Color { + Vertex_const_handle original_vh = tracker.original_vertex(vh); + if(original_vh == Vertex_const_handle() || vertex_map.find(original_vh) == vertex_map.end()) + return CGAL::IO::Color(); + return gso.vertex_color(arr, vertex_map.at(original_vh)); + }; + derived_gso.draw_edge = [&](const Arrangement&, const Halfedge_const_handle& he) { + Halfedge_const_handle original_he = tracker.original_halfedge(he); + if(original_he == Halfedge_const_handle() || halfedge_map.find(original_he) == halfedge_map.end()) return false; + return gso.draw_edge(arr, halfedge_map.at(original_he)); + }; + derived_gso.colored_edge = [&](const Arrangement&, const Halfedge_const_handle& he) { + Halfedge_const_handle original_he = tracker.original_halfedge(he); + if(original_he == Halfedge_const_handle() || halfedge_map.find(original_he) == halfedge_map.end()) return false; + return gso.colored_edge(arr, halfedge_map.at(original_he)); + }; + derived_gso.edge_color = [&](const Arrangement&, const Halfedge_const_handle& he) -> CGAL::IO::Color { + Halfedge_const_handle original_he = tracker.original_halfedge(he); + if(original_he == Halfedge_const_handle() || halfedge_map.find(original_he) == halfedge_map.end()) + return CGAL::IO::Color(); + return gso.edge_color(arr, halfedge_map.at(original_he)); + }; + derived_gso.draw_face = [&](const Arrangement&, const Face_const_handle& fh) { + Face_const_handle original_fh = tracker.original_face(fh); + if(face_map.find(original_fh) == face_map.end()) return false; + return gso.draw_face(arr, face_map.at(original_fh)); + }; + derived_gso.colored_face = [&](const Arrangement&, const Face_const_handle& fh) { + Face_const_handle original_fh = tracker.original_face(fh); + if(face_map.find(original_fh) == face_map.end()) return false; + return gso.draw_face(arr, face_map.at(original_fh)); + }; + derived_gso.face_color = [&](const Arrangement&, const Face_const_handle& fh) -> CGAL::IO::Color { + Face_const_handle original_fh = tracker.original_face(fh); + if(face_map.find(original_fh) == face_map.end()) return CGAL::IO::Color(); + return gso.face_color(arr, face_map.at(original_fh)); + }; + + Arr_viewer viewer(app.activeWindow(), derived_arr, derived_gso, title, initial_bbox); + viewer.show(); + app.exec(); +} + +template +void draw(const Arrangement& arr, const GSOptions& gso, Args&&... args) { + using Geom_traits = typename Arrangement::Geometry_traits_2; + + if constexpr(!has_approximate_traits_v) + return draw_unimplemented(); + else if constexpr(is_or_derived_from_agas_v) + // Arrangements on curved surfaces require special handling. The identification curve must be present to make the + // curved surface homeomorphic to a bounded plane. + return draw_impl_agas(arr, gso, std::forward(args)...); + else + return draw_impl_planar(arr, gso, std::forward(args)...); +} + +} // namespace draw_aos + +/*! + * \brief Draw an arrangement on surface. + * + * \tparam Arrangement + * \tparam GSOptions + * \param arr the arrangement to be drawn + * \param gso graphics scene options + * \param title title of the viewer window + * \param initial_bbox parameter space bounding box to be shown intially. If empty, the approximate bounding box of the + * arrangement is used. For arrangements induced by unbounded curves, the default initial bounding box is computed from + * vertex coordinates. + */ +template +void draw(const Arrangement& arr, + const GSOptions& gso, + const char* title = "2D Arrangement on Surface Viewer", + Bbox_2 initial_bbox = Bbox_2(0, 0, 0, 0)) { + Qt::init_ogl_context(4, 3); + int argc; + QApplication app(argc, nullptr); + draw_aos::draw(arr, gso, title, initial_bbox, app); +} + +/*! + * \brief Draw an arrangement on surface with default graphics scene options. Faces are colored randomly. + * + * \tparam Arrangement + * \param arr the arrangement to be drawn + * \param title title of the viewer window + * \param initial_bbox parameter space bounding box to be shown intially. If empty, the approximate bounding box of the + * arrangement is used. For arrangements induced by unbounded curves, the default initial bounding box is computed from + * vertex coordinates. + */ +template +void draw(const Arrangement& arr, + const char* title = "2D Arrangement on Surface Viewer", + Bbox_2 initial_bbox = Bbox_2(0, 0, 0, 0)) { using Face_const_handle = typename Arrangement::Face_const_handle; using Vertex_const_handle = typename Arrangement::Vertex_const_handle; using Halfedge_const_handle = typename Arrangement::Halfedge_const_handle; @@ -610,32 +775,73 @@ void draw(const CGAL_ARR_TYPE& aos, const char* title = "2D Arrangement on Surfa CGAL::Graphics_scene_options; GSOptions gso; - gso.enable_faces(); - gso.colored_face = [](const Arrangement&, const Face_const_handle&) { return true; }; - gso.face_color = [](const Arrangement&, const Face_const_handle& fh) -> CGAL::IO::Color { - CGAL::Random random((size_t(fh.ptr()))); - return get_random_color(random); - }; - gso.enable_edges(); - gso.colored_edge = [](const Arrangement&, const Halfedge_const_handle&) { return true; }; - gso.edge_color = [](const Arrangement&, const Halfedge_const_handle& heh) -> CGAL::IO::Color { - return CGAL::IO::Color(0, 0, 0); - }; gso.enable_vertices(); + gso.draw_vertex = [](const Arrangement&, const Vertex_const_handle&) { return true; }; gso.colored_vertex = [](const Arrangement&, const Vertex_const_handle&) { return true; }; gso.vertex_color = [](const Arrangement&, const Vertex_const_handle& vh) -> CGAL::IO::Color { return CGAL::IO::Color(255, 0, 0); }; + gso.enable_edges(); + gso.draw_edge = [](const Arrangement&, const Halfedge_const_handle&) { return true; }; + gso.colored_edge = [](const Arrangement&, const Halfedge_const_handle&) { return true; }; + gso.edge_color = [](const Arrangement&, const Halfedge_const_handle& heh) -> CGAL::IO::Color { + return CGAL::IO::Color(0, 0, 0); + }; + gso.enable_faces(); + gso.draw_face = [](const Arrangement&, const Face_const_handle&) { return true; }; + gso.colored_face = [](const Arrangement&, const Face_const_handle&) { return true; }; + gso.face_color = [](const Arrangement&, const Face_const_handle& fh) -> CGAL::IO::Color { + CGAL::Random random(std::size_t(fh.ptr())); + return get_random_color(random); + }; - Qt::init_ogl_context(4, 3); - int argc; - QApplication app(argc, nullptr); - auto viewer = draw_aos::Arr_viewer(app.activeWindow(), aos, gso, title); - viewer.show(); - app.exec(); + draw(arr, gso, title, initial_bbox); } -#undef CGAL_ARR_TYPE +#define CGAL_ARR_TYPE CGAL::Arrangement_on_surface_2 + +/// +template +void add_to_graphics_scene(const CGAL_ARR_TYPE& aos, CGAL::Graphics_scene& graphics_scene, const GSOptions& gso) { + draw_aos::Draw_arr_tool dar(aos, graphics_scene, gso); + dar.add_elements(); +} + +/// +template +void add_to_graphics_scene(const CGAL_ARR_TYPE& aos, CGAL::Graphics_scene& graphics_scene) { + CGAL::Graphics_scene_options + gso; + // colored face? + gso.colored_face = [](const CGAL_ARR_TYPE&, typename CGAL_ARR_TYPE::Face_const_handle) -> bool { return true; }; + + // face color + gso.face_color = [](const CGAL_ARR_TYPE&, typename CGAL_ARR_TYPE::Face_const_handle fh) -> CGAL::IO::Color { + CGAL::Random random((unsigned int)(std::size_t)(&*fh)); + return get_random_color(random); + }; + + add_to_graphics_scene(aos, graphics_scene, gso); +} + +/// Draw an arrangement on surface. +template +void draw_old(const CGAL_ARR_TYPE& aos, + const GSOptions& gso, + const char* title = "2D Arrangement on Surface Basic Viewer") { + CGAL::Graphics_scene graphics_scene; + add_to_graphics_scene(aos, graphics_scene, gso); + draw_graphics_scene(graphics_scene, title); +} + +/// Draw an arrangement on surface. +template +void draw_old(const CGAL_ARR_TYPE& aos, const char* title = "2D Arrangement on Surface Basic Viewer") { + CGAL::Graphics_scene graphics_scene; + add_to_graphics_scene(aos, graphics_scene); + draw_graphics_scene(graphics_scene, title); +} } // namespace CGAL