1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
// Boost.Geometry
// Copyright (c) 2016 Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP
#define BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/arithmetic/cross_product.hpp>
#include <boost/geometry/formulas/gnomonic_spheroid.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/util/math.hpp>
namespace boost { namespace geometry { namespace formula
{
/*!
\brief The intersection of two geodesics using spheroidal gnomonic projection
as proposed by Karney.
\author See
- Charles F.F Karney, Algorithms for geodesics, 2011
https://arxiv.org/pdf/1109.4448.pdf
- GeographicLib forum thread: Intersection between two geodesic lines
https://sourceforge.net/p/geographiclib/discussion/1026621/thread/21aaff9f/
*/
template
<
typename CT,
template <typename, bool, bool, bool, bool, bool> class Inverse,
template <typename, bool, bool, bool, bool> class Direct
>
class gnomonic_intersection
{
public:
template <typename T1, typename T2, typename Spheroid>
static inline bool 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)
{
CT const lon_a1 = lona1;
CT const lat_a1 = lata1;
CT const lon_a2 = lona2;
CT const lat_a2 = lata2;
CT const lon_b1 = lonb1;
CT const lat_b1 = latb1;
CT const lon_b2 = lonb2;
CT const lat_b2 = latb2;
return apply(lon_a1, lat_a1, lon_a2, lat_a2, lon_b1, lat_b1, lon_b2, lat_b2, lon, lat, spheroid);
}
template <typename Spheroid>
static inline bool apply(CT const& lona1, CT const& lata1,
CT const& lona2, CT const& lata2,
CT const& lonb1, CT const& latb1,
CT const& lonb2, CT const& latb2,
CT & lon, CT & lat,
Spheroid const& spheroid)
{
typedef gnomonic_spheroid<CT, Inverse, Direct> gnom_t;
lon = (lona1 + lona2 + lonb1 + lonb2) / 4;
lat = (lata1 + lata2 + latb1 + latb2) / 4;
// TODO: consider normalizing lon
for (int i = 0; i < 10; ++i)
{
CT xa1, ya1, xa2, ya2;
CT xb1, yb1, xb2, yb2;
CT x, y;
double lat1, lon1;
bool ok = gnom_t::forward(lon, lat, lona1, lata1, xa1, ya1, spheroid)
&& gnom_t::forward(lon, lat, lona2, lata2, xa2, ya2, spheroid)
&& gnom_t::forward(lon, lat, lonb1, latb1, xb1, yb1, spheroid)
&& gnom_t::forward(lon, lat, lonb2, latb2, xb2, yb2, spheroid)
&& intersect(xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2, x, y)
&& gnom_t::inverse(lon, lat, x, y, lon1, lat1, spheroid);
if (! ok)
{
return false;
}
if (math::equals(lat1, lat) && math::equals(lon1, lon))
{
break;
}
lat = lat1;
lon = lon1;
}
// NOTE: true is also returned if the number of iterations is too great
// which means that the accuracy of the result is low
return true;
}
private:
static inline bool intersect(CT const& xa1, CT const& ya1, CT const& xa2, CT const& ya2,
CT const& xb1, CT const& yb1, CT const& xb2, CT const& yb2,
CT & x, CT & y)
{
typedef model::point<CT, 3, cs::cartesian> v3d_t;
CT const c0 = 0;
CT const c1 = 1;
v3d_t const va1(xa1, ya1, c1);
v3d_t const va2(xa2, ya2, c1);
v3d_t const vb1(xb1, yb1, c1);
v3d_t const vb2(xb2, yb2, c1);
v3d_t const la = cross_product(va1, va2);
v3d_t const lb = cross_product(vb1, vb2);
v3d_t const p = cross_product(la, lb);
CT const z = get<2>(p);
if (math::equals(z, c0))
{
// degenerated or collinear segments
return false;
}
x = get<0>(p) / z;
y = get<1>(p) / z;
return true;
}
};
}}} // namespace boost::geometry::formula
#endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP
|
