From 21dabff9ee423fee14125b1a17903898408588a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Bartoletti?= Date: Tue, 26 Sep 2023 16:11:39 +0200 Subject: [PATCH] [feature] add Visibility methods based on CGAL Visibility_2 package --- src/algorithm/visibility.cpp | 201 ++++++++++++++++++++++ src/algorithm/visibility.h | 74 ++++++++ test/unit/SFCGAL/algorithm/Visibility.cpp | 156 +++++++++++++++++ 3 files changed, 431 insertions(+) create mode 100644 src/algorithm/visibility.cpp create mode 100644 src/algorithm/visibility.h create mode 100644 test/unit/SFCGAL/algorithm/Visibility.cpp diff --git a/src/algorithm/visibility.cpp b/src/algorithm/visibility.cpp new file mode 100644 index 00000000..64b100e5 --- /dev/null +++ b/src/algorithm/visibility.cpp @@ -0,0 +1,201 @@ +// Copyright (c) 2023-2023, Oslandia. +// SPDX-License-Identifier: LGPL-2.0-or-later +#include +#include + +#include + +#include + +#include +#include +#include + +#include +#include + +namespace SFCGAL { +namespace algorithm { + +using Point_2 = Kernel::Point_2; +using Polygon_2 = CGAL::Polygon_2; + +using Segment_2 = Kernel::Segment_2; +using Traits_2 = CGAL::Arr_segment_traits_2; +using Arrangement_2 = CGAL::Arrangement_2; +using Face_handle = Arrangement_2::Face_handle; +using Halfedge_const_handle = Arrangement_2::Halfedge_const_handle; +using TEV = + CGAL::Triangular_expansion_visibility_2; + +using PolygonWithHoles = CGAL::Polygon_with_holes_2; +/// +/// +/// +auto +visibility(const Geometry &polygon, const Geometry &point) + -> std::unique_ptr +{ + SFCGAL_ASSERT_GEOMETRY_VALIDITY_2D(polygon); + + std::unique_ptr result( + visibility(polygon, point, NoValidityCheck())); + propagateValidityFlag(*result, true); + return result; +} + +auto +visibility(const Geometry &polygon, const Geometry &point, NoValidityCheck) + -> std::unique_ptr +{ + + Point_2 queryPoint{point.as().toPoint_2()}; + std::unique_ptr extRing{new LineString()}; + + // insert geometry into the arrangement + CGAL::Polygon_with_holes_2 pwh{ + polygon.as().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()); + + // Find the face + Arrangement_2::Face_const_handle *face; + CGAL::Arr_naive_point_location pl(arr); + CGAL::Arr_point_location_result::Type obj = + pl.locate(queryPoint); + + // The query point locates in the interior of a face + face = boost::get(&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 visbility that way + if (face != NULL) { + fh = tev.compute_visibility(queryPoint, *face, output_arr); + } else { + // If the point 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()) { + throw(std::exception()); + } + + 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); + } + // 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())); + } + } + + std::unique_ptr result{new Polygon(extRing.release())}; + + return result; +} +/// +/// +/// +auto +visibility(const Geometry &polygon, const Geometry &pointA, + const Geometry &pointB) -> std::unique_ptr +{ + SFCGAL_ASSERT_GEOMETRY_VALIDITY_2D(polygon); + + std::unique_ptr 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 +{ + + Point_2 startPoint{pointA.as().toPoint_2()}; + Point_2 endPoint{pointB.as().toPoint_2()}; + + // insert geometry into the arrangement + CGAL::Polygon_with_holes_2 pwh{ + polygon.as().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 in a boundary segment, find the corresponding half edge + Halfedge_const_handle he = arr.halfedges_begin(); + bool cont = !Segment_2(he->source()->point(), he->target()->point()) + .has_on(startPoint) || + he->source()->point() == startPoint || 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()) { + throw(std::exception()); + } + + cont = !Segment_2(he->source()->point(), he->target()->point()) + .has_on(startPoint) || + he->source()->point() == startPoint || he->face()->is_unbounded(); + } + + // visibility query + Arrangement_2 output_arr; + TEV tev(arr); + Face_handle fh = tev.compute_visibility(startPoint, he, output_arr); + // print out the visibility region. + Arrangement_2::Ccb_halfedge_circulator curr = fh->outer_ccb(); + // print_arrangement(arr); + // print_arrangement(output_arr); + std::unique_ptr extRing{new LineString()}; + extRing->addPoint(Point(curr->target()->point())); + while (++curr != fh->outer_ccb()) { + extRing->addPoint(Point(curr->target()->point())); + } + std::unique_ptr result{new Polygon(extRing.release())}; + + return result; +} + +} // namespace algorithm +} // namespace SFCGAL diff --git a/src/algorithm/visibility.h b/src/algorithm/visibility.h new file mode 100644 index 00000000..c0629f83 --- /dev/null +++ b/src/algorithm/visibility.h @@ -0,0 +1,74 @@ +// Copyright (c) 2023-2023, Oslandia. +// SPDX-License-Identifier: LGPL-2.0-or-later + +#ifndef _SFCGAL_ALGORITHM_VISIBILITY2D_H_ +#define _SFCGAL_ALGORITHM_VISIBILITY2D_H_ + +#include + +#include + +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 + */ +SFCGAL_API auto +visibility(const Geometry &polygon, const Geometry &point) + -> std::unique_ptr; + +/** + * @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 + * @warning No actual validity check is done + */ +SFCGAL_API auto +visibility(const Geometry &polygon, const Geometry &point, NoValidityCheck) + -> std::unique_ptr; + +/** + * @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 and adjacents + */ +SFCGAL_API auto +visibility(const Geometry &polygon, const Geometry &pointA, + const Geometry &pointB) -> std::unique_ptr; + +/** + * @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 and adjacents + */ +SFCGAL_API auto +visibility(const Geometry &polygon, const Geometry &pointA, + const Geometry &pointB, NoValidityCheck) -> std::unique_ptr; + +} // namespace algorithm +} // namespace SFCGAL + +#endif diff --git a/test/unit/SFCGAL/algorithm/Visibility.cpp b/test/unit/SFCGAL/algorithm/Visibility.cpp new file mode 100644 index 00000000..3e616025 --- /dev/null +++ b/test/unit/SFCGAL/algorithm/Visibility.cpp @@ -0,0 +1,156 @@ +/** + * SFCGAL + * + * Copyright (C) 2012-2013 Oslandia + * Copyright (C) 2012-2013 IGN (http://www.ign.fr) + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Library General Public + * License as published by the Free Software Foundation; either + * version 2 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Library General Public License for more details. + + * You should have received a copy of the GNU Library General Public + * License along with this library; if not, see + . + */ +#include + +#include +#include +#include + +using namespace SFCGAL; + +#include "../../../test_config.h" +// always after CGAL +using namespace boost::unit_test; + +BOOST_AUTO_TEST_SUITE(SFCGAL_algorithm_VisibilityTest) + +// algorithm::alphaShapes + +BOOST_AUTO_TEST_CASE(testVisibility_PointInPolygon) +{ + std::vector points; + points.push_back(Point(0.0, 4.0)); + points.push_back(Point(0.0, 0.0)); + points.push_back(Point(3.0, 2.0)); + points.push_back(Point(4.0, 0.0)); + points.push_back(Point(4.0, 4.0)); + points.push_back(Point(1.0, 2.0)); + points.push_back(Point(0.0, 4.0)); + + LineString lineString(points); + Polygon poly(lineString); + + Point queryPoint(0.5, 2.0); + + std::unique_ptr result(algorithm::visibility(poly, queryPoint)); + std::string expectedWkt = "POLYGON((3.0 2.0,1.0 2.0,0.0 4.0,0.0 0.0))"; + BOOST_CHECK_EQUAL(result->asText(1), expectedWkt); +} + +BOOST_AUTO_TEST_CASE(testVisibility_PointInPolygonHole) +{ + std::vector points; + points.push_back(Point(0.0, 4.0)); + points.push_back(Point(0.0, 0.0)); + points.push_back(Point(3.0, 2.0)); + points.push_back(Point(4.0, 0.0)); + points.push_back(Point(4.0, 4.0)); + points.push_back(Point(1.0, 2.0)); + points.push_back(Point(0.0, 4.0)); + + std::vector points_hole; + points_hole.push_back(Point(0.2, 1.75)); + points_hole.push_back(Point(0.9, 1.8)); + points_hole.push_back(Point(0.7, 1.2)); + points_hole.push_back(Point(0.2, 1.75)); + + LineString lineString(points); + LineString hole(points_hole); + + Polygon poly(lineString); + poly.addInteriorRing(hole); + + Point queryPoint(0.5, 2.0); + + std::unique_ptr result(algorithm::visibility(poly, queryPoint)); + std::string expectedWkt = + "POLYGON((0.0 1.6,0.2 1.8,0.9 1.8,1.9 1.3,3.0 2.0,1.0 2.0,0.0 4.0))"; + BOOST_CHECK_EQUAL(result->asText(1), expectedWkt); +} + +BOOST_AUTO_TEST_CASE(testVisibility_SegmentInPolygon) +{ + std::vector points; + points.push_back(Point(0.0, 4.0)); + points.push_back(Point(0.0, 0.0)); + points.push_back(Point(3.0, 2.0)); + points.push_back(Point(4.0, 0.0)); + points.push_back(Point(4.0, 4.0)); + points.push_back(Point(1.0, 2.0)); + points.push_back(Point(0.0, 4.0)); + + LineString lineString(points); + Polygon poly(lineString); + + Point startPoint(1.0, 2.0); + Point endPoint(4.0, 4.0); + + std::unique_ptr result( + algorithm::visibility(poly, startPoint, endPoint)); + // std::string expectedWkt ="POLYGON((3.0 2.0,1.0 2.0,0.0 4.0,0.0 0.0))"; + // BOOST_CHECK_EQUAL(result->asText(1), expectedWkt); + std::cout << poly.asText(1) << "\n"; + std::cout << result->asText(1) << "\n"; +} + +BOOST_AUTO_TEST_CASE(testVisibility_SegmentInPolygonHole) +{ + std::vector points; + points.push_back(Point(1.0, 2.0)); + points.push_back(Point(12.0, 3.0)); + points.push_back(Point(19.0, -2.0)); + points.push_back(Point(12.0, 6.0)); + points.push_back(Point(14.0, 14.0)); + points.push_back(Point(9.0, 5.0)); + points.push_back(Point(1.0, 2.0)); + + std::vector points_hole1; + points_hole1.push_back(Point(8.0, 3.0)); + points_hole1.push_back(Point(8.0, 4.0)); + points_hole1.push_back(Point(10.0, 3.0)); + points_hole1.push_back(Point(8.0, 3.0)); + std::vector points_hole2; + points_hole2.push_back(Point(10.0, 6.0)); + points_hole2.push_back(Point(11.0, 7.0)); + points_hole2.push_back(Point(11.0, 6.0)); + points_hole2.push_back(Point(10.0, 6.0)); + + LineString lineString(points); + LineString hole1(points_hole1); + LineString hole2(points_hole2); + + Polygon poly(lineString); + poly.addInteriorRing(hole1); + poly.addInteriorRing(hole2); + + Point startPoint(19.0, -2.0); + Point endPoint(12.0, 6.0); + + std::cout << poly.asText(1) << "\n"; + + std::unique_ptr result( + algorithm::visibility(poly, startPoint, endPoint)); + // std::string expectedWkt = + // "POLYGON((0.0 1.6,0.2 1.8,0.9 1.8,1.9 1.3,3.0 2.0,1.0 2.0,0.0 4.0))"; + // BOOST_CHECK_EQUAL(result->asText(1), expectedWkt); + std::cout << result->asText(1) << "\n"; +} +BOOST_AUTO_TEST_SUITE_END()