1 // Boost.Geometry
2
3 // Copyright (c) 2016 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_FORMULAS_SJOBERG_INTERSECTION_HPP
12 #define BOOST_GEOMETRY_FORMULAS_SJOBERG_INTERSECTION_HPP
13
14
15 #include <boost/math/constants/constants.hpp>
16
17 #include <boost/geometry/core/radius.hpp>
18 #include <boost/geometry/core/srs.hpp>
19
20 #include <boost/geometry/util/condition.hpp>
21 #include <boost/geometry/util/math.hpp>
22
23 #include <boost/geometry/formulas/flattening.hpp>
24 #include <boost/geometry/formulas/spherical.hpp>
25
26
27 namespace boost { namespace geometry { namespace formula
28 {
29
30 /*!
31 \brief The intersection of two great circles as proposed by Sjoberg.
32 \see See
33 - [Sjoberg02] Lars E. Sjoberg, Intersections on the sphere and ellipsoid, 2002
34 http://link.springer.com/article/10.1007/s00190-001-0230-9
35 */
36 template <typename CT>
37 struct sjoberg_intersection_spherical_02
38 {
39 // TODO: if it will be used as standalone formula
40 // support segments on equator and endpoints on poles
41
applyboost::geometry::formula::sjoberg_intersection_spherical_0242 static inline bool apply(CT const& lon1, CT const& lat1, CT const& lon_a2, CT const& lat_a2,
43 CT const& lon2, CT const& lat2, CT const& lon_b2, CT const& lat_b2,
44 CT & lon, CT & lat)
45 {
46 CT tan_lat = 0;
47 bool res = apply_alt(lon1, lat1, lon_a2, lat_a2,
48 lon2, lat2, lon_b2, lat_b2,
49 lon, tan_lat);
50
51 if (res)
52 {
53 lat = atan(tan_lat);
54 }
55
56 return res;
57 }
58
apply_altboost::geometry::formula::sjoberg_intersection_spherical_0259 static inline bool apply_alt(CT const& lon1, CT const& lat1, CT const& lon_a2, CT const& lat_a2,
60 CT const& lon2, CT const& lat2, CT const& lon_b2, CT const& lat_b2,
61 CT & lon, CT & tan_lat)
62 {
63 CT const cos_lon1 = cos(lon1);
64 CT const sin_lon1 = sin(lon1);
65 CT const cos_lon2 = cos(lon2);
66 CT const sin_lon2 = sin(lon2);
67 CT const sin_lat1 = sin(lat1);
68 CT const sin_lat2 = sin(lat2);
69 CT const cos_lat1 = cos(lat1);
70 CT const cos_lat2 = cos(lat2);
71
72 CT const tan_lat_a2 = tan(lat_a2);
73 CT const tan_lat_b2 = tan(lat_b2);
74
75 return apply(lon1, lon_a2, lon2, lon_b2,
76 sin_lon1, cos_lon1, sin_lat1, cos_lat1,
77 sin_lon2, cos_lon2, sin_lat2, cos_lat2,
78 tan_lat_a2, tan_lat_b2,
79 lon, tan_lat);
80 }
81
82 private:
applyboost::geometry::formula::sjoberg_intersection_spherical_0283 static inline bool apply(CT const& lon1, CT const& lon_a2, CT const& lon2, CT const& lon_b2,
84 CT const& sin_lon1, CT const& cos_lon1, CT const& sin_lat1, CT const& cos_lat1,
85 CT const& sin_lon2, CT const& cos_lon2, CT const& sin_lat2, CT const& cos_lat2,
86 CT const& tan_lat_a2, CT const& tan_lat_b2,
87 CT & lon, CT & tan_lat)
88 {
89 // NOTE:
90 // cos_lat_ = 0 <=> segment on equator
91 // tan_alpha_ = 0 <=> segment vertical
92
93 CT const tan_lat1 = sin_lat1 / cos_lat1; //tan(lat1);
94 CT const tan_lat2 = sin_lat2 / cos_lat2; //tan(lat2);
95
96 CT const dlon1 = lon_a2 - lon1;
97 CT const sin_dlon1 = sin(dlon1);
98 CT const dlon2 = lon_b2 - lon2;
99 CT const sin_dlon2 = sin(dlon2);
100
101 CT const c0 = 0;
102 bool const is_vertical1 = math::equals(sin_dlon1, c0);
103 bool const is_vertical2 = math::equals(sin_dlon2, c0);
104
105 CT tan_alpha1 = 0;
106 CT tan_alpha2 = 0;
107
108 if (is_vertical1 && is_vertical2)
109 {
110 // circles intersect at one of the poles or are collinear
111 return false;
112 }
113 else if (is_vertical1)
114 {
115 CT const cos_dlon2 = cos(dlon2);
116 CT const tan_alpha2_x = cos_lat2 * tan_lat_b2 - sin_lat2 * cos_dlon2;
117 tan_alpha2 = sin_dlon2 / tan_alpha2_x;
118
119 lon = lon1;
120 }
121 else if (is_vertical2)
122 {
123 CT const cos_dlon1 = cos(dlon1);
124 CT const tan_alpha1_x = cos_lat1 * tan_lat_a2 - sin_lat1 * cos_dlon1;
125 tan_alpha1 = sin_dlon1 / tan_alpha1_x;
126
127 lon = lon2;
128 }
129 else
130 {
131 CT const cos_dlon1 = cos(dlon1);
132 CT const cos_dlon2 = cos(dlon2);
133
134 CT const tan_alpha1_x = cos_lat1 * tan_lat_a2 - sin_lat1 * cos_dlon1;
135 CT const tan_alpha2_x = cos_lat2 * tan_lat_b2 - sin_lat2 * cos_dlon2;
136 tan_alpha1 = sin_dlon1 / tan_alpha1_x;
137 tan_alpha2 = sin_dlon2 / tan_alpha2_x;
138
139 CT const T1 = tan_alpha1 * cos_lat1;
140 CT const T2 = tan_alpha2 * cos_lat2;
141 CT const T1T2 = T1*T2;
142 CT const tan_lon_y = T1 * sin_lon2 - T2 * sin_lon1 + T1T2 * (tan_lat1 * cos_lon1 - tan_lat2 * cos_lon2);
143 CT const tan_lon_x = T1 * cos_lon2 - T2 * cos_lon1 - T1T2 * (tan_lat1 * sin_lon1 - tan_lat2 * sin_lon2);
144
145 lon = atan2(tan_lon_y, tan_lon_x);
146 }
147
148 // choose closer result
149 CT const pi = math::pi<CT>();
150 CT const lon_2 = lon > c0 ? lon - pi : lon + pi;
151 CT const lon_dist1 = (std::max)((std::min)(math::longitude_difference<radian>(lon1, lon),
152 math::longitude_difference<radian>(lon_a2, lon)),
153 (std::min)(math::longitude_difference<radian>(lon2, lon),
154 math::longitude_difference<radian>(lon_b2, lon)));
155 CT const lon_dist2 = (std::max)((std::min)(math::longitude_difference<radian>(lon1, lon_2),
156 math::longitude_difference<radian>(lon_a2, lon_2)),
157 (std::min)(math::longitude_difference<radian>(lon2, lon_2),
158 math::longitude_difference<radian>(lon_b2, lon_2)));
159 if (lon_dist2 < lon_dist1)
160 {
161 lon = lon_2;
162 }
163
164 CT const sin_lon = sin(lon);
165 CT const cos_lon = cos(lon);
166
167 if (math::abs(tan_alpha1) >= math::abs(tan_alpha2)) // pick less vertical segment
168 {
169 CT const sin_dlon_1 = sin_lon * cos_lon1 - cos_lon * sin_lon1;
170 CT const cos_dlon_1 = cos_lon * cos_lon1 + sin_lon * sin_lon1;
171 CT const lat_y_1 = sin_dlon_1 + tan_alpha1 * sin_lat1 * cos_dlon_1;
172 CT const lat_x_1 = tan_alpha1 * cos_lat1;
173 tan_lat = lat_y_1 / lat_x_1;
174 }
175 else
176 {
177 CT const sin_dlon_2 = sin_lon * cos_lon2 - cos_lon * sin_lon2;
178 CT const cos_dlon_2 = cos_lon * cos_lon2 + sin_lon * sin_lon2;
179 CT const lat_y_2 = sin_dlon_2 + tan_alpha2 * sin_lat2 * cos_dlon_2;
180 CT const lat_x_2 = tan_alpha2 * cos_lat2;
181 tan_lat = lat_y_2 / lat_x_2;
182 }
183
184 return true;
185 }
186 };
187
188
189 /*! Approximation of dLambda_j [Sjoberg07], expanded into taylor series in e^2
190 Maxima script:
191 dLI_j(c_j, sinB_j, sinB) := integrate(1 / (sqrt(1 - c_j ^ 2 - x ^ 2)*(1 + sqrt(1 - e2*(1 - x ^ 2)))), x, sinB_j, sinB);
192 dL_j(c_j, B_j, B) := -e2 * c_j * dLI_j(c_j, B_j, B);
193 S: taylor(dLI_j(c_j, sinB_j, sinB), e2, 0, 3);
194 assume(c_j < 1);
195 assume(c_j > 0);
196 L1: factor(integrate(sqrt(-x ^ 2 - c_j ^ 2 + 1) / (x ^ 2 + c_j ^ 2 - 1), x));
197 L2: factor(integrate(((x ^ 2 - 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x));
198 L3: factor(integrate(((x ^ 4 - 2 * x ^ 2 + 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x));
199 L4: factor(integrate(((x ^ 6 - 3 * x ^ 4 + 3 * x ^ 2 - 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x));
200
201 \see See
202 - [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007
203 http://link.springer.com/article/10.1007/s00190-007-0204-7
204 */
205 template <unsigned int Order, typename CT>
sjoberg_d_lambda_e_sqr(CT const & sin_betaj,CT const & sin_beta,CT const & Cj,CT const & sqrt_1_Cj_sqr,CT const & e_sqr)206 inline CT sjoberg_d_lambda_e_sqr(CT const& sin_betaj, CT const& sin_beta,
207 CT const& Cj, CT const& sqrt_1_Cj_sqr,
208 CT const& e_sqr)
209 {
210 using math::detail::bounded;
211
212 if (Order == 0)
213 {
214 return 0;
215 }
216
217 CT const c1 = 1;
218 CT const c2 = 2;
219
220 CT const asin_B = asin(bounded(sin_beta / sqrt_1_Cj_sqr, -c1, c1));
221 CT const asin_Bj = asin(sin_betaj / sqrt_1_Cj_sqr);
222 CT const L0 = (asin_B - asin_Bj) / c2;
223
224 if (Order == 1)
225 {
226 return -Cj * e_sqr * L0;
227 }
228
229 CT const c0 = 0;
230 CT const c16 = 16;
231
232 CT const X = sin_beta;
233 CT const Xj = sin_betaj;
234 CT const X_sqr = math::sqr(X);
235 CT const Xj_sqr = math::sqr(Xj);
236 CT const Cj_sqr = math::sqr(Cj);
237 CT const Cj_sqr_plus_one = Cj_sqr + c1;
238 CT const one_minus_Cj_sqr = c1 - Cj_sqr;
239 CT const sqrt_Y = math::sqrt(bounded(-X_sqr + one_minus_Cj_sqr, c0));
240 CT const sqrt_Yj = math::sqrt(-Xj_sqr + one_minus_Cj_sqr);
241 CT const L1 = (Cj_sqr_plus_one * (asin_B - asin_Bj) + X * sqrt_Y - Xj * sqrt_Yj) / c16;
242
243 if (Order == 2)
244 {
245 return -Cj * e_sqr * (L0 + e_sqr * L1);
246 }
247
248 CT const c3 = 3;
249 CT const c5 = 5;
250 CT const c128 = 128;
251
252 CT const E = Cj_sqr * (c3 * Cj_sqr + c2) + c3;
253 CT const F = X * (-c2 * X_sqr + c3 * Cj_sqr + c5);
254 CT const Fj = Xj * (-c2 * Xj_sqr + c3 * Cj_sqr + c5);
255 CT const L2 = (E * (asin_B - asin_Bj) + F * sqrt_Y - Fj * sqrt_Yj) / c128;
256
257 if (Order == 3)
258 {
259 return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * L2));
260 }
261
262 CT const c8 = 8;
263 CT const c9 = 9;
264 CT const c10 = 10;
265 CT const c15 = 15;
266 CT const c24 = 24;
267 CT const c26 = 26;
268 CT const c33 = 33;
269 CT const c6144 = 6144;
270
271 CT const G = Cj_sqr * (Cj_sqr * (Cj_sqr * c15 + c9) + c9) + c15;
272 CT const H = -c10 * Cj_sqr - c26;
273 CT const I = Cj_sqr * (Cj_sqr * c15 + c24) + c33;
274 CT const J = X_sqr * (X * (c8 * X_sqr + H)) + X * I;
275 CT const Jj = Xj_sqr * (Xj * (c8 * Xj_sqr + H)) + Xj * I;
276 CT const L3 = (G * (asin_B - asin_Bj) + J * sqrt_Y - Jj * sqrt_Yj) / c6144;
277
278 // Order 4 and higher
279 return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * (L2 + e_sqr * L3)));
280 }
281
282 /*!
283 \brief The representation of geodesic as proposed by Sjoberg.
284 \see See
285 - [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007
286 http://link.springer.com/article/10.1007/s00190-007-0204-7
287 - [Sjoberg12] Lars E. Sjoberg, Solutions to the ellipsoidal Clairaut constant
288 and the inverse geodetic problem by numerical integration, 2012
289 https://www.degruyter.com/view/j/jogs.2012.2.issue-3/v10156-011-0037-4/v10156-011-0037-4.xml
290 */
291 template <typename CT, unsigned int Order>
292 class sjoberg_geodesic
293 {
sjoberg_geodesic()294 sjoberg_geodesic() {}
295
sign_C(CT const & alphaj)296 static int sign_C(CT const& alphaj)
297 {
298 CT const c0 = 0;
299 CT const c2 = 2;
300 CT const pi = math::pi<CT>();
301 CT const pi_half = pi / c2;
302
303 return (pi_half < alphaj && alphaj < pi) || (-pi_half < alphaj && alphaj < c0) ? -1 : 1;
304 }
305
306 public:
sjoberg_geodesic(CT const & lon,CT const & lat,CT const & alpha,CT const & f)307 sjoberg_geodesic(CT const& lon, CT const& lat, CT const& alpha, CT const& f)
308 : lonj(lon)
309 , latj(lat)
310 , alphaj(alpha)
311 {
312 CT const c0 = 0;
313 CT const c1 = 1;
314 CT const c2 = 2;
315 //CT const pi = math::pi<CT>();
316 //CT const pi_half = pi / c2;
317
318 one_minus_f = c1 - f;
319 e_sqr = f * (c2 - f);
320
321 tan_latj = tan(lat);
322 tan_betaj = one_minus_f * tan_latj;
323 betaj = atan(tan_betaj);
324 sin_betaj = sin(betaj);
325
326 cos_betaj = cos(betaj);
327 sin_alphaj = sin(alphaj);
328 // Clairaut constant (lower-case in the paper)
329 Cj = sign_C(alphaj) * cos_betaj * sin_alphaj;
330 Cj_sqr = math::sqr(Cj);
331 sqrt_1_Cj_sqr = math::sqrt(c1 - Cj_sqr);
332
333 sign_lon_diff = alphaj >= 0 ? 1 : -1; // || alphaj == -pi ?
334 //sign_lon_diff = 1;
335
336 is_on_equator = math::equals(sqrt_1_Cj_sqr, c0);
337 is_Cj_zero = math::equals(Cj, c0);
338
339 t0j = c0;
340 asin_tj_t0j = c0;
341
342 if (! is_Cj_zero)
343 {
344 t0j = sqrt_1_Cj_sqr / Cj;
345 }
346
347 if (! is_on_equator)
348 {
349 //asin_tj_t0j = asin(tan_betaj / t0j);
350 asin_tj_t0j = asin(tan_betaj * Cj / sqrt_1_Cj_sqr);
351 }
352 }
353
354 struct vertex_data
355 {
356 //CT beta0j;
357 CT sin_beta0j;
358 CT dL0j;
359 CT lon0j;
360 };
361
get_vertex_data() const362 vertex_data get_vertex_data() const
363 {
364 CT const c2 = 2;
365 CT const pi = math::pi<CT>();
366 CT const pi_half = pi / c2;
367
368 vertex_data res;
369
370 if (! is_Cj_zero)
371 {
372 //res.beta0j = atan(t0j);
373 //res.sin_beta0j = sin(res.beta0j);
374 res.sin_beta0j = math::sign(t0j) * sqrt_1_Cj_sqr;
375 res.dL0j = d_lambda(res.sin_beta0j);
376 res.lon0j = lonj + sign_lon_diff * (pi_half - asin_tj_t0j + res.dL0j);
377 }
378 else
379 {
380 //res.beta0j = pi_half;
381 //res.sin_beta0j = betaj >= 0 ? 1 : -1;
382 res.sin_beta0j = 1;
383 res.dL0j = 0;
384 res.lon0j = lonj;
385 }
386
387 return res;
388 }
389
is_sin_beta_ok(CT const & sin_beta) const390 bool is_sin_beta_ok(CT const& sin_beta) const
391 {
392 CT const c1 = 1;
393 return math::abs(sin_beta / sqrt_1_Cj_sqr) <= c1;
394 }
395
k_diff(CT const & sin_beta,CT & delta_k) const396 bool k_diff(CT const& sin_beta,
397 CT & delta_k) const
398 {
399 if (is_Cj_zero)
400 {
401 delta_k = 0;
402 return true;
403 }
404
405 // beta out of bounds and not close
406 if (! (is_sin_beta_ok(sin_beta)
407 || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) )
408 {
409 return false;
410 }
411
412 // NOTE: beta may be slightly out of bounds here but d_lambda handles that
413 CT const dLj = d_lambda(sin_beta);
414 delta_k = sign_lon_diff * (/*asin_t_t0j*/ - asin_tj_t0j + dLj);
415
416 return true;
417 }
418
lon_diff(CT const & sin_beta,CT const & t,CT & delta_lon) const419 bool lon_diff(CT const& sin_beta, CT const& t,
420 CT & delta_lon) const
421 {
422 using math::detail::bounded;
423 CT const c1 = 1;
424
425 if (is_Cj_zero)
426 {
427 delta_lon = 0;
428 return true;
429 }
430
431 CT delta_k = 0;
432 if (! k_diff(sin_beta, delta_k))
433 {
434 return false;
435 }
436
437 CT const t_t0j = t / t0j;
438 // NOTE: t may be slightly out of bounds here
439 CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1));
440 delta_lon = sign_lon_diff * asin_t_t0j + delta_k;
441
442 return true;
443 }
444
k_diffs(CT const & sin_beta,vertex_data const & vd,CT & delta_k_before,CT & delta_k_behind,bool check_sin_beta=true) const445 bool k_diffs(CT const& sin_beta, vertex_data const& vd,
446 CT & delta_k_before, CT & delta_k_behind,
447 bool check_sin_beta = true) const
448 {
449 CT const pi = math::pi<CT>();
450
451 if (is_Cj_zero)
452 {
453 delta_k_before = 0;
454 delta_k_behind = sign_lon_diff * pi;
455 return true;
456 }
457
458 // beta out of bounds and not close
459 if (check_sin_beta
460 && ! (is_sin_beta_ok(sin_beta)
461 || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) )
462 {
463 return false;
464 }
465
466 // NOTE: beta may be slightly out of bounds here but d_lambda handles that
467 CT const dLj = d_lambda(sin_beta);
468 delta_k_before = sign_lon_diff * (/*asin_t_t0j*/ - asin_tj_t0j + dLj);
469
470 // This version require no additional dLj calculation
471 delta_k_behind = sign_lon_diff * (pi /*- asin_t_t0j*/ - asin_tj_t0j + vd.dL0j + (vd.dL0j - dLj));
472
473 // [Sjoberg12]
474 //CT const dL101 = d_lambda(sin_betaj, vd.sin_beta0j);
475 // WARNING: the following call might not work if beta was OoB because only the second argument is bounded
476 //CT const dL_01 = d_lambda(sin_beta, vd.sin_beta0j);
477 //delta_k_behind = sign_lon_diff * (pi /*- asin_t_t0j*/ - asin_tj_t0j + dL101 + dL_01);
478
479 return true;
480 }
481
lon_diffs(CT const & sin_beta,CT const & t,vertex_data const & vd,CT & delta_lon_before,CT & delta_lon_behind) const482 bool lon_diffs(CT const& sin_beta, CT const& t, vertex_data const& vd,
483 CT & delta_lon_before, CT & delta_lon_behind) const
484 {
485 using math::detail::bounded;
486 CT const c1 = 1;
487 CT const pi = math::pi<CT>();
488
489 if (is_Cj_zero)
490 {
491 delta_lon_before = 0;
492 delta_lon_behind = sign_lon_diff * pi;
493 return true;
494 }
495
496 CT delta_k_before = 0, delta_k_behind = 0;
497 if (! k_diffs(sin_beta, vd, delta_k_before, delta_k_behind))
498 {
499 return false;
500 }
501
502 CT const t_t0j = t / t0j;
503 // NOTE: t may be slightly out of bounds here
504 CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1));
505 CT const sign_asin_t_t0j = sign_lon_diff * asin_t_t0j;
506 delta_lon_before = sign_asin_t_t0j + delta_k_before;
507 delta_lon_behind = -sign_asin_t_t0j + delta_k_behind;
508
509 return true;
510 }
511
lon(CT const & sin_beta,CT const & t,vertex_data const & vd,CT & lon_before,CT & lon_behind) const512 bool lon(CT const& sin_beta, CT const& t, vertex_data const& vd,
513 CT & lon_before, CT & lon_behind) const
514 {
515 using math::detail::bounded;
516 CT const c1 = 1;
517 CT const pi = math::pi<CT>();
518
519 if (is_Cj_zero)
520 {
521 lon_before = lonj;
522 lon_behind = lonj + sign_lon_diff * pi;
523 return true;
524 }
525
526 if (! (is_sin_beta_ok(sin_beta)
527 || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) )
528 {
529 return false;
530 }
531
532 CT const t_t0j = t / t0j;
533 CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1));
534 CT const dLj = d_lambda(sin_beta);
535 lon_before = lonj + sign_lon_diff * (asin_t_t0j - asin_tj_t0j + dLj);
536 lon_behind = vd.lon0j + (vd.lon0j - lon_before);
537
538 return true;
539 }
540
541
lon(CT const & delta_lon) const542 CT lon(CT const& delta_lon) const
543 {
544 return lonj + delta_lon;
545 }
546
lat(CT const & t) const547 CT lat(CT const& t) const
548 {
549 // t = tan(beta) = (1-f)tan(lat)
550 return atan(t / one_minus_f);
551 }
552
vertex(CT & lon,CT & lat) const553 void vertex(CT & lon, CT & lat) const
554 {
555 lon = get_vertex_data().lon0j;
556 if (! is_Cj_zero)
557 {
558 lat = sjoberg_geodesic::lat(t0j);
559 }
560 else
561 {
562 CT const c2 = 2;
563 lat = math::pi<CT>() / c2;
564 }
565 }
566
lon_of_equator_intersection() const567 CT lon_of_equator_intersection() const
568 {
569 CT const c0 = 0;
570 CT const dLj = d_lambda(c0);
571 CT const asin_tj_t0j = asin(Cj * tan_betaj / sqrt_1_Cj_sqr);
572 return lonj - asin_tj_t0j + dLj;
573 }
574
d_lambda(CT const & sin_beta) const575 CT d_lambda(CT const& sin_beta) const
576 {
577 return sjoberg_d_lambda_e_sqr<Order>(sin_betaj, sin_beta, Cj, sqrt_1_Cj_sqr, e_sqr);
578 }
579
580 // [Sjoberg12]
581 /*CT d_lambda(CT const& sin_beta1, CT const& sin_beta2) const
582 {
583 return sjoberg_d_lambda_e_sqr<Order>(sin_beta1, sin_beta2, Cj, sqrt_1_Cj_sqr, e_sqr);
584 }*/
585
586 CT lonj;
587 CT latj;
588 CT alphaj;
589
590 CT one_minus_f;
591 CT e_sqr;
592
593 CT tan_latj;
594 CT tan_betaj;
595 CT betaj;
596 CT sin_betaj;
597 CT cos_betaj;
598 CT sin_alphaj;
599 CT Cj;
600 CT Cj_sqr;
601 CT sqrt_1_Cj_sqr;
602
603 int sign_lon_diff;
604
605 bool is_on_equator;
606 bool is_Cj_zero;
607
608 CT t0j;
609 CT asin_tj_t0j;
610 };
611
612
613 /*!
614 \brief The intersection of two geodesics as proposed by Sjoberg.
615 \see See
616 - [Sjoberg02] Lars E. Sjoberg, Intersections on the sphere and ellipsoid, 2002
617 http://link.springer.com/article/10.1007/s00190-001-0230-9
618 - [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007
619 http://link.springer.com/article/10.1007/s00190-007-0204-7
620 - [Sjoberg12] Lars E. Sjoberg, Solutions to the ellipsoidal Clairaut constant
621 and the inverse geodetic problem by numerical integration, 2012
622 https://www.degruyter.com/view/j/jogs.2012.2.issue-3/v10156-011-0037-4/v10156-011-0037-4.xml
623 */
624 template
625 <
626 typename CT,
627 template <typename, bool, bool, bool, bool, bool> class Inverse,
628 unsigned int Order = 4
629 >
630 class sjoberg_intersection
631 {
632 typedef sjoberg_geodesic<CT, Order> geodesic_type;
633 typedef Inverse<CT, false, true, false, false, false> inverse_type;
634 typedef typename inverse_type::result_type inverse_result;
635
636 static bool const enable_02 = true;
637 static int const max_iterations_02 = 10;
638 static int const max_iterations_07 = 20;
639
640 public:
641 template <typename T1, typename T2, typename Spheroid>
apply(T1 const & lona1,T1 const & lata1,T1 const & lona2,T1 const & lata2,T2 const & lonb1,T2 const & latb1,T2 const & lonb2,T2 const & latb2,CT & lon,CT & lat,Spheroid const & spheroid)642 static inline bool apply(T1 const& lona1, T1 const& lata1,
643 T1 const& lona2, T1 const& lata2,
644 T2 const& lonb1, T2 const& latb1,
645 T2 const& lonb2, T2 const& latb2,
646 CT & lon, CT & lat,
647 Spheroid const& spheroid)
648 {
649 CT const lon_a1 = lona1;
650 CT const lat_a1 = lata1;
651 CT const lon_a2 = lona2;
652 CT const lat_a2 = lata2;
653 CT const lon_b1 = lonb1;
654 CT const lat_b1 = latb1;
655 CT const lon_b2 = lonb2;
656 CT const lat_b2 = latb2;
657
658 inverse_result const res1 = inverse_type::apply(lon_a1, lat_a1, lon_a2, lat_a2, spheroid);
659 inverse_result const res2 = inverse_type::apply(lon_b1, lat_b1, lon_b2, lat_b2, spheroid);
660
661 return apply(lon_a1, lat_a1, lon_a2, lat_a2, res1.azimuth,
662 lon_b1, lat_b1, lon_b2, lat_b2, res2.azimuth,
663 lon, lat, spheroid);
664 }
665
666 // TODO: Currently may not work correctly if one of the endpoints is the pole
667 template <typename Spheroid>
apply(CT const & lon_a1,CT const & lat_a1,CT const & lon_a2,CT const & lat_a2,CT const & alpha_a1,CT const & lon_b1,CT const & lat_b1,CT const & lon_b2,CT const & lat_b2,CT const & alpha_b1,CT & lon,CT & lat,Spheroid const & spheroid)668 static inline bool apply(CT const& lon_a1, CT const& lat_a1, CT const& lon_a2, CT const& lat_a2, CT const& alpha_a1,
669 CT const& lon_b1, CT const& lat_b1, CT const& lon_b2, CT const& lat_b2, CT const& alpha_b1,
670 CT & lon, CT & lat,
671 Spheroid const& spheroid)
672 {
673 // coordinates in radians
674
675 CT const c0 = 0;
676 CT const c1 = 1;
677
678 CT const f = formula::flattening<CT>(spheroid);
679 CT const one_minus_f = c1 - f;
680
681 geodesic_type geod1(lon_a1, lat_a1, alpha_a1, f);
682 geodesic_type geod2(lon_b1, lat_b1, alpha_b1, f);
683
684 // Cj = 1 if on equator <=> sqrt_1_Cj_sqr = 0
685 // Cj = 0 if vertical <=> sqrt_1_Cj_sqr = 1
686
687 if (geod1.is_on_equator && geod2.is_on_equator)
688 {
689 return false;
690 }
691 else if (geod1.is_on_equator)
692 {
693 lon = geod2.lon_of_equator_intersection();
694 lat = c0;
695 return true;
696 }
697 else if (geod2.is_on_equator)
698 {
699 lon = geod1.lon_of_equator_intersection();
700 lat = c0;
701 return true;
702 }
703
704 // (lon1 - lon2) normalized to (-180, 180]
705 CT const lon1_minus_lon2 = math::longitude_distance_signed<radian>(geod2.lonj, geod1.lonj);
706
707 // vertical segments
708 if (geod1.is_Cj_zero && geod2.is_Cj_zero)
709 {
710 CT const pi = math::pi<CT>();
711
712 // the geodesics are parallel, the intersection point cannot be calculated
713 if ( math::equals(lon1_minus_lon2, c0)
714 || math::equals(lon1_minus_lon2 + (lon1_minus_lon2 < c0 ? pi : -pi), c0) )
715 {
716 return false;
717 }
718
719 lon = c0;
720
721 // the geodesics intersect at one of the poles
722 CT const pi_half = pi / CT(2);
723 CT const abs_lat_a1 = math::abs(lat_a1);
724 CT const abs_lat_a2 = math::abs(lat_a2);
725 if (math::equals(abs_lat_a1, abs_lat_a2))
726 {
727 lat = pi_half;
728 }
729 else
730 {
731 // pick the pole closest to one of the points of the first segment
732 CT const& closer_lat = abs_lat_a1 > abs_lat_a2 ? lat_a1 : lat_a2;
733 lat = closer_lat >= 0 ? pi_half : -pi_half;
734 }
735
736 return true;
737 }
738
739 CT lon_sph = 0;
740
741 // Starting tan(beta)
742 CT t = 0;
743
744 /*if (geod1.is_Cj_zero)
745 {
746 CT const k_base = lon1_minus_lon2 + geod2.sign_lon_diff * geod2.asin_tj_t0j;
747 t = sin(k_base) * geod2.t0j;
748 lon_sph = vertical_intersection_longitude(geod1.lonj, lon_b1, lon_b2);
749 }
750 else if (geod2.is_Cj_zero)
751 {
752 CT const k_base = lon1_minus_lon2 - geod1.sign_lon_diff * geod1.asin_tj_t0j;
753 t = sin(-k_base) * geod1.t0j;
754 lon_sph = vertical_intersection_longitude(geod2.lonj, lon_a1, lon_a2);
755 }
756 else*/
757 {
758 // TODO: Consider using betas instead of latitudes.
759 // Some function calls might be saved this way.
760 CT tan_lat_sph = 0;
761 sjoberg_intersection_spherical_02<CT>::apply_alt(lon_a1, lat_a1, lon_a2, lat_a2,
762 lon_b1, lat_b1, lon_b2, lat_b2,
763 lon_sph, tan_lat_sph);
764 t = one_minus_f * tan_lat_sph; // tan(beta)
765 }
766
767 // TODO: no need to calculate atan here if reduced latitudes were used
768 // instead of latitudes above, in sjoberg_intersection_spherical_02
769 CT const beta = atan(t);
770
771 if (enable_02 && newton_method(geod1, geod2, beta, t, lon1_minus_lon2, lon, lat))
772 {
773 return true;
774 }
775
776 return converge_07(geod1, geod2, beta, t, lon1_minus_lon2, lon_sph, lon, lat);
777 }
778
779 private:
newton_method(geodesic_type const & geod1,geodesic_type const & geod2,CT beta,CT t,CT const & lon1_minus_lon2,CT & lon,CT & lat)780 static inline bool newton_method(geodesic_type const& geod1, geodesic_type const& geod2, // in
781 CT beta, CT t, CT const& lon1_minus_lon2, // in
782 CT & lon, CT & lat) // out
783 {
784 CT const c0 = 0;
785 CT const c1 = 1;
786
787 CT const e_sqr = geod1.e_sqr;
788
789 CT lon1_diff = 0;
790 CT lon2_diff = 0;
791
792 CT abs_dbeta_last = 0;
793
794 // [Sjoberg02] converges faster than solution in [Sjoberg07]
795 // Newton-Raphson method
796 for (int i = 0; i < max_iterations_02; ++i)
797 {
798 CT const sin_beta = sin(beta);
799 CT const cos_beta = cos(beta);
800 CT const cos_beta_sqr = math::sqr(cos_beta);
801 CT const G = c1 - e_sqr * cos_beta_sqr;
802
803 CT f1 = 0;
804 CT f2 = 0;
805
806 if (!geod1.is_Cj_zero)
807 {
808 bool is_beta_ok = geod1.lon_diff(sin_beta, t, lon1_diff);
809
810 if (is_beta_ok)
811 {
812 CT const H = cos_beta_sqr - geod1.Cj_sqr;
813 f1 = geod1.Cj / cos_beta * math::sqrt(G / H);
814 }
815 else
816 {
817 return false;
818 }
819 }
820
821 if (!geod2.is_Cj_zero)
822 {
823 bool is_beta_ok = geod2.lon_diff(sin_beta, t, lon2_diff);
824
825 if (is_beta_ok)
826 {
827 CT const H = cos_beta_sqr - geod2.Cj_sqr;
828 f2 = geod2.Cj / cos_beta * math::sqrt(G / H);
829 }
830 else
831 {
832 return false;
833 }
834 }
835
836 // NOTE: Things may go wrong if the IP is near the vertex
837 // 1. May converge into the wrong direction (from the other way around).
838 // This happens when the starting point is on the other side than the vertex
839 // 2. During converging may "jump" into the other side of the vertex.
840 // In this case sin_beta/sqrt_1_Cj_sqr and t/t0j is not in [-1, 1]
841 // 3. f1-f2 may be 0 which means that the intermediate point is on the vertex
842 // In this case it's not possible to check if this is the correct result
843
844 CT const dbeta_denom = f1 - f2;
845 //CT const dbeta_denom = math::abs(f1) + math::abs(f2);
846
847 if (math::equals(dbeta_denom, c0))
848 {
849 return false;
850 }
851
852 // The sign of dbeta is changed WRT [Sjoberg02]
853 CT const dbeta = (lon1_minus_lon2 + lon1_diff - lon2_diff) / dbeta_denom;
854
855 CT const abs_dbeta = math::abs(dbeta);
856 if (i > 0 && abs_dbeta > abs_dbeta_last)
857 {
858 // The algorithm is not converging
859 // The intersection may be on the other side of the vertex
860 return false;
861 }
862 abs_dbeta_last = abs_dbeta;
863
864 if (math::equals(dbeta, c0))
865 {
866 // Result found
867 break;
868 }
869
870 // Because the sign of dbeta is changed WRT [Sjoberg02] dbeta is subtracted here
871 beta = beta - dbeta;
872
873 t = tan(beta);
874 }
875
876 lat = geod1.lat(t);
877 // NOTE: if Cj is 0 then the result is lonj or lonj+180
878 lon = ! geod1.is_Cj_zero
879 ? geod1.lon(lon1_diff)
880 : geod2.lon(lon2_diff);
881
882 return true;
883 }
884
885 struct geodesics_type
886 {
geodesics_typeboost::geometry::formula::sjoberg_intersection::geodesics_type887 geodesics_type(geodesic_type const& g1, geodesic_type const& g2)
888 : geod1(g1)
889 , geod2(g2)
890 , vertex1(geod1.get_vertex_data())
891 , vertex2(geod2.get_vertex_data())
892 {}
893
894 geodesic_type const& geod1;
895 geodesic_type const& geod2;
896 typename geodesic_type::vertex_data vertex1;
897 typename geodesic_type::vertex_data vertex2;
898 };
899
900 struct converge_07_result
901 {
converge_07_resultboost::geometry::formula::sjoberg_intersection::converge_07_result902 converge_07_result()
903 : lon1(0), lon2(0), k1_diff(0), k2_diff(0), t1(0), t2(0)
904 {}
905
906 CT lon1, lon2;
907 CT k1_diff, k2_diff;
908 CT t1, t2;
909 };
910
converge_07(geodesic_type const & geod1,geodesic_type const & geod2,CT beta,CT t,CT const & lon1_minus_lon2,CT const & lon_sph,CT & lon,CT & lat)911 static inline bool converge_07(geodesic_type const& geod1, geodesic_type const& geod2,
912 CT beta, CT t,
913 CT const& lon1_minus_lon2, CT const& lon_sph,
914 CT & lon, CT & lat)
915 {
916 //CT const c0 = 0;
917 //CT const c1 = 1;
918 //CT const c2 = 2;
919 //CT const pi = math::pi<CT>();
920
921 geodesics_type geodesics(geod1, geod2);
922 converge_07_result result;
923
924 // calculate first pair of longitudes
925 if (!converge_07_step_one(CT(sin(beta)), t, lon1_minus_lon2, geodesics, lon_sph, result, false))
926 {
927 return false;
928 }
929
930 int t_direction = 0;
931
932 CT lon_diff_prev = math::longitude_difference<radian>(result.lon1, result.lon2);
933
934 // [Sjoberg07]
935 for (int i = 2; i < max_iterations_07; ++i)
936 {
937 // pick t candidates from previous result based on dir
938 CT t_cand1 = result.t1;
939 CT t_cand2 = result.t2;
940 // if direction is 0 the closer one is the first
941 if (t_direction < 0)
942 {
943 t_cand1 = (std::min)(result.t1, result.t2);
944 t_cand2 = (std::max)(result.t1, result.t2);
945 }
946 else if (t_direction > 0)
947 {
948 t_cand1 = (std::max)(result.t1, result.t2);
949 t_cand2 = (std::min)(result.t1, result.t2);
950 }
951 else
952 {
953 t_direction = t_cand1 < t_cand2 ? -1 : 1;
954 }
955
956 CT t1 = t;
957 CT beta1 = beta;
958 // check if the further calculation is needed
959 if (converge_07_update(t1, beta1, t_cand1))
960 {
961 break;
962 }
963
964 bool try_t2 = false;
965 converge_07_result result_curr;
966 if (converge_07_step_one(CT(sin(beta1)), t1, lon1_minus_lon2, geodesics, lon_sph, result_curr))
967 {
968 CT const lon_diff1 = math::longitude_difference<radian>(result_curr.lon1, result_curr.lon2);
969 if (lon_diff_prev > lon_diff1)
970 {
971 t = t1;
972 beta = beta1;
973 lon_diff_prev = lon_diff1;
974 result = result_curr;
975 }
976 else if (t_cand1 != t_cand2)
977 {
978 try_t2 = true;
979 }
980 else
981 {
982 // the result is not fully correct but it won't be more accurate
983 break;
984 }
985 }
986 // ! converge_07_step_one
987 else
988 {
989 if (t_cand1 != t_cand2)
990 {
991 try_t2 = true;
992 }
993 else
994 {
995 return false;
996 }
997 }
998
999
1000 if (try_t2)
1001 {
1002 CT t2 = t;
1003 CT beta2 = beta;
1004 // check if the further calculation is needed
1005 if (converge_07_update(t2, beta2, t_cand2))
1006 {
1007 break;
1008 }
1009
1010 if (! converge_07_step_one(CT(sin(beta2)), t2, lon1_minus_lon2, geodesics, lon_sph, result_curr))
1011 {
1012 return false;
1013 }
1014
1015 CT const lon_diff2 = math::longitude_difference<radian>(result_curr.lon1, result_curr.lon2);
1016 if (lon_diff_prev > lon_diff2)
1017 {
1018 t_direction *= -1;
1019 t = t2;
1020 beta = beta2;
1021 lon_diff_prev = lon_diff2;
1022 result = result_curr;
1023 }
1024 else
1025 {
1026 // the result is not fully correct but it won't be more accurate
1027 break;
1028 }
1029 }
1030 }
1031
1032 lat = geod1.lat(t);
1033 lon = ! geod1.is_Cj_zero ? result.lon1 : result.lon2;
1034 math::normalize_longitude<radian>(lon);
1035
1036 return true;
1037 }
1038
converge_07_update(CT & t,CT & beta,CT const & t_new)1039 static inline bool converge_07_update(CT & t, CT & beta, CT const& t_new)
1040 {
1041 CT const c0 = 0;
1042
1043 CT const beta_new = atan(t_new);
1044 CT const dbeta = beta_new - beta;
1045 beta = beta_new;
1046 t = t_new;
1047
1048 return math::equals(dbeta, c0);
1049 }
1050
pick_t(CT const & t1,CT const & t2,int direction)1051 static inline CT const& pick_t(CT const& t1, CT const& t2, int direction)
1052 {
1053 return direction < 0 ? (std::min)(t1, t2) : (std::max)(t1, t2);
1054 }
1055
converge_07_step_one(CT const & sin_beta,CT const & t,CT const & lon1_minus_lon2,geodesics_type const & geodesics,CT const & lon_sph,converge_07_result & result,bool check_sin_beta=true)1056 static inline bool converge_07_step_one(CT const& sin_beta,
1057 CT const& t,
1058 CT const& lon1_minus_lon2,
1059 geodesics_type const& geodesics,
1060 CT const& lon_sph,
1061 converge_07_result & result,
1062 bool check_sin_beta = true)
1063 {
1064 bool ok = converge_07_one_geod(sin_beta, t, geodesics.geod1, geodesics.vertex1, lon_sph,
1065 result.lon1, result.k1_diff, check_sin_beta)
1066 && converge_07_one_geod(sin_beta, t, geodesics.geod2, geodesics.vertex2, lon_sph,
1067 result.lon2, result.k2_diff, check_sin_beta);
1068
1069 if (!ok)
1070 {
1071 return false;
1072 }
1073
1074 CT const k = lon1_minus_lon2 + result.k1_diff - result.k2_diff;
1075
1076 // get 2 possible ts one lesser and one greater than t
1077 // t1 is the closer one
1078 calc_ts(t, k, geodesics.geod1, geodesics.geod2, result.t1, result.t2);
1079
1080 return true;
1081 }
1082
converge_07_one_geod(CT const & sin_beta,CT const & t,geodesic_type const & geod,typename geodesic_type::vertex_data const & vertex,CT const & lon_sph,CT & lon,CT & k_diff,bool check_sin_beta)1083 static inline bool converge_07_one_geod(CT const& sin_beta, CT const& t,
1084 geodesic_type const& geod,
1085 typename geodesic_type::vertex_data const& vertex,
1086 CT const& lon_sph,
1087 CT & lon, CT & k_diff,
1088 bool check_sin_beta)
1089 {
1090 using math::detail::bounded;
1091 CT const c1 = 1;
1092
1093 CT k_diff_before = 0;
1094 CT k_diff_behind = 0;
1095
1096 bool is_beta_ok = geod.k_diffs(sin_beta, vertex, k_diff_before, k_diff_behind, check_sin_beta);
1097
1098 if (! is_beta_ok)
1099 {
1100 return false;
1101 }
1102
1103 CT const asin_t_t0j = ! geod.is_Cj_zero ? asin(bounded(t / geod.t0j, -c1, c1)) : 0;
1104 CT const sign_asin_t_t0j = geod.sign_lon_diff * asin_t_t0j;
1105
1106 CT const lon_before = geod.lonj + sign_asin_t_t0j + k_diff_before;
1107 CT const lon_behind = geod.lonj - sign_asin_t_t0j + k_diff_behind;
1108
1109 CT const lon_dist_before = math::longitude_distance_signed<radian>(lon_before, lon_sph);
1110 CT const lon_dist_behind = math::longitude_distance_signed<radian>(lon_behind, lon_sph);
1111 if (math::abs(lon_dist_before) <= math::abs(lon_dist_behind))
1112 {
1113 k_diff = k_diff_before;
1114 lon = lon_before;
1115 }
1116 else
1117 {
1118 k_diff = k_diff_behind;
1119 lon = lon_behind;
1120 }
1121
1122 return true;
1123 }
1124
calc_ts(CT const & t,CT const & k,geodesic_type const & geod1,geodesic_type const & geod2,CT & t1,CT & t2)1125 static inline void calc_ts(CT const& t, CT const& k,
1126 geodesic_type const& geod1, geodesic_type const& geod2,
1127 CT & t1, CT& t2)
1128 {
1129 CT const c1 = 1;
1130 CT const c2 = 2;
1131
1132 CT const K = sin(k);
1133
1134 BOOST_GEOMETRY_ASSERT(!geod1.is_Cj_zero || !geod2.is_Cj_zero);
1135 if (geod1.is_Cj_zero)
1136 {
1137 t1 = K * geod2.t0j;
1138 t2 = -t1;
1139 }
1140 else if (geod2.is_Cj_zero)
1141 {
1142 t1 = -K * geod1.t0j;
1143 t2 = -t1;
1144 }
1145 else
1146 {
1147 CT const A = math::sqr(geod1.t0j) + math::sqr(geod2.t0j);
1148 CT const B = c2 * geod1.t0j * geod2.t0j * math::sqrt(c1 - math::sqr(K));
1149
1150 CT const K_t01_t02 = K * geod1.t0j * geod2.t0j;
1151 CT const D1 = math::sqrt(A + B);
1152 CT const D2 = math::sqrt(A - B);
1153 CT const t_new1 = K_t01_t02 / D1;
1154 CT const t_new2 = K_t01_t02 / D2;
1155 CT const t_new3 = -t_new1;
1156 CT const t_new4 = -t_new2;
1157
1158 // Pick 2 nearest t_new, one greater and one lesser than current t
1159 CT const abs_t_new1 = math::abs(t_new1);
1160 CT const abs_t_new2 = math::abs(t_new2);
1161 CT const abs_t_max = (std::max)(abs_t_new1, abs_t_new2);
1162 t1 = -abs_t_max; // lesser
1163 t2 = abs_t_max; // greater
1164 if (t1 < t)
1165 {
1166 if (t_new1 < t && t_new1 > t1)
1167 t1 = t_new1;
1168 if (t_new2 < t && t_new2 > t1)
1169 t1 = t_new2;
1170 if (t_new3 < t && t_new3 > t1)
1171 t1 = t_new3;
1172 if (t_new4 < t && t_new4 > t1)
1173 t1 = t_new4;
1174 }
1175 if (t2 > t)
1176 {
1177 if (t_new1 > t && t_new1 < t2)
1178 t2 = t_new1;
1179 if (t_new2 > t && t_new2 < t2)
1180 t2 = t_new2;
1181 if (t_new3 > t && t_new3 < t2)
1182 t2 = t_new3;
1183 if (t_new4 > t && t_new4 < t2)
1184 t2 = t_new4;
1185 }
1186 }
1187
1188 // the first one is the closer one
1189 if (math::abs(t - t2) < math::abs(t - t1))
1190 {
1191 std::swap(t2, t1);
1192 }
1193 }
1194
fj(CT const & cos_beta,CT const & cos2_beta,CT const & Cj,CT const & e_sqr)1195 static inline CT fj(CT const& cos_beta, CT const& cos2_beta, CT const& Cj, CT const& e_sqr)
1196 {
1197 CT const c1 = 1;
1198 CT const Cj_sqr = math::sqr(Cj);
1199 return Cj / cos_beta * math::sqrt((c1 - e_sqr * cos2_beta) / (cos2_beta - Cj_sqr));
1200 }
1201
1202 /*static inline CT vertical_intersection_longitude(CT const& ip_lon, CT const& seg_lon1, CT const& seg_lon2)
1203 {
1204 CT const c0 = 0;
1205 CT const lon_2 = ip_lon > c0 ? ip_lon - pi : ip_lon + pi;
1206
1207 return (std::min)(math::longitude_difference<radian>(ip_lon, seg_lon1),
1208 math::longitude_difference<radian>(ip_lon, seg_lon2))
1209 <=
1210 (std::min)(math::longitude_difference<radian>(lon_2, seg_lon1),
1211 math::longitude_difference<radian>(lon_2, seg_lon2))
1212 ? ip_lon : lon_2;
1213 }*/
1214 };
1215
1216 }}} // namespace boost::geometry::formula
1217
1218
1219 #endif // BOOST_GEOMETRY_FORMULAS_SJOBERG_INTERSECTION_HPP
1220