← Back to branch summary

~ubuntu-branches/ubuntu/saucy/merkaartor/saucy

« back to all changes in this revision

Viewing changes to include/ggl/strategies/spherical/sph_distance.hpp

  • Committer: Bazaar Package Importer
  • Author(s): Bernd Zeimetz
  • Date: 2009-09-13 00:52:12 UTC
  • mto: (1.2.7 upstream) (0.1.3 upstream) (3.1.7 sid)
  • mto: This revision was merged to the branch mainline in revision 10.
  • Revision ID: james.westby@ubuntu.com-20090913005212-pjecal8zxm07x0fj
Import upstream version 0.14+svnfixes~20090912

Show diffs side-by-side

added added

removed removed

Lines of Context:
include/ggl/strategies/spherical/sph_distance.hpp
 
1
// Generic Geometry Library
 
2
//
 
3
// Copyright Barend Gehrels 1995-2009, Geodan Holding B.V. Amsterdam, the Netherlands.
 
4
// Copyright Bruno Lalande 2008, 2009
 
5
// Use, modification and distribution is subject to the Boost Software License,
 
6
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
 
7
// http://www.boost.org/LICENSE_1_0.txt)
 
8
 
 
9
#ifndef GGL_STRATEGY_SPHERICAL_DISTANCE_HPP
 
10
#define GGL_STRATEGY_SPHERICAL_DISTANCE_HPP
 
11
 
 
12
 
 
13
#include <ggl/core/cs.hpp>
 
14
#include <ggl/core/concepts/point_concept.hpp>
 
15
#include <ggl/core/access.hpp>
 
16
#include <ggl/core/radian_access.hpp>
 
17
 
 
18
 
 
19
#include <ggl/strategies/strategy_traits.hpp>
 
20
 
 
21
 
 
22
#include <ggl/strategies/distance_result.hpp>
 
23
 
 
24
// Helper geometry
 
25
#include <ggl/geometries/point_ll.hpp> // to be resolved
 
26
 
 
27
 
 
28
namespace ggl
 
29
{
 
30
namespace strategy
 
31
{
 
32
 
 
33
    namespace distance
 
34
    {
 
35
 
 
36
        /*!
 
37
            \brief Distance calculation for spherical coordinates on a perfect sphere using haversine
 
38
            \ingroup distance
 
39
            \tparam P1 first point type
 
40
            \tparam P2 optional second point type
 
41
            \author Adapted from: http://williams.best.vwh.net/avform.htm
 
42
            \see http://en.wikipedia.org/wiki/Great-circle_distance
 
43
            \note It says: <em>The great circle distance d between two points with coordinates {lat1,lon1} and {lat2,lon2} is given by:
 
44
                        d=acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon1-lon2))
 
45
                    A mathematically equivalent formula, which is less subject to rounding error for short distances is:
 
46
                        d=2*asin(sqrt((sin((lat1-lat2)/2))^2 + cos(lat1)*cos(lat2)*(sin((lon1-lon2)/2))^2))</em>
 
47
        */
 
48
        template <typename P1, typename P2 = P1>
 
49
        class haversine
 
50
        {
 
51
            public :
 
52
                //typedef spherical_distance return_type;
 
53
                typedef double return_type;
 
54
 
 
55
                inline haversine(double r = constants::average_earth_radius)
 
56
                    : m_radius(r)
 
57
                {}
 
58
 
 
59
                inline return_type operator()(const P1& p1, const P2& p2) const
 
60
                {
 
61
                    return calc(get_as_radian<0>(p1), get_as_radian<1>(p1),
 
62
                                    get_as_radian<0>(p2), get_as_radian<1>(p2));
 
63
                }
 
64
 
 
65
            private :
 
66
                double m_radius;
 
67
                typedef typename coordinate_type<P1>::type T1;
 
68
                typedef typename coordinate_type<P2>::type T2;
 
69
 
 
70
                inline return_type calc(const T1& lon1, const T1& lat1, const T2& lon2, const T2& lat2) const
 
71
                {
 
72
                    double a = math::hav(lat2 - lat1) + cos(lat1) * cos(lat2) * math::hav(lon2 - lon1);
 
73
                    double c = 2.0 * asin(sqrt(a));
 
74
                    return return_type(m_radius * c);
 
75
                }
 
76
        };
 
77
 
 
78
 
 
79
 
 
80
        /*!
 
81
            \brief Strategy functor for distance point to segment calculation
 
82
            \ingroup distance
 
83
            \details Class which calculates the distance of a point to a segment, using latlong points
 
84
            \tparam P point type
 
85
            \tparam S segment type
 
86
        */
 
87
        template <typename P, typename S>
 
88
        class ll_point_segment
 
89
        {
 
90
            public :
 
91
                typedef double return_type;
 
92
 
 
93
                inline ll_point_segment(double r = constants::average_earth_radius) : m_radius(r)
 
94
                {}
 
95
 
 
96
                inline return_type operator()(const P& p, const S& s) const
 
97
                {
 
98
                    PR pr, ps1, ps2;
 
99
 
 
100
                    // Select transformation strategy and transform to radians (if necessary)
 
101
                    typename strategy_transform<
 
102
                                typename cs_tag<P>::type,
 
103
                                typename cs_tag<PR>::type,
 
104
                                typename coordinate_system<P>::type,
 
105
                                typename coordinate_system<PR>::type,
 
106
                                dimension<P>::value,
 
107
                                dimension<PR>::value,
 
108
                                P, PR>::type transform_strategy;
 
109
 
 
110
 
 
111
                    // TODO
 
112
                    // ASSUMPTION: segment
 
113
                    // SOLVE THIS USING OTHER FUNCTIONS using get<,>
 
114
                    transform_strategy(p, pr);
 
115
                    transform_strategy(s.first, ps1);
 
116
                    transform_strategy(s.second, ps2);
 
117
                    return calc(pr, ps1, ps2);
 
118
                }
 
119
 
 
120
            private :
 
121
                typedef point_ll
 
122
                    <
 
123
                        typename coordinate_type<P>::type,
 
124
                        cs::geographic<radian>,
 
125
                        ggl::dimension<P>::type::value
 
126
                    > PR;
 
127
                double m_radius;
 
128
 
 
129
                /// Calculate course (bearing) between two points. Might be moved to a "course formula" ...
 
130
                inline double course(const PR& p1, const PR& p2) const
 
131
                {
 
132
                    /***
 
133
                        Course between points
 
134
 
 
135
                        We obtain the initial course, tc1, (at point 1) from point 1 to point 2 by the following. The formula fails if the initial point is a pole. We can special case this with:
 
136
 
 
137
                        IF (cos(lat1) < EPS)   // EPS a small number ~ machine precision
 
138
                          IF (lat1 > 0): tc1= pi        //  starting from N pole
 
139
                          ELSE: tc1= 2*pi         //  starting from S pole
 
140
                          ENDIF
 
141
                        ENDIF
 
142
 
 
143
                        For starting points other than the poles:
 
144
                        IF sin(lon2-lon1)<0: tc1=acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1)))
 
145
                        ELSE: tc1=2*pi-acos((sin(lat2)-sin(lat1)*cos(d))/(sin(d)*cos(lat1)))
 
146
                        ENDIF
 
147
 
 
148
                        An alternative formula, not requiring the pre-computation of d, the distance between the points, is:
 
149
                           tc1=mod(atan2(sin(lon1-lon2)*cos(lat2),
 
150
                                   cos(lat1)*sin(lat2)-sin(lat1)*cos(lat2)*cos(lon1-lon2))
 
151
                                   , 2*pi)
 
152
                     ***/
 
153
                    double dlon = get<0>(p2) - get<0>(p1);
 
154
                    double cos_p2lat = cos(get<1>(p2));
 
155
                    return atan2(sin(dlon) * cos_p2lat,
 
156
                        cos(get<1>(p1)) * sin(get<1>(p2))
 
157
                        - sin(get<1>(p1)) * cos_p2lat * cos(dlon));
 
158
                }
 
159
 
 
160
                inline return_type calc(const PR& p, const PR& sp1, const PR& sp2) const
 
161
                {
 
162
                    /***
 
163
                    Cross track error:
 
164
                    Suppose you are proceeding on a great circle route from A to B (course =crs_AB) and end up at D, perhaps off course.
 
165
                    (We presume that A is ot a pole!) You can calculate the course from A to D (crs_AD) and the distance from A to D (dist_AD)
 
166
                    using the formulae above. In shifteds of these the cross track error, XTD, (distance off course) is given by
 
167
 
 
168
                               XTD =asin(sin(dist_AD)*sin(crs_AD-crs_AB))
 
169
 
 
170
                    (positive XTD means right of course, negative means left)
 
171
                    (If the point A is the N. or S. Pole replace crs_AD-crs_AB with
 
172
                    lon_D-lon_B or lon_B-lon_D, respectively.)
 
173
                     ***/
 
174
 
 
175
                    // Calculate distances, in radians, on the unit sphere
 
176
                    // It seems not useful to let this strategy be templatized, it should be in radians and on the unit sphere
 
177
                    strategy::distance::haversine<PR, PR> strategy(1.0);
 
178
                    double d1 = strategy(sp1, p);
 
179
 
 
180
                    // Actually, calculation of d2 not necessary if we know that the projected point is on the great circle...
 
181
                    double d2 = strategy(sp2, p);
 
182
 
 
183
                    // Source: http://williams.best.vwh.net/avform.htm
 
184
 
 
185
                    double crs_AD = course(sp1, p);
 
186
                    double crs_AB = course(sp1, sp2);
 
187
                    double XTD = fabs(asin(sin(d1) * sin(crs_AD - crs_AB)));
 
188
 
 
189
                    // Return shortest distance, either to projected point on segment sp1-sp2, or to sp1, or to sp2
 
190
                    return return_type(m_radius * (std::min)((std::min)(d1, d2), XTD));
 
191
                }
 
192
        };
 
193
 
 
194
 
 
195
 
 
196
 
 
197
 
 
198
    } // namespace distance
 
199
 
 
200
 
 
201
 
 
202
 
 
203
} // namespace strategy
 
204
 
 
205
 
 
206
#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 
207
template <typename P1, typename P2>
 
208
struct strategy_distance<spherical_tag, spherical_tag, P1, P2>
 
209
{
 
210
    typedef strategy::distance::haversine<P1, P2> type;
 
211
};
 
212
 
 
213
 
 
214
template <typename Point, typename Segment>
 
215
struct strategy_distance_segment<spherical_tag, spherical_tag, Point, Segment>
 
216
{
 
217
    typedef strategy::distance::ll_point_segment<Point, Segment> type;
 
218
};
 
219
 
 
220
 
 
221
// Use this point-segment for geographic as well. For point-point Andoyer is used there, see file "geo_distance.hpp"
 
222
template <typename Point, typename Segment>
 
223
struct strategy_distance_segment<geographic_tag, geographic_tag, Point, Segment>
 
224
{
 
225
    typedef strategy::distance::ll_point_segment<Point, Segment> type;
 
226
};
 
227
 
 
228
 
 
229
 
 
230
 
 
231
template <typename P1, typename P2>
 
232
struct strategy_tag<strategy::distance::haversine<P1, P2> >
 
233
{
 
234
    typedef strategy_tag_distance_point_point type;
 
235
};
 
236
 
 
237
template <typename Point, typename Segment>
 
238
struct strategy_tag<strategy::distance::ll_point_segment<Point, Segment> >
 
239
{
 
240
    typedef strategy_tag_distance_point_segment type;
 
241
};
 
242
 
 
243
 
 
244
 
 
245
#endif
 
246
 
 
247
 
 
248
 
 
249
 
 
250
 
 
251
 
 
252
} // namespace ggl
 
253
 
 
254
 
 
255
#endif // GGL_STRATEGY_SPHERICAL_DISTANCE_HPP