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