feat: draw function and classes for arrangement on surface 2

Now works for segment traits, linear traits and conic traits.

.clang-format

@@ -6,6 +6,7 @@ AllowShortFunctionsOnASingleLine: true
 BinPackParameters: false
 BreakConstructorInitializers: BeforeComma
 BreakBeforeBraces: Custom
+AlwaysBreakTemplateDeclarations: Yes
 BraceWrapping:
   AfterCaseLabel:  false
   AfterClass:      true
@@ -32,4 +33,6 @@ PointerAlignment: Left
 # Control the spaces around conditionals
 SpacesInConditionalStatement: false
 SpaceBeforeParens: false
+# disable sorting of includes
+SortIncludes: false
 ...

Arrangement_on_surface_2/demo/Arrangement_on_surface_2/ArrangementDemoGraphicsView.cpp

@@ -12,44 +12,39 @@
 
 #include "ArrangementDemoGraphicsView.h"
 
-#include <iostream>
-#include <cmath>
-#include <QVarLengthArray>
-#include <QPen>
 #include <QCoreApplication>
-#include <QKeyEvent>
 #include <QFontMetrics>
+#include <QKeyEvent>
+#include <QPen>
+#include <QVarLengthArray>
+#include <cmath>
+#include <iostream>
 
 //! Member function to setup the viewport of the screen
 /*!
   \param parent a Qwidget pointer to the class
 */
-ArrangementDemoGraphicsView::ArrangementDemoGraphicsView( QWidget* parent ) :
+ArrangementDemoGraphicsView::ArrangementDemoGraphicsView(QWidget* parent)
-  QGraphicsView( parent ),
+    : QGraphicsView(parent)
-  maxScale( 500000 ),
+    , maxScale(500000)
-  minScale( 0.0002 )
+    , minScale(0.0002) {
-{
   this->resetTransform();
   this->setResizeAnchor(QGraphicsView::AnchorUnderMouse);
   this->setViewportUpdateMode(QGraphicsView::FullViewportUpdate);
-  this->setMouseTracking( true );
+  this->setMouseTracking(true);
   this->setHorizontalScrollBarPolicy(Qt::ScrollBarAlwaysOff);
   this->setVerticalScrollBarPolicy(Qt::ScrollBarAlwaysOff);
   // TODO: Make options menu work
   this->setRenderHint(QPainter::Antialiasing);
 }
 
-void ArrangementDemoGraphicsView::paintEvent(QPaintEvent* event)
+void ArrangementDemoGraphicsView::paintEvent(QPaintEvent* event) {
-{
   qreal scale = std::sqrt(std::abs(this->transform().determinant()));
-  if (scale > this->maxScale)
+  if(scale > this->maxScale)
     this->scale(this->maxScale / scale, this->maxScale / scale);
-  else if (scale < this->minScale)
+  else if(scale < this->minScale)
     this->scale(this->minScale / scale, this->minScale / scale);
   QGraphicsView::paintEvent(event);
 }
 
-void ArrangementDemoGraphicsView::resetTransform()
+void ArrangementDemoGraphicsView::resetTransform() { this->setTransform({1.0, 0.0, 0.0, -1.0, 0.0, 0.0}); }
-{
-  this->setTransform({1.0, 0.0, 0.0, -1.0, 0.0, 0.0});
-}

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_approximate_point_2.h

@@ -0,0 +1,97 @@
+#ifndef CGAL_DRAW_AOS_ARR_APPROXIMATE_POINT_2_H
+#define CGAL_DRAW_AOS_ARR_APPROXIMATE_POINT_2_H
+
+#include "CGAL/Arr_has.h"
+#include "CGAL/Draw_aos/Arr_approximation_geometry_traits.h"
+#include "CGAL/number_utils.h"
+
+namespace CGAL {
+
+namespace internal {
+
+template <typename Geom_traits, bool Has_approximate_2>
+class Arr_approximate_point_2_impl;
+
+template <typename Geom_traits>
+class Arr_approximate_point_2_impl<Geom_traits, true>
+{
+  using Approx_kernel = Arr_approximation_geometry_traits::Approximation_kernel;
+  using Point_2 = typename Geom_traits::Point_2;
+  using Approx_point = Arr_approximation_geometry_traits::Approx_point;
+
+public:
+  Arr_approximate_point_2_impl(const Geom_traits& traits)
+      : m_approx(traits.approximate_2_object()) {}
+
+  /**
+   * @brief Approximate a point.
+   * TODO: make it work for spherical traits.
+   *
+   * @param pt
+   * @return Point_geom
+   */
+  Approx_point operator()(const Point_2& pt) const { return {m_approx(pt, 0), m_approx(pt, 1)}; }
+
+  /**
+   * @brief Approximate a specific dimension of a point.
+   *
+   * @param pt
+   * @param dim 0 for x, 1 for y, 2 for z if applicable. An exception is thrown if the dimension is invalid.
+   * @return double
+   */
+  double operator()(const Point_2& pt, int dim) const { return m_approx(pt, dim); }
+
+private:
+  typename Geom_traits::Approximate_2 m_approx;
+};
+
+// Fallback to use CGAL::to_double for traits that do not have an approximate_2_object.
+template <typename Geom_traits>
+class Arr_approximate_point_2_impl<Geom_traits, false>
+{
+  using Approx_kernel = Arr_approximation_geometry_traits::Approximation_kernel;
+  using Point_2 = typename Geom_traits::Point_2;
+  using Approx_point = Arr_approximation_geometry_traits::Approx_point;
+
+public:
+  // traits object  is not used in the fallback implementation, but we keep it for consistency.
+  Arr_approximate_point_2_impl(const Geom_traits& traits) {}
+
+  /**
+   * @brief Approximate a point.
+   *
+   * @note this functions does not check if the point is within the bounding box.
+   * @param pt
+   * @return Point_geom
+   */
+  Approx_point operator()(const Point_2& pt) const { return {CGAL::to_double(pt.x()), CGAL::to_double(pt.y())}; }
+
+  /**
+   * @brief Approximate a specific dimension of a point.
+   *
+   * @param pt
+   * @param dim 0 for x, 1 for y, 2 for z if applicable. An exception is thrown if the dimension is invalid.
+   * @return double
+   */
+  double operator()(const Point_2& pt, int dim) const {
+    if(dim == 0) {
+      return CGAL::to_double(pt.x());
+    } else if(dim == 1) {
+      return CGAL::to_double(pt.y());
+    } else if(dim == 2 && pt.dimension() == 3) {
+      return CGAL::to_double(pt.z());
+    } else {
+      throw std::invalid_argument("Invalid dimension for approximate_point");
+    }
+  }
+};
+
+} // namespace internal
+
+template <typename Geom_traits>
+using Arr_approximate_point_2 =
+    internal::Arr_approximate_point_2_impl<Geom_traits, has_approximate_2<Geom_traits>::value>;
+
+} // namespace CGAL
+
+#endif // CGAL_DRAW_AOS_ARR_APPROXIMATE_POINT_2_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_approximation_cache.h

@@ -0,0 +1,105 @@
+#ifndef CGAL_DRAW_AOS_ARR_APPROXIMATION_CACHE_H
+#define CGAL_DRAW_AOS_ARR_APPROXIMATION_CACHE_H
+#include "CGAL/Arr_enums.h"
+#include "CGAL/Draw_aos/Arr_approximation_geometry_traits.h"
+#include "CGAL/Draw_aos/helpers.h"
+#include "CGAL/unordered_flat_map.h"
+#include <boost/range/iterator_range.hpp>
+
+namespace CGAL {
+class Arr_approximation_cache
+{
+  using Approx_geom_traits = Arr_approximation_geometry_traits;
+
+public:
+  using Vertex_cache_obj = Approx_geom_traits::Point_geom;
+  using Halfedge_cache_obj = Approx_geom_traits::Polyline_geom;
+  using Face_cache_obj = Approx_geom_traits::Triangulated_face;
+
+private:
+  using Vertex_const_handle = Arrangement::Vertex_const_handle;
+  using Edge_const_handle = Arrangement::Edge_const_iterator;
+  using Halfedge_const_handle = Arrangement::Halfedge_const_iterator;
+  using Face_const_handle = Arrangement::Face_const_handle;
+  using Vertex_cache = unordered_flat_map<Vertex_const_handle, Vertex_cache_obj>;
+  using Halfedge_cache = unordered_flat_map<Halfedge_const_handle, Halfedge_cache_obj>;
+  using Face_cache = unordered_flat_map<Face_const_handle, Face_cache_obj>;
+
+  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<Vertex_cache_const_iterator>;
+  using Halfedge_cache_range = boost::iterator_range<Halfedge_cache_const_iterator>;
+  using Face_cache_range = boost::iterator_range<Face_cache_const_iterator>;
+
+private:
+  Halfedge_const_handle identify_halfedge(const Halfedge_const_handle& he) const {
+    return he->direction() == ARR_RIGHT_TO_LEFT ? he : he->twin();
+  }
+
+public:
+  std::pair<Vertex_cache_obj&, bool> try_emplace(const Vertex_const_handle& vh) {
+    const auto& [it, inserted] = m_vertex_cache.try_emplace(vh, Vertex_cache_obj());
+    return {it->second, inserted};
+  }
+  std::pair<Halfedge_cache_obj&, bool> try_emplace(const Halfedge_const_handle& he) {
+    const auto& [it, inserted] = m_halfedge_cache.try_emplace(identify_halfedge(he), Halfedge_cache_obj());
+    return {it->second, inserted};
+  }
+  std::pair<Face_cache_obj&, bool> try_emplace(const Face_const_handle& fh) {
+    const auto& [it, inserted] = m_face_cache.try_emplace(fh, Face_cache_obj());
+    return {it->second, inserted};
+  }
+
+  std::pair<const Vertex_cache_obj&, bool> get(const Vertex_const_handle& vh) const {
+    auto it = m_vertex_cache.find(vh);
+    if(it != m_vertex_cache.end()) {
+      return {it->second, true};
+    }
+    return {Vertex_cache_obj(), false};
+  }
+  std::pair<const Halfedge_cache_obj&, bool> get(const Halfedge_const_handle& he) const {
+    auto it = m_halfedge_cache.find(identify_halfedge(he));
+    if(it != m_halfedge_cache.end()) {
+      return {it->second, true};
+    }
+    return {Halfedge_cache_obj(), false};
+  }
+  std::pair<const Face_cache_obj&, bool> get(const Face_const_handle& fh) const {
+    auto it = m_face_cache.find(fh);
+    if(it != m_face_cache.end()) {
+      return {it->second, true};
+    }
+    return {Face_cache_obj(), false};
+  }
+
+  bool has(const Vertex_const_handle& vh) const { return m_vertex_cache.find(vh) != m_vertex_cache.end(); }
+  bool has(const Halfedge_const_handle& he) const {
+    return m_halfedge_cache.find(identify_halfedge(he)) != m_halfedge_cache.end();
+  }
+  bool has(const Face_const_handle& fh) const { return m_face_cache.find(fh) != m_face_cache.end(); }
+
+  Vertex_cache_const_iterator vertex_cache_begin() const { return m_vertex_cache.begin(); }
+  Vertex_cache_const_iterator vertex_cache_end() const { return m_vertex_cache.end(); }
+  Vertex_cache_range vertex_cache() const {
+    return boost::make_iterator_range(vertex_cache_begin(), vertex_cache_end());
+  }
+
+  Halfedge_cache_const_iterator halfedge_cache_begin() const { return m_halfedge_cache.begin(); }
+  Halfedge_cache_const_iterator halfedge_cache_end() const { return m_halfedge_cache.end(); }
+  Halfedge_cache_range halfedge_cache() const {
+    return boost::make_iterator_range(halfedge_cache_begin(), halfedge_cache_end());
+  }
+
+  Face_cache_const_iterator face_cache_begin() const { return m_face_cache.begin(); }
+  Face_cache_const_iterator face_cache_end() const { return m_face_cache.end(); }
+  Face_cache_range face_cache() const { return boost::make_iterator_range(face_cache_begin(), face_cache_end()); }
+
+private:
+  Vertex_cache m_vertex_cache;
+  Halfedge_cache m_halfedge_cache;
+  Face_cache m_face_cache;
+};
+
+} // namespace CGAL
+#endif // CGAL_DRAW_AOS_ARR_APPROXIMATION_CACHE_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_approximation_geometry_traits.h

@@ -0,0 +1,26 @@
+#ifndef CGAL_DRAW_AOS_ARR_APPROXIMATION_GEOMETRY_TRAITS_H
+#define CGAL_DRAW_AOS_ARR_APPROXIMATION_GEOMETRY_TRAITS_H
+
+#include "CGAL/Simple_cartesian.h"
+
+namespace CGAL {
+class Arr_approximation_geometry_traits
+{
+public:
+  using Approximation_kernel = Simple_cartesian<double>;
+  using Approx_point = Approximation_kernel::Point_2;
+  using FT = double;
+  using Point_geom = Approx_point;
+  using Apporx_point_vec = std::vector<Point_geom>;
+  using Polyline_geom = Apporx_point_vec;
+  using Triangle = std::array<std::size_t, 3>;
+  using Triangle_vec = std::vector<Triangle>;
+  struct Triangulated_face
+  {
+    Apporx_point_vec points;
+    Triangle_vec triangles;
+  };
+};
+
+} // namespace CGAL
+#endif // CGAL_DRAW_AOS_ARR_APPROXIMATION_GEOMETRY_TRAITS_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_curve_2.h

@@ -0,0 +1,280 @@
+#ifndef CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_CURVE_2_H
+#define CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_CURVE_2_H
+
+#include "CGAL/Arr_enums.h"
+#include "CGAL/Draw_aos/Arr_bounded_approximate_point_2.h"
+#include "CGAL/Draw_aos/Arr_bounded_compute_y_at_x.h"
+#include "CGAL/Draw_aos/Arr_construct_curve_end.h"
+#include "CGAL/Draw_aos/Arr_construct_segments.h"
+#include "CGAL/Draw_aos/Arr_render_context.h"
+#include "CGAL/Draw_aos/helpers.h"
+#include "CGAL/Draw_aos/Arr_approximation_geometry_traits.h"
+#include "CGAL/basic.h"
+#include <algorithm>
+#include <cstddef>
+#include <iterator>
+#include <optional>
+#include <vector>
+
+namespace CGAL {
+
+/**
+ * @brief Functor to approximate an x-monotone curve within an bounding box.
+ * The bbox here has closed boundary.
+ *
+ * The Approximation is done from xmin to xmax with a given step. For parts outbound the y limits and precedes or
+ * succeeds a part within, the approximation may be skipped but there will be at least one point outside the bbox for
+ * indication.
+ *
+ * @note Bounded approximation is meaningful only when the curve has at least two points within the bbox (boundary
+ * points included).
+ *
+ * TODO: Possible optimizations:
+ * - Specialize for traits that models Approximate_2 on curves.
+ */
+class Arr_bounded_approximate_curve_2
+{
+  using FT = Geom_traits::FT;
+  using Point_2 = Geom_traits::Point_2;
+  using X_monotone_curve_2 = Geom_traits::X_monotone_curve_2;
+  using Halfedge_const_handle = Arrangement::Halfedge_const_iterator;
+  using Approx_point = Arr_approximation_geometry_traits::Approx_point;
+  using Point_geom = Arr_approximation_geometry_traits::Point_geom;
+  using Polyline_geom = Arr_approximation_geometry_traits::Polyline_geom;
+  using Intersections_vector = std::vector<Point_2>;
+
+  struct Execution_context : public Arr_context_delegator<Arr_bounded_render_context>
+  {
+    Execution_context(const Arr_bounded_render_context& ctx,
+                      const X_monotone_curve_2& curve,
+                      const Arr_bounded_approximate_point_2& approx_pt,
+                      const Arr_bounded_compute_y_at_x& compute_y_at_x,
+                      const Intersections_vector& top_inters,
+                      const Intersections_vector& bottom_inters,
+                      Polyline_geom& polyline)
+        : Arr_context_delegator(ctx)
+        , curve(curve)
+        , bounded_compute_y_at_x(compute_y_at_x)
+        , bounded_approx_pt(approx_pt)
+        , top_inters(top_inters)
+        , bottom_inters(bottom_inters)
+        , min_end(Arr_construct_curve_end<Geom_traits>(ctx.traits)(curve, ARR_MIN_END))
+        , max_end(Arr_construct_curve_end<Geom_traits>(ctx.traits)(curve, ARR_MAX_END))
+        , tight_xmin(is_min_end_bounded() ? std::clamp(min_end->x(), FT(ctx.xmin()), FT(ctx.xmax())) : FT(ctx.xmin()))
+        , tight_xmax(is_max_end_bounded() ? std::clamp(max_end->x(), FT(ctx.xmin()), FT(ctx.xmax())) : FT(ctx.xmax()))
+        , tight_ymin(is_min_end_bounded() ? std::clamp(min_end->y(), FT(ctx.ymin()), FT(ctx.ymax())) : FT(ctx.ymin()))
+        , tight_ymax(is_max_end_bounded() ? std::clamp(max_end->y(), FT(ctx.ymin()), FT(ctx.ymax())) : FT(ctx.ymax()))
+        , out_it(std::back_inserter(polyline)) {}
+
+    bool has_y_intersections() const { return !top_inters.empty() || !bottom_inters.empty(); }
+    bool is_min_end_bounded() const { return min_end.has_value(); }
+    bool is_max_end_bounded() const { return max_end.has_value(); }
+    bool is_bounded_curve() const { return is_min_end_bounded() && is_max_end_bounded(); }
+
+    const X_monotone_curve_2& curve;
+    const Arr_bounded_compute_y_at_x& bounded_compute_y_at_x;
+    const Arr_bounded_approximate_point_2& bounded_approx_pt;
+    const Intersections_vector &top_inters, bottom_inters;
+    const std::optional<Point_2> min_end, max_end;
+    const FT tight_xmin, tight_xmax, tight_ymin, tight_ymax;
+    std::back_insert_iterator<Polyline_geom> out_it;
+  };
+
+private:
+  static std::vector<Point_2> compute_intersections(const X_monotone_curve_2& cv1,
+                                                    const X_monotone_curve_2& cv2,
+                                                    const typename Geom_traits::Intersect_2& intersect_2,
+                                                    const Arr_construct_curve_end<Geom_traits>& cst_curve_end) {
+    using Intersect_point = std::pair<Geom_traits::Point_2, Geom_traits::Multiplicity>;
+    using Intersect_curve = Geom_traits::X_monotone_curve_2;
+    using Intersect_type = std::variant<Intersect_point, Intersect_curve>;
+
+    std::vector<Point_2> intersections;
+    auto out_it = std::back_inserter(intersections);
+    intersect_2(cv1, cv2, boost::make_function_output_iterator([&out_it, &cst_curve_end](const Intersect_type& res) {
+                  if(auto* pt = std::get_if<Intersect_point>(&res)) {
+                    *out_it++ = pt->first;
+                    return;
+                  }
+                  if(auto* cv = std::get_if<Intersect_curve>(&res)) {
+                    *out_it++ = cst_curve_end(*cv, ARR_MIN_END).value();
+                    *out_it++ = cst_curve_end(*cv, ARR_MAX_END).value();
+                    return;
+                  }
+                  CGAL_assertion(false && "Unexpected intersection type");
+                }));
+    return intersections;
+  }
+
+  static std::optional<Point_2> first_intersection(Execution_context& ctx) {
+    if(!ctx.top_inters.empty() && !ctx.bottom_inters.empty()) {
+      return ctx->compare_xy_2(ctx.top_inters.front(), ctx.bottom_inters.front()) == CGAL::SMALLER
+                 ? ctx.top_inters.front()
+                 : ctx.bottom_inters.front();
+    }
+    if(!ctx.top_inters.empty()) {
+      return ctx.top_inters.front();
+    }
+    if(!ctx.bottom_inters.empty()) {
+      return ctx.bottom_inters.front();
+    }
+    return std::nullopt;
+  }
+
+  /**
+   * @brief approximate strictly x-monotone curve segment that does not cross the y bounds.
+   *
+   * @precondition: The segment is either inbound or outbound the bbox in the given range.
+   * @param start the x-coordinate of the segment start(exclusive)
+   * @param end the x-coordinate of the segment end(exclusive)
+   * @param step the step to approximate the curve segment, negative values allowed.
+   * @returns true if this part of the curve is within the closed bbox
+   */
+  static void approximate_simple_curve_segment(Execution_context& ctx, const FT& start, const FT& end, double step) {
+    for(FT x = start + step; x < end; x += step) {
+      auto y = ctx.bounded_compute_y_at_x(ctx.curve, x);
+      if(!y.has_value()) {
+        // break as soon as there's no more intersections
+        break;
+      }
+      if(y == ctx->ymin() || y == ctx->ymax()) {
+        // The segment overlaps with the bbox edge. There's no need to insert a dummy point.
+        break;
+      }
+      *ctx.out_it++ = ctx->approx_pt(Point_2(x, y.value()));
+      if(y > ctx->ymax() || y < ctx->ymin()) {
+        // We are outside the bbox. The dummy point is inserted to indicate the curve is outside the bbox.
+        break;
+      }
+    }
+  };
+
+  static void approximate_vertical_curve(Execution_context& ctx) {
+    if(ctx.is_bounded_curve() && ctx.min_end->x() < ctx->xmin() && ctx.min_end->x() > ctx->xmax()) {
+      // The curve is outside the bbox in x direction, no need to approximate
+      return;
+    }
+    if(!ctx.is_bounded_curve() && !ctx.has_y_intersections()) {
+      // The curve has unbounded end and has no intersections with the top or bottom edges,
+      // it must be outbound in x direction.
+      return;
+    }
+    // The vertical curve is now within the x bounds.
+
+    FT tymin = ctx.tight_ymin;
+    FT tymax = ctx.tight_ymax;
+    if(tymax == tymin) {
+      // But the curve is not degenerate. So, either it has only one point within the bbox or it is
+      // entirely outside the bbox in y direction.
+      return;
+    }
+    // Now we gaurantee that the curve has at least two points within the bbox in y direction.
+    // We have to obtain the x coordinate of this vertical curve.
+    FT x = ctx.is_bounded_curve() ? ctx.min_end->x() : first_intersection(ctx)->x();
+    *ctx.out_it++ = ctx->approx_pt(Point_2(x, tymin));
+    *ctx.out_it++ = ctx->approx_pt(Point_2(x, tymax));
+  }
+
+public:
+  Arr_bounded_approximate_curve_2(const Arr_bounded_render_context& ctx,
+                                  const Arr_bounded_approximate_point_2& point_approx)
+      : m_bounded_compute_y_at_x(ctx)
+      , m_approx_pt(point_approx)
+      , m_ctx(ctx)
+      , m_top(ctx.cst_horizontal_segment(ctx.ymax(), ctx.xmin(), ctx.xmax()))
+      , m_bottom(ctx.cst_horizontal_segment(ctx.ymin(), ctx.xmin(), ctx.xmax())) {}
+
+  /**
+   * @brief Approximate an x-monotone curve from left to right within the bounding box.
+   *
+   * @param he non-fictitious halfedge handle
+   * @return const Polyline_geom&
+   */
+  const Polyline_geom& operator()(const Halfedge_const_handle& he) const {
+    CGAL_assertion(!he->is_fictitious());
+
+    auto [polyline, inserted] = m_ctx.cache.try_emplace(he);
+    if(!inserted) {
+      return polyline;
+    }
+
+    if(m_ctx.is_cancelled()) {
+      return polyline;
+    }
+
+    const X_monotone_curve_2& curve = he->curve();
+
+    if(curve.is_degenerate()) {
+      return polyline;
+    }
+
+    auto top_inters = compute_intersections(curve, m_top, m_ctx.intersect_2, m_ctx.cst_curve_end);
+    auto bottom_inters = compute_intersections(curve, m_bottom, m_ctx.intersect_2, m_ctx.cst_curve_end);
+    Execution_context ctx(m_ctx, curve, m_approx_pt, m_bounded_compute_y_at_x, top_inters, bottom_inters, polyline);
+
+    if(ctx->is_vertical_2(curve)) {
+      approximate_vertical_curve(ctx);
+      return polyline;
+    }
+
+    polyline.reserve(top_inters.size() + bottom_inters.size());
+
+    FT txmin = ctx.tight_xmin;
+    FT txmax = ctx.tight_xmax;
+    FT last_x;
+    std::optional<Point_2> first_inter = first_intersection(ctx);
+
+    if(auto y_at_txmin = ctx.bounded_compute_y_at_x(curve, txmin);
+       y_at_txmin.has_value() && y_at_txmin != ctx->ymin() && y_at_txmin != ctx->ymax())
+    {
+      // The tight starting point of the curve is within the bbox and
+      // it's not on the top or bottom edge.
+      *ctx.out_it++ = txmin == ctx->xmin() ? ctx->approx_pt_on_boundary(Point_2(txmin, y_at_txmin.value()))
+                                           : ctx->approx_pt(Point_2(txmin, y_at_txmin.value()));
+      FT segment_end = first_inter.has_value() ? first_inter->x() : txmax;
+      approximate_simple_curve_segment(ctx, txmin, segment_end, ctx->approx_error);
+      last_x = segment_end;
+    } else if(first_inter.has_value()) {
+      last_x = first_inter->x();
+    } else {
+      // We assert that the curve is outbound.
+      // If the min end is bounded, it's obvious.
+      //
+      // If the min end is unbounded, we know that the curve has no intersections with top, bottom or left edge(txmin ==
+      // xmin when the min end is unbounded) and the min end of the curve is outbound (it approaches infinity in one or
+      // both dimension).
+      // Assume that the curve does has one point within the bbox. Note that it's a contiguous
+      // x-monotone curve. So it must cross the top, bottom or left edge to reach the min end from right to left,
+      // which is a contradiction.
+      return polyline;
+    }
+
+    // iterate through the intersections and insert segments in-between.
+    std::merge(top_inters.begin(), top_inters.end(), bottom_inters.begin(), bottom_inters.end(),
+               boost::make_function_output_iterator([&last_x, &ctx](const Point_2& inter) {
+                 approximate_simple_curve_segment(ctx, last_x, inter.x(), ctx->approx_error);
+                 *ctx.out_it++ = ctx->approx_pt_on_boundary(inter);
+                 last_x = inter.x();
+               }),
+               [&ctx](const Point_2& pt1, const Point_2& pt2) { return ctx->compare_xy_2(pt1, pt2) == CGAL::SMALLER; });
+
+    if(auto y_at_txmax = ctx.bounded_compute_y_at_x(curve, txmax);
+       y_at_txmax.has_value() && y_at_txmax != ctx->ymin() && y_at_txmax != ctx->ymax())
+    {
+      approximate_simple_curve_segment(ctx, last_x, txmax, ctx->approx_error);
+      *ctx.out_it++ = txmax == ctx->xmax() ? ctx->approx_pt_on_boundary(Point_2(txmax, y_at_txmax.value()))
+                                           : ctx->approx_pt(Point_2(txmax, y_at_txmax.value()));
+    }
+
+    return polyline;
+  }
+
+private:
+  const Arr_bounded_render_context& m_ctx;
+  const Arr_bounded_approximate_point_2& m_approx_pt;
+  const Arr_bounded_compute_y_at_x m_bounded_compute_y_at_x;
+  const X_monotone_curve_2 m_top;
+  const X_monotone_curve_2 m_bottom;
+};
+} // namespace CGAL
+#endif // CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_CURVE_2_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_face_2.h

@@ -0,0 +1,286 @@
+#ifndef CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_FACE_2_H
+#define CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_FACE_2_H
+
+#include "CGAL/Arr_enums.h"
+#include "CGAL/Bbox_2.h"
+#include "CGAL/Draw_aos/Arr_bounded_approximate_curve_2.h"
+#include "CGAL/Draw_aos/Arr_bounded_approximate_point_2.h"
+#include "CGAL/Draw_aos/Arr_bounded_face_triangulator.h"
+#include "CGAL/Draw_aos/Arr_render_context.h"
+#include "CGAL/basic.h"
+#include <CGAL/Draw_aos/Arr_approximation_geometry_traits.h>
+#include <CGAL/Draw_aos/helpers.h>
+#include <algorithm>
+#include <boost/iterator/function_output_iterator.hpp>
+#include <cstddef>
+#include <functional>
+#include <iostream>
+#include <iterator>
+#include <optional>
+#include <type_traits>
+
+namespace CGAL {
+
+namespace internal {
+
+/**
+ * @brief Patches corners between two boundary points of the bbox
+ * counter-clockwisely.
+ */
+class Patch_boundary
+{
+  using Approx_point = Arr_approximation_geometry_traits::Approx_point;
+
+  enum class Side_of_boundary {
+    Top = 0,
+    Left = 1,
+    Bottom = 2,
+    Right = 3,
+    None = -1,
+  };
+
+private:
+  Side_of_boundary side_of_boundary(const Approx_point& pt) const {
+    if(m_bbox.xmin() <= pt.x() && pt.x() < m_bbox.xmax() && pt.y() == m_bbox.ymax()) {
+      return Side_of_boundary::Top;
+    } else if(pt.x() == m_bbox.xmin() && m_bbox.ymin() <= pt.y() && pt.y() < m_bbox.ymax()) {
+      return Side_of_boundary::Left;
+    } else if(m_bbox.xmin() < pt.x() && pt.x() <= m_bbox.xmax() && pt.y() == m_bbox.ymin()) {
+      return Side_of_boundary::Bottom;
+    } else if(pt.x() == m_bbox.xmax() && m_bbox.ymin() < pt.y() && pt.y() <= m_bbox.ymax()) {
+      return Side_of_boundary::Right;
+    } else {
+      return Side_of_boundary::None;
+    }
+  }
+
+  Approx_point corner_of_side(Side_of_boundary side) const {
+    switch(side) {
+    case Side_of_boundary::Top:
+      // return the top-left corner
+      return Approx_point(m_bbox.xmin(), m_bbox.ymax());
+    case Side_of_boundary::Left:
+      // return the bottom-left corner
+      return Approx_point(m_bbox.xmin(), m_bbox.ymin());
+    case Side_of_boundary::Bottom:
+      // return the bottom-right corner
+      return Approx_point(m_bbox.xmax(), m_bbox.ymin());
+    case Side_of_boundary::Right:
+      // return the top-right corner
+      return Approx_point(m_bbox.xmax(), m_bbox.ymax());
+    default:
+      CGAL_assertion(false && "Invalid side of rectangle");
+      return Approx_point();
+    }
+  }
+
+  Side_of_boundary prev_side(Side_of_boundary side) const {
+    CGAL_assertion(side != Side_of_boundary::None);
+    return static_cast<Side_of_boundary>((static_cast<int>(side) + 3) % 4);
+  }
+
+  Side_of_boundary next_side(Side_of_boundary side) const {
+    CGAL_assertion(side != Side_of_boundary::None);
+    return static_cast<Side_of_boundary>((static_cast<int>(side) + 1) % 4);
+  }
+
+  bool is_on_boundary(const Approx_point& pt) const { return side_of_boundary(pt) != Side_of_boundary::None; }
+
+  double distance_on_same_side(const Approx_point& pt1, const Approx_point& pt2) const {
+    return std::abs(pt1.x() - pt2.x()) + std::abs(pt1.y() - pt2.y());
+  }
+
+public:
+  Patch_boundary(Bbox_2 bbox)
+      : m_bbox(bbox) {}
+
+  template <typename OutputIterator>
+  void operator()(const Approx_point& from, const Approx_point& to, OutputIterator out_it) const {
+    auto from_side = side_of_boundary(from);
+    auto to_side = side_of_boundary(to);
+
+    if(from_side == Side_of_boundary::None || to_side == Side_of_boundary::None) {
+      return;
+    }
+
+    bool from_is_before_to = from_side == to_side && distance_on_same_side(from, corner_of_side(from_side)) <=
+                                                         distance_on_same_side(to, corner_of_side(to_side));
+    // Special case: if both on the same side and the current point is before
+    // the last point in counter-clockwise boundary of the bbox, We have to go
+    // around the bbox. Note that this includes when two points overlaps
+    int num_corners_to_patch =
+        from_is_before_to ? 4 : (static_cast<int>(to_side) - static_cast<int>(from_side) + 4) % 4;
+
+    for(int i = 0; i < num_corners_to_patch; ++i) {
+      Side_of_boundary side = static_cast<Side_of_boundary>((static_cast<int>(from_side) + i) % 4);
+      Approx_point corner = corner_of_side(side);
+      if(corner == from) {
+        continue;
+      }
+      *out_it++ = corner;
+    }
+    return;
+  }
+
+private:
+  const Bbox_2 m_bbox;
+};
+
+} // namespace internal
+
+/**
+ * @brief Bounded face approximation for arrangements.
+ * @note Member functions are not thread-safe.
+ */
+class Arr_bounded_approximate_face_2
+{
+  using Approx_geom_traits = Arr_approximation_geometry_traits;
+  using Face_const_handle = Arrangement::Face_const_handle;
+  using Halfedge_const_handle = Arrangement::Halfedge_const_handle;
+  using Vertex_const_handle = Arrangement::Vertex_const_handle;
+  using Polyline_geom = Approx_geom_traits::Polyline_geom;
+  using Ccb_halfedge_const_circulator = Arrangement::Ccb_halfedge_const_circulator;
+  using Approx_point = Approx_geom_traits::Approx_point;
+  using Patch_boundary = internal::Patch_boundary;
+  using Triangulated_face = Approx_geom_traits::Triangulated_face;
+
+  struct Left_to_right_tag
+  {};
+  struct Right_to_left_tag
+  {};
+
+private:
+  class Execution_context : public Arr_context_delegator<Arr_bounded_render_context>
+  {
+  public:
+    Execution_context(const Arr_bounded_render_context& ctx,
+                      Arr_bounded_face_triangulator& triangulator,
+                      const Arr_bounded_approximate_point_2& bounded_approx_pt,
+                      const Arr_bounded_approximate_curve_2& bounded_approx_curve,
+                      const Patch_boundary& patch_boundary)
+        : Arr_context_delegator(ctx)
+        , triangulator(triangulator)
+        , patch_boundary(patch_boundary)
+        , bounded_approx_pt(bounded_approx_pt)
+        , bounded_approx_curve(bounded_approx_curve) {}
+
+  public:
+    const Arr_bounded_approximate_point_2& bounded_approx_pt;
+    const Arr_bounded_approximate_curve_2& bounded_approx_curve;
+    const Patch_boundary& patch_boundary;
+    Arr_bounded_face_triangulator& triangulator;
+  };
+
+private:
+  static void approximate_vertex(Execution_context& ctx, const Vertex_const_handle& vh) {
+    if(vh->is_at_open_boundary()) {
+      return;
+    }
+    ctx.bounded_approx_pt(vh);
+  }
+
+  template <typename OutputIterator>
+  static void
+  approximate_halfedge_of_ccb(Execution_context& ctx, const Halfedge_const_handle& he, OutputIterator& out_it) {
+    if(he->is_fictitious()) {
+      return;
+    }
+
+    const Polyline_geom& polyline = ctx.bounded_approx_curve(he);
+    if(he->direction() == ARR_LEFT_TO_RIGHT) {
+      std::copy(polyline.begin(), polyline.end(), out_it);
+    } else {
+      std::copy(polyline.rbegin(), polyline.rend(), out_it);
+    }
+  }
+
+  template <typename Ccb_tag>
+  static void approximate_ccb(Execution_context& ctx, const Ccb_halfedge_const_circulator& start_circ) {
+    constexpr bool Is_outer_ccb = std::is_same_v<Ccb_tag, Outer_ccb_tag>;
+    auto ccb_constraint = ctx.triangulator.make_ccb_constraint<Ccb_tag>();
+    auto out_it = ccb_constraint.insert_iterator();
+
+    std::optional<Approx_point> ccb_last_pt, ccb_first_pt;
+    auto counter_clockwise_start_circ = Is_outer_ccb ? start_circ : Ccb_halfedge_const_circulator(start_circ->twin());
+    auto circ = counter_clockwise_start_circ;
+    do {
+      bool is_he_first_pt = true;
+      auto patch_out_it = boost::make_function_output_iterator([&](const Approx_point& pt) {
+        if(ccb_last_pt == pt || !ctx->contains(pt)) {
+          return;
+        }
+        if(is_he_first_pt && ccb_last_pt.has_value()) {
+          ctx.patch_boundary(ccb_last_pt.value(), pt, out_it);
+        }
+
+        *out_it++ = pt;
+        ccb_last_pt = pt;
+        if(!ccb_first_pt.has_value()) {
+          ccb_first_pt = pt;
+        }
+        is_he_first_pt = false;
+      });
+
+      approximate_halfedge_of_ccb(ctx, circ, patch_out_it);
+      approximate_vertex(ctx, circ->target());
+    } while(++circ != counter_clockwise_start_circ);
+
+    if(Is_outer_ccb && !ccb_first_pt.has_value()) {
+      *out_it++ = Approx_point(ctx->xmin(), ctx->ymin());
+      *out_it++ = Approx_point(ctx->xmax(), ctx->ymin());
+      *out_it++ = Approx_point(ctx->xmax(), ctx->ymax());
+      *out_it++ = Approx_point(ctx->xmin(), ctx->ymax());
+    }
+    if(ccb_first_pt.has_value() && ccb_first_pt != ccb_last_pt) {
+      // Close the ccb
+      ctx.patch_boundary(ccb_last_pt.value(), ccb_first_pt.value(), out_it);
+    }
+  }
+
+public:
+  Arr_bounded_approximate_face_2(const Arr_bounded_render_context& ctx,
+                                 const Arr_bounded_approximate_point_2& point_approx,
+                                 const Arr_bounded_approximate_curve_2& curve_approx)
+      : m_ctx(ctx)
+      , m_patch_boundary(ctx.bbox())
+      , m_point_approx(point_approx)
+      , m_curve_approx(curve_approx) {}
+
+  const Triangulated_face& operator()(const Face_const_handle& fh) const {
+    auto [triangulated_face, inserted] = m_ctx.cache.try_emplace(fh);
+    if(!inserted) {
+      return triangulated_face;
+    }
+
+    if(m_ctx.is_cancelled()) {
+      return triangulated_face;
+    }
+
+    CGAL_assertion_msg(!fh->is_fictitious(), "Cannot approximate a fictitious face.");
+    if(!fh->has_outer_ccb()) {
+      // The face is the unbounded face of bounded arrangements
+      return triangulated_face;
+    }
+
+    Arr_bounded_face_triangulator triangulator(m_ctx);
+    Execution_context ctx(m_ctx, triangulator, m_point_approx, m_curve_approx, m_patch_boundary);
+
+    approximate_ccb<Outer_ccb_tag>(ctx, fh->outer_ccb());
+    for(auto inner_ccb = fh->inner_ccbs_begin(); inner_ccb != fh->inner_ccbs_end(); ++inner_ccb) {
+      approximate_ccb<Inner_ccb_tag>(ctx, *inner_ccb);
+    }
+    for(auto isolated_vh = fh->isolated_vertices_begin(); isolated_vh != fh->isolated_vertices_begin(); ++isolated_vh) {
+      approximate_vertex(ctx, isolated_vh);
+    }
+
+    return triangulated_face = (std::move(triangulator));
+  }
+
+private:
+  const Arr_bounded_render_context& m_ctx;
+  const Arr_bounded_approximate_point_2& m_point_approx;
+  const Arr_bounded_approximate_curve_2& m_curve_approx;
+  const Patch_boundary m_patch_boundary;
+};
+} // namespace CGAL
+#endif // CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_FACE_2_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_approximate_point_2.h

@@ -0,0 +1,55 @@
+
+#ifndef CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_POINT_2_H
+#define CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_POINT_2_H
+
+#include "CGAL/Draw_aos/Arr_render_context.h"
+#include <CGAL/Draw_aos/helpers.h>
+#include <CGAL/Draw_aos/Arr_approximate_point_2.h>
+#include <CGAL/Draw_aos/Arr_approximation_geometry_traits.h>
+#include <utility>
+
+namespace CGAL {
+class Arr_bounded_approximate_point_2
+{
+  using Approx_kernel = Arr_approximation_geometry_traits::Approximation_kernel;
+  using Point_2 = Geom_traits::Point_2;
+  using Vertex_const_handle = Arrangement::Vertex_const_handle;
+  using Point_geom = Arr_approximation_geometry_traits::Point_geom;
+
+public:
+  Arr_bounded_approximate_point_2(const Arr_bounded_render_context& ctx)
+      : m_ctx(ctx) {}
+
+  /**
+   * @brief Approximate a vertex within the x-bounded range.
+   *
+   * The function uses cached values if available.
+   * @precondition: The vertex must have an associated point.
+   *
+   * If the point is outside the bounding box, it returns Point_geom{} and false.
+   * Otherwise, it returns the approximated point and true.
+   * @param vh the vertex handle
+   * @return std::pair<const Point_geom&, bool>
+   */
+  std::pair<const Point_geom&, bool> operator()(const Vertex_const_handle& vh) const {
+    const auto& pt = vh->point();
+
+    if(!m_ctx.strictly_contains(pt)) {
+      return {Point_geom(), false};
+    }
+
+    auto [point, inserted] = m_ctx.cache.try_emplace(vh);
+    if(!inserted) {
+      return {point, true};
+    }
+
+    return {point = m_ctx.approx_pt(pt), true};
+  }
+
+private:
+  const Arr_bounded_render_context& m_ctx;
+};
+
+} // namespace CGAL
+
+#endif // CGAL_DRAW_AOS_ARR_BOUNDED_APPROXIMATE_POINT_2_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_compute_y_at_x.h

@@ -0,0 +1,79 @@
+#ifndef CGAL_DRAW_AOS_ARR_BOUNDED_COMPUTE_Y_AT_X_H
+#define CGAL_DRAW_AOS_ARR_BOUNDED_COMPUTE_Y_AT_X_H
+
+#include "CGAL/Arr_enums.h"
+#include "CGAL/Draw_aos/Arr_render_context.h"
+#include "CGAL/basic.h"
+#include <CGAL/Draw_aos/helpers.h>
+#include <CGAL/Draw_aos/Arr_construct_segments.h>
+#include <CGAL/Draw_aos/Arr_construct_curve_end.h>
+#include <boost/iterator/function_output_iterator.hpp>
+#include <optional>
+
+namespace CGAL {
+/**
+ * @brief Functor to compute the y-coordinate at a given x-coordinate for an x-monotone curve within a bounding box.
+ */
+class Arr_bounded_compute_y_at_x
+{
+public:
+  using Point_2 = Geom_traits::Point_2;
+  using X_monotone_curve_2 = Geom_traits::X_monotone_curve_2;
+  using Intersect_2 = Geom_traits::Intersect_2;
+  using Construct_min_vertex_2 = Geom_traits::Construct_min_vertex_2;
+  using FT = Geom_traits::FT;
+  using Approximate_2 = Geom_traits::Approximate_2;
+  using Is_vertical_2 = Geom_traits::Is_vertical_2;
+
+  Arr_bounded_compute_y_at_x(const Arr_bounded_render_context& ctx)
+      : m_ctx(ctx)
+      // TODO: some traits does not have approximate_2_object. we'll need a specialization for them.
+      , m_approx(ctx.traits.approximate_2_object()) {}
+
+  /**
+   * @brief Computes the y-coordinate at a given x-coordinate for an x-monotone curve trimmed
+   * to the bounding box.
+   *
+   * The bounding box here is considered as closed.
+   *
+   * @precondition The curve is not verical
+   * @param curve
+   * @param x
+   * @return true if there is an intersection at given x within the bounding box,
+   * @return false otherwise.
+   */
+  std::optional<FT> operator()(const X_monotone_curve_2& curve, const FT& x) const {
+    CGAL_assertion(!m_ctx.is_vertical_2(curve));
+    if(!m_ctx.contains_x(x)) {
+      return false;
+    }
+
+    auto min_pt = m_ctx.cst_curve_end(curve, ARR_MIN_END);
+    auto max_pt = m_ctx.cst_curve_end(curve, ARR_MAX_END);
+    if(min_pt.has_value() && min_pt->x() == x) {
+      return min_pt->y();
+    }
+    if(max_pt.has_value() && max_pt->x() == x) {
+      return max_pt->y();
+    }
+
+    using Multiplicity = Geom_traits::Multiplicity;
+    using Intersect_point = std::pair<Point_2, Multiplicity>;
+    using Intersect_curve = X_monotone_curve_2;
+    using Intersect_type = std::variant<Intersect_point, Intersect_curve>;
+
+    auto vertical_line = m_ctx.cst_vertical_segment(x, m_ctx.ymin(), m_ctx.ymax());
+    std::optional<FT> y;
+    auto func_out_iter = boost::make_function_output_iterator(
+        [&y, this](const Intersect_type& res) { y = std::get<Intersect_point>(res).first.y(); });
+    m_ctx.intersect_2(curve, vertical_line, func_out_iter);
+    return y;
+  }
+
+private:
+  const Arr_bounded_render_context& m_ctx;
+  Approximate_2 m_approx;
+};
+
+} // namespace CGAL
+#endif // CGAL_DRAW_AOS_ARR_BOUNDED_COMPUTE_Y_AT_X_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_face_triangulator.h

@@ -0,0 +1,297 @@
+#ifndef CGAL_DRAW_AOS_ARR_FACE_TRIANGULATOR_H
+#define CGAL_DRAW_AOS_ARR_FACE_TRIANGULATOR_H
+
+#include "CGAL/Constrained_Delaunay_triangulation_2.h"
+#include "CGAL/Constrained_triangulation_2.h"
+#include "CGAL/Constrained_triangulation_face_base_2.h"
+#include "CGAL/Draw_aos/Arr_approximation_geometry_traits.h"
+#include "CGAL/Draw_aos/Arr_render_context.h"
+#include "CGAL/Triangulation_vertex_base_with_info_2.h"
+#include "CGAL/mark_domain_in_triangulation.h"
+#include "CGAL/unordered_flat_map.h"
+#include <cstddef>
+#include <CGAL/Draw_aos/helpers.h>
+#include <boost/iterator/function_output_iterator.hpp>
+#include <boost/iterator/transform_iterator.hpp>
+#include <fstream>
+#include <iostream>
+#include <utility>
+#include <vector>
+
+namespace CGAL {
+
+/**
+ * @brief Triangulator a bounded face of an arrangement
+ *
+ * @note The face must have an outer CCB.
+ */
+class Arr_bounded_face_triangulator
+{
+  using Approx_geom_traits = Arr_approximation_geometry_traits;
+  using Approx_kernel = Approx_geom_traits::Approximation_kernel;
+  using Approx_point = Approx_geom_traits::Approx_point;
+  using Point_vec = Approx_geom_traits::Apporx_point_vec;
+  using Triangle = Approx_geom_traits::Triangle;
+  using Triangulated_face = Approx_geom_traits::Triangulated_face;
+
+  struct Point_index
+  {
+    constexpr static std::size_t Invalid_index = std::numeric_limits<std::size_t>::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; }
+  };
+  using Epick = CGAL::Exact_predicates_inexact_constructions_kernel;
+  using Vb = CGAL::Triangulation_vertex_base_with_info_2<Point_index, Epick>;
+  using Fb = CGAL::Constrained_triangulation_face_base_2<Epick>;
+  using Tds = CGAL::Triangulation_data_structure_2<Vb, Fb>;
+  using Ct = Constrained_triangulation_2<Epick, Tds, Exact_predicates_tag>;
+  using Point = Epick::Point_2;
+  using Point_with_info = std::pair<Point, Point_index>;
+
+  std::size_t counter{0};
+
+  /**
+   * @brief RAII-style inserter for one CCB in a triangulation.
+   * Collects points and inserts them as a constraint on destruction.
+   * Only one instance per Arr_face_triangulator is allowed at a time.
+   */
+  template <typename Ccb_tag>
+  class Ccb_constraint
+  {
+    constexpr static bool Is_outer_ccb = std::is_same_v<Ccb_tag, Outer_ccb_tag>;
+    friend class Arr_bounded_face_triangulator;
+
+    std::ofstream m_ofs;
+    Ccb_constraint(Arr_bounded_face_triangulator& triangulator)
+        : m_triangulator(&triangulator) {
+      m_ccb_start = m_triangulator->m_points.size();
+      triangulator.m_has_active_constraint = true;
+
+      if(Is_outer_ccb) {
+        std::ofstream ofs_index("/Users/shep/codes/aos_2_js_helper/shapes.txt", std::ios::app);
+        auto name = "outer_ccb_" + std::to_string((*m_triangulator->m_ctx.counter)++) + ".txt";
+        ofs_index << name << std::endl;
+        m_ofs = std::ofstream("/Users/shep/codes/aos_2_js_helper/" + name, std::ios::out | std::ios::trunc);
+      }
+    }
+
+  private:
+    Point offset_boundary_point(const Point& pt) const {
+      constexpr double offset = 1; // It doesn't matter how much we offset the point
+      double x = pt.x(), y = pt.y();
+      const auto& ctx = m_triangulator->m_ctx;
+
+      if(x == ctx.xmin()) {
+        x -= offset;
+      }
+      if(x == ctx.xmax()) {
+        x += offset;
+      }
+      if(y == ctx.ymin()) {
+        y -= offset;
+      }
+      if(y == ctx.ymax()) {
+        y += offset;
+      }
+      return Point(x, y);
+    }
+
+    void insert_ccb() {
+      auto& ct = m_triangulator->m_ct;
+      auto& points = m_triangulator->m_points;
+
+      auto begin = points.begin() + m_ccb_start;
+      auto end = points.end();
+
+      if constexpr(Is_outer_ccb) {
+        auto helpers_iter = m_helper_indices.begin();
+        auto helpers_end = m_helper_indices.end();
+
+        auto concrete_pt_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 = [&points](std::size_t idx) { return points[idx]; };
+
+        auto indexes_begin = boost::make_counting_iterator<std::size_t>(m_ccb_start);
+        auto indexes_end = boost::make_counting_iterator<std::size_t>(points.size());
+        auto filtered_begin = boost::make_filter_iterator(concrete_pt_filter, indexes_begin, indexes_end);
+        auto filtered_end = boost::make_filter_iterator(concrete_pt_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);
+        // {
+        //   std::ofstream ofs_index("/Users/shep/codes/aos_2_js_helper/shapes.txt", std::ios::app);
+        //   auto& ctx = const_cast<Arr_bounded_render_context&>(m_triangulator->m_ctx);
+        //   ofs_index << "ccb_" << ctx.counter << ".txt" << std::endl;
+        //   std::ofstream ofs("/Users/shep/codes/aos_2_js_helper/ccb_" + std::to_string(ctx.counter++) + ".txt",
+        //                     std::ios::out | std::ios::trunc);
+        //   for(auto it = filtered_begin; it != filtered_end; ++it) {
+        //     const auto& pt = it->first;
+        //     ofs << pt.x() << " " << pt.y() << std::endl;
+        //   }
+        // }
+        ct.insert_with_info<Point_with_info>(transformed_begin, transformed_end);
+      } else {
+        ct.insert_with_info<Point_with_info>(begin, end);
+      }
+    }
+
+  public:
+    Ccb_constraint(const Ccb_constraint&) = delete;
+    Ccb_constraint& operator=(const Ccb_constraint&) = delete;
+    Ccb_constraint(Ccb_constraint&& other) noexcept
+        : m_triangulator(other.m_triangulator) {
+      other.m_triangulator = nullptr;
+    }
+    Ccb_constraint& operator=(Ccb_constraint&& other) noexcept {
+      if(this != &other) {
+        m_triangulator = other.m_triangulator;
+        other.m_triangulator = nullptr;
+      }
+      return *this;
+    }
+
+    decltype(auto) insert_iterator() {
+      return boost::make_function_output_iterator([&, this](const Approx_point& pt) {
+        auto& points = m_triangulator->m_points;
+        CGAL_assertion_msg(m_triangulator != nullptr, "Use of destructed or moved Ccb_constraint object.");
+        CGAL_assertion_msg(m_triangulator->m_ctx.contains(pt), "Outbound point in Ccb_constraint.");
+
+        if constexpr(Is_outer_ccb) {
+          bool is_on_boundary = m_triangulator->m_ctx.is_on_boundary(pt);
+          if(is_on_boundary && m_is_last_on_boundary) {
+            m_helper_indices.push_back(points.size());
+            points.emplace_back(offset_boundary_point(points.back().first), Point_index());
+          }
+          m_is_last_on_boundary = is_on_boundary;
+          m_ofs << pt.x() << " " << pt.y() << std::endl;
+        }
+        points.emplace_back(Point(pt.x(), pt.y()), points.size());
+      });
+    }
+
+    ~Ccb_constraint() {
+      if(m_triangulator == nullptr) {
+        return;
+      }
+      if(Is_outer_ccb && m_is_last_on_boundary && m_triangulator->m_points.size() != m_ccb_start) {
+        auto& points = m_triangulator->m_points;
+        const auto& first_pt = points[m_ccb_start].first;
+        const auto& last_pt = points.back().first;
+        bool is_first_on_boundary = m_triangulator->m_ctx.is_on_boundary(first_pt);
+        if(m_is_last_on_boundary && is_first_on_boundary) {
+          m_helper_indices.push_back(points.size());
+          points.emplace_back(offset_boundary_point(last_pt), Point_index());
+        }
+      }
+
+      insert_ccb();
+      m_triangulator->m_has_active_constraint = false;
+      m_triangulator = nullptr;
+    }
+
+  private:
+    Arr_bounded_face_triangulator* m_triangulator;
+    std::size_t m_ccb_start;
+    bool m_is_last_on_boundary = false;
+    std::vector<std::size_t> m_helper_indices; // The offseted point indices when inserting outer ccb constraint
+  };
+
+public:
+  Arr_bounded_face_triangulator(const Arr_bounded_render_context& ctx)
+      : m_ctx(ctx) {}
+  Arr_bounded_face_triangulator(const Arr_bounded_face_triangulator& other) = delete;
+  Arr_bounded_face_triangulator& operator=(const Arr_bounded_face_triangulator& other) = delete;
+
+  operator Triangulated_face() && {
+    CGAL_assertion(!m_has_active_constraint && "There is an active constraint in the triangulator.");
+    CGAL_assertion(m_outer_ccb_processed && "Outer CCB has not been processed yet.");
+
+    // insert_constraint() should be called after insert_with_info(), or info will not be set correctly.
+    auto first_of_pair = [](const Point_with_info& pt_with_info) { return pt_with_info.first; };
+    for(std::size_t i = 0; i < m_ccb_start_indices.size(); ++i) {
+      auto begin = m_points.begin() + m_ccb_start_indices[i];
+      auto end = m_points.begin() + (i + 1 < m_ccb_start_indices.size() ? m_ccb_start_indices[i + 1] : m_points.size());
+      {
+        std::ofstream ofs_index("/Users/shep/codes/aos_2_js_helper/shapes.txt", std::ios::app);
+        auto& ctx = const_cast<Arr_bounded_render_context&>(m_ctx);
+        auto name = "ccb_constraint" + std::to_string((*ctx.counter)++) + ".txt";
+        ofs_index << name << std::endl;
+        std::ofstream ofs("/Users/shep/codes/aos_2_js_helper/" + name, std::ios::out | std::ios::trunc);
+        for(auto it = begin; it != end; ++it) {
+          const auto& pt = it->first;
+          ofs << pt.x() << " " << pt.y() << std::endl;
+        }
+      }
+      m_ct.insert_constraint(boost::make_transform_iterator(begin, first_of_pair),
+                             boost::make_transform_iterator(end, first_of_pair), true);
+    }
+
+    unordered_flat_map<Ct::Face_handle, bool> in_domain_map;
+    boost::associative_property_map<decltype(in_domain_map)> in_domain(in_domain_map);
+    CGAL::mark_domain_in_triangulation(m_ct, in_domain);
+
+    Triangulated_face tf;
+    tf.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();
+      if(!v1.is_valid() || !v2.is_valid() || !v3.is_valid()) {
+        continue;
+      }
+      Triangle tri{v1, v2, v3};
+      tf.triangles.emplace_back(tri);
+    }
+
+    auto transform_first_of_pair = [](const Point_with_info& pt_with_info) {
+      return Approx_point(pt_with_info.first.x(), pt_with_info.first.y());
+    };
+    tf.points = Point_vec(boost::make_transform_iterator(m_points.begin(), transform_first_of_pair),
+                          boost::make_transform_iterator(m_points.end(), transform_first_of_pair));
+    m_ct.clear();
+    return tf;
+  }
+
+  /**
+   * @brief Get an constraint object for a certain ccb.
+   * @note Only one Ccb_constraint can be active per triangulator at a certain time point. When it goes out of scope,
+   * the points collected will be inserted as a constraint into the triangulation.
+   */
+  template <typename Ccb_tag>
+  Ccb_constraint<Ccb_tag> make_ccb_constraint() {
+    if constexpr(std::is_same_v<Ccb_tag, Outer_ccb_tag>) {
+      CGAL_assertion_msg(!m_outer_ccb_processed, "Outer CCB has already been processed.");
+      m_outer_ccb_processed = true;
+    }
+    CGAL_assertion_msg(!m_has_active_constraint, "Only one Ccb_constraint can be active per triangulator at a time.");
+    m_ccb_start_indices.push_back(m_points.size());
+    return Ccb_constraint<Ccb_tag>(*this);
+  }
+
+private:
+  const Arr_bounded_render_context& m_ctx;
+  Ct m_ct;
+  std::vector<Point_with_info> m_points;
+  std::vector<std::size_t> m_ccb_start_indices;
+  bool m_has_active_constraint = false;
+  bool m_outer_ccb_processed = false;
+};
+
+} // namespace CGAL
+
+#endif // CGAL_DRAW_AOS_ARR_FACE_TRIANGULATOR_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_bounded_renderer.h

@@ -0,0 +1,216 @@
+#ifndef CGAL_DRAW_AOS_ARR_BOUNDED_RENDERER_H
+#define CGAL_DRAW_AOS_ARR_BOUNDED_RENDERER_H
+
+#include "CGAL/Arr_trapezoid_ric_point_location.h"
+#include "CGAL/Bbox_2.h"
+#include "CGAL/Draw_aos/Arr_approximation_cache.h"
+#include "CGAL/Draw_aos/Arr_approximation_geometry_traits.h"
+#include "CGAL/Draw_aos/Arr_bounded_approximate_curve_2.h"
+#include "CGAL/Draw_aos/Arr_bounded_approximate_face_2.h"
+#include "CGAL/Draw_aos/Arr_render_context.h"
+#include <CGAL/Draw_aos/helpers.h>
+#include <CGAL/Draw_aos/Arr_portals.h>
+#include <CGAL/IO/Color.h>
+#include <vector>
+
+namespace CGAL {
+/**
+ * @brief Render arrangement on surface within a bounding box.
+ *
+ * @note The class is not thread-safe.
+ */
+class Arr_bounded_renderer
+{
+  using Color = IO::Color;
+  using Point_2 = Geom_traits::Point_2;
+  using Vertex_const_handle = Arrangement::Vertex_const_handle;
+  using Edge_const_handle = Arrangement::Edge_const_iterator;
+  using Halfedge_const_handle = Arrangement::Halfedge_const_handle;
+  using Face_const_handle = Arrangement::Face_const_handle;
+  using Vertex_handle = Arrangement::Vertex_handle;
+  using Halfedge_handle = Arrangement::Halfedge_handle;
+  using Face_handle = Arrangement::Face_handle;
+  using X_monotone_curve_2 = Geom_traits::X_monotone_curve_2;
+  using Point_Location = Arr_trapezoid_ric_point_location<Arrangement>;
+  using Feature_const = std::variant<Vertex_const_handle, Halfedge_const_handle, Face_const_handle>;
+  using Feature_const_vector = std::vector<Feature_const>;
+  using Feature_const_vector_inserter = std::back_insert_iterator<Feature_const_vector>;
+  using Arr_approx_geom_traits = Arr_approximation_geometry_traits;
+  using Point_geom = Arr_approx_geom_traits::Point_geom;
+  using Polyline_geom = Arr_approx_geom_traits::Polyline_geom;
+  using Feature_portal_map = Arr_portals::Feature_portals_map;
+  // QFlags implement this pattern better, but we try not to reply on Qt classes.
+  template <typename E>
+  class Arr_flags
+  {
+  public:
+    Arr_flags(const E& flags)
+        : m_flags(static_cast<int>(flags)) {}
+    Arr_flags(std::initializer_list<E> flags)
+        : m_flags(0) {
+      for(const auto& flag : flags) {
+        m_flags |= static_cast<int>(flag);
+      }
+    }
+
+    bool is_set(E flag) const { return (m_flags & static_cast<int>(flag)) != 0; }
+
+  private:
+    int m_flags;
+  };
+
+  enum class Feature_type {
+    Vertex = 1 << 0,
+    Halfedge = 1 << 1,
+    Face = 1 << 2,
+  };
+
+  struct Execution_context : Arr_context_delegator<Arr_bounded_render_context>
+  {
+    Execution_context(const Arr_bounded_render_context& ctx)
+        : Arr_context_delegator(ctx)
+        , bounded_approx_pt(ctx)
+        , bounded_approx_curve(ctx, bounded_approx_pt)
+        , bounded_approx_face(ctx, bounded_approx_pt, bounded_approx_curve) {}
+
+    const Arr_bounded_approximate_point_2 bounded_approx_pt;
+    const Arr_bounded_approximate_curve_2 bounded_approx_curve;
+    const Arr_bounded_approximate_face_2 bounded_approx_face;
+  };
+
+private:
+  template <typename OutputIterator>
+  static void locate_intersecting_features(const Execution_context& ctx,
+                                           const X_monotone_curve_2& curve,
+                                           OutputIterator out_iter,
+                                           Arr_flags<Feature_type> feats) {
+    using Feature = std::variant<Vertex_handle, Halfedge_handle, Face_handle>;
+    using Feature_vector = std::vector<Feature>;
+
+    auto func_out_iter = boost::make_function_output_iterator([&out_iter, &feats](const Feature& feature) {
+      if(auto* vh = std::get_if<Vertex_handle>(&feature)) {
+        if(!feats.is_set(Feature_type::Vertex) || (*vh)->is_at_open_boundary()) {
+          return;
+        }
+        *out_iter++ = Vertex_const_handle(*vh);
+      } else if(auto* he = std::get_if<Halfedge_handle>(&feature)) {
+        if(!feats.is_set(Feature_type::Halfedge) || (*he)->is_fictitious()) {
+          return;
+        }
+        *out_iter++ = Halfedge_const_handle(*he);
+      } else if(auto* fh = std::get_if<Face_handle>(&feature)) {
+        if(!feats.is_set(Feature_type::Face)) {
+          return;
+        }
+        *out_iter++ = Face_const_handle(*fh);
+      } else {
+        CGAL_assertion(false && "Unexpected feature type");
+      }
+    });
+
+    using Zone_visitor = Arr_compute_zone_visitor<Arrangement, decltype(func_out_iter)>;
+    using Zone_2 = Arrangement_zone_2<Arrangement, Zone_visitor>;
+
+    Zone_visitor zone_visitor(func_out_iter);
+    // use const_cast to fit in the interface but we will not modify the arrangement
+    Zone_2 zone(const_cast<Arrangement&>(ctx->arr), &zone_visitor);
+    zone.init(curve, ctx->point_location);
+    zone.compute_zone();
+  }
+
+  static void discover_faces(Execution_context& ctx, const Face_const_handle& fh) {
+    if(ctx->cache.has(fh)) {
+      // The face has already been discovered and approximated.
+      return;
+    }
+    ctx.bounded_approx_face(fh);
+
+    for(auto inner_ccb = fh->inner_ccbs_begin(); inner_ccb != fh->inner_ccbs_end(); ++inner_ccb) {
+      auto he = *inner_ccb;
+      auto inner_face = he->twin()->face();
+
+      bool is_degenerate = inner_face == fh;
+      bool within_bounds = ctx->strictly_contains(inner_face->outer_ccb()->source()->point());
+      if(is_degenerate || !within_bounds) {
+        continue;
+      }
+      discover_faces(ctx, inner_face);
+    }
+
+    // We don't need to handle isolated vertices.
+    // for(...)
+
+    if(!fh->has_outer_ccb()) {
+      // The unbounded face of bounded arrangements has no outer CCB.
+      return;
+    }
+    // find fully contained adjacent faces
+    auto outer_ccb = fh->outer_ccb();
+    auto circ = outer_ccb;
+    do {
+      auto adj_face = circ->twin()->face();
+      if(adj_face->is_fictitious() || adj_face->is_unbounded()) {
+        // Unbounded faces cannot be contained in out bounding box.
+        continue;
+      }
+
+      bool within_bounds = ctx->strictly_contains(adj_face->outer_ccb()->source()->point());
+      if(!within_bounds) {
+        continue;
+      }
+      // For a face that is not one of the seeding faces,
+      // the face must be contained iff one of its vertices is contained.
+      // For a face that is, we'll touch it sooner or later.
+      discover_faces(ctx, adj_face);
+    } while(++circ != outer_ccb);
+  }
+
+public:
+  Arr_bounded_renderer(Arr_render_context& ctx, Bbox_2 bbox)
+      : m_ctx(ctx)
+      , m_bbox(bbox) {}
+
+  Arr_approximation_cache render() const {
+    Arr_approximation_cache cache;
+
+    if(m_ctx.is_cancelled()) {
+      return cache;
+    }
+
+    Execution_context ctx(Arr_bounded_render_context(m_ctx, m_bbox, cache));
+
+    auto insert_features = boost::make_function_output_iterator([&ctx](const Feature_const& feature) {
+      if(auto* vh = std::get_if<Vertex_const_handle>(&feature)) {
+        ctx.bounded_approx_pt(*vh);
+      } else if(auto* he = std::get_if<Halfedge_const_handle>(&feature)) {
+        ctx.bounded_approx_curve(*he);
+      } else if(auto* fh = std::get_if<Face_const_handle>(&feature)) {
+        discover_faces(ctx, *fh);
+      } else {
+        CGAL_assertion(false && "Unexpected feature type");
+      }
+    });
+
+    auto top = ctx->cst_horizontal_segment(ctx->ymax(), ctx->xmin(), ctx->xmax());
+    auto bottom = ctx->cst_horizontal_segment(ctx->ymin(), ctx->xmin(), ctx->xmax());
+    auto left = ctx->cst_vertical_segment(ctx->xmin(), ctx->ymin(), ctx->ymax());
+    auto right = ctx->cst_vertical_segment(ctx->xmax(), ctx->ymin(), ctx->ymax());
+
+    // top, left are open edges while bottom, right are closed.
+    locate_intersecting_features(ctx, top, insert_features, Feature_type::Face);
+    locate_intersecting_features(ctx, left, insert_features, Feature_type::Face);
+    locate_intersecting_features(ctx, bottom, insert_features, {Feature_type::Face, Feature_type::Halfedge});
+    locate_intersecting_features(ctx, right, insert_features,
+                                 {Feature_type::Face, Feature_type::Halfedge, Feature_type::Vertex});
+
+    return cache;
+  }
+
+private:
+  Arr_render_context m_ctx;
+  const Bbox_2 m_bbox;
+};
+
+} // namespace CGAL
+
+#endif // CGAL_DRAW_AOS_ARR_BOUNDED_RENDERER_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_construct_curve_end.h

@@ -0,0 +1,101 @@
+#ifndef CGAL_ARR_CONSTRUCT_CURVE_END_H
+#define CGAL_ARR_CONSTRUCT_CURVE_END_H
+
+#include "CGAL/Arr_enums.h"
+#include "CGAL/Arr_has.h"
+#include "CGAL/Draw_aos/Arr_approximation_geometry_traits.h"
+#include <optional>
+
+namespace CGAL {
+
+namespace internal {
+
+template <typename Geom_traits>
+class Arr_construct_curve_end_base
+{
+  using Construct_min_vertex_2 = typename Geom_traits::Construct_min_vertex_2;
+  using Construct_max_vertex_2 = typename Geom_traits::Construct_max_vertex_2;
+
+protected:
+  Arr_construct_curve_end_base(const Geom_traits& traits)
+      : m_cst_min_vertex(traits.construct_min_vertex_2_object())
+      , m_cst_max_vertex(traits.construct_max_vertex_2_object()) {}
+
+  Construct_max_vertex_2 m_cst_max_vertex;
+  Construct_min_vertex_2 m_cst_min_vertex;
+};
+
+template <typename Geom_traits, bool Has_parameter_space_in_x>
+class Arr_construct_curve_end_impl;
+
+template <typename Geom_traits>
+class Arr_construct_curve_end_impl<Geom_traits, true> : public Arr_construct_curve_end_base<Geom_traits>
+{
+  using Approx_geom_traits = Arr_approximation_geometry_traits;
+  using Point_2 = typename Geom_traits::Point_2;
+  using X_monotone_curve = typename Geom_traits::X_monotone_curve_2;
+  using Parameter_space_in_x_2 = typename Geom_traits::Parameter_space_in_x_2;
+  using Parameter_space_in_y_2 = typename Geom_traits::Parameter_space_in_y_2;
+
+public:
+  Arr_construct_curve_end_impl(const Geom_traits& traits)
+      : Arr_construct_curve_end_base<Geom_traits>(traits)
+      , m_param_space_in_x(traits.parameter_space_in_x_2_object())
+      , m_param_space_in_y(traits.parameter_space_in_y_2_object()) {}
+
+  /**
+   * @brief Returns an optional type contains the approximated point the given end of the curve is bounded.
+   * Otherwise it's std::nullopt.
+   *
+   * @param curve
+   * @param ce
+   * @return std::optional<Point_2>
+   */
+  std::optional<Point_2> operator()(const X_monotone_curve& curve, Arr_curve_end ce) const {
+    auto ps_x = m_param_space_in_x(curve, ce);
+    auto ps_y = m_param_space_in_y(curve, ce);
+
+    if(ps_x == ARR_INTERIOR && ps_y == ARR_INTERIOR) {
+      return std::make_optional(ce == ARR_MIN_END ? this->m_cst_min_vertex(curve) : this->m_cst_max_vertex(curve));
+    }
+
+    return std::nullopt;
+  }
+
+private:
+  Parameter_space_in_x_2 m_param_space_in_x;
+  Parameter_space_in_y_2 m_param_space_in_y;
+};
+
+template <typename Geom_traits>
+class Arr_construct_curve_end_impl<Geom_traits, false> : public Arr_construct_curve_end_base<Geom_traits>
+{
+  using Approx_geom_traits = Arr_approximation_geometry_traits;
+  using Point_2 = typename Geom_traits::Point_2;
+  using X_monotone_curve = typename Geom_traits::X_monotone_curve_2;
+
+public:
+  Arr_construct_curve_end_impl(const Geom_traits& traits)
+      : Arr_construct_curve_end_base<Geom_traits>(traits) {}
+
+  /**
+   * @brief Construct the min or max vertex of the x-monotone curve based on the curve end enum.
+   *
+   * @param curve
+   * @param ce
+   * @return std::optional<Point_2>
+   */
+  std::optional<Point_2> operator()(const X_monotone_curve& curve, Arr_curve_end ce) const {
+    return ce == ARR_MIN_END ? this->m_cst_min_vertex(curve) : this->m_cst_max_vertex(curve);
+  }
+};
+
+} // namespace internal
+
+template <typename Geom_traits>
+using Arr_construct_curve_end =
+    internal::Arr_construct_curve_end_impl<Geom_traits, has_parameter_space_in_x_2<Geom_traits>::value>;
+
+} // namespace CGAL
+
+#endif // CGAL_ARR_CONSTRUCT_CURVE_END_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_construct_segments.h

@@ -0,0 +1,45 @@
+#ifndef CGAL_DRAW_AOS_ARR_CONSTRUCT_SEGMENTS_H
+#define CGAL_DRAW_AOS_ARR_CONSTRUCT_SEGMENTS_H
+#include "CGAL/Draw_aos/helpers.h"
+
+namespace CGAL {
+class Arr_construct_vertical_segment
+{
+  using Point_2 = Geom_traits::Point_2;
+  using X_monotone_curve_2 = Geom_traits::X_monotone_curve_2;
+  using FT = Geom_traits::FT;
+
+public:
+  Arr_construct_vertical_segment(const Geom_traits& traits)
+      : m_traits(traits) {}
+
+  X_monotone_curve_2 operator()(const FT& x, const FT& ymin, const FT& ymax) const {
+    auto cst_x_curve = m_traits.construct_x_monotone_curve_2_object();
+    return cst_x_curve(Point_2(x, ymin), Point_2(x, ymax));
+  }
+
+private:
+  const Geom_traits& m_traits;
+};
+
+class Arr_construct_horizontal_segment
+{
+  using Point_2 = Geom_traits::Point_2;
+  using X_monotone_curve_2 = Geom_traits::X_monotone_curve_2;
+  using FT = Geom_traits::FT;
+
+public:
+  Arr_construct_horizontal_segment(const Geom_traits& traits)
+      : m_traits(traits) {}
+
+  X_monotone_curve_2 operator()(FT y, FT xmin, FT xmax) const {
+    auto cst_x_curve = m_traits.construct_x_monotone_curve_2_object();
+    return cst_x_curve(Point_2(xmin, y), Point_2(xmax, y));
+  }
+
+private:
+  const Geom_traits& m_traits;
+};
+} // namespace CGAL
+
+#endif

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_graph_conn.h

@@ -0,0 +1,93 @@
+
+#ifndef CGAL_DRAW_AOS_ARR_GRAPH_CONN_H
+#define CGAL_DRAW_AOS_ARR_GRAPH_CONN_H
+#include "CGAL/Union_find.h"
+#include "CGAL/unordered_flat_map.h"
+#include <boost/range/iterator_range.hpp>
+namespace CGAL {
+/**
+ * @brief Arr_graph_conn provides fast connectivity queries for arrangement vertices
+ * based on union-find data structure.
+ */
+template <typename Arr>
+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<Vertex_const_handle>::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());
+    auto [target_it, target_inserted] = m_lookup.try_emplace(target, Union_find_handle());
+    if(source_inserted) {
+      source_it->second = m_uf.make_set(source);
+    }
+    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) {
+    for(const auto& eh : arr.edge_handles()) {
+      insert_halfedge(eh);
+    }
+
+    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<Vertex_const_handle> m_uf;
+  unordered_flat_map<Vertex_const_handle, Union_find_handle> m_lookup;
+};
+} // namespace CGAL
+#endif // CGAL_DRAW_AOS_ARR_GRAPH_CONN_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_portals.h

@@ -0,0 +1,133 @@
+#ifndef CGAL_DRAW_AOS_ARR_CREATE_PORTALS_H
+#define CGAL_DRAW_AOS_ARR_CREATE_PORTALS_H
+
+#include "CGAL/Arr_vertical_decomposition_2.h"
+#include "CGAL/Draw_aos/Arr_approximate_point_2.h"
+#include "CGAL/Draw_aos/Arr_construct_segments.h"
+#include "CGAL/Draw_aos/helpers.h"
+#include "CGAL/Draw_aos/Arr_graph_conn.h"
+#include "CGAL/Object.h"
+#include "CGAL/basic.h"
+#include "CGAL/unordered_flat_map.h"
+#include <boost/iterator/function_output_iterator.hpp>
+#include <CGAL/Draw_aos/helpers.h>
+#include <CGAL/Draw_aos/Arr_approximation_geometry_traits.h>
+#include <limits>
+#include <utility>
+
+namespace CGAL {
+/**
+ * @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.
+ */
+class Arr_portals
+{
+
+  using Vertex_const_handle = Arrangement::Vertex_const_handle;
+  using Halfedge_const_handle = Arrangement::Halfedge_const_handle;
+  using Face_const_handle = Arrangement::Face_const_handle;
+  using Point_2 = Geom_traits::Point_2;
+  using Approx_point = Arr_approximation_geometry_traits::Approx_point;
+  using X_monotone_curve_2 = Geom_traits::X_monotone_curve_2;
+  using Feature_const = std::variant<Vertex_const_handle, Halfedge_const_handle, Face_const_handle>;
+
+public:
+  // Pair composed of a portal entry point and the portaled vertex on the inner ccb.
+  using Portal = std::pair<std::optional<Approx_point>, Vertex_const_handle>;
+  using Portal_vector = std::vector<Portal>;
+  // Map from a feature to its portals sorted by the x coordinate of the virtual vertical segments.
+  using Feature_portals_map = unordered_flat_map<Feature_const, Portal_vector>;
+
+private:
+  // Use this function to locate the intersection point of a vertical ray shooted from a point
+  // to it's upper x-monotone curve(it can't be vertical, or we should've got its min vertex from vertical
+  // decomposition).
+  static Point_2 upper_intersection(const Geom_traits& traits, const Point_2& pt, const X_monotone_curve_2& curve) {
+    Arr_construct_vertical_segment cst_vertical_segment(traits);
+    Point_2 intersection_point;
+    auto intersect = traits.intersect_2_object();
+    auto vertical_line = cst_vertical_segment(pt.x(), pt.y(), std::numeric_limits<double>::max());
+    bool found_intersection = false;
+
+    using Multiplicity = Geom_traits::Multiplicity;
+    using Intersect_point = std::pair<Point_2, Multiplicity>;
+    using Intersect_curve = X_monotone_curve_2;
+    using Intersect_type = std::variant<Intersect_point, Intersect_curve>;
+
+    intersect(curve, vertical_line, boost::make_function_output_iterator([&](const Intersect_type& res) {
+                found_intersection = true;
+                if(std::holds_alternative<Intersect_point>(res)) {
+                  intersection_point = std::get<Intersect_point>(res).first;
+                  return;
+                }
+                CGAL_assertion(false && "Unexpected intersection type");
+              }));
+
+    return found_intersection ? intersection_point : Point_2(pt.x(), std::numeric_limits<double>::max());
+  }
+
+public:
+  Feature_portals_map create(const Arrangement& arr) const {
+    using Object_pair = std::pair<Object, Object>;
+    using Vert_decomp_entry = std::pair<Vertex_const_handle, Object_pair>;
+
+    Arr_graph_conn conn(arr);
+    auto visited_ccbs = std::unordered_set<Vertex_const_handle>();
+    Feature_portals_map feature_portals;
+    auto intersect = arr.traits()->intersect_2_object();
+    auto approx_pt = Arr_approximate_point_2<Geom_traits>(*arr.traits());
+
+    auto func_out_iter = boost::make_function_output_iterator([&](const Vert_decomp_entry& entry) {
+      const auto& [vh, obj_pair] = entry;
+
+      const auto& ccb_main_vertex = conn.ccb_representative_vertex(vh);
+      if(visited_ccbs.find(ccb_main_vertex) != visited_ccbs.end()) {
+        // This ccb has already been processed
+        return;
+      }
+
+      const auto& above_feat = obj_pair.second;
+      if(Vertex_const_handle above_vh; CGAL::assign(above_vh, above_feat)) {
+        if(conn.is_connected(above_vh, vh)) {
+          // This upper vertex is connected to vh, skip it
+          return;
+        }
+
+        const auto& [it, _] = feature_portals.try_emplace(above_vh, std::vector<Portal>{});
+        if((above_vh)->is_at_open_boundary()) {
+          it->second.emplace_back(std::nullopt, vh);
+        } else {
+          it->second.emplace_back(approx_pt((above_vh)->point()), vh);
+        }
+      } else if(Halfedge_const_handle above_he; CGAL::assign(above_he, above_feat)) {
+        if(conn.is_connected((above_he)->source(), vh)) {
+          return;
+        }
+
+        const auto& [it, _] = feature_portals.try_emplace(above_he, std::vector<Portal>{});
+        if((above_he)->is_fictitious()) {
+          it->second.emplace_back(std::nullopt, vh);
+          return;
+        }
+        it->second.emplace_back(approx_pt(upper_intersection(*arr.traits(), vh->point(), (above_he)->curve())), vh);
+      } else if(Face_const_handle above_fh; CGAL::assign(above_fh, above_feat)) {
+        // We don't create portals for the unbounded face in bounded arrangements.
+        CGAL_assertion(above_fh->is_unbounded() && !above_fh->has_outer_ccb());
+      } else {
+        return;
+      }
+
+      visited_ccbs.insert(ccb_main_vertex);
+    });
+
+    decompose(arr, func_out_iter);
+    return feature_portals;
+  }
+};
+} // namespace CGAL
+#endif // CGAL_DRAW_AOS_ARR_CREATE_PORTALS_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_render_context.h

@@ -0,0 +1,218 @@
+#ifndef CGAL_DRAW_AOS_ARR_RENDER_CONTEXT_H
+#define CGAL_DRAW_AOS_ARR_RENDER_CONTEXT_H
+#include "CGAL/Arr_point_location_result.h"
+#include "CGAL/Arr_trapezoid_ric_point_location.h"
+#include "CGAL/Bbox_2.h"
+#include "CGAL/Draw_aos/Arr_approximate_point_2.h"
+#include "CGAL/Draw_aos/Arr_approximation_cache.h"
+#include "CGAL/Draw_aos/Arr_construct_curve_end.h"
+#include "CGAL/Draw_aos/Arr_construct_segments.h"
+#include <algorithm>
+#include <atomic>
+#include <chrono>
+#include <CGAL/Arrangement_2.h>
+#include <CGAL/Draw_aos/helpers.h>
+#include <cstdlib>
+#include <fstream>
+#include <iterator>
+#include <limits>
+#include <memory>
+
+namespace CGAL {
+
+class Arr_cancellable_context_mixin
+{
+  using Clock = std::chrono::steady_clock;
+  using Duration = Clock::duration;
+  using Time_point = std::chrono::time_point<Clock, Duration>;
+
+protected:
+  Arr_cancellable_context_mixin()
+      : m_start_time(Clock::now())
+      , m_done(std::make_shared<std::atomic<bool>>(false)) {}
+
+public:
+  Time_point start_time() const { return m_start_time; }
+
+  Time_point end_time() const { return m_end_time; }
+
+  Duration elapsed_time() const { return Clock::now() - m_start_time; }
+
+  bool is_cancelled() const { return m_done->load(); }
+
+  void cancel() {
+    m_done->store(true, std::memory_order_relaxed);
+    m_end_time = Clock::now();
+  }
+
+private:
+  Time_point m_start_time, m_end_time;
+  std::shared_ptr<std::atomic<bool>> m_done;
+};
+
+class Arr_bounds_context_mixin
+{
+  using Approx_point = Arr_approximation_geometry_traits::Approx_point;
+
+protected:
+  Arr_bounds_context_mixin(const Bbox_2& bbox)
+      : m_bbox(bbox) {}
+
+public:
+  double xmin() const { return m_bbox.xmin(); }
+  double xmax() const { return m_bbox.xmax(); }
+  double ymin() const { return m_bbox.ymin(); }
+  double ymax() const { return m_bbox.ymax(); }
+
+  const Bbox_2& bbox() const { return m_bbox; }
+
+  template <typename FT>
+  bool strictly_contains_x(FT x) const {
+    return xmin() < x && x <= xmax();
+  }
+
+  template <typename FT>
+  bool strictly_contains_y(FT y) const {
+    return ymin() < y && y <= ymax();
+  }
+
+  template <typename Point>
+  bool strictly_contains(const Point& pt) const {
+    return strictly_contains_x(pt.x()) && strictly_contains_y(pt.y());
+  }
+
+  template <typename FT>
+  bool contains_x(FT x) const {
+    return xmin() <= x && x <= xmax();
+  }
+
+  template <typename FT>
+  bool contains_y(FT y) const {
+    return ymin() <= y && y <= ymax();
+  }
+
+  template <typename Point>
+  bool contains(const Point& pt) const {
+    return contains_x(pt.x()) && contains_y(pt.y());
+  }
+
+  template <typename Point>
+  bool is_on_boundary(const Point& pt) const {
+    return (pt.x() == xmin() || pt.x() == xmax()) && contains_y(pt.y()) ||
+           (pt.y() == ymin() || pt.y() == ymax()) && contains_x(pt.x());
+  }
+
+private:
+  const Bbox_2 m_bbox;
+};
+
+class Arr_geom_traits_context_mixin
+{
+public:
+  Arr_geom_traits_context_mixin(const Geom_traits& _traits)
+      : traits(_traits)
+      , cst_curve_end(_traits)
+      , cst_horizontal_segment(_traits)
+      , cst_vertical_segment(_traits)
+      , intersect_2(_traits.intersect_2_object())
+      , compare_xy_2(_traits.compare_xy_2_object())
+      , is_vertical_2(_traits.is_vertical_2_object())
+      , approx_pt(_traits) {}
+
+  const Geom_traits& traits;
+  const Arr_construct_curve_end<Geom_traits> cst_curve_end;
+  const Arr_construct_vertical_segment cst_vertical_segment;
+  const Arr_construct_horizontal_segment cst_horizontal_segment;
+  const Geom_traits::Intersect_2 intersect_2;
+  const Geom_traits::Compare_xy_2 compare_xy_2;
+  const Geom_traits::Is_vertical_2 is_vertical_2;
+  const Arr_approximate_point_2<Geom_traits> approx_pt;
+};
+
+class Arr_render_context : public Arr_cancellable_context_mixin, public Arr_geom_traits_context_mixin
+{
+  using Point_location = Arr_trapezoid_ric_point_location<Arrangement>;
+
+public:
+  Arr_render_context(const Arrangement& arr, const Point_location& pl, double approx_error)
+      : Arr_cancellable_context_mixin()
+      , Arr_geom_traits_context_mixin(*arr.traits())
+      , arr(arr)
+      , point_location(pl)
+      , counter(std::make_shared<std::size_t>(0)) // TODO: remove this after debugging
+      , approx_error(approx_error) {}
+
+public:
+  const double approx_error;
+  const Arrangement& arr;
+  std::shared_ptr<std::size_t> counter; // TODO: remove this after debugging
+  const Point_location& point_location;
+};
+
+class Arr_bounded_render_context : public Arr_render_context, public Arr_bounds_context_mixin
+{
+  using Approx_point = Arr_approximation_geometry_traits::Approx_point;
+  using Point_2 = Geom_traits::Point_2;
+  constexpr static double ep_base = std::numeric_limits<double>::epsilon();
+
+public:
+  Arr_bounded_render_context(const Arr_render_context& ctx, const Bbox_2& bbox, Arr_approximation_cache& cache)
+      : Arr_render_context(ctx)
+      , ep_xmin(std::max(std::abs(ep_base * bbox.xmin()), ep_base))
+      , ep_xmax(std::max(std::abs(ep_base * bbox.xmax()), ep_base))
+      , ep_ymin(std::max(std::abs(ep_base * bbox.ymin()), ep_base))
+      , ep_ymax(std::max(std::abs(ep_base * bbox.ymax()), ep_base))
+      , Arr_bounds_context_mixin(bbox)
+      , cache(cache) {
+
+    // TODO: remove this after debugging
+    std::ofstream ofs_index("/Users/shep/codes/aos_2_js_helper/shapes.txt", std::ios::out | std::ios::trunc);
+  }
+
+  Approx_point approx_pt_on_boundary(const Point_2& pt) const {
+    double x = this->approx_pt(pt, 0);
+    double y = this->approx_pt(pt, 1);
+
+    if(std::abs(x - xmin()) < ep_xmin) {
+      x = xmin();
+    } else if(std::abs(x - xmax()) < ep_xmax) {
+      x = xmax();
+    } else if(std::abs(y - ymin()) < ep_ymin) {
+      y = ymin();
+    } else if(std::abs(y - ymax()) < ep_ymax) {
+      y = ymax();
+    } else {
+      // We shall not call this function if not approximated from a boundary point.
+      CGAL_assertion(false && "Failed to match point to the boundary");
+    }
+
+    return Approx_point(x, y);
+  }
+
+public:
+  Arr_approximation_cache& cache;
+
+private:
+  // floating point epsilon at boundary coordinates
+  double ep_xmin, ep_xmax, ep_ymin, ep_ymax;
+};
+
+template <typename Context>
+class Arr_context_delegator
+{
+public:
+  Arr_context_delegator(const Context& ctx)
+      : m_ctx(ctx) {}
+
+  const Context* operator->() const { return &m_ctx; }
+  const Context& operator*() const { return m_ctx; }
+
+  // implicit conversion to the context type
+  operator const Context&() const { return m_ctx; }
+
+private:
+  const Context& m_ctx;
+};
+
+} // namespace CGAL
+#endif // CGAL_DRAW_AOS_ARR_RENDER_CONTEXT_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/Arr_viewer.h

@@ -0,0 +1,234 @@
+
+#ifndef ARR_VIEWER_H
+#define ARR_VIEWER_H
+#include <CGAL/Draw_aos/helpers.h>
+#include <array>
+#include <boost/iterator/function_output_iterator.hpp>
+#include "CGAL/Arr_trapezoid_ric_point_location.h"
+#include "CGAL/Bbox_2.h"
+#include "CGAL/Draw_aos/Arr_bounded_renderer.h"
+#include "CGAL/Draw_aos/Arr_render_context.h"
+#include "CGAL/Graphics_scene.h"
+#include "CGAL/Graphics_scene_options.h"
+#include "CGAL/Qt/camera.h"
+#include <QtWidgets/QApplication>
+#include <QtWidgets/QWidget>
+#include <QtOpenGLWidgets/QtOpenGLWidgets>
+#include <QtGui/QOpenGLFunctions>
+#include <QtGui/QMouseEvent>
+#include <QtGui/QKeyEvent>
+#include <CGAL/Basic_viewer.h>
+#include <cstddef>
+#include <boost/range/iterator_range.hpp>
+
+namespace CGAL {
+
+class Arr_viewer : public Qt::Basic_viewer
+{
+  using Basic_viewer = Qt::Basic_viewer;
+  using Vertex_const_handle = Arrangement::Vertex_const_handle;
+  using Halfedge_const_handle = Arrangement::Halfedge_const_handle;
+  using Face_const_handle = Arrangement::Face_const_handle;
+  using Graphics_scene_options =
+      Graphics_scene_options<Arrangement, Vertex_const_handle, Halfedge_const_handle, Face_const_handle>;
+
+private:
+  // Function to check if the camera's state has changed
+  bool is_camera_changed() {
+    QMatrix4x4 proj_mat, mv_mat;
+    this->camera_->computeProjectionMatrix();
+    this->camera_->computeModelViewMatrix();
+    this->camera_->getProjectionMatrix(proj_mat.data());
+    this->camera_->getModelViewMatrix(mv_mat.data());
+    if(proj_mat == m_last_proj_matrix && mv_mat == m_last_modelview_matrix) {
+      return false;
+    }
+    m_last_proj_matrix = proj_mat;
+    m_last_modelview_matrix = mv_mat;
+    return true;
+  }
+
+  Bbox_2 initial_bbox() const {
+    Bbox_2 bbox;
+    for(const auto& vh : m_arr.vertex_handles()) {
+      bbox += vh->point().bbox();
+    }
+    if(bbox.x_span() == 0 || bbox.y_span() == 0) {
+      // make a default bbox around the degenrate rect
+      bbox = Bbox_2(bbox.xmin() - 1, bbox.ymin() - 1, bbox.xmax() + 1, bbox.ymax() + 1);
+    }
+    return bbox;
+  }
+
+  /**
+   * @brief Computes the bounding box of the view from orthogonal camera.
+   *
+   * Currently we reverse the entire projection to find out the view bounding box.
+   * @return Bbox_2
+   */
+  Bbox_2 view_bbox_from_camera() const {
+    QMatrix4x4 modelview_matrix, projection_matrix;
+    camera_->getModelViewMatrix(modelview_matrix.data());
+    camera_->getProjectionMatrix(projection_matrix.data());
+    QMatrix4x4 inverse_mvp = (projection_matrix * modelview_matrix).inverted();
+
+    // Define 4 corners of the near plane in NDC (-1 to 1 in x and y)
+    QVector4D clip_space_corners[] = {QVector4D(-1.0, -1.0, 0.0, 1.0), QVector4D(-1.0, 1.0, 0.0, 1.0),
+                                      QVector4D(1.0, -1.0, 0.0, 1.0), QVector4D(1.0, 1.0, 0.0, 1.0)};
+
+    double xmin = std::numeric_limits<double>::max();
+    double xmax = std::numeric_limits<double>::lowest();
+    double ymin = std::numeric_limits<double>::max();
+    double ymax = std::numeric_limits<double>::lowest();
+
+    for(const QVector4D& corner : clip_space_corners) {
+      QVector4D world = inverse_mvp * corner;
+      if(world.w() != 0.0)
+        world /= world.w();
+      double x = world.x();
+      double y = world.y();
+
+      xmin = std::min(xmin, x);
+      xmax = std::max(xmax, x);
+      ymin = std::min(ymin, y);
+      ymax = std::max(ymax, y);
+    }
+
+    return Bbox_2(xmin, ymin, xmax, ymax);
+  }
+
+  double get_approx_error(const Bbox_2& bbox) const {
+    std::array<GLint, 4> viewport;
+    camera_->getViewport(viewport.data());
+    double width = static_cast<double>(viewport[2]);
+    return 0.5;
+  }
+
+public:
+  Arr_viewer(QWidget* parent,
+             const Arrangement& arr,
+             Graphics_scene_options options,
+             const char* title = "Arrangement Viewer")
+      : Basic_viewer(parent, m_scene, title)
+      , m_scene_options(options)
+      , m_pl(arr)
+      , m_arr(arr) {}
+
+  void rerender(Bbox_2 bbox) {
+    Arr_render_context ctx(m_arr, m_pl, get_approx_error(bbox));
+    Arr_bounded_renderer renderer(ctx, bbox);
+    const auto& cache = renderer.render();
+
+    m_scene.clear();
+
+    int counter = 0;
+
+    // add faces
+    for(const auto& [fh, face_tris] : cache.face_cache()) {
+      const auto& points = face_tris.points;
+      const auto& tris = face_tris.triangles;
+      bool draw_face = m_scene_options.colored_face(m_arr, fh);
+      for(const auto& t : tris) {
+        if(draw_face) {
+          m_scene.face_begin(m_scene_options.face_color(m_arr, fh));
+        } else {
+          m_scene.face_begin();
+        }
+        for(const auto idx : t) {
+          m_scene.add_point_in_face(points[idx]);
+        }
+        m_scene.face_end();
+      }
+    }
+
+    // add edges
+    counter = 0;
+    for(const auto& [he, polyline] : cache.halfedge_cache()) {
+      if(polyline.size() < 2) {
+        continue; // skip degenerate edges
+      }
+      // ofs_index << "polyline_" << counter << ".txt" << std::endl;
+      // std::ofstream ofs("/Users/shep/codes/aos_2_js_helper/polyline_" + std::to_string(counter++) + ".txt");
+      bool draw_colored_edge = m_scene_options.colored_edge(m_arr, he);
+      auto color = draw_colored_edge ? m_scene_options.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(mid_pt.x() <= bbox.xmin() || mid_pt.x() > bbox.xmax() || mid_pt.y() <= bbox.ymin() ||
+           mid_pt.y() > bbox.ymax())
+        {
+          continue;
+        }
+        // ofs << cur_pt << "\n";
+        if(draw_colored_edge) {
+          m_scene.add_segment(cur_pt, next_pt, color);
+        } else {
+          m_scene.add_segment(cur_pt, next_pt);
+        }
+      }
+      // ofs << polyline.back() << "\n";
+      // ofs << std::endl;
+    }
+    // add vertices
+    for(const auto& [vh, pt] : cache.vertex_cache()) {
+      // std::cout << pt << std::endl;
+      if(m_scene_options.colored_vertex(m_arr, vh)) {
+        m_scene.add_point(pt, m_scene_options.vertex_color(m_arr, vh));
+      } else {
+        m_scene.add_point(pt);
+      }
+    }
+    Basic_viewer::redraw();
+  }
+
+  virtual void draw() override {
+    if(is_camera_changed()) {
+      Bbox_2 bbox = view_bbox_from_camera();
+      // shrink the bbox by 10% for testing
+      double dx = (bbox.xmax() - bbox.xmin()) * 0.1;
+      double dy = (bbox.ymax() - bbox.ymin()) * 0.1;
+      bbox = Bbox_2(bbox.xmin() + dx, bbox.ymin() + dy, bbox.xmax() - dx, bbox.ymax() - dy);
+      std::cout << "Camera changed, recomputing arrangement bounding box: " << bbox << std::endl;
+      rerender(bbox);
+    }
+    Basic_viewer::draw();
+  }
+
+  virtual ~Arr_viewer() {}
+
+private:
+  Graphics_scene m_scene;
+  Graphics_scene_options m_scene_options;
+  const Arrangement& m_arr;
+  Arr_trapezoid_ric_point_location<Arrangement> m_pl;
+  QMatrix4x4 m_last_proj_matrix;
+  QMatrix4x4 m_last_modelview_matrix;
+};
+
+void draw_viewer(const Arrangement& arr) {
+  Qt::init_ogl_context(4, 3);
+  int argc;
+  QApplication app(argc, nullptr);
+  Graphics_scene_options<Arrangement, Arrangement::Vertex_const_handle, Arrangement::Halfedge_const_handle,
+                         Arrangement::Face_const_handle>
+      gso;
+  gso.enable_faces();
+  gso.enable_edges();
+  gso.enable_vertices();
+  gso.face_color = [](const Arrangement&, const Arrangement::Face_const_handle& fh) -> CGAL::IO::Color {
+    CGAL::Random random((size_t(&*fh)));
+    return get_random_color(random);
+  };
+  gso.colored_face = [](const Arrangement&, const Arrangement::Face_const_handle&) { return true; };
+  gso.vertex_color = [](const Arrangement&, const Arrangement::Vertex_const_handle& vh) -> CGAL::IO::Color {
+    CGAL::Random random((size_t(&*vh)));
+    return get_random_color(random);
+  };
+  Arr_viewer viewer(app.activeWindow(), arr, gso, "Arrangement Viewer");
+  viewer.show();
+  app.exec();
+}
+
+} // namespace CGAL
+#endif // ARR_VIEWER_H

Arrangement_on_surface_2/include/CGAL/Draw_aos/helpers.h

@@ -0,0 +1,20 @@
+#ifndef CGAL_DRAW_AOS_HELPERS_H
+#define CGAL_DRAW_AOS_HELPERS_H
+#include "CGAL/Arr_linear_traits_2.h"
+#include <CGAL/Arr_segment_traits_2.h>
+#include <CGAL/Arrangement_2.h>
+namespace CGAL {
+
+using Exact_kernel = CGAL::Exact_predicates_exact_constructions_kernel;
+// using Geom_traits = Arr_linear_traits_2<Exact_kernel>;
+using Geom_traits = CGAL::Arr_linear_traits_2<Exact_kernel>;
+using Arrangement = Arrangement_2<Geom_traits>;
+
+struct Inner_ccb_tag
+{};
+struct Outer_ccb_tag
+{};
+
+} // namespace CGAL
+
+#endif // CGAL_DRAW_AOS_HELPERS_H

Arrangement_on_surface_2/include/CGAL/draw_arrangement_2.h

@@ -18,21 +18,31 @@
 
 #include <CGAL/license/Arrangement_on_surface_2.h>
 
+#include <CGAL/Draw_aos/helpers.h>
+#include <CGAL/Draw_aos/Arr_viewer.h>
+#include <CGAL/Draw_aos/Arr_bounded_renderer.h>
+#include <CGAL/Draw_aos/Arr_bounded_approximate_curve_2.h>
+#include <CGAL/Draw_aos/Arr_bounded_approximate_face_2.h>
+#include <CGAL/Draw_aos/Arr_bounded_approximate_point_2.h>
+#include <CGAL/Draw_aos/Arr_approximation_cache.h>
+#include <CGAL/Draw_aos/Arr_bounded_compute_y_at_x.h>
+#include <CGAL/Draw_aos/Arr_render_context.h>
+#include <CGAL/Draw_aos/Arr_bounded_face_triangulator.h>
+#include <CGAL/Draw_aos/Arr_approximation_geometry_traits.h>
 
 #include <CGAL/config.h>
 
-#include <unordered_map>
 #include <cstdlib>
-#include <random>
 #include <type_traits>
+#include <unordered_map>
 
+#include <CGAL/Arr_geodesic_arc_on_sphere_traits_2.h>
+#include <CGAL/Arrangement_2.h>
+#include <CGAL/Arrangement_on_surface_2.h>
 #include <CGAL/Basic_viewer.h>
 #include <CGAL/Graphics_scene.h>
 #include <CGAL/Graphics_scene_options.h>
 #include <CGAL/Random.h>
-#include <CGAL/Arrangement_on_surface_2.h>
-#include <CGAL/Arrangement_2.h>
-#include <CGAL/Arr_geodesic_arc_on_sphere_traits_2.h>
 
 namespace CGAL {
 
@@ -47,12 +57,12 @@ namespace draw_function_for_arrangement_2 {
 // Primary templates: detection fails by default
 // Does the traits have approximate_2_object()?
 template <typename, typename = std::void_t<>>
-struct has_approximate_2_object : std::false_type {};
+struct has_approximate_2_object : std::false_type
+{};
 
 // Specialization: detection succeeds if decltype(T::approximate_2_object()) is valid
 template <typename T>
-struct has_approximate_2_object<T, std::void_t<decltype(std::declval<T>().approximate_2_object())>> :
+struct has_approximate_2_object<T, std::void_t<decltype(std::declval<T>().approximate_2_object())>> : std::true_type
-    std::true_type
 {};
 
 // Convenience variable
@@ -64,50 +74,51 @@ inline constexpr bool has_approximate_2_object_v = has_approximate_2_object<T>::
 // Primary templates: detection fails by default
 // Does a class have operator()(const Point&)?
 template <typename, typename, typename = std::void_t<>>
-struct has_operator_point : std::false_type {};
+struct has_operator_point : std::false_type
+{};
 
 // Specialization: detection succeeds if decltype works out
-  template <typename T, typename A>
+template <typename T, typename A>
-  struct has_operator_point<T, A, std::void_t<decltype(std::declval<A>()(std::declval<const typename T::Point_2&>()))>> :
+struct has_operator_point<T, A, std::void_t<decltype(std::declval<A>()(std::declval<const typename T::Point_2&>()))>>
-    std::true_type
+    : std::true_type
 {};
 
 // Convenience variable
-  template <typename T, typename A>
+template <typename T, typename A>
-  inline constexpr bool has_operator_point_v = has_operator_point<T, A>::value;
+inline constexpr bool has_operator_point_v = has_operator_point<T, A>::value;
 
 // ========
 
 // Primary templates: detection fails by default
 // Does a class have operator()(const X_monotone_curve&)?
 template <typename, typename, typename = std::void_t<>>
-struct has_operator_xcv : std::false_type {};
+struct has_operator_xcv : std::false_type
+{};
 
 // Specialization: detection succeeds if decltype works out
-struct Dummy_output {};
+struct Dummy_output
+{};
 using Concrete_output_iterator = Dummy_output*;
 
 template <typename T, typename A>
-struct has_operator_xcv<T, A, std::void_t<decltype(
+struct has_operator_xcv<T,
-                                            std::declval<A>()(
+                        A,
-                                              std::declval<const typename T::X_monotone_curve_2&>(),
+                        std::void_t<decltype(std::declval<A>()(std::declval<const typename T::X_monotone_curve_2&>(),
-                                              std::declval<double>,
+                                                               std::declval<double>,
-                                              std::declval<Concrete_output_iterator>(),
+                                                               std::declval<Concrete_output_iterator>(),
-                                              std::declval<bool>()
+                                                               std::declval<bool>()))>> : std::true_type
-                                            )
-                                          )
-                      >> : std::true_type
 {};
 
 // Convenience variable
-  template <typename T, typename A>
+template <typename T, typename A>
-  inline constexpr bool has_operator_xcv_v = has_operator_xcv<T, A>::value;
+inline constexpr bool has_operator_xcv_v = has_operator_xcv<T, A>::value;
 
 // ========
 
 // Helper: detect whether T is or derives from Arr_geodesic_arc_on_sphere_traits_2<*, *, *>
 template <typename T>
-struct is_or_derived_from_agas {
+struct is_or_derived_from_agas
+{
 private:
   template <typename Kernel_, int AtanX, int AtanY>
   static std::true_type test(const Arr_geodesic_arc_on_sphere_traits_2<Kernel_, AtanX, AtanY>*);
@@ -124,8 +135,9 @@ inline constexpr bool is_or_derived_from_agas_v = is_or_derived_from_agas<T>::va
 // ========
 
 ///
-template<typename Arr, typename GSOptions>
+template <typename Arr, typename GSOptions>
-class Draw_arr_tool {
+class Draw_arr_tool
+{
 public:
   using Halfedge_const_handle = typename Arr::Halfedge_const_handle;
   using Vertex_const_handle = typename Arr::Vertex_const_handle;
@@ -139,20 +151,18 @@ public:
 
   /*! Construct
    */
-  Draw_arr_tool(Arr& a_aos, CGAL::Graphics_scene& a_gs, const GSOptions& a_gso) :
+  Draw_arr_tool(Arr& a_aos, CGAL::Graphics_scene& a_gs, const GSOptions& a_gso)
-    m_aos(a_aos), m_gs(a_gs),
+      : m_aos(a_aos)
-    m_gso(a_gso)
+      , m_gs(a_gs)
-  {}
+      , m_gso(a_gso) {}
 
   /// Add a face.
   void add_face(Face_const_handle face) {
     // std::cout << "add_face()\n";
-    for (Inner_ccb_const_iterator it = face->inner_ccbs_begin();
+    for(Inner_ccb_const_iterator it = face->inner_ccbs_begin(); it != face->inner_ccbs_end(); ++it)
-         it != face->inner_ccbs_end(); ++it)
+      add_ccb(*it);
-     add_ccb(*it);
 
-    for (Outer_ccb_const_iterator it = face->outer_ccbs_begin();
+    for(Outer_ccb_const_iterator it = face->outer_ccbs_begin(); it != face->outer_ccbs_end(); ++it) {
-         it != face->outer_ccbs_end(); ++it) {
       add_ccb(*it);
       draw_region(*it);
     }
@@ -164,10 +174,11 @@ public:
     auto curr = circ;
     do {
       auto new_face = curr->twin()->face();
-      if (m_visited.find(new_face) != m_visited.end()) continue;
+      if(m_visited.find(new_face) != m_visited.end())
+        continue;
       m_visited[new_face] = true;
       add_face(new_face);
-    } while (++curr != circ);
+    } while(++curr != circ);
   }
 
   /// Draw a region.
@@ -188,9 +199,10 @@ public:
      *
      * For now we use C++14 features.
      */
-    if (m_gso.colored_face(m_aos, circ->face()))
+    if(m_gso.colored_face(m_aos, circ->face()))
       m_gs.face_begin(m_gso.face_color(m_aos, circ->face()));
-    else m_gs.face_begin();
+    else
+      m_gs.face_begin();
 
     const auto* traits = this->m_aos.geometry_traits();
     auto ext = find_smallest(circ, *traits);
@@ -198,10 +210,11 @@ public:
 
     do {
       // Skip halfedges that are "antenas":
-      while (curr->face() == curr->twin()->face()) curr = curr->twin()->next();
+      while(curr->face() == curr->twin()->face())
+        curr = curr->twin()->next();
       draw_region_impl1(*traits, curr);
       curr = curr->next();
-    } while (curr != ext);
+    } while(curr != ext);
 
     m_gs.face_end();
   }
@@ -209,9 +222,10 @@ public:
   /// Compile time dispatching
 
   ///
-  template <typename T, typename A, std::enable_if_t<! has_operator_point_v<T, A>, int> = 0>
+  template <typename T, typename A, std::enable_if_t<!has_operator_point_v<T, A>, int> = 0>
-  void draw_region_impl2(const T& /* traits */, const A& /* approximate */, Halfedge_const_handle curr)
+  void draw_region_impl2(const T& /* traits */, const A& /* approximate */, Halfedge_const_handle curr) {
-  { draw_exact_region(curr); }
+    draw_exact_region(curr);
+  }
 
   ///
   template <typename T, typename A, std::enable_if_t<has_operator_point_v<T, A>, int> = 0>
@@ -221,12 +235,13 @@ public:
 
   /*! Draw a region, where the traits does not has approximate_2_object.
    */
-  template <typename T, std::enable_if_t<! has_approximate_2_object_v<T> && ! is_or_derived_from_agas_v<T>, int> = 0>
+  template <typename T, std::enable_if_t<!has_approximate_2_object_v<T> && !is_or_derived_from_agas_v<T>, int> = 0>
-  void draw_region_impl1(const T& /* traits */, Halfedge_const_handle curr)
+  void draw_region_impl1(const T& /* traits */, Halfedge_const_handle curr) {
-  { draw_exact_region(curr); }
+    draw_exact_region(curr);
+  }
 
   ///
-  template <typename T, std::enable_if_t<has_approximate_2_object_v<T> && ! is_or_derived_from_agas_v<T>, int> = 0>
+  template <typename T, std::enable_if_t<has_approximate_2_object_v<T> && !is_or_derived_from_agas_v<T>, int> = 0>
   auto draw_region_impl1(const T& traits, Halfedge_const_handle curr) {
     using Approximate = typename Gt::Approximate_2;
     draw_region_impl2(traits, traits.approximate_2_object(), curr);
@@ -253,35 +268,37 @@ public:
     double error(0.01); // TODO? (this->pixel_ratio());
     bool l2r = curr->direction() == ARR_LEFT_TO_RIGHT;
     approx(curr->curve(), error, std::back_inserter(polyline), l2r);
-    if (polyline.empty()) return;
+    if(polyline.empty())
+      return;
     auto it = polyline.begin();
     auto prev = it++;
-    for (; it != polyline.end(); prev = it++) m_gs.add_point_in_face(*prev);
+    for(; it != polyline.end(); prev = it++)
+      m_gs.add_point_in_face(*prev);
   }
 
   /*! Draw an exact curve.
    */
   template <typename XMonotoneCurve>
-  void draw_exact_curve(const XMonotoneCurve& curve,
+  void draw_exact_curve(const XMonotoneCurve& curve, bool colored, const CGAL::IO::Color& c) {
-                        bool colored, const CGAL::IO::Color& c) {
     const auto* traits = this->m_aos.geometry_traits();
     auto ctr_min = traits->construct_min_vertex_2_object();
     auto ctr_max = traits->construct_max_vertex_2_object();
     m_gs.add_segment(ctr_min(curve), ctr_max(curve));
-    if (colored) m_gs.add_segment(ctr_min(curve), ctr_max(curve), c);
+    if(colored)
-    else m_gs.add_segment(ctr_min(curve), ctr_max(curve));
+      m_gs.add_segment(ctr_min(curve), ctr_max(curve), c);
+    else
+      m_gs.add_segment(ctr_min(curve), ctr_max(curve));
   }
 
   /*! Draw a region in an exact manner.
    *  This fallback simply draws the curve in an exact manner (and even this is not guaranteed).
    */
-  void draw_exact_region(Halfedge_const_handle curr)
+  void draw_exact_region(Halfedge_const_handle curr) { draw_exact_curve(curr->curve(), false, CGAL::IO::Color()); }
-  { draw_exact_curve(curr->curve(), false, CGAL::IO::Color()); }
 
   /// Add all faces.
   template <typename Traits>
   void add_faces(const Traits&) {
-    for (auto it = m_aos.unbounded_faces_begin(); it != m_aos.unbounded_faces_end(); ++it)
+    for(auto it = m_aos.unbounded_faces_begin(); it != m_aos.unbounded_faces_end(); ++it)
       add_face(it);
   }
 
@@ -290,49 +307,52 @@ public:
   /*! Draw a point using approximate coordinates.
    */
   template <typename Approximate>
-  void draw_approximate_point(const Point& p, const Approximate& approx,
+  void draw_approximate_point(const Point& p, const Approximate& approx, bool colored, const CGAL::IO::Color& color) {
-                              bool colored, const CGAL::IO::Color& color) {
+    if(colored)
-    if (colored) m_gs.add_point(approx(p), color);
+      m_gs.add_point(approx(p), color);
-    else m_gs.add_point(approx(p));
+    else
+      m_gs.add_point(approx(p));
   }
 
   ///
   void draw_exact_point(const Point& p, bool colored, const CGAL::IO::Color& color) {
-    if (colored) m_gs.add_point(p, color);
+    if(colored)
-    else m_gs.add_point(p);
+      m_gs.add_point(p, color);
+    else
+      m_gs.add_point(p);
   }
 
   ///
-  template <typename T, typename A, std::enable_if_t<! has_operator_point_v<T, A>, int> = 0>
+  template <typename T, typename A, std::enable_if_t<!has_operator_point_v<T, A>, int> = 0>
-  void draw_point_impl2(const T& /* traits */, const A& /* approximate */, const Point& p,
+  void draw_point_impl2(
-                        bool colored, const CGAL::IO::Color& c)
+      const T& /* traits */, const A& /* approximate */, const Point& p, bool colored, const CGAL::IO::Color& c) {
-  { draw_exact_point(p, colored, c); }
+    draw_exact_point(p, colored, c);
+  }
 
   ///
   template <typename T, typename A, std::enable_if_t<has_operator_point_v<T, A>, int> = 0>
-  auto draw_point_impl2(const T& /* traits */, const A& approx, const Point& p,
+  auto
-                        bool colored, const CGAL::IO::Color& c) {
+  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);
   }
 
   /*! Draw a point, where the traits does not has approximate_2_object.
    */
-  template <typename T, std::enable_if_t<! has_approximate_2_object_v<T> && ! is_or_derived_from_agas_v<T>, int> = 0>
+  template <typename T, std::enable_if_t<!has_approximate_2_object_v<T> && !is_or_derived_from_agas_v<T>, int> = 0>
-  void draw_point_impl1(const T& /* traits */, const Point& p,
+  void draw_point_impl1(const T& /* traits */, const Point& p, bool colored, const CGAL::IO::Color& c) {
-                        bool colored, const CGAL::IO::Color& c)
+    draw_exact_point(p, colored, c);
-  { draw_exact_point(p, colored, c); }
+  }
 
   /*! Draw a point, where the traits does have approximate_2_object.
    */
-  template <typename T, std::enable_if_t<has_approximate_2_object_v<T> && ! is_or_derived_from_agas_v<T>, int> = 0>
+  template <typename T, std::enable_if_t<has_approximate_2_object_v<T> && !is_or_derived_from_agas_v<T>, int> = 0>
-  auto draw_point_impl1(const T& traits, const Point& p,
+  auto draw_point_impl1(const T& traits, const Point& p, bool colored, const CGAL::IO::Color& c) {
-                        bool colored, const CGAL::IO::Color& c) {
     draw_point_impl2(traits, traits.approximate_2_object(), p, colored, c);
   }
 
   /*! Draw a geodesic point.
    */
-    template <typename T, std::enable_if_t<is_or_derived_from_agas_v<T>, int> = 0>
+  template <typename T, std::enable_if_t<is_or_derived_from_agas_v<T>, int> = 0>
   void draw_point_impl1(const T& traits, const Point& p, bool colored, const CGAL::IO::Color& color) {
     using Traits = T;
     using Ak = typename Traits::Approximate_kernel;
@@ -343,10 +363,12 @@ public:
     auto x = ap.dx();
     auto y = ap.dy();
     auto z = ap.dz();
-    auto l = std::sqrt(x*x + y*y + z*z);
+    auto l = std::sqrt(x * x + y * y + z * z);
-    Approx_point_3 p3(x/l, y/l, z/l);
+    Approx_point_3 p3(x / l, y / l, z / l);
-    if (colored) m_gs.add_point(p3, color);
+    if(colored)
-    else m_gs.add_point(p3);
+      m_gs.add_point(p3, color);
+    else
+      m_gs.add_point(p3);
   }
 
   /// Draw a point.
@@ -357,10 +379,10 @@ public:
 
   ///
   template <typename Kernel, int AtanX, int AtanY>
-  Halfedge_const_handle
+  Halfedge_const_handle find_smallest(Ccb_halfedge_const_circulator circ,
-  find_smallest(Ccb_halfedge_const_circulator circ,
+                                      Arr_geodesic_arc_on_sphere_traits_2<Kernel, AtanX, AtanY> const&) {
-                Arr_geodesic_arc_on_sphere_traits_2<Kernel, AtanX, AtanY> const&)
+    return circ;
-  { return circ; }
+  }
 
   /*! Find the halfedge incident to the lexicographically smallest vertex
    * along the CCB, such that there is no other halfedge underneath.
@@ -374,32 +396,36 @@ public:
 
     // Find the first halfedge directed from left to right
     auto curr = circ;
-    do if (curr->direction() == CGAL::ARR_LEFT_TO_RIGHT) break;
+    do
-    while (++curr != circ);
+      if(curr->direction() == CGAL::ARR_LEFT_TO_RIGHT)
+        break;
+    while(++curr != circ);
     Halfedge_const_handle ext = curr;
 
     // Find the halfedge incident to the lexicographically smallest vertex,
     //  such that there is no other halfedge underneath.
     do {
       // Discard edges not directed from left to right:
-      if (curr->direction() != CGAL::ARR_LEFT_TO_RIGHT) continue;
+      if(curr->direction() != CGAL::ARR_LEFT_TO_RIGHT)
+        continue;
 
       auto res = cmp_xy(curr->source()->point(), ext->source()->point());
 
       // Discard the edges inciden to a point strictly larger than the point
       // incident to the stored extreme halfedge:
-      if (res == LARGER) continue;
+      if(res == LARGER)
+        continue;
 
       // Store the edge inciden to a point strictly smaller:
-      if (res == SMALLER) {
+      if(res == SMALLER) {
         ext = curr;
         continue;
       }
 
       // The incident points are equal; compare the halfedges themselves:
-      if (cmp_y(curr->curve(), ext->curve(), curr->source()->point()) == SMALLER)
+      if(cmp_y(curr->curve(), ext->curve(), curr->source()->point()) == SMALLER)
         ext = curr;
-    } while (++curr != circ);
+    } while(++curr != circ);
 
     return ext;
   }
@@ -410,29 +436,33 @@ public:
     // std::cout << "ratio: " << this->pixel_ratio() << std::endl;
     m_visited.clear();
 
-    if (m_aos.is_empty()) return;
+    if(m_aos.is_empty())
+      return;
 
-    if (m_gso.are_faces_enabled()) add_faces(*(this->m_aos.geometry_traits()));
+    if(m_gso.are_faces_enabled())
+      add_faces(*(this->m_aos.geometry_traits()));
 
     // Add edges that do not separate faces.
-    if (m_gso.are_edges_enabled()) {
+    if(m_gso.are_edges_enabled()) {
-      for (auto it = m_aos.edges_begin(); it != m_aos.edges_end(); ++it) {
+      for(auto it = m_aos.edges_begin(); it != m_aos.edges_end(); ++it) {
-        if (it->face() != it->twin()->face()) {
+        if(it->face() != it->twin()->face()) {
-          if (m_gso.draw_edge(m_aos, it)) {
+          if(m_gso.draw_edge(m_aos, it)) {
-            if (m_gso.colored_edge(m_aos, it))
+            if(m_gso.colored_edge(m_aos, it))
               draw_curve(it->curve(), true, m_gso.edge_color(m_aos, it));
-            else draw_curve(it->curve(), false, CGAL::IO::Color());
+            else
+              draw_curve(it->curve(), false, CGAL::IO::Color());
           }
         }
       }
     }
 
     // Add all points
-    if (m_gso.are_vertices_enabled()) {
+    if(m_gso.are_vertices_enabled()) {
-      for (auto it = m_aos.vertices_begin(); it != m_aos.vertices_end(); ++it) {
+      for(auto it = m_aos.vertices_begin(); it != m_aos.vertices_end(); ++it) {
-        if (m_gso.colored_vertex(m_aos, it))
+        if(m_gso.colored_vertex(m_aos, it))
           draw_point(it->point(), true, m_gso.vertex_color(m_aos, it));
-        else draw_point(it->point(), false, CGAL::IO::Color());
+        else
+          draw_point(it->point(), false, CGAL::IO::Color());
       }
     }
 
@@ -448,45 +478,52 @@ public:
   template <typename XMonotoneCurve, typename Approximate>
   void draw_approximate_curve(const XMonotoneCurve& curve,
                               const Approximate& approx,
-                              bool colored, const CGAL::IO::Color& c) {
+                              bool colored,
+                              const CGAL::IO::Color& c) {
     // std::cout << "draw_approximate_curve\n";
     std::vector<typename Gt::Approximate_point_2> polyline;
     double error(0.01); // TODO? (this->pixel_ratio());
     approx(curve, error, std::back_inserter(polyline));
-    if (polyline.empty()) return;
+    if(polyline.empty())
+      return;
     auto it = polyline.begin();
     auto prev = it++;
-    for (; it != polyline.end(); prev = it++) {
+    for(; it != polyline.end(); prev = it++) {
-      if (colored) m_gs.add_segment(*prev, *it, c);
+      if(colored)
-      else m_gs.add_segment(*prev, *it);
+        m_gs.add_segment(*prev, *it, c);
+      else
+        m_gs.add_segment(*prev, *it);
     }
   }
 
   ///
-  template <typename T, typename A, std::enable_if_t<! has_operator_point_v<T, A>, int> = 0>
+  template <typename T, typename A, std::enable_if_t<!has_operator_point_v<T, A>, int> = 0>
-  void draw_curve_impl2(const T& /* traits */, const A& /* approximate */, const X_monotone_curve& xcv,
+  void draw_curve_impl2(const T& /* traits */,
-                        bool colored, const CGAL::IO::Color& c)
+                        const A& /* approximate */,
-  { draw_exact_curve(xcv, colored, c); }
+                        const X_monotone_curve& xcv,
+                        bool colored,
+                        const CGAL::IO::Color& c) {
+    draw_exact_curve(xcv, colored, c);
+  }
 
   ///
   template <typename T, typename A, std::enable_if_t<has_operator_point_v<T, A>, int> = 0>
-  auto draw_curve_impl2(const T& /* traits */, const A& approx, const X_monotone_curve& xcv,
+  auto draw_curve_impl2(
-                        bool colored, const CGAL::IO::Color& c) {
+      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);
   }
 
   /*! Draw a curve, where the traits does not has approximate_2_object.
    */
-  template <typename T, std::enable_if_t<! has_approximate_2_object_v<T> && ! is_or_derived_from_agas_v<T>, int> = 0>
+  template <typename T, std::enable_if_t<!has_approximate_2_object_v<T> && !is_or_derived_from_agas_v<T>, int> = 0>
-  void draw_curve_impl1(const T& /* traits */, const X_monotone_curve& xcv,
+  void draw_curve_impl1(const T& /* traits */, const X_monotone_curve& xcv, bool colored, const CGAL::IO::Color& c) {
-                        bool colored, const CGAL::IO::Color& c)
+    draw_exact_curve(xcv, colored, c);
-  { draw_exact_curve(xcv, colored, c); }
+  }
 
   /*! Draw a curve, where the traits does have approximate_2_object.
    */
-  template <typename T, std::enable_if_t<has_approximate_2_object_v<T> && ! is_or_derived_from_agas_v<T>, int> = 0>
+  template <typename T, std::enable_if_t<has_approximate_2_object_v<T> && !is_or_derived_from_agas_v<T>, int> = 0>
-  auto draw_curve_impl1(const T& traits, const X_monotone_curve& xcv,
+  auto draw_curve_impl1(const T& traits, const X_monotone_curve& xcv, bool colored, const CGAL::IO::Color& c) {
-                        bool colored, const CGAL::IO::Color& c) {
     using Approximate = typename Gt::Approximate_2;
     draw_curve_impl2(traits, traits.approximate_2_object(), xcv, colored, c);
   }
@@ -494,8 +531,7 @@ public:
   /*! Draw a geodesic curve
    */
   template <typename T, std::enable_if_t<is_or_derived_from_agas_v<T>, int> = 0>
-  void draw_curve_impl1(const T& traits, const X_monotone_curve& xcv,
+  void draw_curve_impl1(const T& traits, const X_monotone_curve& xcv, bool colored, const CGAL::IO::Color& c) {
-                        bool colored, const CGAL::IO::Color& c) {
     // std::cout << "draw_curve (geodesic)\n";
     using Traits = T;
     using Kernel = typename Traits::Kernel;
@@ -511,16 +547,18 @@ public:
     auto x = it->dx();
     auto y = it->dy();
     auto z = it->dz();
-    auto l = std::sqrt(x*x + y*y + z*z);
+    auto l = std::sqrt(x * x + y * y + z * z);
-    Approx_point_3 prev(x/l, y/l, z/l);
+    Approx_point_3 prev(x / l, y / l, z / l);
-    for (++it; it != apoints.end(); ++it) {
+    for(++it; it != apoints.end(); ++it) {
       auto x = it->dx();
       auto y = it->dy();
       auto z = it->dz();
-      auto l = std::sqrt(x*x + y*y + z*z);
+      auto l = std::sqrt(x * x + y * y + z * z);
-      Approx_point_3 next(x/l, y/l, z/l);
+      Approx_point_3 next(x / l, y / l, z / l);
-      if (colored) m_gs.add_segment(prev, next, c);
+      if(colored)
-      else m_gs.add_segment(prev, next);
+        m_gs.add_segment(prev, next, c);
+      else
+        m_gs.add_segment(prev, next);
       prev = next;
     }
   }
@@ -570,31 +608,22 @@ protected:
 
 ///
 template <typename GeometryTraits_2, typename TopologyTraits, class GSOptions>
-void add_to_graphics_scene(const CGAL_ARR_TYPE& aos,
+void add_to_graphics_scene(const CGAL_ARR_TYPE& aos, CGAL::Graphics_scene& graphics_scene, const GSOptions& gso) {
-                           CGAL::Graphics_scene& graphics_scene,
-                           const GSOptions& gso) {
   draw_function_for_arrangement_2::Draw_arr_tool dar(aos, graphics_scene, gso);
   dar.add_elements();
 }
 
 ///
 template <typename GeometryTraits_2, typename TopologyTraits>
-void add_to_graphics_scene(const CGAL_ARR_TYPE& aos,
+void add_to_graphics_scene(const CGAL_ARR_TYPE& aos, CGAL::Graphics_scene& graphics_scene) {
-                           CGAL::Graphics_scene& graphics_scene) {
+  CGAL::Graphics_scene_options<CGAL_ARR_TYPE, typename CGAL_ARR_TYPE::Vertex_const_handle,
-  CGAL::Graphics_scene_options<CGAL_ARR_TYPE,
+                               typename CGAL_ARR_TYPE::Halfedge_const_handle, typename CGAL_ARR_TYPE::Face_const_handle>
-                               typename CGAL_ARR_TYPE::Vertex_const_handle,
+      gso;
-                               typename CGAL_ARR_TYPE::Halfedge_const_handle,
-                               typename CGAL_ARR_TYPE::Face_const_handle>
-    gso;
   // colored face?
-  gso.colored_face=[](const CGAL_ARR_TYPE&,
+  gso.colored_face = [](const CGAL_ARR_TYPE&, typename CGAL_ARR_TYPE::Face_const_handle) -> bool { return true; };
-                      typename CGAL_ARR_TYPE::Face_const_handle) -> bool
-  { return true; };
 
   // face color
-  gso.face_color=[](const CGAL_ARR_TYPE&,
+  gso.face_color = [](const CGAL_ARR_TYPE&, typename CGAL_ARR_TYPE::Face_const_handle fh) -> CGAL::IO::Color {
-                    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);
   };
@@ -604,18 +633,17 @@ void add_to_graphics_scene(const CGAL_ARR_TYPE& aos,
 
 /// Draw an arrangement on surface.
 template <typename GeometryTraits_2, typename TopologyTraits, class GSOptions>
-void draw(const CGAL_ARR_TYPE& aos, const GSOptions& gso,
+void draw(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 <typename GeometryTraits_2, typename TopologyTraits>
-void draw(const CGAL_ARR_TYPE& aos,
+void draw(const CGAL_ARR_TYPE& aos, const char* title = "2D Arrangement on Surface Basic Viewer") {
-          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);

Arrangement_on_surface_2/test/Arrangement_on_surface_2/CMakeLists.txt

@@ -6,6 +6,6 @@ project(Arrangement_on_surface_2_Tests)
 
 enable_testing()
 
-find_package(CGAL REQUIRED COMPONENTS Core)
+find_package(CGAL REQUIRED COMPONENTS Core Qt6)
 
 include(${CMAKE_CURRENT_SOURCE_DIR}/cgal_test.cmake)

Arrangement_on_surface_2/test/Arrangement_on_surface_2/approx_y_at_x.h

@@ -0,0 +1,11 @@
+#ifndef APPROX_Y_AT_X_H
+#define APPROX_Y_AT_X_H
+#include <CGAL/Arr_segment_traits_2.h>
+#include <CGAL/Arrangement_2.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+
+using Exact_kernel = CGAL::Exact_predicates_exact_constructions_kernel;
+using Segment_traits_2 = CGAL::Arr_segment_traits_2<Exact_kernel>;
+using Arr = CGAL::Arrangement_2<Segment_traits_2>;
+double approx_y_at_x(const Arr::Traits_2& traits, const Arr::X_monotone_curve_2& curve, double x) {}
+#endif // APPROX_Y_AT_X_H

Arrangement_on_surface_2/test/Arrangement_on_surface_2/cgal_test.cmake

@@ -1404,6 +1404,39 @@ function(test_algebraic_traits_core)
     VERTEX)
 endfunction()
 
+#---------------------------------------------------------------------#
+# Draw_aos tests
+#---------------------------------------------------------------------#
+function (test_draw_aos_arr_bounded_approximate_face_2)
+  set(nt ${CGAL_GMPQ_NT})
+  set(kernel ${CARTESIAN_KERNEL})
+  set(geom_traits ${GEODESIC_ARC_ON_SPHERE_GEOM_TRAITS})
+  set(topol_traits ${SPHERICAL_TOPOL_TRAITS})
+  set(flags "-DTEST_NT=${nt} -DTEST_KERNEL=${kernel} -DTEST_GEOM_TRAITS=${geom_traits} -DTEST_TOPOL_TRAITS=${topol_traits}")
+
+  compile_test_with_flags(test_draw_aos_arr_bounded_approximate_face_2 draw_aos "${flags}")
+  # target_link_libraries(test_draw_aos_arr_bounded_approximate_face_2)
+endfunction()
+
+function (test_drawing_planar)
+  set(nt ${CGAL_GMPQ_NT})
+  set(kernel ${CARTESIAN_KERNEL})
+  set(geom_traits ${POLYLINE_GEOM_TRAITS})
+  set(flags "-DCGAL_USE_BASIC_VIEWER")
+  compile_test_with_flags(test_drawing_planar drawing "${flags}")
+  # target_link_libraries(test_drawing_planar PRIVATE Qt6::Core Qt6::Gui Qt6::Widgets)
+  target_link_libraries(test_drawing_planar PRIVATE CGAL::CGAL_Basic_viewer CGAL::CGAL_Qt6)
+endfunction()
+
+function (test_drawing_spherical)
+  set(nt ${CGAL_GMPQ_NT})
+  set(kernel ${CARTESIAN_KERNEL})
+  set(geom_traits ${POLYLINE_GEOM_TRAITS})
+  set(flags "-DTEST_NT=${nt} -DTEST_KERNEL=${kernel} -DTEST_GEOM_TRAITS=${geom_traits}")
+
+  compile_test_with_flags(test_drawing_spherical drawing "${flags}")
+  target_link_libraries(test_drawing_spherical PRIVATE CGAL::CGAL_Basic_viewer CGAL::CGAL_Qt6)
+endfunction()
 
 configure("")
 compile_and_run(construction_test_suite_generator)
@@ -1486,6 +1519,11 @@ compile_and_run(test_io)
 compile_and_run(test_sgm)
 
 compile_and_run(test_polycurve_intersection)
+
+test_draw_aos_arr_bounded_approximate_face_2()
+test_drawing_planar()
+test_drawing_spherical()  
+
 if(CGAL_DISABLE_GMP)
   get_directory_property(LIST_OF_TESTS TESTS)
   foreach(_test ${LIST_OF_TESTS})

Arrangement_on_surface_2/test/Arrangement_on_surface_2/outer_ccb_0.txt

@@ -0,0 +1,4 @@
+0.200004 -1.24026
+3.73824 -1.24026
+3.73824 1.94415
+0.200004 1.94415

Arrangement_on_surface_2/test/Arrangement_on_surface_2/print_arr.h

@@ -0,0 +1,122 @@
+#ifndef _PRINT_ARR_H_
+#define _PRINT_ARR_H_
+#include <CGAL/Arrangement_2.h>
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <iostream>
+
+//-----------------------------------------------------------------------------
+// Print all neighboring vertices to a given arrangement vertex.
+//
+template <class Arrangement>
+void print_neighboring_vertices(typename Arrangement::Vertex_const_handle v) {
+  if(v->is_isolated()) {
+    std::cout << "The vertex (" << v->point() << ") is isolated\n";
+    return;
+  }
+
+  std::cout << "The neighbors of the vertex (" << v->point() << ") are:";
+  typename Arrangement::Halfedge_around_vertex_const_circulator first, curr;
+  first = curr = v->incident_halfedges();
+  do
+    std::cout << " (" << curr->source()->point() << ")";
+  while(++curr != first);
+  std::cout << std::endl;
+}
+
+//-----------------------------------------------------------------------------
+// Print all vertices (points) and edges (curves) along a connected component
+// boundary.
+//
+template <typename Arrangement>
+void print_ccb(typename Arrangement::Ccb_halfedge_const_circulator circ) {
+  std::cout << "(" << circ->source()->point() << ")";
+  typename Arrangement::Ccb_halfedge_const_circulator curr = circ;
+  do {
+    typename Arrangement::Halfedge_const_handle e = curr;
+    std::cout << "   [" << e->curve() << "]   "
+              << "(" << e->target()->point() << ")";
+  } while(++curr != circ);
+  std::cout << std::endl;
+}
+
+//-----------------------------------------------------------------------------
+// Print the boundary description of an arrangement face.
+//
+template <typename Arrangement>
+void print_face(typename Arrangement::Face_const_handle f) {
+  // Print the outer boundary.
+  if(f->is_unbounded()) {
+    std::cout << "Unbounded face.\n";
+  } else {
+    std::cout << "Outer boundary: ";
+    print_ccb<Arrangement>(f->outer_ccb());
+  }
+
+  // Print the boundary of each of the holes.
+  size_t index = 1;
+  for(auto hole = f->holes_begin(); hole != f->holes_end(); ++hole, ++index) {
+    std::cout << "    Hole #" << index << ": ";
+    // The following statement pacifies msvc.
+    typename Arrangement::Ccb_halfedge_const_circulator circ = *hole;
+    print_ccb<Arrangement>(circ);
+  }
+
+  // Print the isolated vertices.
+  index = 1;
+  for(auto iv = f->isolated_vertices_begin(); iv != f->isolated_vertices_end(); ++iv, ++index) {
+    std::cout << "    Isolated vertex #" << index << ": "
+              << "(" << iv->point() << ")" << std::endl;
+  }
+}
+
+//-----------------------------------------------------------------------------
+// Print the given arrangement.
+//
+template <typename Arrangement>
+void print_arrangement(const Arrangement& arr) {
+  CGAL_precondition(arr.is_valid());
+
+  // Print the arrangement vertices.
+  std::cout << arr.number_of_vertices() << " vertices:\n";
+  for(auto vit = arr.vertices_begin(); vit != arr.vertices_end(); ++vit) {
+    std::cout << "(" << vit->point() << ")";
+    if(vit->is_isolated())
+      std::cout << " - Isolated.\n";
+    else
+      std::cout << " - degree " << vit->degree() << std::endl;
+  }
+
+  // Print the arrangement edges.
+  std::cout << arr.number_of_edges() << " edges:\n";
+  for(auto eit = arr.edges_begin(); eit != arr.edges_end(); ++eit)
+    std::cout << "[" << eit->curve() << "]\n";
+
+  // Print the arrangement faces.
+  std::cout << arr.number_of_faces() << " faces:\n";
+  for(auto fit = arr.faces_begin(); fit != arr.faces_end(); ++fit)
+    print_face<Arrangement>(fit);
+}
+
+//-----------------------------------------------------------------------------
+// Print the size of the given arrangement.
+//
+template <typename Arrangement>
+void print_arrangement_size(const Arrangement& arr) {
+  std::cout << "The arrangement size:\n"
+            << "   |V| = " << arr.number_of_vertices() << ",  |E| = " << arr.number_of_edges()
+            << ",  |F| = " << arr.number_of_faces() << std::endl;
+}
+
+//-----------------------------------------------------------------------------
+// Print the size of the given unbounded arrangement.
+//
+template <typename Arrangement>
+void print_unbounded_arrangement_size(const Arrangement& arr) {
+  std::cout << "The arrangement size:\n"
+            << "   |V| = " << arr.number_of_vertices() << " (plus " << arr.number_of_vertices_at_infinity()
+            << " at infinity)"
+            << ",  |E| = " << arr.number_of_edges() << ",  |F| = " << arr.number_of_faces() << " ("
+            << arr.number_of_unbounded_faces() << " unbounded)\n\n";
+}
+
+#endif

Arrangement_on_surface_2/test/Arrangement_on_surface_2/test_draw_aos_arr_bounded_approximate_face_2.cpp

@@ -0,0 +1,194 @@
+#include "CGAL/Bbox_2.h"
+#include "CGAL/Draw_aos/Arr_approximation_geometry_traits.h"
+#include <CGAL/Draw_aos/Arr_bounded_approximate_face_2.h>
+#include <cstddef>
+#include <iterator>
+#include <vector>
+
+using Approx_point = CGAL::Arr_approximation_geometry_traits::Approx_point;
+
+struct Test_case
+{
+  std::vector<Approx_point> points_to_insert;
+  std::vector<Approx_point> expected_points;
+  CGAL::Bbox_2 bbox;
+
+  Test_case(std::vector<Approx_point> points, std::vector<Approx_point> expected, CGAL::Bbox_2 b)
+      : points_to_insert(std::move(points))
+      , expected_points(std::move(expected))
+      , bbox(b) {}
+};
+
+void test_geom_fill_corners_output_iterator(const Test_case& test_case) {
+  std::vector<Approx_point> points{};
+  auto inserter = std::back_inserter(points);
+  CGAL::Geom_fill_corners_output_iterator<std::back_insert_iterator<decltype(points)>> geom_inserter(
+      inserter, CGAL::Bbox_2(test_case.bbox));
+
+  for(const auto& pt : test_case.points_to_insert) {
+    geom_inserter++ = pt;
+  }
+
+  if(points != test_case.expected_points) {
+    std::cerr << "Test failed for bbox: " << test_case.bbox << std::endl;
+    std::cerr << "Points inserted: \n";
+    for(const auto& pt : test_case.points_to_insert) {
+      std::cerr << pt << "\n";
+    }
+    std::cerr << "\nExpected: \n";
+    for(const auto& pt : test_case.expected_points) {
+      std::cerr << pt << "\n";
+    }
+    std::cerr << "\nGot:\n";
+    for(const auto& pt : points) {
+      std::cerr << pt << "\n";
+    }
+    std::cerr << std::endl;
+    assert(false);
+  }
+}
+
+int main() {
+  CGAL::Bbox_2 bbox(0, 0, 10, 10);
+  std::vector<Test_case> cases{
+      /**
+       *  --------------------
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  1-------------------
+       */
+      Test_case({Approx_point(0, 0)}, {Approx_point(0, 0)}, bbox),
+
+      /**
+       *  2-------------------
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  1-------------------
+       *
+       *  Expected:
+       *
+       *  4------------------3
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  1------------------2
+       */
+      Test_case({Approx_point(0, 0), Approx_point(0, 10)},
+                {Approx_point(0, 0), Approx_point(10, 0), Approx_point(10, 10), Approx_point(0, 10)}, bbox),
+      /**
+       *  --------------------
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  1------------------2
+       *
+       *  Expected:
+       *
+       *  2-------------------
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  1------------------2
+       */
+      Test_case({Approx_point(0, 0), Approx_point(10, 0)}, {Approx_point(0, 0), Approx_point(10, 0)}, bbox),
+      /**
+       *  --------------------
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  1---------2---------
+       *
+       *  Expected:
+       *
+       *  --------------------
+       *  |                  |
+       *  |                  |
+       *  2                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  1-------------------
+       */
+      Test_case({Approx_point(0, 0), Approx_point(5, 0)}, {Approx_point(0, 0), Approx_point(5, 0)}, bbox),
+      /**
+       *  --------------------
+       *  |                  |
+       *  |                  |
+       *  |                  1
+       *  |                  |
+       *  |                  2
+       *  |                  |
+       *  --------------------
+       *
+       *  Expected:
+       *
+       *  3------------------2
+       *  |                  |
+       *  |                  |
+       *  |                  1
+       *  |                  |
+       *  |                  6
+       *  |                  |
+       *  4------------------5
+       */
+      Test_case({Approx_point(10, 5), Approx_point(10, 3)},
+                {Approx_point(10, 5), Approx_point(10, 10), Approx_point(0, 10), Approx_point(0, 0),
+                 Approx_point(10, 0), Approx_point(10, 3)},
+                bbox),
+      /**
+       *  -------------------1
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  2
+       *  |                  |
+       *  -------------------|
+       *
+       *  Expected:
+       *
+       *  2------------------1
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  |
+       *  |                  5
+       *  |                  |
+       *  3------------------4
+       */
+      Test_case({Approx_point(10, 10), Approx_point(10, 3)},
+                {
+                    Approx_point(10, 10),
+                    Approx_point(0, 10),
+                    Approx_point(0, 0),
+                    Approx_point(10, 0),
+                    Approx_point(10, 3),
+                },
+                bbox),
+  };
+
+  for(size_t i = 0; i < cases.size(); ++i) {
+    const auto& test_case = cases[i];
+    test_geom_fill_corners_output_iterator(test_case);
+  }
+}

Arrangement_on_surface_2/test/Arrangement_on_surface_2/test_drawing_planar.cpp

@@ -0,0 +1,322 @@
+
+#include "CGAL/Arr_algebraic_segment_traits_2.h"
+#include "CGAL/Arr_circle_segment_traits_2.h"
+#include "CGAL/Arr_conic_traits_2.h"
+#include "CGAL/Arr_default_dcel.h"
+#include "CGAL/Arr_enums.h"
+#include "CGAL/Arr_linear_traits_2.h"
+#include "CGAL/Arr_rational_function_traits_2.h"
+#include "CGAL/Arr_segment_traits_2.h"
+#include "CGAL/Arr_spherical_topology_traits_2.h"
+#include "CGAL/Arrangement_2.h"
+#include "CGAL/Arrangement_on_surface_2.h"
+#include "CGAL/CORE_algebraic_number_traits.h"
+#include "CGAL/Draw_aos/Arr_viewer.h"
+#include <array>
+#include <fstream>
+#include <iostream>
+#include "CGAL/Exact_predicates_exact_constructions_kernel.h"
+#include "CGAL/Exact_predicates_inexact_constructions_kernel.h"
+#include <CGAL/draw_arrangement_2.h>
+
+// void draw_segments_arr_1() {
+//   using Exact_kernel = CGAL::Exact_predicates_exact_constructions_kernel;
+//   using Segment_traits = CGAL::Arr_segment_traits_2<Exact_kernel>;
+//   using Point_2 = Segment_traits::Point_2;
+//   using Arrangement = CGAL::Arrangement_2<Segment_traits>;
+//   // Make a square
+//   Arrangement arr;
+//   auto traits = arr.traits();
+//   auto cst_x_curve = traits->construct_x_monotone_curve_2_object();
+//   auto square = {cst_x_curve({0, 0}, {5, 0}), cst_x_curve({5, 0}, {5, 5}), cst_x_curve({5, 5}, {0, 5}),
+//                  cst_x_curve({0, 5}, {0, 0})};
+//   insert(arr, square.begin(), square.end());
+//   auto hole_triangle = {
+//       cst_x_curve({1, 1}, {2, 1}),
+//       cst_x_curve({2, 1}, {2, 2}),
+//       cst_x_curve({2, 2}, {1, 1}),
+//   };
+//   insert(arr, hole_triangle.begin(), hole_triangle.end());
+//   // CGAL::draw_viewer(arr);
+// }
+
+// void draw_segments_arr_2() {
+//   using Exact_kernel = CGAL::Exact_predicates_exact_constructions_kernel;
+//   using Segment_traits = CGAL::Arr_segment_traits_2<Exact_kernel>;
+//   using Point_2 = Segment_traits::Point_2;
+//   using X_monotone_curve_2 = Segment_traits::X_monotone_curve_2;
+//   using Arrangement = CGAL::Arrangement_2<Segment_traits>;
+//   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<Point_2, 6> 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<X_monotone_curve_2, 6> 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_2{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_viewer(arr);
+// }
+
+// void draw_segments_arr_3() {
+//   // generate random segments and draw them
+//   using Exact_kernel = CGAL::Exact_predicates_exact_constructions_kernel;
+//   using Segment_traits = CGAL::Arr_segment_traits_2<Exact_kernel>;
+//   using Point_2 = Segment_traits::Point_2;
+//   using Arrangement = CGAL::Arrangement_2<Segment_traits>;
+//   using X_monotone_curve_2 = Segment_traits::X_monotone_curve_2;
+//   using Random = CGAL::Random;
+
+//   Arrangement arr;
+//   auto traits = arr.traits();
+//   auto cst_x_curve = traits->construct_x_monotone_curve_2_object();
+//   Random random;
+//   std::vector<X_monotone_curve_2> segments;
+//   std::cout << "Generating random segments..." << std::endl;
+//   for(int i = 0; i < 100; ++i) {
+//     // generate random points
+//     Point_2 p1(random.get_double(-100, 100), random.get_double(-100, 100));
+//     Point_2 p2(random.get_double(-100, 100), random.get_double(-100, 100));
+//     // create a segment
+//     X_monotone_curve_2 seg = cst_x_curve(p1, p2);
+//     segments.push_back(seg);
+//   }
+
+//   std::cout << "Inserting segments into the arrangement..." << std::endl;
+//   // insert segments into the arrangement
+//   CGAL::insert(arr, segments.begin(), segments.end());
+
+//   // draw the arrangement
+//   // CGAL::draw_viewer(arr);
+// }
+
+void draw_linear_arr_1() {
+  using Exact_kernel = CGAL::Exact_predicates_exact_constructions_kernel;
+  using Traits = CGAL::Arr_linear_traits_2<Exact_kernel>;
+  using Point_2 = Traits::Point_2;
+  using Line_2 = Traits::Line_2;
+  using Ray_2 = Traits::Ray_2;
+  using Curve_2 = Traits::Curve_2;
+  using Arrangement = CGAL::Arrangement_2<Traits>;
+  using Face_const_handle = Arrangement::Face_const_handle;
+  using Halfedge_const_handle = Arrangement::Halfedge_const_iterator;
+  using X_monotone_curve_2 = Traits::X_monotone_curve_2;
+
+  Arrangement arr;
+  auto x_axis = X_monotone_curve_2(Ray_2(Point_2(0, 0), Point_2(1, 0)));
+  auto y_axis = X_monotone_curve_2(Ray_2(Point_2(0, 0), Point_2(0, 1)));
+  CGAL::insert(arr, x_axis);
+  CGAL::insert(arr, y_axis);
+
+  CGAL::draw_viewer(arr);
+}
+
+void draw_linear_arr_2() {
+  using Exact_kernel = CGAL::Exact_predicates_exact_constructions_kernel;
+  using Traits = CGAL::Arr_linear_traits_2<Exact_kernel>;
+  using Point_2 = Traits::Point_2;
+  using Line_2 = Traits::Line_2;
+  using Segment_2 = Traits::Segment_2;
+  using Ray_2 = Traits::Ray_2;
+  using Curve_2 = Traits::Curve_2;
+  using Arrangement = CGAL::Arrangement_2<Traits>;
+  using Face_const_handle = Arrangement::Face_const_handle;
+  using Halfedge_const_handle = Arrangement::Halfedge_const_iterator;
+  using X_monotone_curve_2 = Traits::X_monotone_curve_2;
+
+  Arrangement arr;
+  auto& traits = *arr.traits();
+  // Insert a n*n grid, each cell is a square of size 5
+  int n = 1;
+  // for(int i = 0; i < n; ++i) {
+  //   Point_2 p1(i * 5, 0);
+  //   Point_2 p2(i * 5, 1);
+  //   CGAL::insert(arr, Curve_2(Line_2(p1, p2)));
+  // }
+  // for(int i = 0; i < n; ++i) {
+  //   Point_2 p1(0, i * 5);
+  //   Point_2 p2(1, i * 5);
+  //   CGAL::insert(arr, Curve_2(Line_2(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_2 p1(i * 5 + 1, j * 5 + 1);
+      Point_2 p2(i * 5 + 4, j * 5 + 4);
+      CGAL::insert(arr, Curve_2(Segment_2(p1, Point_2(p2.x(), p1.y()))));
+      CGAL::insert(arr, Curve_2(Segment_2(Point_2(p1.x(), p2.y()), p2)));
+      CGAL::insert(arr, Curve_2(Segment_2(p1, Point_2(p1.x(), p2.y()))));
+      CGAL::insert(arr, Curve_2(Segment_2(Point_2(p2.x(), p1.y()), p2)));
+
+      // Insert a triangle inside the square
+      Point_2 tri_p1(i * 5 + 2, j * 5 + 2);
+      Point_2 tri_p2(i * 5 + 3, j * 5 + 2);
+      Point_2 tri_p3(i * 5 + 2.5, j * 5 + 3);
+      CGAL::insert(arr, Curve_2(Segment_2(tri_p1, tri_p2)));
+      CGAL::insert(arr, Curve_2(Segment_2(tri_p2, tri_p3)));
+      CGAL::insert(arr, Curve_2(Segment_2(tri_p3, tri_p1)));
+
+      // Connect the triangle to the square
+      Point_2 top(i * 5 + 2.5, j * 5 + 4);
+      CGAL::insert(arr, Curve_2(Segment_2(tri_p1, top)));
+    }
+  }
+
+  auto arr2 = arr;
+
+  CGAL::draw_viewer(arr);
+}
+
+// // supports segments
+// void draw_circle_segs_arr() {
+//   using Exact_kernel = CGAL::Exact_predicates_exact_constructions_kernel;
+//   using Traits = CGAL::Arr_circle_segment_traits_2<Exact_kernel>;
+//   using Point_2 = Traits::Point_2;
+//   using Curve_2 = Traits::Curve_2;
+//   using Arrangement = CGAL::Arrangement_2<Traits>;
+
+//   auto traits = Traits();
+//   Arrangement arr;
+//   auto cv1 = Curve_2(Exact_kernel::Circle_2({0, 0}, 10));
+//   CGAL::insert(arr, cv1);
+//   CGAL::draw(arr);
+// }
+
+// void draw_conic_arcs_arr() {
+//   using Nt_traits = CGAL::CORE_algebraic_number_traits;
+//   using Rational = Nt_traits::Rational;
+//   using Rat_kernel = CGAL::Cartesian<Rational>;
+//   using Rat_point = Rat_kernel::Point_2;
+//   using Rat_segment = Rat_kernel::Segment_2;
+//   using Rat_circle = Rat_kernel::Circle_2;
+//   using Algebraic = Nt_traits::Algebraic;
+//   using Alg_kernel = CGAL::Cartesian<Algebraic>;
+//   using Traits = CGAL::Arr_conic_traits_2<Rat_kernel, Alg_kernel, Nt_traits>;
+//   using Point = Traits::Point_2;
+//   using Conic_arc = Traits::Curve_2;
+//   using X_monotone_conic_arc = Traits::X_monotone_curve_2;
+//   using Arrangement = CGAL::Arrangement_2<Traits>;
+
+//   Arrangement arr;
+//   auto traits = Traits();
+//   auto cst_x_curve = traits.construct_curve_2_object();
+
+//   auto vert_seg = cst_x_curve(Rat_segment(Rat_point(0, 0), Rat_point(1, 0)));
+//   auto hor_seg = cst_x_curve(Rat_segment(Rat_point(0, 0), Rat_point(0, 1)));
+
+//   CGAL::insert(arr, vert_seg);
+//   CGAL::insert(arr, hor_seg);
+
+//   CGAL::draw(arr);
+// }
+
+// void draw_algebraic_arr() {
+// #if CGAL_USE_GMP && CGAL_USE_MPFI
+// #include <CGAL/Gmpz.h>
+//   using Integer = CGAL::Gmpz;
+// #elif CGAL_USE_CORE
+// #include <CGAL/CORE_BigInt.h>
+//   using Integer = CORE::BigInt;
+// #else
+// #include <CGAL/leda_integer.h>
+//   using Integer = LEDA::integer;
+// #endif
+//   using Traits = CGAL::Arr_algebraic_segment_traits_2<Integer>;
+//   using Arrangement = CGAL::Arrangement_2<Traits>;
+//   using Polynomial = Traits::Polynomial_2;
+//   using X_monotone_curve_2 = Traits::X_monotone_curve_2;
+//   using Parameter_space_in_x_2 = Traits::Parameter_space_in_x_2;
+
+//   Arrangement arr;
+//   auto traits = arr.traits();
+//   X_monotone_curve_2 cv;
+//   auto param_space_in_x = traits->parameter_space_in_x_2_object();
+//   auto ctr_cv = traits->construct_curve_2_object();
+//   Polynomial x = CGAL::shift(Polynomial(1), 1, 0);
+//   Polynomial y = CGAL::shift(Polynomial(1), 1, 1);
+//   auto cst_x_curve = traits->construct_x_monotone_segment_2_object();
+//   auto curve = ctr_cv(CGAL::ipower(x, 4) + CGAL::ipower(y, 3) - 1);
+//   CGAL::insert(arr, curve);
+//   // CGAL::draw(arr);
+// }
+
+// void draw_rational_arr() {
+//   using AK1 = CGAL::Algebraic_kernel_d_1<CORE::BigInt>;
+//   using Traits = CGAL::Arr_rational_function_traits_2<AK1>;
+//   using Arrangement = CGAL::Arrangement_2<Traits>;
+//   using Polynomial = Traits::Polynomial_1;
+//   using Alg_real = Traits::Algebraic_real_1;
+//   using Bound = Traits::Bound;
+
+//   auto traits = Traits();
+//   auto approx = traits.approximate_2_object();
+//   auto cst_x_curve = traits.construct_x_monotone_curve_2_object();
+//   Arrangement arr;
+//   Polynomial x = CGAL::shift(Polynomial(1), 1);
+//   Polynomial P1 = CGAL::ipower(x, 4) - 6 * x * x + 8;
+//   Alg_real l(Bound(-2.1)), r(Bound(2.1));
+//   auto cv1 = cst_x_curve(P1, l, r);
+//   CGAL::insert(arr, cv1);
+//   // CGAL::draw(arr);
+// }
+
+// void draw_spherical_arr() {
+//   using Exact_kernel = CGAL::Exact_predicates_exact_constructions_kernel;
+//   using Direction_3 = Exact_kernel::Direction_3;
+//   using Geodesic_traits = CGAL::Arr_geodesic_arc_on_sphere_traits_2<Exact_kernel>;
+//   using Spherical_topo_traits = CGAL::Arr_spherical_topology_traits_2<Geodesic_traits>;
+//   using Arrangement = CGAL::Arrangement_on_surface_2<Geodesic_traits, Spherical_topo_traits>;
+//   using Point_2 = Geodesic_traits::Point_2;
+
+//   Arrangement arr;
+//   auto traits = arr.geometry_traits();
+//   auto cst_pt = traits->construct_point_2_object();
+//   auto cst_param = traits->parameter_space_in_x_2_object();
+
+//   Point_2 p1 = cst_pt(Direction_3(1, 0, 0));
+// }
+
+using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
+using Point = Kernel::Point_2;
+// void write_polyline(std::string filename, const std::vector<Point>& points) {
+//   std::ofstream ofs_index("/Users/shep/codes/aos_2_js_helper/shapes.txt");
+//   ofs_index << filename << std::endl;
+//   std::ofstream ofs("/Users/shep/codes/aos_2_js_helper/" + filename);
+//   for(const auto& pt : points) {
+//     ofs << pt << "\n";
+//   }
+//   ofs << std::endl;
+// }
+
+int main() {
+  // draw_segments_arr_2();
+  draw_linear_arr_2();
+  // test_zone();
+  // draw_conic_arcs_arr();
+  // draw_algebraic_arr();
+  // draw_rational_arr();
+  // draw_circle_segs_arr();
+  return 0;
+}

Arrangement_on_surface_2/test/Arrangement_on_surface_2/test_drawing_spherical.cpp

@@ -0,0 +1,99 @@
+#include "CGAL/Constrained_triangulation_2.h"
+#include "CGAL/Graphics_scene.h"
+#include "CGAL/Triangulation_data_structure_2.h"
+#include "CGAL/Triangulation_vertex_base_with_info_2.h"
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <algorithm>
+#include <cstddef>
+#include <utility>
+#include <vector>
+#include "CGAL/Random.h"
+#include "CGAL/IO/Color.h"
+#include "CGAL/draw_triangulation_2.h"
+#include "CGAL/mark_domain_in_triangulation.h"
+#include <CGAL/Graphics_scene_options.h>
+
+template <class PT>
+struct Polygon_triangulation_gs_options : public CGAL::Graphics_scene_options<typename PT::Triangulation,
+                                                                              typename PT::Vertex_handle,
+                                                                              typename PT::Finite_edges_iterator,
+                                                                              typename PT::Finite_faces_iterator>
+{
+  using T2 = typename PT::Triangulation;
+  template <class IPM>
+  Polygon_triangulation_gs_options(IPM ipm) {
+    this->colored_face = [](const T2&, const typename PT::Finite_faces_iterator) -> bool { return true; };
+
+    this->face_color = [ipm](const T2&, const typename PT::Finite_faces_iterator fh) -> CGAL::IO::Color {
+      if(!get(ipm, fh)) {
+        std::cout << "Face not in domain, returning black color." << std::endl;
+        return CGAL::IO::Color(0, 0, 0);
+      }
+      CGAL::Random random((unsigned int)(std::size_t)(&*fh));
+      return CGAL::get_random_color(random);
+    };
+
+    this->draw_face = [ipm](const T2&, const typename PT::Finite_faces_iterator fh) -> bool { return true; };
+
+    this->draw_edge = [ipm](const T2& pt, const typename PT::Finite_edges_iterator eh) -> bool {
+      typename PT::Face_handle fh1 = eh->first;
+      typename PT::Face_handle fh2 = pt.mirror_edge(*eh).first;
+      return get(ipm, fh1) || get(ipm, fh2);
+    };
+  }
+};
+
+int main() {
+  using K = CGAL::Exact_predicates_inexact_constructions_kernel;
+  using Vb = CGAL::Triangulation_vertex_base_with_info_2<size_t, K>;
+  using Fb = CGAL::Constrained_triangulation_face_base_2<K>;
+  using Tds = CGAL::Triangulation_data_structure_2<Vb, Fb>;
+  using Cdt = CGAL::Constrained_triangulation_2<K, Tds, CGAL::Exact_predicates_tag>;
+  using Point = K::Point_2;
+
+  Cdt cdt;
+  std::vector<Point> outer = {Point(0, 0), Point(5, 0), Point(5, 5), Point(0, 5), Point(0, 0)};
+  std::vector<Point> inner = {Point(0, 0), Point(0.5, 0), Point(0.5, 0.5), Point(0.5, 0),
+                              Point(1, 0), Point(1, 1),   Point(0, 1)};
+  std::vector<Point> outer_constraint = {Point(0, 0), Point(-1, -1), Point(5, 0), Point(6, -1),
+                                         Point(5, 5), Point(6, 6),   Point(0, 5), Point(-1, 6)};
+
+  auto add_info = [](const Point& p) { return std::make_pair(p, 1); };
+
+  cdt.insert_with_info<std::pair<Point, size_t>>(boost::make_transform_iterator(outer.begin(), add_info),
+                                                 boost::make_transform_iterator(outer.end(), add_info));
+  // cdt.insert_with_info<std::pair<Point, size_t>>(boost::make_transform_iterator(inner.begin(), add_info),
+  //                                                boost::make_transform_iterator(inner.end(), add_info));
+  cdt.insert_constraint(outer_constraint.begin(), outer_constraint.end(), true);
+  // cdt.insert_constraint(inner.begin(), inner.end(), true);
+  using In_domain_map = CGAL::unordered_flat_map<typename Cdt::Face_handle, bool>;
+  In_domain_map in_domain_map;
+  boost::associative_property_map<In_domain_map> in_domain(in_domain_map);
+
+  CGAL::mark_domain_in_triangulation(cdt, in_domain);
+
+  std::cout << "Number of faces in triangulation: " << cdt.number_of_faces() << std::endl;
+
+  CGAL::Graphics_scene_options<Cdt::Triangulation, Cdt::Vertex_handle, Cdt::Finite_edges_iterator,
+                               Cdt::Finite_faces_iterator>
+      gso;
+  gso.face_color = [&](const Cdt::Triangulation&, const typename Cdt::Finite_faces_iterator fh) -> CGAL::IO::Color {
+    if(!in_domain_map[fh]) {
+      return CGAL::IO::Color(255, 255, 255); // black for faces not in domain
+    }
+    std::array<std::size_t, 3> vertices;
+    for(int i = 0; i < 3; ++i) {
+      vertices[i] = fh->vertex(i)->info();
+    }
+    if(std::any_of(vertices.begin(), vertices.end(), [](auto idx) { return idx == 0; })) {
+      return CGAL::IO::Color(255, 255, 255);
+    }
+    CGAL::Random rand((std::size_t)(&*fh));
+    return CGAL::get_random_color(rand);
+  };
+  gso.colored_face = [&](const Cdt::Triangulation&, const typename Cdt::Finite_faces_iterator) -> bool {
+    return true; // always color faces
+  };
+
+  CGAL::draw(cdt, gso, "Polygon Triangulation Viewer");
+}

Arrangement_on_surface_2/test/Arrangement_on_surface_2/zone_exp.h

@@ -0,0 +1,71 @@
+#ifndef ZONE_EXP_H
+#define ZONE_EXP_H
+#include "CGAL/Arr_trapezoid_ric_point_location.h"
+#include <CGAL/Arr_algebraic_segment_traits_2.h>
+#include <CGAL/Arr_linear_traits_2.h>
+#include <CGAL/Arr_segment_traits_2.h>
+#include <CGAL/Arrangement_2.h>
+#include <CGAL/Arrangement_2/Arr_compute_zone_visitor.h>
+#include <CGAL/Arrangement_zone_2.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+#include <CGAL/draw_arrangement_2.h>
+
+void test_zone() {
+  using Exact_kernel = CGAL::Exact_predicates_exact_constructions_kernel;
+  using Segment_traits = CGAL::Arr_segment_traits_2<Exact_kernel>;
+  using Point_2 = Segment_traits::Point_2;
+  using Arrangement = CGAL::Arrangement_2<Segment_traits>;
+  using X_monotone_curve_2 = Segment_traits::X_monotone_curve_2;
+  using Traits = CGAL::Arr_linear_traits_2<Exact_kernel>;
+  using Face_const_handle = Arrangement::Face_const_handle;
+  using Halfedge_const_handle = Arrangement::Halfedge_const_handle;
+  using Vertex_const_handle = Arrangement::Vertex_const_handle;
+  using Face_handle = Arrangement::Face_handle;
+  using Halfedge_handle = Arrangement::Halfedge_handle;
+  using Vertex_handle = Arrangement::Vertex_handle;
+  using Halfedge_const_handle = Arrangement::Halfedge_const_handle;
+
+  Arrangement arr;
+  auto traits = arr.traits();
+
+  auto cst_x_curve = traits->construct_x_monotone_curve_2_object();
+  auto tri = {
+      cst_x_curve({-1, 0}, {1, 0}),
+      cst_x_curve({1, 0}, {0, 3}),
+      cst_x_curve({0, 3}, {-1, 0}),
+  };
+  CGAL::insert(arr, tri.begin(), tri.end());
+
+  // point location
+  auto pl = CGAL::Arr_trapezoid_ric_point_location<Arrangement>(arr);
+
+  // test zone construction
+  using Feature = std::variant<Vertex_handle, Halfedge_handle, Face_handle>;
+  using Feature_vector = std::vector<Feature>;
+  using Feature_vector_out_iter = std::back_insert_iterator<Feature_vector>;
+  using Zone_visitor = CGAL::Arr_compute_zone_visitor<Arrangement, Feature_vector_out_iter>;
+
+  Feature_vector zone_features;
+  Feature_vector_out_iter out_iter(zone_features);
+  Zone_visitor zone_visitor(out_iter);
+  auto cv = cst_x_curve({-1, 0}, {1, 0});
+  CGAL::Arrangement_zone_2<Arrangement, Zone_visitor> zone(arr, &zone_visitor);
+
+  zone.init(cv, pl);
+  zone.compute_zone();
+
+  for(const auto& feature : zone_features) {
+    if(auto* hf = std::get_if<Halfedge_handle>(&feature)) {
+      std::cout << "Halfedge: " << (*hf)->curve() << std::endl;
+      Halfedge_const_handle che(*hf);
+      std::cout << "  const curve" << che->curve() << std::endl;
+    } else if(auto* vh = std::get_if<Vertex_handle>(&feature)) {
+      std::cout << "Vertex: " << (*vh)->point() << std::endl;
+    } else if(auto* fh = std::get_if<Face_handle>(&feature)) {
+      std::cout << "Face: " << std::endl;
+    }
+  }
+
+  CGAL::draw(arr);
+}
+#endif // ZONE_EXP_H

Basic_viewer/include/CGAL/Basic_viewer.h

@@ -13,8 +13,8 @@
 #ifndef CGAL_BASIC_VIEWER_H
 #define CGAL_BASIC_VIEWER_H
 
-#include <CGAL/license/GraphicsView.h>
 #include <CGAL/Graphics_scene.h>
+#include <CGAL/license/GraphicsView.h>
 
 // compatibility
 #if defined(CGAL_USE_BASIC_VIEWER) && !defined(CGAL_USE_BASIC_VIEWER_QT)
@@ -30,14 +30,10 @@
 // #elif defined(CGAL_USE_BASIC_VIEWER_GLFW)
 // #include <CGAL/GLFW/Basic_viewer.h>
 #else
-namespace CGAL
+namespace CGAL {
-{
+inline void draw_graphics_scene(const Graphics_scene&, const char* = "CGAL Basic Viewer") {
-  inline
+  std::cerr << "Impossible to draw, CGAL_USE_BASIC_VIEWER is not defined." << std::endl;
-  void draw_graphics_scene(const Graphics_scene&,
+}
-                           const char* ="CGAL Basic Viewer")
-  {
-    std::cerr<<"Impossible to draw, CGAL_USE_BASIC_VIEWER is not defined."<<std::endl;
-  }
 } // End namespace CGAL
 #endif