diff options
Diffstat (limited to 'boost/geometry/strategies/spherical/distance_cross_track.hpp')
-rw-r--r-- | boost/geometry/strategies/spherical/distance_cross_track.hpp | 531 |
1 files changed, 485 insertions, 46 deletions
diff --git a/boost/geometry/strategies/spherical/distance_cross_track.hpp b/boost/geometry/strategies/spherical/distance_cross_track.hpp index a40f03dbaf..486c7effbd 100644 --- a/boost/geometry/strategies/spherical/distance_cross_track.hpp +++ b/boost/geometry/strategies/spherical/distance_cross_track.hpp @@ -1,6 +1,11 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2014. +// Modifications copyright (c) 2014, Oracle and/or its affiliates. + +// Contributed and/or modified by Menelaos Karavelas, 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 @@ -9,15 +14,17 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP +#include <algorithm> #include <boost/config.hpp> #include <boost/concept_check.hpp> #include <boost/mpl/if.hpp> -#include <boost/type_traits.hpp> +#include <boost/type_traits/is_void.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/core/access.hpp> #include <boost/geometry/core/radian_access.hpp> +#include <boost/geometry/core/tags.hpp> #include <boost/geometry/algorithms/detail/course.hpp> @@ -25,39 +32,297 @@ #include <boost/geometry/strategies/concepts/distance_concept.hpp> #include <boost/geometry/strategies/spherical/distance_haversine.hpp> -#include <boost/geometry/util/promote_floating_point.hpp> #include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/promote_floating_point.hpp> +#include <boost/geometry/util/select_calculation_type.hpp> #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK # include <boost/geometry/io/dsv/write.hpp> #endif - namespace boost { namespace geometry { namespace strategy { namespace distance { -/*! -\brief Strategy functor for distance point to segment calculation -\ingroup strategies -\details Class which calculates the distance of a point to a segment, for points on a sphere or globe -\see http://williams.best.vwh.net/avform.htm -\tparam CalculationType \tparam_calculation -\tparam Strategy underlying point-point distance strategy, defaults to haversine -\qbk{ -[heading See also] -[link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)] -} +namespace comparable +{ +/* + Given a spherical segment AB and a point D, we are interested in + computing the distance of D from AB. This is usually known as the + cross track distance. + + If the projection (along great circles) of the point D lies inside + the segment AB, then the distance (cross track error) XTD is given + by the formula (see http://williams.best.vwh.net/avform.htm#XTE): + + XTD = asin( sin(dist_AD) * sin(crs_AD-crs_AB) ) + + where dist_AD is the great circle distance between the points A and + B, and crs_AD, crs_AB is the course (bearing) between the points A, + D and A, B, respectively. + + If the point D does not project inside the arc AB, then the distance + of D from AB is the minimum of the two distances dist_AD and dist_BD. + + Our reference implementation for this procedure is listed below + (this was the old Boost.Geometry implementation of the cross track distance), + where: + * The member variable m_strategy is the underlying haversine strategy. + * p stands for the point D. + * sp1 stands for the segment endpoint A. + * sp2 stands for the segment endpoint B. + + ================= reference implementation -- start ================= + + return_type d1 = m_strategy.apply(sp1, p); + return_type d3 = m_strategy.apply(sp1, sp2); + + if (geometry::math::equals(d3, 0.0)) + { + // "Degenerate" segment, return either d1 or d2 + return d1; + } + + return_type d2 = m_strategy.apply(sp2, p); + + return_type crs_AD = geometry::detail::course<return_type>(sp1, p); + return_type crs_AB = geometry::detail::course<return_type>(sp1, sp2); + return_type crs_BA = crs_AB - geometry::math::pi<return_type>(); + return_type crs_BD = geometry::detail::course<return_type>(sp2, p); + return_type d_crs1 = crs_AD - crs_AB; + return_type d_crs2 = crs_BD - crs_BA; + + // d1, d2, d3 are in principle not needed, only the sign matters + return_type projection1 = cos( d_crs1 ) * d1 / d3; + return_type projection2 = cos( d_crs2 ) * d2 / d3; + + if (projection1 > 0.0 && projection2 > 0.0) + { + return_type XTD + = radius() * math::abs( asin( sin( d1 / radius() ) * sin( d_crs1 ) )); + + // Return shortest distance, projected point on segment sp1-sp2 + return return_type(XTD); + } + else + { + // Return shortest distance, project either on point sp1 or sp2 + return return_type( (std::min)( d1 , d2 ) ); + } + + ================= reference implementation -- end ================= + + + Motivation + ---------- + In what follows we develop a comparable version of the cross track + distance strategy, that meets the following goals: + * It is more efficient than the original cross track strategy (less + operations and less calls to mathematical functions). + * Distances using the comparable cross track strategy can not only + be compared with other distances using the same strategy, but also with + distances computed with the comparable version of the haversine strategy. + * It can serve as the basis for the computation of the cross track distance, + as it is more efficient to compute its comparable version and + transform that to the actual cross track distance, rather than + follow/use the reference implementation listed above. + + Major idea + ---------- + The idea here is to use the comparable haversine strategy to compute + the distances d1, d2 and d3 in the above listing. Once we have done + that we need also to make sure that instead of returning XTD (as + computed above) that we return a distance CXTD that is compatible + with the comparable haversine distance. To achieve this CXTD must satisfy + the relation: + XTD = 2 * R * asin( sqrt(XTD) ) + where R is the sphere's radius. + + Below we perform the mathematical analysis that show how to compute CXTD. + + + Mathematical analysis + --------------------- + Below we use the following trigonometric identities: + sin(2 * x) = 2 * sin(x) * cos(x) + cos(asin(x)) = sqrt(1 - x^2) + + Observation: + The distance d1 needed when the projection of the point D is within the + segment must be the true distance. However, comparable::haversine<> + returns a comparable distance instead of the one needed. + To remedy this, we implicitly compute what is needed. + More precisely, we need to compute sin(true_d1): + + sin(true_d1) = sin(2 * asin(sqrt(d1))) + = 2 * sin(asin(sqrt(d1)) * cos(asin(sqrt(d1))) + = 2 * sqrt(d1) * sqrt(1-(sqrt(d1))^2) + = 2 * sqrt(d1 - d1 * d1) + This relation is used below. + + As we mentioned above the goal is to find CXTD (named "a" below for + brevity) such that ("b" below stands for "d1", and "c" for "d_crs1"): + + 2 * R * asin(sqrt(a)) == R * asin(2 * sqrt(b-b^2) * sin(c)) + + Analysis: + 2 * R * asin(sqrt(a)) == R * asin(2 * sqrt(b-b^2) * sin(c)) + <=> 2 * asin(sqrt(a)) == asin(sqrt(b-b^2) * sin(c)) + <=> sin(2 * asin(sqrt(a))) == 2 * sqrt(b-b^2) * sin(c) + <=> 2 * sin(asin(sqrt(a))) * cos(asin(sqrt(a))) == 2 * sqrt(b-b^2) * sin(c) + <=> 2 * sqrt(a) * sqrt(1-a) == 2 * sqrt(b-b^2) * sin(c) + <=> sqrt(a) * sqrt(1-a) == sqrt(b-b^2) * sin(c) + <=> sqrt(a-a^2) == sqrt(b-b^2) * sin(c) + <=> a-a^2 == (b-b^2) * (sin(c))^2 + + Consider the quadratic equation: x^2-x+p^2 == 0, + where p = sqrt(b-b^2) * sin(c); its discriminant is: + d = 1 - 4 * p^2 = 1 - 4 * (b-b^2) * (sin(c))^2 + + The two solutions are: + a_1 = (1 - sqrt(d)) / 2 + a_2 = (1 + sqrt(d)) / 2 + + Which one to choose? + "a" refers to the distance (on the unit sphere) of D from the + supporting great circle Circ(A,B) of the segment AB. + The two different values for "a" correspond to the lengths of the two + arcs delimited D and the points of intersection of Circ(A,B) and the + great circle perperdicular to Circ(A,B) passing through D. + Clearly, the value we want is the smallest among these two distances, + hence the root we must choose is the smallest root among the two. + + So the answer is: + CXTD = ( 1 - sqrt(1 - 4 * (b-b^2) * (sin(c))^2) ) / 2 + + Therefore, in order to implement the comparable version of the cross + track strategy we need to: + (1) Use the comparable version of the haversine strategy instead of + the non-comparable one. + (2) Instead of return XTD when D projects inside the segment AB, we + need to return CXTD, given by the following formula: + CXTD = ( 1 - sqrt(1 - 4 * (d1-d1^2) * (sin(d_crs1))^2) ) / 2; + + + Complexity Analysis + ------------------- + In the analysis that follows we refer to the actual implementation below. + In particular, instead of computing CXTD as above, we use the more + efficient (operation-wise) computation of CXTD shown here: + + return_type sin_d_crs1 = sin(d_crs1); + return_type d1_x_sin = d1 * sin_d_crs1; + return_type d = d1_x_sin * (sin_d_crs1 - d1_x_sin); + return d / (0.5 + math::sqrt(0.25 - d)); + + Notice that instead of computing: + 0.5 - 0.5 * sqrt(1 - 4 * d) = 0.5 - sqrt(0.25 - d) + we use the following formula instead: + d / (0.5 + sqrt(0.25 - d)). + This is done for numerical robustness. The expression 0.5 - sqrt(0.25 - x) + has large numerical errors for values of x close to 0 (if using doubles + the error start to become large even when d is as large as 0.001). + To remedy that, we re-write 0.5 - sqrt(0.25 - x) as: + 0.5 - sqrt(0.25 - d) + = (0.5 - sqrt(0.25 - d) * (0.5 - sqrt(0.25 - d)) / (0.5 + sqrt(0.25 - d)). + The numerator is the difference of two squares: + (0.5 - sqrt(0.25 - d) * (0.5 - sqrt(0.25 - d)) + = 0.5^2 - (sqrt(0.25 - d))^ = 0.25 - (0.25 - d) = d, + which gives the expression we use. + + For the complexity analysis, we distinguish between two cases: + (A) The distance is realized between the point D and an + endpoint of the segment AB + + Gains: + Since we are using comparable::haversine<> which is called + 3 times, we gain: + -> 3 calls to sqrt + -> 3 calls to asin + -> 6 multiplications + + Loses: None + + So the net gain is: + -> 6 function calls (sqrt/asin) + -> 6 arithmetic operations + + If we use comparable::cross_track<> to compute + cross_track<> we need to account for a call to sqrt, a call + to asin and 2 multiplications. In this case the net gain is: + -> 4 function calls (sqrt/asin) + -> 4 arithmetic operations + + + (B) The distance is realized between the point D and an + interior point of the segment AB + + Gains: + Since we are using comparable::haversine<> which is called + 3 times, we gain: + -> 3 calls to sqrt + -> 3 calls to asin + -> 6 multiplications + Also we gain the operations used to compute XTD: + -> 2 calls to sin + -> 1 call to asin + -> 1 call to abs + -> 2 multiplications + -> 1 division + So the total gains are: + -> 9 calls to sqrt/sin/asin + -> 1 call to abs + -> 8 multiplications + -> 1 division + + Loses: + To compute a distance compatible with comparable::haversine<> + we need to perform a few more operations, namely: + -> 1 call to sin + -> 1 call to sqrt + -> 2 multiplications + -> 1 division + -> 1 addition + -> 2 subtractions + + So roughly speaking the net gain is: + -> 8 fewer function calls and 3 fewer arithmetic operations + + If we were to implement cross_track directly from the + comparable version (much like what haversine<> does using + comparable::haversine<>) we need additionally + -> 2 function calls (asin/sqrt) + -> 2 multiplications + + So it pays off to re-implement cross_track<> to use + comparable::cross_track<>; in this case the net gain would be: + -> 6 function calls + -> 1 arithmetic operation + + Summary/Conclusion + ------------------ + Following the mathematical and complexity analysis above, the + comparable cross track strategy (as implemented below) satisfies + all the goal mentioned in the beginning: + * It is more efficient than its non-comparable counter-part. + * Comparable distances using this new strategy can also be compared + with comparable distances computed with the comparable haversine + strategy. + * It turns out to be more efficient to compute the actual cross + track distance XTD by first computing CXTD, and then computing + XTD by means of the formula: + XTD = 2 * R * asin( sqrt(CXTD) ) */ + template < typename CalculationType = void, - typename Strategy = haversine<double, CalculationType> + typename Strategy = comparable::haversine<double, CalculationType> > class cross_track { @@ -106,6 +371,14 @@ public : typedef typename return_type<Point, PointOfSegment>::type return_type; +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK + std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " " << crs_AD * geometry::math::r2d << std::endl; + std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " " << crs_AB * geometry::math::r2d << std::endl; + std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " " << crs_BD * geometry::math::r2d << std::endl; + std::cout << "Projection AD-AB " << projection1 << " : " << d_crs1 * geometry::math::r2d << std::endl; + std::cout << "Projection BD-BA " << projection2 << " : " << d_crs2 * geometry::math::r2d << std::endl; +#endif + // http://williams.best.vwh.net/avform.htm#XTE return_type d1 = m_strategy.apply(sp1, p); return_type d3 = m_strategy.apply(sp1, sp2); @@ -129,25 +402,34 @@ public : return_type projection1 = cos( d_crs1 ) * d1 / d3; return_type projection2 = cos( d_crs2 ) * d2 / d3; -#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK - std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " " << crs_AD * geometry::math::r2d << std::endl; - std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " " << crs_AB * geometry::math::r2d << std::endl; - std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " " << crs_BD * geometry::math::r2d << std::endl; - std::cout << "Projection AD-AB " << projection1 << " : " << d_crs1 * geometry::math::r2d << std::endl; - std::cout << "Projection BD-BA " << projection2 << " : " << d_crs2 * geometry::math::r2d << std::endl; -#endif - - if(projection1 > 0.0 && projection2 > 0.0) + if (projection1 > 0.0 && projection2 > 0.0) { - return_type XTD = radius() * geometry::math::abs( asin( sin( d1 / radius() ) * sin( d_crs1 ) )); +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK + return_type XTD = radius() * geometry::math::abs( asin( sin( d1 ) * sin( d_crs1 ) )); - #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK std::cout << "Projection ON the segment" << std::endl; - std::cout << "XTD: " << XTD << " d1: " << d1 << " d2: " << d2 << std::endl; + std::cout << "XTD: " << XTD + << " d1: " << (d1 * radius()) + << " d2: " << (d2 * radius()) + << std::endl; #endif - - // Return shortest distance, projected point on segment sp1-sp2 - return return_type(XTD); + return_type const half(0.5); + return_type const quarter(0.25); + + return_type sin_d_crs1 = sin(d_crs1); + /* + This is the straightforward obvious way to continue: + + return_type discriminant + = 1.0 - 4.0 * (d1 - d1 * d1) * sin_d_crs1 * sin_d_crs1; + return 0.5 - 0.5 * math::sqrt(discriminant); + + Below we optimize the number of arithmetic operations + and account for numerical robustness: + */ + return_type d1_x_sin = d1 * sin_d_crs1; + return_type d = d1_x_sin * (sin_d_crs1 - d1_x_sin); + return d / (half + math::sqrt(quarter - d)); } else { @@ -164,6 +446,95 @@ public : { return m_strategy.radius(); } private : + Strategy m_strategy; +}; + +} // namespace comparable + + +/*! +\brief Strategy functor for distance point to segment calculation +\ingroup strategies +\details Class which calculates the distance of a point to a segment, for points on a sphere or globe +\see http://williams.best.vwh.net/avform.htm +\tparam CalculationType \tparam_calculation +\tparam Strategy underlying point-point distance strategy, defaults to haversine + +\qbk{ +[heading See also] +[link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)] +} + +*/ +template +< + typename CalculationType = void, + typename Strategy = haversine<double, CalculationType> +> +class cross_track +{ +public : + template <typename Point, typename PointOfSegment> + struct return_type + : promote_floating_point + < + typename select_calculation_type + < + Point, + PointOfSegment, + CalculationType + >::type + > + {}; + + inline cross_track() + {} + + explicit inline cross_track(typename Strategy::radius_type const& r) + : m_strategy(r) + {} + + inline cross_track(Strategy const& s) + : m_strategy(s) + {} + + + // It might be useful in the future + // to overload constructor with strategy info. + // crosstrack(...) {} + + + template <typename Point, typename PointOfSegment> + inline typename return_type<Point, PointOfSegment>::type + apply(Point const& p, PointOfSegment const& sp1, PointOfSegment const& sp2) const + { + +#if !defined(BOOST_MSVC) + BOOST_CONCEPT_ASSERT + ( + (concept::PointDistanceStrategy<Strategy, Point, PointOfSegment>) + ); +#endif + typedef typename return_type<Point, PointOfSegment>::type return_type; + typedef cross_track<CalculationType, Strategy> this_type; + + typedef typename services::comparable_type + < + this_type + >::type comparable_type; + + comparable_type cstrategy + = services::get_comparable<this_type>::apply(m_strategy); + + return_type const a = cstrategy.apply(p, sp1, sp2); + return_type const c = return_type(2.0) * asin(math::sqrt(a)); + return c * radius(); + } + + inline typename Strategy::radius_type radius() const + { return m_strategy.radius(); } + +private : Strategy m_strategy; }; @@ -190,8 +561,10 @@ struct return_type<cross_track<CalculationType, Strategy>, P, PS> template <typename CalculationType, typename Strategy> struct comparable_type<cross_track<CalculationType, Strategy> > { - // There is no shortcut, so the strategy itself is its comparable type - typedef cross_track<CalculationType, Strategy> type; + typedef comparable::cross_track + < + CalculationType, typename comparable_type<Strategy>::type + > type; }; @@ -207,7 +580,8 @@ struct get_comparable<cross_track<CalculationType, Strategy> > cross_track<CalculationType, Strategy> >::type comparable_type; public : - static inline comparable_type apply(cross_track<CalculationType, Strategy> const& strategy) + static inline comparable_type + apply(cross_track<CalculationType, Strategy> const& strategy) { return comparable_type(strategy.radius()); } @@ -218,21 +592,95 @@ template < typename CalculationType, typename Strategy, - typename P, typename PS + typename P, + typename PS > struct result_from_distance<cross_track<CalculationType, Strategy>, P, PS> { private : - typedef typename cross_track<CalculationType, Strategy>::template return_type<P, PS> return_type; + typedef typename cross_track + < + CalculationType, Strategy + >::template return_type<P, PS>::type return_type; public : template <typename T> - static inline return_type apply(cross_track<CalculationType, Strategy> const& , T const& distance) + static inline return_type + apply(cross_track<CalculationType, Strategy> const& , T const& distance) { return distance; } }; +// Specializations for comparable::cross_track +template <typename RadiusType, typename CalculationType> +struct tag<comparable::cross_track<RadiusType, CalculationType> > +{ + typedef strategy_tag_distance_point_segment type; +}; + + +template +< + typename RadiusType, + typename CalculationType, + typename P, + typename PS +> +struct return_type<comparable::cross_track<RadiusType, CalculationType>, P, PS> + : comparable::cross_track + < + RadiusType, CalculationType + >::template return_type<P, PS> +{}; + + +template <typename RadiusType, typename CalculationType> +struct comparable_type<comparable::cross_track<RadiusType, CalculationType> > +{ + typedef comparable::cross_track<RadiusType, CalculationType> type; +}; + + +template <typename RadiusType, typename CalculationType> +struct get_comparable<comparable::cross_track<RadiusType, CalculationType> > +{ +private : + typedef comparable::cross_track<RadiusType, CalculationType> this_type; +public : + static inline this_type apply(this_type const& input) + { + return input; + } +}; + + +template +< + typename RadiusType, + typename CalculationType, + typename P, + typename PS +> +struct result_from_distance + < + comparable::cross_track<RadiusType, CalculationType>, P, PS + > +{ +private : + typedef comparable::cross_track<RadiusType, CalculationType> strategy_type; + typedef typename return_type<strategy_type, P, PS>::type return_type; +public : + template <typename T> + static inline return_type apply(strategy_type const& strategy, + T const& distance) + { + return_type const s + = sin( (distance / strategy.radius()) / return_type(2.0) ); + return s * s; + } +}; + /* @@ -309,17 +757,8 @@ struct default_strategy } // namespace services #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - }} // namespace strategy::distance - -#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - - -#endif - - }} // namespace boost::geometry - #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_HPP |