diff options
Diffstat (limited to 'boost/geometry/strategies')
27 files changed, 1153 insertions, 440 deletions
diff --git a/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp b/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp index 1d59e13cf6..0cd0cdc4db 100644 --- a/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp +++ b/boost/geometry/strategies/agnostic/hull_graham_andrew.hpp @@ -24,6 +24,7 @@ #include <boost/range.hpp> +#include <boost/geometry/core/assert.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/core/point_type.hpp> #include <boost/geometry/strategies/convex_hull.hpp> @@ -354,7 +355,7 @@ private: { std::copy(boost::begin(first), boost::end(first), out); - BOOST_ASSERT(closed ? !boost::empty(second) : boost::size(second) > 1); + BOOST_GEOMETRY_ASSERT(closed ? !boost::empty(second) : boost::size(second) > 1); std::copy(++boost::rbegin(second), // skip the first Point closed ? boost::rend(second) : --boost::rend(second), // skip the last Point if open out); diff --git a/boost/geometry/strategies/agnostic/relate.hpp b/boost/geometry/strategies/agnostic/relate.hpp index 9e8753251d..676207694f 100644 --- a/boost/geometry/strategies/agnostic/relate.hpp +++ b/boost/geometry/strategies/agnostic/relate.hpp @@ -1,17 +1,17 @@ // 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 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) -// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle - #ifndef BOOST_GEOMETRY_STRATEGY_AGNOSTIC_RELATE_HPP #define BOOST_GEOMETRY_STRATEGY_AGNOSTIC_RELATE_HPP -#include <boost/geometry/algorithms/detail/relate/relate.hpp> +#include <boost/geometry/algorithms/relate.hpp> namespace boost { namespace geometry @@ -25,7 +25,7 @@ struct relate { static inline bool apply(Geometry1 const& geometry1, Geometry2 const& geometry2) { - return detail::relate::relate<StaticMask>(geometry1, geometry2); + return geometry::relate(geometry1, geometry2, StaticMask()); } }; @@ -47,7 +47,10 @@ struct default_strategy<AnyTag1, AnyTag2, AnyTag1, AnyTag2, AnyCS, AnyCS, Geomet < Geometry1, Geometry2, - detail::relate::static_mask_within + typename detail::de9im::static_mask_within_type + < + Geometry1, Geometry2 + >::type > type; }; @@ -58,7 +61,10 @@ struct default_strategy<AnyTag1, AnyTag2, AnyTag1, areal_tag, AnyCS, AnyCS, Geom < Geometry1, Geometry2, - detail::relate::static_mask_within + typename detail::de9im::static_mask_within_type + < + Geometry1, Geometry2 + >::type > type; }; @@ -84,7 +90,10 @@ struct default_strategy<AnyTag1, AnyTag2, AnyTag1, AnyTag2, AnyCS, AnyCS, Geomet < Geometry1, Geometry2, - detail::relate::static_mask_covered_by + typename detail::de9im::static_mask_covered_by_type + < + Geometry1, Geometry2 + >::type > type; }; @@ -95,7 +104,10 @@ struct default_strategy<AnyTag1, AnyTag2, AnyTag1, areal_tag, AnyCS, AnyCS, Geom < Geometry1, Geometry2, - detail::relate::static_mask_covered_by + typename detail::de9im::static_mask_covered_by_type + < + Geometry1, Geometry2 + >::type > type; }; diff --git a/boost/geometry/strategies/agnostic/side_by_azimuth.hpp b/boost/geometry/strategies/agnostic/side_by_azimuth.hpp deleted file mode 100644 index 14c69a0597..0000000000 --- a/boost/geometry/strategies/agnostic/side_by_azimuth.hpp +++ /dev/null @@ -1,87 +0,0 @@ -// 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/buffer.hpp b/boost/geometry/strategies/buffer.hpp index 7dbe03b4a9..86bdcc37bf 100644 --- a/boost/geometry/strategies/buffer.hpp +++ b/boost/geometry/strategies/buffer.hpp @@ -66,7 +66,8 @@ enum piece_type buffered_round_end, buffered_flat_end, buffered_point, - buffered_concave // always on the inside + buffered_concave, // always on the inside + piece_type_unknown }; @@ -82,6 +83,17 @@ enum join_selector join_spike // collinear, with overlap, next segment goes back }; +/*! +\brief Enumerates types of result codes from buffer strategies +\ingroup enum +*/ +enum result_code +{ + result_normal, + result_error_numerical, + result_no_output +}; + }} // namespace strategy::buffer diff --git a/boost/geometry/strategies/cartesian/buffer_end_round.hpp b/boost/geometry/strategies/cartesian/buffer_end_round.hpp index a233f1c4be..3d7d5bb467 100644 --- a/boost/geometry/strategies/cartesian/buffer_end_round.hpp +++ b/boost/geometry/strategies/cartesian/buffer_end_round.hpp @@ -1,6 +1,11 @@ // 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 @@ -62,8 +67,7 @@ private : DistanceType const& buffer_distance, RangeOut& range_out) const { - PromotedType const two = 2.0; - PromotedType const two_pi = two * geometry::math::pi<PromotedType>(); + PromotedType const two_pi = geometry::math::two_pi<PromotedType>(); std::size_t point_buffer_count = m_points_per_circle; diff --git a/boost/geometry/strategies/cartesian/buffer_join_miter.hpp b/boost/geometry/strategies/cartesian/buffer_join_miter.hpp index 99ec80527f..5358156f6d 100644 --- a/boost/geometry/strategies/cartesian/buffer_join_miter.hpp +++ b/boost/geometry/strategies/cartesian/buffer_join_miter.hpp @@ -9,7 +9,7 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_MITER_HPP #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_MITER_HPP -#include <boost/assert.hpp> +#include <boost/geometry/core/assert.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/policies/compare.hpp> #include <boost/geometry/util/math.hpp> @@ -99,7 +99,7 @@ public: if (distance > max_distance) { - BOOST_ASSERT(distance != 0.0); + BOOST_GEOMETRY_ASSERT(distance != 0.0); promoted_type const proportion = max_distance / distance; set<0>(p, get<0>(vertex) + dx * proportion); diff --git a/boost/geometry/strategies/cartesian/buffer_join_round.hpp b/boost/geometry/strategies/cartesian/buffer_join_round.hpp index 9ec51cd1ec..e1a433e1ee 100644 --- a/boost/geometry/strategies/cartesian/buffer_join_round.hpp +++ b/boost/geometry/strategies/cartesian/buffer_join_round.hpp @@ -2,6 +2,11 @@ // 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) @@ -11,7 +16,6 @@ #include <algorithm> -#include <boost/assert.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/policies/compare.hpp> #include <boost/geometry/strategies/buffer.hpp> @@ -19,6 +23,7 @@ #include <boost/geometry/util/select_most_precise.hpp> #ifdef BOOST_GEOMETRY_DEBUG_BUFFER_WARN +#include <iostream> #include <boost/geometry/io/wkt/wkt.hpp> #endif @@ -76,8 +81,7 @@ private : PromotedType const dx2 = get<0>(perp2) - get<0>(vertex); PromotedType const dy2 = get<1>(perp2) - get<1>(vertex); - PromotedType const two = 2.0; - PromotedType const two_pi = two * geometry::math::pi<PromotedType>(); + PromotedType const two_pi = geometry::math::two_pi<PromotedType>(); PromotedType const angle1 = atan2(dy1, dx1); PromotedType angle2 = atan2(dy2, dx2); @@ -96,14 +100,14 @@ private : // - 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 - int const n = (std::max)(static_cast<int>( - ceil(m_points_per_circle * angle_diff / two_pi)), 1); + std::size_t const n = (std::max)(static_cast<std::size_t>( + ceil(m_points_per_circle * angle_diff / two_pi)), std::size_t(1)); 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) + for (std::size_t i = 0; i < n - 1; i++, a -= diff) { Point p; set<0>(p, get<0>(vertex) + buffer_distance * cos(a)); diff --git a/boost/geometry/strategies/cartesian/buffer_join_round_by_divide.hpp b/boost/geometry/strategies/cartesian/buffer_join_round_by_divide.hpp index 1444c795af..ed3a010cd5 100644 --- a/boost/geometry/strategies/cartesian/buffer_join_round_by_divide.hpp +++ b/boost/geometry/strategies/cartesian/buffer_join_round_by_divide.hpp @@ -9,7 +9,6 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_ROUND_BY_DIVIDE_HPP #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_BUFFER_JOIN_ROUND_BY_DIVIDE_HPP -#include <boost/assert.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/policies/compare.hpp> #include <boost/geometry/strategies/buffer.hpp> diff --git a/boost/geometry/strategies/cartesian/buffer_point_circle.hpp b/boost/geometry/strategies/cartesian/buffer_point_circle.hpp index 86ebc43c9c..f289857177 100644 --- a/boost/geometry/strategies/cartesian/buffer_point_circle.hpp +++ b/boost/geometry/strategies/cartesian/buffer_point_circle.hpp @@ -85,8 +85,7 @@ public : promoted_type const buffer_distance = distance_strategy.apply(point, point, strategy::buffer::buffer_side_left); - promoted_type const two = 2.0; - promoted_type const two_pi = two * geometry::math::pi<promoted_type>(); + promoted_type const two_pi = geometry::math::two_pi<promoted_type>(); promoted_type const diff = two_pi / promoted_type(m_count); promoted_type a = 0; diff --git a/boost/geometry/strategies/cartesian/buffer_side_straight.hpp b/boost/geometry/strategies/cartesian/buffer_side_straight.hpp index 24655ab3d7..75247377fa 100644 --- a/boost/geometry/strategies/cartesian/buffer_side_straight.hpp +++ b/boost/geometry/strategies/cartesian/buffer_side_straight.hpp @@ -9,6 +9,8 @@ #include <cstddef> +#include <boost/math/special_functions/fpclassify.hpp> + #include <boost/geometry/core/coordinate_type.hpp> #include <boost/geometry/core/access.hpp> #include <boost/geometry/util/math.hpp> @@ -51,9 +53,9 @@ public : typename OutputRange, typename DistanceStrategy > - static inline void apply( + static inline result_code apply( Point const& input_p1, Point const& input_p2, - strategy::buffer::buffer_side_selector side, + buffer_side_selector side, DistanceStrategy const& distance, OutputRange& output_range) { @@ -73,19 +75,46 @@ public : // For normalization [0,1] (=dot product d.d, sqrt) promoted_type const length = geometry::math::sqrt(dx * dx + dy * dy); + if (! boost::math::isfinite(length)) + { + // In case of coordinates differences of e.g. 1e300, length + // will overflow and we should not generate output +#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_WARN + std::cout << "Error in length calculation for points " + << geometry::wkt(input_p1) << " " << geometry::wkt(input_p2) + << " length: " << length << std::endl; +#endif + return result_error_numerical; + } + if (geometry::math::equals(length, 0)) { // Coordinates are simplified and therefore most often not equal. // But if simplify is skipped, or for lines with two // equal points, length is 0 and we cannot generate output. - return; + return result_no_output; } + promoted_type const d = distance.apply(input_p1, input_p2, side); + // Generate the normalized perpendicular p, to the left (ccw) promoted_type const px = -dy / length; promoted_type const py = dx / length; - promoted_type const d = distance.apply(input_p1, input_p2, side); + if (geometry::math::equals(px, 0) + && geometry::math::equals(py, 0)) + { + // This basically should not occur - because of the checks above. + // There are no unit tests triggering this condition +#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_WARN + std::cout << "Error in perpendicular calculation for points " + << geometry::wkt(input_p1) << " " << geometry::wkt(input_p2) + << " length: " << length + << " distance: " << d + << std::endl; +#endif + return result_no_output; + } output_range.resize(2); @@ -93,6 +122,8 @@ public : set<1>(output_range.front(), get<1>(input_p1) + py * d); set<0>(output_range.back(), get<0>(input_p2) + px * d); set<1>(output_range.back(), get<1>(input_p2) + py * d); + + return result_normal; } #endif // DOXYGEN_SHOULD_SKIP_THIS }; diff --git a/boost/geometry/strategies/cartesian/cart_intersect.hpp b/boost/geometry/strategies/cartesian/cart_intersect.hpp index a7bd385226..8a9857376a 100644 --- a/boost/geometry/strategies/cartesian/cart_intersect.hpp +++ b/boost/geometry/strategies/cartesian/cart_intersect.hpp @@ -119,15 +119,10 @@ struct relate_cartesian_segments boost::ignore_unused_variable_warning(robust_policy); - typedef typename select_calculation_type - <Segment1, Segment2, CalculationType>::type coordinate_type; - using geometry::detail::equals::equals_point_point; bool const a_is_point = equals_point_point(robust_a1, robust_a2); bool const b_is_point = equals_point_point(robust_b1, robust_b2); - typedef side::side_by_triangle<coordinate_type> side; - if(a_is_point && b_is_point) { return equals_point_point(robust_a1, robust_b2) @@ -136,20 +131,32 @@ struct relate_cartesian_segments ; } + typedef typename select_calculation_type + <Segment1, Segment2, CalculationType>::type coordinate_type; + + typedef side::side_by_triangle<coordinate_type> side; + side_info sides; sides.set<0>(side::apply(robust_b1, robust_b2, robust_a1), side::apply(robust_b1, robust_b2, robust_a2)); - sides.set<1>(side::apply(robust_a1, robust_a2, robust_b1), - side::apply(robust_a1, robust_a2, robust_b2)); - bool collinear = sides.collinear(); + if (sides.same<0>()) + { + // Both points are at same side of other segment, we can leave + return Policy::disjoint(); + } - if (sides.same<0>() || sides.same<1>()) + sides.set<1>(side::apply(robust_a1, robust_a2, robust_b1), + side::apply(robust_a1, robust_a2, robust_b2)); + + if (sides.same<1>()) { // Both points are at same side of other segment, we can leave return Policy::disjoint(); } + bool collinear = sides.collinear(); + typedef typename select_most_precise < coordinate_type, double @@ -264,7 +271,7 @@ private: 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"); + //BOOST_GEOMETRY_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 diff --git a/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp b/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp index 0357f17e7a..47763a9f30 100644 --- a/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp +++ b/boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp @@ -1,8 +1,8 @@ // 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. @@ -22,6 +22,7 @@ #include <cstddef> +#include <boost/math/special_functions/fpclassify.hpp> #include <boost/mpl/if.hpp> #include <boost/numeric/conversion/cast.hpp> #include <boost/type_traits.hpp> @@ -206,11 +207,18 @@ public : Point >::type coordinate_type; - set<0>(centroid, - boost::numeric_cast<coordinate_type>(state.sum_x / a3)); - set<1>(centroid, - boost::numeric_cast<coordinate_type>(state.sum_y / a3)); - return true; + // Prevent NaN centroid coordinates + if (boost::math::isfinite(a3)) + { + // NOTE: above calculation_type is checked, not the centroid coordinate_type + // which means that the centroid can still be filled with INF + // if e.g. calculation_type is double and centroid contains floats + set<0>(centroid, + boost::numeric_cast<coordinate_type>(state.sum_x / a3)); + set<1>(centroid, + boost::numeric_cast<coordinate_type>(state.sum_y / a3)); + return true; + } } return false; diff --git a/boost/geometry/strategies/cartesian/centroid_weighted_length.hpp b/boost/geometry/strategies/cartesian/centroid_weighted_length.hpp index b788738d15..0735e925b6 100644 --- a/boost/geometry/strategies/cartesian/centroid_weighted_length.hpp +++ b/boost/geometry/strategies/cartesian/centroid_weighted_length.hpp @@ -1,7 +1,12 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2009-2012 Mateusz Loskot, London, UK. -// Copyright (c) 2009-2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2009-2015 Mateusz Loskot, London, UK. +// Copyright (c) 2009-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 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. @@ -13,9 +18,13 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_WEIGHTED_LENGTH_HPP #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CENTROID_WEIGHTED_LENGTH_HPP +#include <boost/math/special_functions/fpclassify.hpp> +#include <boost/numeric/conversion/cast.hpp> + #include <boost/geometry/algorithms/detail/distance/interface.hpp> #include <boost/geometry/algorithms/detail/distance/point_to_geometry.hpp> #include <boost/geometry/arithmetic/arithmetic.hpp> +#include <boost/geometry/util/for_each_coordinate.hpp> #include <boost/geometry/util/select_most_precise.hpp> #include <boost/geometry/strategies/centroid.hpp> #include <boost/geometry/strategies/default_distance_result.hpp> @@ -91,17 +100,37 @@ public : static inline bool result(state_type const& state, Point& centroid) { distance_type const zero = distance_type(); - if (! geometry::math::equals(state.length, zero)) + if (! geometry::math::equals(state.length, zero) + && boost::math::isfinite(state.length)) // Prevent NaN centroid coordinates { - assign_zero(centroid); - add_point(centroid, state.average_sum); - divide_value(centroid, state.length); + // NOTE: above distance_type is checked, not the centroid coordinate_type + // which means that the centroid can still be filled with INF + // if e.g. distance_type is double and centroid contains floats + geometry::for_each_coordinate(centroid, set_sum_div_length(state)); return true; } return false; } + struct set_sum_div_length + { + state_type const& m_state; + set_sum_div_length(state_type const& state) + : m_state(state) + {} + template <typename Pt, std::size_t Dimension> + void apply(Pt & centroid) const + { + typedef typename geometry::coordinate_type<Pt>::type coordinate_type; + geometry::set<Dimension>( + centroid, + boost::numeric_cast<coordinate_type>( + geometry::get<Dimension>(m_state.average_sum) / m_state.length + ) + ); + } + }; }; #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS diff --git a/boost/geometry/strategies/cartesian/side_of_intersection.hpp b/boost/geometry/strategies/cartesian/side_of_intersection.hpp index 39487676c1..db57644ad5 100644 --- a/boost/geometry/strategies/cartesian/side_of_intersection.hpp +++ b/boost/geometry/strategies/cartesian/side_of_intersection.hpp @@ -2,6 +2,11 @@ // Copyright (c) 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 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) @@ -10,11 +15,25 @@ #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_OF_INTERSECTION_HPP +#include <limits> + +#include <boost/core/ignore_unused.hpp> +#include <boost/type_traits/is_integral.hpp> +#include <boost/type_traits/make_unsigned.hpp> + #include <boost/geometry/arithmetic/determinant.hpp> #include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/assert.hpp> #include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/algorithms/detail/assign_indexed_point.hpp> +#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp> #include <boost/geometry/util/math.hpp> +#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG +#include <boost/math/common_factor_ct.hpp> +#include <boost/math/common_factor_rt.hpp> +#include <boost/multiprecision/cpp_int.hpp> +#endif namespace boost { namespace geometry { @@ -22,6 +41,94 @@ namespace boost { namespace geometry namespace strategy { namespace side { +namespace detail +{ + +// A tool for multiplication of integers avoiding overflow +// It's a temporary workaround until we can use Multiprecision +// The algorithm is based on Karatsuba algorithm +// see: http://en.wikipedia.org/wiki/Karatsuba_algorithm +template <typename T> +struct multiplicable_integral +{ + // Currently this tool can't be used with non-integral coordinate types. + // Also side_of_intersection strategy sign_of_product() and sign_of_compare() + // functions would have to be modified to properly support floating-point + // types (comparisons and multiplication). + BOOST_STATIC_ASSERT(boost::is_integral<T>::value); + + static const std::size_t bits = CHAR_BIT * sizeof(T); + static const std::size_t half_bits = bits / 2; + typedef typename boost::make_unsigned<T>::type unsigned_type; + static const unsigned_type base = unsigned_type(1) << half_bits; // 2^half_bits + + int m_sign; + unsigned_type m_ms; + unsigned_type m_ls; + + multiplicable_integral(int sign, unsigned_type ms, unsigned_type ls) + : m_sign(sign), m_ms(ms), m_ls(ls) + {} + + explicit multiplicable_integral(T const& val) + { + unsigned_type val_u = val > 0 ? + unsigned_type(val) + : val == (std::numeric_limits<T>::min)() ? + unsigned_type((std::numeric_limits<T>::max)()) + 1 + : unsigned_type(-val); + // MMLL -> S 00MM 00LL + m_sign = math::sign(val); + m_ms = val_u >> half_bits; // val_u / base + m_ls = val_u - m_ms * base; + } + + friend multiplicable_integral operator*(multiplicable_integral const& a, + multiplicable_integral const& b) + { + // (S 00MM 00LL) * (S 00MM 00LL) -> (S Z2MM 00LL) + unsigned_type z2 = a.m_ms * b.m_ms; + unsigned_type z0 = a.m_ls * b.m_ls; + unsigned_type z1 = (a.m_ms + a.m_ls) * (b.m_ms + b.m_ls) - z2 - z0; + // z0 may be >= base so it must be normalized to allow comparison + unsigned_type z0_ms = z0 >> half_bits; // z0 / base + return multiplicable_integral(a.m_sign * b.m_sign, + z2 * base + z1 + z0_ms, + z0 - base * z0_ms); + } + + friend bool operator<(multiplicable_integral const& a, + multiplicable_integral const& b) + { + if ( a.m_sign == b.m_sign ) + { + bool u_less = a.m_ms < b.m_ms + || (a.m_ms == b.m_ms && a.m_ls < b.m_ls); + return a.m_sign > 0 ? u_less : (! u_less); + } + else + { + return a.m_sign < b.m_sign; + } + } + + friend bool operator>(multiplicable_integral const& a, + multiplicable_integral const& b) + { + return b < a; + } + + template <typename CmpVal> + void check_value(CmpVal const& cmp_val) const + { + unsigned_type b = base; // a workaround for MinGW - undefined reference base + CmpVal val = CmpVal(m_sign) * (CmpVal(m_ms) * CmpVal(b) + CmpVal(m_ls)); + BOOST_GEOMETRY_ASSERT(cmp_val == val); + } +}; + +} // namespace detail + // 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 @@ -29,15 +136,93 @@ namespace strategy { namespace side // It can be used for either integer (rescaled) points, and also for FP class side_of_intersection { +private : + template <typename T, typename U> + static inline + int sign_of_product(T const& a, U const& b) + { + return a == 0 || b == 0 ? 0 + : a > 0 && b > 0 ? 1 + : a < 0 && b < 0 ? 1 + : -1; + } + + template <typename T> + static inline + int sign_of_compare(T const& a, T const& b, T const& c, T const& d) + { + // Both a*b and c*d are positive + // We have to judge if a*b > c*d + + using side::detail::multiplicable_integral; + multiplicable_integral<T> ab = multiplicable_integral<T>(a) + * multiplicable_integral<T>(b); + multiplicable_integral<T> cd = multiplicable_integral<T>(c) + * multiplicable_integral<T>(d); + + int result = ab > cd ? 1 + : ab < cd ? -1 + : 0 + ; + +#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG + using namespace boost::multiprecision; + cpp_int const lab = cpp_int(a) * cpp_int(b); + cpp_int const lcd = cpp_int(c) * cpp_int(d); + + ab.check_value(lab); + cd.check_value(lcd); + + int result2 = lab > lcd ? 1 + : lab < lcd ? -1 + : 0 + ; + BOOST_GEOMETRY_ASSERT(result == result2); +#endif + + return result; + } + + template <typename T> + static inline + int sign_of_addition_of_two_products(T const& a, T const& b, T const& c, T const& d) + { + // sign of a*b+c*d, 1 if positive, -1 if negative, else 0 + int const ab = sign_of_product(a, b); + int const cd = sign_of_product(c, d); + if (ab == 0) + { + return cd; + } + if (cd == 0) + { + return ab; + } + + if (ab == cd) + { + // Both positive or both negative + return ab; + } + + // One is positive, one is negative, both are non zero + // If ab is positive, we have to judge if a*b > -c*d (then 1 because sum is positive) + // If ab is negative, we have to judge if c*d > -a*b (idem) + return ab == 1 + ? sign_of_compare(a, b, -c, d) + : sign_of_compare(c, d, -a, b); + } + + 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> + template <typename T, typename Segment, typename Point> static inline T side_value(Segment const& a, Segment const& b, - Segment const& c) + Segment const& c, Point const& fallback_point) { // The first point of the three segments is reused several times T const ax = get<0, 0>(a); @@ -67,9 +252,15 @@ public : 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; + // Assuming they intersect (this method should be called for + // segments known to intersect), they are collinear and overlap. + // They have one or two intersection points - we don't know and + // have to rely on the fallback intersection point + + Point c1, c2; + geometry::detail::assign_point_from_index<0>(c, c1); + geometry::detail::assign_point_from_index<1>(c, c2); + return side_by_triangle<>::apply(c1, c2, fallback_point); } // Cramer's rule: da (see cart_intersect.hpp) @@ -82,7 +273,9 @@ public : // 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> + +#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG + T const result1 = geometry::detail::determinant<T> ( dx_c * d, dy_c * d, d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da @@ -93,15 +286,52 @@ public : // Therefore, the sign is always the same as that result, and the // resulting side (left,right,collinear) is the same + // The first row we divide again by d because of determinant multiply rule + T const result2 = d * geometry::detail::determinant<T> + ( + dx_c, dy_c, + d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da + ); + // Write out: + T const result3 = d * (dx_c * (d * (ay - cy) + dy_a * da) + - dy_c * (d * (ax - cx) + dx_a * da)); + // Write out in braces: + T const result4 = d * (dx_c * d * (ay - cy) + dx_c * dy_a * da + - dy_c * d * (ax - cx) - dy_c * dx_a * da); + // Write in terms of d * XX + da * YY + T const result5 = d * (d * (dx_c * (ay - cy) - dy_c * (ax - cx)) + + da * (dx_c * dy_a - dy_c * dx_a)); + + boost::ignore_unused(result1, result2, result3, result4, result5); + //return result; +#endif + + // We consider the results separately + // (in the end we only have to return the side-value 1,0 or -1) + + // To avoid multiplications we judge the product (easy, avoids *d) + // and the sign of p*q+r*s (more elaborate) + T const result = sign_of_product + ( + d, + sign_of_addition_of_two_products + ( + d, dx_c * (ay - cy) - dy_c * (ax - cx), + da, dx_c * dy_a - dy_c * dx_a + ) + ); return result; + } - template <typename Segment> - static inline int apply(Segment const& a, Segment const& b, Segment const& c) + template <typename Segment, typename Point> + static inline int apply(Segment const& a, Segment const& b, + Segment const& c, + Point const& fallback_point) { typedef typename geometry::coordinate_type<Segment>::type coordinate_type; - coordinate_type const s = side_value<coordinate_type>(a, b, c); + coordinate_type const s = side_value<coordinate_type>(a, b, c, fallback_point); coordinate_type const zero = coordinate_type(); return math::equals(s, zero) ? 0 : s > zero ? 1 diff --git a/boost/geometry/strategies/geographic/distance_thomas.hpp b/boost/geometry/strategies/geographic/distance_thomas.hpp new file mode 100644 index 0000000000..7252b723dd --- /dev/null +++ b/boost/geometry/strategies/geographic/distance_thomas.hpp @@ -0,0 +1,156 @@ +// Boost.Geometry + +// Copyright (c) 2007-2012 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 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_THOMAS_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_THOMAS_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/thomas_inverse.hpp> + +namespace boost { namespace geometry +{ + +namespace strategy { namespace distance +{ + +/*! +\brief The solution of the inverse problem of geodesics on latlong coordinates, + Forsyth-Andoyer-Lambert type approximation with second order terms. +\ingroup distance +\tparam Spheroid The reference spheroid model +\tparam CalculationType \tparam_calculation +\author See + - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965 + http://www.dtic.mil/docs/citations/AD0627893 + - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970 + http://www.dtic.mil/docs/citations/AD703541 +*/ +template +< + typename Spheroid, + typename CalculationType = void +> +class thomas +{ +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 thomas() + : m_spheroid() + {} + + explicit inline thomas(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::thomas_inverse + < + typename calculation_type<Point1, Point2>::type, + true, false + >::apply(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<thomas<Spheroid, CalculationType> > +{ + typedef strategy_tag_distance_point_point type; +}; + + +template <typename Spheroid, typename CalculationType, typename P1, typename P2> +struct return_type<thomas<Spheroid, CalculationType>, P1, P2> + : thomas<Spheroid, CalculationType>::template calculation_type<P1, P2> +{}; + + +template <typename Spheroid, typename CalculationType> +struct comparable_type<thomas<Spheroid, CalculationType> > +{ + typedef thomas<Spheroid, CalculationType> type; +}; + + +template <typename Spheroid, typename CalculationType> +struct get_comparable<thomas<Spheroid, CalculationType> > +{ + static inline thomas<Spheroid, CalculationType> apply(thomas<Spheroid, CalculationType> const& input) + { + return input; + } +}; + +template <typename Spheroid, typename CalculationType, typename P1, typename P2> +struct result_from_distance<thomas<Spheroid, CalculationType>, P1, P2 > +{ + template <typename T> + static inline typename return_type<thomas<Spheroid, CalculationType>, P1, P2>::type + apply(thomas<Spheroid, CalculationType> const& , T const& value) + { + return value; + } +}; + + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::distance + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_THOMAS_HPP diff --git a/boost/geometry/strategies/geographic/distance_vincenty.hpp b/boost/geometry/strategies/geographic/distance_vincenty.hpp index 56dd14bbdc..65bfa19939 100644 --- a/boost/geometry/strategies/geographic/distance_vincenty.hpp +++ b/boost/geometry/strategies/geographic/distance_vincenty.hpp @@ -82,12 +82,13 @@ public : { 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(); + typename calculation_type<Point1, Point2>::type, + true, false + >::apply(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 diff --git a/boost/geometry/strategies/geographic/side_andoyer.hpp b/boost/geometry/strategies/geographic/side_andoyer.hpp new file mode 100644 index 0000000000..e0f0c04067 --- /dev/null +++ b/boost/geometry/strategies/geographic/side_andoyer.hpp @@ -0,0 +1,55 @@ +// Boost.Geometry + +// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2014, 2015. +// Modifications copyright (c) 2014-2015 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_SIDE_ANDOYER_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_ANDOYER_HPP + + +#include <boost/geometry/algorithms/detail/andoyer_inverse.hpp> + +#include <boost/geometry/strategies/geographic/side_detail.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 andoyer + : public detail::by_azimuth<geometry::detail::andoyer_inverse, Model, CalculationType> +{ + typedef detail::by_azimuth<geometry::detail::andoyer_inverse, Model, CalculationType> base_t; + +public: + andoyer(Model const& model = Model()) + : base_t(model) + {} +}; + +}} // namespace strategy::side + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_ANDOYER_HPP diff --git a/boost/geometry/strategies/geographic/side_detail.hpp b/boost/geometry/strategies/geographic/side_detail.hpp new file mode 100644 index 0000000000..c00213d072 --- /dev/null +++ b/boost/geometry/strategies/geographic/side_detail.hpp @@ -0,0 +1,139 @@ +// Boost.Geometry + +// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2014, 2015. +// Modifications copyright (c) 2014-2015 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_SIDE_DETAIL_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_DETAIL_HPP + +#include <boost/geometry/core/cs.hpp> +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/radian_access.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/concepts/side_concept.hpp> + + +namespace boost { namespace geometry +{ + + +namespace strategy { namespace side +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +/*! +\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 InverseFormula Geodesic inverse solution formula. +\tparam Model Reference model of coordinate system. +\tparam CalculationType \tparam_calculation + */ +template <template<typename, bool, bool> class InverseFormula, + typename Model, + typename CalculationType = void> +class by_azimuth +{ +public: + 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; + + typedef InverseFormula<calc_t, false, true> inverse_formula; + + calc_t a1p = azimuth<calc_t, inverse_formula>(p1, p, m_model); + calc_t a12 = azimuth<calc_t, inverse_formula>(p1, p2, m_model); + + calc_t const pi = math::pi<calc_t>(); + + // instead of the formula from XTD + //calc_t a_diff = asin(sin(a1p - a12)); + + calc_t a_diff = a1p - a12; + // normalize, angle in [-pi, pi] + while ( a_diff > pi ) + a_diff -= calc_t(2) * pi; + while ( a_diff < -pi ) + a_diff += calc_t(2) * pi; + + // NOTE: in general it shouldn't be required to support the pi/-pi case + // because in non-cartesian systems it makes sense to check the side + // only "between" the endpoints. + // However currently the winding strategy calls the side strategy + // for vertical segments to check if the point is "between the endpoints. + // This could be avoided since the side strategy is not required for that + // because meridian is the shortest path. So a difference of + // longitudes would be sufficient (of course normalized to [-pi, pi]). + + // NOTE: with the above said, the pi/-pi check is temporary + // however in case if this was required + // the geodesics on ellipsoid aren't "symmetrical" + // therefore instead of comparing a_diff to pi and -pi + // one should probably use inverse azimuths and compare + // the difference to 0 as well + + // positive azimuth is on the right side + return math::equals(a_diff, 0) + || math::equals(a_diff, pi) + || math::equals(a_diff, -pi) ? 0 + : a_diff > 0 ? -1 // right + : 1; // left + } + +private: + template <typename ResultType, + typename InverseFormulaType, + typename Point1, + typename Point2, + typename ModelT> + static inline ResultType azimuth(Point1 const& point1, Point2 const& point2, ModelT const& model) + { + return InverseFormulaType::apply(get_as_radian<0>(point1), + get_as_radian<1>(point1), + get_as_radian<0>(point2), + get_as_radian<1>(point2), + model).azimuth; + } + + Model m_model; +}; + +} // detail +#endif // DOXYGEN_NO_DETAIL + +}} // namespace strategy::side + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_DETAIL_HPP diff --git a/boost/geometry/strategies/geographic/side_thomas.hpp b/boost/geometry/strategies/geographic/side_thomas.hpp new file mode 100644 index 0000000000..a96cabdaab --- /dev/null +++ b/boost/geometry/strategies/geographic/side_thomas.hpp @@ -0,0 +1,55 @@ +// Boost.Geometry + +// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2014, 2015. +// Modifications copyright (c) 2014-2015 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_SIDE_THOMAS_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_THOMAS_HPP + + +#include <boost/geometry/algorithms/detail/thomas_inverse.hpp> + +#include <boost/geometry/strategies/geographic/side_detail.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 thomas + : public detail::by_azimuth<geometry::detail::thomas_inverse, Model, CalculationType> +{ + typedef detail::by_azimuth<geometry::detail::thomas_inverse, Model, CalculationType> base_t; + +public: + thomas(Model const& model = Model()) + : base_t(model) + {} +}; + +}} // namespace strategy::side + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_THOMAS_HPP diff --git a/boost/geometry/strategies/geographic/side_vincenty.hpp b/boost/geometry/strategies/geographic/side_vincenty.hpp new file mode 100644 index 0000000000..a8ef9550d6 --- /dev/null +++ b/boost/geometry/strategies/geographic/side_vincenty.hpp @@ -0,0 +1,55 @@ +// Boost.Geometry + +// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2014, 2015. +// Modifications copyright (c) 2014-2015 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_SIDE_VINCENTY_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_VINCENTY_HPP + + +#include <boost/geometry/algorithms/detail/vincenty_inverse.hpp> + +#include <boost/geometry/strategies/geographic/side_detail.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 vincenty + : public detail::by_azimuth<geometry::detail::vincenty_inverse, Model, CalculationType> +{ + typedef detail::by_azimuth<geometry::detail::vincenty_inverse, Model, CalculationType> base_t; + +public: + vincenty(Model const& model = Model()) + : base_t(model) + {} +}; + +}} // namespace strategy::side + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_VINCENTY_HPP diff --git a/boost/geometry/strategies/intersection_result.hpp b/boost/geometry/strategies/intersection_result.hpp index b467b92a7f..9c6d17ec5e 100644 --- a/boost/geometry/strategies/intersection_result.hpp +++ b/boost/geometry/strategies/intersection_result.hpp @@ -1,6 +1,11 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2007-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 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 @@ -20,140 +25,6 @@ namespace boost { namespace geometry { -/*! - \brief Dimensionally Extended 9 Intersection Matrix - \details - \ingroup overlay - \see http://gis.hsr.ch/wiki/images/3/3d/9dem_springer.pdf -*/ -struct de9im -{ - int ii, ib, ie, - bi, bb, be, - ei, eb, ee; - - inline de9im() - : ii(-1), ib(-1), ie(-1) - , bi(-1), bb(-1), be(-1) - , ei(-1), eb(-1), ee(-1) - { - } - - inline de9im(int ii0, int ib0, int ie0, - int bi0, int bb0, int be0, - int ei0, int eb0, int ee0) - : ii(ii0), ib(ib0), ie(ie0) - , bi(bi0), bb(bb0), be(be0) - , ei(ei0), eb(eb0), ee(ee0) - {} - - inline bool equals() const - { - return ii >= 0 && ie < 0 && be < 0 && ei < 0 && eb < 0; - } - - inline bool disjoint() const - { - return ii < 0 && ib < 0 && bi < 0 && bb < 0; - } - - inline bool intersects() const - { - return ii >= 0 || bb >= 0 || bi >= 0 || ib >= 0; - } - - inline bool touches() const - { - return ii < 0 && (bb >= 0 || bi >= 0 || ib >= 0); - } - - inline bool crosses() const - { - return (ii >= 0 && ie >= 0) || (ii == 0); - } - - inline bool overlaps() const - { - return ii >= 0 && ie >= 0 && ei >= 0; - } - - inline bool within() const - { - return ii >= 0 && ie < 0 && be < 0; - } - - inline bool contains() const - { - return ii >= 0 && ei < 0 && eb < 0; - } - - - static inline char as_char(int v) - { - return v >= 0 && v < 10 ? static_cast<char>('0' + v) : '-'; - } - -#if defined(HAVE_MATRIX_AS_STRING) - inline std::string matrix_as_string(std::string const& tab, std::string const& nl) const - { - std::string ret; - ret.reserve(9 + tab.length() * 3 + nl.length() * 3); - ret += tab; ret += as_char(ii); ret += as_char(ib); ret += as_char(ie); ret += nl; - ret += tab; ret += as_char(bi); ret += as_char(bb); ret += as_char(be); ret += nl; - ret += tab; ret += as_char(ei); ret += as_char(eb); ret += as_char(ee); - return ret; - } - - inline std::string matrix_as_string() const - { - return matrix_as_string("", ""); - } -#endif - -}; - -struct de9im_segment : public de9im -{ - bool collinear; // true if segments are aligned (for equal,overlap,touch) - bool opposite; // true if direction is reversed (for equal,overlap,touch) - bool parallel; // true if disjoint but parallel - bool degenerate; // true for segment(s) of zero length - - double ra, rb; // temp - - inline de9im_segment() - : de9im() - , collinear(false) - , opposite(false) - , parallel(false) - , degenerate(false) - {} - - inline de9im_segment(double a, double b, - int ii0, int ib0, int ie0, - int bi0, int bb0, int be0, - int ei0, int eb0, int ee0, - bool c = false, bool o = false, bool p = false, bool d = false) - : de9im(ii0, ib0, ie0, bi0, bb0, be0, ei0, eb0, ee0) - , collinear(c) - , opposite(o) - , parallel(p) - , degenerate(d) - , ra(a), rb(b) - {} - - -#if defined(HAVE_MATRIX_AS_STRING) - inline std::string as_string() const - { - std::string ret = matrix_as_string(); - ret += collinear ? "c" : "-"; - ret += opposite ? "o" : "-"; - return ret; - } -#endif -}; - template <typename SegmentRatio> struct fraction_type { @@ -209,7 +80,7 @@ struct segment_intersection_info typedef PromotedType promoted_type; CoordinateType dx_a, dy_a; - CoordinateType dx_b, dy_b; // TODO b can be removed + CoordinateType dx_b, dy_b; SegmentRatio robust_ra; SegmentRatio robust_rb; }; diff --git a/boost/geometry/strategies/spherical/area_huiller.hpp b/boost/geometry/strategies/spherical/area_huiller.hpp index e55fdb2b0e..37d8d20124 100644 --- a/boost/geometry/strategies/spherical/area_huiller.hpp +++ b/boost/geometry/strategies/spherical/area_huiller.hpp @@ -1,6 +1,12 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2007-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 +// 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 @@ -31,18 +37,22 @@ namespace strategy { namespace area \tparam PointOfSegment point type of segments of rings/polygons \tparam CalculationType \tparam_calculation \author Barend Gehrels. Adapted from: -- http://www.soe.ucsc.edu/~pang/160/f98/Gems/GemsIV/sph_poly.c +- http://webdocs.cs.ualberta.ca/~graphics/books/GraphicsGems/gemsiv/sph_poly.c +- http://tog.acm.org/resources/GraphicsGems/gemsiv/sph_poly.c - http://williams.best.vwh.net/avform.htm -\note The version in Gems didn't account for polygons crossing the 180 meridian. +\note The version in Graphics Gems IV (page 132-137) didn't account for +polygons crossing the 0 and 180 meridians. The fix for this algorithm +can be found in Graphics Gems V (pages 45-46). See: +- http://kysmykseka.net/koti/wizardry/Game%20Development/Programming/Graphics%20Gems%204.pdf +- http://kysmykseka.net/koti/wizardry/Game%20Development/Programming/Graphics%20Gems%205.pdf \note This version works for convex and non-convex polygons, for 180 meridian crossing polygons and for polygons with holes. However, some cases (especially 180 meridian cases) must still be checked. \note The version which sums angles, which is often seen, doesn't handle non-convex polygons correctly. -\note The version which sums longitudes, see -http://trs-new.jpl.nasa.gov/dspace/bitstream/2014/40409/1/07-03.pdf, is simple -and works well in most cases but not in 180 meridian crossing cases. This probably -could be solved. +\note The version which sums longitudes, see http://hdl.handle.net/2014/40409, +is simple and works well in most cases but not in 180 meridian crossing cases. +This probably could be solved. \note This version is made for spherical equatorial coordinate systems @@ -113,8 +123,12 @@ public : calculation_type const half = 0.5; calculation_type const two = 2.0; calculation_type const four = 4.0; - calculation_type const two_pi = two * geometry::math::pi<calculation_type>(); - calculation_type const half_pi = half * geometry::math::pi<calculation_type>(); + calculation_type const pi + = geometry::math::pi<calculation_type>(); + calculation_type const two_pi + = geometry::math::two_pi<calculation_type>(); + calculation_type const half_pi + = geometry::math::half_pi<calculation_type>(); // Distance p1 p2 calculation_type a = state.distance_over_unit_sphere.apply(p1, p2); @@ -128,33 +142,31 @@ public : // E: spherical excess, using l'Huiller's formula // [tg(e / 4)]2 = tg[s / 2] tg[(s-a) / 2] tg[(s-b) / 2] tg[(s-c) / 2] - calculation_type E = four + calculation_type excess = four * atan(geometry::math::sqrt(geometry::math::abs(tan(s / two) * tan((s - a) / two) * tan((s - b) / two) * tan((s - c) / two)))); - E = geometry::math::abs(E); + excess = geometry::math::abs(excess); // In right direction: positive, add area. In left direction: negative, subtract area. - // Longitude comparisons are not so obvious. If one is negative, other is positive, + // Longitude comparisons are not so obvious. If one is negative and other is positive, // we have to take the dateline into account. - // TODO: check this / enhance this, should be more robust. See also the "grow" for ll - // TODO: use minmax or "smaller"/"compare" strategy for this - calculation_type lon1 = geometry::get_as_radian<0>(p1) < 0 - ? geometry::get_as_radian<0>(p1) + two_pi - : geometry::get_as_radian<0>(p1); - calculation_type lon2 = geometry::get_as_radian<0>(p2) < 0 - ? geometry::get_as_radian<0>(p2) + two_pi - : geometry::get_as_radian<0>(p2); + calculation_type lon_diff = geometry::get_as_radian<0>(p2) + - geometry::get_as_radian<0>(p1); + if (lon_diff <= 0) + { + lon_diff += two_pi; + } - if (lon2 < lon1) + if (lon_diff > pi) { - E = -E; + excess = -excess; } - state.sum += E; + state.sum += excess; } } diff --git a/boost/geometry/strategies/spherical/distance_cross_track.hpp b/boost/geometry/strategies/spherical/distance_cross_track.hpp index 486c7effbd..31b59e77ff 100644 --- a/boost/geometry/strategies/spherical/distance_cross_track.hpp +++ b/boost/geometry/strategies/spherical/distance_cross_track.hpp @@ -340,6 +340,8 @@ public : > {}; + typedef typename Strategy::radius_type radius_type; + inline cross_track() {} @@ -372,11 +374,16 @@ 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; + std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " " + << crs_AD * geometry::math::r2d<return_type>() << std::endl; + std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " " + << crs_AB * geometry::math::r2d<return_type>() << 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<return_type>() << std::endl; + std::cout << "Projection BD-BA " << projection2 << " : " + << d_crs2 * geometry::math::r2d<return_type>() << std::endl; #endif // http://williams.best.vwh.net/avform.htm#XTE @@ -487,6 +494,8 @@ public : > {}; + typedef typename Strategy::radius_type radius_type; + inline cross_track() {} diff --git a/boost/geometry/strategies/spherical/distance_cross_track_point_box.hpp b/boost/geometry/strategies/spherical/distance_cross_track_point_box.hpp index fe32f77236..14c5bd4eda 100644 --- a/boost/geometry/strategies/spherical/distance_cross_track_point_box.hpp +++ b/boost/geometry/strategies/spherical/distance_cross_track_point_box.hpp @@ -1,13 +1,14 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2008-2014 Bruno Lalande, Paris, France. -// Copyright (c) 2008-2014 Barend Gehrels, Amsterdam, the Netherlands. -// Copyright (c) 2009-2014 Mateusz Loskot, London, UK. +// Copyright (c) 2008-2015 Bruno Lalande, Paris, France. +// Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2009-2015 Mateusz Loskot, London, UK. -// 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 Adam Wulkiewicz, on behalf of Oracle +// 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 @@ -16,14 +17,23 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_POINT_BOX_HPP #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_POINT_BOX_HPP +#include <boost/config.hpp> +#include <boost/concept_check.hpp> +#include <boost/mpl/if.hpp> +#include <boost/type_traits/is_void.hpp> #include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/assert.hpp> +#include <boost/geometry/core/point_type.hpp> +#include <boost/geometry/core/radian_access.hpp> +#include <boost/geometry/core/tags.hpp> #include <boost/geometry/strategies/distance.hpp> +#include <boost/geometry/strategies/concepts/distance_concept.hpp> +#include <boost/geometry/strategies/spherical/distance_cross_track.hpp> #include <boost/geometry/util/math.hpp> -#include <boost/geometry/util/calculation_type.hpp> - +#include <boost/geometry/algorithms/detail/assign_box_corners.hpp> namespace boost { namespace geometry @@ -32,131 +42,197 @@ namespace boost { namespace geometry namespace strategy { namespace distance { + +/*! +\brief Strategy functor for distance point to box calculation +\ingroup strategies +\details Class which calculates the distance of a point to a box, for +points and boxes on a sphere or globe +\tparam CalculationType \tparam_calculation +\tparam Strategy underlying point-segment distance strategy, defaults +to cross track + +\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> + typename Strategy = cross_track<CalculationType> > class cross_track_point_box { public: template <typename Point, typename Box> struct return_type - : promote_floating_point - < - typename select_calculation_type - < - Point, - typename point_type<Box>::type, - CalculationType - >::type - > + : services::return_type<Strategy, Point, typename point_type<Box>::type> {}; + typedef typename Strategy::radius_type radius_type; + inline cross_track_point_box() {} explicit inline cross_track_point_box(typename Strategy::radius_type const& r) - : m_pp_strategy(r) + : m_ps_strategy(r) {} inline cross_track_point_box(Strategy const& s) - : m_pp_strategy(s) + : m_ps_strategy(s) {} - + + + // It might be useful in the future + // to overload constructor with strategy info. + // crosstrack(...) {} + template <typename Point, typename Box> inline typename return_type<Point, Box>::type apply(Point const& point, Box const& box) const { - #if !defined(BOOST_MSVC) BOOST_CONCEPT_ASSERT ( - (concept::PointDistanceStrategy + (concept::PointSegmentDistanceStrategy < - Strategy, Point, - typename point_type<Box>::type + Strategy, Point, typename point_type<Box>::type >) ); #endif + // this method assumes that the coordinates of the point and + // the box are normalized + typedef typename return_type<Point, Box>::type return_type; - typedef typename point_type<Box>::type box_point_t; - - // Create (counterclockwise) array of points, the fifth one closes it - // If every point is on the LEFT side (=1) or ON the border (=0) - // the distance should be equal to 0. + typedef typename point_type<Box>::type box_point_type; // TODO: This strategy as well as other cross-track strategies // and therefore e.g. spherical within(Point, Box) may not work // properly for a Box degenerated to a Segment or Point - boost::array<box_point_t, 5> bp; - geometry::detail::assign_box_corners_oriented<true>(box, bp); - bp[4] = bp[0]; + box_point_type bottom_left, bottom_right, top_left, top_right; + geometry::detail::assign_box_corners(box, + bottom_left, bottom_right, + top_left, top_right); + + return_type const plon = geometry::get_as_radian<0>(point); + return_type const plat = geometry::get_as_radian<1>(point); + + return_type const lon_min = geometry::get_as_radian<0>(bottom_left); + return_type const lat_min = geometry::get_as_radian<1>(bottom_left); + return_type const lon_max = geometry::get_as_radian<0>(top_right); + return_type const lat_max = geometry::get_as_radian<1>(top_right); + + return_type const pi = math::pi<return_type>(); + return_type const two_pi = math::two_pi<return_type>(); + + // First check if the point is within the band defined by the + // minimum and maximum longitude of the box; if yes, determine + // if the point is above, below or inside the box and compute + // the distance (easy in this case) + // + // Notice that the point may not be inside the longitude range + // of the box, but the shifted point may be inside the + // longitude range of the box; in this case the point is still + // considered as inside the longitude range band of the box + if ((plon >= lon_min && plon <= lon_max) || plon + two_pi <= lon_max) + { + if (plat > lat_max) + { + return services::result_from_distance + < + Strategy, Point, box_point_type + >::apply(m_ps_strategy, radius() * (plat - lat_max)); + } + else if (plat < lat_min) + { + return services::result_from_distance + < + Strategy, Point, box_point_type + >::apply(m_ps_strategy, radius() * (lat_min - plat)); + } + else + { + BOOST_GEOMETRY_ASSERT(plat >= lat_min && plat <= lat_max); + return return_type(0); + } + } - for (int i = 1; i < 5; i++) + // Otherwise determine which αμονγ the two medirian segments of the + // box the point is closest to, and compute the distance of + // the point to this closest segment + + // Below lon_midway is the longitude of the meridian that: + // (1) is midway between the meridians of the left and right + // meridians of the box, and + // (2) does not intersect the box + return_type const two = 2.0; + bool use_left_segment; + if (lon_max > pi) { - box_point_t const& p1 = bp[i - 1]; - box_point_t const& p2 = bp[i]; + // the box crosses the antimeridian - return_type const crs_AD = geometry::detail::course<return_type>(p1, point); - return_type const crs_AB = geometry::detail::course<return_type>(p1, p2); - return_type const d_crs1 = crs_AD - crs_AB; - return_type const sin_d_crs1 = sin(d_crs1); + // midway longitude = lon_min - (lon_min + (lon_max - 2 * pi)) / 2; + return_type const lon_midway = (lon_min - lon_max) / two + pi; + BOOST_GEOMETRY_ASSERT(lon_midway >= -pi && lon_midway <= pi); - // this constant sin() is here to be consistent with the side strategy - return_type const sigXTD = asin(sin(0.001) * sin_d_crs1); + use_left_segment = plon > lon_midway; + } + else + { + // the box does not cross the antimeridian - // If the point is on the right side of the edge - if ( sigXTD > 0 ) + return_type const lon_sum = lon_min + lon_max; + if (math::equals(lon_sum, return_type(0))) { - return_type const crs_BA = crs_AB - geometry::math::pi<return_type>(); - return_type const crs_BD = geometry::detail::course<return_type>(p2, point); - return_type const d_crs2 = crs_BD - crs_BA; + // special case: the box is symmetric with respect to + // the prime meridian; the midway meridian is the antimeridian - return_type const projection1 = cos( d_crs1 ); - return_type const projection2 = cos( d_crs2 ); + use_left_segment = plon < lon_min; + } + else + { + // midway long. = lon_min - (2 * pi - (lon_max - lon_min)) / 2; + return_type lon_midway = (lon_min + lon_max) / two - pi; - if(projection1 > 0.0 && projection2 > 0.0) + // normalize the midway longitude + if (lon_midway > pi) { - return_type const d1 = m_pp_strategy.apply(p1, point); - return_type const - XTD = radius() - * geometry::math::abs( - asin( sin( d1 / radius() ) * sin_d_crs1 ) - ); - - return return_type(XTD); + lon_midway -= two_pi; } - else + else if (lon_midway < -pi) { - // OPTIMIZATION - // Return d1 if projection1 <= 0 and d2 if projection2 <= 0 - // if both == 0 then return d1 or d2 - // both shouldn't be < 0 - - return_type const d1 = m_pp_strategy.apply(p1, point); - return_type const d2 = m_pp_strategy.apply(p2, point); - - return return_type((std::min)( d1 , d2 )); + lon_midway += two_pi; } + BOOST_GEOMETRY_ASSERT(lon_midway >= -pi && lon_midway <= pi); + + // if lon_sum is positive the midway meridian is left + // of the box, or right of the box otherwise + use_left_segment = lon_sum > 0 + ? (plon < lon_min && plon >= lon_midway) + : (plon <= lon_max || plon > lon_midway); } } - // Return 0 if the point isn't on the right side of any segment - return return_type(0); + return use_left_segment + ? m_ps_strategy.apply(point, bottom_left, top_left) + : m_ps_strategy.apply(point, bottom_right, top_right); } inline typename Strategy::radius_type radius() const - { return m_pp_strategy.radius(); } - -private : + { + return m_ps_strategy.radius(); + } - Strategy m_pp_strategy; +private: + Strategy m_ps_strategy; }; + #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS namespace services { @@ -170,67 +246,68 @@ struct tag<cross_track_point_box<CalculationType, Strategy> > template <typename CalculationType, typename Strategy, typename P, typename Box> struct return_type<cross_track_point_box<CalculationType, Strategy>, P, Box> - : cross_track_point_box<CalculationType, Strategy>::template return_type<P, Box> + : cross_track_point_box + < + CalculationType, Strategy + >::template return_type<P, Box> {}; template <typename CalculationType, typename Strategy> struct comparable_type<cross_track_point_box<CalculationType, Strategy> > { - // There is no shortcut, so the strategy itself is its comparable type - typedef cross_track_point_box<CalculationType, Strategy> type; + typedef cross_track_point_box + < + CalculationType, typename comparable_type<Strategy>::type + > type; }; -template -< - typename CalculationType, - typename Strategy -> +template <typename CalculationType, typename Strategy> struct get_comparable<cross_track_point_box<CalculationType, Strategy> > { - typedef typename comparable_type - < - cross_track_point_box<CalculationType, Strategy> - >::type comparable_type; -public : - static inline comparable_type apply( - cross_track_point_box<CalculationType, Strategy> const& strategy) + typedef cross_track_point_box<CalculationType, Strategy> this_strategy; + typedef typename comparable_type<this_strategy>::type comparable_type; + +public: + static inline comparable_type apply(this_strategy const& strategy) { - return cross_track_point_box<CalculationType, Strategy>(strategy.radius()); + return comparable_type(strategy.radius()); } }; -template -< - typename CalculationType, - typename Strategy, - typename P, typename Box -> +template <typename CalculationType, typename Strategy, typename P, typename Box> struct result_from_distance < - cross_track_point_box<CalculationType, Strategy>, - P, - Box + cross_track_point_box<CalculationType, Strategy>, P, Box > { -private : - typedef typename cross_track_point_box +private: + typedef cross_track_point_box<CalculationType, Strategy> this_strategy; + + typedef typename this_strategy::template return_type < - CalculationType, Strategy - >::template return_type<P, Box> return_type; -public : + P, Box + >::type return_type; + +public: template <typename T> - static inline return_type apply( - cross_track_point_box<CalculationType, Strategy> const& , - T const& distance) + static inline return_type apply(this_strategy const& strategy, + T const& distance) { - return distance; + Strategy s(strategy.radius()); + + return result_from_distance + < + Strategy, P, typename point_type<Box>::type + >::apply(s, distance); } }; +// define cross_track_point_box<default_point_segment_strategy> as +// default point-box strategy for the spherical equatorial coordinate system template <typename Point, typename Box, typename Strategy> struct default_strategy < @@ -247,7 +324,7 @@ struct default_strategy boost::is_void<Strategy>, typename default_strategy < - point_tag, point_tag, + point_tag, segment_tag, Point, typename point_type<Box>::type, spherical_equatorial_tag, spherical_equatorial_tag >::type, diff --git a/boost/geometry/strategies/strategies.hpp b/boost/geometry/strategies/strategies.hpp index bf436dba2b..28850020af 100644 --- a/boost/geometry/strategies/strategies.hpp +++ b/boost/geometry/strategies/strategies.hpp @@ -4,8 +4,10 @@ // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. // Copyright (c) 2009-2012 Mateusz Loskot, London, UK. -// 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 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. @@ -14,8 +16,6 @@ // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) -// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle - #ifndef BOOST_GEOMETRY_STRATEGIES_STRATEGIES_HPP #define BOOST_GEOMETRY_STRATEGIES_STRATEGIES_HPP @@ -64,7 +64,11 @@ #include <boost/geometry/strategies/spherical/ssf.hpp> #include <boost/geometry/strategies/geographic/distance_andoyer.hpp> +#include <boost/geometry/strategies/geographic/distance_thomas.hpp> #include <boost/geometry/strategies/geographic/distance_vincenty.hpp> +//#include <boost/geometry/strategies/geographic/side_andoyer.hpp> +//#include <boost/geometry/strategies/geographic/side_thomas.hpp> +//#include <boost/geometry/strategies/geographic/side_vincenty.hpp> #include <boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp> #include <boost/geometry/strategies/agnostic/buffer_distance_asymmetric.hpp> diff --git a/boost/geometry/strategies/strategy_transform.hpp b/boost/geometry/strategies/strategy_transform.hpp index 9d3b1d910c..99fcb720c6 100644 --- a/boost/geometry/strategies/strategy_transform.hpp +++ b/boost/geometry/strategies/strategy_transform.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. @@ -28,6 +33,7 @@ #include <boost/geometry/strategies/transform.hpp> #include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/promote_floating_point.hpp> #include <boost/geometry/util/select_coordinate_type.hpp> namespace boost { namespace geometry @@ -130,7 +136,15 @@ struct degree_radian_vv assert_dimension<P1, 2>(); assert_dimension<P2, 2>(); - detail::transform_coordinates<P1, P2, 0, 2, F>::transform(p1, p2, math::d2r); + typedef typename promote_floating_point + < + typename select_coordinate_type<P1, P2>::type + >::type calculation_type; + + detail::transform_coordinates + < + P1, P2, 0, 2, F + >::transform(p1, p2, math::d2r<calculation_type>()); return true; } }; @@ -143,7 +157,16 @@ struct degree_radian_vv_3 assert_dimension<P1, 3>(); assert_dimension<P2, 3>(); - detail::transform_coordinates<P1, P2, 0, 2, F>::transform(p1, p2, math::d2r); + typedef typename promote_floating_point + < + typename select_coordinate_type<P1, P2>::type + >::type calculation_type; + + detail::transform_coordinates + < + P1, P2, 0, 2, F + >::transform(p1, p2, math::d2r<calculation_type>()); + // Copy height or other third dimension set<2>(p2, get<2>(p1)); return true; diff --git a/boost/geometry/strategies/transform/matrix_transformers.hpp b/boost/geometry/strategies/transform/matrix_transformers.hpp index 27a3a2ae80..699b91b3aa 100644 --- a/boost/geometry/strategies/transform/matrix_transformers.hpp +++ b/boost/geometry/strategies/transform/matrix_transformers.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. @@ -40,6 +45,7 @@ #include <boost/geometry/core/coordinate_dimension.hpp> #include <boost/geometry/core/cs.hpp> #include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/promote_floating_point.hpp> #include <boost/geometry/util/select_coordinate_type.hpp> #include <boost/geometry/util/select_most_precise.hpp> @@ -340,7 +346,8 @@ struct as_radian<degree> template <typename T> static inline T get(T const& value) { - return value * math::d2r; + typedef typename promote_floating_point<T>::type promoted_type; + return value * math::d2r<promoted_type>(); } }; |