EIC code displayed by LXR master/​include/​boost/​geometry/​algorithms/​detail/​is_valid/​polyhedral_surface.hpp	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-08-20 08:44:11

 // Boost.Geometry
 
 // Copyright (c) 2025 Oracle and/or its affiliates.
 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
 
 // Use, modification and distribution is subject to the Boost Software License,
 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
 // http://www.boost.org/LICENSE_1_0.txt)
 
 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_IS_VALID_POLYHEDRAL_SURFACE_HPP
 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_IS_VALID_POLYHEDRAL_SURFACE_HPP
 
 #include <boost/geometry/algorithms/buffer.hpp>
 #include <boost/geometry/algorithms/correct.hpp>
 #include <boost/geometry/algorithms/detail/is_valid/interface.hpp>
 #include <boost/geometry/algorithms/detail/is_valid/polygon.hpp>
 #include <boost/geometry/algorithms/detail/calculate_point_order.hpp>
 #include <boost/geometry/algorithms/distance.hpp> //needed for rtree
 #include <boost/geometry/index/predicates.hpp>
 #include <boost/geometry/index/rtree.hpp>
 #include <boost/geometry/io/wkt/wkt.hpp>
 #include <boost/geometry/iterators/point_reverse_iterator.hpp>
 
 #include <vector>
 #include <stack>
 
 namespace boost { namespace geometry
 {
 
 #ifndef DOXYGEN_NO_DETAIL
 namespace detail { namespace is_valid
 {
 
 struct polygon_frontier_type
 {
     polygon_frontier_type(std::size_t size)
         : m_stack(), m_visited(size, false)
     {}
 
     void push(std::size_t index)
     {
         if (!m_visited[index])
         {
             m_stack.push(index);
             m_visited[index] = true;
         }
     }
     bool empty()
     {
         return m_stack.empty();
     }
     std::size_t top()
     {
         return m_stack.top();
     }
     void pop()
     {
         m_stack.pop();
     }
 private :
     std::stack<std::size_t> m_stack;
     std::vector<bool> m_visited;
 };
 
 // Implementation of the polygon segment intersection algorithm
 // This function checks if a segment (s1, s2) intersects with a 2D polygon in 3D space.
 // proj1 and proj2 are the indices of the coordinate axes used to determine the sign of three points
 // of the polygon after projection.
 // It returns:
 // - 0 if the segment does not intersect the polygon
 // - 1 if the segment intersects the polygon in the interior
 // - 2 if the segment intersects the polygon at a vertex
 // - 3 if the segment intersects the polygon at an edge
 // Reference: R.J.Segura, F.R.Feito - An algorithm for determining intersection segment-polygon in
 // 3D, Computers & Graphics 1998
 template <typename Point, typename Polygon, typename Strategy>
 int polygon_segment_intersection(Point const& s1,
                                  Point const& s2,
                                  Polygon const& polygon,
                                  int proj1,
                                  int proj2,
                                  Strategy const& strategy)
 {
     BOOST_GEOMETRY_ASSERT(geometry::dimension<Point>::value == 3);
     BOOST_GEOMETRY_ASSERT(geometry::dimension<Polygon>::value == 3);
 
     auto const pit_a = points_begin(polygon);
     auto const& p_a = *pit_a;
 
     int cut = 0;
 
     for (auto pit = boost::next(pit_a); pit != boost::prior(points_end(polygon)); ++pit)
     {
         auto const& p_b = *pit;
         auto const& p_c = *boost::next(pit);
 
         int const sign = strategy.side().apply(p_a, p_b, p_c, proj1, proj2);
 
         int const sign1 = strategy.side3().apply(s1, p_a, s2, p_b);
         int const sign2 = strategy.side3().apply(s1, p_c, p_b, s2);
         int const sign3 = strategy.side3().apply(s1, p_a, p_c, s2);
 
         // Insects at p_b or p_a or p_c -> vertex
         if ((sign1 == 0 && sign2 == 0) || (sign1 == 0 && sign3 == 0) || (sign2 == 0 && sign3 == 0))
         {
             return 2;
         }
 
         // Intersects at edge p_a, p_c -> internal edge
         if (sign1 == 0 && sign2 == sign3) cut += sign;
 
         // Intersects at edge p_b, p_c -> external edge
         if (sign2 == 0 && sign1 == sign3) return 3;
 
         // Intersects at edge p_a, p_b -> internal edge
         if (sign3 == 0 && sign1 == sign2) cut += sign;
 
         // Intersects inside p_a, p_b, p_c -> inside
         if (sign1 == sign2 && sign2 == sign3) cut += 2 * sign;
     }
 
     return (cut == 2) ? 1 : 0;
 }
 
 template <typename PolyhedralSurface>
 struct valid_intersection
 {
     valid_intersection(std::size_t num_of_polygons)
         : m_polygon_projections(num_of_polygons)
     {}
 
     using coordinate_type = typename coordinate_type<PolyhedralSurface>::type;
     using polygon_3d = typename boost::range_value<PolyhedralSurface>::type;
     using ring_3d = typename boost::geometry::ring_type<polygon_3d>::type;
     using point_3d = typename boost::geometry::point_type<ring_3d>::type;
     using segment_3d = typename boost::geometry::model::segment<point_3d>;
     using box_3d = typename boost::geometry::model::box<point_3d>;
 
     using point_2d = boost::geometry::model::point
         <coordinate_type, 2, boost::geometry::cs::cartesian>;
     using linestring_2d = typename boost::geometry::model::linestring<point_2d>;
     using polygon_2d = typename std::conditional
         <
             boost::geometry::point_order<polygon_3d>::value == boost::geometry::counterclockwise,
             typename boost::geometry::model::polygon<point_2d, true>,
             typename boost::geometry::model::polygon<point_2d, false>
         >::type;
     using ring_2d = typename boost::geometry::ring_type<polygon_2d>::type;
 
 private:
 
     enum class intersection_type
     {
         none,
         valid,
         invalid
     };
 
     struct polygon_projection
     {
         polygon_2d projected_polygon;
         int p1_idx;
         int p2_idx;
         int p3_idx;
         int proj1;
         int proj2;
         int dir;
     };
 
     template <typename Strategy>
     static std::tuple<int, int, int> affinely_independent(point_3d const& p1,
                                                           point_3d const& p2,
                                                           point_3d const& p3,
                                                           Strategy const& strategy)
     {
         int side = strategy.side().apply(p1, p2, p3, 0, 1);
         if (side != 0) return {0, 1, side};
 
         side = strategy.side().apply(p1, p2, p3, 0, 2);
         if (side != 0) return {0, 2, side};
 
         side = strategy.side().apply(p1, p2, p3, 1, 2);
         if (side != 0) return {1, 2, side};
 
         return {-1, -1, 0}; // All points are collinear
     }
 
     template <typename Strategy>
     static std::tuple<int, int, int> find_non_collinear_points(polygon_3d const& polygon,
                                                                int& p1_index,
                                                                int& p2_index,
                                                                int& p3_index,
                                                                Strategy const& strategy)
     {
         for (auto pit1 = points_begin(polygon); pit1 != points_end(polygon); ++pit1)
         {
             for (auto pit2 = boost::next(pit1); pit2 != points_end(polygon); ++pit2)
             {
                 for (auto pit3 = boost::next(pit2); pit3 != points_end(polygon); ++pit3)
                 {
                     // A set of 3 independent points is sufficient to construct a plane.
                     auto res = affinely_independent(*pit1, *pit2, *pit3, strategy);
                     if (std::get<0>(res) != -1)
                     {
                         p1_index = std::distance(points_begin(polygon), pit1);
                         p2_index = std::distance(points_begin(polygon), pit2);
                         p3_index = std::distance(points_begin(polygon), pit3);
                         return res;
                     }
                 }
             }
         }
         return {-1, -1, 0}; // All points are collinear
     }
 
     template <typename VisitPolicy>
     static bool has_correct_intersection_orientation(polygon_2d const& polygon,
                                                      point_2d const& segment_point1,
                                                      point_2d const& segment_point2,
                                                      VisitPolicy& visitor)
     {
         if (boost::geometry::equals(segment_point1, segment_point2))
         {
             visitor.template apply<failure_duplicate_points>();
             return false;
         }
 
         boost::geometry::closing_iterator<ring_2d> it(boost::geometry::exterior_ring(polygon));
         auto begin = it;
         boost::geometry::closing_iterator<ring_2d> end(boost::geometry::exterior_ring(polygon), true);
 
         for (; it != end - 1; ++it)
         {
             if (boost::geometry::equals(*it, segment_point1))
             {
                 auto prev = (it == begin) ? end - 1 : boost::prior(it);
 
                 if (boost::geometry::equals(segment_point2, *prev))
                 {
                     // valid intersection and consistent orientation
                     return true;
                 }
 
                 auto next = (it == end) ? begin : boost::next(it);
 
                 if (boost::geometry::equals(segment_point2, *next))
                 {
                     // valid intersection but inconsistent orientation
                     visitor.template apply<failure_inconsistent_orientation>();
                     return false;
                 }
             }
         }
 
         visitor.template apply<failure_invalid_intersection>();
 
         return false;
     }
 
     point_2d project_point_to_coordinate_plane(point_3d const& point, int proj1, int proj2)
     {
         if (proj1 == 1 && proj2 == 2) return point_2d{geometry::get<1>(point), geometry::get<2>(point)};
         if (proj1 == 0 && proj2 == 2) return point_2d{geometry::get<0>(point), geometry::get<2>(point)};
         if (proj1 == 0 && proj2 == 1) return point_2d{geometry::get<0>(point), geometry::get<1>(point)};
         return point_2d{};
     }
 
     // Check if a segment intersects a polygon in 3D space
     template <typename VisitPolicy, typename Strategy>
     intersection_type check_polygon_segment_intersection(point_3d const& segment_point1,
                                                          point_3d const& segment_point2,
                                                          polygon_3d const& polygon,
                                                          std::size_t const polygon_index,
                                                          VisitPolicy& visitor,
                                                          Strategy const& strategy)
     {
         auto const& projection = m_polygon_projections[polygon_index];
         auto const& projected_polygon = projection.projected_polygon;
 
         auto const& p1 = polygon.outer()[projection.p1_idx];
         auto const& p2 = polygon.outer()[projection.p2_idx];
         auto const& p3 = polygon.outer()[projection.p3_idx];
 
         int proj1 = projection.proj1;
         int proj2 = projection.proj2;
         int dir = projection.dir;
 
         // Check if the segment intersects the plane of the polygon
         auto sign1 = strategy.side3().apply(p1, p2, p3, segment_point1);
         auto sign2 = strategy.side3().apply(p1, p2, p3, segment_point2);
 
         if ((sign1 > 0 && sign2 > 0) || (sign1 < 0 && sign2 < 0))
         {
             return intersection_type::none; // No intersection, this is valid
         }
 
         // Check if the segment is on the same plane as the polygon
         if (sign1 == 0 && sign2 == 0)
         {
             auto const& projected_p1 = project_point_to_coordinate_plane(segment_point1, proj1, proj2);
             auto const& projected_p2 = project_point_to_coordinate_plane(segment_point2, proj1, proj2);
 
             linestring_2d projected_segment = {projected_p1, projected_p2};
 
             // Check if the projected segment intersects the polygon
             std::tuple
             <
                 boost::geometry::model::multi_point<point_2d>,
                 boost::geometry::model::multi_linestring<linestring_2d>> intersections;
             boost::geometry::intersection(projected_segment, projected_polygon, intersections);
 
             auto const& points = std::get<0>(intersections);
             auto const& linestrings = std::get<1>(intersections);
 
             if (boost::geometry::is_empty(points) && boost::geometry::is_empty(linestrings))
             {
                 return intersection_type::none; // No intersection, this is valid
             }
 
             if (!points.empty())
             {
                 if (points.size() == 1 &&
                     (boost::geometry::equals(points.front(), projected_p1) ||
                     boost::geometry::equals(points.front(), projected_p2)))
                 {
                     return intersection_type::none; // Valid empty intersection
                 }
                 visitor.template apply<failure_invalid_intersection>();
                 return intersection_type::invalid; // Invalid intersection
             }
 
             if (linestrings.size() == 1)
             {
                 if (boost::geometry::equals(linestrings.front(), projected_segment))
                 {
                     auto const& pp1 = dir < 0 ? projected_p2 : projected_p1;
                     auto const& pp2 = dir < 0 ? projected_p1 : projected_p2;
                     if (!has_correct_intersection_orientation(projected_polygon, pp1, pp2, visitor))
                     {
                         return intersection_type::invalid;
                     }
                     return intersection_type::valid; // The only case of valid non empty intersection
                 }
             }
             visitor.template apply<failure_invalid_intersection>();
             return intersection_type::invalid;
         }
 
         int inter = polygon_segment_intersection(segment_point1, segment_point2, polygon, proj1,
             proj2, strategy);
 
         if ((inter == 1 || inter == 3) || (inter == 2 && sign1 !=0 && sign2 != 0))
         {
             visitor.template apply<failure_invalid_intersection>();
             return intersection_type::invalid;
         }
 
         return intersection_type::none;
     }
 
     // Helper function to check if any segment of a ring intersects a polygon
     template <typename VisitPolicy, typename Strategy>
     intersection_type check_polygon_polygon_intersection(polygon_3d const& polygon1,
                                                          polygon_3d const& polygon2,
                                                          std::size_t const& index2,
                                                          VisitPolicy& visitor,
                                                          Strategy const& strategy)
     {
         bool valid_intersection = false;
         auto const& ring1 = boost::geometry::exterior_ring(polygon1);
 
         for (auto it = boost::begin(ring1); it + 1 != boost::end(ring1); ++it)
         {
             auto result = check_polygon_segment_intersection(*it, *(it + 1), polygon2, index2,
                 visitor, strategy);
 
             if (result == intersection_type::invalid)
             {
                 return intersection_type::invalid; // Found invalid intersection
             }
 
             if (result == intersection_type::valid)
             {
                 valid_intersection = true; // Found a valid intersection
             }
         }
         // Either no intersection or valid intersection
         return valid_intersection ? intersection_type::valid : intersection_type::none;
     }
 
  public:
 
     template <typename VisitPolicy, typename Strategy>
     intersection_type apply(polygon_3d const& polygon1,
                             polygon_3d const& polygon2,
                             std::size_t const& index1,
                             std::size_t const& index2,
                             VisitPolicy& visitor,
                             Strategy const& strategy)
     {
         auto intersection1 = check_polygon_polygon_intersection(polygon1, polygon2, index2,
             visitor, strategy);
 
         if (intersection1 == intersection_type::invalid)
         {
             return intersection_type::invalid; // Found invalid intersection
         }
 
         auto intersection2 = check_polygon_polygon_intersection(polygon2, polygon1, index1,
             visitor, strategy);
 
         if (intersection2 == intersection_type::invalid)
         {
             return intersection_type::invalid; // Found invalid intersection
         }
 
         if (intersection1 == intersection_type::none && intersection2 == intersection_type::none)
         {
             return intersection_type::none; // No intersection found
         }
 
         return intersection_type::valid;
     }
 
     template <typename VisitPolicy, typename Strategy>
     intersection_type initialize(polygon_3d const& polygon,
                                  std::size_t const polygon_index,
                                  VisitPolicy& visitor,
                                  Strategy const& strategy)
     {
         // TODO: Support interior rings
         ring_3d const& ring = boost::geometry::exterior_ring(polygon);
 
         if (boost::size(ring) < detail::minimum_ring_size<ring_3d>::value)
         {
             visitor.template apply<failure_few_points_on_face>();
             return intersection_type::invalid;
         }
 
         // Get three non-collinear points from the polygon to define the plane
         int p1_idx = 0;
         int p2_idx = 0;
         int p3_idx = 0;
 
         auto res = find_non_collinear_points(polygon, p1_idx, p2_idx, p3_idx, strategy);
 
         if (std::get<0>(res) == -1)
         {
             visitor.template apply<failure_collinear_points_on_face>();
             return intersection_type::invalid;
         }
 
         int proj1 = std::get<0>(res);
         int proj2 = std::get<1>(res);
         int dir = std::get<2>(res);
 
         auto const& p1 = polygon.outer()[p1_idx];
         auto const& p2 = polygon.outer()[p2_idx];
         auto const& p3 = polygon.outer()[p3_idx];
 
         // Check for non-coplanar points
         for (auto pit = points_begin(polygon); pit != points_end(polygon); ++pit)
         {
             if (strategy.side3().apply(p1, p2, p3, *pit) != 0)
             {
                 visitor.template apply<failure_non_coplanar_points_on_face>();
                 return intersection_type::invalid;
             }
         }
 
         // Project the polygon
         polygon_2d projected_polygon;
 
         auto add_projected_points = [&](auto begin_it, auto end_it)
         {
             for (auto it = begin_it; it != end_it; ++it)
             {
                 projected_polygon.outer().emplace_back(
                     project_point_to_coordinate_plane(*it, proj1, proj2));
             }
         };
 
         if (dir > 0) {
             add_projected_points(points_begin(ring), points_end(ring));
         } else {
             add_projected_points(points_rbegin(ring), points_rend(ring));
         }
 
         // Check if the projected polygon is valid
         if (!resolve_dynamic::is_valid<polygon_2d>::apply(projected_polygon, visitor, strategy))
         {
             return intersection_type::invalid;
         }
 
         // Cache the projection for future use
         m_polygon_projections[polygon_index] = polygon_projection{std::move(projected_polygon),
             p1_idx, p2_idx, p3_idx, proj1, proj2, dir};
 
         return intersection_type::valid;
     }
 
 private:
     std::vector<polygon_projection> m_polygon_projections;
 };
 
 template <typename PolyhedralSurface>
 class is_valid_polyhedral_surface
 {
 public:
     template <typename VisitPolicy, typename Strategy>
     static inline bool apply(PolyhedralSurface const& surface,
                              VisitPolicy& visitor,
                              Strategy const& strategy)
     {
         boost::ignore_unused(strategy);
 
         valid_intersection<PolyhedralSurface> valid_intersection(surface.size());
 
         // Check for colinear/coplanar points and compute projections
         for (std::size_t i = 0; i < surface.size(); ++i)
         {
             auto intersection = valid_intersection.initialize(surface[i], i, visitor, strategy);
 
             using intersection_type = decltype(intersection);
 
             if (intersection == intersection_type::invalid)
             {
                 return false;
             }
         }
 
         // Create an rtree to use it for polygon pair intersections
         using box_3d = typename boost::geometry::detail::is_valid::valid_intersection
             <PolyhedralSurface>::box_3d;
         using value = std::pair<box_3d, std::size_t>;
 
         boost::geometry::index::rtree<value, boost::geometry::index::rstar<4>> rtree;
         std::vector<std::unique_ptr<box_3d>> boxes(surface.size());
 
         for (std::size_t i = 0; i < surface.size(); ++i)
         {
             auto b = boost::geometry::return_envelope<box_3d>(surface[i], strategy);
             boxes[i] = std::make_unique<box_3d>(b);
             rtree.insert(std::make_pair(b, i)); // use balancing algorithm
         }
 
         // We check pairs of polygons for having a valid intersection.
         // We iterate in a BFS manner to check connectivity of the surface.
 
         auto it1 = boost::begin(surface);
 
         polygon_frontier_type polygon_frontier(surface.size());
         polygon_frontier.push(0);
 
         std::vector<std::size_t> remaining_polygons(surface.size());
         std::vector<std::size_t> positions(surface.size());
 
         for (std::size_t i = 0; i < surface.size(); ++i)
         {
             remaining_polygons[i] = i;
             positions[i] = i;
         }
 
         int remaining_count = surface.size();
 
         while (!polygon_frontier.empty())
         {
             auto index1 = polygon_frontier.top();
             polygon_frontier.pop();
 
             it1 = boost::begin(surface) + index1;
 
             // Remove index from remaining polygons and update positions
             auto moving_index = remaining_polygons[remaining_count - 1];
             remaining_polygons[positions[index1]] = moving_index;
             positions[moving_index] = positions[index1];
 
             // Query all polygons that intersect with the current polygon.
             // We use a buffer to avoid precision issues with the rtree query.
             std::vector<value> intersecting_boxes;
             auto const& b_query = *boxes[index1];
             auto const& b_query_buf = boost::geometry::return_buffer<box_3d>(b_query, 0.0001);
 
             rtree.query(boost::geometry::index::intersects(b_query_buf),
                         std::back_inserter(intersecting_boxes));
 
             for (int i = 0; i < remaining_count - 1; ++i)
             {
                 auto index2 = remaining_polygons[i];
 
                 // Skip if index2 is not in the intersecting boxes
                 if (std::find_if(intersecting_boxes.begin(), intersecting_boxes.end(),
                     [index2](const value& v) { return v.second == index2; }
                 ) == intersecting_boxes.end())
                 {
                     continue;
                 }
 
                 auto intersection = valid_intersection.apply(*it1, surface[index2], index1, index2,
                     visitor, strategy);
 
                 using intersection_type = decltype(intersection);
 
                 if (intersection == intersection_type::invalid)
                 {
                     return false; //  Invalid intersection found (including inconsistent orientation)
                 }
                 else if (intersection == intersection_type::valid)
                 {
                     polygon_frontier.push(index2);
                 }
             }
             remaining_count--;
         }
 
         if (remaining_count > 0)
         {
             visitor.template apply<failure_disconnected_surface>();
             return false; // Not all polygons were visited, hence disconnected
         }
 
         return true;
     }
 };
 
 }} // namespace detail::is_valid
 #endif // DOXYGEN_NO_DETAIL
 
 #ifndef DOXYGEN_NO_DISPATCH
 namespace dispatch
 {
 
 
 // Reference (for validity of Polyhedral Surfaces): OGC
 // OpenGIS® Implementation Standard for Geographic
 // information - Simple feature access - Part 1: Common
 // architecture: 6.1.12
 template <typename PolyhedralSurface, bool AllowEmptyMultiGeometries>
 struct is_valid
     <
         PolyhedralSurface, polyhedral_surface_tag, AllowEmptyMultiGeometries
     > : detail::is_valid::is_valid_polyhedral_surface<PolyhedralSurface>
 {};
 
 
 } // namespace dispatch
 #endif // DOXYGEN_NO_DISPATCH
 
 }} // namespace boost::geometry
 
 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_IS_VALID_POLYHEDRAL_SURFACE_HPP

This page was automatically generated by the 2.3.7 LXR engine. The LXR team