1 // Boost.Geometry
2 
3 // Copyright (c) 2017 Oracle and/or its affiliates.
4 
5 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
6 
7 // Use, modification and distribution is subject to the Boost Software License,
8 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
9 // http://www.boost.org/LICENSE_1_0.txt)
10 
11 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_WITHIN_MULTI_POINT_HPP
12 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_WITHIN_MULTI_POINT_HPP
13 
14 
15 #include <algorithm>
16 #include <vector>
17 
18 #include <boost/range.hpp>
19 #include <boost/type_traits/is_same.hpp>
20 
21 #include <boost/geometry/algorithms/detail/disjoint/box_box.hpp>
22 #include <boost/geometry/algorithms/detail/disjoint/point_box.hpp>
23 #include <boost/geometry/algorithms/detail/expand_by_epsilon.hpp>
24 #include <boost/geometry/algorithms/detail/within/point_in_geometry.hpp>
25 #include <boost/geometry/algorithms/envelope.hpp>
26 #include <boost/geometry/algorithms/detail/partition.hpp>
27 #include <boost/geometry/core/tag.hpp>
28 #include <boost/geometry/core/tag_cast.hpp>
29 #include <boost/geometry/core/tags.hpp>
30 
31 #include <boost/geometry/geometries/box.hpp>
32 
33 #include <boost/geometry/index/rtree.hpp>
34 
35 #include <boost/geometry/policies/compare.hpp>
36 
37 #include <boost/geometry/strategies/covered_by.hpp>
38 #include <boost/geometry/strategies/disjoint.hpp>
39 
40 
41 namespace boost { namespace geometry {
42 
43 #ifndef DOXYGEN_NO_DETAIL
44 namespace detail { namespace within {
45 
46 struct multi_point_point
47 {
48     template <typename MultiPoint, typename Point, typename Strategy>
49     static inline bool apply(MultiPoint const& multi_point,
50                              Point const& point,
51                              Strategy const& strategy)
52     {
53         typedef typename boost::range_const_iterator<MultiPoint>::type iterator;
54         for ( iterator it = boost::begin(multi_point) ; it != boost::end(multi_point) ; ++it )
55         {
56             if (! strategy.apply(*it, point))
57             {
58                 return false;
59             }
60         }
61 
62         // all points of MultiPoint inside Point
63         return true;
64     }
65 };
66 
67 // NOTE: currently the strategy is ignored, math::equals() is used inside geometry::less<>
68 struct multi_point_multi_point
69 {
70     template <typename MultiPoint1, typename MultiPoint2, typename Strategy>
71     static inline bool apply(MultiPoint1 const& multi_point1,
72                              MultiPoint2 const& multi_point2,
73                              Strategy const& /*strategy*/)
74     {
75         typedef typename boost::range_value<MultiPoint2>::type point2_type;
76 
77         geometry::less<> const less = geometry::less<>();
78 
79         std::vector<point2_type> points2(boost::begin(multi_point2), boost::end(multi_point2));
80         std::sort(points2.begin(), points2.end(), less);
81 
82         bool result = false;
83 
84         typedef typename boost::range_const_iterator<MultiPoint1>::type iterator;
85         for ( iterator it = boost::begin(multi_point1) ; it != boost::end(multi_point1) ; ++it )
86         {
87             if (! std::binary_search(points2.begin(), points2.end(), *it, less))
88             {
89                 return false;
90             }
91             else
92             {
93                 result = true;
94             }
95         }
96 
97         return result;
98     }
99 };
100 
101 
102 // TODO: the complexity could be lesser
103 //   the second geometry could be "prepared"/sorted
104 // For Linear geometries partition could be used
105 // For Areal geometries point_in_geometry() would have to call the winding
106 //   strategy differently, currently it linearly calls the strategy for each
107 //   segment. So the segments would have to be sorted in a way consistent with
108 //   the strategy and then the strategy called only for the segments in range.
109 template <bool Within>
110 struct multi_point_single_geometry
111 {
112     template <typename MultiPoint, typename LinearOrAreal, typename Strategy>
113     static inline bool apply(MultiPoint const& multi_point,
114                              LinearOrAreal const& linear_or_areal,
115                              Strategy const& strategy)
116     {
117         typedef typename boost::range_value<MultiPoint>::type point1_type;
118         typedef typename point_type<LinearOrAreal>::type point2_type;
119         typedef model::box<point2_type> box2_type;
120 
121         // Create envelope of geometry
122         box2_type box;
123         geometry::envelope(linear_or_areal, box, strategy.get_envelope_strategy());
124         geometry::detail::expand_by_epsilon(box);
125 
126         typedef typename strategy::covered_by::services::default_strategy
127             <
128                 point1_type, box2_type
129             >::type point_in_box_type;
130 
131         // Test each Point with envelope and then geometry if needed
132         // If in the exterior, break
133         bool result = false;
134 
135         typedef typename boost::range_const_iterator<MultiPoint>::type iterator;
136         for ( iterator it = boost::begin(multi_point) ; it != boost::end(multi_point) ; ++it )
137         {
138             int in_val = 0;
139 
140             // exterior of box and of geometry
141             if (! point_in_box_type::apply(*it, box)
142                 || (in_val = point_in_geometry(*it, linear_or_areal, strategy)) < 0)
143             {
144                 result = false;
145                 break;
146             }
147 
148             // interior : interior/boundary
149             if (Within ? in_val > 0 : in_val >= 0)
150             {
151                 result = true;
152             }
153         }
154 
155         return result;
156     }
157 };
158 
159 
160 // TODO: same here, probably the complexity could be lesser
161 template <bool Within>
162 struct multi_point_multi_geometry
163 {
164     template <typename MultiPoint, typename LinearOrAreal, typename Strategy>
165     static inline bool apply(MultiPoint const& multi_point,
166                              LinearOrAreal const& linear_or_areal,
167                              Strategy const& strategy)
168     {
169         typedef typename point_type<LinearOrAreal>::type point2_type;
170         typedef model::box<point2_type> box2_type;
171         static const bool is_linear = is_same
172             <
173                 typename tag_cast
174                     <
175                         typename tag<LinearOrAreal>::type,
176                         linear_tag
177                     >::type,
178                 linear_tag
179             >::value;
180 
181         typename Strategy::envelope_strategy_type const
182             envelope_strategy = strategy.get_envelope_strategy();
183 
184         // TODO: box pairs could be constructed on the fly, inside the rtree
185 
186         // Prepare range of envelopes and ids
187         std::size_t count2 = boost::size(linear_or_areal);
188         typedef std::pair<box2_type, std::size_t> box_pair_type;
189         typedef std::vector<box_pair_type> box_pair_vector;
190         box_pair_vector boxes(count2);
191         for (std::size_t i = 0 ; i < count2 ; ++i)
192         {
193             geometry::envelope(linear_or_areal, boxes[i].first, envelope_strategy);
194             geometry::detail::expand_by_epsilon(boxes[i].first);
195             boxes[i].second = i;
196         }
197 
198         // Create R-tree
199         index::rtree<box_pair_type, index::rstar<4> > rtree(boxes.begin(), boxes.end());
200 
201         // For each point find overlapping envelopes and test corresponding single geometries
202         // If a point is in the exterior break
203         bool result = false;
204 
205         typedef typename boost::range_const_iterator<MultiPoint>::type iterator;
206         for ( iterator it = boost::begin(multi_point) ; it != boost::end(multi_point) ; ++it )
207         {
208             // TODO: investigate the possibility of using satisfies
209             // TODO: investigate the possibility of using iterative queries (optimization below)
210             box_pair_vector inters_boxes;
211             rtree.query(index::intersects(*it), std::back_inserter(inters_boxes));
212 
213             bool found_interior = false;
214             bool found_boundary = false;
215             int boundaries = 0;
216 
217             typedef typename box_pair_vector::const_iterator iterator;
218             for ( iterator box_it = inters_boxes.begin() ; box_it != inters_boxes.end() ; ++box_it )
219             {
220                 int in_val = point_in_geometry(*it, range::at(linear_or_areal, box_it->second), strategy);
221 
222                 if (in_val > 0)
223                     found_interior = true;
224                 else if (in_val == 0)
225                     ++boundaries;
226 
227                 // If the result was set previously (interior or
228                 // interior/boundary found) the only thing that needs to be
229                 // done for other points is to make sure they're not
230                 // overlapping the exterior no need to analyse boundaries.
231                 if (result && in_val >= 0)
232                 {
233                     break;
234                 }
235             }
236 
237             if ( boundaries > 0)
238             {
239                 if (is_linear && boundaries % 2 == 0)
240                     found_interior = true;
241                 else
242                     found_boundary = true;
243             }
244 
245             // exterior
246             if (! found_interior && ! found_boundary)
247             {
248                 result = false;
249                 break;
250             }
251 
252             // interior : interior/boundary
253             if (Within ? found_interior : (found_interior || found_boundary))
254             {
255                 result = true;
256             }
257         }
258 
259         return result;
260     }
261 };
262 
263 }} // namespace detail::within
264 #endif // DOXYGEN_NO_DETAIL
265 
266 }} // namespace boost::geometry
267 
268 
269 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_WITHIN_MULTI_POINT_HPP
270