From 43d7f824917a2bc5f720b67162eb202ac2fa5e29 Mon Sep 17 00:00:00 2001 From: Menelaos Karavelas Date: Wed, 19 Nov 2014 23:52:59 +0200 Subject: [PATCH] [strategies][distance][spherical][cross track] update comments' section and make it a stand-alone explanatory document; polish code; remove old code; fix header includes; fix long lines; --- .../spherical/distance_cross_track.hpp | 323 ++++++++++++------ 1 file changed, 212 insertions(+), 111 deletions(-) diff --git a/include/boost/geometry/strategies/spherical/distance_cross_track.hpp b/include/boost/geometry/strategies/spherical/distance_cross_track.hpp index 7a8e435a2..11dd7d9a9 100644 --- a/include/boost/geometry/strategies/spherical/distance_cross_track.hpp +++ b/include/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 #include #include #include -#include +#include #include #include #include +#include #include @@ -25,15 +32,15 @@ #include #include -#include #include +#include +#include #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK # include #endif - namespace boost { namespace geometry { @@ -45,14 +52,108 @@ 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(sp1, p); + return_type crs_AB = geometry::detail::course(sp1, sp2); + return_type crs_BA = crs_AB - geometry::math::pi(); + return_type crs_BD = geometry::detail::course(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 achive 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 --------------------- - Identities used: - sin(2*x) = 2*sin(x)*cos(x) - cos(asin(x)) = sqrt(1-x^2) + 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 need when the projection of the point is is within the + 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. @@ -64,8 +165,10 @@ namespace comparable = 2 * sqrt(d1 - d1 * d1) This relation is used below. - Goal: find an a such that: - 2 * R * asin(sqrt(a)) == R * asin(2 * sqrt(b-b^2) * sin(c)) + 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)) @@ -77,28 +180,48 @@ namespace comparable <=> 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 - Its discriminant is: d = 1 - 4 * p^2 = 1 - 4 * (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 - So solutions are: - x_{1,2} = (1 + sqrt(d))/2, (1 - sqrt(d))/2 + The two solutions are: + a_1 = (1 - sqrt(d)) / 2 + a_2 = (1 + sqrt(d)) / 2 Which one to choose? - My understanding is that "a" is the distance of the third - point to its closest point on the segment. The two - different values for "a" correspond to the lengths of the two - arcs delimited by these to points. - So, my guess is the value we want is the smallest root (unit - tests confirm this choice). + "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: - a = ( 1 - sqrt(1 - 4 * (b-b^2) * (sin(c))^2) ) / 2; + 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 - ------------------ - (A) The distance is realized between the point and an - endpoint of the segment + ------------------- + 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 0.5 - math::sqrt(0.25 + d1_x_sin * (d1_x_sin - sin_d_crs1)); + + 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 @@ -120,8 +243,8 @@ namespace comparable -> 4 arithmetic operations - (B) The distance is realized between the point and an - interior point of the segment + (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 @@ -159,11 +282,24 @@ namespace comparable -> 2 function calls (asin/sqrt) -> 2 multiplications - So it seems like it pays off to re-implement cross_track<> to - use comparable::cross_track<>; in this case the net gain - would be: + 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 -> 2 arithmetic operations + + Summary/Conlusion + ----------------- + 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 @@ -218,6 +354,14 @@ public : typedef typename return_type::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); @@ -244,7 +388,13 @@ public : if (projection1 > 0.0 && projection2 > 0.0) { #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK + return_type XTD = radius() * geometry::math::abs( asin( sin( d1 ) * sin( d_crs1 ) )); + std::cout << "Projection ON the segment" << std::endl; + std::cout << "XTD: " << XTD + << " d1: " << (d1 * radius()) + << " d2: " << (d2 * radius()) + << std::endl; #endif return_type const half(0.5); return_type const quarter(0.25); @@ -283,6 +433,7 @@ private : } // namespace comparable + /*! \brief Strategy functor for distance point to segment calculation \ingroup strategies @@ -346,64 +497,7 @@ public : (concept::PointDistanceStrategy) ); #endif - typedef typename return_type::type return_type; - -#ifdef BOOST_GEOMETRY_USE_OLD_CROSS_TRACK_STRATEGY - // http://williams.best.vwh.net/avform.htm#XTE - 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(sp1, p); - return_type crs_AB = geometry::detail::course(sp1, sp2); - return_type crs_BA = crs_AB - geometry::math::pi(); - return_type crs_BD = geometry::detail::course(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; - -#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) - { - return_type XTD = radius() * geometry::math::abs( asin( sin( d1 / radius() ) * 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; -#endif - - // Return shortest distance, projected point on segment sp1-sp2 - return return_type(XTD); - } - else - { -#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK - std::cout << "Projection OUTSIDE the segment" << std::endl; -#endif - - // Return shortest distance, project either on point sp1 or sp2 - return return_type( (std::min)( d1 , d2 ) ); - } -#else // BOOST_GEOMETRY_USE_OLD_CROSS_TRACK_STRATEGY - typedef cross_track this_type; typedef typename services::comparable_type @@ -417,7 +511,6 @@ public : return_type const a = cstrategy.apply(p, sp1, sp2); return_type const c = return_type(2.0) * asin(math::sqrt(a)); return c * radius(); -#endif // BOOST_GEOMETRY_USE_OLD_CROSS_TRACK_STRATEGY } inline typename Strategy::radius_type radius() const @@ -450,15 +543,10 @@ struct return_type, P, PS> template struct comparable_type > { -#if 0 - // There is no shortcut, so the strategy itself is its comparable type - typedef cross_track type; -#else typedef comparable::cross_track < CalculationType, typename comparable_type::type > type; -#endif }; @@ -474,7 +562,8 @@ struct get_comparable > cross_track >::type comparable_type; public : - static inline comparable_type apply(cross_track const& strategy) + static inline comparable_type + apply(cross_track const& strategy) { return comparable_type(strategy.radius()); } @@ -485,7 +574,8 @@ template < typename CalculationType, typename Strategy, - typename P, typename PS + typename P, + typename PS > struct result_from_distance, P, PS> { @@ -496,7 +586,8 @@ private : >::template return_type::type return_type; public : template - static inline return_type apply(cross_track const& , T const& distance) + static inline return_type + apply(cross_track const& , T const& distance) { return distance; } @@ -511,9 +602,18 @@ struct tag > }; -template -struct return_type, P1, P2> - : comparable::cross_track::template return_type +template +< + typename RadiusType, + typename CalculationType, + typename P, + typename PS +> +struct return_type, P, PS> + : comparable::cross_track + < + RadiusType, CalculationType + >::template return_type {}; @@ -537,15 +637,25 @@ public : }; -template -struct result_from_distance, P1, P2> +template +< + typename RadiusType, + typename CalculationType, + typename P, + typename PS +> +struct result_from_distance + < + comparable::cross_track, P, PS + > { private : typedef comparable::cross_track strategy_type; - typedef typename return_type::type return_type; + typedef typename return_type::type return_type; public : template - static inline return_type apply(strategy_type const& strategy, T const& distance) + static inline return_type apply(strategy_type const& strategy, + T const& distance) { return_type const s = sin( (distance / strategy.radius()) / return_type(2.0) ); @@ -629,17 +739,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