From 1cdab4d669769b07f27cd99f10c4a690e6857da5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lo=C3=AFc=20Bartoletti?= Date: Fri, 23 Aug 2024 15:20:21 +0200 Subject: [PATCH] Add Buffer3D Algorithm --- src/PolyhedralSurface.cpp | 17 ++ src/PolyhedralSurface.h | 6 + src/algorithm/buffer3D.cpp | 293 ++++++++++++++++++++ src/algorithm/buffer3D.h | 85 ++++++ test/unit/SFCGAL/algorithm/Buffer3DTest.cpp | 115 ++++++++ 5 files changed, 516 insertions(+) create mode 100644 src/algorithm/buffer3D.cpp create mode 100644 src/algorithm/buffer3D.h create mode 100644 test/unit/SFCGAL/algorithm/Buffer3DTest.cpp diff --git a/src/PolyhedralSurface.cpp b/src/PolyhedralSurface.cpp index 815b8672..039a26ec 100644 --- a/src/PolyhedralSurface.cpp +++ b/src/PolyhedralSurface.cpp @@ -25,6 +25,23 @@ PolyhedralSurface::PolyhedralSurface(const std::vector &polygons) } } +PolyhedralSurface::PolyhedralSurface(const std::unique_ptr &geometry) +{ + if (geometry->is()) { + *this = static_cast(*geometry); + } else if (geometry->is()) { + const TriangulatedSurface &triangulatedSurface = + geometry->as(); + for (size_t i = 0; i < triangulatedSurface.numTriangles(); ++i) { + this->addPolygon(triangulatedSurface.triangleN(i)); + } + } else if (geometry->is()) { + this->addPolygon(geometry->as()); + } else { + throw std::invalid_argument("Cannot convert geometry to PolyhedralSurface"); + } +} + /// /// /// diff --git a/src/PolyhedralSurface.h b/src/PolyhedralSurface.h index 8649b926..8698fa4f 100644 --- a/src/PolyhedralSurface.h +++ b/src/PolyhedralSurface.h @@ -42,6 +42,12 @@ class SFCGAL_API PolyhedralSurface : public Surface { * Constructor with a vector of polygons */ PolyhedralSurface(const std::vector &polygons); + + /** + * Constructor with a Geometry + */ + PolyhedralSurface(const std::unique_ptr &geometry); + /** * Constructor from a Polyhedron (detail::MarkedPolyhedron or * CGAL::Polyhedron_3) diff --git a/src/algorithm/buffer3D.cpp b/src/algorithm/buffer3D.cpp new file mode 100644 index 00000000..795fedbf --- /dev/null +++ b/src/algorithm/buffer3D.cpp @@ -0,0 +1,293 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace SFCGAL { +namespace algorithm { + +Buffer3D::Buffer3D(const Geometry &inputGeometry, double radius, int segments) + : _radius(radius), _segments(segments) +{ + if (inputGeometry.is()) { + _inputPoints.push_back(inputGeometry.as().toPoint_3()); + } else if (inputGeometry.is()) { + const LineString &ls = inputGeometry.as(); + for (size_t i = 0; i < ls.numPoints(); ++i) { + _inputPoints.push_back(ls.pointN(i).toPoint_3()); + } + } else { + throw std::invalid_argument("Input geometry must be a Point or LineString"); + } +} + +PolyhedralSurface +Buffer3D::compute(BufferType type) const +{ + if (_inputPoints.size() == 1) { + return computePointBuffer(); + } + + switch (type) { + case ROUND: + return computeRoundBuffer(); + case CYLSPHERE: + return computeCylSphereBuffer(); + case FLAT: + return computeFlatBuffer(); + default: + throw std::invalid_argument("Invalid buffer type"); + } +} + +PolyhedralSurface +Buffer3D::computePointBuffer() const +{ + Kernel::Point_3 center(_inputPoints[0].x(), _inputPoints[0].y(), + _inputPoints[0].z()); + Sphere sphere(_radius, center, _segments / 2, _segments); + return PolyhedralSurface(sphere.generatePolyhedron()); +} + +PolyhedralSurface +Buffer3D::computeRoundBuffer() const +{ + typedef Kernel::Point_3 Point_3; + typedef Point_3 *point_iterator; + typedef std::pair point_range; + typedef std::list polyline; + typedef CGAL::Nef_polyhedron_3 Nef_polyhedron; + + // Create sphere + Point_3 center(0, 0, 0); + SFCGAL::Sphere sphere(_radius, center, _segments / 2, _segments); + + // Generate polyhedron from sphere + CGAL::Polyhedron_3 spherePolyhedron = sphere.generatePolyhedron(); + + // Convert Polyhedron to a Nef_polyhedron + Nef_polyhedron N0(spherePolyhedron); + + // Create a polyline from _inputPoints + polyline poly; + if (!_inputPoints.empty()) { + // Create a non-const copy of the points + std::vector points_copy(_inputPoints.begin(), _inputPoints.end()); + poly.push_back(point_range(&points_copy.front(), &points_copy.back() + 1)); + + // Create a Nef_polyhedron from the polyline + Nef_polyhedron N1(poly.begin(), poly.end(), + Nef_polyhedron::Polylines_tag()); + + // Perform Minkowski sum + Nef_polyhedron result = CGAL::minkowski_sum_3(N0, N1); + + // Convert result to SFCGAL::PolyhedralSurface + SFCGAL::detail::MarkedPolyhedron out; + result.convert_to_polyhedron(out); + return PolyhedralSurface(out); + } else { + // If _inputPoints is empty, return an empty PolyhedralSurface + return PolyhedralSurface(); + } +} + +PolyhedralSurface +Buffer3D::computeCylSphereBuffer() const +{ + std::vector> meshes; + + for (size_t i = 0; i < _inputPoints.size() - 1; ++i) { + // Create and add cylinder + Kernel::Vector_3 axis(_inputPoints[i + 1].x() - _inputPoints[i].x(), + _inputPoints[i + 1].y() - _inputPoints[i].y(), + _inputPoints[i + 1].z() - _inputPoints[i].z()); + Kernel::FT height = CGAL::sqrt(CGAL::to_double(axis.squared_length())); + Kernel::Point_3 base(_inputPoints[i].x(), _inputPoints[i].y(), + _inputPoints[i].z()); + Cylinder cyl(base, axis, _radius, height, _segments); + meshes.push_back(cyl.generatePolyhedron()); + + // Add sphere at joints + if (i > 0) { + Kernel::Point_3 sphereCenter(_inputPoints[i].x(), _inputPoints[i].y(), + _inputPoints[i].z()); + Sphere sphere(_radius, sphereCenter, _segments / 2, _segments); + meshes.push_back(sphere.generatePolyhedron()); + } + } + + // Add sphere at the end point + Kernel::Point_3 endSphereCenter(_inputPoints.back().x(), + _inputPoints.back().y(), + _inputPoints.back().z()); + Sphere endSphere(_radius, endSphereCenter, _segments / 2, _segments); + meshes.push_back(endSphere.generatePolyhedron()); + + // Combine all meshes into a single PolyhedralSurface + // TODO: Use Union3D to merge interior faces? + PolyhedralSurface poly_surface; + for (const auto &mesh : meshes) { + poly_surface.addPolygons(mesh); + } + + return poly_surface; +} + +PolyhedralSurface +Buffer3D::computeFlatBuffer() const +{ + std::vector line_points; + for (const auto &p : _inputPoints) { + line_points.emplace_back(p.x(), p.y(), p.z()); + } + + Cylinder::Surface_mesh buffer; + std::vector> rings; + + std::vector bisector_planes; + for (size_t i = 1; i < line_points.size() - 1; ++i) { + bisector_planes.push_back(compute_bisector_plane( + line_points[i - 1], line_points[i], line_points[i + 1])); + } + + for (size_t i = 0; i < line_points.size() - 1; ++i) { + Kernel::Vector_3 axis = normalize(line_points[i + 1] - line_points[i]); + Kernel::Point_3 extended_start = + extend_point(line_points[i], -axis, _radius); + Kernel::Point_3 extended_end = + extend_point(line_points[i + 1], axis, _radius); + + std::vector start_circle = + create_circle_points(extended_start, axis, _radius, _segments); + std::vector end_circle = + create_circle_points(extended_end, axis, _radius, _segments); + + std::vector start_ring, end_ring; + + for (size_t j = 0; j < _segments; ++j) { + Kernel::Point_3 start_point = start_circle[j]; + Kernel::Point_3 end_point = end_circle[j]; + + if (i > 0) { + start_point = intersect_segment_plane(start_point, end_point, + bisector_planes[i - 1]); + } + if (i < line_points.size() - 2) { + end_point = + intersect_segment_plane(start_point, end_point, bisector_planes[i]); + } + + Cylinder::Surface_mesh::Vertex_index v1 = buffer.add_vertex(start_point); + Cylinder::Surface_mesh::Vertex_index v2 = buffer.add_vertex(end_point); + start_ring.push_back(v1); + end_ring.push_back(v2); + + if (j > 0) { + buffer.add_face(start_ring[j - 1], start_ring[j], end_ring[j], + end_ring[j - 1]); + } + } + buffer.add_face(start_ring.back(), start_ring.front(), end_ring.front(), + end_ring.back()); + + rings.push_back(start_ring); + if (i == line_points.size() - 2) + rings.push_back(end_ring); + } + + // Add caps + Kernel::Point_3 start_center = + extend_point(line_points.front(), + normalize(line_points[1] - line_points.front()), -_radius); + Kernel::Point_3 end_center = extend_point( + line_points.back(), + normalize(line_points.back() - line_points[line_points.size() - 2]), + _radius); + Cylinder::Surface_mesh::Vertex_index start_center_index = + buffer.add_vertex(start_center); + Cylinder::Surface_mesh::Vertex_index end_center_index = + buffer.add_vertex(end_center); + + for (size_t i = 0; i < _segments; ++i) { + buffer.add_face(start_center_index, rings.front()[(i + 1) % _segments], + rings.front()[i]); + buffer.add_face(end_center_index, rings.back()[i], + rings.back()[(i + 1) % _segments]); + } + + return PolyhedralSurface(buffer); +} + +Kernel::Vector_3 +Buffer3D::normalize(const Kernel::Vector_3 &v) const +{ + double length = std::sqrt(CGAL::to_double(v.squared_length())); + return v / length; +} + +Kernel::Point_3 +Buffer3D::extend_point(const Kernel::Point_3 &point, + const Kernel::Vector_3 &direction, double distance) const +{ + return point + direction * distance; +} + +std::vector +Buffer3D::create_circle_points(const Kernel::Point_3 ¢er, + const Kernel::Vector_3 &axis, double radius, + int segments) const +{ + std::vector points; + Kernel::Vector_3 perpendicular = + CGAL::cross_product(axis, Kernel::Vector_3(0, 0, 1)); + if (perpendicular == CGAL::NULL_VECTOR) { + perpendicular = CGAL::cross_product(axis, Kernel::Vector_3(0, 1, 0)); + } + perpendicular = normalize(perpendicular); + Kernel::Vector_3 perpendicular2 = + normalize(CGAL::cross_product(axis, perpendicular)); + + for (int i = 0; i < segments; ++i) { + double angle = 2.0 * M_PI * i / segments; + Kernel::Vector_3 offset = radius * (std::cos(angle) * perpendicular + + std::sin(angle) * perpendicular2); + points.push_back(center + offset); + } + return points; +} + +Kernel::Plane_3 +Buffer3D::compute_bisector_plane(const Kernel::Point_3 &p1, + const Kernel::Point_3 &p2, + const Kernel::Point_3 &p3) const +{ + Kernel::Vector_3 v1 = normalize(p2 - p1); + Kernel::Vector_3 v2 = normalize(p3 - p2); + Kernel::Vector_3 bisector = v1 + v2; + return Kernel::Plane_3(p2, bisector); +} + +Kernel::Point_3 +Buffer3D::intersect_segment_plane(const Kernel::Point_3 &p1, + const Kernel::Point_3 &p2, + const Kernel::Plane_3 &plane) const +{ + Kernel::Vector_3 v = p2 - p1; + Kernel::FT t = + -plane.a() * p1.x() - plane.b() * p1.y() - plane.c() * p1.z() - plane.d(); + t /= plane.a() * v.x() + plane.b() * v.y() + plane.c() * v.z(); + return p1 + t * v; +} + +} // namespace algorithm +} // namespace SFCGAL diff --git a/src/algorithm/buffer3D.h b/src/algorithm/buffer3D.h new file mode 100644 index 00000000..bd10d9c4 --- /dev/null +++ b/src/algorithm/buffer3D.h @@ -0,0 +1,85 @@ +#include "SFCGAL//io/OBJ.h" +#include "SFCGAL/Cylinder.h" +#include "SFCGAL/Geometry.h" +#include "SFCGAL/LineString.h" +#include "SFCGAL/Point.h" +#include "SFCGAL/PolyhedralSurface.h" +#include "SFCGAL/Sphere.h" +#include "SFCGAL/algorithm/minkowskiSum3D.h" +#include +#include +#include + +namespace SFCGAL { +namespace algorithm { + +/** + * @brief Computes a 3D buffer around a Point or LineString. + * @ingroup public_api + */ +class Buffer3D { +public: + /** + * @brief Buffer type enumeration + */ + enum BufferType { + ROUND, ///< Minkowski sum with a sphere + CYLSPHERE, ///< Union of cylinders and spheres + FLAT ///< Construction of a disk on the bisector plane + }; + + /** + * @brief Constructs a Buffer3D object + * @param inputGeometry The input geometry (must be a Point or LineString) + * @param radius The buffer radius + * @param segments The number of segments used to approximate curved surfaces + * @throws std::invalid_argument if the input geometry is not a Point or + * LineString + */ + Buffer3D(const Geometry &inputGeometry, double radius, int segments); + + /** + * @brief Computes the 3D buffer + * @param type The type of buffer to compute + * @return A PolyhedralSurface representing the 3D buffer + * @throws std::invalid_argument if an invalid buffer type is provided + */ + PolyhedralSurface + compute(BufferType type) const; + +private: + std::vector _inputPoints; + double _radius; + int _segments; + + PolyhedralSurface + computePointBuffer() const; + PolyhedralSurface + computeRoundBuffer() const; + PolyhedralSurface + computeCylSphereBuffer() const; + PolyhedralSurface + computeFlatBuffer() const; + + // Helper functions for FLAT buffer + CGAL::Vector_3 + normalize(const CGAL::Vector_3 &v) const; + CGAL::Point_3 + extend_point(const CGAL::Point_3 &point, + const CGAL::Vector_3 &direction, double distance) const; + std::vector> + create_circle_points(const CGAL::Point_3 ¢er, + const CGAL::Vector_3 &axis, double radius, + int segments) const; + CGAL::Plane_3 + compute_bisector_plane(const CGAL::Point_3 &p1, + const CGAL::Point_3 &p2, + const CGAL::Point_3 &p3) const; + CGAL::Point_3 + intersect_segment_plane(const CGAL::Point_3 &p1, + const CGAL::Point_3 &p2, + const CGAL::Plane_3 &plane) const; +}; + +} // namespace algorithm +} // namespace SFCGAL diff --git a/test/unit/SFCGAL/algorithm/Buffer3DTest.cpp b/test/unit/SFCGAL/algorithm/Buffer3DTest.cpp new file mode 100644 index 00000000..998c363a --- /dev/null +++ b/test/unit/SFCGAL/algorithm/Buffer3DTest.cpp @@ -0,0 +1,115 @@ +#include + +#include "SFCGAL/GeometryCollection.h" +#include "SFCGAL/LineString.h" +#include "SFCGAL/MultiLineString.h" +#include "SFCGAL/MultiPoint.h" +#include "SFCGAL/MultiPolygon.h" +#include "SFCGAL/MultiSolid.h" +#include "SFCGAL/Point.h" +#include "SFCGAL/Polygon.h" +#include "SFCGAL/PolyhedralSurface.h" +#include "SFCGAL/Solid.h" +#include "SFCGAL/Triangle.h" +#include "SFCGAL/TriangulatedSurface.h" +#include "SFCGAL/algorithm/buffer3D.h" +#include "SFCGAL/io/wkt.h" + +using namespace SFCGAL; + +#include "../../../test_config.h" +// always after CGAL +using namespace boost::unit_test; + +BOOST_AUTO_TEST_SUITE(SFCGAL_algorithm_Buffer3DTest) + +// algorithm::buffer3D + +BOOST_AUTO_TEST_CASE(testBuffer3D_Point) +{ + double radius = 10.0; + int segments = 16; + + Point point(0, 0, 0); + algorithm::Buffer3D buffer3d(point, radius, segments); + + std::vector bufferTypes = { + algorithm::Buffer3D::ROUND, algorithm::Buffer3D::CYLSPHERE, + algorithm::Buffer3D::FLAT}; + + for (auto bufferType : bufferTypes) { + PolyhedralSurface buffer = buffer3d.compute(bufferType); + BOOST_CHECK(buffer.is3D()); + BOOST_CHECK(buffer.numGeometries() > 0); + + std::string typeName; + switch (bufferType) { + case algorithm::Buffer3D::ROUND: + typeName = "ROUND"; + break; + case algorithm::Buffer3D::CYLSPHERE: + typeName = "CYLSPHERE"; + break; + case algorithm::Buffer3D::FLAT: + typeName = "FLAT"; + break; + } + + std::cout << "Point " << typeName << " Buffer 3D WKT: " << buffer.asText(1) + << std::endl; + SFCGAL::io::OBJ::save(buffer, "/tmp/point_" + typeName + "_buffer_3d.obj"); + } +} + +BOOST_AUTO_TEST_CASE(testBuffer3D_LineString) +{ + double radius = 10.0; + int segments = 16; + + std::vector points = {Point(-100, 0, 0), Point(40, -70, 0), + Point(40, 50, 40), Point(-90, -60, 60), + Point(0, 0, -100), Point(30, 0, 150)}; + LineString lineString(points); + + algorithm::Buffer3D buffer3d(lineString, radius, segments); + + std::vector bufferTypes = { + algorithm::Buffer3D::ROUND, algorithm::Buffer3D::CYLSPHERE, + algorithm::Buffer3D::FLAT}; + + for (auto bufferType : bufferTypes) { + PolyhedralSurface buffer = buffer3d.compute(bufferType); + BOOST_CHECK(buffer.is3D()); + BOOST_CHECK(buffer.numGeometries() > 0); + + std::string typeName; + switch (bufferType) { + case algorithm::Buffer3D::ROUND: + typeName = "ROUND"; + break; + case algorithm::Buffer3D::CYLSPHERE: + typeName = "CYLSPHERE"; + break; + case algorithm::Buffer3D::FLAT: + typeName = "FLAT"; + break; + } + + std::cout << "LineString " << typeName + << " Buffer 3D WKT: " << buffer.asText(1) << std::endl; + SFCGAL::io::OBJ::save(buffer, + "/tmp/linestring_" + typeName + "_buffer_3d.obj"); + } +} + +BOOST_AUTO_TEST_CASE(testBuffer3D_InvalidGeometry) +{ + double radius = 10.0; + int segments = 16; + + Triangle triangle(Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0)); + BOOST_CHECK_THROW(algorithm::Buffer3D(triangle, radius, segments), + std::invalid_argument); +} + +BOOST_AUTO_TEST_SUITE_END()