diff options
Diffstat (limited to 'boost/geometry/strategies')
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 7b7cd1890f..446d2f02cd 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 bc0c46f644..73bd21ac73 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 a960a6f1f9..1d59e13cf6 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 f4ed7a634f..641533fc6a 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 318047fadb..9e8753251d 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 0000000000..14c69a0597 --- /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 8ad3bbc50d..99e7d9b50f 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 74780d6165..a233f1c4be 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 8fcf3b996c..99ec80527f 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 9e467c85a0..9ec51cd1ec 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 f64a82d8fc..86ebc43c9c 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 66af2d2e9c..a7bd385226 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 76e2f7144c..c12f6e2024 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 f199fb80e5..0357f17e7a 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 5d589ffc86..77443d46a9 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 0000000000..39487676c1 --- /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 a258ddb9b4..5ba9b1603d 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 a0cbbd21ed..6e75fa95a6 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 185a511464..e4f326d3ee 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 0000000000..64de8c1a41 --- /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 0000000000..56dd14bbdc --- /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 0000000000..3beedc7809 --- /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 a40f03dbaf..486c7effbd 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 8db29c5157..8b32056f3e 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 c4c5f24eeb..3f7be05557 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 41562950fd..81f3205e90 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 afc5d7046f..bf436dba2b 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> |