Subject: [Boost-commit] svn:boost r85063 - in trunk: boost/geometry/extensions/gis/geographic/strategies libs/geometry/extensions/test/gis/latlong
From: barend.gehrels_at_[hidden]
Date: 2013-07-17 17:08:39
Author: barendgehrels
Date: 2013-07-17 17:08:38 EDT (Wed, 17 Jul 2013)
New Revision: 85063
URL: http://svn.boost.org/trac/boost/changeset/85063
Log:
[geometry][extension] adapted vincenty to new strategy structure
Text files modified:
trunk/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp | 202 ++++++++++++++++++---------------------
trunk/libs/geometry/extensions/test/gis/latlong/vincenty.cpp | 17 ++
2 files changed, 108 insertions(+), 111 deletions(-)
Modified: trunk/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp
==============================================================================
--- trunk/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp Wed Jul 17 16:58:28 2013 (r85062)
+++ trunk/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp 2013-07-17 17:08:38 EDT (Wed, 17 Jul 2013) (r85063)
@@ -33,8 +33,7 @@
/*!
\brief Distance calculation formulae on latlong coordinates, after Vincenty, 1975
\ingroup distance
-\tparam Point1 \tparam_first_point
-\tparam Point2 \tparam_second_point
+\tparam RadiusType type of radius (of the Earth)
\tparam CalculationType \tparam_calculation
\author See http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
\author Adapted from various implementations to get it close to the original document
@@ -45,114 +44,118 @@
*/
template
<
- typename Point1,
- typename Point2 = Point1,
+ typename RadiusType,
typename CalculationType = void
>
class vincenty
{
public :
- typedef typename promote_floating_point
- <
- typename select_most_precise
- <
- typename select_calculation_type
- <
- Point1,
- Point2,
- CalculationType
- >::type,
- double // to avoid bad results in float
- >::type
- >::type calculation_type;
+ template <typename Point1, typename Point2>
+ struct calculation_type
+ : promote_floating_point
+ <
+ typename select_calculation_type
+ <
+ Point1,
+ Point2,
+ CalculationType
+ >::type
+ >
+ {};
inline vincenty()
{}
- explicit inline vincenty(geometry::detail::ellipsoid<calculation_type> const& e)
+ explicit inline vincenty(geometry::detail::ellipsoid<RadiusType> const& e)
: m_ellipsoid(e)
{}
- inline calculation_type apply(Point1 const& p1, Point2 const& p2) const
- {
- return calculate(get_as_radian<0>(p1), get_as_radian<1>(p1),
- get_as_radian<0>(p2), get_as_radian<1>(p2));
+ template <typename Point1, typename Point2>
+ inline typename calculation_type<Point1, Point2>::type
+ apply(Point1 const& point1, Point2 const& point2) const
+ {
+ return calculate<typename calculation_type<Point1, Point2>::type>
+ (
+ get_as_radian<0>(point1), get_as_radian<1>(point1),
+ get_as_radian<0>(point2), get_as_radian<1>(point2)
+ );
}
- inline geometry::detail::ellipsoid<calculation_type> ellipsoid() const
+ inline geometry::detail::ellipsoid<RadiusType> ellipsoid() const
{
return m_ellipsoid;
}
private :
- geometry::detail::ellipsoid<calculation_type> m_ellipsoid;
+ geometry::detail::ellipsoid<RadiusType> m_ellipsoid;
- inline calculation_type calculate(calculation_type const& lon1,
- calculation_type const& lat1,
- calculation_type const& lon2,
- calculation_type const& lat2) const
- {
- calculation_type const c2 = 2;
- calculation_type const pi = geometry::math::pi<calculation_type>();
- calculation_type const two_pi = c2 * pi;
+ template <typename CT, typename T>
+ inline CT calculate(T const& lon1,
+ T const& lat1,
+ T const& lon2,
+ T const& lat2) const
+ {
+ CT const c2 = 2;
+ CT const pi = geometry::math::pi<CT>();
+ CT const two_pi = c2 * pi;
// lambda: difference in longitude on an auxiliary sphere
- calculation_type L = lon2 - lon1;
- calculation_type lambda = L;
+ CT L = lon2 - lon1;
+ CT lambda = L;
if (L < -pi) L += two_pi;
if (L > pi) L -= two_pi;
if (math::equals(lat1, lat2) && math::equals(lon1, lon2))
{
- return calculation_type(0);
+ return CT(0);
}
// U: reduced latitude, defined by tan U = (1-f) tan phi
- calculation_type const c1 = 1;
- calculation_type const one_min_f = c1 - m_ellipsoid.f();
+ CT const c1 = 1;
+ CT const one_min_f = c1 - m_ellipsoid.f();
- calculation_type const U1 = atan(one_min_f * tan(lat1)); // above (1)
- calculation_type const U2 = atan(one_min_f * tan(lat2)); // above (1)
+ CT const U1 = atan(one_min_f * tan(lat1)); // above (1)
+ CT const U2 = atan(one_min_f * tan(lat2)); // above (1)
- calculation_type const cos_U1 = cos(U1);
- calculation_type const cos_U2 = cos(U2);
- calculation_type const sin_U1 = sin(U1);
- calculation_type const sin_U2 = sin(U2);
+ CT const cos_U1 = cos(U1);
+ CT const cos_U2 = cos(U2);
+ CT const sin_U1 = sin(U1);
+ CT const sin_U2 = sin(U2);
// alpha: azimuth of the geodesic at the equator
- calculation_type cos2_alpha;
- calculation_type sin_alpha;
+ CT cos2_alpha;
+ CT sin_alpha;
// sigma: angular distance p1,p2 on the sphere
// sigma1: angular distance on the sphere from the equator to p1
// sigma_m: angular distance on the sphere from the equator to the midpoint of the line
- calculation_type sigma;
- calculation_type sin_sigma;
- calculation_type cos2_sigma_m;
-
- calculation_type previous_lambda;
-
- calculation_type const c3 = 3;
- calculation_type const c4 = 4;
- calculation_type const c6 = 6;
- calculation_type const c16 = 16;
+ CT sigma;
+ CT sin_sigma;
+ CT cos2_sigma_m;
+
+ CT previous_lambda;
+
+ CT const c3 = 3;
+ CT const c4 = 4;
+ CT const c6 = 6;
+ CT const c16 = 16;
- calculation_type const c_e_12 = 1e-12;
+ CT const c_e_12 = 1e-12;
do
{
previous_lambda = lambda; // (13)
- calculation_type sin_lambda = sin(lambda);
- calculation_type cos_lambda = cos(lambda);
+ CT sin_lambda = sin(lambda);
+ CT cos_lambda = cos(lambda);
sin_sigma = sqrt(math::sqr(cos_U2 * sin_lambda) + math::sqr(cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda)); // (14)
- calculation_type cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
+ CT cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17)
cos2_alpha = c1 - math::sqr(sin_alpha);
cos2_sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
- calculation_type C = m_ellipsoid.f()/c16 * cos2_alpha * (c4 + m_ellipsoid.f() * (c4 - c3 * cos2_alpha)); // (10)
+ CT C = m_ellipsoid.f()/c16 * cos2_alpha * (c4 + m_ellipsoid.f() * (c4 - c3 * cos2_alpha)); // (10)
sigma = atan2(sin_sigma, cos_sigma); // (16)
lambda = L + (c1 - C) * m_ellipsoid.f() * sin_alpha *
(sigma + C * sin_sigma * ( cos2_sigma_m + C * cos_sigma * (-c1 + c2 * math::sqr(cos2_sigma_m)))); // (11)
@@ -160,24 +163,24 @@
} while (geometry::math::abs(previous_lambda - lambda) > c_e_12
&& geometry::math::abs(lambda) < pi);
- calculation_type sqr_u = cos2_alpha * (math::sqr(m_ellipsoid.a()) - math::sqr(m_ellipsoid.b())) / math::sqr(m_ellipsoid.b()); // above (1)
+ CT sqr_u = cos2_alpha * (math::sqr(m_ellipsoid.a()) - math::sqr(m_ellipsoid.b())) / math::sqr(m_ellipsoid.b()); // above (1)
// Oops getting hard here
// (again, problem is that ttmath cannot divide by doubles, which is OK)
- calculation_type const c47 = 47;
- calculation_type const c74 = 74;
- calculation_type const c128 = 128;
- calculation_type const c256 = 256;
- calculation_type const c175 = 175;
- calculation_type const c320 = 320;
- calculation_type const c768 = 768;
- calculation_type const c1024 = 1024;
- calculation_type const c4096 = 4096;
- calculation_type const c16384 = 16384;
-
- calculation_type A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3)
- calculation_type B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4)
- calculation_type delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m)
+ CT const c47 = 47;
+ CT const c74 = 74;
+ CT const c128 = 128;
+ CT const c256 = 256;
+ CT const c175 = 175;
+ CT const c320 = 320;
+ CT const c768 = 768;
+ CT const c1024 = 1024;
+ CT const c4096 = 4096;
+ CT const c16384 = 16384;
+
+ CT A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3)
+ CT B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4)
+ CT delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m)
- (B/c6) * cos2_sigma_m * (-c3 + c4 * math::sqr(sin_sigma)) * (-c3 + c4 * cos2_sigma_m))); // (6)
return m_ellipsoid.b() * A * (sigma - delta_sigma); // (19)
@@ -188,57 +191,41 @@
namespace services
{
-template <typename Point1, typename Point2>
-struct tag<strategy::distance::vincenty<Point1, Point2> >
+template <typename RadiusType, typename CalculationType>
+struct tag<vincenty<RadiusType, CalculationType> >
{
typedef strategy_tag_distance_point_point type;
};
-template <typename Point1, typename Point2>
-struct return_type<strategy::distance::vincenty<Point1, Point2> >
-{
- typedef typename strategy::distance::vincenty<Point1, Point2>::calculation_type type;
-};
-
-
-template <typename Point1, typename Point2, typename P1, typename P2>
-struct similar_type<vincenty<Point1, Point2>, P1, P2>
-{
- typedef vincenty<P1, P2> type;
-};
-
+template <typename RadiusType, typename CalculationType, typename P1, typename P2>
+struct return_type<vincenty<RadiusType, CalculationType>, P1, P2>
+ : vincenty<RadiusType, CalculationType>::template calculation_type<P1, P2>
+{};
-template <typename Point1, typename Point2, typename P1, typename P2>
-struct get_similar<vincenty<Point1, Point2>, P1, P2>
-{
- static inline vincenty<P1, P2> apply(vincenty<Point1, Point2> const& input)
- {
- return vincenty<P1, P2>(input.ellipsoid());
- }
-};
-template <typename Point1, typename Point2>
-struct comparable_type<vincenty<Point1, Point2> >
+template <typename RadiusType, typename CalculationType>
+struct comparable_type<vincenty<RadiusType, CalculationType> >
{
- typedef vincenty<Point1, Point2> type;
+ typedef vincenty<RadiusType, CalculationType> type;
};
-template <typename Point1, typename Point2>
-struct get_comparable<vincenty<Point1, Point2> >
+template <typename RadiusType, typename CalculationType>
+struct get_comparable<vincenty<RadiusType, CalculationType> >
{
- static inline vincenty<Point1, Point2> apply(vincenty<Point1, Point2> const& input)
+ static inline vincenty<RadiusType, CalculationType> apply(vincenty<RadiusType, CalculationType> const& input)
{
return input;
}
};
-template <typename Point1, typename Point2>
-struct result_from_distance<vincenty<Point1, Point2> >
+template <typename RadiusType, typename CalculationType, typename P1, typename P2>
+struct result_from_distance<vincenty<RadiusType, CalculationType>, P1, P2 >
{
template <typename T>
- static inline typename return_type<vincenty<Point1, Point2> >::type apply(vincenty<Point1, Point2> const& , T const& value)
+ static inline typename return_type<vincenty<RadiusType, CalculationType>, P1, P2>::type
+ apply(vincenty<RadiusType, CalculationType> const& , T const& value)
{
return value;
}
@@ -260,3 +247,4 @@
#endif // BOOST_GEOMETRY_EXTENSIONS_GIS_GEOGRAPHIC_STRATEGIES_VINCENTY_HPP
+
Modified: trunk/libs/geometry/extensions/test/gis/latlong/vincenty.cpp
==============================================================================
--- trunk/libs/geometry/extensions/test/gis/latlong/vincenty.cpp Wed Jul 17 16:58:28 2013 (r85062)
+++ trunk/libs/geometry/extensions/test/gis/latlong/vincenty.cpp 2013-07-17 17:08:38 EDT (Wed, 17 Jul 2013) (r85063)
@@ -33,12 +33,21 @@
template <typename P1, typename P2>
void test_vincenty(double lon1, double lat1, double lon2, double lat2, double expected_km)
{
- typedef bg::strategy::distance::vincenty<P1, P2> vincenty_type;
-
- BOOST_CONCEPT_ASSERT( (bg::concept::PointDistanceStrategy<vincenty_type>) );
+ // Set radius type, but for integer coordinates we want to have floating point radius type
+ typedef typename bg::promote_floating_point
+ <
+ typename bg::coordinate_type<P1>::type
+ >::type rtype;
+
+ typedef bg::strategy::distance::vincenty<rtype> vincenty_type;
+
+ BOOST_CONCEPT_ASSERT(
+ (
+ bg::concept::PointDistanceStrategy<vincenty_type, P1, P2>)
+ );
vincenty_type vincenty;
- typedef typename bg::strategy::distance::services::return_type<vincenty_type>::type return_type;
+ typedef typename bg::strategy::distance::services::return_type<vincenty_type, P1, P2>::type return_type;
P1 p1, p2;