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