]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Boost.Geometry |
2 | ||
3 | // Copyright (c) 2016, Oracle and/or its affiliates. | |
4 | // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle | |
5 | ||
6 | // Use, modification and distribution is subject to the Boost Software License, | |
7 | // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at | |
8 | // http://www.boost.org/LICENSE_1_0.txt) | |
9 | ||
10 | #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_INTERSECTION_HPP | |
11 | #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_INTERSECTION_HPP | |
12 | ||
13 | #include <algorithm> | |
14 | ||
15 | #include <boost/geometry/core/cs.hpp> | |
16 | #include <boost/geometry/core/access.hpp> | |
17 | #include <boost/geometry/core/radian_access.hpp> | |
18 | #include <boost/geometry/core/tags.hpp> | |
19 | ||
20 | #include <boost/geometry/algorithms/detail/assign_values.hpp> | |
21 | #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp> | |
22 | #include <boost/geometry/algorithms/detail/equals/point_point.hpp> | |
23 | #include <boost/geometry/algorithms/detail/recalculate.hpp> | |
24 | ||
25 | #include <boost/geometry/arithmetic/arithmetic.hpp> | |
26 | #include <boost/geometry/arithmetic/cross_product.hpp> | |
27 | #include <boost/geometry/arithmetic/dot_product.hpp> | |
28 | #include <boost/geometry/formulas/spherical.hpp> | |
29 | ||
30 | #include <boost/geometry/geometries/concepts/point_concept.hpp> | |
31 | #include <boost/geometry/geometries/concepts/segment_concept.hpp> | |
32 | ||
33 | #include <boost/geometry/policies/robustness/segment_ratio.hpp> | |
34 | ||
35 | #include <boost/geometry/strategies/side_info.hpp> | |
36 | #include <boost/geometry/strategies/intersection.hpp> | |
37 | #include <boost/geometry/strategies/intersection_result.hpp> | |
38 | ||
39 | #include <boost/geometry/util/math.hpp> | |
40 | #include <boost/geometry/util/select_calculation_type.hpp> | |
41 | ||
42 | ||
43 | namespace boost { namespace geometry | |
44 | { | |
45 | ||
46 | namespace strategy { namespace intersection | |
47 | { | |
48 | ||
49 | // NOTE: | |
50 | // The coordinates of crossing IP may be calculated with small precision in some cases. | |
51 | // For double, near the equator noticed error ~1e-9 so far greater than | |
52 | // machine epsilon which is ~1e-16. This error is ~0.04m. | |
53 | // E.g. consider two cases, one near the origin and the second one rotated by 90 deg around Z or SN axis. | |
54 | // After the conversion from spherical degrees to cartesian 3d the following coordinates | |
55 | // are calculated: | |
56 | // for sph (-1 -1, 1 1) deg cart3d ys are -0.017449748351250485 and 0.017449748351250485 | |
57 | // for sph (89 -1, 91 1) deg cart3d xs are 0.017449748351250571 and -0.017449748351250450 | |
58 | // During the conversion degrees must first be converted to radians and then radians | |
59 | // are passed into trigonometric functions. The error may have several causes: | |
60 | // 1. Radians cannot represent exactly the same angles as degrees. | |
61 | // 2. Different longitudes are passed into sin() for x, corresponding to cos() for y, | |
62 | // and for different angle the error of the result may be different. | |
63 | // 3. These non-corresponding cartesian coordinates are used in calculation, | |
64 | // e.g. multiplied several times in cross and dot products. | |
65 | // If it was a problem this strategy could e.g. "normalize" longitudes before the conversion using the source units | |
66 | // by rotating the globe around Z axis, so moving longitudes always the same way towards the origin, | |
67 | // assuming this could help which is not clear. | |
68 | // For now, intersection points near the endpoints are checked explicitly if needed (if the IP is near the endpoint) | |
69 | // to generate precise result for them. Only the crossing (i) case may suffer from lower precision. | |
70 | ||
71 | template <typename Policy, typename CalculationType = void> | |
72 | struct relate_spherical_segments | |
73 | { | |
74 | typedef typename Policy::return_type return_type; | |
75 | ||
76 | enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 }; | |
77 | ||
78 | template <typename CoordinateType, typename SegmentRatio, typename Vector3d> | |
79 | struct segment_intersection_info | |
80 | { | |
81 | typedef typename select_most_precise | |
82 | < | |
83 | CoordinateType, double | |
84 | >::type promoted_type; | |
85 | ||
86 | promoted_type comparable_length_a() const | |
87 | { | |
88 | return robust_ra.denominator(); | |
89 | } | |
90 | ||
91 | promoted_type comparable_length_b() const | |
92 | { | |
93 | return robust_rb.denominator(); | |
94 | } | |
95 | ||
96 | template <typename Point, typename Segment1, typename Segment2> | |
97 | void assign_a(Point& point, Segment1 const& a, Segment2 const& b) const | |
98 | { | |
99 | assign(point, a, b); | |
100 | } | |
101 | template <typename Point, typename Segment1, typename Segment2> | |
102 | void assign_b(Point& point, Segment1 const& a, Segment2 const& b) const | |
103 | { | |
104 | assign(point, a, b); | |
105 | } | |
106 | ||
107 | template <typename Point, typename Segment1, typename Segment2> | |
108 | void assign(Point& point, Segment1 const& a, Segment2 const& b) const | |
109 | { | |
110 | if (ip_flag == ipi_inters) | |
111 | { | |
112 | // TODO: assign the rest of coordinates | |
113 | point = formula::cart3d_to_sph<Point>(intersection_point); | |
114 | } | |
115 | else if (ip_flag == ipi_at_a1) | |
116 | { | |
117 | detail::assign_point_from_index<0>(a, point); | |
118 | } | |
119 | else if (ip_flag == ipi_at_a2) | |
120 | { | |
121 | detail::assign_point_from_index<1>(a, point); | |
122 | } | |
123 | else if (ip_flag == ipi_at_b1) | |
124 | { | |
125 | detail::assign_point_from_index<0>(b, point); | |
126 | } | |
127 | else // ip_flag == ipi_at_b2 | |
128 | { | |
129 | detail::assign_point_from_index<1>(b, point); | |
130 | } | |
131 | } | |
132 | ||
133 | Vector3d intersection_point; | |
134 | SegmentRatio robust_ra; | |
135 | SegmentRatio robust_rb; | |
136 | intersection_point_flag ip_flag; | |
137 | }; | |
138 | ||
139 | // Relate segments a and b | |
140 | template <typename Segment1, typename Segment2, typename RobustPolicy> | |
141 | static inline return_type apply(Segment1 const& a, Segment2 const& b, | |
142 | RobustPolicy const& robust_policy) | |
143 | { | |
144 | typedef typename point_type<Segment1>::type point1_t; | |
145 | typedef typename point_type<Segment2>::type point2_t; | |
146 | point1_t a1, a2; | |
147 | point2_t b1, b2; | |
148 | ||
149 | // TODO: use indexed_point_view if possible? | |
150 | detail::assign_point_from_index<0>(a, a1); | |
151 | detail::assign_point_from_index<1>(a, a2); | |
152 | detail::assign_point_from_index<0>(b, b1); | |
153 | detail::assign_point_from_index<1>(b, b2); | |
154 | ||
155 | return apply(a, b, robust_policy, a1, a2, b1, b2); | |
156 | } | |
157 | ||
158 | // Relate segments a and b | |
159 | template <typename Segment1, typename Segment2, typename RobustPolicy, typename Point1, typename Point2> | |
160 | static inline return_type apply(Segment1 const& a, Segment2 const& b, | |
161 | RobustPolicy const&, | |
162 | Point1 const& a1, Point1 const& a2, Point2 const& b1, Point2 const& b2) | |
163 | { | |
164 | BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) ); | |
165 | BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) ); | |
166 | ||
167 | // TODO: check only 2 first coordinates here? | |
168 | using geometry::detail::equals::equals_point_point; | |
169 | bool a_is_point = equals_point_point(a1, a2); | |
170 | bool b_is_point = equals_point_point(b1, b2); | |
171 | ||
172 | if(a_is_point && b_is_point) | |
173 | { | |
174 | return equals_point_point(a1, b2) | |
175 | ? Policy::degenerate(a, true) | |
176 | : Policy::disjoint() | |
177 | ; | |
178 | } | |
179 | ||
180 | typedef typename select_calculation_type | |
181 | <Segment1, Segment2, CalculationType>::type calc_t; | |
182 | ||
183 | calc_t const c0 = 0; | |
184 | calc_t const c1 = 1; | |
185 | ||
186 | typedef model::point<calc_t, 3, cs::cartesian> vec3d_t; | |
187 | ||
188 | using namespace formula; | |
189 | vec3d_t const a1v = sph_to_cart3d<vec3d_t>(a1); | |
190 | vec3d_t const a2v = sph_to_cart3d<vec3d_t>(a2); | |
191 | vec3d_t const b1v = sph_to_cart3d<vec3d_t>(b1); | |
192 | vec3d_t const b2v = sph_to_cart3d<vec3d_t>(b2); | |
193 | ||
194 | vec3d_t norm1 = cross_product(a1v, a2v); | |
195 | vec3d_t norm2 = cross_product(b1v, b2v); | |
196 | ||
197 | side_info sides; | |
198 | // not normalized normals, the same as in SSF | |
199 | sides.set<0>(sph_side_value(norm2, a1v), sph_side_value(norm2, a2v)); | |
200 | if (sides.same<0>()) | |
201 | { | |
202 | // Both points are at same side of other segment, we can leave | |
203 | return Policy::disjoint(); | |
204 | } | |
205 | // not normalized normals, the same as in SSF | |
206 | sides.set<1>(sph_side_value(norm1, b1v), sph_side_value(norm1, b2v)); | |
207 | if (sides.same<1>()) | |
208 | { | |
209 | // Both points are at same side of other segment, we can leave | |
210 | return Policy::disjoint(); | |
211 | } | |
212 | ||
213 | // NOTE: at this point the segments may still be disjoint | |
214 | ||
215 | bool collinear = sides.collinear(); | |
216 | ||
217 | calc_t const len1 = math::sqrt(dot_product(norm1, norm1)); | |
218 | calc_t const len2 = math::sqrt(dot_product(norm2, norm2)); | |
219 | ||
220 | // point or opposite sides of a sphere, assume point | |
221 | if (math::equals(len1, c0)) | |
222 | { | |
223 | a_is_point = true; | |
224 | if (sides.get<0, 0>() == 0 || sides.get<0, 1>() == 0) | |
225 | { | |
226 | sides.set<0>(0, 0); | |
227 | } | |
228 | } | |
229 | else | |
230 | { | |
231 | // normalize | |
232 | divide_value(norm1, len1); | |
233 | } | |
234 | ||
235 | if (math::equals(len2, c0)) | |
236 | { | |
237 | b_is_point = true; | |
238 | if (sides.get<1, 0>() == 0 || sides.get<1, 1>() == 0) | |
239 | { | |
240 | sides.set<1>(0, 0); | |
241 | } | |
242 | } | |
243 | else | |
244 | { | |
245 | // normalize | |
246 | divide_value(norm2, len2); | |
247 | } | |
248 | ||
249 | // check both degenerated once more | |
250 | if (a_is_point && b_is_point) | |
251 | { | |
252 | return equals_point_point(a1, b2) | |
253 | ? Policy::degenerate(a, true) | |
254 | : Policy::disjoint() | |
255 | ; | |
256 | } | |
257 | ||
258 | // NOTE: at this point one of the segments may be degenerated | |
259 | // and the segments may still be disjoint | |
260 | ||
261 | calc_t dot_n1n2 = dot_product(norm1, norm2); | |
262 | ||
263 | // NOTE: this is technically not needed since theoretically above sides | |
264 | // are calculated, but just in case check the normals. | |
265 | // Have in mind that SSF side strategy doesn't check this. | |
266 | // collinear if normals are equal or opposite: cos(a) in {-1, 1} | |
267 | if (!collinear && math::equals(math::abs(dot_n1n2), c1)) | |
268 | { | |
269 | collinear = true; | |
270 | sides.set<0>(0, 0); | |
271 | sides.set<1>(0, 0); | |
272 | } | |
273 | ||
274 | if (collinear) | |
275 | { | |
276 | if (a_is_point) | |
277 | { | |
278 | return collinear_one_degenerted<calc_t>(a, true, b1, b2, a1, a2, b1v, b2v, norm2, a1v); | |
279 | } | |
280 | else if (b_is_point) | |
281 | { | |
282 | // b2 used to be consistent with (degenerated) checks above (is it needed?) | |
283 | return collinear_one_degenerted<calc_t>(b, false, a1, a2, b1, b2, a1v, a2v, norm1, b1v); | |
284 | } | |
285 | else | |
286 | { | |
287 | calc_t dist_a1_a2, dist_a1_b1, dist_a1_b2; | |
288 | calc_t dist_b1_b2, dist_b1_a1, dist_b1_a2; | |
289 | // use shorter segment | |
290 | if (len1 <= len2) | |
291 | { | |
292 | calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, norm1, b1v, dist_a1_a2, dist_a1_b1); | |
293 | calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, norm1, b2v, dist_a1_a2, dist_a1_b2); | |
294 | dist_b1_b2 = dist_a1_b2 - dist_a1_b1; | |
295 | dist_b1_a1 = -dist_a1_b1; | |
296 | dist_b1_a2 = dist_a1_a2 - dist_a1_b1; | |
297 | } | |
298 | else | |
299 | { | |
300 | calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, norm2, a1v, dist_b1_b2, dist_b1_a1); | |
301 | calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, norm2, a2v, dist_b1_b2, dist_b1_a2); | |
302 | dist_a1_a2 = dist_b1_a2 - dist_b1_a1; | |
303 | dist_a1_b1 = -dist_b1_a1; | |
304 | dist_a1_b2 = dist_b1_b2 - dist_b1_a1; | |
305 | } | |
306 | ||
307 | segment_ratio<calc_t> ra_from(dist_b1_a1, dist_b1_b2); | |
308 | segment_ratio<calc_t> ra_to(dist_b1_a2, dist_b1_b2); | |
309 | segment_ratio<calc_t> rb_from(dist_a1_b1, dist_a1_a2); | |
310 | segment_ratio<calc_t> rb_to(dist_a1_b2, dist_a1_a2); | |
311 | ||
312 | // NOTE: this is probably not needed | |
313 | int const a1_wrt_b = position_value(c0, dist_a1_b1, dist_a1_b2); | |
314 | int const a2_wrt_b = position_value(dist_a1_a2, dist_a1_b1, dist_a1_b2); | |
315 | int const b1_wrt_a = position_value(c0, dist_b1_a1, dist_b1_a2); | |
316 | int const b2_wrt_a = position_value(dist_b1_b2, dist_b1_a1, dist_b1_a2); | |
317 | ||
318 | if (a1_wrt_b == 1) | |
319 | { | |
320 | ra_from.assign(0, dist_b1_b2); | |
321 | rb_from.assign(0, dist_a1_a2); | |
322 | } | |
323 | else if (a1_wrt_b == 3) | |
324 | { | |
325 | ra_from.assign(dist_b1_b2, dist_b1_b2); | |
326 | rb_to.assign(0, dist_a1_a2); | |
327 | } | |
328 | ||
329 | if (a2_wrt_b == 1) | |
330 | { | |
331 | ra_to.assign(0, dist_b1_b2); | |
332 | rb_from.assign(dist_a1_a2, dist_a1_a2); | |
333 | } | |
334 | else if (a2_wrt_b == 3) | |
335 | { | |
336 | ra_to.assign(dist_b1_b2, dist_b1_b2); | |
337 | rb_to.assign(dist_a1_a2, dist_a1_a2); | |
338 | } | |
339 | ||
340 | if ((a1_wrt_b < 1 && a2_wrt_b < 1) || (a1_wrt_b > 3 && a2_wrt_b > 3)) | |
341 | { | |
342 | return Policy::disjoint(); | |
343 | } | |
344 | ||
345 | bool const opposite = dot_n1n2 < c0; | |
346 | ||
347 | return Policy::segments_collinear(a, b, opposite, | |
348 | a1_wrt_b, a2_wrt_b, b1_wrt_a, b2_wrt_a, | |
349 | ra_from, ra_to, rb_from, rb_to); | |
350 | } | |
351 | } | |
352 | else // crossing | |
353 | { | |
354 | if (a_is_point || b_is_point) | |
355 | { | |
356 | return Policy::disjoint(); | |
357 | } | |
358 | ||
359 | vec3d_t i1; | |
360 | intersection_point_flag ip_flag; | |
361 | calc_t dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1; | |
362 | if (calculate_ip_data(a1, a2, b1, b2, a1v, a2v, b1v, b2v, norm1, norm2, sides, | |
363 | i1, dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1, ip_flag)) | |
364 | { | |
365 | // intersects | |
366 | segment_intersection_info | |
367 | < | |
368 | calc_t, | |
369 | segment_ratio<calc_t>, | |
370 | vec3d_t | |
371 | > sinfo; | |
372 | ||
373 | sinfo.robust_ra.assign(dist_a1_i1, dist_a1_a2); | |
374 | sinfo.robust_rb.assign(dist_b1_i1, dist_b1_b2); | |
375 | sinfo.intersection_point = i1; | |
376 | sinfo.ip_flag = ip_flag; | |
377 | ||
378 | return Policy::segments_crosses(sides, sinfo, a, b); | |
379 | } | |
380 | else | |
381 | { | |
382 | return Policy::disjoint(); | |
383 | } | |
384 | } | |
385 | } | |
386 | ||
387 | private: | |
388 | template <typename CalcT, typename Segment, typename Point1, typename Point2, typename Vec3d> | |
389 | static inline return_type collinear_one_degenerted(Segment const& segment, bool degenerated_a, | |
390 | Point1 const& a1, Point1 const& a2, | |
391 | Point2 const& b1, Point2 const& b2, | |
392 | Vec3d const& v1, Vec3d const& v2, Vec3d const& norm, | |
393 | Vec3d const& vother) | |
394 | { | |
395 | CalcT dist_1_2, dist_1_o; | |
396 | return ! calculate_collinear_data(a1, a2, b1, b2, v1, v2, norm, vother, dist_1_2, dist_1_o) | |
397 | ? Policy::disjoint() | |
398 | : Policy::one_degenerate(segment, segment_ratio<CalcT>(dist_1_o, dist_1_2), degenerated_a); | |
399 | } | |
400 | ||
401 | template <typename Point1, typename Point2, typename Vec3d, typename CalcT> | |
402 | static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, | |
403 | Point2 const& b1, Point2 const& b2, | |
404 | Vec3d const& a1v, // in | |
405 | Vec3d const& a2v, // in | |
406 | Vec3d const& norm1, // in | |
407 | Vec3d const& b1v_or_b2v, // in | |
408 | CalcT& dist_a1_a2, CalcT& dist_a1_i1) // out | |
409 | { | |
410 | // calculate dist_a1_a2 and dist_a1_i1 | |
411 | calculate_dists(a1v, a2v, norm1, b1v_or_b2v, dist_a1_a2, dist_a1_i1); | |
412 | ||
413 | // if i1 is close to a1 and b1 or b2 is equal to a1 | |
414 | if (is_endpoint_equal(dist_a1_i1, a1, b1, b2)) | |
415 | { | |
416 | dist_a1_i1 = 0; | |
417 | return true; | |
418 | } | |
419 | // or i1 is close to a2 and b1 or b2 is equal to a2 | |
420 | else if (is_endpoint_equal(dist_a1_a2 - dist_a1_i1, a2, b1, b2)) | |
421 | { | |
422 | dist_a1_i1 = dist_a1_a2; | |
423 | return true; | |
424 | } | |
425 | ||
426 | // or i1 is on b | |
427 | return segment_ratio<CalcT>(dist_a1_i1, dist_a1_a2).on_segment(); | |
428 | } | |
429 | ||
430 | template <typename Point1, typename Point2, typename Vec3d, typename CalcT> | |
431 | static inline bool calculate_ip_data(Point1 const& a1, Point1 const& a2, // in | |
432 | Point2 const& b1, Point2 const& b2, // in | |
433 | Vec3d const& a1v, Vec3d const& a2v, // in | |
434 | Vec3d const& b1v, Vec3d const& b2v, // in | |
435 | Vec3d const& norm1, Vec3d const& norm2, // in | |
436 | side_info const& sides, // in | |
437 | Vec3d & i1, // in-out | |
438 | CalcT& dist_a1_a2, CalcT& dist_a1_i1, // out | |
439 | CalcT& dist_b1_b2, CalcT& dist_b1_i1, // out | |
440 | intersection_point_flag& ip_flag) // out | |
441 | { | |
442 | // great circles intersections | |
443 | i1 = cross_product(norm1, norm2); | |
444 | // NOTE: the length should be greater than 0 at this point | |
445 | // if the normals were not normalized and their dot product | |
446 | // not checked before this function is called the length | |
447 | // should be checked here (math::equals(len, c0)) | |
448 | CalcT const len = math::sqrt(dot_product(i1, i1)); | |
449 | divide_value(i1, len); // normalize i1 | |
450 | ||
451 | calculate_dists(a1v, a2v, norm1, i1, dist_a1_a2, dist_a1_i1); | |
452 | ||
453 | // choose the opposite side of the globe if the distance is shorter | |
454 | { | |
455 | CalcT const d = abs_distance(dist_a1_a2, dist_a1_i1); | |
456 | if (d > CalcT(0)) | |
457 | { | |
458 | CalcT const dist_a1_i2 = dist_of_i2(dist_a1_i1); | |
459 | CalcT const d2 = abs_distance(dist_a1_a2, dist_a1_i2); | |
460 | if (d2 < d) | |
461 | { | |
462 | dist_a1_i1 = dist_a1_i2; | |
463 | multiply_value(i1, CalcT(-1)); // the opposite intersection | |
464 | } | |
465 | } | |
466 | } | |
467 | ||
468 | bool is_on_a = false, is_near_a1 = false, is_near_a2 = false; | |
469 | if (! is_potentially_crossing(dist_a1_a2, dist_a1_i1, is_on_a, is_near_a1, is_near_a2)) | |
470 | { | |
471 | return false; | |
472 | } | |
473 | ||
474 | calculate_dists(b1v, b2v, norm2, i1, dist_b1_b2, dist_b1_i1); | |
475 | ||
476 | bool is_on_b = false, is_near_b1 = false, is_near_b2 = false; | |
477 | if (! is_potentially_crossing(dist_b1_b2, dist_b1_i1, is_on_b, is_near_b1, is_near_b2)) | |
478 | { | |
479 | return false; | |
480 | } | |
481 | ||
482 | // reassign the IP if some endpoints overlap | |
483 | using geometry::detail::equals::equals_point_point; | |
484 | if (is_near_a1) | |
485 | { | |
486 | if (is_near_b1 && equals_point_point(a1, b1)) | |
487 | { | |
488 | dist_a1_i1 = 0; | |
489 | dist_b1_i1 = 0; | |
490 | //i1 = a1v; | |
491 | ip_flag = ipi_at_a1; | |
492 | return true; | |
493 | } | |
494 | ||
495 | if (is_near_b2 && equals_point_point(a1, b2)) | |
496 | { | |
497 | dist_a1_i1 = 0; | |
498 | dist_b1_i1 = dist_b1_b2; | |
499 | //i1 = a1v; | |
500 | ip_flag = ipi_at_a1; | |
501 | return true; | |
502 | } | |
503 | } | |
504 | ||
505 | if (is_near_a2) | |
506 | { | |
507 | if (is_near_b1 && equals_point_point(a2, b1)) | |
508 | { | |
509 | dist_a1_i1 = dist_a1_a2; | |
510 | dist_b1_i1 = 0; | |
511 | //i1 = a2v; | |
512 | ip_flag = ipi_at_a2; | |
513 | return true; | |
514 | } | |
515 | ||
516 | if (is_near_b2 && equals_point_point(a2, b2)) | |
517 | { | |
518 | dist_a1_i1 = dist_a1_a2; | |
519 | dist_b1_i1 = dist_b1_b2; | |
520 | //i1 = a2v; | |
521 | ip_flag = ipi_at_a2; | |
522 | return true; | |
523 | } | |
524 | } | |
525 | ||
526 | // at this point we know that the endpoints doesn't overlap | |
527 | // reassign IP and distance if the IP is on a segment and one of | |
528 | // the endpoints of the other segment lies on the former segment | |
529 | if (is_on_a) | |
530 | { | |
531 | if (is_near_b1 && sides.template get<1, 0>() == 0) // b1 wrt a | |
532 | { | |
533 | dist_b1_i1 = 0; | |
534 | //i1 = b1v; | |
535 | ip_flag = ipi_at_b1; | |
536 | return true; | |
537 | } | |
538 | ||
539 | if (is_near_b2 && sides.template get<1, 1>() == 0) // b2 wrt a | |
540 | { | |
541 | dist_b1_i1 = dist_b1_b2; | |
542 | //i1 = b2v; | |
543 | ip_flag = ipi_at_b2; | |
544 | return true; | |
545 | } | |
546 | } | |
547 | ||
548 | if (is_on_b) | |
549 | { | |
550 | if (is_near_a1 && sides.template get<0, 0>() == 0) // a1 wrt b | |
551 | { | |
552 | dist_a1_i1 = 0; | |
553 | //i1 = a1v; | |
554 | ip_flag = ipi_at_a1; | |
555 | return true; | |
556 | } | |
557 | ||
558 | if (is_near_a2 && sides.template get<0, 1>() == 0) // a2 wrt b | |
559 | { | |
560 | dist_a1_i1 = dist_a1_a2; | |
561 | //i1 = a2v; | |
562 | ip_flag = ipi_at_a2; | |
563 | return true; | |
564 | } | |
565 | } | |
566 | ||
567 | ip_flag = ipi_inters; | |
568 | ||
569 | return is_on_a && is_on_b; | |
570 | } | |
571 | ||
572 | template <typename Vec3d, typename CalcT> | |
573 | static inline void calculate_dists(Vec3d const& a1v, // in | |
574 | Vec3d const& a2v, // in | |
575 | Vec3d const& norm1, // in | |
576 | Vec3d const& i1, // in | |
577 | CalcT& dist_a1_a2, CalcT& dist_a1_i1) // out | |
578 | { | |
579 | CalcT const c0 = 0; | |
580 | CalcT const c1 = 1; | |
581 | CalcT const c2 = 2; | |
582 | CalcT const c4 = 4; | |
583 | ||
584 | CalcT cos_a1_a2 = dot_product(a1v, a2v); | |
585 | dist_a1_a2 = -cos_a1_a2 + c1; // [1, -1] -> [0, 2] representing [0, pi] | |
586 | ||
587 | CalcT cos_a1_i1 = dot_product(a1v, i1); | |
588 | dist_a1_i1 = -cos_a1_i1 + c1; // [0, 2] representing [0, pi] | |
589 | if (dot_product(norm1, cross_product(a1v, i1)) < c0) // left or right of a1 on a | |
590 | { | |
591 | dist_a1_i1 = -dist_a1_i1; // [0, 2] -> [0, -2] representing [0, -pi] | |
592 | } | |
593 | if (dist_a1_i1 <= -c2) // <= -pi | |
594 | { | |
595 | dist_a1_i1 += c4; // += 2pi | |
596 | } | |
597 | } | |
598 | ||
599 | // the dist of the ip on the other side of the sphere | |
600 | template <typename CalcT> | |
601 | static inline CalcT dist_of_i2(CalcT const& dist_a1_i1) | |
602 | { | |
603 | CalcT const c2 = 2; | |
604 | CalcT const c4 = 4; | |
605 | ||
606 | CalcT dist_a1_i2 = dist_a1_i1 - c2; // dist_a1_i2 = dist_a1_i1 - pi; | |
607 | if (dist_a1_i2 <= -c2) // <= -pi | |
608 | { | |
609 | dist_a1_i2 += c4; // += 2pi; | |
610 | } | |
611 | return dist_a1_i2; | |
612 | } | |
613 | ||
614 | template <typename CalcT> | |
615 | static inline CalcT abs_distance(CalcT const& dist_a1_a2, CalcT const& dist_a1_i1) | |
616 | { | |
617 | if (dist_a1_i1 < CalcT(0)) | |
618 | return -dist_a1_i1; | |
619 | else if (dist_a1_i1 > dist_a1_a2) | |
620 | return dist_a1_i1 - dist_a1_a2; | |
621 | else | |
622 | return CalcT(0); | |
623 | } | |
624 | ||
625 | template <typename CalcT> | |
626 | static inline bool is_potentially_crossing(CalcT const& dist_a1_a2, CalcT const& dist_a1_i1, // in | |
627 | bool& is_on_a, bool& is_near_a1, bool& is_near_a2) // out | |
628 | { | |
629 | is_on_a = segment_ratio<CalcT>(dist_a1_i1, dist_a1_a2).on_segment(); | |
630 | is_near_a1 = is_near(dist_a1_i1); | |
631 | is_near_a2 = is_near(dist_a1_a2 - dist_a1_i1); | |
632 | return is_on_a || is_near_a1 || is_near_a2; | |
633 | } | |
634 | ||
635 | template <typename CalcT, typename P1, typename P2> | |
636 | static inline bool is_endpoint_equal(CalcT const& dist, | |
637 | P1 const& ai, P2 const& b1, P2 const& b2) | |
638 | { | |
639 | using geometry::detail::equals::equals_point_point; | |
640 | return is_near(dist) && (equals_point_point(ai, b1) || equals_point_point(ai, b2)); | |
641 | } | |
642 | ||
643 | template <typename CalcT> | |
644 | static inline bool is_near(CalcT const& dist) | |
645 | { | |
646 | CalcT const small_number = CalcT(boost::is_same<CalcT, float>::value ? 0.0001 : 0.00000001); | |
647 | return math::abs(dist) <= small_number; | |
648 | } | |
649 | ||
650 | template <typename ProjCoord1, typename ProjCoord2> | |
651 | static inline int position_value(ProjCoord1 const& ca1, | |
652 | ProjCoord2 const& cb1, | |
653 | ProjCoord2 const& cb2) | |
654 | { | |
655 | // S1x 0 1 2 3 4 | |
656 | // S2 |----------> | |
657 | return math::equals(ca1, cb1) ? 1 | |
658 | : math::equals(ca1, cb2) ? 3 | |
659 | : cb1 < cb2 ? | |
660 | ( ca1 < cb1 ? 0 | |
661 | : ca1 > cb2 ? 4 | |
662 | : 2 ) | |
663 | : ( ca1 > cb1 ? 0 | |
664 | : ca1 < cb2 ? 4 | |
665 | : 2 ); | |
666 | } | |
667 | }; | |
668 | ||
669 | ||
670 | #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS | |
671 | namespace services | |
672 | { | |
673 | ||
674 | /*template <typename Policy, typename CalculationType> | |
675 | struct default_strategy<spherical_polar_tag, Policy, CalculationType> | |
676 | { | |
677 | typedef relate_spherical_segments<Policy, CalculationType> type; | |
678 | };*/ | |
679 | ||
680 | template <typename Policy, typename CalculationType> | |
681 | struct default_strategy<spherical_equatorial_tag, Policy, CalculationType> | |
682 | { | |
683 | typedef relate_spherical_segments<Policy, CalculationType> type; | |
684 | }; | |
685 | ||
686 | template <typename Policy, typename CalculationType> | |
687 | struct default_strategy<geographic_tag, Policy, CalculationType> | |
688 | { | |
689 | typedef relate_spherical_segments<Policy, CalculationType> type; | |
690 | }; | |
691 | ||
692 | } // namespace services | |
693 | #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS | |
694 | ||
695 | ||
696 | }} // namespace strategy::intersection | |
697 | ||
698 | }} // namespace boost::geometry | |
699 | ||
700 | ||
701 | #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_INTERSECTION_HPP |