boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-02-11 11:52:11 +00:00
Code Issues Projects Releases Wiki Activity

[strategies][distance][spherical][cross track] update comments' section and make it

Browse Source 
a stand-alone explanatory document; polish code; remove old code; fix header includes;
fix long lines;
This commit is contained in:
Menelaos Karavelas
2014-11-19 23:52:59 +02:00
parent 445fa3fd8e
commit 43d7f82491
1 changed files with 212 additions and 111 deletions
Download Patch File Download Diff File Expand all files Collapse all files

323
include/boost/geometry/strategies/spherical/distance_cross_track.hpp
View File

@@ -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,15 +32,15 @@
#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
{
@@ -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<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 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<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);
@@ -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<Strategy, Point, PointOfSegment>)
);
#endif
typedef typename return_type<Point, PointOfSegment>::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<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;
#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<CalculationType, Strategy> 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<cross_track<CalculationType, Strategy>, P, PS>
template <typename CalculationType, typename Strategy>
struct comparable_type<cross_track<CalculationType, Strategy> >
{
#if 0
// There is no shortcut, so the strategy itself is its comparable type
typedef cross_track<CalculationType, Strategy> type;
#else
typedef comparable::cross_track
<
CalculationType, typename comparable_type<Strategy>::type
> type;
#endif
};
@@ -474,7 +562,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());
}
@@ -485,7 +574,8 @@ template
<
typename CalculationType,
typename Strategy,
typename P, typename PS
typename P,
typename PS
>
struct result_from_distance<cross_track<CalculationType, Strategy>, P, PS>
{
@@ -496,7 +586,8 @@ private :
>::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;
}
@@ -511,9 +602,18 @@ struct tag<comparable::cross_track<RadiusType, CalculationType> >
};
template <typename RadiusType, typename CalculationType, typename P1, typename P2>
struct return_type<comparable::cross_track<RadiusType, CalculationType>, P1, P2>
: comparable::cross_track<RadiusType, CalculationType>::template return_type<P1, P2>
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>
{};
@@ -537,15 +637,25 @@ public :
};
template <typename RadiusType, typename CalculationType, typename P1, typename P2>
struct result_from_distance<comparable::cross_track<RadiusType, CalculationType>, P1, P2>
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, P1, P2>::type return_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)
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