Chris@16
|
1 // Boost.Geometry (aka GGL, Generic Geometry Library)
|
Chris@16
|
2
|
Chris@16
|
3 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
|
Chris@16
|
4
|
Chris@16
|
5 // Use, modification and distribution is subject to the Boost Software License,
|
Chris@16
|
6 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
|
Chris@16
|
7 // http://www.boost.org/LICENSE_1_0.txt)
|
Chris@16
|
8
|
Chris@16
|
9 #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP
|
Chris@16
|
10 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP
|
Chris@16
|
11
|
Chris@16
|
12
|
Chris@16
|
13 #include <boost/config.hpp>
|
Chris@16
|
14 #include <boost/concept_check.hpp>
|
Chris@16
|
15 #include <boost/mpl/if.hpp>
|
Chris@16
|
16 #include <boost/type_traits.hpp>
|
Chris@16
|
17
|
Chris@16
|
18 #include <boost/geometry/core/cs.hpp>
|
Chris@16
|
19 #include <boost/geometry/core/access.hpp>
|
Chris@16
|
20 #include <boost/geometry/core/radian_access.hpp>
|
Chris@16
|
21
|
Chris@16
|
22
|
Chris@16
|
23 #include <boost/geometry/strategies/distance.hpp>
|
Chris@16
|
24 #include <boost/geometry/strategies/concepts/distance_concept.hpp>
|
Chris@16
|
25 #include <boost/geometry/strategies/spherical/distance_haversine.hpp>
|
Chris@16
|
26
|
Chris@16
|
27 #include <boost/geometry/util/promote_floating_point.hpp>
|
Chris@16
|
28 #include <boost/geometry/util/math.hpp>
|
Chris@16
|
29
|
Chris@16
|
30 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
|
Chris@16
|
31 # include <boost/geometry/io/dsv/write.hpp>
|
Chris@16
|
32 #endif
|
Chris@16
|
33
|
Chris@16
|
34
|
Chris@16
|
35
|
Chris@16
|
36 namespace boost { namespace geometry
|
Chris@16
|
37 {
|
Chris@16
|
38
|
Chris@16
|
39 namespace strategy { namespace distance
|
Chris@16
|
40 {
|
Chris@16
|
41
|
Chris@16
|
42 /*!
|
Chris@16
|
43 \brief Strategy functor for distance point to segment calculation
|
Chris@16
|
44 \ingroup strategies
|
Chris@16
|
45 \details Class which calculates the distance of a point to a segment, for points on a sphere or globe
|
Chris@16
|
46 \see http://williams.best.vwh.net/avform.htm
|
Chris@16
|
47 \tparam CalculationType \tparam_calculation
|
Chris@16
|
48 \tparam Strategy underlying point-point distance strategy, defaults to haversine
|
Chris@16
|
49
|
Chris@16
|
50 \qbk{
|
Chris@16
|
51 [heading See also]
|
Chris@16
|
52 [link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)]
|
Chris@16
|
53 }
|
Chris@16
|
54
|
Chris@16
|
55 */
|
Chris@16
|
56 template
|
Chris@16
|
57 <
|
Chris@16
|
58 typename CalculationType = void,
|
Chris@16
|
59 typename Strategy = haversine<double, CalculationType>
|
Chris@16
|
60 >
|
Chris@16
|
61 class cross_track
|
Chris@16
|
62 {
|
Chris@16
|
63 public :
|
Chris@16
|
64 template <typename Point, typename PointOfSegment>
|
Chris@16
|
65 struct return_type
|
Chris@16
|
66 : promote_floating_point
|
Chris@16
|
67 <
|
Chris@16
|
68 typename select_calculation_type
|
Chris@16
|
69 <
|
Chris@16
|
70 Point,
|
Chris@16
|
71 PointOfSegment,
|
Chris@16
|
72 CalculationType
|
Chris@16
|
73 >::type
|
Chris@16
|
74 >
|
Chris@16
|
75 {};
|
Chris@16
|
76
|
Chris@16
|
77 inline cross_track()
|
Chris@16
|
78 {}
|
Chris@16
|
79
|
Chris@16
|
80 explicit inline cross_track(typename Strategy::radius_type const& r)
|
Chris@16
|
81 : m_strategy(r)
|
Chris@16
|
82 {}
|
Chris@16
|
83
|
Chris@16
|
84 inline cross_track(Strategy const& s)
|
Chris@16
|
85 : m_strategy(s)
|
Chris@16
|
86 {}
|
Chris@16
|
87
|
Chris@16
|
88
|
Chris@16
|
89 // It might be useful in the future
|
Chris@16
|
90 // to overload constructor with strategy info.
|
Chris@16
|
91 // crosstrack(...) {}
|
Chris@16
|
92
|
Chris@16
|
93
|
Chris@16
|
94 template <typename Point, typename PointOfSegment>
|
Chris@16
|
95 inline typename return_type<Point, PointOfSegment>::type
|
Chris@16
|
96 apply(Point const& p, PointOfSegment const& sp1, PointOfSegment const& sp2) const
|
Chris@16
|
97 {
|
Chris@16
|
98
|
Chris@16
|
99 #if !defined(BOOST_MSVC)
|
Chris@16
|
100 BOOST_CONCEPT_ASSERT
|
Chris@16
|
101 (
|
Chris@16
|
102 (concept::PointDistanceStrategy<Strategy, Point, PointOfSegment>)
|
Chris@16
|
103 );
|
Chris@16
|
104 #endif
|
Chris@16
|
105
|
Chris@16
|
106 typedef typename return_type<Point, PointOfSegment>::type return_type;
|
Chris@16
|
107
|
Chris@16
|
108 // http://williams.best.vwh.net/avform.htm#XTE
|
Chris@16
|
109 return_type d1 = m_strategy.apply(sp1, p);
|
Chris@16
|
110 return_type d3 = m_strategy.apply(sp1, sp2);
|
Chris@16
|
111
|
Chris@16
|
112 if (geometry::math::equals(d3, 0.0))
|
Chris@16
|
113 {
|
Chris@16
|
114 // "Degenerate" segment, return either d1 or d2
|
Chris@16
|
115 return d1;
|
Chris@16
|
116 }
|
Chris@16
|
117
|
Chris@16
|
118 return_type d2 = m_strategy.apply(sp2, p);
|
Chris@16
|
119
|
Chris@16
|
120 return_type crs_AD = course(sp1, p);
|
Chris@16
|
121 return_type crs_AB = course(sp1, sp2);
|
Chris@16
|
122 return_type crs_BA = crs_AB - geometry::math::pi<return_type>();
|
Chris@16
|
123 return_type crs_BD = course(sp2, p);
|
Chris@16
|
124 return_type d_crs1 = crs_AD - crs_AB;
|
Chris@16
|
125 return_type d_crs2 = crs_BD - crs_BA;
|
Chris@16
|
126
|
Chris@16
|
127 // d1, d2, d3 are in principle not needed, only the sign matters
|
Chris@16
|
128 return_type projection1 = cos( d_crs1 ) * d1 / d3;
|
Chris@16
|
129 return_type projection2 = cos( d_crs2 ) * d2 / d3;
|
Chris@16
|
130
|
Chris@16
|
131 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
|
Chris@16
|
132 std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " " << crs_AD * geometry::math::r2d << std::endl;
|
Chris@16
|
133 std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " " << crs_AB * geometry::math::r2d << std::endl;
|
Chris@16
|
134 std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " " << crs_BD * geometry::math::r2d << std::endl;
|
Chris@16
|
135 std::cout << "Projection AD-AB " << projection1 << " : " << d_crs1 * geometry::math::r2d << std::endl;
|
Chris@16
|
136 std::cout << "Projection BD-BA " << projection2 << " : " << d_crs2 * geometry::math::r2d << std::endl;
|
Chris@16
|
137 #endif
|
Chris@16
|
138
|
Chris@16
|
139 if(projection1 > 0.0 && projection2 > 0.0)
|
Chris@16
|
140 {
|
Chris@16
|
141 return_type XTD = radius() * geometry::math::abs( asin( sin( d1 / radius() ) * sin( d_crs1 ) ));
|
Chris@16
|
142
|
Chris@16
|
143 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
|
Chris@16
|
144 std::cout << "Projection ON the segment" << std::endl;
|
Chris@16
|
145 std::cout << "XTD: " << XTD << " d1: " << d1 << " d2: " << d2 << std::endl;
|
Chris@16
|
146 #endif
|
Chris@16
|
147
|
Chris@16
|
148 // Return shortest distance, projected point on segment sp1-sp2
|
Chris@16
|
149 return return_type(XTD);
|
Chris@16
|
150 }
|
Chris@16
|
151 else
|
Chris@16
|
152 {
|
Chris@16
|
153 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
|
Chris@16
|
154 std::cout << "Projection OUTSIDE the segment" << std::endl;
|
Chris@16
|
155 #endif
|
Chris@16
|
156
|
Chris@16
|
157 // Return shortest distance, project either on point sp1 or sp2
|
Chris@16
|
158 return return_type( (std::min)( d1 , d2 ) );
|
Chris@16
|
159 }
|
Chris@16
|
160 }
|
Chris@16
|
161
|
Chris@16
|
162 inline typename Strategy::radius_type radius() const
|
Chris@16
|
163 { return m_strategy.radius(); }
|
Chris@16
|
164
|
Chris@16
|
165 private :
|
Chris@16
|
166
|
Chris@16
|
167 Strategy m_strategy;
|
Chris@16
|
168
|
Chris@16
|
169 /// Calculate course (bearing) between two points. Might be moved to a "course formula" ...
|
Chris@16
|
170 template <typename Point1, typename Point2>
|
Chris@16
|
171 inline typename return_type<Point1, Point2>::type
|
Chris@16
|
172 course(Point1 const& p1, Point2 const& p2) const
|
Chris@16
|
173 {
|
Chris@16
|
174 typedef typename return_type<Point1, Point2>::type return_type;
|
Chris@16
|
175
|
Chris@16
|
176 // http://williams.best.vwh.net/avform.htm#Crs
|
Chris@16
|
177 return_type dlon = get_as_radian<0>(p2) - get_as_radian<0>(p1);
|
Chris@16
|
178 return_type cos_p2lat = cos(get_as_radian<1>(p2));
|
Chris@16
|
179
|
Chris@16
|
180 // "An alternative formula, not requiring the pre-computation of d"
|
Chris@16
|
181 return atan2(sin(dlon) * cos_p2lat,
|
Chris@16
|
182 cos(get_as_radian<1>(p1)) * sin(get_as_radian<1>(p2))
|
Chris@16
|
183 - sin(get_as_radian<1>(p1)) * cos_p2lat * cos(dlon));
|
Chris@16
|
184 }
|
Chris@16
|
185
|
Chris@16
|
186 };
|
Chris@16
|
187
|
Chris@16
|
188
|
Chris@16
|
189
|
Chris@16
|
190 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
|
Chris@16
|
191 namespace services
|
Chris@16
|
192 {
|
Chris@16
|
193
|
Chris@16
|
194 template <typename CalculationType, typename Strategy>
|
Chris@16
|
195 struct tag<cross_track<CalculationType, Strategy> >
|
Chris@16
|
196 {
|
Chris@16
|
197 typedef strategy_tag_distance_point_segment type;
|
Chris@16
|
198 };
|
Chris@16
|
199
|
Chris@16
|
200
|
Chris@16
|
201 template <typename CalculationType, typename Strategy, typename P, typename PS>
|
Chris@16
|
202 struct return_type<cross_track<CalculationType, Strategy>, P, PS>
|
Chris@16
|
203 : cross_track<CalculationType, Strategy>::template return_type<P, PS>
|
Chris@16
|
204 {};
|
Chris@16
|
205
|
Chris@16
|
206
|
Chris@16
|
207 template <typename CalculationType, typename Strategy>
|
Chris@16
|
208 struct comparable_type<cross_track<CalculationType, Strategy> >
|
Chris@16
|
209 {
|
Chris@16
|
210 // There is no shortcut, so the strategy itself is its comparable type
|
Chris@16
|
211 typedef cross_track<CalculationType, Strategy> type;
|
Chris@16
|
212 };
|
Chris@16
|
213
|
Chris@16
|
214
|
Chris@16
|
215 template
|
Chris@16
|
216 <
|
Chris@16
|
217 typename CalculationType,
|
Chris@16
|
218 typename Strategy
|
Chris@16
|
219 >
|
Chris@16
|
220 struct get_comparable<cross_track<CalculationType, Strategy> >
|
Chris@16
|
221 {
|
Chris@16
|
222 typedef typename comparable_type
|
Chris@16
|
223 <
|
Chris@16
|
224 cross_track<CalculationType, Strategy>
|
Chris@16
|
225 >::type comparable_type;
|
Chris@16
|
226 public :
|
Chris@16
|
227 static inline comparable_type apply(cross_track<CalculationType, Strategy> const& strategy)
|
Chris@16
|
228 {
|
Chris@16
|
229 return comparable_type(strategy.radius());
|
Chris@16
|
230 }
|
Chris@16
|
231 };
|
Chris@16
|
232
|
Chris@16
|
233
|
Chris@16
|
234 template
|
Chris@16
|
235 <
|
Chris@16
|
236 typename CalculationType,
|
Chris@16
|
237 typename Strategy,
|
Chris@16
|
238 typename P, typename PS
|
Chris@16
|
239 >
|
Chris@16
|
240 struct result_from_distance<cross_track<CalculationType, Strategy>, P, PS>
|
Chris@16
|
241 {
|
Chris@16
|
242 private :
|
Chris@16
|
243 typedef typename cross_track<CalculationType, Strategy>::template return_type<P, PS> return_type;
|
Chris@16
|
244 public :
|
Chris@16
|
245 template <typename T>
|
Chris@16
|
246 static inline return_type apply(cross_track<CalculationType, Strategy> const& , T const& distance)
|
Chris@16
|
247 {
|
Chris@16
|
248 return distance;
|
Chris@16
|
249 }
|
Chris@16
|
250 };
|
Chris@16
|
251
|
Chris@16
|
252
|
Chris@16
|
253 template
|
Chris@16
|
254 <
|
Chris@16
|
255 typename CalculationType,
|
Chris@16
|
256 typename Strategy
|
Chris@16
|
257 >
|
Chris@16
|
258 struct strategy_point_point<cross_track<CalculationType, Strategy> >
|
Chris@16
|
259 {
|
Chris@16
|
260 typedef Strategy type;
|
Chris@16
|
261 };
|
Chris@16
|
262
|
Chris@16
|
263
|
Chris@16
|
264
|
Chris@16
|
265 /*
|
Chris@16
|
266
|
Chris@16
|
267 TODO: spherical polar coordinate system requires "get_as_radian_equatorial<>"
|
Chris@16
|
268
|
Chris@16
|
269 template <typename Point, typename PointOfSegment, typename Strategy>
|
Chris@16
|
270 struct default_strategy
|
Chris@16
|
271 <
|
Chris@16
|
272 segment_tag, Point, PointOfSegment,
|
Chris@16
|
273 spherical_polar_tag, spherical_polar_tag,
|
Chris@16
|
274 Strategy
|
Chris@16
|
275 >
|
Chris@16
|
276 {
|
Chris@16
|
277 typedef cross_track
|
Chris@16
|
278 <
|
Chris@16
|
279 void,
|
Chris@16
|
280 typename boost::mpl::if_
|
Chris@16
|
281 <
|
Chris@16
|
282 boost::is_void<Strategy>,
|
Chris@16
|
283 typename default_strategy
|
Chris@16
|
284 <
|
Chris@16
|
285 point_tag, Point, PointOfSegment,
|
Chris@16
|
286 spherical_polar_tag, spherical_polar_tag
|
Chris@16
|
287 >::type,
|
Chris@16
|
288 Strategy
|
Chris@16
|
289 >::type
|
Chris@16
|
290 > type;
|
Chris@16
|
291 };
|
Chris@16
|
292 */
|
Chris@16
|
293
|
Chris@16
|
294 template <typename Point, typename PointOfSegment, typename Strategy>
|
Chris@16
|
295 struct default_strategy
|
Chris@16
|
296 <
|
Chris@16
|
297 segment_tag, Point, PointOfSegment,
|
Chris@16
|
298 spherical_equatorial_tag, spherical_equatorial_tag,
|
Chris@16
|
299 Strategy
|
Chris@16
|
300 >
|
Chris@16
|
301 {
|
Chris@16
|
302 typedef cross_track
|
Chris@16
|
303 <
|
Chris@16
|
304 void,
|
Chris@16
|
305 typename boost::mpl::if_
|
Chris@16
|
306 <
|
Chris@16
|
307 boost::is_void<Strategy>,
|
Chris@16
|
308 typename default_strategy
|
Chris@16
|
309 <
|
Chris@16
|
310 point_tag, Point, PointOfSegment,
|
Chris@16
|
311 spherical_equatorial_tag, spherical_equatorial_tag
|
Chris@16
|
312 >::type,
|
Chris@16
|
313 Strategy
|
Chris@16
|
314 >::type
|
Chris@16
|
315 > type;
|
Chris@16
|
316 };
|
Chris@16
|
317
|
Chris@16
|
318
|
Chris@16
|
319
|
Chris@16
|
320 } // namespace services
|
Chris@16
|
321 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
|
Chris@16
|
322
|
Chris@16
|
323
|
Chris@16
|
324 }} // namespace strategy::distance
|
Chris@16
|
325
|
Chris@16
|
326
|
Chris@16
|
327 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
|
Chris@16
|
328
|
Chris@16
|
329
|
Chris@16
|
330 #endif
|
Chris@16
|
331
|
Chris@16
|
332
|
Chris@16
|
333 }} // namespace boost::geometry
|
Chris@16
|
334
|
Chris@16
|
335
|
Chris@16
|
336 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP
|