File alphaShapes.cpp
File List > algorithm > alphaShapes.cpp
Go to the documentation of this file
// Copyright (c) 2022-2022, Oslandia.
// SPDX-License-Identifier: LGPL-2.0-or-later
#include "SFCGAL/algorithm/alphaShapes.h"
#include "SFCGAL/GeometryCollection.h"
#include "SFCGAL/MultiPolygon.h"
#include "SFCGAL/Polygon.h"
#include "SFCGAL/PolyhedralSurface.h"
#include "SFCGAL/Triangle.h"
#include "SFCGAL/TriangulatedSurface.h"
#include "SFCGAL/detail/GetPointsVisitor.h"
#include "SFCGAL/Exception.h"
#include "SFCGAL/Kernel.h"
#include <boost/format.hpp>
#include <CGAL/Alpha_shape_2.h>
#include <CGAL/Alpha_shape_face_base_2.h>
#include <CGAL/Alpha_shape_vertex_base_2.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/algorithm.h>
#include <CGAL/Arr_non_caching_segment_basic_traits_2.h>
#include <CGAL/Arrangement_2.h>
#include <vector>
namespace SFCGAL::algorithm {
using Vb = CGAL::Alpha_shape_vertex_base_2<Kernel>;
using Fb = CGAL::Alpha_shape_face_base_2<Kernel>;
using Tds = CGAL::Triangulation_data_structure_2<Vb, Fb>;
using Triangulation_2 = CGAL::Delaunay_triangulation_2<Kernel, Tds>;
using Point_2 = CGAL::Point_2<Kernel>;
using Segment_2 = CGAL::Segment_2<Kernel>;
using Alpha_shape_2 = CGAL::Alpha_shape_2<Triangulation_2>;
using Alpha_shape_edges_iterator = Alpha_shape_2::Alpha_shape_edges_iterator;
using Traits_2 = CGAL::Arr_non_caching_segment_basic_traits_2<Kernel>;
using Arrangement = CGAL::Arrangement_2<Traits_2>;
template <class OutputIterator>
void
alpha_edges(const Alpha_shape_2 &A, OutputIterator out)
{
auto it = A.alpha_shape_edges_begin();
auto end = A.alpha_shape_edges_end();
for (; it != end; ++it) {
if (A.classify(*it) == 2) {
*out++ = A.segment(*it);
}
}
}
static auto
computeAlpha(const Geometry &g, Alpha_shape_2 &alphaShape, double alpha = 0,
size_t nb_components = 1) -> double
{
using CGAL::object_cast;
double result = -1.0;
if (g.isEmpty()) {
return result;
}
SFCGAL::detail::GetPointsVisitor getPointVisitor;
const_cast<Geometry &>(g).accept(getPointVisitor);
// collect points
if (getPointVisitor.points.size() < 4) {
return result;
}
std::vector<Point_2> points;
points.reserve(getPointVisitor.points.size());
for (auto &point : getPointVisitor.points) {
points.push_back(point->toPoint_2());
}
std::vector<Segment_2> segments;
alphaShape.make_alpha_shape(points.begin(), points.end());
alphaShape.set_alpha(Kernel::FT(alpha));
alpha_edges(alphaShape, std::back_inserter(segments));
result = CGAL::to_double(*alphaShape.find_optimal_alpha(nb_components));
return result;
}
static auto
alpha_to_geometry(const Alpha_shape_2 &A, bool allow_holes)
-> std::unique_ptr<Geometry>
{
std::vector<Segment_2> segments;
alpha_edges(A, std::back_inserter(segments));
Arrangement arr;
CGAL::insert_non_intersecting_curves(arr, segments.begin(), segments.end());
auto poly{std::make_unique<Polygon>()};
for (auto f = arr.faces_begin(); f != arr.faces_end(); f++) {
auto ring{std::make_unique<LineString>()};
for (auto h = f->holes_begin(); h != f->holes_end(); h++) {
auto he = *h;
do {
ring->addPoint(he->source()->point());
} while (++he != *h);
}
if (ring->numPoints() > 3) {
ring->addPoint(ring->startPoint());
if (f->is_unbounded()) {
poly->setExteriorRing(ring.release());
} else if (allow_holes) {
poly->addInteriorRing(ring.release());
}
}
}
std::unique_ptr<Geometry> result = std::move(poly);
return result;
}
auto
optimal_alpha_shapes(const Geometry &g, bool allow_holes, size_t nb_components)
-> std::unique_ptr<Geometry>
{
Alpha_shape_2 A;
const double optimalAlpha{computeAlpha(g, A, 10000, nb_components)};
if (optimalAlpha < 0) {
return std::unique_ptr<Geometry>(new GeometryCollection());
}
A.set_alpha(optimalAlpha);
return alpha_to_geometry(A, allow_holes);
}
auto
alphaShapes(const Geometry &g, double alpha, bool allow_holes)
-> std::unique_ptr<Geometry>
{
using CGAL::object_cast;
Alpha_shape_2 A;
const double optimalAlpha{computeAlpha(g, A, alpha)};
if (optimalAlpha < 0) {
return std::unique_ptr<Geometry>(new GeometryCollection());
}
return alpha_to_geometry(A, allow_holes);
}
} // namespace SFCGAL::algorithm