]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Boost.Geometry (aka GGL, Generic Geometry Library) |
2 | ||
3 | // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. | |
4 | // Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland. | |
5 | ||
6 | // This file was modified by Oracle on 2013, 2014, 2016. | |
7 | // Modifications copyright (c) 2013-2016 Oracle and/or its affiliates. | |
8 | // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle | |
9 | ||
10 | // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library | |
11 | // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. | |
12 | ||
13 | // Use, modification and distribution is subject to the Boost Software License, | |
14 | // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at | |
15 | // http://www.boost.org/LICENSE_1_0.txt) | |
16 | ||
17 | #ifndef BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP | |
18 | #define BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP | |
19 | ||
20 | ||
21 | #include <boost/core/ignore_unused.hpp> | |
22 | ||
23 | #include <boost/geometry/util/math.hpp> | |
24 | #include <boost/geometry/util/select_calculation_type.hpp> | |
25 | ||
26 | #include <boost/geometry/strategies/side.hpp> | |
27 | #include <boost/geometry/strategies/covered_by.hpp> | |
28 | #include <boost/geometry/strategies/within.hpp> | |
29 | ||
30 | ||
31 | namespace boost { namespace geometry | |
32 | { | |
33 | ||
34 | namespace strategy { namespace within | |
35 | { | |
36 | ||
37 | // 1 deg or pi/180 rad | |
38 | template <typename Point, | |
39 | typename CalculationType = typename coordinate_type<Point>::type> | |
40 | struct winding_small_angle | |
41 | { | |
42 | typedef typename coordinate_system<Point>::type cs_t; | |
43 | typedef math::detail::constants_on_spheroid | |
44 | < | |
45 | CalculationType, | |
46 | typename cs_t::units | |
47 | > constants; | |
48 | ||
49 | static inline CalculationType apply() | |
50 | { | |
51 | return constants::half_period() / CalculationType(180); | |
52 | } | |
53 | }; | |
54 | ||
55 | ||
56 | // Fix for https://svn.boost.org/trac/boost/ticket/9628 | |
57 | // For floating point coordinates, the <D> coordinate of a point is compared | |
58 | // with the segment's points using some EPS. If the coordinates are "equal" | |
59 | // the sides are calculated. Therefore we can treat a segment as a long areal | |
60 | // geometry having some width. There is a small ~triangular area somewhere | |
61 | // between the segment's effective area and a segment's line used in sides | |
62 | // calculation where the segment is on the one side of the line but on the | |
63 | // other side of a segment (due to the width). | |
64 | // Below picture assuming D = 1, if D = 0 horiz<->vert, E<->N, RIGHT<->UP. | |
65 | // For the s1 of a segment going NE the real side is RIGHT but the point may | |
66 | // be detected as LEFT, like this: | |
67 | // RIGHT | |
68 | // ___-----> | |
69 | // ^ O Pt __ __ | |
70 | // EPS __ __ | |
71 | // v__ __ BUT DETECTED AS LEFT OF THIS LINE | |
72 | // _____7 | |
73 | // _____/ | |
74 | // _____/ | |
75 | // In the code below actually D = 0, so segments are nearly-vertical | |
76 | // Called when the point is on the same level as one of the segment's points | |
77 | // but the point is not aligned with a vertical segment | |
78 | template <typename CSTag> | |
79 | struct winding_side_equal | |
80 | { | |
81 | typedef typename strategy::side::services::default_strategy | |
82 | < | |
83 | CSTag | |
84 | >::type strategy_side_type; | |
85 | ||
86 | template <typename Point, typename PointOfSegment> | |
87 | static inline int apply(Point const& point, | |
88 | PointOfSegment const& se, | |
89 | int count) | |
90 | { | |
91 | typedef typename coordinate_type<PointOfSegment>::type scoord_t; | |
92 | typedef typename coordinate_system<PointOfSegment>::type::units units_t; | |
93 | ||
94 | if (math::equals(get<1>(point), get<1>(se))) | |
95 | return 0; | |
96 | ||
97 | // Create a horizontal segment intersecting the original segment's endpoint | |
98 | // equal to the point, with the derived direction (E/W). | |
99 | PointOfSegment ss1, ss2; | |
100 | set<1>(ss1, get<1>(se)); | |
101 | set<0>(ss1, get<0>(se)); | |
102 | set<1>(ss2, get<1>(se)); | |
103 | scoord_t ss20 = get<0>(se); | |
104 | if (count > 0) | |
105 | { | |
106 | ss20 += winding_small_angle<PointOfSegment>::apply(); | |
107 | } | |
108 | else | |
109 | { | |
110 | ss20 -= winding_small_angle<PointOfSegment>::apply(); | |
111 | } | |
112 | math::normalize_longitude<units_t>(ss20); | |
113 | set<0>(ss2, ss20); | |
114 | ||
115 | // Check the side using this vertical segment | |
116 | return strategy_side_type::apply(ss1, ss2, point); | |
117 | } | |
118 | }; | |
119 | // The optimization for cartesian | |
120 | template <> | |
121 | struct winding_side_equal<cartesian_tag> | |
122 | { | |
123 | template <typename Point, typename PointOfSegment> | |
124 | static inline int apply(Point const& point, | |
125 | PointOfSegment const& se, | |
126 | int count) | |
127 | { | |
128 | // NOTE: for D=0 the signs would be reversed | |
129 | return math::equals(get<1>(point), get<1>(se)) ? | |
130 | 0 : | |
131 | get<1>(point) < get<1>(se) ? | |
132 | // assuming count is equal to 1 or -1 | |
133 | -count : // ( count > 0 ? -1 : 1) : | |
134 | count; // ( count > 0 ? 1 : -1) ; | |
135 | } | |
136 | }; | |
137 | ||
138 | ||
139 | template <typename Point, | |
140 | typename CalculationType, | |
141 | typename CSTag = typename cs_tag<Point>::type> | |
142 | struct winding_check_touch | |
143 | { | |
144 | typedef CalculationType calc_t; | |
145 | typedef typename coordinate_system<Point>::type::units units_t; | |
146 | typedef math::detail::constants_on_spheroid<CalculationType, units_t> constants; | |
147 | ||
148 | template <typename PointOfSegment, typename State> | |
149 | static inline int apply(Point const& point, | |
150 | PointOfSegment const& seg1, | |
151 | PointOfSegment const& seg2, | |
152 | State& state, | |
153 | bool& eq1, | |
154 | bool& eq2) | |
155 | { | |
156 | calc_t const pi = constants::half_period(); | |
157 | calc_t const pi2 = pi / calc_t(2); | |
158 | ||
159 | calc_t const px = get<0>(point); | |
160 | calc_t const s1x = get<0>(seg1); | |
161 | calc_t const s2x = get<0>(seg2); | |
162 | calc_t const py = get<1>(point); | |
163 | calc_t const s1y = get<1>(seg1); | |
164 | calc_t const s2y = get<1>(seg2); | |
165 | ||
166 | // NOTE: lat in {-90, 90} and arbitrary lon | |
167 | // it doesn't matter what lon it is if it's a pole | |
168 | // so e.g. if one of the segment endpoints is a pole | |
169 | // then only the other lon matters | |
170 | ||
171 | bool eq1_strict = math::equals(s1x, px); | |
172 | bool eq2_strict = math::equals(s2x, px); | |
173 | ||
174 | eq1 = eq1_strict // lon strictly equal to s1 | |
175 | || math::equals(s1y, pi2) || math::equals(s1y, -pi2); // s1 is pole | |
176 | eq2 = eq2_strict // lon strictly equal to s2 | |
177 | || math::equals(s2y, pi2) || math::equals(s2y, -pi2); // s2 is pole | |
178 | ||
179 | // segment overlapping pole | |
180 | calc_t s1x_anti = s1x + constants::half_period(); | |
181 | math::normalize_longitude<units_t, calc_t>(s1x_anti); | |
182 | bool antipodal = math::equals(s2x, s1x_anti); | |
183 | if (antipodal) | |
184 | { | |
185 | eq1 = eq2 = eq1 || eq2; | |
186 | ||
187 | // segment overlapping pole and point is pole | |
188 | if (math::equals(py, pi2) || math::equals(py, -pi2)) | |
189 | { | |
190 | eq1 = eq2 = true; | |
191 | } | |
192 | } | |
193 | ||
194 | // Both equal p -> segment vertical | |
195 | // The only thing which has to be done is check if point is ON segment | |
196 | if (eq1 && eq2) | |
197 | { | |
198 | // segment endpoints on the same sides of the globe | |
199 | if (! antipodal | |
200 | // p's lat between segment endpoints' lats | |
201 | ? (s1y <= py && s2y >= py) || (s2y <= py && s1y >= py) | |
202 | // going through north or south pole? | |
203 | : (pi - s1y - s2y <= pi | |
204 | ? (eq1_strict && s1y <= py) || (eq2_strict && s2y <= py) // north | |
205 | || math::equals(py, pi2) // point on north pole | |
206 | : (eq1_strict && s1y >= py) || (eq2_strict && s2y >= py)) // south | |
207 | || math::equals(py, -pi2) // point on south pole | |
208 | ) | |
209 | { | |
210 | state.m_touches = true; | |
211 | } | |
212 | return true; | |
213 | } | |
214 | return false; | |
215 | } | |
216 | }; | |
217 | // The optimization for cartesian | |
218 | template <typename Point, typename CalculationType> | |
219 | struct winding_check_touch<Point, CalculationType, cartesian_tag> | |
220 | { | |
221 | typedef CalculationType calc_t; | |
222 | ||
223 | template <typename PointOfSegment, typename State> | |
224 | static inline bool apply(Point const& point, | |
225 | PointOfSegment const& seg1, | |
226 | PointOfSegment const& seg2, | |
227 | State& state, | |
228 | bool& eq1, | |
229 | bool& eq2) | |
230 | { | |
231 | calc_t const px = get<0>(point); | |
232 | calc_t const s1x = get<0>(seg1); | |
233 | calc_t const s2x = get<0>(seg2); | |
234 | ||
235 | eq1 = math::equals(s1x, px); | |
236 | eq2 = math::equals(s2x, px); | |
237 | ||
238 | // Both equal p -> segment vertical | |
239 | // The only thing which has to be done is check if point is ON segment | |
240 | if (eq1 && eq2) | |
241 | { | |
242 | calc_t const py = get<1>(point); | |
243 | calc_t const s1y = get<1>(seg1); | |
244 | calc_t const s2y = get<1>(seg2); | |
245 | if ((s1y <= py && s2y >= py) || (s2y <= py && s1y >= py)) | |
246 | { | |
247 | state.m_touches = true; | |
248 | } | |
249 | return true; | |
250 | } | |
251 | return false; | |
252 | } | |
253 | }; | |
254 | ||
255 | ||
256 | // Called if point is not aligned with a vertical segment | |
257 | template <typename Point, | |
258 | typename CalculationType, | |
259 | typename CSTag = typename cs_tag<Point>::type> | |
260 | struct winding_calculate_count | |
261 | { | |
262 | typedef CalculationType calc_t; | |
263 | typedef typename coordinate_system<Point>::type::units units_t; | |
264 | ||
265 | static inline bool greater(calc_t const& l, calc_t const& r) | |
266 | { | |
267 | calc_t diff = l - r; | |
268 | math::normalize_longitude<units_t, calc_t>(diff); | |
269 | return diff > calc_t(0); | |
270 | } | |
271 | ||
272 | static inline int apply(calc_t const& p, | |
273 | calc_t const& s1, calc_t const& s2, | |
274 | bool eq1, bool eq2) | |
275 | { | |
276 | // Probably could be optimized by avoiding normalization for some comparisons | |
277 | // e.g. s1 > p could be calculated from p > s1 | |
278 | ||
279 | // If both segment endpoints were poles below checks wouldn't be enough | |
280 | // but this means that either both are the same or that they are N/S poles | |
281 | // and therefore the segment is not valid. | |
282 | // If needed (eq1 && eq2 ? 0) could be returned | |
283 | ||
284 | return | |
285 | eq1 ? (greater(s2, p) ? 1 : -1) // Point on level s1, E/W depending on s2 | |
286 | : eq2 ? (greater(s1, p) ? -1 : 1) // idem | |
287 | : greater(p, s1) && greater(s2, p) ? 2 // Point between s1 -> s2 --> E | |
288 | : greater(p, s2) && greater(s1, p) ? -2 // Point between s2 -> s1 --> W | |
289 | : 0; | |
290 | } | |
291 | }; | |
292 | // The optimization for cartesian | |
293 | template <typename Point, typename CalculationType> | |
294 | struct winding_calculate_count<Point, CalculationType, cartesian_tag> | |
295 | { | |
296 | typedef CalculationType calc_t; | |
297 | ||
298 | static inline int apply(calc_t const& p, | |
299 | calc_t const& s1, calc_t const& s2, | |
300 | bool eq1, bool eq2) | |
301 | { | |
302 | return | |
303 | eq1 ? (s2 > p ? 1 : -1) // Point on level s1, E/W depending on s2 | |
304 | : eq2 ? (s1 > p ? -1 : 1) // idem | |
305 | : s1 < p && s2 > p ? 2 // Point between s1 -> s2 --> E | |
306 | : s2 < p && s1 > p ? -2 // Point between s2 -> s1 --> W | |
307 | : 0; | |
308 | } | |
309 | }; | |
310 | ||
311 | ||
312 | /*! | |
313 | \brief Within detection using winding rule | |
314 | \ingroup strategies | |
315 | \tparam Point \tparam_point | |
316 | \tparam PointOfSegment \tparam_segment_point | |
317 | \tparam CalculationType \tparam_calculation | |
318 | \author Barend Gehrels | |
319 | \note The implementation is inspired by terralib http://www.terralib.org (LGPL) | |
320 | \note but totally revised afterwards, especially for cases on segments | |
321 | \note Only dependant on "side", -> agnostic, suitable for spherical/latlong | |
322 | ||
323 | \qbk{ | |
324 | [heading See also] | |
325 | [link geometry.reference.algorithms.within.within_3_with_strategy within (with strategy)] | |
326 | } | |
327 | */ | |
328 | template | |
329 | < | |
330 | typename Point, | |
331 | typename PointOfSegment = Point, | |
332 | typename CalculationType = void | |
333 | > | |
334 | class winding | |
335 | { | |
336 | typedef typename select_calculation_type | |
337 | < | |
338 | Point, | |
339 | PointOfSegment, | |
340 | CalculationType | |
341 | >::type calculation_type; | |
342 | ||
343 | ||
344 | typedef typename strategy::side::services::default_strategy | |
345 | < | |
346 | typename cs_tag<Point>::type | |
347 | >::type strategy_side_type; | |
348 | ||
349 | ||
350 | /*! subclass to keep state */ | |
351 | class counter | |
352 | { | |
353 | int m_count; | |
354 | bool m_touches; | |
355 | ||
356 | inline int code() const | |
357 | { | |
358 | return m_touches ? 0 : m_count == 0 ? -1 : 1; | |
359 | } | |
360 | ||
361 | public : | |
362 | friend class winding; | |
363 | ||
364 | template <typename P, typename CT, typename CST> | |
365 | friend struct winding_check_touch; | |
366 | ||
367 | inline counter() | |
368 | : m_count(0) | |
369 | , m_touches(false) | |
370 | {} | |
371 | ||
372 | }; | |
373 | ||
374 | ||
375 | static inline int check_segment(Point const& point, | |
376 | PointOfSegment const& seg1, PointOfSegment const& seg2, | |
377 | counter& state, bool& eq1, bool& eq2) | |
378 | { | |
379 | if (winding_check_touch<Point, calculation_type> | |
380 | ::apply(point, seg1, seg2, state, eq1, eq2)) | |
381 | { | |
382 | return 0; | |
383 | } | |
384 | ||
385 | calculation_type const p = get<0>(point); | |
386 | calculation_type const s1 = get<0>(seg1); | |
387 | calculation_type const s2 = get<0>(seg2); | |
388 | return winding_calculate_count<Point, calculation_type> | |
389 | ::apply(p, s1, s2, eq1, eq2); | |
390 | } | |
391 | ||
392 | ||
393 | public : | |
394 | ||
395 | // Typedefs and static methods to fulfill the concept | |
396 | typedef Point point_type; | |
397 | typedef PointOfSegment segment_point_type; | |
398 | typedef counter state_type; | |
399 | ||
400 | static inline bool apply(Point const& point, | |
401 | PointOfSegment const& s1, PointOfSegment const& s2, | |
402 | counter& state) | |
403 | { | |
404 | typedef typename cs_tag<Point>::type cs_t; | |
405 | ||
406 | bool eq1 = false; | |
407 | bool eq2 = false; | |
408 | boost::ignore_unused(eq2); | |
409 | ||
410 | int count = check_segment(point, s1, s2, state, eq1, eq2); | |
411 | if (count != 0) | |
412 | { | |
413 | int side = 0; | |
414 | if (count == 1 || count == -1) | |
415 | { | |
416 | side = winding_side_equal<cs_t>::apply(point, eq1 ? s1 : s2, count); | |
417 | } | |
418 | else // count == 2 || count == -2 | |
419 | { | |
420 | // 1 left, -1 right | |
421 | side = strategy_side_type::apply(s1, s2, point); | |
422 | } | |
423 | ||
424 | if (side == 0) | |
425 | { | |
426 | // Point is lying on segment | |
427 | state.m_touches = true; | |
428 | state.m_count = 0; | |
429 | return false; | |
430 | } | |
431 | ||
432 | // Side is NEG for right, POS for left. | |
433 | // The count is -2 for down, 2 for up (or -1/1) | |
434 | // Side positive thus means UP and LEFTSIDE or DOWN and RIGHTSIDE | |
435 | // See accompagnying figure (TODO) | |
436 | if (side * count > 0) | |
437 | { | |
438 | state.m_count += count; | |
439 | } | |
440 | } | |
441 | return ! state.m_touches; | |
442 | } | |
443 | ||
444 | static inline int result(counter const& state) | |
445 | { | |
446 | return state.code(); | |
447 | } | |
448 | }; | |
449 | ||
450 | ||
451 | #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS | |
452 | ||
453 | namespace services | |
454 | { | |
455 | ||
456 | // Register using "areal_tag" for ring, polygon, multi-polygon | |
457 | template <typename AnyTag, typename Point, typename Geometry> | |
458 | struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, cartesian_tag, cartesian_tag, Point, Geometry> | |
459 | { | |
460 | typedef winding<Point, typename geometry::point_type<Geometry>::type> type; | |
461 | }; | |
462 | ||
463 | template <typename AnyTag, typename Point, typename Geometry> | |
464 | struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, spherical_tag, spherical_tag, Point, Geometry> | |
465 | { | |
466 | typedef winding<Point, typename geometry::point_type<Geometry>::type> type; | |
467 | }; | |
468 | ||
469 | // TODO: use linear_tag and pointlike_tag the same way how areal_tag is used | |
470 | ||
471 | template <typename Point, typename Geometry, typename AnyTag> | |
472 | struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, cartesian_tag, cartesian_tag, Point, Geometry> | |
473 | { | |
474 | typedef winding<Point, typename geometry::point_type<Geometry>::type> type; | |
475 | }; | |
476 | ||
477 | template <typename Point, typename Geometry, typename AnyTag> | |
478 | struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, spherical_tag, spherical_tag, Point, Geometry> | |
479 | { | |
480 | typedef winding<Point, typename geometry::point_type<Geometry>::type> type; | |
481 | }; | |
482 | ||
483 | } // namespace services | |
484 | ||
485 | #endif | |
486 | ||
487 | ||
488 | }} // namespace strategy::within | |
489 | ||
490 | ||
491 | ||
492 | #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS | |
493 | namespace strategy { namespace covered_by { namespace services | |
494 | { | |
495 | ||
496 | // Register using "areal_tag" for ring, polygon, multi-polygon | |
497 | template <typename AnyTag, typename Point, typename Geometry> | |
498 | struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, cartesian_tag, cartesian_tag, Point, Geometry> | |
499 | { | |
500 | typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type; | |
501 | }; | |
502 | ||
503 | template <typename AnyTag, typename Point, typename Geometry> | |
504 | struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, spherical_tag, spherical_tag, Point, Geometry> | |
505 | { | |
506 | typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type; | |
507 | }; | |
508 | ||
509 | // TODO: use linear_tag and pointlike_tag the same way how areal_tag is used | |
510 | ||
511 | template <typename Point, typename Geometry, typename AnyTag> | |
512 | struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, cartesian_tag, cartesian_tag, Point, Geometry> | |
513 | { | |
514 | typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type; | |
515 | }; | |
516 | ||
517 | template <typename Point, typename Geometry, typename AnyTag> | |
518 | struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, spherical_tag, spherical_tag, Point, Geometry> | |
519 | { | |
520 | typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type; | |
521 | }; | |
522 | ||
523 | }}} // namespace strategy::covered_by::services | |
524 | #endif | |
525 | ||
526 | ||
527 | }} // namespace boost::geometry | |
528 | ||
529 | ||
530 | #endif // BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POLY_WINDING_HPP |