Merge branch 'visibility2D' into 'master'

[feature] add Visibility methods based on CGAL Visibility_2 package

See merge request Oslandia/SFCGAL!287

.cirrus.yml

@@ -36,7 +36,7 @@ task:
     - brew install cmake boost cgal gmp mpfr
     - cmake -DSFCGAL_BUILD_TESTS=ON -S . -B build
     - cmake --build build
-    - ctest --test-dir build
+    - ctest -VV --test-dir build
   matrix:
     macos_instance:
       image: ghcr.io/cirruslabs/macos-ventura-base:latest

NEWS

@@ -2,6 +2,7 @@
     * Add polygon partition (Raphaël Delhome, Loïc Bartoletti)
     * WKT: Fix triangle code (Loïc Bartoletti)
     * Straight Skeleton: Add a version with extrusion
+    * Add visibility algorithms (Loïc Bartoletti)
 1.4.1 (2022-01-27):
     * Add alpha-shapes algorithm (Loïc Bartoletti, Hugo Mercier)
     * Fix build and tests for MSYS2/MinGW (Loïc Bartoletti)

src/LineString.h

@@ -186,6 +186,15 @@ class SFCGAL_API LineString : public Geometry {
   bool
   isClosed() const;
 
+  /**
+   * closes the LineString
+   */
+  inline void
+  closes()
+  {
+    _points.push_back(_points.front().clone());
+  }
+
   //-- iterators
 
   inline iterator

src/algorithm/visibility.cpp

@@ -0,0 +1,195 @@
+// Copyright (c) 2023-2023, Oslandia.
+// SPDX-License-Identifier: LGPL-2.0-or-later
+#include <SFCGAL/Polygon.h>
+#include <SFCGAL/algorithm/visibility.h>
+
+#include <SFCGAL/algorithm/isValid.h>
+
+#include <SFCGAL/Exception.h>
+
+#include <CGAL/Arr_segment_traits_2.h>
+#include <CGAL/Arrangement_2.h>
+#include <CGAL/Triangular_expansion_visibility_2.h>
+
+#include <CGAL/Arr_naive_point_location.h>
+#include <memory>
+
+namespace SFCGAL {
+namespace algorithm {
+
+using Point_2               = Kernel::Point_2;
+using Polygon_2             = CGAL::Polygon_2<Kernel>;
+using Segment_2             = Kernel::Segment_2;
+using Traits_2              = CGAL::Arr_segment_traits_2<Kernel>;
+using Arrangement_2         = CGAL::Arrangement_2<Traits_2>;
+using Face_handle           = Arrangement_2::Face_handle;
+using Halfedge_const_handle = Arrangement_2::Halfedge_const_handle;
+using TEV =
+    CGAL::Triangular_expansion_visibility_2<Arrangement_2, CGAL::Tag_true>;
+using PolygonWithHoles = CGAL::Polygon_with_holes_2<Kernel>;
+
+static auto
+query_visibility(Face_handle fh, Halfedge_const_handle he)
+    -> std::unique_ptr<Polygon>
+{
+
+  std::unique_ptr<LineString> extRing{new LineString()};
+  // Make sure the visibility polygon we find has an outer boundary
+  if (fh->has_outer_ccb()) {
+    Arrangement_2::Ccb_halfedge_circulator curr = fh->outer_ccb();
+
+    // find the right halfedge first
+    if (he != Halfedge_const_handle())
+      while (++curr != fh->outer_ccb())
+        if (curr->source()->point() == he->source()->point())
+          break;
+
+    Arrangement_2::Ccb_halfedge_circulator first = curr;
+    extRing->addPoint(Point(curr->source()->point()));
+
+    // Save the points from the visibility polygon
+    while (++curr != first) {
+      extRing->addPoint(Point(curr->source()->point()));
+    }
+  }
+
+  extRing->closes();
+  std::unique_ptr<Polygon> result{new Polygon(extRing.release())};
+
+  return result;
+}
+
+///
+///
+///
+auto
+visibility(const Geometry &polygon, const Geometry &point)
+    -> std::unique_ptr<Polygon>
+{
+  SFCGAL_ASSERT_GEOMETRY_VALIDITY_2D(polygon);
+
+  std::unique_ptr<Polygon> result(
+      visibility(polygon, point, NoValidityCheck()));
+  propagateValidityFlag(*result, true);
+  return result;
+}
+
+auto
+visibility(const Geometry &polygon, const Geometry &point, NoValidityCheck)
+    -> std::unique_ptr<Polygon>
+{
+
+  Point_2                     queryPoint{point.as<Point>().toPoint_2()};
+  std::unique_ptr<LineString> extRing{new LineString()};
+
+  // insert geometry into the arrangement
+  CGAL::Polygon_with_holes_2 pwh{
+      polygon.as<Polygon>().toPolygon_with_holes_2()};
+  Arrangement_2 arr;
+
+  CGAL::insert(arr, pwh.outer_boundary().edges_begin(),
+               pwh.outer_boundary().edges_end());
+  // the holes
+  for (PolygonWithHoles::Hole_const_iterator hit = pwh.holes_begin();
+       hit != pwh.holes_end(); ++hit)
+    CGAL::insert(arr, hit->edges_begin(), hit->edges_end());
+
+  // Find the face
+  Arrangement_2::Face_const_handle                    *face;
+  CGAL::Arr_naive_point_location<Arrangement_2>        pl(arr);
+  CGAL::Arr_point_location_result<Arrangement_2>::Type obj =
+      pl.locate(queryPoint);
+
+  // The query point locates in the interior of a face
+  face = boost::get<Arrangement_2::Face_const_handle>(&obj);
+  Arrangement_2         output_arr;
+  Face_handle           fh;
+  Halfedge_const_handle he = Halfedge_const_handle();
+
+  // Create Triangular Expansion Visibility object.
+  TEV tev(arr);
+
+  // If the point is within a face, we can compute the visibility that way
+  if (face != nullptr) {
+    fh = tev.compute_visibility(queryPoint, *face, output_arr);
+  } else {
+    // If the point is in a boundary segment, find the corresponding half edge
+    he        = arr.halfedges_begin();
+    bool cont = !Segment_2(he->source()->point(), he->target()->point())
+                     .has_on(queryPoint) ||
+                he->source()->point() == queryPoint ||
+                he->face()->is_unbounded();
+    // While we are not in the right half edge, or while q is the source,
+    // continue
+    while (cont) {
+      he++;
+      if (he == arr.halfedges_end()) {
+        BOOST_THROW_EXCEPTION(
+            Exception("Can not find corresponding half edge."));
+      }
+
+      cont = !Segment_2(he->source()->point(), he->target()->point())
+                  .has_on(queryPoint) ||
+             he->source()->point() == queryPoint || he->face()->is_unbounded();
+    }
+
+    // Use the half edge to compute the visibility
+    fh = tev.compute_visibility(queryPoint, he, output_arr);
+  }
+
+  return query_visibility(fh, he);
+}
+///
+///
+///
+auto
+visibility(const Geometry &polygon, const Geometry &pointA,
+           const Geometry &pointB) -> std::unique_ptr<Polygon>
+{
+  SFCGAL_ASSERT_GEOMETRY_VALIDITY_2D(polygon);
+
+  std::unique_ptr<Polygon> result(
+      visibility(polygon, pointA, pointB, NoValidityCheck()));
+  propagateValidityFlag(*result, true);
+  return result;
+}
+
+auto
+visibility(const Geometry &polygon, const Geometry &pointA,
+           const Geometry &pointB, NoValidityCheck) -> std::unique_ptr<Polygon>
+{
+
+  Point_2 startPoint{pointA.as<Point>().toPoint_2()};
+  Point_2 endPoint{pointB.as<Point>().toPoint_2()};
+
+  // insert geometry into the arrangement
+  CGAL::Polygon_with_holes_2 pwh{
+      polygon.as<Polygon>().toPolygon_with_holes_2()};
+  Arrangement_2 arr;
+
+  CGAL::insert(arr, pwh.outer_boundary().edges_begin(),
+               pwh.outer_boundary().edges_end());
+  for (PolygonWithHoles::Hole_const_iterator hit = pwh.holes_begin();
+       hit != pwh.holes_end(); ++hit)
+    CGAL::insert(arr, hit->edges_begin(), hit->edges_end());
+
+  // If the point is in a boundary segment, find the corresponding half edge
+  Halfedge_const_handle he = arr.halfedges_begin();
+  while (he->source()->point() != startPoint ||
+         he->target()->point() != endPoint) {
+    he++;
+    if (he == arr.halfedges_end()) {
+      BOOST_THROW_EXCEPTION(Exception("Can not find corresponding half edge."));
+    }
+  }
+
+  // visibility query
+  Arrangement_2 output_arr;
+  TEV           tev(arr);
+  Face_handle   fh = tev.compute_visibility(endPoint, he, output_arr);
+
+  return query_visibility(fh, he);
+}
+
+} // namespace algorithm
+} // namespace SFCGAL

src/algorithm/visibility.h

@@ -0,0 +1,78 @@
+// Copyright (c) 2023-2023, Oslandia.
+// SPDX-License-Identifier: LGPL-2.0-or-later
+
+#ifndef _SFCGAL_ALGORITHM_VISIBILITY2D_H_
+#define _SFCGAL_ALGORITHM_VISIBILITY2D_H_
+
+#include <SFCGAL/config.h>
+
+#include <memory>
+
+namespace SFCGAL {
+class Geometry;
+class Polygon;
+} // namespace SFCGAL
+
+namespace SFCGAL {
+namespace algorithm {
+struct NoValidityCheck;
+
+/**
+ * @brief build the visibility polygon of a Point inside a Polygon
+ * @param polygon input geometry
+ * @param point input geometry
+ * @ingroup public_api
+ * @pre polygon is a valid geometry
+ * @pre point must be inside polygon or on the boundary
+ */
+SFCGAL_API auto
+visibility(const Geometry &polygon, const Geometry &point)
+    -> std::unique_ptr<Polygon>;
+
+/**
+ * @brief build the visibility polygon of a Point inside a Polygon
+ * @param polygon input geometry
+ * @param point input geometry
+ * @ingroup public_api
+ * @pre polygon is a valid geometry
+ * @pre point must be inside polygon or on the boundary
+ * @warning No actual validity check is done
+ */
+SFCGAL_API auto
+visibility(const Geometry &polygon, const Geometry &point, NoValidityCheck)
+    -> std::unique_ptr<Polygon>;
+
+/**
+ * @brief build the visibility polygon of the segment [pointA ; pointB] on a
+ * Polygon
+ * @param polygon input geometry
+ * @param pointA input geometry
+ * @param pointB input geometry
+ * @ingroup public_api
+ * @pre polygon is a valid geometry
+ * @pre pointA and pointB must be vertices of poly, adjacents and respect the
+ * direction
+ */
+SFCGAL_API auto
+visibility(const Geometry &polygon, const Geometry &pointA,
+           const Geometry &pointB) -> std::unique_ptr<Polygon>;
+
+/**
+ * @brief build the visibility polygon of a Point inside a Polygon
+ * @param polygon input geometry
+ * @param pointA input geometry
+ * @param pointB input geometry
+ * @ingroup public_api
+ * @pre polygon is a valid geometry
+ * @warning No actual validity check is done
+ * @pre pointA and pointB must be vertices of poly, adjacents and respect the
+ * direction
+ */
+SFCGAL_API auto
+visibility(const Geometry &polygon, const Geometry &pointA,
+           const Geometry &pointB, NoValidityCheck) -> std::unique_ptr<Polygon>;
+
+} // namespace algorithm
+} // namespace SFCGAL
+
+#endif

src/capi/sfcgal_c.cpp

@@ -42,6 +42,7 @@
 #include <SFCGAL/algorithm/straightSkeleton.h>
 #include <SFCGAL/algorithm/tesselate.h>
 #include <SFCGAL/algorithm/union.h>
+#include <SFCGAL/algorithm/visibility.h>
 #include <SFCGAL/algorithm/volume.h>
 #include <SFCGAL/triangulate/triangulate2DZ.h>
 
@@ -1381,3 +1382,73 @@ sfcgal_optimal_convex_partition_2(const sfcgal_geometry_t *geom)
 
   return result.release();
 }
+
+extern "C" sfcgal_geometry_t *
+sfcgal_geometry_visibility_point(const sfcgal_geometry_t *polygon,
+                                 const sfcgal_geometry_t *point)
+{
+
+  const auto *poly = reinterpret_cast<const SFCGAL::Geometry *>(polygon);
+  const auto *pt   = reinterpret_cast<const SFCGAL::Geometry *>(point);
+  std::unique_ptr<SFCGAL::Geometry> result;
+
+  if (poly->geometryTypeId() != SFCGAL::TYPE_POLYGON) {
+    SFCGAL_ERROR("visibility() only applies to polygons");
+    return result.release();
+  }
+
+  if (pt->geometryTypeId() != SFCGAL::TYPE_POINT) {
+    SFCGAL_ERROR("second argument must be a point");
+    return result.release();
+  }
+
+  try {
+    result = SFCGAL::algorithm::visibility(poly->as<const SFCGAL::Polygon>(),
+                                           pt->as<const SFCGAL::Point>());
+  } catch (std::exception &e) {
+    SFCGAL_WARNING("During visibility(A, B) :");
+    SFCGAL_WARNING("  with A: %s", poly->asText().c_str());
+    SFCGAL_WARNING("  and B: %s", pt->asText().c_str());
+    SFCGAL_ERROR("%s", e.what());
+    return result.release();
+  }
+
+  return result.release();
+}
+
+extern "C" sfcgal_geometry_t *
+sfcgal_geometry_visibility_segment(const sfcgal_geometry_t *polygon,
+                                   const sfcgal_geometry_t *pointA,
+                                   const sfcgal_geometry_t *pointB)
+{
+  const auto *poly = reinterpret_cast<const SFCGAL::Geometry *>(polygon);
+  const auto *ptA  = reinterpret_cast<const SFCGAL::Geometry *>(pointA);
+  const auto *ptB  = reinterpret_cast<const SFCGAL::Geometry *>(pointB);
+  std::unique_ptr<SFCGAL::Geometry> result;
+
+  if (poly->geometryTypeId() != SFCGAL::TYPE_POLYGON) {
+    SFCGAL_ERROR("visibility() only applies to polygons");
+    return result.release();
+  }
+
+  if ((ptA->geometryTypeId() != SFCGAL::TYPE_POINT) ||
+      (ptB->geometryTypeId() != SFCGAL::TYPE_POINT)) {
+    SFCGAL_ERROR("second and third argument must be a point");
+    return result.release();
+  }
+
+  try {
+    result = SFCGAL::algorithm::visibility(poly->as<const SFCGAL::Polygon>(),
+                                           ptA->as<const SFCGAL::Point>(),
+                                           ptB->as<const SFCGAL::Point>());
+  } catch (std::exception &e) {
+    SFCGAL_WARNING("During visibility(A, B, C) :");
+    SFCGAL_WARNING("  with A: %s", poly->asText().c_str());
+    SFCGAL_WARNING("  and B: %s", ptA->asText().c_str());
+    SFCGAL_WARNING("  and C: %s", ptB->asText().c_str());
+    SFCGAL_ERROR("%s", e.what());
+    return result.release();
+  }
+
+  return result.release();
+}