diff options
Diffstat (limited to 'boost/geometry/strategies/spherical')
9 files changed, 975 insertions, 199 deletions
diff --git a/boost/geometry/strategies/spherical/area.hpp b/boost/geometry/strategies/spherical/area.hpp index 1ab61b644f..d12ed9498f 100644 --- a/boost/geometry/strategies/spherical/area.hpp +++ b/boost/geometry/strategies/spherical/area.hpp @@ -1,6 +1,8 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2016-2017 Oracle and/or its affiliates. +// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland. + +// Copyright (c) 2016-2018 Oracle and/or its affiliates. // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -13,9 +15,9 @@ #include <boost/geometry/formulas/area_formulas.hpp> -#include <boost/geometry/core/radius.hpp> -#include <boost/geometry/core/srs.hpp> +#include <boost/geometry/srs/sphere.hpp> #include <boost/geometry/strategies/area.hpp> +#include <boost/geometry/strategies/spherical/get_radius.hpp> namespace boost { namespace geometry @@ -24,11 +26,12 @@ namespace boost { namespace geometry namespace strategy { namespace area { + /*! \brief Spherical area calculation \ingroup strategies \details Calculates area on the surface of a sphere using the trapezoidal rule -\tparam PointOfSegment \tparam_segment_point +\tparam RadiusTypeOrSphere \tparam_radius_or_sphere \tparam CalculationType \tparam_calculation \qbk{ @@ -38,7 +41,7 @@ namespace strategy { namespace area */ template < - typename PointOfSegment, + typename RadiusTypeOrSphere = double, typename CalculationType = void > class spherical @@ -46,34 +49,35 @@ class spherical // Enables special handling of long segments static const bool LongSegment = false; -typedef typename boost::mpl::if_c - < - boost::is_void<CalculationType>::type::value, - typename select_most_precise +public: + template <typename Geometry> + struct result_type + : strategy::area::detail::result_type < - typename coordinate_type<PointOfSegment>::type, - double - >::type, - CalculationType - >::type CT; - -protected : - struct excess_sum + Geometry, + CalculationType + > + {}; + + template <typename Geometry> + class state { - CT m_sum; + friend class spherical; - // Keep track if encircles some pole - size_t m_crosses_prime_meridian; + typedef typename result_type<Geometry>::type return_type; - inline excess_sum() + public: + inline state() : m_sum(0) , m_crosses_prime_meridian(0) {} - template <typename SphereType> - inline CT area(SphereType sphere) const + + private: + template <typename RadiusType> + inline return_type area(RadiusType const& r) const { - CT result; - CT radius = geometry::get_radius<0>(sphere); + return_type result; + return_type radius = r; // Encircles pole if(m_crosses_prime_meridian % 2 == 1) @@ -81,12 +85,12 @@ protected : size_t times_crosses_prime_meridian = 1 + (m_crosses_prime_meridian / 2); - result = CT(2) - * geometry::math::pi<CT>() + result = return_type(2) + * geometry::math::pi<return_type>() * times_crosses_prime_meridian - geometry::math::abs(m_sum); - if(geometry::math::sign<CT>(m_sum) == 1) + if(geometry::math::sign<return_type>(m_sum) == 1) { result = - result; } @@ -99,54 +103,62 @@ protected : return result; } + + return_type m_sum; + + // Keep track if encircles some pole + size_t m_crosses_prime_meridian; }; public : - typedef CT return_type; - typedef PointOfSegment segment_point_type; - typedef excess_sum state_type; - typedef geometry::srs::sphere<CT> sphere_type; // For backward compatibility reasons the radius is set to 1 inline spherical() - : m_sphere(1.0) - {} - - template <typename T> - explicit inline spherical(geometry::srs::sphere<T> const& sphere) - : m_sphere(geometry::get_radius<0>(sphere)) + : m_radius(1.0) {} - explicit inline spherical(CT const& radius) - : m_sphere(radius) + template <typename RadiusOrSphere> + explicit inline spherical(RadiusOrSphere const& radius_or_sphere) + : m_radius(strategy_detail::get_radius + < + RadiusOrSphere + >::apply(radius_or_sphere)) {} + template <typename PointOfSegment, typename Geometry> inline void apply(PointOfSegment const& p1, PointOfSegment const& p2, - excess_sum& state) const + state<Geometry>& st) const { if (! geometry::math::equals(get<0>(p1), get<0>(p2))) { - typedef geometry::formula::area_formulas<CT> area_formulas; + typedef geometry::formula::area_formulas + < + typename result_type<Geometry>::type + > area_formulas; - state.m_sum += area_formulas::template spherical<LongSegment>(p1, p2); + st.m_sum += area_formulas::template spherical<LongSegment>(p1, p2); // Keep track whenever a segment crosses the prime meridian if (area_formulas::crosses_prime_meridian(p1, p2)) { - state.m_crosses_prime_meridian++; + st.m_crosses_prime_meridian++; } } } - inline return_type result(excess_sum const& state) const + template <typename Geometry> + inline typename result_type<Geometry>::type + result(state<Geometry> const& st) const { - return state.area(m_sphere); + return st.area(m_radius); } private : - /// srs Sphere - sphere_type m_sphere; + typename strategy_detail::get_radius + < + RadiusTypeOrSphere + >::type m_radius; }; #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS @@ -155,17 +167,17 @@ namespace services { -template <typename Point> -struct default_strategy<spherical_equatorial_tag, Point> +template <> +struct default_strategy<spherical_equatorial_tag> { - typedef strategy::area::spherical<Point> type; + typedef strategy::area::spherical<> type; }; // Note: spherical polar coordinate system requires "get_as_radian_equatorial" -template <typename Point> -struct default_strategy<spherical_polar_tag, Point> +template <> +struct default_strategy<spherical_polar_tag> { - typedef strategy::area::spherical<Point> type; + typedef strategy::area::spherical<> type; }; } // namespace services diff --git a/boost/geometry/strategies/spherical/densify.hpp b/boost/geometry/strategies/spherical/densify.hpp new file mode 100644 index 0000000000..97f4605a98 --- /dev/null +++ b/boost/geometry/strategies/spherical/densify.hpp @@ -0,0 +1,189 @@ +// Boost.Geometry + +// Copyright (c) 2017-2018, Oracle and/or its affiliates. + +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// Licensed under the Boost Software License version 1.0. +// http://www.boost.org/users/license.html + +#ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DENSIFY_HPP +#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DENSIFY_HPP + + +#include <boost/geometry/algorithms/detail/convert_point_to_point.hpp> +#include <boost/geometry/algorithms/detail/signed_size_type.hpp> +#include <boost/geometry/arithmetic/arithmetic.hpp> +#include <boost/geometry/arithmetic/cross_product.hpp> +#include <boost/geometry/arithmetic/dot_product.hpp> +#include <boost/geometry/arithmetic/normalize.hpp> +#include <boost/geometry/core/assert.hpp> +#include <boost/geometry/core/coordinate_dimension.hpp> +#include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/core/radian_access.hpp> +#include <boost/geometry/formulas/spherical.hpp> +#include <boost/geometry/srs/sphere.hpp> +#include <boost/geometry/strategies/densify.hpp> +#include <boost/geometry/strategies/spherical/get_radius.hpp> +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/select_most_precise.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace densify +{ + + +/*! +\brief Densification of spherical segment. +\ingroup strategies +\tparam RadiusTypeOrSphere \tparam_radius_or_sphere +\tparam CalculationType \tparam_calculation + +\qbk{ +[heading See also] +[link geometry.reference.algorithms.densify.densify_4_with_strategy densify (with strategy)] +} + */ +template +< + typename RadiusTypeOrSphere = double, + typename CalculationType = void +> +class spherical +{ +public: + // For consistency with area strategy the radius is set to 1 + inline spherical() + : m_radius(1.0) + {} + + template <typename RadiusOrSphere> + explicit inline spherical(RadiusOrSphere const& radius_or_sphere) + : m_radius(strategy_detail::get_radius + < + RadiusOrSphere + >::apply(radius_or_sphere)) + {} + + template <typename Point, typename AssignPolicy, typename T> + inline void apply(Point const& p0, Point const& p1, AssignPolicy & policy, T const& length_threshold) const + { + typedef typename AssignPolicy::point_type out_point_t; + typedef typename select_most_precise + < + typename coordinate_type<Point>::type, + typename coordinate_type<out_point_t>::type, + CalculationType + >::type calc_t; + + calc_t const c0 = 0; + calc_t const c1 = 1; + calc_t const pi = math::pi<calc_t>(); + + typedef model::point<calc_t, 3, cs::cartesian> point3d_t; + point3d_t const xyz0 = formula::sph_to_cart3d<point3d_t>(p0); + point3d_t const xyz1 = formula::sph_to_cart3d<point3d_t>(p1); + calc_t const dot01 = geometry::dot_product(xyz0, xyz1); + calc_t const angle01 = acos(dot01); + + BOOST_GEOMETRY_ASSERT(length_threshold > T(0)); + + signed_size_type n = signed_size_type(angle01 * m_radius / length_threshold); + if (n <= 0) + return; + + point3d_t axis; + if (! math::equals(angle01, pi)) + { + axis = geometry::cross_product(xyz0, xyz1); + geometry::detail::vec_normalize(axis); + } + else // antipodal + { + calc_t const half_pi = math::half_pi<calc_t>(); + calc_t const lat = geometry::get_as_radian<1>(p0); + + if (math::equals(lat, half_pi)) + { + // pointing east, segment lies on prime meridian, going south + axis = point3d_t(c0, c1, c0); + } + else if (math::equals(lat, -half_pi)) + { + // pointing west, segment lies on prime meridian, going north + axis = point3d_t(c0, -c1, c0); + } + else + { + // lon rotated west by pi/2 at equator + calc_t const lon = geometry::get_as_radian<0>(p0); + axis = point3d_t(sin(lon), -cos(lon), c0); + } + } + + calc_t step = angle01 / (n + 1); + + calc_t a = step; + for (signed_size_type i = 0 ; i < n ; ++i, a += step) + { + // Axis-Angle rotation + // see: https://en.wikipedia.org/wiki/Axis-angle_representation + calc_t const cos_a = cos(a); + calc_t const sin_a = sin(a); + // cos_a * v + point3d_t s1 = xyz0; + geometry::multiply_value(s1, cos_a); + // sin_a * (n x v) + point3d_t s2 = geometry::cross_product(axis, xyz0); + geometry::multiply_value(s2, sin_a); + // (1 - cos_a)(n.v) * n + point3d_t s3 = axis; + geometry::multiply_value(s3, (c1 - cos_a) * geometry::dot_product(axis, xyz0)); + // v_rot = cos_a * v + sin_a * (n x v) + (1 - cos_a)(n.v) * e + point3d_t v_rot = s1; + geometry::add_point(v_rot, s2); + geometry::add_point(v_rot, s3); + + out_point_t p = formula::cart3d_to_sph<out_point_t>(v_rot); + geometry::detail::conversion::point_to_point + < + Point, out_point_t, + 2, dimension<out_point_t>::value + >::apply(p0, p); + + policy.apply(p); + } + } + +private: + typename strategy_detail::get_radius + < + RadiusTypeOrSphere + >::type m_radius; +}; + + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +template <> +struct default_strategy<spherical_equatorial_tag> +{ + typedef strategy::densify::spherical<> type; +}; + + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::densify + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DENSIFY_HPP diff --git a/boost/geometry/strategies/spherical/distance_cross_track.hpp b/boost/geometry/strategies/spherical/distance_cross_track.hpp index 82a2fb3d6f..db029dff90 100644 --- a/boost/geometry/strategies/spherical/distance_cross_track.hpp +++ b/boost/geometry/strategies/spherical/distance_cross_track.hpp @@ -352,6 +352,11 @@ public : : m_strategy(s) {} + //TODO: apply a more general strategy getter + inline Strategy get_distance_strategy() const + { + return m_strategy; + } // It might be useful in the future // to overload constructor with strategy info. @@ -526,6 +531,11 @@ public : : m_strategy(s) {} + //TODO: apply a more general strategy getter + inline Strategy get_distance_strategy() const + { + return m_strategy; + } // It might be useful in the future // to overload constructor with strategy info. diff --git a/boost/geometry/strategies/spherical/distance_cross_track_box_box.hpp b/boost/geometry/strategies/spherical/distance_cross_track_box_box.hpp new file mode 100644 index 0000000000..efa7a728a6 --- /dev/null +++ b/boost/geometry/strategies/spherical/distance_cross_track_box_box.hpp @@ -0,0 +1,440 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2016-2018 Oracle and/or its affiliates. +// Contributed and/or modified by Vissarion Fisikopoulos, 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 +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_BOX_BOX_HPP +#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_BOX_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/algorithms/detail/assign_box_corners.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace distance +{ + +namespace details +{ + +template <typename ReturnType> +class cross_track_box_box_generic +{ +public : + + template <typename Point, typename Strategy> + ReturnType static inline diagonal_case(Point topA, + Point topB, + Point bottomA, + Point bottomB, + bool north_shortest, + bool non_overlap, + Strategy ps_strategy) + { + if (north_shortest && non_overlap) + { + return ps_strategy.get_distance_strategy().apply(topA, bottomB); + } + if (north_shortest && !non_overlap) + { + return ps_strategy.apply(topA, topB, bottomB); + } + if (!north_shortest && non_overlap) + { + return ps_strategy.get_distance_strategy().apply(bottomA, topB); + } + return ps_strategy.apply(bottomA, topB, bottomB); + } + + + template + < + typename Box1, + typename Box2, + typename Strategy + > + ReturnType static inline apply (Box1 const& box1, + Box2 const& box2, + Strategy ps_strategy) + { + + // this method assumes that the coordinates of the point and + // the box are normalized + + typedef typename point_type<Box1>::type box_point_type1; + typedef typename point_type<Box2>::type box_point_type2; + + box_point_type1 bottom_left1, bottom_right1, top_left1, top_right1; + geometry::detail::assign_box_corners(box1, + bottom_left1, bottom_right1, + top_left1, top_right1); + + box_point_type2 bottom_left2, bottom_right2, top_left2, top_right2; + geometry::detail::assign_box_corners(box2, + bottom_left2, bottom_right2, + top_left2, top_right2); + + ReturnType lon_min1 = geometry::get_as_radian<0>(bottom_left1); + ReturnType const lat_min1 = geometry::get_as_radian<1>(bottom_left1); + ReturnType lon_max1 = geometry::get_as_radian<0>(top_right1); + ReturnType const lat_max1 = geometry::get_as_radian<1>(top_right1); + + ReturnType lon_min2 = geometry::get_as_radian<0>(bottom_left2); + ReturnType const lat_min2 = geometry::get_as_radian<1>(bottom_left2); + ReturnType lon_max2 = geometry::get_as_radian<0>(top_right2); + ReturnType const lat_max2 = geometry::get_as_radian<1>(top_right2); + + ReturnType const two_pi = math::two_pi<ReturnType>(); + + // Test which sides of the boxes are closer and if boxes cross + // antimeridian + bool right_wrap; + + if (lon_min2 > 0 && lon_max2 < 0) // box2 crosses antimeridian + { +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK_BOX_BOX + std::cout << "(box2 crosses antimeridian)"; +#endif + right_wrap = lon_min2 - lon_max1 < lon_min1 - lon_max2; + lon_max2 += two_pi; + if (lon_min1 > 0 && lon_max1 < 0) // both boxes crosses antimeridian + { + lon_max1 += two_pi; + } + } + else if (lon_min1 > 0 && lon_max1 < 0) // only box1 crosses antimeridian + { +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK_BOX_BOX + std::cout << "(box1 crosses antimeridian)"; +#endif + return apply(box2, box1, ps_strategy); + } + else + { + right_wrap = lon_max1 <= lon_min2 + ? lon_min2 - lon_max1 < two_pi - (lon_max2 - lon_min1) + : lon_min1 - lon_max2 > two_pi - (lon_max1 - lon_min2); + + } + + // Check1: if box2 crosses the band defined by the + // minimum and maximum longitude of box1; if yes, determine + // if the box2 is above, below or intersects/is inside box1 and compute + // the distance (easy in this case) + + bool lon_min12 = lon_min1 <= lon_min2; + bool right = lon_max1 <= lon_min2; + bool left = lon_min1 >= lon_max2; + bool lon_max12 = lon_max1 <= lon_max2; + + if ((lon_min12 && !right) + || (!left && !lon_max12) + || (!lon_min12 && lon_max12)) + { +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK_BOX_BOX + std::cout << "(up-down)\n"; +#endif + if (lat_min1 > lat_max2) + { + return geometry::strategy::distance::services::result_from_distance + < + Strategy, box_point_type1, box_point_type2 + >::apply(ps_strategy, ps_strategy.get_distance_strategy() + .meridian(lat_min1, lat_max2)); + } + else if (lat_max1 < lat_min2) + { + return geometry::strategy::distance::services::result_from_distance + < + Strategy, box_point_type1, box_point_type2 + >::apply(ps_strategy, ps_strategy.get_distance_strategy(). + meridian(lat_min2, lat_max1)); + } + else + { + //BOOST_GEOMETRY_ASSERT(plat >= lat_min && plat <= lat_max); + return ReturnType(0); + } + } + + // Check2: if box2 is right/left of box1 + // the max lat of box2 should be less than the max lat of box1 + bool bottom_max; + + ReturnType top_common = (std::min)(lat_max1, lat_max2); + ReturnType bottom_common = (std::max)(lat_min1, lat_min2); + + // true if the closest points are on northern hemisphere + bool north_shortest = top_common + bottom_common > 0; + // true if box bands do not overlap + bool non_overlap = top_common < bottom_common; + + if (north_shortest) + { + bottom_max = lat_max1 >= lat_max2; + } + else + { + bottom_max = lat_min1 <= lat_min2; + } + +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK_BOX_BOX + std::cout << "(diagonal)"; +#endif + if (bottom_max && !right_wrap) + { +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK_BOX_BOX + std::cout << "(bottom left)"; +#endif + return diagonal_case(top_right2, top_left1, + bottom_right2, bottom_left1, + north_shortest, non_overlap, + ps_strategy); + } + if (bottom_max && right_wrap) + { +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK_BOX_BOX + std::cout << "(bottom right)"; +#endif + return diagonal_case(top_left2, top_right1, + bottom_left2, bottom_right1, + north_shortest, non_overlap, + ps_strategy); + } + if (!bottom_max && !right_wrap) + { +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK_BOX_BOX + std::cout << "(top left)"; +#endif + return diagonal_case(top_left1, top_right2, + bottom_left1, bottom_right2, + north_shortest, non_overlap, + ps_strategy); + } + if (!bottom_max && right_wrap) + { +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK_BOX_BOX + std::cout << "(top right)"; +#endif + return diagonal_case(top_right1, top_left2, + bottom_right1, bottom_left2, + north_shortest, non_overlap, + ps_strategy); + } + return ReturnType(0); + } +}; + +} //namespace details + +/*! +\brief Strategy functor for distance box to box calculation +\ingroup strategies +\details Class which calculates the distance of a box to a box, for +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 = cross_track<CalculationType> +> +class cross_track_box_box +{ +public: + template <typename Box1, typename Box2> + struct return_type + : services::return_type<Strategy, + typename point_type<Box1>::type, + typename point_type<Box2>::type> + {}; + + typedef typename Strategy::radius_type radius_type; + + inline cross_track_box_box() + {} + + explicit inline cross_track_box_box(typename Strategy::radius_type const& r) + : m_ps_strategy(r) + {} + + inline cross_track_box_box(Strategy const& s) + : m_ps_strategy(s) + {} + + + // It might be useful in the future + // to overload constructor with strategy info. + // crosstrack(...) {} + + template <typename Box1, typename Box2> + inline typename return_type<Box1, Box2>::type + apply(Box1 const& box1, Box2 const& box2) const + { +#if !defined(BOOST_MSVC) + BOOST_CONCEPT_ASSERT + ( + (concepts::PointSegmentDistanceStrategy + < + Strategy, + typename point_type<Box1>::type, + typename point_type<Box2>::type + >) + ); +#endif + typedef typename return_type<Box1, Box2>::type return_type; + return details::cross_track_box_box_generic + <return_type>::apply(box1, box2, m_ps_strategy); + } + + inline typename Strategy::radius_type radius() const + { + return m_ps_strategy.radius(); + } + +private: + Strategy m_ps_strategy; +}; + + + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +template <typename CalculationType, typename Strategy> +struct tag<cross_track_box_box<CalculationType, Strategy> > +{ + typedef strategy_tag_distance_box_box type; +}; + + +template <typename CalculationType, typename Strategy, typename Box1, typename Box2> +struct return_type<cross_track_box_box<CalculationType, Strategy>, Box1, Box2> + : cross_track_box_box + < + CalculationType, Strategy + >::template return_type<Box1, Box2> +{}; + + +template <typename CalculationType, typename Strategy> +struct comparable_type<cross_track_box_box<CalculationType, Strategy> > +{ + typedef cross_track_box_box + < + CalculationType, typename comparable_type<Strategy>::type + > type; +}; + + +template <typename CalculationType, typename Strategy> +struct get_comparable<cross_track_box_box<CalculationType, Strategy> > +{ + typedef cross_track_box_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 comparable_type(strategy.radius()); + } +}; + + +template <typename CalculationType, typename Strategy, typename Box1, typename Box2> +struct result_from_distance + < + cross_track_box_box<CalculationType, Strategy>, Box1, Box2 + > +{ +private: + typedef cross_track_box_box<CalculationType, Strategy> this_strategy; + + typedef typename this_strategy::template return_type + < + Box1, Box2 + >::type return_type; + +public: + template <typename T> + static inline return_type apply(this_strategy const& strategy, + T const& distance) + { + Strategy s(strategy.radius()); + + return result_from_distance + < + Strategy, + typename point_type<Box1>::type, + typename point_type<Box2>::type + >::apply(s, distance); + } +}; + + +// define cross_track_box_box<default_point_segment_strategy> as +// default box-box strategy for the spherical equatorial coordinate system +template <typename Box1, typename Box2, typename Strategy> +struct default_strategy + < + box_tag, box_tag, Box1, Box2, + spherical_equatorial_tag, spherical_equatorial_tag, + Strategy + > +{ + typedef cross_track_box_box + < + void, + typename boost::mpl::if_ + < + boost::is_void<Strategy>, + typename default_strategy + < + point_tag, segment_tag, + typename point_type<Box1>::type, typename point_type<Box2>::type, + spherical_equatorial_tag, spherical_equatorial_tag + >::type, + Strategy + >::type + > type; +}; + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::distance + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_CROSS_TRACK_BOX_BOX_HPP 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 ee805c36d1..9c8e7f1a3e 100644 --- a/boost/geometry/strategies/spherical/distance_cross_track_point_box.hpp +++ b/boost/geometry/strategies/spherical/distance_cross_track_point_box.hpp @@ -4,11 +4,12 @@ // 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, 2015. -// Modifications copyright (c) 2014-2015, Oracle and/or its affiliates. +// This file was modified by Oracle on 2014-2017. +// Modifications copyright (c) 2014-2017, Oracle and/or its affiliates. -// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle +// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle // 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 @@ -42,92 +43,46 @@ 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 = cross_track<CalculationType> -> -class cross_track_point_box +namespace details { -public: - template <typename Point, typename Box> - struct return_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_ps_strategy(r) - {} - - inline cross_track_point_box(Strategy const& s) - : m_ps_strategy(s) - {} +template <typename ReturnType> +class cross_track_point_box_generic +{ +public : - // 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 + template + < + typename Point, + typename Box, + typename Strategy + > + ReturnType static inline apply (Point const& point, + Box const& box, + Strategy ps_strategy) { -#if !defined(BOOST_MSVC) - BOOST_CONCEPT_ASSERT - ( - (concepts::PointSegmentDistanceStrategy - < - 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_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 - 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); + ReturnType const plon = geometry::get_as_radian<0>(point); + ReturnType 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); + ReturnType const lon_min = geometry::get_as_radian<0>(bottom_left); + ReturnType const lat_min = geometry::get_as_radian<1>(bottom_left); + ReturnType const lon_max = geometry::get_as_radian<0>(top_right); + ReturnType 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>(); + ReturnType const pi = math::pi<ReturnType>(); + ReturnType const two_pi = math::two_pi<ReturnType>(); + + typedef typename point_type<Box>::type box_point_type; // First check if the point is within the band defined by the // minimum and maximum longitude of the box; if yes, determine @@ -142,22 +97,22 @@ public: { if (plat > lat_max) { - return services::result_from_distance - < - Strategy, Point, box_point_type - >::apply(m_ps_strategy, radius() * (plat - lat_max)); + return geometry::strategy::distance::services::result_from_distance + < + Strategy, Point, box_point_type + >::apply(ps_strategy, ps_strategy.get_distance_strategy().meridian(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)); + return geometry::strategy::distance::services::result_from_distance + < + Strategy, Point, box_point_type + >::apply(ps_strategy, ps_strategy.get_distance_strategy().meridian(lat_min, plat)); } else { BOOST_GEOMETRY_ASSERT(plat >= lat_min && plat <= lat_max); - return return_type(0); + return ReturnType(0); } } @@ -169,14 +124,14 @@ public: // (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; + ReturnType const two = 2.0; bool use_left_segment; if (lon_max > pi) { // the box crosses the antimeridian // midway longitude = lon_min - (lon_min + (lon_max - 2 * pi)) / 2; - return_type const lon_midway = (lon_min - lon_max) / two + pi; + ReturnType const lon_midway = (lon_min - lon_max) / two + pi; BOOST_GEOMETRY_ASSERT(lon_midway >= -pi && lon_midway <= pi); use_left_segment = plon > lon_midway; @@ -185,8 +140,8 @@ public: { // the box does not cross the antimeridian - return_type const lon_sum = lon_min + lon_max; - if (math::equals(lon_sum, return_type(0))) + ReturnType const lon_sum = lon_min + lon_max; + if (math::equals(lon_sum, ReturnType(0))) { // special case: the box is symmetric with respect to // the prime meridian; the midway meridian is the antimeridian @@ -196,7 +151,7 @@ public: else { // midway long. = lon_min - (2 * pi - (lon_max - lon_min)) / 2; - return_type lon_midway = (lon_min + lon_max) / two - pi; + ReturnType lon_midway = (lon_min + lon_max) / two - pi; // normalize the midway longitude if (lon_midway > pi) @@ -212,14 +167,79 @@ public: // 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); + ? (plon < lon_min && plon >= lon_midway) + : (plon <= lon_max || plon > lon_midway); } } return use_left_segment - ? m_ps_strategy.apply(point, bottom_left, top_left) - : m_ps_strategy.apply(point, bottom_right, top_right); + ? ps_strategy.apply(point, bottom_left, top_left) + : ps_strategy.apply(point, bottom_right, top_right); + } +}; + +} //namespace details + +/*! +\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 = cross_track<CalculationType> +> +class cross_track_point_box +{ +public: + template <typename Point, typename Box> + struct return_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_ps_strategy(r) + {} + + inline cross_track_point_box(Strategy const& 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 + ( + (concepts::PointSegmentDistanceStrategy + < + Strategy, Point, typename point_type<Box>::type + >) + ); +#endif + typedef typename return_type<Point, Box>::type return_type; + return details::cross_track_point_box_generic + <return_type>::apply(point, box, m_ps_strategy); } inline typename Strategy::radius_type radius() const diff --git a/boost/geometry/strategies/spherical/distance_haversine.hpp b/boost/geometry/strategies/spherical/distance_haversine.hpp index 8b32056f3e..1a4cbdf36b 100644 --- a/boost/geometry/strategies/spherical/distance_haversine.hpp +++ b/boost/geometry/strategies/spherical/distance_haversine.hpp @@ -2,6 +2,12 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. +// This file was modified by Oracle on 2017, 2018. +// Modifications copyright (c) 2017-2018, Oracle and/or its affiliates. + +// Contributed and/or modified by Vissarion Fysikopoulos, 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 // http://www.boost.org/LICENSE_1_0.txt) @@ -10,16 +16,18 @@ #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_DISTANCE_HAVERSINE_HPP -#include <boost/geometry/core/cs.hpp> #include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/cs.hpp> #include <boost/geometry/core/radian_access.hpp> -#include <boost/geometry/util/math.hpp> -#include <boost/geometry/util/select_calculation_type.hpp> -#include <boost/geometry/util/promote_floating_point.hpp> +#include <boost/geometry/srs/sphere.hpp> #include <boost/geometry/strategies/distance.hpp> +#include <boost/geometry/strategies/spherical/get_radius.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 @@ -40,7 +48,7 @@ namespace comparable // - applying asin (which is strictly (monotone) increasing) template < - typename RadiusType, + typename RadiusTypeOrSphere = double, typename CalculationType = void > class haversine @@ -59,10 +67,21 @@ public : > {}; - typedef RadiusType radius_type; + typedef typename strategy_detail::get_radius + < + RadiusTypeOrSphere + >::type radius_type; + + inline haversine() + : m_radius(1.0) + {} - explicit inline haversine(RadiusType const& r = 1.0) - : m_radius(r) + template <typename RadiusOrSphere> + explicit inline haversine(RadiusOrSphere const& radius_or_sphere) + : m_radius(strategy_detail::get_radius + < + RadiusOrSphere + >::apply(radius_or_sphere)) {} template <typename Point1, typename Point2> @@ -75,7 +94,13 @@ public : ); } - inline RadiusType radius() const + template <typename T1, typename T2> + inline radius_type meridian(T1 lat1, T2 lat2) const + { + return m_radius * (lat1 - lat2); + } + + inline radius_type radius() const { return m_radius; } @@ -90,7 +115,7 @@ private : + cos(lat1) * cos(lat2) * math::hav(lon2 - lon1); } - RadiusType m_radius; + radius_type m_radius; }; @@ -101,7 +126,7 @@ private : \brief Distance calculation for spherical coordinates on a perfect sphere using haversine \ingroup strategies -\tparam RadiusType \tparam_radius +\tparam RadiusTypeOrSphere \tparam_radius_or_sphere \tparam CalculationType \tparam_calculation \author Adapted from: http://williams.best.vwh.net/avform.htm \see http://en.wikipedia.org/wiki/Great-circle_distance @@ -112,22 +137,19 @@ A mathematically equivalent formula, which is less subject to rounding error for short distances is: d=2*asin(sqrt((sin((lat1-lat2) / 2))^2 + cos(lat1)*cos(lat2)*(sin((lon1-lon2) / 2))^2)) - - \qbk{ [heading See also] [link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)] } - */ template < - typename RadiusType, + typename RadiusTypeOrSphere = double, typename CalculationType = void > class haversine { - typedef comparable::haversine<RadiusType, CalculationType> comparable_type; + typedef comparable::haversine<RadiusTypeOrSphere, CalculationType> comparable_type; public : template <typename Point1, typename Point2> @@ -135,14 +157,28 @@ public : : services::return_type<comparable_type, Point1, Point2> {}; - typedef RadiusType radius_type; + typedef typename strategy_detail::get_radius + < + RadiusTypeOrSphere + >::type radius_type; + + /*! + \brief Default constructor, radius set to 1.0 for the unit sphere + */ + inline haversine() + : m_radius(1.0) + {} /*! \brief Constructor - \param radius radius of the sphere, defaults to 1.0 for the unit sphere + \param radius_or_sphere radius of the sphere or sphere model */ - inline haversine(RadiusType const& radius = 1.0) - : m_radius(radius) + template <typename RadiusOrSphere> + explicit inline haversine(RadiusOrSphere const& radius_or_sphere) + : m_radius(strategy_detail::get_radius + < + RadiusOrSphere + >::apply(radius_or_sphere)) {} /*! @@ -162,16 +198,29 @@ public : } /*! + \brief meridian distance calculation + \return the calculated distance (including multiplying with radius) + \param p1 first point + \param p2 second point + */ + + template <typename T1, typename T2> + inline radius_type meridian(T1 lat1, T2 lat2) const + { + return m_radius * (lat1 - lat2); + } + + /*! \brief access to radius value \return the radius */ - inline RadiusType radius() const + inline radius_type radius() const { return m_radius; } private : - RadiusType m_radius; + radius_type m_radius; }; diff --git a/boost/geometry/strategies/spherical/get_radius.hpp b/boost/geometry/strategies/spherical/get_radius.hpp new file mode 100644 index 0000000000..411642b13a --- /dev/null +++ b/boost/geometry/strategies/spherical/get_radius.hpp @@ -0,0 +1,81 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland. + +// Copyright (c) 2016-2018 Oracle and/or its affiliates. +// Contributed and/or modified by Vissarion Fisikopoulos, 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 +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_GET_RADIUS_HPP +#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_GET_RADIUS_HPP + + +#include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/core/radius.hpp> +#include <boost/geometry/core/tag.hpp> +#include <boost/geometry/core/tags.hpp> +#include <boost/geometry/util/select_most_precise.hpp> + + +namespace boost { namespace geometry +{ + + +#ifndef DOXYGEN_NO_DETAIL +namespace strategy_detail +{ + +template +< + typename RadiusTypeOrSphere, + typename Tag = typename tag<RadiusTypeOrSphere>::type +> +struct get_radius +{ + typedef typename geometry::radius_type<RadiusTypeOrSphere>::type type; + static type apply(RadiusTypeOrSphere const& sphere) + { + return geometry::get_radius<0>(sphere); + } +}; + +template <typename RadiusTypeOrSphere> +struct get_radius<RadiusTypeOrSphere, void> +{ + typedef RadiusTypeOrSphere type; + static type apply(RadiusTypeOrSphere const& radius) + { + return radius; + } +}; + +// For backward compatibility +template <typename Point> +struct get_radius<Point, point_tag> +{ + typedef typename select_most_precise + < + typename coordinate_type<Point>::type, + double + >::type type; + + template <typename RadiusOrSphere> + static typename get_radius<RadiusOrSphere>::type + apply(RadiusOrSphere const& radius_or_sphere) + { + return get_radius<RadiusOrSphere>::apply(radius_or_sphere); + } +}; + + +} // namespace strategy_detail +#endif // DOXYGEN_NO_DETAIL + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_GET_RADIUS_HPP diff --git a/boost/geometry/strategies/spherical/intersection.hpp b/boost/geometry/strategies/spherical/intersection.hpp index b5ae878b85..1cc7b20764 100644 --- a/boost/geometry/strategies/spherical/intersection.hpp +++ b/boost/geometry/strategies/spherical/intersection.hpp @@ -1,5 +1,7 @@ // Boost.Geometry +// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland. + // Copyright (c) 2016-2017, Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -118,7 +120,7 @@ struct ecef_segments { typedef area::spherical < - typename point_type<Geometry>::type, + typename coordinate_type<Geometry>::type, CalculationType > type; }; @@ -161,38 +163,12 @@ struct ecef_segments template <typename CoordinateType, typename SegmentRatio, typename Vector3d> struct segment_intersection_info { - typedef typename select_most_precise - < - CoordinateType, double - >::type promoted_type; - segment_intersection_info(CalcPolicy const& calc) : calc_policy(calc) {} - promoted_type comparable_length_a() const - { - return robust_ra.denominator(); - } - - promoted_type comparable_length_b() const - { - return robust_rb.denominator(); - } - - template <typename Point, typename Segment1, typename Segment2> - void assign_a(Point& point, Segment1 const& a, Segment2 const& b) const - { - assign(point, a, b); - } - template <typename Point, typename Segment1, typename Segment2> - void assign_b(Point& point, Segment1 const& a, Segment2 const& b) const - { - assign(point, a, b); - } - template <typename Point, typename Segment1, typename Segment2> - void assign(Point& point, Segment1 const& a, Segment2 const& b) const + void calculate(Point& point, Segment1 const& a, Segment2 const& b) const { if (ip_flag == ipi_inters) { diff --git a/boost/geometry/strategies/spherical/point_in_poly_winding.hpp b/boost/geometry/strategies/spherical/point_in_poly_winding.hpp index 0f1a901d13..c942cbe460 100644 --- a/boost/geometry/strategies/spherical/point_in_poly_winding.hpp +++ b/boost/geometry/strategies/spherical/point_in_poly_winding.hpp @@ -1,7 +1,7 @@ // 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) 2013-2017 Adam Wulkiewicz, Lodz, Poland. // This file was modified by Oracle on 2013, 2014, 2016, 2017. // Modifications copyright (c) 2013-2017 Oracle and/or its affiliates. @@ -169,7 +169,7 @@ public: int side = 0; if (ci.count == 1 || ci.count == -1) { - side = side_equal(point, eq1 ? s1 : s2, ci, s1, s2); + side = side_equal(point, eq1 ? s1 : s2, ci); } else // count == 2 || count == -2 { @@ -431,8 +431,7 @@ private: // but the point is not aligned with a vertical segment inline int side_equal(Point const& point, PointOfSegment const& se, - count_info const& ci, - PointOfSegment const& s1, PointOfSegment const& s2) const + count_info const& ci) const { typedef typename coordinate_type<PointOfSegment>::type scoord_t; typedef typename coordinate_system<PointOfSegment>::type::units units_t; |