diff options
Diffstat (limited to 'boost/geometry/strategies')
33 files changed, 2401 insertions, 484 deletions
diff --git a/boost/geometry/strategies/area.hpp b/boost/geometry/strategies/area.hpp index e192d9b28b..866f37e846 100644 --- a/boost/geometry/strategies/area.hpp +++ b/boost/geometry/strategies/area.hpp @@ -3,6 +3,7 @@ // 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) 2017 Adam Wulkiewicz, Lodz, Poland. // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -14,35 +15,76 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_AREA_HPP #define BOOST_GEOMETRY_STRATEGIES_AREA_HPP + +#include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/util/select_most_precise.hpp> + #include <boost/mpl/assert.hpp> -#include <boost/geometry/strategies/tags.hpp> namespace boost { namespace geometry { -namespace strategy { namespace area { namespace services +namespace strategy { namespace area +{ + + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +// If user specified a CalculationType, use that type, whatever it is +// and whatever the Geometry is. +// Else, use Geometry's coordinate-type promoted to double if needed. +template +< + typename Geometry, + typename CalculationType +> +struct result_type +{ + typedef CalculationType type; +}; + +template +< + typename Geometry +> +struct result_type<Geometry, void> + : select_most_precise + < + typename coordinate_type<Geometry>::type, + double + > +{}; + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + + +namespace services { /*! \brief Traits class binding a default area strategy to a coordinate system \ingroup area \tparam Tag tag of coordinate system - \tparam PointOfSegment point-type */ -template <typename Tag, typename PointOfSegment> +template <typename Tag> struct default_strategy { BOOST_MPL_ASSERT_MSG ( - false, NOT_IMPLEMENTED_FOR_THIS_POINT_TYPE - , (types<PointOfSegment>) + false, NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM + , (types<Tag>) ); }; -}}} // namespace strategy::area::services +} // namespace services + +}} // namespace strategy::area }} // namespace boost::geometry diff --git a/boost/geometry/strategies/area_result.hpp b/boost/geometry/strategies/area_result.hpp new file mode 100644 index 0000000000..b8a6c4c85e --- /dev/null +++ b/boost/geometry/strategies/area_result.hpp @@ -0,0 +1,116 @@ +// Boost.Geometry + +// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland. + +// 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_AREA_RESULT_HPP +#define BOOST_GEOMETRY_STRATEGIES_AREA_RESULT_HPP + + +#include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/core/cs.hpp> + +#include <boost/geometry/strategies/area.hpp> +#include <boost/geometry/strategies/default_strategy.hpp> + +#include <boost/geometry/util/select_most_precise.hpp> +#include <boost/geometry/util/select_sequence_element.hpp> + + +#include <boost/mpl/if.hpp> +#include <boost/type_traits/is_same.hpp> +#include <boost/variant/variant_fwd.hpp> + + +namespace boost { namespace geometry +{ + + +/*! +\brief Meta-function defining return type of area function +\ingroup area +\note The return-type is defined by Geometry and Strategy + */ +template +< + typename Geometry, + typename Strategy = default_strategy +> +struct area_result + : Strategy::template result_type<Geometry> +{}; + +template <BOOST_VARIANT_ENUM_PARAMS(typename T), typename Strategy> +struct area_result<boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)>, Strategy> + : geometry::area_result + < + typename geometry::util::select_sequence_element + < + typename boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)>::types + >::type, + Strategy + > +{}; + +template <typename Geometry> +struct area_result<Geometry, default_strategy> + : geometry::area_result + < + Geometry, + typename geometry::strategy::area::services::default_strategy + < + typename cs_tag<Geometry>::type + >::type + > +{}; + +#ifndef DOXYGEN_NO_DETAIL +namespace detail { namespace area +{ + +template <typename Curr, typename Next> +struct pred_more_precise_default_area_result +{ + typedef typename geometry::area_result<Curr, default_strategy>::type curr_result_t; + typedef typename geometry::area_result<Next, default_strategy>::type next_result_t; + + typedef typename boost::mpl::if_c + < + boost::is_same + < + curr_result_t, + typename geometry::select_most_precise + < + curr_result_t, + next_result_t + >::type + >::value, + Curr, + Next + >::type type; +}; + +}} // namespace detail::area +#endif //DOXYGEN_NO_DETAIL + +template <BOOST_VARIANT_ENUM_PARAMS(typename T)> +struct area_result<boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)>, default_strategy> + : geometry::area_result + < + typename geometry::util::select_sequence_element + < + typename boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)>::types, + geometry::detail::area::pred_more_precise_default_area_result + >::type, + default_strategy + > +{}; + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_AREA_RESULT_HPP diff --git a/boost/geometry/strategies/cartesian/area.hpp b/boost/geometry/strategies/cartesian/area.hpp new file mode 100644 index 0000000000..27708424c2 --- /dev/null +++ b/boost/geometry/strategies/cartesian/area.hpp @@ -0,0 +1,146 @@ +// 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) 2017 Adam Wulkiewicz, Lodz, Poland. + +// This file was modified by Oracle on 2016, 2017. +// Modifications copyright (c) 2016-2017, 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. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AREA_HPP +#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AREA_HPP + + +#include <boost/mpl/if.hpp> + +//#include <boost/geometry/arithmetic/determinant.hpp> +#include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/core/coordinate_dimension.hpp> +#include <boost/geometry/strategies/area.hpp> +#include <boost/geometry/util/select_most_precise.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace area +{ + +/*! +\brief Cartesian area calculation +\ingroup strategies +\details Calculates cartesian area using the trapezoidal rule +\tparam CalculationType \tparam_calculation + +\qbk{ +[heading See also] +[link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)] +} + +*/ +template +< + typename CalculationType = void +> +class cartesian +{ +public : + template <typename Geometry> + struct result_type + : strategy::area::detail::result_type + < + Geometry, + CalculationType + > + {}; + + template <typename Geometry> + class state + { + friend class cartesian; + + typedef typename result_type<Geometry>::type return_type; + + public: + inline state() + : sum(0) + { + // Strategy supports only 2D areas + assert_dimension<Geometry, 2>(); + } + + private: + inline return_type area() const + { + return_type const two = 2; + return sum / two; + } + + return_type sum; + }; + + template <typename PointOfSegment, typename Geometry> + static inline void apply(PointOfSegment const& p1, + PointOfSegment const& p2, + state<Geometry>& st) + { + typedef typename state<Geometry>::return_type return_type; + + // Below formulas are equivalent, however the two lower ones + // suffer less from accuracy loss for great values of coordinates. + // See: https://svn.boost.org/trac/boost/ticket/11928 + + // SUM += x2 * y1 - x1 * y2; + // state.sum += detail::determinant<return_type>(p2, p1); + + // SUM += (x2 - x1) * (y2 + y1) + //state.sum += (return_type(get<0>(p2)) - return_type(get<0>(p1))) + // * (return_type(get<1>(p2)) + return_type(get<1>(p1))); + + // SUM += (x1 + x2) * (y1 - y2) + st.sum += (return_type(get<0>(p1)) + return_type(get<0>(p2))) + * (return_type(get<1>(p1)) - return_type(get<1>(p2))); + } + + template <typename Geometry> + static inline typename result_type<Geometry>::type + result(state<Geometry>& st) + { + return st.area(); + } + +}; + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +namespace services +{ + template <> + struct default_strategy<cartesian_tag> + { + typedef strategy::area::cartesian<> type; + }; + +} // namespace services + +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::area + + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AREA_HPP diff --git a/boost/geometry/strategies/cartesian/area_surveyor.hpp b/boost/geometry/strategies/cartesian/area_surveyor.hpp index b3f19b1b1e..91df31745a 100644 --- a/boost/geometry/strategies/cartesian/area_surveyor.hpp +++ b/boost/geometry/strategies/cartesian/area_surveyor.hpp @@ -3,6 +3,7 @@ // 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) 2017 Adam Wulkiewicz, Lodz, Poland. // This file was modified by Oracle on 2016, 2017. // Modifications copyright (c) 2016-2017, Oracle and/or its affiliates. @@ -20,133 +21,9 @@ #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AREA_SURVEYOR_HPP -#include <boost/mpl/if.hpp> - -//#include <boost/geometry/arithmetic/determinant.hpp> -#include <boost/geometry/core/coordinate_type.hpp> -#include <boost/geometry/core/coordinate_dimension.hpp> -#include <boost/geometry/strategies/area.hpp> -#include <boost/geometry/util/select_most_precise.hpp> - - -namespace boost { namespace geometry -{ - -namespace strategy { namespace area -{ - -/*! -\brief Area calculation for cartesian points -\ingroup strategies -\details Calculates area using the Surveyor's formula, a well-known - triangulation algorithm -\tparam PointOfSegment \tparam_segment_point -\tparam CalculationType \tparam_calculation - -\qbk{ -[heading See also] -[link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)] -} - -*/ -template -< - typename PointOfSegment, - typename CalculationType = void -> -class surveyor -{ -public : - // If user specified a calculation type, use that type, - // whatever it is and whatever the point-type is. - // Else, use the pointtype, but at least double - typedef typename - boost::mpl::if_c - < - boost::is_void<CalculationType>::type::value, - typename select_most_precise - < - typename coordinate_type<PointOfSegment>::type, - double - >::type, - CalculationType - >::type return_type; - - -private : - - class summation - { - friend class surveyor; - - return_type sum; - public : - - inline summation() : sum(return_type()) - { - // Strategy supports only 2D areas - assert_dimension<PointOfSegment, 2>(); - } - inline return_type area() const - { - return_type result = sum; - return_type const two = 2; - result /= two; - return result; - } - }; - -public : - typedef summation state_type; - typedef PointOfSegment segment_point_type; - - static inline void apply(PointOfSegment const& p1, - PointOfSegment const& p2, - summation& state) - { - // Below formulas are equivalent, however the two lower ones - // suffer less from accuracy loss for great values of coordinates. - // See: https://svn.boost.org/trac/boost/ticket/11928 - - // SUM += x2 * y1 - x1 * y2; - // state.sum += detail::determinant<return_type>(p2, p1); - - // SUM += (x2 - x1) * (y2 + y1) - //state.sum += (return_type(get<0>(p2)) - return_type(get<0>(p1))) - // * (return_type(get<1>(p2)) + return_type(get<1>(p1))); - - // SUM += (x1 + x2) * (y1 - y2) - state.sum += (return_type(get<0>(p1)) + return_type(get<0>(p2))) - * (return_type(get<1>(p1)) - return_type(get<1>(p2))); - } - - static inline return_type result(summation const& state) - { - return state.area(); - } - -}; - -#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - -namespace services -{ - template <typename Point> - struct default_strategy<cartesian_tag, Point> - { - typedef strategy::area::surveyor<Point> type; - }; - -} // namespace services - -#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - - -}} // namespace strategy::area - - - -}} // namespace boost::geometry +// keep this file for now, for backward compatibility +// functionality-wise, make it equivalent to boost/geometry/strategies/cartesian/area.hpp +#include <boost/geometry/strategies/cartesian/area.hpp> #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AREA_SURVEYOR_HPP diff --git a/boost/geometry/strategies/cartesian/centroid_average.hpp b/boost/geometry/strategies/cartesian/centroid_average.hpp index c12f6e2024..32f551cf95 100644 --- a/boost/geometry/strategies/cartesian/centroid_average.hpp +++ b/boost/geometry/strategies/cartesian/centroid_average.hpp @@ -3,6 +3,7 @@ // 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) 2017 Adam Wulkiewicz, Lodz, Poland. // This file was modified by Oracle on 2015. // Modifications copyright (c) 2015 Oracle and/or its affiliates. @@ -23,6 +24,7 @@ #include <cstddef> #include <boost/geometry/algorithms/assign.hpp> +#include <boost/geometry/algorithms/detail/signed_size_type.hpp> #include <boost/geometry/arithmetic/arithmetic.hpp> #include <boost/geometry/core/coordinate_type.hpp> #include <boost/geometry/core/point_type.hpp> @@ -53,7 +55,7 @@ private : class sum { friend class average; - std::size_t count; + signed_size_type count; PointCentroid centroid; public : diff --git a/boost/geometry/strategies/cartesian/densify.hpp b/boost/geometry/strategies/cartesian/densify.hpp new file mode 100644 index 0000000000..23651637b9 --- /dev/null +++ b/boost/geometry/strategies/cartesian/densify.hpp @@ -0,0 +1,133 @@ +// 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_CARTESIAN_DENSIFY_HPP +#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_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/dot_product.hpp> +#include <boost/geometry/core/assert.hpp> +#include <boost/geometry/core/coordinate_dimension.hpp> +#include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/strategies/densify.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 cartesian segment. +\ingroup strategies +\tparam CalculationType \tparam_calculation + +\qbk{ +[heading See also] +[link geometry.reference.algorithms.densify.densify_4_with_strategy densify (with strategy)] +} + */ +template +< + typename CalculationType = void +> +class cartesian +{ +public: + template <typename Point, typename AssignPolicy, typename T> + static inline void apply(Point const& p0, Point const& p1, AssignPolicy & policy, T const& length_threshold) + { + 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; + + typedef model::point<calc_t, geometry::dimension<Point>::value, cs::cartesian> calc_point_t; + + calc_point_t cp0, cp1; + geometry::detail::conversion::convert_point_to_point(p0, cp0); + geometry::detail::conversion::convert_point_to_point(p1, cp1); + + // dir01 = xy1 - xy0 + calc_point_t dir01 = cp1; + geometry::subtract_point(dir01, cp0); + calc_t const dot01 = geometry::dot_product(dir01, dir01); + calc_t const len = math::sqrt(dot01); + + BOOST_GEOMETRY_ASSERT(length_threshold > T(0)); + + signed_size_type n = signed_size_type(len / length_threshold); + if (n <= 0) + { + return; + } + + // NOTE: Normalization will not work for integral coordinates + // normalize + //geometry::divide_value(dir01, len); + + calc_t step = len / (n + 1); + + calc_t d = step; + for (signed_size_type i = 0 ; i < n ; ++i, d += step) + { + // pd = xy0 + d * dir01 + calc_point_t pd = dir01; + + // without normalization + geometry::multiply_value(pd, calc_t(i + 1)); + geometry::divide_value(pd, calc_t(n + 1)); + // with normalization + //geometry::multiply_value(pd, d); + + geometry::add_point(pd, cp0); + + // NOTE: Only needed if types calc_point_t and out_point_t are different + // otherwise pd could simply be passed into policy + out_point_t p; + assert_dimension_equal<calc_point_t, out_point_t>(); + geometry::detail::conversion::convert_point_to_point(pd, p); + + policy.apply(p); + } + } +}; + + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +template <> +struct default_strategy<cartesian_tag> +{ + typedef strategy::densify::cartesian<> type; +}; + + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::densify + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_DENSIFY_HPP diff --git a/boost/geometry/strategies/cartesian/intersection.hpp b/boost/geometry/strategies/cartesian/intersection.hpp index 233bb50b64..10fd73ba54 100644 --- a/boost/geometry/strategies/cartesian/intersection.hpp +++ b/boost/geometry/strategies/cartesian/intersection.hpp @@ -1,7 +1,7 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. -// Copyright (c) 2013-2014 Adam Wulkiewicz, Lodz, Poland. +// Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland. // This file was modified by Oracle on 2014, 2016, 2017. // Modifications copyright (c) 2014-2017, Oracle and/or its affiliates. @@ -33,7 +33,7 @@ #include <boost/geometry/util/promote_integral.hpp> #include <boost/geometry/util/select_calculation_type.hpp> -#include <boost/geometry/strategies/cartesian/area_surveyor.hpp> +#include <boost/geometry/strategies/cartesian/area.hpp> #include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp> #include <boost/geometry/strategies/cartesian/envelope_segment.hpp> #include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp> @@ -103,9 +103,8 @@ struct cartesian_segments template <typename Geometry> struct area_strategy { - typedef area::surveyor + typedef area::cartesian < - typename point_type<Geometry>::type, CalculationType > type; }; @@ -144,6 +143,7 @@ struct cartesian_segments template <typename CoordinateType, typename SegmentRatio> struct segment_intersection_info { + private : typedef typename select_most_precise < CoordinateType, double @@ -198,6 +198,45 @@ struct cartesian_segments >(numerator * dy_promoted / denominator)); } + public : + template <typename Point, typename Segment1, typename Segment2> + void calculate(Point& point, Segment1 const& a, Segment2 const& b) const + { + bool use_a = true; + + // Prefer one segment if one is on or near an endpoint + bool const a_near_end = robust_ra.near_end(); + bool const b_near_end = robust_rb.near_end(); + if (a_near_end && ! b_near_end) + { + use_a = true; + } + else if (b_near_end && ! a_near_end) + { + use_a = false; + } + else + { + // Prefer shorter segment + promoted_type const len_a = comparable_length_a(); + promoted_type const len_b = comparable_length_b(); + if (len_b < len_a) + { + use_a = false; + } + // else use_a is true but was already assigned like that + } + + if (use_a) + { + assign_a(point, a, b); + } + else + { + assign_b(point, a, b); + } + } + CoordinateType dx_a, dy_a; CoordinateType dx_b, dy_b; SegmentRatio robust_ra; diff --git a/boost/geometry/strategies/concepts/area_concept.hpp b/boost/geometry/strategies/concepts/area_concept.hpp index 4eec6d1fc3..66982e3f27 100644 --- a/boost/geometry/strategies/concepts/area_concept.hpp +++ b/boost/geometry/strategies/concepts/area_concept.hpp @@ -3,6 +3,7 @@ // 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) 2017 Adam Wulkiewicz, Lodz, Poland. // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -26,19 +27,16 @@ namespace boost { namespace geometry { namespace concepts \brief Checks strategy for area \ingroup area */ -template <typename Strategy> +template <typename Geometry, typename Strategy> class AreaStrategy { #ifndef DOXYGEN_NO_CONCEPT_MEMBERS - // 1) must define state_type, - typedef typename Strategy::state_type state_type; + // 1) must define state template, + typedef typename Strategy::template state<Geometry> state_type; - // 2) must define return_type, - typedef typename Strategy::return_type return_type; - - // 3) must define point_type, of polygon (segments) - typedef typename Strategy::segment_point_type spoint_type; + // 2) must define result_type template, + typedef typename Strategy::template result_type<Geometry>::type return_type; struct check_methods { @@ -47,11 +45,11 @@ class AreaStrategy Strategy const* str = 0; state_type *st = 0; - // 4) must implement a method apply with the following signature - spoint_type const* sp = 0; + // 3) must implement a method apply with the following signature + typename geometry::point_type<Geometry>::type const* sp = 0; str->apply(*sp, *sp, *st); - // 5) must implement a static method result with the following signature + // 4) must implement a static method result with the following signature return_type r = str->result(*st); boost::ignore_unused_variable_warning(r); diff --git a/boost/geometry/strategies/default_area_result.hpp b/boost/geometry/strategies/default_area_result.hpp index 8adfa5d6ea..65818ea91a 100644 --- a/boost/geometry/strategies/default_area_result.hpp +++ b/boost/geometry/strategies/default_area_result.hpp @@ -3,6 +3,7 @@ // 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) 2017 Adam Wulkiewicz, Lodz, Poland. // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -15,10 +16,7 @@ #define BOOST_GEOMETRY_STRATEGIES_DEFAULT_AREA_RESULT_HPP -#include <boost/geometry/core/cs.hpp> -#include <boost/geometry/core/coordinate_type.hpp> -#include <boost/geometry/strategies/area.hpp> -#include <boost/geometry/util/select_most_precise.hpp> +#include <boost/geometry/strategies/area_result.hpp> namespace boost { namespace geometry @@ -33,16 +31,8 @@ namespace boost { namespace geometry template <typename Geometry> struct default_area_result -{ - typedef typename point_type<Geometry>::type point_type; - typedef typename strategy::area::services::default_strategy - < - typename cs_tag<point_type>::type, - point_type - >::type strategy_type; - - typedef typename strategy_type::return_type type; -}; + : area_result<Geometry> +{}; }} // namespace boost::geometry diff --git a/boost/geometry/strategies/densify.hpp b/boost/geometry/strategies/densify.hpp new file mode 100644 index 0000000000..cde061af92 --- /dev/null +++ b/boost/geometry/strategies/densify.hpp @@ -0,0 +1,42 @@ +// Boost.Geometry + +// Copyright (c) 2017, 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_DENSIFY_HPP +#define BOOST_GEOMETRY_STRATEGIES_DENSIFY_HPP + + +#include <boost/mpl/assert.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace densify +{ + +namespace services +{ + +template <typename CSTag> +struct default_strategy +{ + BOOST_MPL_ASSERT_MSG + ( + false, NOT_IMPLEMENTED_FOR_THIS_CS + , (types<CSTag>) + ); +}; + +} // namespace services + +}} // namespace strategy::densify + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_DENSIFY_HPP diff --git a/boost/geometry/strategies/geographic/area.hpp b/boost/geometry/strategies/geographic/area.hpp index 40d6c8243e..ac7d933ce7 100644 --- a/boost/geometry/strategies/geographic/area.hpp +++ b/boost/geometry/strategies/geographic/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 @@ -12,7 +14,7 @@ #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP -#include <boost/geometry/core/srs.hpp> +#include <boost/geometry/srs/spheroid.hpp> #include <boost/geometry/formulas/area_formulas.hpp> #include <boost/geometry/formulas/authalic_radius_sqr.hpp> @@ -32,7 +34,6 @@ namespace strategy { namespace area \ingroup strategies \details Geographic area calculation by trapezoidal rule plus integral approximation that gives the ellipsoidal correction -\tparam PointOfSegment \tparam_segment_point \tparam FormulaPolicy Formula used to calculate azimuths \tparam SeriesOrder The order of approximation of the geodesic integral \tparam Spheroid The spheroid model @@ -43,12 +44,12 @@ namespace strategy { namespace area \qbk{ [heading See also] -[link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)] +\* [link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)] +\* [link geometry.reference.srs.srs_spheroid srs::spheroid] } */ template < - typename PointOfSegment, typename FormulaPolicy = strategy::andoyer, std::size_t SeriesOrder = strategy::default_order<FormulaPolicy>::value, typename Spheroid = srs::spheroid<double>, @@ -63,58 +64,67 @@ class geographic //Select default types in case they are not set - 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; + Geometry, + CalculationType + > + {}; protected : struct spheroid_constants { + typedef typename boost::mpl::if_c + < + boost::is_void<CalculationType>::value, + typename geometry::radius_type<Spheroid>::type, + CalculationType + >::type calc_t; + Spheroid m_spheroid; - CT const m_a2; // squared equatorial radius - CT const m_e2; // squared eccentricity - CT const m_ep2; // squared second eccentricity - CT const m_ep; // second eccentricity - CT const m_c2; // squared authalic radius + calc_t const m_a2; // squared equatorial radius + calc_t const m_e2; // squared eccentricity + calc_t const m_ep2; // squared second eccentricity + calc_t const m_ep; // second eccentricity + calc_t const m_c2; // squared authalic radius inline spheroid_constants(Spheroid const& spheroid) : m_spheroid(spheroid) , m_a2(math::sqr(get_radius<0>(spheroid))) - , m_e2(formula::eccentricity_sqr<CT>(spheroid)) - , m_ep2(m_e2 / (CT(1.0) - m_e2)) + , m_e2(formula::eccentricity_sqr<calc_t>(spheroid)) + , m_ep2(m_e2 / (calc_t(1.0) - m_e2)) , m_ep(math::sqrt(m_ep2)) , m_c2(formula_dispatch::authalic_radius_sqr < - CT, Spheroid, srs_spheroid_tag + calc_t, Spheroid, srs_spheroid_tag >::apply(m_a2, m_e2)) {} }; - struct area_sums +public: + template <typename Geometry> + class state { - CT m_excess_sum; - CT m_correction_sum; + friend class geographic; - // Keep track if encircles some pole - std::size_t m_crosses_prime_meridian; + typedef typename result_type<Geometry>::type return_type; - inline area_sums() + public: + inline state() : m_excess_sum(0) , m_correction_sum(0) , m_crosses_prime_meridian(0) {} - inline CT area(spheroid_constants spheroid_const) const + + private: + inline return_type area(spheroid_constants const& spheroid_const) const { - CT result; + return_type result; - CT sum = spheroid_const.m_c2 * m_excess_sum + return_type sum = spheroid_const.m_c2 * m_excess_sum + spheroid_const.m_e2 * spheroid_const.m_a2 * m_correction_sum; // If encircles some pole @@ -123,13 +133,13 @@ protected : std::size_t times_crosses_prime_meridian = 1 + (m_crosses_prime_meridian / 2); - result = CT(2.0) - * geometry::math::pi<CT>() + result = return_type(2.0) + * geometry::math::pi<return_type>() * spheroid_const.m_c2 - * CT(times_crosses_prime_meridian) + * return_type(times_crosses_prime_meridian) - geometry::math::abs(sum); - if (geometry::math::sign<CT>(sum) == 1) + if (geometry::math::sign<return_type>(sum) == 1) { result = - result; } @@ -142,48 +152,52 @@ protected : return result; } + + return_type m_excess_sum; + return_type m_correction_sum; + + // Keep track if encircles some pole + std::size_t m_crosses_prime_meridian; }; public : - typedef CT return_type; - typedef PointOfSegment segment_point_type; - typedef area_sums state_type; - explicit inline geographic(Spheroid const& spheroid = Spheroid()) : m_spheroid_constants(spheroid) {} + template <typename PointOfSegment, typename Geometry> inline void apply(PointOfSegment const& p1, PointOfSegment const& p2, - area_sums& state) const + state<Geometry>& st) const { - if (! geometry::math::equals(get<0>(p1), get<0>(p2))) { - typedef geometry::formula::area_formulas < - CT, SeriesOrder, ExpandEpsN + typename result_type<Geometry>::type, + SeriesOrder, ExpandEpsN > area_formulas; typename area_formulas::return_type_ellipsoidal result = area_formulas::template ellipsoidal<FormulaPolicy::template inverse> (p1, p2, m_spheroid_constants); - state.m_excess_sum += result.spherical_term; - state.m_correction_sum += result.ellipsoidal_term; + st.m_excess_sum += result.spherical_term; + st.m_correction_sum += result.ellipsoidal_term; // 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(area_sums const& state) const + template <typename Geometry> + inline typename result_type<Geometry>::type + result(state<Geometry> const& st) const { - return state.area(m_spheroid_constants); + return st.area(m_spheroid_constants); } private: @@ -197,10 +211,10 @@ namespace services { -template <typename Point> -struct default_strategy<geographic_tag, Point> +template <> +struct default_strategy<geographic_tag> { - typedef strategy::area::geographic<Point> type; + typedef strategy::area::geographic<> type; }; #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS diff --git a/boost/geometry/strategies/geographic/azimuth.hpp b/boost/geometry/strategies/geographic/azimuth.hpp index 79c49750fb..b918caccea 100644 --- a/boost/geometry/strategies/geographic/azimuth.hpp +++ b/boost/geometry/strategies/geographic/azimuth.hpp @@ -12,7 +12,7 @@ #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AZIMUTH_HPP -#include <boost/geometry/core/srs.hpp> +#include <boost/geometry/srs/spheroid.hpp> #include <boost/geometry/strategies/azimuth.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp> diff --git a/boost/geometry/strategies/geographic/densify.hpp b/boost/geometry/strategies/geographic/densify.hpp new file mode 100644 index 0000000000..a31ba72200 --- /dev/null +++ b/boost/geometry/strategies/geographic/densify.hpp @@ -0,0 +1,136 @@ +// 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_GEOGRAPHIC_DENSIFY_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DENSIFY_HPP + + +#include <boost/geometry/algorithms/detail/convert_point_to_point.hpp> +#include <boost/geometry/algorithms/detail/signed_size_type.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/srs/spheroid.hpp> +#include <boost/geometry/strategies/densify.hpp> +#include <boost/geometry/strategies/geographic/parameters.hpp> +#include <boost/geometry/util/select_most_precise.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace densify +{ + + +/*! +\brief Densification of geographic segment. +\ingroup strategies +\tparam FormulaPolicy The geodesic formulas used internally. +\tparam Spheroid The spheroid model. +\tparam CalculationType \tparam_calculation + +\qbk{ +[heading See also] +\* [link geometry.reference.algorithms.densify.densify_4_with_strategy densify (with strategy)] +\* [link geometry.reference.srs.srs_spheroid srs::spheroid] +} + */ +template +< + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +class geographic +{ +public: + geographic() + : m_spheroid() + {} + + explicit geographic(Spheroid const& spheroid) + : m_spheroid(spheroid) + {} + + 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; + + typedef typename FormulaPolicy::template direct<calc_t, true, false, false, false> direct_t; + typedef typename FormulaPolicy::template inverse<calc_t, true, true, false, false, false> inverse_t; + + typename inverse_t::result_type + inv_r = inverse_t::apply(get_as_radian<0>(p0), get_as_radian<1>(p0), + get_as_radian<0>(p1), get_as_radian<1>(p1), + m_spheroid); + + BOOST_GEOMETRY_ASSERT(length_threshold > T(0)); + + signed_size_type n = signed_size_type(inv_r.distance / length_threshold); + if (n <= 0) + return; + + calc_t step = inv_r.distance / (n + 1); + + calc_t current = step; + for (signed_size_type i = 0 ; i < n ; ++i, current += step) + { + typename direct_t::result_type + dir_r = direct_t::apply(get_as_radian<0>(p0), get_as_radian<1>(p0), + current, inv_r.azimuth, + m_spheroid); + + out_point_t p; + set_from_radian<0>(p, dir_r.lon2); + set_from_radian<1>(p, dir_r.lat2); + geometry::detail::conversion::point_to_point + < + Point, out_point_t, + 2, dimension<out_point_t>::value + >::apply(p0, p); + + policy.apply(p); + } + } + +private: + Spheroid m_spheroid; +}; + + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +template <> +struct default_strategy<geographic_tag> +{ + typedef strategy::densify::geographic<> 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/geographic/distance.hpp b/boost/geometry/strategies/geographic/distance.hpp index 98ee46a369..01f766c10e 100644 --- a/boost/geometry/strategies/geographic/distance.hpp +++ b/boost/geometry/strategies/geographic/distance.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2016 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014-2017. -// Modifications copyright (c) 2014-2017 Oracle and/or its affiliates. +// This file was modified by Oracle on 2014-2018. +// Modifications copyright (c) 2014-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 @@ -20,12 +20,13 @@ #include <boost/geometry/core/coordinate_type.hpp> #include <boost/geometry/core/radian_access.hpp> #include <boost/geometry/core/radius.hpp> -#include <boost/geometry/core/srs.hpp> #include <boost/geometry/formulas/andoyer_inverse.hpp> #include <boost/geometry/formulas/elliptic_arc_length.hpp> #include <boost/geometry/formulas/flattening.hpp> +#include <boost/geometry/srs/spheroid.hpp> + #include <boost/geometry/strategies/distance.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp> @@ -34,6 +35,7 @@ #include <boost/geometry/util/promote_floating_point.hpp> #include <boost/geometry/util/select_calculation_type.hpp> +#include <boost/geometry/geometries/point_xy.hpp> namespace boost { namespace geometry { @@ -41,6 +43,19 @@ namespace boost { namespace geometry namespace strategy { namespace distance { +/*! +\brief Distance calculation for geographic coordinates on a spheroid +\ingroup strategies +\tparam FormulaPolicy Formula used to calculate azimuths +\tparam Spheroid The spheroid model +\tparam CalculationType \tparam_calculation + +\qbk{ +[heading See also] +\* [link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)] +\* [link geometry.reference.srs.srs_spheroid srs::spheroid] +} +*/ template < typename FormulaPolicy = strategy::andoyer, @@ -110,6 +125,19 @@ public : return apply(lon1, lat1, lon2, lat2, m_spheroid); } + // points on a meridian not crossing poles + template <typename CT> + inline CT meridian(CT lat1, CT lat2) const + { + typedef typename formula::elliptic_arc_length + < + CT, strategy::default_order<FormulaPolicy>::value + > elliptic_arc_length; + + return elliptic_arc_length::meridian_not_crossing_pole_dist(lat1, lat2, + m_spheroid); + } + inline Spheroid const& model() const { return m_spheroid; diff --git a/boost/geometry/strategies/geographic/distance_cross_track.hpp b/boost/geometry/strategies/geographic/distance_cross_track.hpp index 799208a96d..be930a3fd4 100644 --- a/boost/geometry/strategies/geographic/distance_cross_track.hpp +++ b/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -40,7 +40,7 @@ #include <boost/geometry/formulas/mean_radius.hpp> #ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK -# include <boost/geometry/io/dsv/write.hpp> +#include <boost/geometry/io/dsv/write.hpp> #endif #ifndef BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS @@ -165,19 +165,29 @@ private : } template <typename CT> - CT static inline normalize(CT g4) + CT static inline normalize(CT g4, CT& der) { CT const pi = math::pi<CT>(); - if (g4 < 0 && g4 < -pi)//close to -270 + if (g4 < -1.25*pi)//close to -270 { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "g4=" << g4 << ", close to -270" << std::endl; +#endif return g4 + 1.5 * pi; } - else if (g4 > 0 && g4 > pi)//close to 270 + else if (g4 > 1.25*pi)//close to 270 { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "g4=" << g4 << ", close to 270" << std::endl; +#endif return - g4 + 1.5 * pi; } - else if (g4 < 0 && g4 > -pi)//close to -90 + else if (g4 < 0 && g4 > -0.75*pi)//close to -90 { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "g4=" << g4 << ", close to -90" << std::endl; +#endif + der = -der; return -g4 - pi/2; } return g4 - pi/2; @@ -190,8 +200,8 @@ private : CT lon3, CT lat3, //query point p3 Spheroid const& spheroid) { - typedef typename FormulaPolicy::template inverse<CT, true, false, false, true, true> - inverse_distance_quantities_type; + typedef typename FormulaPolicy::template inverse<CT, true, true, false, true, true> + inverse_distance_azimuth_quantities_type; typedef typename FormulaPolicy::template inverse<CT, false, true, false, false, false> inverse_azimuth_type; typedef typename FormulaPolicy::template inverse<CT, false, true, true, false, false> @@ -204,7 +214,7 @@ private : result_distance_point_segment<CT> result; // Constants - CT const f = geometry::formula::flattening<CT>(spheroid); + //CT const f = geometry::formula::flattening<CT>(spheroid); CT const pi = math::pi<CT>(); CT const half_pi = pi / CT(2); CT const c0 = CT(0); @@ -227,11 +237,28 @@ private : //Note: antipodal points on equator does not define segment on equator //but pass by the pole CT diff = geometry::math::longitude_distance_signed<geometry::radian>(lon1, lon2); - if (math::equals(lat1, c0) && math::equals(lat2, c0) - && !math::equals(math::abs(diff), pi)) + + typedef typename formula::elliptic_arc_length<CT> elliptic_arc_length; + + bool meridian_not_crossing_pole = + elliptic_arc_length::meridian_not_crossing_pole(lat1, lat2, diff); + + bool meridian_crossing_pole = + elliptic_arc_length::meridian_crossing_pole(diff); + + //bool meridian_crossing_pole = math::equals(math::abs(diff), pi); + //bool meridian_not_crossing_pole = math::equals(math::abs(diff), c0); + + if (math::equals(lat1, c0) && math::equals(lat2, c0) && !meridian_crossing_pole) { #ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK std::cout << "Equatorial segment" << std::endl; + std::cout << "segment=(" << lon1 * math::r2d<CT>(); + std::cout << "," << lat1 * math::r2d<CT>(); + std::cout << "),(" << lon2 * math::r2d<CT>(); + std::cout << "," << lat2 * math::r2d<CT>(); + std::cout << ")\np=(" << lon3 * math::r2d<CT>(); + std::cout << "," << lat3 * math::r2d<CT>() << ")\n"; #endif if (lon3 <= lon1) { @@ -244,7 +271,12 @@ private : return non_iterative_case(lon3, lat1, lon3, lat3, spheroid); } - if (math::equals(math::abs(diff), pi)) + if ( (meridian_not_crossing_pole || meridian_crossing_pole ) && lat1 > lat2) + { + std::swap(lat1,lat2); + } + + if (meridian_crossing_pole) { #ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK std::cout << "Meridian segment" << std::endl; @@ -349,12 +381,9 @@ private : geometry::cs::spherical_equatorial<geometry::radian> > point; - CT bet1 = atan((1 - f) * tan(lon1)); - CT bet2 = atan((1 - f) * tan(lon2)); - CT bet3 = atan((1 - f) * tan(lon3)); - point p1 = point(bet1, lat1); - point p2 = point(bet2, lat2); - point p3 = point(bet3, lat3); + point p1 = point(lon1, lat1); + point p2 = point(lon2, lat2); + point p3 = point(lon3, lat3); geometry::strategy::distance::cross_track<CT> cross_track(earth_radius); CT s34 = cross_track.apply(p3, p1, p2); @@ -390,22 +419,28 @@ private : CT a4 = inverse_azimuth_type::apply(res14.lon2, res14.lat2, lon2, lat2, spheroid).azimuth; - res34 = inverse_distance_quantities_type::apply(res14.lon2, res14.lat2, - lon3, lat3, spheroid); + res34 = inverse_distance_azimuth_quantities_type::apply(res14.lon2, res14.lat2, + lon3, lat3, spheroid); g4 = res34.azimuth - a4; - delta_g4 = normalize(g4); + CT M43 = res34.geodesic_scale; // cos(s14/earth_radius) is the spherical limit CT m34 = res34.reduced_length; CT der = (M43 / m34) * sin(g4); + + // normalize (g4 - pi/2) + delta_g4 = normalize(g4, der); + s14 = s14 - delta_g4 / der; #ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK std::cout << "p4=" << res14.lon2 * math::r2d<CT>() << "," << res14.lat2 * math::r2d<CT>() << std::endl; - std::cout << "delta_g4=" << delta_g4 << std::endl; + std::cout << "a34=" << res34.azimuth * math::r2d<CT>() << std::endl; + std::cout << "a4=" << a4 * math::r2d<CT>() << std::endl; std::cout << "g4=" << g4 * math::r2d<CT>() << std::endl; + std::cout << "delta_g4=" << delta_g4 * math::r2d<CT>() << std::endl; std::cout << "der=" << der << std::endl; std::cout << "M43=" << M43 << std::endl; std::cout << "spherical limit=" << cos(s14/earth_radius) << std::endl; @@ -448,20 +483,20 @@ private : std::cout << "s34(sph) =" << s34_sph << std::endl; std::cout << "s34(geo) =" - << inverse_distance_quantities_type::apply(get<0>(p4), get<1>(p4), lon3, lat3, spheroid).distance + << inverse_distance_azimuth_quantities_type::apply(get<0>(p4), get<1>(p4), lon3, lat3, spheroid).distance << ", p4=(" << get<0>(p4) * math::r2d<double>() << "," << get<1>(p4) * math::r2d<double>() << ")" << std::endl; - CT s31 = inverse_distance_quantities_type::apply(lon3, lat3, lon1, lat1, spheroid).distance; - CT s32 = inverse_distance_quantities_type::apply(lon3, lat3, lon2, lat2, spheroid).distance; + CT s31 = inverse_distance_azimuth_quantities_type::apply(lon3, lat3, lon1, lat1, spheroid).distance; + CT s32 = inverse_distance_azimuth_quantities_type::apply(lon3, lat3, lon2, lat2, spheroid).distance; CT a4 = inverse_azimuth_type::apply(get<0>(p4), get<1>(p4), lon2, lat2, spheroid).azimuth; geometry::formula::result_direct<CT> res4 = direct_distance_type::apply(get<0>(p4), get<1>(p4), .04, a4, spheroid); - CT p4_plus = inverse_distance_quantities_type::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; + CT p4_plus = inverse_distance_azimuth_quantities_type::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; geometry::formula::result_direct<CT> res1 = direct_distance_type::apply(lon1, lat1, s14-.04, a12, spheroid); - CT p4_minus = inverse_distance_quantities_type::apply(res1.lon2, res1.lat2, lon3, lat3, spheroid).distance; + CT p4_minus = inverse_distance_azimuth_quantities_type::apply(res1.lon2, res1.lat2, lon3, lat3, spheroid).distance; std::cout << "s31=" << s31 << "\ns32=" << s32 << "\np4_plus=" << p4_plus << ", p4=(" << res4.lon2 * math::r2d<double>() << "," << res4.lat2 * math::r2d<double>() << ")" @@ -606,6 +641,30 @@ public : } }; +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + typename P, + typename PS +> +struct result_from_distance<geographic_cross_track<FormulaPolicy, Spheroid, CalculationType>, P, PS> +{ +private : + typedef typename geographic_cross_track + < + FormulaPolicy, Spheroid, CalculationType + >::template return_type<P, PS>::type return_type; +public : + template <typename T> + static inline return_type + apply(geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> const& , T const& distance) + { + return distance; + } +}; + template <typename Point, typename PointOfSegment> struct default_strategy diff --git a/boost/geometry/strategies/geographic/distance_cross_track_box_box.hpp b/boost/geometry/strategies/geographic/distance_cross_track_box_box.hpp new file mode 100644 index 0000000000..4f6b3b45b7 --- /dev/null +++ b/boost/geometry/strategies/geographic/distance_cross_track_box_box.hpp @@ -0,0 +1,207 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// 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) + +#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_CROSS_TRACK_BOX_BOX_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_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/strategies/spherical/distance_cross_track_box_box.hpp> + +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/algorithms/detail/assign_box_corners.hpp> + + +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 FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +class geographic_cross_track_box_box +{ +public: + typedef geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> Strategy; + + template <typename Box1, typename Box2> + struct return_type + : services::return_type<Strategy, typename point_type<Box1>::type, typename point_type<Box2>::type> + {}; + + inline geographic_cross_track_box_box() + {} + + 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, Strategy()); + } +}; + + + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +template <typename Strategy, typename Spheroid, typename CalculationType> +struct tag<geographic_cross_track_box_box<Strategy, Spheroid, CalculationType> > +{ + typedef strategy_tag_distance_box_box type; +}; + + +template <typename Strategy, typename Spheroid, typename CalculationType, typename Box1, typename Box2> +struct return_type<geographic_cross_track_box_box<Strategy, Spheroid, CalculationType>, Box1, Box2> + : geographic_cross_track_box_box + < + Strategy, Spheroid, CalculationType + >::template return_type<Box1, Box2> +{}; + +template <typename Strategy, typename Spheroid, typename Box1, typename Box2> +struct return_type<geographic_cross_track_box_box<Strategy, Spheroid>, Box1, Box2> + : geographic_cross_track_box_box + < + Strategy, Spheroid + >::template return_type<Box1, Box2> +{}; + +template <typename Strategy, typename Box1, typename Box2> +struct return_type<geographic_cross_track_box_box<Strategy>, Box1, Box2> + : geographic_cross_track_box_box + < + Strategy + >::template return_type<Box1, Box2> +{}; + +template <typename Strategy, typename Spheroid, typename CalculationType> +struct comparable_type<geographic_cross_track_box_box<Strategy, Spheroid, CalculationType> > +{ + typedef geographic_cross_track_box_box + < + typename comparable_type<Strategy>::type, Spheroid, CalculationType + > type; +}; + + +template <typename Strategy, typename Spheroid, typename CalculationType> +struct get_comparable<geographic_cross_track_box_box<Strategy, Spheroid, CalculationType> > +{ + typedef geographic_cross_track_box_box<Strategy, Spheroid, CalculationType> this_strategy; + typedef typename comparable_type<this_strategy>::type comparable_type; + +public: + static inline comparable_type apply(this_strategy const& /*strategy*/) + { + return comparable_type(); + } +}; + + +template <typename Strategy, typename Spheroid, typename CalculationType, typename Box1, typename Box2> +struct result_from_distance + < + geographic_cross_track_box_box<Strategy, Spheroid, CalculationType>, Box1, Box2 + > +{ +private: + typedef geographic_cross_track_box_box<Strategy, Spheroid, CalculationType> 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) + { + result_from_distance + < + Strategy, + typename point_type<Box1>::type, + typename point_type<Box2>::type + >::apply(strategy, distance); + } +}; + +template <typename Box1, typename Box2> +struct default_strategy + < + box_tag, box_tag, Box1, Box2, + geographic_tag, geographic_tag + > +{ + typedef geographic_cross_track_box_box<> type; +}; + + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::distance + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_CROSS_TRACK_BOX_BOX_HPP diff --git a/boost/geometry/strategies/geographic/distance_cross_track_point_box.hpp b/boost/geometry/strategies/geographic/distance_cross_track_point_box.hpp new file mode 100644 index 0000000000..016427428c --- /dev/null +++ b/boost/geometry/strategies/geographic/distance_cross_track_point_box.hpp @@ -0,0 +1,218 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// 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) + +#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_CROSS_TRACK_POINT_BOX_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_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/strategies/geographic/distance_cross_track.hpp> +#include <boost/geometry/strategies/spherical/distance_cross_track_point_box.hpp> + +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/algorithms/detail/assign_box_corners.hpp> + + +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 FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +class geographic_cross_track_point_box +{ +public: + typedef geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> Strategy; + + template <typename Point, typename Box> + struct return_type + : services::return_type<Strategy, Point, typename point_type<Box>::type> + {}; + + inline geographic_cross_track_point_box() + {} + + 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, Strategy()); + } +}; + + + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +template <typename Strategy, typename Spheroid, typename CalculationType> +struct tag<geographic_cross_track_point_box<Strategy, Spheroid, CalculationType> > +{ + typedef strategy_tag_distance_point_box type; +}; + + +template <typename Strategy, typename Spheroid, typename CalculationType, typename P, typename Box> +struct return_type<geographic_cross_track_point_box<Strategy, Spheroid, CalculationType>, P, Box> + : geographic_cross_track_point_box + < + Strategy, Spheroid, CalculationType + >::template return_type<P, Box> +{}; + +template <typename Strategy, typename Spheroid, typename P, typename Box> +struct return_type<geographic_cross_track_point_box<Strategy, Spheroid>, P, Box> + : geographic_cross_track_point_box + < + Strategy, Spheroid + >::template return_type<P, Box> +{}; + +template <typename Strategy, typename P, typename Box> +struct return_type<geographic_cross_track_point_box<Strategy>, P, Box> + : geographic_cross_track_point_box + < + Strategy + >::template return_type<P, Box> +{}; + +template <typename Strategy, typename Spheroid, typename CalculationType> +struct comparable_type<geographic_cross_track_point_box<Strategy, Spheroid, CalculationType> > +{ + typedef geographic_cross_track_point_box + < + Strategy, Spheroid, CalculationType + > type; +}; + + +template <typename Strategy, typename Spheroid, typename CalculationType> +struct get_comparable<geographic_cross_track_point_box<Strategy, Spheroid, CalculationType> > +{ + typedef geographic_cross_track_point_box<Strategy, Spheroid, CalculationType> this_strategy; + typedef typename comparable_type<this_strategy>::type comparable_type; + +public: + static inline comparable_type apply(this_strategy const&) + { + return comparable_type(); + } +}; + + +template <typename Strategy, typename Spheroid, typename CalculationType, typename P, typename Box> +struct result_from_distance + < + geographic_cross_track_point_box<Strategy, Spheroid, CalculationType>, P, Box + > +{ +private: + typedef geographic_cross_track_point_box<Strategy, Spheroid, CalculationType> this_strategy; + + typedef typename this_strategy::template return_type + < + P, Box + >::type return_type; + +public: + template <typename T> + static inline return_type apply(this_strategy const& strategy, + T const& distance) + { + result_from_distance + < + Strategy, P, typename point_type<Box>::type + >::apply(strategy, distance); + } +}; + +template <typename Point, typename Box> +struct default_strategy + < + point_tag, box_tag, Point, Box, + geographic_tag, geographic_tag + > +{ + typedef geographic_cross_track_point_box<> type; +}; + +template <typename Box, typename Point> +struct default_strategy + < + box_tag, point_tag, Box, Point, + geographic_tag, geographic_tag + > +{ + typedef typename default_strategy + < + point_tag, box_tag, Point, Box, + geographic_tag, geographic_tag + >::type type; +}; + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::distance + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_CROSS_TRACK_POINT_BOX_HPP diff --git a/boost/geometry/strategies/geographic/envelope_segment.hpp b/boost/geometry/strategies/geographic/envelope_segment.hpp index 3641b39428..5bada63f63 100644 --- a/boost/geometry/strategies/geographic/envelope_segment.hpp +++ b/boost/geometry/strategies/geographic/envelope_segment.hpp @@ -14,7 +14,9 @@ #include <boost/geometry/algorithms/detail/envelope/segment.hpp> #include <boost/geometry/algorithms/detail/normalize.hpp> -#include <boost/geometry/core/srs.hpp> + +#include <boost/geometry/srs/spheroid.hpp> + #include <boost/geometry/strategies/envelope.hpp> #include <boost/geometry/strategies/geographic/azimuth.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp> diff --git a/boost/geometry/strategies/geographic/intersection.hpp b/boost/geometry/strategies/geographic/intersection.hpp index e91659d40e..b017097f3e 100644 --- a/boost/geometry/strategies/geographic/intersection.hpp +++ b/boost/geometry/strategies/geographic/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 @@ -15,7 +17,6 @@ #include <boost/geometry/core/cs.hpp> #include <boost/geometry/core/access.hpp> #include <boost/geometry/core/radian_access.hpp> -#include <boost/geometry/core/srs.hpp> #include <boost/geometry/core/tags.hpp> #include <boost/geometry/algorithms/detail/assign_values.hpp> @@ -33,6 +34,8 @@ #include <boost/geometry/policies/robustness/segment_ratio.hpp> +#include <boost/geometry/srs/spheroid.hpp> + #include <boost/geometry/strategies/geographic/area.hpp> #include <boost/geometry/strategies/geographic/distance.hpp> #include <boost/geometry/strategies/geographic/envelope_segment.hpp> @@ -106,7 +109,6 @@ struct geographic_segments { typedef area::geographic < - typename point_type<Geometry>::type, FormulaPolicy, Order, Spheroid, @@ -152,34 +154,8 @@ struct geographic_segments template <typename CoordinateType, typename SegmentRatio> struct segment_intersection_info { - typedef typename select_most_precise - < - CoordinateType, double - >::type promoted_type; - - 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/geographic/intersection_elliptic.hpp b/boost/geometry/strategies/geographic/intersection_elliptic.hpp index 76e9940fe3..fce9735255 100644 --- a/boost/geometry/strategies/geographic/intersection_elliptic.hpp +++ b/boost/geometry/strategies/geographic/intersection_elliptic.hpp @@ -11,7 +11,7 @@ #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_ELLIPTIC_HPP -#include <boost/geometry/core/srs.hpp> +#include <boost/geometry/srs/spheroid.hpp> #include <boost/geometry/formulas/geographic.hpp> diff --git a/boost/geometry/strategies/geographic/side.hpp b/boost/geometry/strategies/geographic/side.hpp index 9276965a97..1201dc2f6d 100644 --- a/boost/geometry/strategies/geographic/side.hpp +++ b/boost/geometry/strategies/geographic/side.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014-2017. -// Modifications copyright (c) 2014-2017 Oracle and/or its affiliates. +// This file was modified by Oracle on 2014-2018. +// Modifications copyright (c) 2014-2018 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -18,10 +18,11 @@ #include <boost/geometry/core/access.hpp> #include <boost/geometry/core/radian_access.hpp> #include <boost/geometry/core/radius.hpp> -#include <boost/geometry/core/srs.hpp> #include <boost/geometry/formulas/spherical.hpp> +#include <boost/geometry/srs/spheroid.hpp> + #include <boost/geometry/util/math.hpp> #include <boost/geometry/util/promote_floating_point.hpp> #include <boost/geometry/util/select_calculation_type.hpp> @@ -48,6 +49,11 @@ namespace strategy { namespace side \tparam FormulaPolicy Geodesic solution formula policy. \tparam Spheroid Reference model of coordinate system. \tparam CalculationType \tparam_calculation + +\qbk{ +[heading See also] +[link geometry.reference.srs.srs_spheroid srs::spheroid] +} */ template < 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; diff --git a/boost/geometry/strategies/strategies.hpp b/boost/geometry/strategies/strategies.hpp index e7f3604abf..7d6cb618c6 100644 --- a/boost/geometry/strategies/strategies.hpp +++ b/boost/geometry/strategies/strategies.hpp @@ -3,9 +3,10 @@ // 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) 2017 Adam Wulkiewicz, Lodz, Poland. -// This file was modified by Oracle on 2014-2017. -// Modifications copyright (c) 2014-2017 Oracle and/or its affiliates. +// This file was modified by Oracle on 2014-2018. +// Modifications copyright (c) 2014-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 @@ -30,6 +31,7 @@ #include <boost/geometry/strategies/compare.hpp> #include <boost/geometry/strategies/convex_hull.hpp> #include <boost/geometry/strategies/covered_by.hpp> +#include <boost/geometry/strategies/densify.hpp> #include <boost/geometry/strategies/disjoint.hpp> #include <boost/geometry/strategies/distance.hpp> #include <boost/geometry/strategies/envelope.hpp> @@ -40,7 +42,7 @@ #include <boost/geometry/strategies/transform.hpp> #include <boost/geometry/strategies/within.hpp> -#include <boost/geometry/strategies/cartesian/area_surveyor.hpp> +#include <boost/geometry/strategies/cartesian/area.hpp> #include <boost/geometry/strategies/cartesian/azimuth.hpp> #include <boost/geometry/strategies/cartesian/box_in_box.hpp> #include <boost/geometry/strategies/cartesian/buffer_end_flat.hpp> @@ -54,6 +56,7 @@ #include <boost/geometry/strategies/cartesian/centroid_average.hpp> #include <boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp> #include <boost/geometry/strategies/cartesian/centroid_weighted_length.hpp> +#include <boost/geometry/strategies/cartesian/densify.hpp> #include <boost/geometry/strategies/cartesian/disjoint_segment_box.hpp> #include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp> #include <boost/geometry/strategies/cartesian/distance_pythagoras_point_box.hpp> @@ -70,9 +73,11 @@ #include <boost/geometry/strategies/spherical/area.hpp> #include <boost/geometry/strategies/spherical/azimuth.hpp> +#include <boost/geometry/strategies/spherical/densify.hpp> #include <boost/geometry/strategies/spherical/disjoint_segment_box.hpp> #include <boost/geometry/strategies/spherical/distance_haversine.hpp> #include <boost/geometry/strategies/spherical/distance_cross_track.hpp> +#include <boost/geometry/strategies/spherical/distance_cross_track_box_box.hpp> #include <boost/geometry/strategies/spherical/distance_cross_track_point_box.hpp> #include <boost/geometry/strategies/spherical/compare.hpp> #include <boost/geometry/strategies/spherical/envelope_segment.hpp> @@ -82,10 +87,13 @@ #include <boost/geometry/strategies/geographic/area.hpp> #include <boost/geometry/strategies/geographic/azimuth.hpp> +#include <boost/geometry/strategies/geographic/densify.hpp> #include <boost/geometry/strategies/geographic/disjoint_segment_box.hpp> #include <boost/geometry/strategies/geographic/distance.hpp> #include <boost/geometry/strategies/geographic/distance_andoyer.hpp> #include <boost/geometry/strategies/geographic/distance_cross_track.hpp> +#include <boost/geometry/strategies/geographic/distance_cross_track_box_box.hpp> +#include <boost/geometry/strategies/geographic/distance_cross_track_point_box.hpp> #include <boost/geometry/strategies/geographic/distance_thomas.hpp> #include <boost/geometry/strategies/geographic/distance_vincenty.hpp> #include <boost/geometry/strategies/geographic/envelope_segment.hpp> diff --git a/boost/geometry/strategies/transform/matrix_transformers.hpp b/boost/geometry/strategies/transform/matrix_transformers.hpp index e0ac6496f3..a3f3223c95 100644 --- a/boost/geometry/strategies/transform/matrix_transformers.hpp +++ b/boost/geometry/strategies/transform/matrix_transformers.hpp @@ -169,11 +169,11 @@ public : typedef typename geometry::coordinate_type<P2>::type ct2; set<0>(p2, boost::numeric_cast<ct2>( - c1 * m_matrix(0,0) + c2 * m_matrix(0,1) + c3 * m_matrix(0,2) + m_matrix(0,3))); + c1 * qvm::A<0,0>(m_matrix) + c2 * qvm::A<0,1>(m_matrix) + c3 * qvm::A<0,2>(m_matrix) + qvm::A<0,3>(m_matrix))); set<1>(p2, boost::numeric_cast<ct2>( - c1 * m_matrix(1,0) + c2 * m_matrix(1,1) + c3 * m_matrix(1,2) + m_matrix(1,3))); + c1 * qvm::A<1,0>(m_matrix) + c2 * qvm::A<1,1>(m_matrix) + c3 * qvm::A<1,2>(m_matrix) + qvm::A<1,3>(m_matrix))); set<2>(p2, boost::numeric_cast<ct2>( - c1 * m_matrix(2,0) + c2 * m_matrix(2,1) + c3 * m_matrix(2,2) + m_matrix(2,3))); + c1 * qvm::A<2,0>(m_matrix) + c2 * qvm::A<2,1>(m_matrix) + c3 * qvm::A<2,2>(m_matrix) + qvm::A<2,3>(m_matrix))); return true; } diff --git a/boost/geometry/strategies/transform/srs_transformer.hpp b/boost/geometry/strategies/transform/srs_transformer.hpp new file mode 100644 index 0000000000..e1815e2362 --- /dev/null +++ b/boost/geometry/strategies/transform/srs_transformer.hpp @@ -0,0 +1,102 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2008-2012 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 2017, 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_TRANSFORM_SRS_TRANSFORMER_HPP +#define BOOST_GEOMETRY_STRATEGIES_TRANSFORM_SRS_TRANSFORMER_HPP + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace transform +{ + +/*! + \brief Transformation strategy to do transform using a forward + Map Projection or SRS transformation. + \ingroup transform + \tparam ProjectionOrTransformation SRS projection or transformation type + */ +template +< + typename ProjectionOrTransformation +> +class srs_forward_transformer +{ +public: + inline srs_forward_transformer() + {} + + template <typename Parameters> + inline srs_forward_transformer(Parameters const& parameters) + : m_proj_or_transform(parameters) + {} + + template <typename Parameters1, typename Parameters2> + inline srs_forward_transformer(Parameters1 const& parameters1, Parameters2 const& parameters2) + : m_proj_or_transform(parameters1, parameters2) + {} + + template <typename Geometry1, typename Geometry2> + inline bool apply(Geometry1 const& g1, Geometry2 & g2) const + { + return m_proj_or_transform.forward(g1, g2); + } + +private: + ProjectionOrTransformation m_proj_or_transform; +}; + + +/*! + \brief Transformation strategy to do transform using an inverse + Map Projection or SRS transformation. + \ingroup transform + \tparam ProjectionOrTransformation SRS projection or transformation type + */ +template +< + typename ProjectionOrTransformation +> +class srs_inverse_transformer +{ +public: + inline srs_inverse_transformer() + {} + + template <typename Parameters> + inline srs_inverse_transformer(Parameters const& parameters) + : m_proj_or_transform(parameters) + {} + + template <typename Parameters1, typename Parameters2> + inline srs_inverse_transformer(Parameters1 const& parameters1, Parameters2 const& parameters2) + : m_proj_or_transform(parameters1, parameters2) + {} + + template <typename Geometry1, typename Geometry2> + inline bool apply(Geometry1 const& g1, Geometry2 & g2) const + { + return m_proj_or_transform.inverse(g1, g2); + } + +private: + ProjectionOrTransformation m_proj_or_transform; +}; + + +}} // namespace strategy::transform + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_TRANSFORM_SRS_TRANSFORMER_HPP |