summaryrefslogtreecommitdiff
path: root/boost/geometry/strategies
diff options
context:
space:
mode:
Diffstat (limited to 'boost/geometry/strategies')
-rw-r--r--boost/geometry/strategies/agnostic/buffer_distance_asymmetric.hpp4
-rw-r--r--boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp5
-rw-r--r--boost/geometry/strategies/agnostic/hull_graham_andrew.hpp2
-rw-r--r--boost/geometry/strategies/agnostic/point_in_poly_winding.hpp12
-rw-r--r--boost/geometry/strategies/agnostic/relate.hpp31
-rw-r--r--boost/geometry/strategies/agnostic/side_by_azimuth.hpp87
-rw-r--r--boost/geometry/strategies/agnostic/simplify_douglas_peucker.hpp51
-rw-r--r--boost/geometry/strategies/cartesian/buffer_end_round.hpp12
-rw-r--r--boost/geometry/strategies/cartesian/buffer_join_miter.hpp2
-rw-r--r--boost/geometry/strategies/cartesian/buffer_join_round.hpp48
-rw-r--r--boost/geometry/strategies/cartesian/buffer_point_circle.hpp14
-rw-r--r--boost/geometry/strategies/cartesian/cart_intersect.hpp169
-rw-r--r--boost/geometry/strategies/cartesian/centroid_average.hpp18
-rw-r--r--boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp9
-rw-r--r--boost/geometry/strategies/cartesian/side_by_triangle.hpp179
-rw-r--r--boost/geometry/strategies/cartesian/side_of_intersection.hpp119
-rw-r--r--boost/geometry/strategies/comparable_distance_result.hpp20
-rw-r--r--boost/geometry/strategies/concepts/distance_concept.hpp66
-rw-r--r--boost/geometry/strategies/distance_result.hpp36
-rw-r--r--boost/geometry/strategies/geographic/distance_andoyer.hpp224
-rw-r--r--boost/geometry/strategies/geographic/distance_vincenty.hpp161
-rw-r--r--boost/geometry/strategies/geographic/mapping_ssf.hpp185
-rw-r--r--boost/geometry/strategies/spherical/distance_cross_track.hpp531
-rw-r--r--boost/geometry/strategies/spherical/distance_haversine.hpp2
-rw-r--r--boost/geometry/strategies/spherical/side_by_cross_track.hpp42
-rw-r--r--boost/geometry/strategies/spherical/ssf.hpp109
-rw-r--r--boost/geometry/strategies/strategies.hpp3
27 files changed, 1870 insertions, 271 deletions
diff --git a/boost/geometry/strategies/agnostic/buffer_distance_asymmetric.hpp b/boost/geometry/strategies/agnostic/buffer_distance_asymmetric.hpp
index 7b7cd18..446d2f0 100644
--- a/boost/geometry/strategies/agnostic/buffer_distance_asymmetric.hpp
+++ b/boost/geometry/strategies/agnostic/buffer_distance_asymmetric.hpp
@@ -9,6 +9,8 @@
#ifndef BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_BUFFER_DISTANCE_ASYMMETRIC_HPP
#define BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_BUFFER_DISTANCE_ASYMMETRIC_HPP
+#include <boost/core/ignore_unused.hpp>
+
#include <boost/geometry/strategies/buffer.hpp>
#include <boost/geometry/util/math.hpp>
@@ -79,6 +81,8 @@ public :
inline NumericType max_distance(JoinStrategy const& join_strategy,
EndStrategy const& end_strategy) const
{
+ boost::ignore_unused(join_strategy, end_strategy);
+
NumericType const left = geometry::math::abs(m_left);
NumericType const right = geometry::math::abs(m_right);
NumericType const dist = (std::max)(left, right);
diff --git a/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp b/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp
index bc0c46f..73bd21a 100644
--- a/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp
+++ b/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp
@@ -9,6 +9,9 @@
#ifndef BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_BUFFER_DISTANCE_SYMMETRIC_HPP
#define BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_BUFFER_DISTANCE_SYMMETRIC_HPP
+
+#include <boost/core/ignore_unused.hpp>
+
#include <boost/geometry/strategies/buffer.hpp>
#include <boost/geometry/util/math.hpp>
@@ -76,6 +79,8 @@ public :
inline NumericType max_distance(JoinStrategy const& join_strategy,
EndStrategy const& end_strategy) const
{
+ boost::ignore_unused(join_strategy, end_strategy);
+
NumericType const dist = geometry::math::abs(m_distance);
return (std::max)(join_strategy.max_distance(dist),
end_strategy.max_distance(dist));
diff --git a/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp b/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp
index a960a6f..1d59e13 100644
--- a/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp
+++ b/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp
@@ -364,7 +364,7 @@ private:
// count describes a closed case but comparison with min size of closed
// gives the result compatible also with open
// here core_detail::closure::minimum_ring_size<closed> could be used
- if ( count < 4 )
+ if (count < 4)
{
// there should be only one missing
*out++ = *boost::begin(first);
diff --git a/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp b/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp
index f4ed7a6..641533f 100644
--- a/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp
+++ b/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp
@@ -73,7 +73,7 @@ struct winding_side_equal
PointOfSegment ss1, ss2;
set<1-D>(ss1, get<1-D>(se));
set<1-D>(ss2, get<1-D>(se));
- if ( count > 0 ) // UP
+ if (count > 0) // UP
{
set<D>(ss1, 0);
set<D>(ss2, 1);
@@ -127,7 +127,7 @@ struct winding_side_between
set<1-D>(ss1, get<1-D>(s1));
set<1-D>(ss2, get<1-D>(s1));
- if ( count > 0 ) // UP
+ if (count > 0) // UP
{
set<D>(ss1, 0);
set<D>(ss2, 1);
@@ -140,9 +140,9 @@ struct winding_side_between
int const seg_side = strategy_side_type::apply(ss1, ss2, s2);
- if ( seg_side != 0 ) // segment not vertical
+ if (seg_side != 0) // segment not vertical
{
- if ( strategy_side_type::apply(ss1, ss2, point) == -seg_side ) // point on the opposite side than s2
+ if (strategy_side_type::apply(ss1, ss2, point) == -seg_side) // point on the opposite side than s2
{
return -seg_side;
}
@@ -151,7 +151,7 @@ struct winding_side_between
set<1-D>(ss1, get<1-D>(s2));
set<1-D>(ss2, get<1-D>(s2));
- if ( strategy_side_type::apply(ss1, ss2, point) == seg_side ) // point behind s2
+ if (strategy_side_type::apply(ss1, ss2, point) == seg_side) // point behind s2
{
return seg_side;
}
@@ -308,7 +308,7 @@ public :
if (count != 0)
{
int side = 0;
- if ( count == 1 || count == -1 )
+ if (count == 1 || count == -1)
{
side = winding_side_equal<cs_t>
::template apply<1>(point, eq1 ? s1 : s2, count);
diff --git a/boost/geometry/strategies/agnostic/relate.hpp b/boost/geometry/strategies/agnostic/relate.hpp
index 318047f..9e87532 100644
--- a/boost/geometry/strategies/agnostic/relate.hpp
+++ b/boost/geometry/strategies/agnostic/relate.hpp
@@ -20,10 +20,9 @@ namespace boost { namespace geometry
namespace strategy { namespace relate
{
-template <typename StaticMask>
+template <typename Geometry1, typename Geometry2, typename StaticMask>
struct relate
{
- template <typename Geometry1, typename Geometry2>
static inline bool apply(Geometry1 const& geometry1, Geometry2 const& geometry2)
{
return detail::relate::relate<StaticMask>(geometry1, geometry2);
@@ -44,13 +43,23 @@ namespace services
template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2, typename AnyCS>
struct default_strategy<AnyTag1, AnyTag2, AnyTag1, AnyTag2, AnyCS, AnyCS, Geometry1, Geometry2>
{
- typedef strategy::relate::relate<detail::relate::static_mask_within> type;
+ typedef strategy::relate::relate
+ <
+ Geometry1,
+ Geometry2,
+ detail::relate::static_mask_within
+ > type;
};
template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2, typename AnyCS>
struct default_strategy<AnyTag1, AnyTag2, AnyTag1, areal_tag, AnyCS, AnyCS, Geometry1, Geometry2>
{
- typedef strategy::relate::relate<detail::relate::static_mask_within> type;
+ typedef strategy::relate::relate
+ <
+ Geometry1,
+ Geometry2,
+ detail::relate::static_mask_within
+ > type;
};
@@ -71,13 +80,23 @@ namespace strategy { namespace covered_by { namespace services
template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2, typename AnyCS>
struct default_strategy<AnyTag1, AnyTag2, AnyTag1, AnyTag2, AnyCS, AnyCS, Geometry1, Geometry2>
{
- typedef strategy::relate::relate<detail::relate::static_mask_covered_by> type;
+ typedef strategy::relate::relate
+ <
+ Geometry1,
+ Geometry2,
+ detail::relate::static_mask_covered_by
+ > type;
};
template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2, typename AnyCS>
struct default_strategy<AnyTag1, AnyTag2, AnyTag1, areal_tag, AnyCS, AnyCS, Geometry1, Geometry2>
{
- typedef strategy::relate::relate<detail::relate::static_mask_covered_by> type;
+ typedef strategy::relate::relate
+ <
+ Geometry1,
+ Geometry2,
+ detail::relate::static_mask_covered_by
+ > type;
};
diff --git a/boost/geometry/strategies/agnostic/side_by_azimuth.hpp b/boost/geometry/strategies/agnostic/side_by_azimuth.hpp
new file mode 100644
index 0000000..14c69a0
--- /dev/null
+++ b/boost/geometry/strategies/agnostic/side_by_azimuth.hpp
@@ -0,0 +1,87 @@
+// Boost.Geometry
+
+// Copyright (c) 2007-2012 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 Adam Wulkiewicz, 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
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_SIDE_BY_AZIMUTH_HPP
+#define BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_SIDE_BY_AZIMUTH_HPP
+
+#include <boost/mpl/if.hpp>
+#include <boost/type_traits.hpp>
+#include <boost/core/ignore_unused.hpp>
+
+#include <boost/geometry/core/cs.hpp>
+#include <boost/geometry/core/access.hpp>
+#include <boost/geometry/core/radian_access.hpp>
+
+#include <boost/geometry/algorithms/detail/azimuth.hpp>
+
+#include <boost/geometry/util/math.hpp>
+#include <boost/geometry/util/promote_floating_point.hpp>
+#include <boost/geometry/util/select_calculation_type.hpp>
+
+#include <boost/geometry/strategies/side.hpp>
+//#include <boost/geometry/strategies/concepts/side_concept.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+
+namespace strategy { namespace side
+{
+
+/*!
+\brief Check at which side of a segment a point lies
+ left of segment (> 0), right of segment (< 0), on segment (0)
+\ingroup strategies
+\tparam Model Reference model of coordinate system.
+\tparam CalculationType \tparam_calculation
+ */
+template <typename Model, typename CalculationType = void>
+class side_by_azimuth
+{
+public:
+ side_by_azimuth(Model const& model = Model())
+ : m_model(model)
+ {}
+
+ template <typename P1, typename P2, typename P>
+ inline int apply(P1 const& p1, P2 const& p2, P const& p)
+ {
+ typedef typename promote_floating_point
+ <
+ typename select_calculation_type_alt
+ <
+ CalculationType,
+ P1, P2, P
+ >::type
+ >::type calc_t;
+
+ calc_t d1 = 0.001;
+ calc_t crs_AD = geometry::detail::azimuth<calc_t>(p1, p, m_model);
+ calc_t crs_AB = geometry::detail::azimuth<calc_t>(p1, p2, m_model);
+ calc_t XTD = asin(sin(d1) * sin(crs_AD - crs_AB));
+
+ return math::equals(XTD, 0) ? 0 : XTD < 0 ? 1 : -1;
+ }
+
+private:
+ Model m_model;
+};
+
+}} // namespace strategy::side
+
+
+}} // namespace boost::geometry
+
+
+#endif // BOOST_GEOMETRY_STRATEGIES_AGNOSTIC_SIDE_BY_AZIMUTH_HPP
diff --git a/boost/geometry/strategies/agnostic/simplify_douglas_peucker.hpp b/boost/geometry/strategies/agnostic/simplify_douglas_peucker.hpp
index 8ad3bbc..99e7d9b 100644
--- a/boost/geometry/strategies/agnostic/simplify_douglas_peucker.hpp
+++ b/boost/geometry/strategies/agnostic/simplify_douglas_peucker.hpp
@@ -1,8 +1,13 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
-// Copyright (c) 1995, 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 1995, 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 1995 Maarten Hilferink, Amsterdam, the Netherlands
+// This file was modified by Oracle on 2015.
+// Modifications copyright (c) 2015, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
+
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
@@ -15,6 +20,9 @@
#include <cstddef>
+#ifdef BOOST_GEOMETRY_DEBUG_DOUGLAS_PEUCKER
+#include <iostream>
+#endif
#include <vector>
#include <boost/range.hpp>
@@ -23,10 +31,7 @@
#include <boost/geometry/strategies/distance.hpp>
-
-//#define GL_DEBUG_DOUGLAS_PEUCKER
-
-#ifdef GL_DEBUG_DOUGLAS_PEUCKER
+#ifdef BOOST_GEOMETRY_DEBUG_DOUGLAS_PEUCKER
#include <boost/geometry/io/dsv/write.hpp>
#endif
@@ -126,7 +131,7 @@ namespace detail
// because we want to consider a candidate point in between
if (size <= 2)
{
-#ifdef GL_DEBUG_DOUGLAS_PEUCKER
+#ifdef BOOST_GEOMETRY_DEBUG_DOUGLAS_PEUCKER
if (begin != end)
{
std::cout << "ignore between " << dsv(begin->p)
@@ -140,7 +145,7 @@ namespace detail
iterator_type last = end - 1;
-#ifdef GL_DEBUG_DOUGLAS_PEUCKER
+#ifdef BOOST_GEOMETRY_DEBUG_DOUGLAS_PEUCKER
std::cout << "find between " << dsv(begin->p)
<< " and " << dsv(last->p)
<< " size=" << size << std::endl;
@@ -155,7 +160,7 @@ namespace detail
{
distance_type dist = ps_distance_strategy.apply(it->p, begin->p, last->p);
-#ifdef GL_DEBUG_DOUGLAS_PEUCKER
+#ifdef BOOST_GEOMETRY_DEBUG_DOUGLAS_PEUCKER
std::cout << "consider " << dsv(it->p)
<< " at " << double(dist)
<< ((dist > max_dist) ? " maybe" : " no")
@@ -173,7 +178,7 @@ namespace detail
// and handle segments in between recursively
if ( less()(max_dist, md) )
{
-#ifdef GL_DEBUG_DOUGLAS_PEUCKER
+#ifdef BOOST_GEOMETRY_DEBUG_DOUGLAS_PEUCKER
std::cout << "use " << dsv(candidate->p) << std::endl;
#endif
@@ -193,6 +198,10 @@ namespace detail
OutputIterator out,
distance_type max_distance) const
{
+#ifdef BOOST_GEOMETRY_DEBUG_DOUGLAS_PEUCKER
+ std::cout << "max distance: " << max_distance
+ << std::endl << std::endl;
+#endif
distance_strategy_type strategy;
// Copy coordinates, a vector of references to all points
@@ -228,8 +237,6 @@ namespace detail
}
};
-
-
}
#endif // DOXYGEN_NO_DETAIL
@@ -269,18 +276,28 @@ public :
PointDistanceStrategy
>::distance_type distance_type;
- typedef distance_type return_type;
-
template <typename Range, typename OutputIterator>
static inline OutputIterator apply(Range const& range,
OutputIterator out,
- distance_type max_distance)
+ distance_type const& max_distance)
{
- return detail::douglas_peucker
+ namespace services = strategy::distance::services;
+
+ typedef typename services::comparable_type
<
- Point,
PointDistanceStrategy
- >().apply(range, out, max_distance);
+ >::type comparable_distance_strategy_type;
+
+ return detail::douglas_peucker
+ <
+ Point, comparable_distance_strategy_type
+ >().apply(range, out,
+ services::result_from_distance
+ <
+ comparable_distance_strategy_type, Point, Point
+ >::apply(comparable_distance_strategy_type(),
+ max_distance)
+ );
}
};
diff --git a/boost/geometry/strategies/cartesian/buffer_end_round.hpp b/boost/geometry/strategies/cartesian/buffer_end_round.hpp
index 74780d6..a233f1c 100644
--- a/boost/geometry/strategies/cartesian/buffer_end_round.hpp
+++ b/boost/geometry/strategies/cartesian/buffer_end_round.hpp
@@ -95,8 +95,9 @@ public :
//! \brief Constructs the strategy
//! \param points_per_circle points which would be used for a full circle
+ //! (if points_per_circle is smaller than 4, it is internally set to 4)
explicit inline end_round(std::size_t points_per_circle = 90)
- : m_points_per_circle(points_per_circle)
+ : m_points_per_circle((points_per_circle < 4u) ? 4u : points_per_circle)
{}
#ifndef DOXYGEN_SHOULD_SKIP_THIS
@@ -106,7 +107,7 @@ public :
inline void apply(Point const& penultimate_point,
Point const& perp_left_point,
Point const& ultimate_point,
- Point const& ,
+ Point const& perp_right_point,
buffer_side_selector side,
DistanceStrategy const& distance,
RangeOut& range_out) const
@@ -142,6 +143,13 @@ public :
set<1>(shifted_point, get<1>(ultimate_point) + dist_half_diff * sin(alpha));
generate_points(shifted_point, alpha, (dist_left + dist_right) / two, range_out);
}
+
+ if (m_points_per_circle % 2 == 1)
+ {
+ // For a half circle, if the number of points is not even,
+ // we should insert the end point too, to generate a full cap
+ range_out.push_back(perp_right_point);
+ }
}
template <typename NumericType>
diff --git a/boost/geometry/strategies/cartesian/buffer_join_miter.hpp b/boost/geometry/strategies/cartesian/buffer_join_miter.hpp
index 8fcf3b9..99ec805 100644
--- a/boost/geometry/strategies/cartesian/buffer_join_miter.hpp
+++ b/boost/geometry/strategies/cartesian/buffer_join_miter.hpp
@@ -35,6 +35,8 @@ namespace strategy { namespace buffer
their length. The miter is not changed to a bevel form (as done in some
other software), it is just adapted to the specified miter_limit but keeps
its miter form.
+ If the buffer distance is 5.0, and the miter limit is 2.0, generated points
+ will be located at a distance of at most 10.0 (2*5) units.
This strategy is only applicable for Cartesian coordinate systems.
\qbk{
diff --git a/boost/geometry/strategies/cartesian/buffer_join_round.hpp b/boost/geometry/strategies/cartesian/buffer_join_round.hpp
index 9e467c8..9ec51cd 100644
--- a/boost/geometry/strategies/cartesian/buffer_join_round.hpp
+++ b/boost/geometry/strategies/cartesian/buffer_join_round.hpp
@@ -1,6 +1,6 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
-// Copyright (c) 2012-2014 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2012-2015 Barend Gehrels, Amsterdam, the Netherlands.
// 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,6 +9,8 @@
#ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_ROUND_HPP
#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_ROUND_HPP
+#include <algorithm>
+
#include <boost/assert.hpp>
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/policies/compare.hpp>
@@ -69,34 +71,38 @@ private :
DistanceType const& buffer_distance,
RangeOut& range_out) const
{
- PromotedType dx1 = get<0>(perp1) - get<0>(vertex);
- PromotedType dy1 = get<1>(perp1) - get<1>(vertex);
- PromotedType dx2 = get<0>(perp2) - get<0>(vertex);
- PromotedType dy2 = get<1>(perp2) - get<1>(vertex);
+ PromotedType const dx1 = get<0>(perp1) - get<0>(vertex);
+ PromotedType const dy1 = get<1>(perp1) - get<1>(vertex);
+ PromotedType const dx2 = get<0>(perp2) - get<0>(vertex);
+ PromotedType const dy2 = get<1>(perp2) - get<1>(vertex);
- BOOST_ASSERT(buffer_distance != 0);
+ PromotedType const two = 2.0;
+ PromotedType const two_pi = two * geometry::math::pi<PromotedType>();
- dx1 /= buffer_distance;
- dy1 /= buffer_distance;
- dx2 /= buffer_distance;
- dy2 /= buffer_distance;
+ PromotedType const angle1 = atan2(dy1, dx1);
+ PromotedType angle2 = atan2(dy2, dx2);
+ while (angle2 > angle1)
+ {
+ angle2 -= two_pi;
+ }
+ PromotedType const angle_diff = angle1 - angle2;
- PromotedType angle_diff = acos(dx1 * dx2 + dy1 * dy2);
+ // Divide the angle into an integer amount of steps to make it
+ // visually correct also for a low number of points / circle
- PromotedType two = 2.0;
- PromotedType steps = m_points_per_circle;
- int n = boost::numeric_cast<int>(steps * angle_diff
- / (two * geometry::math::pi<PromotedType>()));
+ // If a full circle is divided into 3 parts (e.g. angle is 125),
+ // the one point in between must still be generated
+ // The calculation below:
+ // - generates 1 point in between for an angle of 125 based on 3 points
+ // - generates 0 points in between for an angle of 90 based on 4 points
- if (n <= 1)
- {
- return;
- }
+ int const n = (std::max)(static_cast<int>(
+ ceil(m_points_per_circle * angle_diff / two_pi)), 1);
- PromotedType const angle1 = atan2(dy1, dx1);
- PromotedType diff = angle_diff / PromotedType(n);
+ PromotedType const diff = angle_diff / static_cast<PromotedType>(n);
PromotedType a = angle1 - diff;
+ // Walk to n - 1 to avoid generating the last point
for (int i = 0; i < n - 1; i++, a -= diff)
{
Point p;
diff --git a/boost/geometry/strategies/cartesian/buffer_point_circle.hpp b/boost/geometry/strategies/cartesian/buffer_point_circle.hpp
index f64a82d..86ebc43 100644
--- a/boost/geometry/strategies/cartesian/buffer_point_circle.hpp
+++ b/boost/geometry/strategies/cartesian/buffer_point_circle.hpp
@@ -1,5 +1,12 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
-// Copyright (c) 2012-2014 Barend Gehrels, Amsterdam, the Netherlands.
+
+// Copyright (c) 2012-2015 Barend Gehrels, Amsterdam, the Netherlands.
+
+// This file was modified by Oracle on 2015.
+// Modifications copyright (c) 2015, 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
// http://www.boost.org/LICENSE_1_0.txt)
@@ -45,9 +52,10 @@ class point_circle
{
public :
//! \brief Constructs the strategy
- //! \param count number of points for the created circle
+ //! \param count number of points for the created circle (if count
+ //! is smaller than 3, count is internally set to 3)
explicit point_circle(std::size_t count = 90)
- : m_count(count)
+ : m_count((count < 3u) ? 3u : count)
{}
#ifndef DOXYGEN_SHOULD_SKIP_THIS
diff --git a/boost/geometry/strategies/cartesian/cart_intersect.hpp b/boost/geometry/strategies/cartesian/cart_intersect.hpp
index 66af2d2..a7bd385 100644
--- a/boost/geometry/strategies/cartesian/cart_intersect.hpp
+++ b/boost/geometry/strategies/cartesian/cart_intersect.hpp
@@ -1,7 +1,13 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
-// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
-// Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland.
+// Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2013-2014 Adam Wulkiewicz, Lodz, Poland.
+
+// 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
+// Contributed and/or modified by Adam Wulkiewicz, 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
@@ -194,7 +200,11 @@ struct relate_cartesian_segments
get<1>(robust_b1) - get<1>(robust_a1),
robust_db0, robust_db);
- if (robust_da0 == 0)
+ math::detail::equals_factor_policy<robust_coordinate_type>
+ policy(robust_dx_a, robust_dy_a, robust_dx_b, robust_dy_b);
+ robust_coordinate_type const zero = 0;
+ if (math::detail::equals_by_policy(robust_da0, zero, policy)
+ || math::detail::equals_by_policy(robust_db0, zero, policy))
{
// If this is the case, no rescaling is done for FP precision.
// We set it to collinear, but it indicates a robustness issue.
@@ -211,25 +221,31 @@ struct relate_cartesian_segments
if (collinear)
{
- bool const collinear_use_first
- = geometry::math::abs(robust_dx_a) + geometry::math::abs(robust_dx_b)
- >= geometry::math::abs(robust_dy_a) + geometry::math::abs(robust_dy_b);
-
- // Degenerate cases: segments of single point, lying on other segment, are not disjoint
- // This situation is collinear too
-
- if (collinear_use_first)
+ std::pair<bool, bool> const collinear_use_first
+ = is_x_more_significant(geometry::math::abs(robust_dx_a),
+ geometry::math::abs(robust_dy_a),
+ geometry::math::abs(robust_dx_b),
+ geometry::math::abs(robust_dy_b),
+ a_is_point, b_is_point);
+
+ if (collinear_use_first.second)
{
- return relate_collinear<0, ratio_type>(a, b,
- robust_a1, robust_a2, robust_b1, robust_b2,
- a_is_point, b_is_point);
- }
- else
- {
- // Y direction contains larger segments (maybe dx is zero)
- return relate_collinear<1, ratio_type>(a, b,
- robust_a1, robust_a2, robust_b1, robust_b2,
- a_is_point, b_is_point);
+ // Degenerate cases: segments of single point, lying on other segment, are not disjoint
+ // This situation is collinear too
+
+ if (collinear_use_first.first)
+ {
+ return relate_collinear<0, ratio_type>(a, b,
+ robust_a1, robust_a2, robust_b1, robust_b2,
+ a_is_point, b_is_point);
+ }
+ else
+ {
+ // Y direction contains larger segments (maybe dx is zero)
+ return relate_collinear<1, ratio_type>(a, b,
+ robust_a1, robust_a2, robust_b1, robust_b2,
+ a_is_point, b_is_point);
+ }
}
}
@@ -237,6 +253,40 @@ struct relate_cartesian_segments
}
private:
+ // first is true if x is more significant
+ // second is true if the more significant difference is not 0
+ template <typename RobustCoordinateType>
+ static inline std::pair<bool, bool>
+ is_x_more_significant(RobustCoordinateType const& abs_robust_dx_a,
+ RobustCoordinateType const& abs_robust_dy_a,
+ RobustCoordinateType const& abs_robust_dx_b,
+ RobustCoordinateType const& abs_robust_dy_b,
+ bool const a_is_point,
+ bool const b_is_point)
+ {
+ //BOOST_ASSERT_MSG(!(a_is_point && b_is_point), "both segments shouldn't be degenerated");
+
+ // for degenerated segments the second is always true because this function
+ // shouldn't be called if both segments were degenerated
+
+ if (a_is_point)
+ {
+ return std::make_pair(abs_robust_dx_b >= abs_robust_dy_b, true);
+ }
+ else if (b_is_point)
+ {
+ return std::make_pair(abs_robust_dx_a >= abs_robust_dy_a, true);
+ }
+ else
+ {
+ RobustCoordinateType const min_dx = (std::min)(abs_robust_dx_a, abs_robust_dx_b);
+ RobustCoordinateType const min_dy = (std::min)(abs_robust_dy_a, abs_robust_dy_b);
+ return min_dx == min_dy ?
+ std::make_pair(true, min_dx > RobustCoordinateType(0)) :
+ std::make_pair(min_dx > min_dy, true);
+ }
+ }
+
template
<
std::size_t Dimension,
@@ -319,17 +369,58 @@ private:
RobustType const length_a = oa_2 - oa_1; // no abs, see above
RobustType const length_b = ob_2 - ob_1;
- RatioType const ra_from(oa_1 - ob_1, length_b);
- RatioType const ra_to(oa_2 - ob_1, length_b);
- RatioType const rb_from(ob_1 - oa_1, length_a);
- RatioType const rb_to(ob_2 - oa_1, length_a);
+ RatioType ra_from(oa_1 - ob_1, length_b);
+ RatioType ra_to(oa_2 - ob_1, length_b);
+ RatioType rb_from(ob_1 - oa_1, length_a);
+ RatioType rb_to(ob_2 - oa_1, length_a);
+
+ // use absolute measure to detect endpoints intersection
+ // NOTE: it'd be possible to calculate bx_wrt_a using ax_wrt_b values
+ int const a1_wrt_b = position_value(oa_1, ob_1, ob_2);
+ int const a2_wrt_b = position_value(oa_2, ob_1, ob_2);
+ int const b1_wrt_a = position_value(ob_1, oa_1, oa_2);
+ int const b2_wrt_a = position_value(ob_2, oa_1, oa_2);
+
+ // fix the ratios if necessary
+ // CONSIDER: fixing ratios also in other cases, if they're inconsistent
+ // e.g. if ratio == 1 or 0 (so IP at the endpoint)
+ // but position value indicates that the IP is in the middle of the segment
+ // because one of the segments is very long
+ // In such case the ratios could be moved into the middle direction
+ // by some small value (e.g. EPS+1ULP)
+ if (a1_wrt_b == 1)
+ {
+ ra_from.assign(0, 1);
+ rb_from.assign(0, 1);
+ }
+ else if (a1_wrt_b == 3)
+ {
+ ra_from.assign(1, 1);
+ rb_to.assign(0, 1);
+ }
+
+ if (a2_wrt_b == 1)
+ {
+ ra_to.assign(0, 1);
+ rb_from.assign(1, 1);
+ }
+ else if (a2_wrt_b == 3)
+ {
+ ra_to.assign(1, 1);
+ rb_to.assign(1, 1);
+ }
- if ((ra_from.left() && ra_to.left()) || (ra_from.right() && ra_to.right()))
+ if ((a1_wrt_b < 1 && a2_wrt_b < 1) || (a1_wrt_b > 3 && a2_wrt_b > 3))
+ //if ((ra_from.left() && ra_to.left()) || (ra_from.right() && ra_to.right()))
{
return Policy::disjoint();
}
- return Policy::segments_collinear(a, b, ra_from, ra_to, rb_from, rb_to);
+ bool const opposite = math::sign(length_a) != math::sign(length_b);
+
+ return Policy::segments_collinear(a, b, opposite,
+ a1_wrt_b, a2_wrt_b, b1_wrt_a, b2_wrt_a,
+ ra_from, ra_to, rb_from, rb_to);
}
/// Relate segments where one is degenerate
@@ -351,8 +442,32 @@ private:
// b1/b2 (4..4)
// Ratio: (4-2)/(6-2)
RatioType const ratio(d - s1, s2 - s1);
+
+ if (!ratio.on_segment())
+ {
+ return Policy::disjoint();
+ }
+
return Policy::one_degenerate(degenerate_segment, ratio, a_degenerate);
}
+
+ template <typename ProjCoord1, typename ProjCoord2>
+ static inline int position_value(ProjCoord1 const& ca1,
+ ProjCoord2 const& cb1,
+ ProjCoord2 const& cb2)
+ {
+ // S1x 0 1 2 3 4
+ // S2 |---------->
+ return math::equals(ca1, cb1) ? 1
+ : math::equals(ca1, cb2) ? 3
+ : cb1 < cb2 ?
+ ( ca1 < cb1 ? 0
+ : ca1 > cb2 ? 4
+ : 2 )
+ : ( ca1 > cb1 ? 0
+ : ca1 < cb2 ? 4
+ : 2 );
+ }
};
diff --git a/boost/geometry/strategies/cartesian/centroid_average.hpp b/boost/geometry/strategies/cartesian/centroid_average.hpp
index 76e2f71..c12f6e2 100644
--- a/boost/geometry/strategies/cartesian/centroid_average.hpp
+++ b/boost/geometry/strategies/cartesian/centroid_average.hpp
@@ -4,6 +4,11 @@
// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
+// This file was modified by Oracle on 2015.
+// Modifications copyright (c) 2015 Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
@@ -15,6 +20,8 @@
#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_AVERAGE_HPP
+#include <cstddef>
+
#include <boost/geometry/algorithms/assign.hpp>
#include <boost/geometry/arithmetic/arithmetic.hpp>
#include <boost/geometry/core/coordinate_type.hpp>
@@ -46,7 +53,7 @@ private :
class sum
{
friend class average;
- int count;
+ std::size_t count;
PointCentroid centroid;
public :
@@ -68,10 +75,15 @@ public :
state.count++;
}
- static inline void result(sum const& state, PointCentroid& centroid)
+ static inline bool result(sum const& state, PointCentroid& centroid)
{
centroid = state.centroid;
- divide_value(centroid, state.count);
+ if ( state.count > 0 )
+ {
+ divide_value(centroid, state.count);
+ return true;
+ }
+ return false;
}
};
diff --git a/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp b/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp
index f199fb8..0357f17 100644
--- a/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp
+++ b/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp
@@ -4,6 +4,11 @@
// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
+// This file was modified by Oracle on 2015.
+// Modifications copyright (c) 2015 Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
@@ -15,6 +20,8 @@
#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_BASHEIN_DETMER_HPP
+#include <cstddef>
+
#include <boost/mpl/if.hpp>
#include <boost/numeric/conversion/cast.hpp>
#include <boost/type_traits.hpp>
@@ -143,7 +150,7 @@ private :
class sums
{
friend class bashein_detmer;
- int count;
+ std::size_t count;
calculation_type sum_a2;
calculation_type sum_x;
calculation_type sum_y;
diff --git a/boost/geometry/strategies/cartesian/side_by_triangle.hpp b/boost/geometry/strategies/cartesian/side_by_triangle.hpp
index 5d589ff..77443d4 100644
--- a/boost/geometry/strategies/cartesian/side_by_triangle.hpp
+++ b/boost/geometry/strategies/cartesian/side_by_triangle.hpp
@@ -1,8 +1,13 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
-// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
-// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
-// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
+// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
+// Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
+
+// This file was modified by Oracle on 2015.
+// Modifications copyright (c) 2015, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
@@ -22,6 +27,9 @@
#include <boost/geometry/util/select_coordinate_type.hpp>
#include <boost/geometry/strategies/side.hpp>
+#include <boost/geometry/algorithms/detail/relate/less.hpp>
+#include <boost/geometry/algorithms/detail/equals/point_point.hpp>
+
namespace boost { namespace geometry
{
@@ -38,6 +46,24 @@ namespace strategy { namespace side
template <typename CalculationType = void>
class side_by_triangle
{
+ template <typename Policy>
+ struct eps_policy
+ {
+ eps_policy() {}
+ template <typename Type>
+ eps_policy(Type const& a, Type const& b, Type const& c, Type const& d)
+ : policy(a, b, c, d)
+ {}
+ Policy policy;
+ };
+
+ struct eps_empty
+ {
+ eps_empty() {}
+ template <typename Type>
+ eps_empty(Type const&, Type const&, Type const&, Type const&) {}
+ };
+
public :
// Template member function, because it is not always trivial
@@ -47,23 +73,34 @@ public :
// Types can be all three different. Therefore it is
// not implemented (anymore) as "segment"
- template <typename coordinate_type, typename promoted_type, typename P1, typename P2, typename P>
- static inline promoted_type side_value(P1 const& p1, P2 const& p2, P const& p)
+ template
+ <
+ typename CoordinateType,
+ typename PromotedType,
+ typename P1,
+ typename P2,
+ typename P,
+ typename EpsPolicy
+ >
+ static inline
+ PromotedType side_value(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & eps_policy)
{
- coordinate_type const x = get<0>(p);
- coordinate_type const y = get<1>(p);
+ CoordinateType const x = get<0>(p);
+ CoordinateType const y = get<1>(p);
- coordinate_type const sx1 = get<0>(p1);
- coordinate_type const sy1 = get<1>(p1);
- coordinate_type const sx2 = get<0>(p2);
- coordinate_type const sy2 = get<1>(p2);
+ CoordinateType const sx1 = get<0>(p1);
+ CoordinateType const sy1 = get<1>(p1);
+ CoordinateType const sx2 = get<0>(p2);
+ CoordinateType const sy2 = get<1>(p2);
- promoted_type const dx = sx2 - sx1;
- promoted_type const dy = sy2 - sy1;
- promoted_type const dpx = x - sx1;
- promoted_type const dpy = y - sy1;
+ PromotedType const dx = sx2 - sx1;
+ PromotedType const dy = sy2 - sy1;
+ PromotedType const dpx = x - sx1;
+ PromotedType const dpy = y - sy1;
- return geometry::detail::determinant<promoted_type>
+ eps_policy = EpsPolicy(dx, dy, dpx, dpy);
+
+ return geometry::detail::determinant<PromotedType>
(
dx, dy,
dpx, dpy
@@ -71,9 +108,99 @@ public :
}
+ template
+ <
+ typename CoordinateType,
+ typename PromotedType,
+ typename P1,
+ typename P2,
+ typename P
+ >
+ static inline
+ PromotedType side_value(P1 const& p1, P2 const& p2, P const& p)
+ {
+ eps_empty dummy;
+ return side_value<CoordinateType, PromotedType>(p1, p2, p, dummy);
+ }
+
+
+ template
+ <
+ typename CoordinateType,
+ typename PromotedType,
+ bool AreAllIntegralCoordinates
+ >
+ struct compute_side_value
+ {
+ template <typename P1, typename P2, typename P, typename EpsPolicy>
+ static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
+ {
+ return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
+ }
+ };
+
+ template <typename CoordinateType, typename PromotedType>
+ struct compute_side_value<CoordinateType, PromotedType, false>
+ {
+ template <typename P1, typename P2, typename P, typename EpsPolicy>
+ static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
+ {
+ // For robustness purposes, first check if any two points are
+ // the same; in this case simply return that the points are
+ // collinear
+ if (geometry::detail::equals::equals_point_point(p1, p2)
+ || geometry::detail::equals::equals_point_point(p1, p)
+ || geometry::detail::equals::equals_point_point(p2, p))
+ {
+ return PromotedType(0);
+ }
+
+ // The side_by_triangle strategy computes the signed area of
+ // the point triplet (p1, p2, p); as such it is (in theory)
+ // invariant under cyclic permutations of its three arguments.
+ //
+ // In the context of numerical errors that arise in
+ // floating-point computations, and in order to make the strategy
+ // consistent with respect to cyclic permutations of its three
+ // arguments, we cyclically permute them so that the first
+ // argument is always the lexicographically smallest point.
+
+ geometry::detail::relate::less less;
+ if (less(p, p1))
+ {
+ if (less(p, p2))
+ {
+ // p is the lexicographically smallest
+ return side_value<CoordinateType, PromotedType>(p, p1, p2, epsp);
+ }
+ else
+ {
+ // p2 is the lexicographically smallest
+ return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
+ }
+ }
+
+ if (less(p1, p2))
+ {
+ // p1 is the lexicographically smallest
+ return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
+ }
+ else
+ {
+ // p2 is the lexicographically smallest
+ return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
+ }
+ }
+ };
+
+
template <typename P1, typename P2, typename P>
static inline int apply(P1 const& p1, P2 const& p2, P const& p)
{
+ typedef typename coordinate_type<P1>::type coordinate_type1;
+ typedef typename coordinate_type<P2>::type coordinate_type2;
+ typedef typename coordinate_type<P>::type coordinate_type3;
+
typedef typename boost::mpl::if_c
<
boost::is_void<CalculationType>::type::value,
@@ -81,10 +208,9 @@ public :
<
typename select_most_precise
<
- typename coordinate_type<P1>::type,
- typename coordinate_type<P2>::type
+ coordinate_type1, coordinate_type2
>::type,
- typename coordinate_type<P>::type
+ coordinate_type3
>::type,
CalculationType
>::type coordinate_type;
@@ -96,10 +222,19 @@ public :
double
>::type promoted_type;
- promoted_type const s = side_value<coordinate_type, promoted_type>(p1, p2, p);
- promoted_type const zero = promoted_type();
+ bool const are_all_integral_coordinates =
+ boost::is_integral<coordinate_type1>::value
+ && boost::is_integral<coordinate_type2>::value
+ && boost::is_integral<coordinate_type3>::value;
- return math::equals(s, zero) ? 0
+ eps_policy< math::detail::equals_factor_policy<promoted_type> > epsp;
+ promoted_type s = compute_side_value
+ <
+ coordinate_type, promoted_type, are_all_integral_coordinates
+ >::apply(p1, p2, p, epsp);
+
+ promoted_type const zero = promoted_type();
+ return math::detail::equals_by_policy(s, zero, epsp.policy) ? 0
: s > zero ? 1
: -1;
}
diff --git a/boost/geometry/strategies/cartesian/side_of_intersection.hpp b/boost/geometry/strategies/cartesian/side_of_intersection.hpp
new file mode 100644
index 0000000..3948767
--- /dev/null
+++ b/boost/geometry/strategies/cartesian/side_of_intersection.hpp
@@ -0,0 +1,119 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
+#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
+
+
+#include <boost/geometry/arithmetic/determinant.hpp>
+#include <boost/geometry/core/access.hpp>
+#include <boost/geometry/core/coordinate_type.hpp>
+#include <boost/geometry/util/math.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+namespace strategy { namespace side
+{
+
+// Calculates the side of the intersection-point (if any) of
+// of segment a//b w.r.t. segment c
+// This is calculated without (re)calculating the IP itself again and fully
+// based on integer mathematics; there are no divisions
+// It can be used for either integer (rescaled) points, and also for FP
+class side_of_intersection
+{
+public :
+
+ // Calculates the side of the intersection-point (if any) of
+ // of segment a//b w.r.t. segment c
+ // This is calculated without (re)calculating the IP itself again and fully
+ // based on integer mathematics
+ template <typename T, typename Segment>
+ static inline T side_value(Segment const& a, Segment const& b,
+ Segment const& c)
+ {
+ // The first point of the three segments is reused several times
+ T const ax = get<0, 0>(a);
+ T const ay = get<0, 1>(a);
+ T const bx = get<0, 0>(b);
+ T const by = get<0, 1>(b);
+ T const cx = get<0, 0>(c);
+ T const cy = get<0, 1>(c);
+
+ T const dx_a = get<1, 0>(a) - ax;
+ T const dy_a = get<1, 1>(a) - ay;
+
+ T const dx_b = get<1, 0>(b) - bx;
+ T const dy_b = get<1, 1>(b) - by;
+
+ T const dx_c = get<1, 0>(c) - cx;
+ T const dy_c = get<1, 1>(c) - cy;
+
+ // Cramer's rule: d (see cart_intersect.hpp)
+ T const d = geometry::detail::determinant<T>
+ (
+ dx_a, dy_a,
+ dx_b, dy_b
+ );
+
+ T const zero = T();
+ if (d == zero)
+ {
+ // There is no IP of a//b, they are collinear or parallel
+ // We don't have to divide but we can already conclude the side-value
+ // is meaningless and the resulting determinant will be 0
+ return zero;
+ }
+
+ // Cramer's rule: da (see cart_intersect.hpp)
+ T const da = geometry::detail::determinant<T>
+ (
+ dx_b, dy_b,
+ ax - bx, ay - by
+ );
+
+ // IP is at (ax + (da/d) * dx_a, ay + (da/d) * dy_a)
+ // Side of IP is w.r.t. c is: determinant(dx_c, dy_c, ipx-cx, ipy-cy)
+ // We replace ipx by expression above and multiply each term by d
+ T const result = geometry::detail::determinant<T>
+ (
+ dx_c * d, dy_c * d,
+ d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da
+ );
+
+ // Note: result / (d * d)
+ // is identical to the side_value of side_by_triangle
+ // Therefore, the sign is always the same as that result, and the
+ // resulting side (left,right,collinear) is the same
+
+ return result;
+
+ }
+
+ template <typename Segment>
+ static inline int apply(Segment const& a, Segment const& b, Segment const& c)
+ {
+ typedef typename geometry::coordinate_type<Segment>::type coordinate_type;
+ coordinate_type const s = side_value<coordinate_type>(a, b, c);
+ coordinate_type const zero = coordinate_type();
+ return math::equals(s, zero) ? 0
+ : s > zero ? 1
+ : -1;
+ }
+
+};
+
+
+}} // namespace strategy::side
+
+}} // namespace boost::geometry
+
+
+#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP
diff --git a/boost/geometry/strategies/comparable_distance_result.hpp b/boost/geometry/strategies/comparable_distance_result.hpp
index a258ddb..5ba9b16 100644
--- a/boost/geometry/strategies/comparable_distance_result.hpp
+++ b/boost/geometry/strategies/comparable_distance_result.hpp
@@ -1,6 +1,6 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
-// Copyright (c) 2014, Oracle and/or its affiliates.
+// Copyright (c) 2014-2015, Oracle and/or its affiliates.
// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
@@ -91,9 +91,9 @@ struct comparable_distance_result
// A set of all variant type combinations that are compatible and
// implemented
typedef typename util::combine_if<
- typename mpl::vector1<Geometry1>,
+ typename boost::mpl::vector1<Geometry1>,
typename boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)>::types,
- mpl::always<mpl::true_>
+ boost::mpl::always<boost::mpl::true_>
>::type possible_input_types;
// The (possibly variant) result type resulting from these combinations
@@ -101,11 +101,11 @@ struct comparable_distance_result
typename transform_variant<
possible_input_types,
resolve_strategy::comparable_distance_result<
- mpl::first<mpl::_>,
- mpl::second<mpl::_>,
+ boost::mpl::first<boost::mpl::_>,
+ boost::mpl::second<boost::mpl::_>,
Strategy
>,
- mpl::back_inserter<mpl::vector0<> >
+ boost::mpl::back_inserter<boost::mpl::vector0<> >
>::type
>::type type;
};
@@ -144,7 +144,7 @@ struct comparable_distance_result
<
typename boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)>::types,
typename boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)>::types,
- mpl::always<mpl::true_>
+ boost::mpl::always<boost::mpl::true_>
>::type possible_input_types;
// The (possibly variant) result type resulting from these combinations
@@ -152,11 +152,11 @@ struct comparable_distance_result
typename transform_variant<
possible_input_types,
resolve_strategy::comparable_distance_result<
- mpl::first<mpl::_>,
- mpl::second<mpl::_>,
+ boost::mpl::first<boost::mpl::_>,
+ boost::mpl::second<boost::mpl::_>,
Strategy
>,
- mpl::back_inserter<mpl::vector0<> >
+ boost::mpl::back_inserter<boost::mpl::vector0<> >
>::type
>::type type;
};
diff --git a/boost/geometry/strategies/concepts/distance_concept.hpp b/boost/geometry/strategies/concepts/distance_concept.hpp
index a0cbbd2..6e75fa9 100644
--- a/boost/geometry/strategies/concepts/distance_concept.hpp
+++ b/boost/geometry/strategies/concepts/distance_concept.hpp
@@ -1,8 +1,13 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
-// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
-// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
-// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
+// Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2008-2014 Bruno Lalande, Paris, France.
+// Copyright (c) 2009-2014 Mateusz Loskot, London, UK.
+
+// 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
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
@@ -19,6 +24,8 @@
#include <boost/concept_check.hpp>
#include <boost/core/ignore_unused.hpp>
+#include <boost/mpl/assert.hpp>
+#include <boost/type_traits/is_same.hpp>
#include <boost/geometry/util/parameter_type_of.hpp>
@@ -26,13 +33,15 @@
#include <boost/geometry/geometries/segment.hpp>
#include <boost/geometry/geometries/point.hpp>
+#include <boost/geometry/strategies/tags.hpp>
+
namespace boost { namespace geometry { namespace concept
{
/*!
- \brief Checks strategy for point-segment-distance
+ \brief Checks strategy for point-point or point-box or box-box distance
\ingroup distance
*/
template <typename Strategy, typename Point1, typename Point2>
@@ -57,7 +66,7 @@ private :
ApplyMethod, 1
>::type ptype2;
- // 2) must define meta-function return_type
+ // 2) must define meta-function "return_type"
typedef typename strategy::distance::services::return_type
<
Strategy, ptype1, ptype2
@@ -75,6 +84,16 @@ private :
Strategy
>::type tag;
+ static const bool is_correct_strategy_tag =
+ boost::is_same<tag, strategy_tag_distance_point_point>::value
+ || boost::is_same<tag, strategy_tag_distance_point_box>::value
+ || boost::is_same<tag, strategy_tag_distance_box_box>::value;
+
+ BOOST_MPL_ASSERT_MSG
+ ((is_correct_strategy_tag),
+ INCORRECT_STRATEGY_TAG,
+ (types<tag>));
+
// 5) must implement apply with arguments
Strategy* str = 0;
ptype1 *p1 = 0;
@@ -111,7 +130,7 @@ public :
/*!
- \brief Checks strategy for point-segment-distance
+ \brief Checks strategy for point-segment distance
\ingroup strategy_concepts
*/
template <typename Strategy, typename Point, typename PointOfSegment>
@@ -125,6 +144,7 @@ private :
template <typename ApplyMethod>
static void apply(ApplyMethod)
{
+ // 1) inspect and define both arguments of apply
typedef typename parameter_type_of
<
ApplyMethod, 0
@@ -135,10 +155,28 @@ private :
ApplyMethod, 1
>::type sptype;
- // must define meta-function return_type
- typedef typename strategy::distance::services::return_type<Strategy, ptype, sptype>::type rtype;
+ namespace services = strategy::distance::services;
+ // 2) must define meta-function "tag"
+ typedef typename services::tag<Strategy>::type tag;
+ BOOST_MPL_ASSERT_MSG
+ ((boost::is_same
+ <
+ tag, strategy_tag_distance_point_segment
+ >::value),
+ INCORRECT_STRATEGY_TAG,
+ (types<tag>));
+ // 3) must define meta-function "return_type"
+ typedef typename services::return_type
+ <
+ Strategy, ptype, sptype
+ >::type rtype;
+
+ // 4) must define meta-function "comparable_type"
+ typedef typename services::comparable_type<Strategy>::type ctype;
+
+ // 5) must implement apply with arguments
Strategy *str = 0;
ptype *p = 0;
sptype *sp1 = 0;
@@ -146,8 +184,16 @@ private :
rtype r = str->apply(*p, *sp1, *sp2);
- boost::ignore_unused_variable_warning(str);
- boost::ignore_unused_variable_warning(r);
+ // 6) must define (meta-)struct "get_comparable" with apply
+ ctype cstrategy = services::get_comparable<Strategy>::apply(*str);
+
+ // 7) must define (meta-)struct "result_from_distance" with apply
+ r = services::result_from_distance
+ <
+ Strategy, ptype, sptype
+ >::apply(*str, rtype(1.0));
+
+ boost::ignore_unused(str, r, cstrategy);
}
};
diff --git a/boost/geometry/strategies/distance_result.hpp b/boost/geometry/strategies/distance_result.hpp
index 185a511..e4f326d 100644
--- a/boost/geometry/strategies/distance_result.hpp
+++ b/boost/geometry/strategies/distance_result.hpp
@@ -1,13 +1,13 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
-// Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands.
-// Copyright (c) 2008-2014 Bruno Lalande, Paris, France.
-// Copyright (c) 2009-2014 Mateusz Loskot, London, UK.
-// Copyright (c) 2013-2014 Adam Wulkiewicz, Lodz, Poland.
-// Copyright (c) 2014 Samuel Debionne, Grenoble, France.
+// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
+// Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
+// Copyright (c) 2013-2015 Adam Wulkiewicz, Lodz, Poland.
+// Copyright (c) 2014-2015 Samuel Debionne, Grenoble, France.
-// This file was modified by Oracle on 2014.
-// Modifications copyright (c) 2014, Oracle and/or its affiliates.
+// This file was modified by Oracle on 2014, 2015.
+// Modifications copyright (c) 2014-2015, Oracle and/or its affiliates.
// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
@@ -100,14 +100,14 @@ struct distance_result
// A set of all variant type combinations that are compatible and
// implemented
typedef typename util::combine_if<
- typename mpl::vector1<Geometry1>,
+ typename boost::mpl::vector1<Geometry1>,
typename boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)>::types,
// Here we want should remove most of the combinations that
// are not valid, mostly to limit the size of the resulting MPL set.
// But is_implementedn is not ready for prime time
//
- // util::is_implemented2<mpl::_1, mpl::_2, dispatch::distance<mpl::_1, mpl::_2> >
- mpl::always<mpl::true_>
+ // util::is_implemented2<boost::mpl::_1, boost::mpl::_2, dispatch::distance<boost::mpl::_1, boost::mpl::_2> >
+ boost::mpl::always<boost::mpl::true_>
>::type possible_input_types;
// The (possibly variant) result type resulting from these combinations
@@ -115,11 +115,11 @@ struct distance_result
typename transform_variant<
possible_input_types,
resolve_strategy::distance_result<
- mpl::first<mpl::_>,
- mpl::second<mpl::_>,
+ boost::mpl::first<boost::mpl::_>,
+ boost::mpl::second<boost::mpl::_>,
Strategy
>,
- mpl::back_inserter<mpl::vector0<> >
+ boost::mpl::back_inserter<boost::mpl::vector0<> >
>::type
>::type type;
};
@@ -163,8 +163,8 @@ struct distance_result
// resulting MPL vector.
// But is_implemented is not ready for prime time
//
- // util::is_implemented2<mpl::_1, mpl::_2, dispatch::distance<mpl::_1, mpl::_2> >
- mpl::always<mpl::true_>
+ // util::is_implemented2<boost::mpl::_1, boost::mpl::_2, dispatch::distance<boost::mpl::_1, boost::mpl::_2> >
+ boost::mpl::always<boost::mpl::true_>
>::type possible_input_types;
// The (possibly variant) result type resulting from these combinations
@@ -172,11 +172,11 @@ struct distance_result
typename transform_variant<
possible_input_types,
resolve_strategy::distance_result<
- mpl::first<mpl::_>,
- mpl::second<mpl::_>,
+ boost::mpl::first<boost::mpl::_>,
+ boost::mpl::second<boost::mpl::_>,
Strategy
>,
- mpl::back_inserter<mpl::vector0<> >
+ boost::mpl::back_inserter<boost::mpl::vector0<> >
>::type
>::type type;
};
diff --git a/boost/geometry/strategies/geographic/distance_andoyer.hpp b/boost/geometry/strategies/geographic/distance_andoyer.hpp
new file mode 100644
index 0000000..64de8c1
--- /dev/null
+++ b/boost/geometry/strategies/geographic/distance_andoyer.hpp
@@ -0,0 +1,224 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2007-2012 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 Adam Wulkiewicz, 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
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ANDOYER_HPP
+#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ANDOYER_HPP
+
+
+#include <boost/geometry/core/coordinate_type.hpp>
+#include <boost/geometry/core/radian_access.hpp>
+#include <boost/geometry/core/radius.hpp>
+#include <boost/geometry/core/srs.hpp>
+
+#include <boost/geometry/algorithms/detail/flattening.hpp>
+
+#include <boost/geometry/strategies/distance.hpp>
+
+#include <boost/geometry/util/math.hpp>
+#include <boost/geometry/util/promote_floating_point.hpp>
+#include <boost/geometry/util/select_calculation_type.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+namespace strategy { namespace distance
+{
+
+
+/*!
+\brief Point-point distance approximation taking flattening into account
+\ingroup distance
+\tparam Spheroid The reference spheroid model
+\tparam CalculationType \tparam_calculation
+\author After Andoyer, 19xx, republished 1950, republished by Meeus, 1999
+\note Although not so well-known, the approximation is very good: in all cases the results
+are about the same as Vincenty. In my (Barend's) testcases the results didn't differ more than 6 m
+\see http://nacc.upc.es/tierra/node16.html
+\see http://sci.tech-archive.net/Archive/sci.geo.satellite-nav/2004-12/2724.html
+\see http://home.att.net/~srschmitt/great_circle_route.html (implementation)
+\see http://www.codeguru.com/Cpp/Cpp/algorithms/article.php/c5115 (implementation)
+\see http://futureboy.homeip.net/frinksamp/navigation.frink (implementation)
+\see http://www.voidware.com/earthdist.htm (implementation)
+*/
+template
+<
+ typename Spheroid,
+ typename CalculationType = void
+>
+class andoyer
+{
+public :
+ template <typename Point1, typename Point2>
+ struct calculation_type
+ : promote_floating_point
+ <
+ typename select_calculation_type
+ <
+ Point1,
+ Point2,
+ CalculationType
+ >::type
+ >
+ {};
+
+ typedef Spheroid model_type;
+
+ inline andoyer()
+ : m_spheroid()
+ {}
+
+ explicit inline andoyer(Spheroid const& spheroid)
+ : m_spheroid(spheroid)
+ {}
+
+
+ template <typename Point1, typename Point2>
+ inline typename calculation_type<Point1, Point2>::type
+ apply(Point1 const& point1, Point2 const& point2) const
+ {
+ return calc<typename calculation_type<Point1, Point2>::type>
+ (
+ get_as_radian<0>(point1), get_as_radian<1>(point1),
+ get_as_radian<0>(point2), get_as_radian<1>(point2)
+ );
+ }
+
+ inline Spheroid const& model() const
+ {
+ return m_spheroid;
+ }
+
+private :
+ template <typename CT, typename T>
+ inline CT calc(T const& lon1,
+ T const& lat1,
+ T const& lon2,
+ T const& lat2) const
+ {
+ CT const G = (lat1 - lat2) / 2.0;
+ CT const lambda = (lon1 - lon2) / 2.0;
+
+ if (geometry::math::equals(lambda, 0.0)
+ && geometry::math::equals(G, 0.0))
+ {
+ return 0.0;
+ }
+
+ CT const F = (lat1 + lat2) / 2.0;
+
+ CT const sinG2 = math::sqr(sin(G));
+ CT const cosG2 = math::sqr(cos(G));
+ CT const sinF2 = math::sqr(sin(F));
+ CT const cosF2 = math::sqr(cos(F));
+ CT const sinL2 = math::sqr(sin(lambda));
+ CT const cosL2 = math::sqr(cos(lambda));
+
+ CT const S = sinG2 * cosL2 + cosF2 * sinL2;
+ CT const C = cosG2 * cosL2 + sinF2 * sinL2;
+
+ CT const c0 = 0;
+ CT const c1 = 1;
+ CT const c2 = 2;
+ CT const c3 = 3;
+
+ if (geometry::math::equals(S, c0) || geometry::math::equals(C, c0))
+ {
+ return c0;
+ }
+
+ CT const radius_a = CT(get_radius<0>(m_spheroid));
+ CT const flattening = geometry::detail::flattening<CT>(m_spheroid);
+
+ CT const omega = atan(math::sqrt(S / C));
+ CT const r3 = c3 * math::sqrt(S * C) / omega; // not sure if this is r or greek nu
+ CT const D = c2 * omega * radius_a;
+ CT const H1 = (r3 - c1) / (c2 * C);
+ CT const H2 = (r3 + c1) / (c2 * S);
+
+ return D * (c1 + flattening * (H1 * sinF2 * cosG2 - H2 * cosF2 * sinG2) );
+ }
+
+ Spheroid m_spheroid;
+};
+
+
+#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
+namespace services
+{
+
+template <typename Spheroid, typename CalculationType>
+struct tag<andoyer<Spheroid, CalculationType> >
+{
+ typedef strategy_tag_distance_point_point type;
+};
+
+
+template <typename Spheroid, typename CalculationType, typename P1, typename P2>
+struct return_type<andoyer<Spheroid, CalculationType>, P1, P2>
+ : andoyer<Spheroid, CalculationType>::template calculation_type<P1, P2>
+{};
+
+
+template <typename Spheroid, typename CalculationType>
+struct comparable_type<andoyer<Spheroid, CalculationType> >
+{
+ typedef andoyer<Spheroid, CalculationType> type;
+};
+
+
+template <typename Spheroid, typename CalculationType>
+struct get_comparable<andoyer<Spheroid, CalculationType> >
+{
+ static inline andoyer<Spheroid, CalculationType> apply(andoyer<Spheroid, CalculationType> const& input)
+ {
+ return input;
+ }
+};
+
+template <typename Spheroid, typename CalculationType, typename P1, typename P2>
+struct result_from_distance<andoyer<Spheroid, CalculationType>, P1, P2>
+{
+ template <typename T>
+ static inline typename return_type<andoyer<Spheroid, CalculationType>, P1, P2>::type
+ apply(andoyer<Spheroid, CalculationType> const& , T const& value)
+ {
+ return value;
+ }
+};
+
+
+template <typename Point1, typename Point2>
+struct default_strategy<point_tag, point_tag, Point1, Point2, geographic_tag, geographic_tag>
+{
+ typedef strategy::distance::andoyer
+ <
+ srs::spheroid
+ <
+ typename select_coordinate_type<Point1, Point2>::type
+ >
+ > type;
+};
+
+
+} // namespace services
+#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
+
+
+}} // namespace strategy::distance
+
+
+}} // namespace boost::geometry
+
+
+#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ANDOYER_HPP
diff --git a/boost/geometry/strategies/geographic/distance_vincenty.hpp b/boost/geometry/strategies/geographic/distance_vincenty.hpp
new file mode 100644
index 0000000..56dd14b
--- /dev/null
+++ b/boost/geometry/strategies/geographic/distance_vincenty.hpp
@@ -0,0 +1,161 @@
+// Boost.Geometry
+
+// Copyright (c) 2007-2012 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 Adam Wulkiewicz, 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
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_VINCENTY_HPP
+#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_VINCENTY_HPP
+
+
+#include <boost/geometry/core/coordinate_type.hpp>
+#include <boost/geometry/core/radian_access.hpp>
+
+#include <boost/geometry/strategies/distance.hpp>
+
+#include <boost/geometry/util/promote_floating_point.hpp>
+#include <boost/geometry/util/select_calculation_type.hpp>
+
+#include <boost/geometry/algorithms/detail/vincenty_inverse.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace strategy { namespace distance
+{
+
+/*!
+\brief Distance calculation formulae on latlong coordinates, after Vincenty, 1975
+\ingroup distance
+\tparam Spheroid The reference spheroid model
+\tparam CalculationType \tparam_calculation
+\author See
+ - http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
+ - http://www.icsm.gov.au/gda/gdav2.3.pdf
+\author Adapted from various implementations to get it close to the original document
+ - http://www.movable-type.co.uk/scripts/LatLongVincenty.html
+ - http://exogen.case.edu/projects/geopy/source/geopy.distance.html
+ - http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink
+
+*/
+template
+<
+ typename Spheroid,
+ typename CalculationType = void
+>
+class vincenty
+{
+public :
+ template <typename Point1, typename Point2>
+ struct calculation_type
+ : promote_floating_point
+ <
+ typename select_calculation_type
+ <
+ Point1,
+ Point2,
+ CalculationType
+ >::type
+ >
+ {};
+
+ typedef Spheroid model_type;
+
+ inline vincenty()
+ : m_spheroid()
+ {}
+
+ explicit inline vincenty(Spheroid const& spheroid)
+ : m_spheroid(spheroid)
+ {}
+
+ template <typename Point1, typename Point2>
+ inline typename calculation_type<Point1, Point2>::type
+ apply(Point1 const& point1, Point2 const& point2) const
+ {
+ return geometry::detail::vincenty_inverse
+ <
+ typename calculation_type<Point1, Point2>::type
+ >(get_as_radian<0>(point1),
+ get_as_radian<1>(point1),
+ get_as_radian<0>(point2),
+ get_as_radian<1>(point2),
+ m_spheroid).distance();
+ }
+
+ inline Spheroid const& model() const
+ {
+ return m_spheroid;
+ }
+
+private :
+ Spheroid m_spheroid;
+};
+
+#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
+namespace services
+{
+
+template <typename Spheroid, typename CalculationType>
+struct tag<vincenty<Spheroid, CalculationType> >
+{
+ typedef strategy_tag_distance_point_point type;
+};
+
+
+template <typename Spheroid, typename CalculationType, typename P1, typename P2>
+struct return_type<vincenty<Spheroid, CalculationType>, P1, P2>
+ : vincenty<Spheroid, CalculationType>::template calculation_type<P1, P2>
+{};
+
+
+template <typename Spheroid, typename CalculationType>
+struct comparable_type<vincenty<Spheroid, CalculationType> >
+{
+ typedef vincenty<Spheroid, CalculationType> type;
+};
+
+
+template <typename Spheroid, typename CalculationType>
+struct get_comparable<vincenty<Spheroid, CalculationType> >
+{
+ static inline vincenty<Spheroid, CalculationType> apply(vincenty<Spheroid, CalculationType> const& input)
+ {
+ return input;
+ }
+};
+
+template <typename Spheroid, typename CalculationType, typename P1, typename P2>
+struct result_from_distance<vincenty<Spheroid, CalculationType>, P1, P2 >
+{
+ template <typename T>
+ static inline typename return_type<vincenty<Spheroid, CalculationType>, P1, P2>::type
+ apply(vincenty<Spheroid, CalculationType> const& , T const& value)
+ {
+ return value;
+ }
+};
+
+
+} // namespace services
+#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
+
+
+// We might add a vincenty-like strategy also for point-segment distance, but to calculate the projected point is not trivial
+
+
+
+}} // namespace strategy::distance
+
+
+}} // namespace boost::geometry
+
+
+#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_VINCENTY_HPP
diff --git a/boost/geometry/strategies/geographic/mapping_ssf.hpp b/boost/geometry/strategies/geographic/mapping_ssf.hpp
new file mode 100644
index 0000000..3beedc7
--- /dev/null
+++ b/boost/geometry/strategies/geographic/mapping_ssf.hpp
@@ -0,0 +1,185 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2011-2012 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 Adam Wulkiewicz, 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
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_MAPPING_SSF_HPP
+#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_MAPPING_SSF_HPP
+
+
+#include <boost/core/ignore_unused.hpp>
+
+#include <boost/geometry/core/radius.hpp>
+
+#include <boost/geometry/util/math.hpp>
+#include <boost/geometry/util/promote_floating_point.hpp>
+#include <boost/geometry/util/select_calculation_type.hpp>
+
+#include <boost/geometry/strategies/side.hpp>
+#include <boost/geometry/strategies/spherical/ssf.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+namespace strategy { namespace side
+{
+
+
+// An enumeration type defining types of mapping of geographical
+// latitude to spherical latitude.
+// See: http://en.wikipedia.org/wiki/Great_ellipse
+// http://en.wikipedia.org/wiki/Latitude#Auxiliary_latitudes
+enum mapping_type { mapping_geodetic, mapping_reduced, mapping_geocentric };
+
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace detail
+{
+
+template <typename Spheroid, mapping_type Mapping>
+struct mapper
+{
+ explicit inline mapper(Spheroid const& /*spheroid*/) {}
+
+ template <typename CalculationType>
+ static inline CalculationType const& apply(CalculationType const& lat)
+ {
+ return lat;
+ }
+};
+
+template <typename Spheroid>
+struct mapper<Spheroid, mapping_reduced>
+{
+ typedef typename promote_floating_point
+ <
+ typename radius_type<Spheroid>::type
+ >::type fraction_type;
+
+ explicit inline mapper(Spheroid const& spheroid)
+ {
+ fraction_type const a = geometry::get_radius<0>(spheroid);
+ fraction_type const b = geometry::get_radius<2>(spheroid);
+ b_div_a = b / a;
+ }
+
+ template <typename CalculationType>
+ inline CalculationType apply(CalculationType const& lat) const
+ {
+ return atan(static_cast<CalculationType>(b_div_a) * tan(lat));
+ }
+
+ fraction_type b_div_a;
+};
+
+template <typename Spheroid>
+struct mapper<Spheroid, mapping_geocentric>
+{
+ typedef typename promote_floating_point
+ <
+ typename radius_type<Spheroid>::type
+ >::type fraction_type;
+
+ explicit inline mapper(Spheroid const& spheroid)
+ {
+ fraction_type const a = geometry::get_radius<0>(spheroid);
+ fraction_type const b = geometry::get_radius<2>(spheroid);
+ sqr_b_div_a = b / a;
+ sqr_b_div_a *= sqr_b_div_a;
+ }
+
+ template <typename CalculationType>
+ inline CalculationType apply(CalculationType const& lat) const
+ {
+ return atan(static_cast<CalculationType>(sqr_b_div_a) * tan(lat));
+ }
+
+ fraction_type sqr_b_div_a;
+};
+
+}
+#endif // DOXYGEN_NO_DETAIL
+
+
+/*!
+\brief Check at which side of a geographical segment a point lies
+ left of segment (> 0), right of segment (< 0), on segment (0).
+ The check is performed by mapping the geographical coordinates
+ to spherical coordinates and using spherical_side_formula.
+\ingroup strategies
+\tparam Spheroid The reference spheroid model
+\tparam Mapping The type of mapping of geographical to spherical latitude
+\tparam CalculationType \tparam_calculation
+ */
+template <typename Spheroid,
+ mapping_type Mapping = mapping_geodetic,
+ typename CalculationType = void>
+class mapping_spherical_side_formula
+{
+
+public :
+ inline mapping_spherical_side_formula()
+ : m_mapper(Spheroid())
+ {}
+
+ explicit inline mapping_spherical_side_formula(Spheroid const& spheroid)
+ : m_mapper(spheroid)
+ {}
+
+ template <typename P1, typename P2, typename P>
+ inline int apply(P1 const& p1, P2 const& p2, P const& p)
+ {
+ typedef typename promote_floating_point
+ <
+ typename select_calculation_type_alt
+ <
+ CalculationType,
+ P1, P2, P
+ >::type
+ >::type calculation_type;
+
+ calculation_type lon1 = get_as_radian<0>(p1);
+ calculation_type lat1 = m_mapper.template apply<calculation_type>(get_as_radian<1>(p1));
+ calculation_type lon2 = get_as_radian<0>(p2);
+ calculation_type lat2 = m_mapper.template apply<calculation_type>(get_as_radian<1>(p2));
+ calculation_type lon = get_as_radian<0>(p);
+ calculation_type lat = m_mapper.template apply<calculation_type>(get_as_radian<1>(p));
+
+ return detail::spherical_side_formula(lon1, lat1, lon2, lat2, lon, lat);
+ }
+
+private:
+ side::detail::mapper<Spheroid, Mapping> const m_mapper;
+};
+
+// The specialization for geodetic latitude which can be used directly
+template <typename Spheroid,
+ typename CalculationType>
+class mapping_spherical_side_formula<Spheroid, mapping_geodetic, CalculationType>
+{
+
+public :
+ inline mapping_spherical_side_formula() {}
+ explicit inline mapping_spherical_side_formula(Spheroid const& /*spheroid*/) {}
+
+ template <typename P1, typename P2, typename P>
+ static inline int apply(P1 const& p1, P2 const& p2, P const& p)
+ {
+ return spherical_side_formula<CalculationType>::apply(p1, p2, p);
+ }
+};
+
+}} // namespace strategy::side
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_MAPPING_SSF_HPP
diff --git a/boost/geometry/strategies/spherical/distance_cross_track.hpp b/boost/geometry/strategies/spherical/distance_cross_track.hpp
index a40f03d..486c7ef 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
diff --git a/boost/geometry/strategies/spherical/distance_haversine.hpp b/boost/geometry/strategies/spherical/distance_haversine.hpp
index 8db29c5..8b32056 100644
--- a/boost/geometry/strategies/spherical/distance_haversine.hpp
+++ b/boost/geometry/strategies/spherical/distance_haversine.hpp
@@ -158,7 +158,7 @@ public :
typedef typename calculation_type<Point1, Point2>::type calculation_type;
calculation_type const a = comparable_type::apply(p1, p2);
calculation_type const c = calculation_type(2.0) * asin(math::sqrt(a));
- return m_radius * c;
+ return calculation_type(m_radius) * c;
}
/*!
diff --git a/boost/geometry/strategies/spherical/side_by_cross_track.hpp b/boost/geometry/strategies/spherical/side_by_cross_track.hpp
index c4c5f24..3f7be05 100644
--- a/boost/geometry/strategies/spherical/side_by_cross_track.hpp
+++ b/boost/geometry/strategies/spherical/side_by_cross_track.hpp
@@ -2,6 +2,11 @@
// Copyright (c) 2007-2012 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 Adam Wulkiewicz, 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
// http://www.boost.org/LICENSE_1_0.txt)
@@ -9,18 +14,15 @@
#ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SIDE_BY_CROSS_TRACK_HPP
#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SIDE_BY_CROSS_TRACK_HPP
-#include <boost/mpl/if.hpp>
-#include <boost/type_traits.hpp>
-#include <boost/core/ignore_unused.hpp>
-
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/algorithms/detail/course.hpp>
-#include <boost/geometry/util/select_coordinate_type.hpp>
#include <boost/geometry/util/math.hpp>
+#include <boost/geometry/util/promote_floating_point.hpp>
+#include <boost/geometry/util/select_calculation_type.hpp>
#include <boost/geometry/strategies/side.hpp>
//#include <boost/geometry/strategies/concepts/side_concept.hpp>
@@ -47,27 +49,19 @@ public :
template <typename P1, typename P2, typename P>
static inline int apply(P1 const& p1, P2 const& p2, P const& p)
{
- typedef typename boost::mpl::if_c
+ typedef typename promote_floating_point
<
- boost::is_void<CalculationType>::type::value,
- typename select_most_precise
+ typename select_calculation_type_alt
<
- typename select_most_precise
- <
- typename coordinate_type<P1>::type,
- typename coordinate_type<P2>::type
- >::type,
- typename coordinate_type<P>::type
- >::type,
- CalculationType
- >::type coordinate_type;
-
- boost::ignore_unused<coordinate_type>();
-
- double d1 = 0.001; // m_strategy.apply(sp1, p);
- double crs_AD = geometry::detail::course<double>(p1, p);
- double crs_AB = geometry::detail::course<double>(p1, p2);
- double XTD = asin(sin(d1) * sin(crs_AD - crs_AB));
+ CalculationType,
+ P1, P2, P
+ >::type
+ >::type calc_t;
+
+ calc_t d1 = 0.001; // m_strategy.apply(sp1, p);
+ calc_t crs_AD = geometry::detail::course<calc_t>(p1, p);
+ calc_t crs_AB = geometry::detail::course<calc_t>(p1, p2);
+ calc_t XTD = asin(sin(d1) * sin(crs_AD - crs_AB));
return math::equals(XTD, 0) ? 0 : XTD < 0 ? 1 : -1;
}
diff --git a/boost/geometry/strategies/spherical/ssf.hpp b/boost/geometry/strategies/spherical/ssf.hpp
index 4156295..81f3205 100644
--- a/boost/geometry/strategies/spherical/ssf.hpp
+++ b/boost/geometry/strategies/spherical/ssf.hpp
@@ -16,8 +16,9 @@
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/radian_access.hpp>
-#include <boost/geometry/util/select_coordinate_type.hpp>
#include <boost/geometry/util/math.hpp>
+#include <boost/geometry/util/promote_floating_point.hpp>
+#include <boost/geometry/util/select_calculation_type.hpp>
#include <boost/geometry/strategies/side.hpp>
//#include <boost/geometry/strategies/concepts/side_concept.hpp>
@@ -30,6 +31,43 @@ namespace boost { namespace geometry
namespace strategy { namespace side
{
+#ifndef DOXYGEN_NO_DETAIL
+namespace detail
+{
+
+template <typename T>
+int spherical_side_formula(T const& lambda1, T const& delta1,
+ T const& lambda2, T const& delta2,
+ T const& lambda, T const& delta)
+{
+ // Create temporary points (vectors) on unit a sphere
+ T const cos_delta1 = cos(delta1);
+ T const c1x = cos_delta1 * cos(lambda1);
+ T const c1y = cos_delta1 * sin(lambda1);
+ T const c1z = sin(delta1);
+
+ T const cos_delta2 = cos(delta2);
+ T const c2x = cos_delta2 * cos(lambda2);
+ T const c2y = cos_delta2 * sin(lambda2);
+ T const c2z = sin(delta2);
+
+ // (Third point is converted directly)
+ T const cos_delta = cos(delta);
+
+ // Apply the "Spherical Side Formula" as presented on my blog
+ T const dist
+ = (c1y * c2z - c1z * c2y) * cos_delta * cos(lambda)
+ + (c1z * c2x - c1x * c2z) * cos_delta * sin(lambda)
+ + (c1x * c2y - c1y * c2x) * sin(delta);
+
+ T zero = T();
+ return dist > zero ? 1
+ : dist < zero ? -1
+ : 0;
+}
+
+}
+#endif // DOXYGEN_NO_DETAIL
/*!
\brief Check at which side of a Great Circle segment a point lies
@@ -45,60 +83,25 @@ public :
template <typename P1, typename P2, typename P>
static inline int apply(P1 const& p1, P2 const& p2, P const& p)
{
- typedef typename boost::mpl::if_c
+ typedef typename promote_floating_point
<
- boost::is_void<CalculationType>::type::value,
-
- // Select at least a double...
- typename select_most_precise
+ typename select_calculation_type_alt
<
- typename select_most_precise
- <
- typename select_most_precise
- <
- typename coordinate_type<P1>::type,
- typename coordinate_type<P2>::type
- >::type,
- typename coordinate_type<P>::type
- >::type,
- double
- >::type,
- CalculationType
- >::type coordinate_type;
-
- // Convenient shortcuts
- typedef coordinate_type ct;
- ct const lambda1 = get_as_radian<0>(p1);
- ct const delta1 = get_as_radian<1>(p1);
- ct const lambda2 = get_as_radian<0>(p2);
- ct const delta2 = get_as_radian<1>(p2);
- ct const lambda = get_as_radian<0>(p);
- ct const delta = get_as_radian<1>(p);
-
- // Create temporary points (vectors) on unit a sphere
- ct const cos_delta1 = cos(delta1);
- ct const c1x = cos_delta1 * cos(lambda1);
- ct const c1y = cos_delta1 * sin(lambda1);
- ct const c1z = sin(delta1);
-
- ct const cos_delta2 = cos(delta2);
- ct const c2x = cos_delta2 * cos(lambda2);
- ct const c2y = cos_delta2 * sin(lambda2);
- ct const c2z = sin(delta2);
-
- // (Third point is converted directly)
- ct const cos_delta = cos(delta);
-
- // Apply the "Spherical Side Formula" as presented on my blog
- ct const dist
- = (c1y * c2z - c1z * c2y) * cos_delta * cos(lambda)
- + (c1z * c2x - c1x * c2z) * cos_delta * sin(lambda)
- + (c1x * c2y - c1y * c2x) * sin(delta);
-
- ct zero = ct();
- return dist > zero ? 1
- : dist < zero ? -1
- : 0;
+ CalculationType,
+ P1, P2, P
+ >::type
+ >::type calculation_type;
+
+ calculation_type const lambda1 = get_as_radian<0>(p1);
+ calculation_type const delta1 = get_as_radian<1>(p1);
+ calculation_type const lambda2 = get_as_radian<0>(p2);
+ calculation_type const delta2 = get_as_radian<1>(p2);
+ calculation_type const lambda = get_as_radian<0>(p);
+ calculation_type const delta = get_as_radian<1>(p);
+
+ return detail::spherical_side_formula(lambda1, delta1,
+ lambda2, delta2,
+ lambda, delta);
}
};
diff --git a/boost/geometry/strategies/strategies.hpp b/boost/geometry/strategies/strategies.hpp
index afc5d70..bf436db 100644
--- a/boost/geometry/strategies/strategies.hpp
+++ b/boost/geometry/strategies/strategies.hpp
@@ -63,6 +63,9 @@
#include <boost/geometry/strategies/spherical/compare_circular.hpp>
#include <boost/geometry/strategies/spherical/ssf.hpp>
+#include <boost/geometry/strategies/geographic/distance_andoyer.hpp>
+#include <boost/geometry/strategies/geographic/distance_vincenty.hpp>
+
#include <boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp>
#include <boost/geometry/strategies/agnostic/buffer_distance_asymmetric.hpp>
#include <boost/geometry/strategies/agnostic/hull_graham_andrew.hpp>