detail/within/multi_point.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2017-2023, Oracle and/or its affiliates.
0006
0007 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0008 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0009
0010 // Use, modification and distribution is subject to the Boost Software License,
0011 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0012 // http://www.boost.org/LICENSE_1_0.txt)
0013
0014 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_WITHIN_MULTI_POINT_HPP
0015 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_WITHIN_MULTI_POINT_HPP
0016
0017
0018 #include <algorithm>
0019 #include <vector>
0020
0021 #include <boost/range/begin.hpp>
0022 #include <boost/range/end.hpp>
0023 #include <boost/range/size.hpp>
0024 #include <boost/range/value_type.hpp>
0025
0026 #include <boost/geometry/algorithms/detail/disjoint/box_box.hpp>
0027 #include <boost/geometry/algorithms/detail/disjoint/point_box.hpp>
0028 #include <boost/geometry/algorithms/detail/expand_by_epsilon.hpp>
0029 #include <boost/geometry/algorithms/detail/within/point_in_geometry.hpp>
0030 #include <boost/geometry/algorithms/envelope.hpp>
0031 #include <boost/geometry/algorithms/detail/partition.hpp>
0032 #include <boost/geometry/core/tag.hpp>
0033 #include <boost/geometry/core/tags.hpp>
0034
0035 #include <boost/geometry/geometries/box.hpp>
0036
0037 #include <boost/geometry/index/rtree.hpp>
0038
0039 #include <boost/geometry/policies/compare.hpp>
0040
0041 #include <boost/geometry/strategies/covered_by.hpp>
0042 #include <boost/geometry/strategies/disjoint.hpp>
0043
0044 #include <boost/geometry/util/constexpr.hpp>
0045 #include <boost/geometry/util/type_traits.hpp>
0046
0047
0048 namespace boost { namespace geometry {
0049
0050 #ifndef DOXYGEN_NO_DETAIL
0051 namespace detail { namespace within {
0052
0053 struct multi_point_point
0054 {
0055     template <typename MultiPoint, typename Point, typename Strategy>
0056     static inline bool apply(MultiPoint const& multi_point,
0057                              Point const& point,
0058                              Strategy const& strategy)
0059     {
0060         auto const s = strategy.relate(multi_point, point);
0061
0062         for (auto it = boost::begin(multi_point); it != boost::end(multi_point); ++it)
0063         {
0064             if (! s.apply(*it, point))
0065             {
0066                 return false;
0067             }
0068         }
0069
0070         // all points of MultiPoint inside Point
0071         return true;
0072     }
0073 };
0074
0075 // NOTE: currently the strategy is ignored, math::equals() is used inside geometry::less<>
0076 struct multi_point_multi_point
0077 {
0078     template <typename MultiPoint1, typename MultiPoint2, typename Strategy>
0079     static inline bool apply(MultiPoint1 const& multi_point1,
0080                              MultiPoint2 const& multi_point2,
0081                              Strategy const& /*strategy*/)
0082     {
0083         typedef typename boost::range_value<MultiPoint2>::type point2_type;
0084         typedef geometry::less<void, -1, Strategy> less_type;
0085
0086         less_type const less = less_type();
0087
0088         std::vector<point2_type> points2(boost::begin(multi_point2), boost::end(multi_point2));
0089         std::sort(points2.begin(), points2.end(), less);
0090
0091         bool result = false;
0092
0093         for (auto it = boost::begin(multi_point1); it != boost::end(multi_point1); ++it)
0094         {
0095             if (! std::binary_search(points2.begin(), points2.end(), *it, less))
0096             {
0097                 return false;
0098             }
0099             else
0100             {
0101                 result = true;
0102             }
0103         }
0104
0105         return result;
0106     }
0107 };
0108
0109
0110 // TODO: the complexity could be lesser
0111 //   the second geometry could be "prepared"/sorted
0112 // For Linear geometries partition could be used
0113 // For Areal geometries point_in_geometry() would have to call the winding
0114 //   strategy differently, currently it linearly calls the strategy for each
0115 //   segment. So the segments would have to be sorted in a way consistent with
0116 //   the strategy and then the strategy called only for the segments in range.
0117 template <bool Within>
0118 struct multi_point_single_geometry
0119 {
0120     template <typename MultiPoint, typename LinearOrAreal, typename Strategy>
0121     static inline bool apply(MultiPoint const& multi_point,
0122                              LinearOrAreal const& linear_or_areal,
0123                              Strategy const& strategy)
0124     {
0125         // Create envelope of geometry
0126         model::box<point_type_t<LinearOrAreal>> box;
0127         geometry::envelope(linear_or_areal, box, strategy);
0128         geometry::detail::expand_by_epsilon(box);
0129
0130         // Test each Point with envelope and then geometry if needed
0131         // If in the exterior, break
0132         bool result = false;
0133
0134         for (auto it = boost::begin(multi_point); it != boost::end(multi_point); ++it )
0135         {
0136             typedef decltype(strategy.covered_by(*it, box)) point_in_box_type;
0137
0138             int in_val = 0;
0139
0140             // exterior of box and of geometry
0141             if (! point_in_box_type::apply(*it, box)
0142                 || (in_val = point_in_geometry(*it, linear_or_areal, strategy)) < 0)
0143             {
0144                 result = false;
0145                 break;
0146             }
0147
0148             // interior : interior/boundary
0149             if (Within ? in_val > 0 : in_val >= 0)
0150             {
0151                 result = true;
0152             }
0153         }
0154
0155         return result;
0156     }
0157 };
0158
0159
0160 // TODO: same here, probably the complexity could be lesser
0161 template <bool Within>
0162 struct multi_point_multi_geometry
0163 {
0164     template <typename MultiPoint, typename LinearOrAreal, typename Strategy>
0165     static inline bool apply(MultiPoint const& multi_point,
0166                              LinearOrAreal const& linear_or_areal,
0167                              Strategy const& strategy)
0168     {
0169         using point2_type = point_type_t<LinearOrAreal>;
0170         using box2_type = model::box<point2_type>;
0171         using box_pair_type = std::pair<box2_type, std::size_t>;
0172         using box_pair_vector = std::vector<box_pair_type>;
0173
0174         static const bool is_linear = util::is_linear<LinearOrAreal>::value;
0175
0176         // TODO: box pairs could be constructed on the fly, inside the rtree
0177
0178         // Prepare range of envelopes and ids
0179         std::size_t count2 = boost::size(linear_or_areal);
0180         box_pair_vector boxes(count2);
0181         for (std::size_t i = 0 ; i < count2 ; ++i)
0182         {
0183             geometry::envelope(linear_or_areal, boxes[i].first, strategy);
0184             geometry::detail::expand_by_epsilon(boxes[i].first);
0185             boxes[i].second = i;
0186         }
0187
0188         // Create R-tree
0189         typedef index::parameters<index::rstar<4>, Strategy> index_parameters_type;
0190         index::rtree<box_pair_type, index_parameters_type>
0191             rtree(boxes.begin(), boxes.end(),
0192                   index_parameters_type(index::rstar<4>(), strategy));
0193
0194         // For each point find overlapping envelopes and test corresponding single geometries
0195         // If a point is in the exterior break
0196         bool result = false;
0197
0198         for (auto it = boost::begin(multi_point); it != boost::end(multi_point); ++it)
0199         {
0200             // TODO: investigate the possibility of using satisfies
0201             // TODO: investigate the possibility of using iterative queries (optimization below)
0202             box_pair_vector inters_boxes;
0203             rtree.query(index::intersects(*it), std::back_inserter(inters_boxes));
0204
0205             bool found_interior = false;
0206             bool found_boundary = false;
0207             int boundaries = 0;
0208
0209             typedef typename box_pair_vector::const_iterator box_iterator;
0210             for (box_iterator box_it = inters_boxes.begin() ;
0211                  box_it != inters_boxes.end() ; ++box_it )
0212             {
0213                 int const in_val = point_in_geometry(*it,
0214                     range::at(linear_or_areal, box_it->second), strategy);
0215
0216                 if (in_val > 0)
0217                 {
0218                     found_interior = true;
0219                 }
0220                 else if (in_val == 0)
0221                 {
0222                     ++boundaries;
0223                 }
0224
0225                 // If the result was set previously (interior or
0226                 // interior/boundary found) the only thing that needs to be
0227                 // done for other points is to make sure they're not
0228                 // overlapping the exterior no need to analyse boundaries.
0229                 if (result && in_val >= 0)
0230                 {
0231                     break;
0232                 }
0233             }
0234
0235             if (boundaries > 0)
0236             {
0237                 if BOOST_GEOMETRY_CONSTEXPR (is_linear)
0238                 {
0239                     if (boundaries % 2 == 0)
0240                     {
0241                         found_interior = true;
0242                     }
0243                     else
0244                     {
0245                         found_boundary = true;
0246                     }
0247                 }
0248                 else
0249                 {
0250                     found_boundary = true;
0251                 }
0252             }
0253
0254             // exterior
0255             if (! found_interior && ! found_boundary)
0256             {
0257                 result = false;
0258                 break;
0259             }
0260
0261             // interior : interior/boundary
0262             if (Within ? found_interior : (found_interior || found_boundary))
0263             {
0264                 result = true;
0265             }
0266         }
0267
0268         return result;
0269     }
0270 };
0271
0272 }} // namespace detail::within
0273 #endif // DOXYGEN_NO_DETAIL
0274
0275 }} // namespace boost::geometry
0276
0277
0278 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_WITHIN_MULTI_POINT_HPP