boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-01-31 08:12:13 +00:00
Code Issues Projects Releases Wiki Activity

[formulas] [strategies] Distance point-segment use the meridian formula. Use new static version of distance strategy.

Browse Source
This commit is contained in:
Vissarion Fysikopoulos
2017-11-01 17:42:48 +02:00
parent b17ad43f7f
commit f458d8d28e
5 changed files with 112 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files

45
include/boost/geometry/formulas/distance_point_segment.hpp
View File

@@ -18,6 +18,10 @@
#define BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS 100
#endif
#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG
#include <iostream>
#endif
/*!
\brief Algorithm to compute the distance between a segment and a point using
direct and inverse geodesic problems as subroutines. The algorithm
@@ -33,19 +37,20 @@ template
<
typename CT,
typename Units,
template <typename, bool, bool, bool, bool ,bool> class Inverse,
template <typename, bool, bool, bool, bool> class Direct,
typename FormulaPolicy,
//template <typename, bool, bool, bool, bool ,bool> class Inverse,
//template <typename, bool, bool, bool, bool> class Direct,
bool EnableClosestPoint = false
>
class distance_point_segment{
public:
typedef Inverse<CT, true, false, false, true, true> inverse_distance_quantities_type;
typedef Inverse<CT, true, false, false, false, false> inverse_distance_type;
typedef Inverse<CT, false, true, false, false, false> inverse_azimuth_type;
typedef Inverse<CT, false, true, true, false, false> inverse_azimuth_reverse_type;
typedef Direct<CT, true, false, false, false> direct_distance_type;
typedef typename FormulaPolicy::template inverse<CT, true, false, false, true, true> inverse_distance_quantities_type;
typedef typename FormulaPolicy::template inverse<CT, true, false, false, false, false> inverse_distance_type;
typedef typename FormulaPolicy::template inverse<CT, false, true, false, false, false> inverse_azimuth_type;
typedef typename FormulaPolicy::template inverse<CT, false, true, true, false, false> inverse_azimuth_reverse_type;
typedef typename FormulaPolicy::template direct<CT, true, false, false, false> direct_distance_type;
struct result_distance_point_segment
{
@@ -80,9 +85,16 @@ public:
CT lon2, CT lat2, //p2
Spheroid const& spheroid)
{
CT distance = inverse_distance_quantities_type::apply(lon1, lat1,
lon2, lat2,
spheroid).distance;
//CT distance = inverse_distance_type::apply(lon1, lat1,
//lon2, lat2,
//spheroid).distance;
//geometry::model::point<CT, 2, geometry::cs::geographic<geometry::radian> > point;
// CT distance = geometry::distance(point(lon1, lat1), point(lon2, lat2), geometry::strategy::distance::andoyer<>());
//geometry::strategy::distance::geographic<> str(spheroid);
CT distance = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT>::apply(lon1, lat1, lon2, lat2, spheroid);
return non_iterative_case(lon1, lat1, distance);
}
@@ -153,14 +165,25 @@ public:
}
return non_iterative_case(lon3, lat1, lon3, lat3, spheroid);
}
/*
CT d1 = inverse_distance_type::apply(lon1, lat1,
lon3, lat3, spheroid).distance;
CT d3 = inverse_distance_type::apply(lon1, lat1,
lon2, lat2, spheroid).distance;
*/
CT d1 = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT>
::apply(lon1, lat1, lon3, lat3, spheroid);
CT d3 = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT>
::apply(lon1, lat1, lon2, lat2, spheroid);
if (geometry::math::equals(d3, c0))
{
#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG
std::cout << "Degenerate segment" << std::endl;
std::cout << "d1=" << d1 << std::endl;
#endif
return non_iterative_case(lon1, lat2, d1);
}

50
include/boost/geometry/formulas/elliptic_arc_length.hpp
View File

@@ -22,7 +22,6 @@
#include <boost/geometry/formulas/flattening.hpp>
namespace boost { namespace geometry { namespace formula
{
@@ -36,6 +35,53 @@ class elliptic_arc_length
public :
struct result
{
result()
: distance(0)
, meridian(false)
{}
CT distance;
bool meridian;
};
template <typename T, typename Spheroid>
static result apply(T lon1, T lat1, T lon2, T lat2, Spheroid const& spheroid)
{
result res;
CT c0 = 0;
CT pi = math::pi<CT>();
CT half_pi = pi/CT(2);
CT diff = math::longitude_distance_signed<geometry::radian>(lon1, lon2);
if (math::equals(diff, c0))
{
// single meridian not crossing pole
if (lat1 > lat2)
{
std::swap(lat1, lat2);
}
res.distance = apply(lat2, spheroid) - apply(lat1, spheroid);
res.meridian = true;
}
if (math::equals(math::abs(diff), pi))
{
// meridian crosses pole
CT lat_sign = 1;
if (lat1+lat2 < c0)
{
lat_sign = CT(-1);
}
res.distance = math::abs(lat_sign * CT(2) * apply(half_pi, spheroid)
- apply(lat1, spheroid) - apply(lat2, spheroid));
res.meridian = true;
}
return res;
}
// Distance computation on meridians using series approximations
// to elliptic integrals. Formula to compute distance from lattitude 0 to lat
// https://en.wikipedia.org/wiki/Meridian_arc
@@ -48,7 +94,7 @@ public :
CT n = f / (CT(2) - f);
CT M = a/(1+n);
CT C0 = 1;
if (Order == 0)
{
return M * C0 * lat;

63
include/boost/geometry/strategies/geographic/distance.hpp
View File

@@ -73,6 +73,29 @@ public :
: m_spheroid(spheroid)
{}
template <typename CT>
static inline CT apply(CT lon1, CT lat1, CT lon2, CT lat2,
Spheroid const& spheroid)
{
typedef typename formula::elliptic_arc_length
<
CT, strategy::default_order<FormulaPolicy>::value
> elliptic_arc_length;
typename elliptic_arc_length::result res =
elliptic_arc_length::apply(lon1, lat1, lon2, lat2, spheroid);
if (res.meridian)
{
return res.distance;
}
return FormulaPolicy::template inverse
<
CT, true, false, false, false, false
>::apply(lon1, lat1, lon2, lat2, spheroid).distance;
}
template <typename Point1, typename Point2>
inline typename calculation_type<Point1, Point2>::type
apply(Point1 const& point1, Point2 const& point2) const
@@ -84,45 +107,7 @@ public :
CT lon2 = get_as_radian<0>(point2);
CT lat2 = get_as_radian<1>(point2);
CT c0 = 0;
CT pi = math::pi<CT>();
CT half_pi = pi/CT(2);
CT diff = math::longitude_distance_signed<geometry::radian>(lon1, lon2);
typedef typename formula::elliptic_arc_length
<
CT, strategy::default_order<FormulaPolicy>::value
> elliptic_arc_length;
if (math::equals(diff, c0))
{
// single meridian not crossing pole
if (lat1 > lat2)
{
std::swap(lat1, lat2);
}
return elliptic_arc_length::apply(lat2, m_spheroid)
- elliptic_arc_length::apply(lat1, m_spheroid);
}
if (math::equals(math::abs(diff), pi))
{
// meridian crosses pole
CT lat_sign = 1;
if (lat1+lat2 < c0)
{
lat_sign = CT(-1);
}
return math::abs(lat_sign * CT(2) *
elliptic_arc_length::apply(half_pi, m_spheroid)
- elliptic_arc_length::apply(lat1, m_spheroid)
- elliptic_arc_length::apply(lat2, m_spheroid));
}
return FormulaPolicy::template inverse
<
CT, true, false, false, false, false
>::apply(lon1, lat1, lon2, lat2, m_spheroid).distance;
return apply(lon1, lat1, lon2, lat2, m_spheroid);
}
inline Spheroid const& model() const

3
include/boost/geometry/strategies/geographic/distance_cross_track.hpp
View File

@@ -106,8 +106,7 @@ public :
<
CT,
units_type,
FormulaPolicy::template inverse,
FormulaPolicy::template direct
FormulaPolicy
>::apply(get<0>(sp1), get<1>(sp1),
get<0>(sp2), get<1>(sp2),
get<0>(p), get<1>(p),

8
test/algorithms/distance/distance_geo_pl_l.cpp
View File

@@ -8,6 +8,9 @@
// Licensed under the Boost Software License version 1.0.
// http://www.boost.org/users/license.html
#define BOOST_GEOMETRY_TEST_DEBUG
#define BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG
#include <iostream>
#ifndef BOOST_TEST_MODULE
@@ -318,13 +321,12 @@ void test_distance_point_segment(Strategy_pp const& strategy_pp,
pp_distance("POINT(90 -89)", "POINT(90 90)", strategy_pp),
strategy_ps);
// degenerate segment
tester::apply("p-s-14",
tester::apply("p-s-15",
"POINT(90 90)",
"SEGMENT(0.5 -90,175.5 -90)",
pp_distance("POINT(90 -90)", "POINT(90 90)", strategy_pp),
pp_distance("POINT(0.5 -90)", "POINT(90 90)", strategy_pp),
strategy_ps);
}
//===========================================================================