diff options
Diffstat (limited to 'boost/geometry/strategies/geographic')
16 files changed, 2087 insertions, 345 deletions
diff --git a/boost/geometry/strategies/geographic/area.hpp b/boost/geometry/strategies/geographic/area.hpp new file mode 100644 index 0000000000..e1d3b09b5a --- /dev/null +++ b/boost/geometry/strategies/geographic/area.hpp @@ -0,0 +1,216 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2016-2017 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_GEOGRAPHIC_AREA_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP + + +#include <boost/geometry/core/srs.hpp> + +#include <boost/geometry/formulas/area_formulas.hpp> +#include <boost/geometry/formulas/flattening.hpp> + +#include <boost/geometry/strategies/geographic/parameters.hpp> + +#include <boost/math/special_functions/atanh.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace area +{ + +/*! +\brief Geographic area calculation +\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 +\tparam CalculationType \tparam_calculation +\author See +- Danielsen JS, The area under the geodesic. Surv Rev 30(232): 61–66, 1989 +- Charles F.F Karney, Algorithms for geodesics, 2011 https://arxiv.org/pdf/1109.4448.pdf + +\qbk{ +[heading See also] +[link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)] +} +*/ +template +< + typename PointOfSegment, + typename FormulaPolicy = strategy::andoyer, + std::size_t SeriesOrder = strategy::default_order<FormulaPolicy>::value, + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +class geographic +{ + // Switch between two kinds of approximation(series in eps and n v.s.series in k ^ 2 and e'^2) + static const bool ExpandEpsN = true; + // LongSegment Enables special handling of long segments + static const bool LongSegment = false; + + //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 + < + typename coordinate_type<PointOfSegment>::type, + double + >::type, + CalculationType + >::type CT; + +protected : + struct spheroid_constants + { + 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; // authalic radius + + inline spheroid_constants(Spheroid const& spheroid) + : m_spheroid(spheroid) + , m_a2(math::sqr(get_radius<0>(spheroid))) + , m_e2(formula::flattening<CT>(spheroid) + * (CT(2.0) - CT(formula::flattening<CT>(spheroid)))) + , m_ep2(m_e2 / (CT(1.0) - m_e2)) + , m_ep(math::sqrt(m_ep2)) + , m_c2((m_a2 / CT(2.0)) + + ((math::sqr(get_radius<2>(spheroid)) * boost::math::atanh(math::sqrt(m_e2))) + / (CT(2.0) * math::sqrt(m_e2)))) + {} + }; + + struct area_sums + { + CT m_excess_sum; + CT m_correction_sum; + + // Keep track if encircles some pole + std::size_t m_crosses_prime_meridian; + + inline area_sums() + : m_excess_sum(0) + , m_correction_sum(0) + , m_crosses_prime_meridian(0) + {} + inline CT area(spheroid_constants spheroid_const) const + { + CT result; + + CT sum = spheroid_const.m_c2 * m_excess_sum + + spheroid_const.m_e2 * spheroid_const.m_a2 * m_correction_sum; + + // If encircles some pole + if (m_crosses_prime_meridian % 2 == 1) + { + std::size_t times_crosses_prime_meridian + = 1 + (m_crosses_prime_meridian / 2); + + result = CT(2.0) + * geometry::math::pi<CT>() + * spheroid_const.m_c2 + * CT(times_crosses_prime_meridian) + - geometry::math::abs(sum); + + if (geometry::math::sign<CT>(sum) == 1) + { + result = - result; + } + + } + else + { + result = sum; + } + + return result; + } + }; + +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) + {} + + inline void apply(PointOfSegment const& p1, + PointOfSegment const& p2, + area_sums& state) const + { + + if (! geometry::math::equals(get<0>(p1), get<0>(p2))) + { + + typedef geometry::formula::area_formulas + < + CT, 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; + + // Keep track whenever a segment crosses the prime meridian + geometry::formula::area_formulas<CT> + ::crosses_prime_meridian(p1, p2, state); + } + } + + inline return_type result(area_sums const& state) const + { + return state.area(m_spheroid_constants); + } + +private: + spheroid_constants m_spheroid_constants; + +}; + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +namespace services +{ + + +template <typename Point> +struct default_strategy<geographic_tag, Point> +{ + typedef strategy::area::geographic<Point> type; +}; + +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +} + +}} // namespace strategy::area + + + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP diff --git a/boost/geometry/strategies/geographic/azimuth.hpp b/boost/geometry/strategies/geographic/azimuth.hpp new file mode 100644 index 0000000000..47f59d1033 --- /dev/null +++ b/boost/geometry/strategies/geographic/azimuth.hpp @@ -0,0 +1,103 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2016-2017 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_GEOGRAPHIC_AZIMUTH_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AZIMUTH_HPP + + +#include <boost/geometry/core/srs.hpp> + +#include <boost/geometry/strategies/azimuth.hpp> +#include <boost/geometry/strategies/geographic/parameters.hpp> + +#include <boost/mpl/if.hpp> +#include <boost/type_traits/is_void.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace azimuth +{ + +template +< + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +class geographic +{ +public : + + typedef Spheroid model_type; + + inline geographic() + : m_spheroid() + {} + + explicit inline geographic(Spheroid const& spheroid) + : m_spheroid(spheroid) + {} + + inline model_type const& model() const + { + return m_spheroid; + } + + template <typename T> + inline void apply(T const& lon1_rad, T const& lat1_rad, + T const& lon2_rad, T const& lat2_rad, + T& a1, T& a2) const + { + typedef typename boost::mpl::if_ + < + boost::is_void<CalculationType>, T, CalculationType + >::type calc_t; + + typedef typename FormulaPolicy::template inverse<calc_t, false, true, true, false, false> inverse_type; + typedef typename inverse_type::result_type inverse_result; + inverse_result i_res = inverse_type::apply(calc_t(lon1_rad), calc_t(lat1_rad), + calc_t(lon2_rad), calc_t(lat2_rad), + m_spheroid); + a1 = i_res.azimuth; + a2 = i_res.reverse_azimuth; + } + +private : + Spheroid m_spheroid; +}; + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +namespace services +{ + +template <typename CalculationType> +struct default_strategy<geographic_tag, CalculationType> +{ + typedef strategy::azimuth::geographic + < + strategy::andoyer, + srs::spheroid<double>, + CalculationType + > type; +}; + +} + +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +}} // namespace strategy::azimuth + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AZIMUTH_HPP diff --git a/boost/geometry/strategies/geographic/distance.hpp b/boost/geometry/strategies/geographic/distance.hpp new file mode 100644 index 0000000000..d3656f449c --- /dev/null +++ b/boost/geometry/strategies/geographic/distance.hpp @@ -0,0 +1,195 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// 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. + +// 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_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_HPP + + +#include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/core/radian_access.hpp> +#include <boost/geometry/core/radius.hpp> +#include <boost/geometry/core/srs.hpp> + +#include <boost/geometry/formulas/andoyer_inverse.hpp> +#include <boost/geometry/formulas/flattening.hpp> + +#include <boost/geometry/strategies/distance.hpp> +#include <boost/geometry/strategies/geographic/parameters.hpp> + +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/promote_floating_point.hpp> +#include <boost/geometry/util/select_calculation_type.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace distance +{ + +template +< + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +class geographic +{ +public : + template <typename Point1, typename Point2> + struct calculation_type + : promote_floating_point + < + typename select_calculation_type + < + Point1, + Point2, + CalculationType + >::type + > + {}; + + typedef Spheroid model_type; + + inline geographic() + : m_spheroid() + {} + + explicit inline geographic(Spheroid const& spheroid) + : m_spheroid(spheroid) + {} + + template <typename Point1, typename Point2> + inline typename calculation_type<Point1, Point2>::type + apply(Point1 const& point1, Point2 const& point2) const + { + return FormulaPolicy::template inverse + < + typename calculation_type<Point1, Point2>::type, + true, false, false, false, false + >::apply(get_as_radian<0>(point1), get_as_radian<1>(point1), + get_as_radian<0>(point2), get_as_radian<1>(point2), + m_spheroid).distance; + } + + inline Spheroid const& model() const + { + return m_spheroid; + } + +private : + Spheroid m_spheroid; +}; + + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct tag<geographic<FormulaPolicy, Spheroid, CalculationType> > +{ + typedef strategy_tag_distance_point_point type; +}; + + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + typename P1, + typename P2 +> +struct return_type<geographic<FormulaPolicy, Spheroid, CalculationType>, P1, P2> + : geographic<FormulaPolicy, Spheroid, CalculationType>::template calculation_type<P1, P2> +{}; + + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct comparable_type<geographic<FormulaPolicy, Spheroid, CalculationType> > +{ + typedef geographic<FormulaPolicy, Spheroid, CalculationType> type; +}; + + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct get_comparable<geographic<FormulaPolicy, Spheroid, CalculationType> > +{ + static inline geographic<FormulaPolicy, Spheroid, CalculationType> + apply(geographic<FormulaPolicy, Spheroid, CalculationType> const& input) + { + return input; + } +}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + typename P1, + typename P2 +> +struct result_from_distance<geographic<FormulaPolicy, Spheroid, CalculationType>, P1, P2> +{ + template <typename T> + static inline typename return_type<geographic<FormulaPolicy, Spheroid, CalculationType>, P1, P2>::type + apply(geographic<FormulaPolicy, Spheroid, CalculationType> const& , T const& value) + { + return value; + } +}; + + +template <typename Point1, typename Point2> +struct default_strategy<point_tag, point_tag, Point1, Point2, geographic_tag, geographic_tag> +{ + typedef strategy::distance::geographic + < + strategy::andoyer, + srs::spheroid + < + typename select_coordinate_type<Point1, Point2>::type + > + > type; +}; + + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::distance + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_HPP diff --git a/boost/geometry/strategies/geographic/distance_andoyer.hpp b/boost/geometry/strategies/geographic/distance_andoyer.hpp index 1946cd1090..d732951642 100644 --- a/boost/geometry/strategies/geographic/distance_andoyer.hpp +++ b/boost/geometry/strategies/geographic/distance_andoyer.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2016 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014, 2016. -// Modifications copyright (c) 2014-2016 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 @@ -11,24 +11,12 @@ // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) -#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ANDOYER_HPP -#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ANDOYER_HPP +#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_DETAIL_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_DETAIL_HPP -#include <boost/geometry/core/coordinate_type.hpp> -#include <boost/geometry/core/radian_access.hpp> -#include <boost/geometry/core/radius.hpp> -#include <boost/geometry/core/srs.hpp> - -#include <boost/geometry/algorithms/detail/flattening.hpp> - -#include <boost/geometry/formulas/andoyer_inverse.hpp> - -#include <boost/geometry/strategies/distance.hpp> - -#include <boost/geometry/util/math.hpp> -#include <boost/geometry/util/promote_floating_point.hpp> -#include <boost/geometry/util/select_calculation_type.hpp> +#include <boost/geometry/strategies/geographic/distance.hpp> +#include <boost/geometry/strategies/geographic/parameters.hpp> namespace boost { namespace geometry @@ -57,55 +45,28 @@ are about the same as Vincenty. In my (Barend's) testcases the results didn't di */ template < - typename Spheroid, + typename Spheroid = srs::spheroid<double>, typename CalculationType = void > class andoyer + : public strategy::distance::geographic + < + strategy::andoyer, Spheroid, CalculationType + > { -public : - template <typename Point1, typename Point2> - struct calculation_type - : promote_floating_point - < - typename select_calculation_type - < - Point1, - Point2, - CalculationType - >::type - > - {}; - - typedef Spheroid model_type; + typedef strategy::distance::geographic + < + strategy::andoyer, Spheroid, CalculationType + > base_type; +public : inline andoyer() - : m_spheroid() + : base_type() {} explicit inline andoyer(Spheroid const& spheroid) - : m_spheroid(spheroid) + : base_type(spheroid) {} - - template <typename Point1, typename Point2> - inline typename calculation_type<Point1, Point2>::type - apply(Point1 const& point1, Point2 const& point2) const - { - return geometry::formula::andoyer_inverse - < - typename calculation_type<Point1, Point2>::type, - true, false - >::apply(get_as_radian<0>(point1), get_as_radian<1>(point1), - get_as_radian<0>(point2), get_as_radian<1>(point2), - m_spheroid).distance; - } - - inline Spheroid const& model() const - { - return m_spheroid; - } - -private : - Spheroid m_spheroid; }; @@ -154,19 +115,6 @@ struct result_from_distance<andoyer<Spheroid, CalculationType>, P1, P2> }; -template <typename Point1, typename Point2> -struct default_strategy<point_tag, point_tag, Point1, Point2, geographic_tag, geographic_tag> -{ - typedef strategy::distance::andoyer - < - srs::spheroid - < - typename select_coordinate_type<Point1, Point2>::type - > - > type; -}; - - } // namespace services #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS @@ -177,4 +125,4 @@ struct default_strategy<point_tag, point_tag, Point1, Point2, geographic_tag, ge }} // namespace boost::geometry -#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ANDOYER_HPP +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_DETAIL_HPP diff --git a/boost/geometry/strategies/geographic/distance_thomas.hpp b/boost/geometry/strategies/geographic/distance_thomas.hpp index 39e0ecfa6f..490920c778 100644 --- a/boost/geometry/strategies/geographic/distance_thomas.hpp +++ b/boost/geometry/strategies/geographic/distance_thomas.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2015, 2016. -// Modifications copyright (c) 2015-2016 Oracle and/or its affiliates. +// This file was modified by Oracle on 2015-2017. +// Modifications copyright (c) 2015-2017 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -15,15 +15,9 @@ #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_THOMAS_HPP -#include <boost/geometry/core/coordinate_type.hpp> -#include <boost/geometry/core/radian_access.hpp> +#include <boost/geometry/strategies/geographic/distance.hpp> +#include <boost/geometry/strategies/geographic/parameters.hpp> -#include <boost/geometry/strategies/distance.hpp> - -#include <boost/geometry/util/promote_floating_point.hpp> -#include <boost/geometry/util/select_calculation_type.hpp> - -#include <boost/geometry/formulas/thomas_inverse.hpp> namespace boost { namespace geometry { @@ -45,57 +39,28 @@ namespace strategy { namespace distance */ template < - typename Spheroid, + typename Spheroid = srs::spheroid<double>, typename CalculationType = void > class thomas + : public strategy::distance::geographic + < + strategy::thomas, Spheroid, CalculationType + > { -public : - template <typename Point1, typename Point2> - struct calculation_type - : promote_floating_point - < - typename select_calculation_type - < - Point1, - Point2, - CalculationType - >::type - > - {}; - - typedef Spheroid model_type; + typedef strategy::distance::geographic + < + strategy::thomas, Spheroid, CalculationType + > base_type; +public : inline thomas() - : m_spheroid() + : base_type() {} explicit inline thomas(Spheroid const& spheroid) - : m_spheroid(spheroid) + : base_type(spheroid) {} - - template <typename Point1, typename Point2> - inline typename calculation_type<Point1, Point2>::type - apply(Point1 const& point1, Point2 const& point2) const - { - return geometry::formula::thomas_inverse - < - typename calculation_type<Point1, Point2>::type, - true, false - >::apply(get_as_radian<0>(point1), - get_as_radian<1>(point1), - get_as_radian<0>(point2), - get_as_radian<1>(point2), - m_spheroid).distance; - } - - inline Spheroid const& model() const - { - return m_spheroid; - } - -private : - Spheroid m_spheroid; }; #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS diff --git a/boost/geometry/strategies/geographic/distance_vincenty.hpp b/boost/geometry/strategies/geographic/distance_vincenty.hpp index e79e9aeb46..41146db9ff 100644 --- a/boost/geometry/strategies/geographic/distance_vincenty.hpp +++ b/boost/geometry/strategies/geographic/distance_vincenty.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014, 2016. -// Modifications copyright (c) 2014-2016 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 @@ -15,15 +15,9 @@ #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_VINCENTY_HPP -#include <boost/geometry/core/coordinate_type.hpp> -#include <boost/geometry/core/radian_access.hpp> +#include <boost/geometry/strategies/geographic/distance.hpp> +#include <boost/geometry/strategies/geographic/parameters.hpp> -#include <boost/geometry/strategies/distance.hpp> - -#include <boost/geometry/util/promote_floating_point.hpp> -#include <boost/geometry/util/select_calculation_type.hpp> - -#include <boost/geometry/formulas/vincenty_inverse.hpp> namespace boost { namespace geometry { @@ -47,57 +41,28 @@ namespace strategy { namespace distance */ template < - typename Spheroid, + typename Spheroid = srs::spheroid<double>, typename CalculationType = void > class vincenty + : public strategy::distance::geographic + < + strategy::vincenty, Spheroid, CalculationType + > { -public : - template <typename Point1, typename Point2> - struct calculation_type - : promote_floating_point - < - typename select_calculation_type - < - Point1, - Point2, - CalculationType - >::type - > - {}; - - typedef Spheroid model_type; + typedef strategy::distance::geographic + < + strategy::vincenty, Spheroid, CalculationType + > base_type; +public: inline vincenty() - : m_spheroid() + : base_type() {} explicit inline vincenty(Spheroid const& spheroid) - : m_spheroid(spheroid) + : base_type(spheroid) {} - - template <typename Point1, typename Point2> - inline typename calculation_type<Point1, Point2>::type - apply(Point1 const& point1, Point2 const& point2) const - { - return geometry::formula::vincenty_inverse - < - typename calculation_type<Point1, Point2>::type, - true, false - >::apply(get_as_radian<0>(point1), - get_as_radian<1>(point1), - get_as_radian<0>(point2), - get_as_radian<1>(point2), - m_spheroid).distance; - } - - inline Spheroid const& model() const - { - return m_spheroid; - } - -private : - Spheroid m_spheroid; }; #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS diff --git a/boost/geometry/strategies/geographic/envelope_segment.hpp b/boost/geometry/strategies/geographic/envelope_segment.hpp new file mode 100644 index 0000000000..3641b39428 --- /dev/null +++ b/boost/geometry/strategies/geographic/envelope_segment.hpp @@ -0,0 +1,104 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2017 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_GEOGRAPHIC_ENVELOPE_SEGMENT_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ENVELOPE_SEGMENT_HPP + + +#include <boost/geometry/algorithms/detail/envelope/segment.hpp> +#include <boost/geometry/algorithms/detail/normalize.hpp> +#include <boost/geometry/core/srs.hpp> +#include <boost/geometry/strategies/envelope.hpp> +#include <boost/geometry/strategies/geographic/azimuth.hpp> +#include <boost/geometry/strategies/geographic/parameters.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace envelope +{ + +template +< + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = geometry::srs::spheroid<double>, + typename CalculationType = void +> +class geographic_segment +{ +public: + typedef Spheroid model_type; + + inline geographic_segment() + : m_spheroid() + {} + + explicit inline geographic_segment(Spheroid const& spheroid) + : m_spheroid(spheroid) + {} + + template <typename Point1, typename Point2, typename Box> + inline void apply(Point1 const& point1, Point2 const& point2, Box& box) const + { + Point1 p1_normalized = detail::return_normalized<Point1>(point1); + Point2 p2_normalized = detail::return_normalized<Point2>(point2); + + geometry::strategy::azimuth::geographic + < + FormulaPolicy, + Spheroid, + CalculationType + > azimuth_geographic(m_spheroid); + + typedef typename coordinate_system<Point1>::type::units units_type; + + detail::envelope::envelope_segment_impl + < + geographic_tag + >::template apply<units_type>(geometry::get<0>(p1_normalized), + geometry::get<1>(p1_normalized), + geometry::get<0>(p2_normalized), + geometry::get<1>(p2_normalized), + box, + azimuth_geographic); + + } + +private: + Spheroid m_spheroid; +}; + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +namespace services +{ + +template <typename CalculationType> +struct default_strategy<geographic_tag, CalculationType> +{ + typedef strategy::envelope::geographic_segment + < + strategy::andoyer, + srs::spheroid<double>, + CalculationType + > type; +}; + +} + +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::envelope + +}} //namepsace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ENVELOPE_SEGMENT_HPP diff --git a/boost/geometry/strategies/geographic/intersection.hpp b/boost/geometry/strategies/geographic/intersection.hpp new file mode 100644 index 0000000000..1708c274c0 --- /dev/null +++ b/boost/geometry/strategies/geographic/intersection.hpp @@ -0,0 +1,897 @@ +// Boost.Geometry + +// Copyright (c) 2016-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_GEOGRAPHIC_INTERSECTION_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP + +#include <algorithm> + +#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> +#include <boost/geometry/algorithms/detail/assign_indexed_point.hpp> +#include <boost/geometry/algorithms/detail/equals/point_point.hpp> +#include <boost/geometry/algorithms/detail/recalculate.hpp> + +#include <boost/geometry/formulas/andoyer_inverse.hpp> +#include <boost/geometry/formulas/sjoberg_intersection.hpp> +#include <boost/geometry/formulas/spherical.hpp> + +#include <boost/geometry/geometries/concepts/point_concept.hpp> +#include <boost/geometry/geometries/concepts/segment_concept.hpp> + +#include <boost/geometry/policies/robustness/segment_ratio.hpp> + +#include <boost/geometry/strategies/geographic/area.hpp> +#include <boost/geometry/strategies/geographic/distance.hpp> +#include <boost/geometry/strategies/geographic/parameters.hpp> +#include <boost/geometry/strategies/geographic/side.hpp> +#include <boost/geometry/strategies/intersection.hpp> +#include <boost/geometry/strategies/intersection_result.hpp> +#include <boost/geometry/strategies/side_info.hpp> + +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/select_calculation_type.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace intersection +{ + +// CONSIDER: Improvement of the robustness/accuracy/repeatability by +// moving all segments to 0 longitude +// picking latitudes closer to 0 +// etc. + +template +< + typename FormulaPolicy = strategy::andoyer, + unsigned int Order = strategy::default_order<FormulaPolicy>::value, + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +struct geographic_segments +{ + typedef side::geographic + < + FormulaPolicy, Spheroid, CalculationType + > side_strategy_type; + + inline side_strategy_type get_side_strategy() const + { + return side_strategy_type(m_spheroid); + } + + template <typename Geometry1, typename Geometry2> + struct point_in_geometry_strategy + { + typedef strategy::within::winding + < + typename point_type<Geometry1>::type, + typename point_type<Geometry2>::type, + side_strategy_type, + CalculationType + > type; + }; + + template <typename Geometry1, typename Geometry2> + inline typename point_in_geometry_strategy<Geometry1, Geometry2>::type + get_point_in_geometry_strategy() const + { + typedef typename point_in_geometry_strategy + < + Geometry1, Geometry2 + >::type strategy_type; + return strategy_type(get_side_strategy()); + } + + template <typename Geometry> + struct area_strategy + { + typedef area::geographic + < + typename point_type<Geometry>::type, + FormulaPolicy, + Order, + Spheroid, + CalculationType + > type; + }; + + template <typename Geometry> + inline typename area_strategy<Geometry>::type get_area_strategy() const + { + typedef typename area_strategy<Geometry>::type strategy_type; + return strategy_type(m_spheroid); + } + + template <typename Geometry> + struct distance_strategy + { + typedef distance::geographic + < + FormulaPolicy, + Spheroid, + CalculationType + > type; + }; + + template <typename Geometry> + inline typename distance_strategy<Geometry>::type get_distance_strategy() const + { + typedef typename distance_strategy<Geometry>::type strategy_type; + return strategy_type(m_spheroid); + } + + enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 }; + + 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 + { + if (ip_flag == ipi_inters) + { + // TODO: assign the rest of coordinates + set_from_radian<0>(point, lon); + set_from_radian<1>(point, lat); + } + else if (ip_flag == ipi_at_a1) + { + detail::assign_point_from_index<0>(a, point); + } + else if (ip_flag == ipi_at_a2) + { + detail::assign_point_from_index<1>(a, point); + } + else if (ip_flag == ipi_at_b1) + { + detail::assign_point_from_index<0>(b, point); + } + else // ip_flag == ipi_at_b2 + { + detail::assign_point_from_index<1>(b, point); + } + } + + CoordinateType lon; + CoordinateType lat; + SegmentRatio robust_ra; + SegmentRatio robust_rb; + intersection_point_flag ip_flag; + }; + + explicit geographic_segments(Spheroid const& spheroid = Spheroid()) + : m_spheroid(spheroid) + {} + + // Relate segments a and b + template + < + typename Segment1, + typename Segment2, + typename Policy, + typename RobustPolicy + > + inline typename Policy::return_type apply(Segment1 const& a, Segment2 const& b, + Policy const& policy, + RobustPolicy const& robust_policy) const + { + typedef typename point_type<Segment1>::type point1_t; + typedef typename point_type<Segment2>::type point2_t; + point1_t a1, a2; + point2_t b1, b2; + + detail::assign_point_from_index<0>(a, a1); + detail::assign_point_from_index<1>(a, a2); + detail::assign_point_from_index<0>(b, b1); + detail::assign_point_from_index<1>(b, b2); + + return apply(a, b, policy, robust_policy, a1, a2, b1, b2); + } + + // Relate segments a and b + template + < + typename Segment1, + typename Segment2, + typename Policy, + typename RobustPolicy, + typename Point1, + typename Point2 + > + inline typename Policy::return_type apply(Segment1 const& a, Segment2 const& b, + Policy const&, RobustPolicy const&, + Point1 a1, Point1 a2, Point2 b1, Point2 b2) const + { + bool is_a_reversed = get<1>(a1) > get<1>(a2); + bool is_b_reversed = get<1>(b1) > get<1>(b2); + + if (is_a_reversed) + { + std::swap(a1, a2); + } + + if (is_b_reversed) + { + std::swap(b1, b2); + } + + return apply<Policy>(a, b, a1, a2, b1, b2, is_a_reversed, is_b_reversed); + } + +private: + // Relate segments a and b + template + < + typename Policy, + typename Segment1, + typename Segment2, + typename Point1, + typename Point2 + > + inline typename Policy::return_type apply(Segment1 const& a, Segment2 const& b, + Point1 const& a1, Point1 const& a2, + Point2 const& b1, Point2 const& b2, + bool is_a_reversed, bool is_b_reversed) const + { + BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) ); + BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) ); + + typedef typename select_calculation_type + <Segment1, Segment2, CalculationType>::type calc_t; + + // normalized spheroid + srs::spheroid<calc_t> spheroid = normalized_spheroid<calc_t>(m_spheroid); + + // TODO: check only 2 first coordinates here? + using geometry::detail::equals::equals_point_point; + bool a_is_point = equals_point_point(a1, a2); + bool b_is_point = equals_point_point(b1, b2); + + if(a_is_point && b_is_point) + { + return equals_point_point(a1, b2) + ? Policy::degenerate(a, true) + : Policy::disjoint() + ; + } + + calc_t const a1_lon = get_as_radian<0>(a1); + calc_t const a1_lat = get_as_radian<1>(a1); + calc_t const a2_lon = get_as_radian<0>(a2); + calc_t const a2_lat = get_as_radian<1>(a2); + calc_t const b1_lon = get_as_radian<0>(b1); + calc_t const b1_lat = get_as_radian<1>(b1); + calc_t const b2_lon = get_as_radian<0>(b2); + calc_t const b2_lat = get_as_radian<1>(b2); + + side_info sides; + + // NOTE: potential optimization, don't calculate distance at this point + // this would require to reimplement inverse strategy to allow + // calculation of distance if needed, probably also storing intermediate + // results somehow inside an object. + typedef typename FormulaPolicy::template inverse<calc_t, true, true, false, false, false> inverse_dist_azi; + typedef typename inverse_dist_azi::result_type inverse_result; + + // TODO: no need to call inverse formula if we know that the points are equal + // distance can be set to 0 in this case and azimuth may be not calculated + bool const is_equal_a1_b1 = equals_point_point(a1, b1); + bool const is_equal_a2_b1 = equals_point_point(a2, b1); + + inverse_result res_b1_b2 = inverse_dist_azi::apply(b1_lon, b1_lat, b2_lon, b2_lat, spheroid); + inverse_result res_b1_a1 = inverse_dist_azi::apply(b1_lon, b1_lat, a1_lon, a1_lat, spheroid); + inverse_result res_b1_a2 = inverse_dist_azi::apply(b1_lon, b1_lat, a2_lon, a2_lat, spheroid); + sides.set<0>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_b1_a1.azimuth, res_b1_b2.azimuth), + is_equal_a2_b1 ? 0 : formula::azimuth_side_value(res_b1_a2.azimuth, res_b1_b2.azimuth)); + if (sides.same<0>()) + { + // Both points are at the same side of other segment, we can leave + return Policy::disjoint(); + } + + bool const is_equal_a1_b2 = equals_point_point(a1, b2); + + inverse_result res_a1_a2 = inverse_dist_azi::apply(a1_lon, a1_lat, a2_lon, a2_lat, spheroid); + inverse_result res_a1_b1 = inverse_dist_azi::apply(a1_lon, a1_lat, b1_lon, b1_lat, spheroid); + inverse_result res_a1_b2 = inverse_dist_azi::apply(a1_lon, a1_lat, b2_lon, b2_lat, spheroid); + sides.set<1>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_a1_b1.azimuth, res_a1_a2.azimuth), + is_equal_a1_b2 ? 0 : formula::azimuth_side_value(res_a1_b2.azimuth, res_a1_a2.azimuth)); + if (sides.same<1>()) + { + // Both points are at the same side of other segment, we can leave + return Policy::disjoint(); + } + + // NOTE: at this point the segments may still be disjoint + // NOTE: at this point one of the segments may be degenerated + + bool collinear = sides.collinear(); + + if (! collinear) + { + // WARNING: the side strategy doesn't have the info about the other + // segment so it may return results inconsistent with this intersection + // strategy, as it checks both segments for consistency + + if (sides.get<0, 0>() == 0 && sides.get<0, 1>() == 0) + { + collinear = true; + sides.set<1>(0, 0); + } + else if (sides.get<1, 0>() == 0 && sides.get<1, 1>() == 0) + { + collinear = true; + sides.set<0>(0, 0); + } + } + + if (collinear) + { + if (a_is_point) + { + return collinear_one_degenerated<Policy, calc_t>(a, true, b1, b2, a1, a2, res_b1_b2, res_b1_a1, is_b_reversed); + } + else if (b_is_point) + { + return collinear_one_degenerated<Policy, calc_t>(b, false, a1, a2, b1, b2, res_a1_a2, res_a1_b1, is_a_reversed); + } + else + { + calc_t dist_a1_a2, dist_a1_b1, dist_a1_b2; + calc_t dist_b1_b2, dist_b1_a1, dist_b1_a2; + // use shorter segment + if (res_a1_a2.distance <= res_b1_b2.distance) + { + calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, dist_a1_a2, dist_a1_b1); + calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b2, dist_a1_a2, dist_a1_b2); + dist_b1_b2 = dist_a1_b2 - dist_a1_b1; + dist_b1_a1 = -dist_a1_b1; + dist_b1_a2 = dist_a1_a2 - dist_a1_b1; + } + else + { + calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a1, dist_b1_b2, dist_b1_a1); + calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a2, dist_b1_b2, dist_b1_a2); + dist_a1_a2 = dist_b1_a2 - dist_b1_a1; + dist_a1_b1 = -dist_b1_a1; + dist_a1_b2 = dist_b1_b2 - dist_b1_a1; + } + + // NOTE: this is probably not needed + calc_t const c0 = 0; + int a1_on_b = position_value(c0, dist_a1_b1, dist_a1_b2); + int a2_on_b = position_value(dist_a1_a2, dist_a1_b1, dist_a1_b2); + int b1_on_a = position_value(c0, dist_b1_a1, dist_b1_a2); + int b2_on_a = position_value(dist_b1_b2, dist_b1_a1, dist_b1_a2); + + if ((a1_on_b < 1 && a2_on_b < 1) || (a1_on_b > 3 && a2_on_b > 3)) + { + return Policy::disjoint(); + } + + if (a1_on_b == 1) + { + dist_b1_a1 = 0; + dist_a1_b1 = 0; + } + else if (a1_on_b == 3) + { + dist_b1_a1 = dist_b1_b2; + dist_a1_b2 = 0; + } + + if (a2_on_b == 1) + { + dist_b1_a2 = 0; + dist_a1_b1 = dist_a1_a2; + } + else if (a2_on_b == 3) + { + dist_b1_a2 = dist_b1_b2; + dist_a1_b2 = dist_a1_a2; + } + + bool opposite = ! same_direction(res_a1_a2.azimuth, res_b1_b2.azimuth); + + // NOTE: If segment was reversed opposite, positions and segment ratios has to be altered + if (is_a_reversed) + { + // opposite + opposite = ! opposite; + // positions + std::swap(a1_on_b, a2_on_b); + b1_on_a = 4 - b1_on_a; + b2_on_a = 4 - b2_on_a; + // distances for ratios + std::swap(dist_b1_a1, dist_b1_a2); + dist_a1_b1 = dist_a1_a2 - dist_a1_b1; + dist_a1_b2 = dist_a1_a2 - dist_a1_b2; + } + if (is_b_reversed) + { + // opposite + opposite = ! opposite; + // positions + a1_on_b = 4 - a1_on_b; + a2_on_b = 4 - a2_on_b; + std::swap(b1_on_a, b2_on_a); + // distances for ratios + dist_b1_a1 = dist_b1_b2 - dist_b1_a1; + dist_b1_a2 = dist_b1_b2 - dist_b1_a2; + std::swap(dist_a1_b1, dist_a1_b2); + } + + segment_ratio<calc_t> ra_from(dist_b1_a1, dist_b1_b2); + segment_ratio<calc_t> ra_to(dist_b1_a2, dist_b1_b2); + segment_ratio<calc_t> rb_from(dist_a1_b1, dist_a1_a2); + segment_ratio<calc_t> rb_to(dist_a1_b2, dist_a1_a2); + + return Policy::segments_collinear(a, b, opposite, + a1_on_b, a2_on_b, b1_on_a, b2_on_a, + ra_from, ra_to, rb_from, rb_to); + } + } + else // crossing or touching + { + if (a_is_point || b_is_point) + { + return Policy::disjoint(); + } + + calc_t lon = 0, lat = 0; + intersection_point_flag ip_flag; + calc_t dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1; + if (calculate_ip_data(a1, a2, b1, b2, + a1_lon, a1_lat, a2_lon, a2_lat, + b1_lon, b1_lat, b2_lon, b2_lat, + res_a1_a2, res_a1_b1, res_a1_b2, + res_b1_b2, res_b1_a1, res_b1_a2, + sides, spheroid, + lon, lat, + dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1, + ip_flag)) + { + // NOTE: If segment was reversed sides and segment ratios has to be altered + if (is_a_reversed) + { + // sides + sides_reverse_segment<0>(sides); + // distance for ratio + dist_a1_i1 = dist_a1_a2 - dist_a1_i1; + // ip flag + ip_flag_reverse_segment(ip_flag, ipi_at_a1, ipi_at_a2); + } + if (is_b_reversed) + { + // sides + sides_reverse_segment<1>(sides); + // distance for ratio + dist_b1_i1 = dist_b1_b2 - dist_b1_i1; + // ip flag + ip_flag_reverse_segment(ip_flag, ipi_at_b1, ipi_at_b2); + } + + // intersects + segment_intersection_info + < + calc_t, + segment_ratio<calc_t> + > sinfo; + + sinfo.lon = lon; + sinfo.lat = lat; + sinfo.robust_ra.assign(dist_a1_i1, dist_a1_a2); + sinfo.robust_rb.assign(dist_b1_i1, dist_b1_b2); + sinfo.ip_flag = ip_flag; + + return Policy::segments_crosses(sides, sinfo, a, b); + } + else + { + return Policy::disjoint(); + } + } + } + + template <typename Policy, typename CalcT, typename Segment, typename Point1, typename Point2, typename ResultInverse> + static inline typename Policy::return_type + collinear_one_degenerated(Segment const& segment, bool degenerated_a, + Point1 const& a1, Point1 const& a2, + Point2 const& b1, Point2 const& b2, + ResultInverse const& res_a1_a2, + ResultInverse const& res_a1_bi, + bool is_other_reversed) + { + CalcT dist_1_2, dist_1_o; + if (! calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_bi, dist_1_2, dist_1_o)) + { + return Policy::disjoint(); + } + + // NOTE: If segment was reversed segment ratio has to be altered + if (is_other_reversed) + { + // distance for ratio + dist_1_o = dist_1_2 - dist_1_o; + } + + return Policy::one_degenerate(segment, segment_ratio<CalcT>(dist_1_o, dist_1_2), degenerated_a); + } + + // TODO: instead of checks below test bi against a1 and a2 here? + // in order to make this independent from is_near() + template <typename Point1, typename Point2, typename ResultInverse, typename CalcT> + static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in + Point2 const& b1, Point2 const& b2, // in + ResultInverse const& res_a1_a2, // in + ResultInverse const& res_a1_bi, // in + CalcT& dist_a1_a2, CalcT& dist_a1_bi) // out + { + dist_a1_a2 = res_a1_a2.distance; + + dist_a1_bi = res_a1_bi.distance; + if (! same_direction(res_a1_bi.azimuth, res_a1_a2.azimuth)) + { + dist_a1_bi = -dist_a1_bi; + } + + // if i1 is close to a1 and b1 or b2 is equal to a1 + if (is_endpoint_equal(dist_a1_bi, a1, b1, b2)) + { + dist_a1_bi = 0; + return true; + } + // or i1 is close to a2 and b1 or b2 is equal to a2 + else if (is_endpoint_equal(dist_a1_a2 - dist_a1_bi, a2, b1, b2)) + { + dist_a1_bi = dist_a1_a2; + return true; + } + + // or i1 is on b + return segment_ratio<CalcT>(dist_a1_bi, dist_a1_a2).on_segment(); + } + + template <typename Point1, typename Point2, typename CalcT, typename ResultInverse, typename Spheroid_> + static inline bool calculate_ip_data(Point1 const& a1, Point1 const& a2, // in + Point2 const& b1, Point2 const& b2, // in + CalcT const& a1_lon, CalcT const& a1_lat, // in + CalcT const& a2_lon, CalcT const& a2_lat, // in + CalcT const& b1_lon, CalcT const& b1_lat, // in + CalcT const& b2_lon, CalcT const& b2_lat, // in + ResultInverse const& res_a1_a2, // in + ResultInverse const& res_a1_b1, // in + ResultInverse const& res_a1_b2, // in + ResultInverse const& res_b1_b2, // in + ResultInverse const& res_b1_a1, // in + ResultInverse const& res_b1_a2, // in + side_info const& sides, // in + Spheroid_ const& spheroid, // in + CalcT & lon, CalcT & lat, // out + CalcT& dist_a1_a2, CalcT& dist_a1_ip, // out + CalcT& dist_b1_b2, CalcT& dist_b1_ip, // out + intersection_point_flag& ip_flag) // out + { + dist_a1_a2 = res_a1_a2.distance; + dist_b1_b2 = res_b1_b2.distance; + + // assign the IP if some endpoints overlap + using geometry::detail::equals::equals_point_point; + if (equals_point_point(a1, b1)) + { + lon = a1_lon; + lat = a1_lat; + dist_a1_ip = 0; + dist_b1_ip = 0; + ip_flag = ipi_at_a1; + return true; + } + else if (equals_point_point(a1, b2)) + { + lon = a1_lon; + lat = a1_lat; + dist_a1_ip = 0; + dist_b1_ip = dist_b1_b2; + ip_flag = ipi_at_a1; + return true; + } + else if (equals_point_point(a2, b1)) + { + lon = a2_lon; + lat = a2_lat; + dist_a1_ip = dist_a1_a2; + dist_b1_ip = 0; + ip_flag = ipi_at_a2; + return true; + } + else if (equals_point_point(a2, b2)) + { + lon = a2_lon; + lat = a2_lat; + dist_a1_ip = dist_a1_a2; + dist_b1_ip = dist_b1_b2; + ip_flag = ipi_at_a2; + return true; + } + + // at this point we know that the endpoints doesn't overlap + // check cases when an endpoint lies on the other geodesic + if (sides.template get<0, 0>() == 0) // a1 wrt b + { + if (res_b1_a1.distance <= res_b1_b2.distance + && same_direction(res_b1_a1.azimuth, res_b1_b2.azimuth)) + { + lon = a1_lon; + lat = a1_lat; + dist_a1_ip = 0; + dist_b1_ip = res_b1_a1.distance; + ip_flag = ipi_at_a1; + return true; + } + else + { + return false; + } + } + else if (sides.template get<0, 1>() == 0) // a2 wrt b + { + if (res_b1_a2.distance <= res_b1_b2.distance + && same_direction(res_b1_a2.azimuth, res_b1_b2.azimuth)) + { + lon = a2_lon; + lat = a2_lat; + dist_a1_ip = res_a1_a2.distance; + dist_b1_ip = res_b1_a2.distance; + ip_flag = ipi_at_a2; + return true; + } + else + { + return false; + } + } + else if (sides.template get<1, 0>() == 0) // b1 wrt a + { + if (res_a1_b1.distance <= res_a1_a2.distance + && same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth)) + { + lon = b1_lon; + lat = b1_lat; + dist_a1_ip = res_a1_b1.distance; + dist_b1_ip = 0; + ip_flag = ipi_at_b1; + return true; + } + else + { + return false; + } + } + else if (sides.template get<1, 1>() == 0) // b2 wrt a + { + if (res_a1_b2.distance <= res_a1_a2.distance + && same_direction(res_a1_b2.azimuth, res_a1_a2.azimuth)) + { + lon = b2_lon; + lat = b2_lat; + dist_a1_ip = res_a1_b2.distance; + dist_b1_ip = res_b1_b2.distance; + ip_flag = ipi_at_b2; + return true; + } + else + { + return false; + } + } + + // At this point neither the endpoints overlaps + // nor any andpoint lies on the other geodesic + // So the endpoints should lie on the opposite sides of both geodesics + + bool const ok = formula::sjoberg_intersection<CalcT, FormulaPolicy::template inverse, Order> + ::apply(a1_lon, a1_lat, a2_lon, a2_lat, res_a1_a2.azimuth, + b1_lon, b1_lat, b2_lon, b2_lat, res_b1_b2.azimuth, + lon, lat, spheroid); + + if (! ok) + { + return false; + } + + typedef typename FormulaPolicy::template inverse<CalcT, true, true, false, false, false> inverse_dist_azi; + typedef typename inverse_dist_azi::result_type inverse_result; + + inverse_result const res_a1_ip = inverse_dist_azi::apply(a1_lon, a1_lat, lon, lat, spheroid); + dist_a1_ip = res_a1_ip.distance; + if (! same_direction(res_a1_ip.azimuth, res_a1_a2.azimuth)) + { + dist_a1_ip = -dist_a1_ip; + } + + bool is_on_a = segment_ratio<CalcT>(dist_a1_ip, dist_a1_a2).on_segment(); + // NOTE: not fully consistent with equals_point_point() since radians are always used. + bool is_on_a1 = math::equals(lon, a1_lon) && math::equals(lat, a1_lat); + bool is_on_a2 = math::equals(lon, a2_lon) && math::equals(lat, a2_lat); + + if (! (is_on_a || is_on_a1 || is_on_a2)) + { + return false; + } + + inverse_result const res_b1_ip = inverse_dist_azi::apply(b1_lon, b1_lat, lon, lat, spheroid); + dist_b1_ip = res_b1_ip.distance; + if (! same_direction(res_b1_ip.azimuth, res_b1_b2.azimuth)) + { + dist_b1_ip = -dist_b1_ip; + } + + bool is_on_b = segment_ratio<CalcT>(dist_b1_ip, dist_b1_b2).on_segment(); + // NOTE: not fully consistent with equals_point_point() since radians are always used. + bool is_on_b1 = math::equals(lon, b1_lon) && math::equals(lat, b1_lat); + bool is_on_b2 = math::equals(lon, b2_lon) && math::equals(lat, b2_lat); + + if (! (is_on_b || is_on_b1 || is_on_b2)) + { + return false; + } + + ip_flag = ipi_inters; + + if (is_on_b1) + { + lon = b1_lon; + lat = b1_lat; + dist_b1_ip = 0; + ip_flag = ipi_at_b1; + } + else if (is_on_b2) + { + lon = b2_lon; + lat = b2_lat; + dist_b1_ip = res_b1_b2.distance; + ip_flag = ipi_at_b2; + } + + if (is_on_a1) + { + lon = a1_lon; + lat = a1_lat; + dist_a1_ip = 0; + ip_flag = ipi_at_a1; + } + else if (is_on_a2) + { + lon = a2_lon; + lat = a2_lat; + dist_a1_ip = res_a1_a2.distance; + ip_flag = ipi_at_a2; + } + + return true; + } + + template <typename CalcT, typename P1, typename P2> + static inline bool is_endpoint_equal(CalcT const& dist, + P1 const& ai, P2 const& b1, P2 const& b2) + { + using geometry::detail::equals::equals_point_point; + return is_near(dist) && (equals_point_point(ai, b1) || equals_point_point(ai, b2)); + } + + template <typename CalcT> + static inline bool is_near(CalcT const& dist) + { + // NOTE: This strongly depends on the Inverse method + CalcT const small_number = CalcT(boost::is_same<CalcT, float>::value ? 0.0001 : 0.00000001); + return math::abs(dist) <= small_number; + } + + template <typename ProjCoord1, typename ProjCoord2> + static inline int position_value(ProjCoord1 const& ca1, + ProjCoord2 const& cb1, + ProjCoord2 const& cb2) + { + // S1x 0 1 2 3 4 + // S2 |----------> + return math::equals(ca1, cb1) ? 1 + : math::equals(ca1, cb2) ? 3 + : cb1 < cb2 ? + ( ca1 < cb1 ? 0 + : ca1 > cb2 ? 4 + : 2 ) + : ( ca1 > cb1 ? 0 + : ca1 < cb2 ? 4 + : 2 ); + } + + template <typename CalcT> + static inline bool same_direction(CalcT const& azimuth1, CalcT const& azimuth2) + { + // distance between two angles normalized to (-180, 180] + CalcT const angle_diff = math::longitude_distance_signed<radian>(azimuth1, azimuth2); + return math::abs(angle_diff) <= math::half_pi<CalcT>(); + } + + template <int Which> + static inline void sides_reverse_segment(side_info & sides) + { + // names assuming segment A is reversed (Which == 0) + int a1_wrt_b = sides.template get<Which, 0>(); + int a2_wrt_b = sides.template get<Which, 1>(); + std::swap(a1_wrt_b, a2_wrt_b); + sides.template set<Which>(a1_wrt_b, a2_wrt_b); + int b1_wrt_a = sides.template get<1 - Which, 0>(); + int b2_wrt_a = sides.template get<1 - Which, 1>(); + sides.template set<1 - Which>(-b1_wrt_a, -b2_wrt_a); + } + + static inline void ip_flag_reverse_segment(intersection_point_flag & ip_flag, + intersection_point_flag const& ipi_at_p1, + intersection_point_flag const& ipi_at_p2) + { + ip_flag = ip_flag == ipi_at_p1 ? ipi_at_p2 : + ip_flag == ipi_at_p2 ? ipi_at_p1 : + ip_flag; + } + + template <typename CalcT, typename SpheroidT> + static inline srs::spheroid<CalcT> normalized_spheroid(SpheroidT const& spheroid) + { + return srs::spheroid<CalcT>(CalcT(1), + CalcT(get_radius<2>(spheroid)) // b/a + / CalcT(get_radius<0>(spheroid))); + } + +private: + Spheroid m_spheroid; +}; + + +}} // namespace strategy::intersection + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_HPP diff --git a/boost/geometry/strategies/geographic/intersection_elliptic.hpp b/boost/geometry/strategies/geographic/intersection_elliptic.hpp new file mode 100644 index 0000000000..76e9940fe3 --- /dev/null +++ b/boost/geometry/strategies/geographic/intersection_elliptic.hpp @@ -0,0 +1,243 @@ +// Boost.Geometry + +// Copyright (c) 2016-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_GEOGRAPHIC_INTERSECTION_ELLIPTIC_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_ELLIPTIC_HPP + + +#include <boost/geometry/core/srs.hpp> + +#include <boost/geometry/formulas/geographic.hpp> + +#include <boost/geometry/strategies/spherical/intersection.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace intersection +{ + +template <typename Spheroid> +struct great_elliptic_segments_calc_policy + : spherical_segments_calc_policy +{ + explicit great_elliptic_segments_calc_policy(Spheroid const& spheroid = Spheroid()) + : m_spheroid(spheroid) + {} + + template <typename Point, typename Point3d> + Point from_cart3d(Point3d const& point_3d) const + { + return formula::cart3d_to_geo<Point>(point_3d, m_spheroid); + } + + template <typename Point3d, typename Point> + Point3d to_cart3d(Point const& point) const + { + return formula::geo_to_cart3d<Point3d>(point, m_spheroid); + } + + // relate_xxx_calc_policy must live londer than plane because it contains + // Spheroid object and plane keeps the reference to that object. + template <typename Point3d> + struct plane + { + typedef typename coordinate_type<Point3d>::type coord_t; + + // not normalized + plane(Point3d const& p1, Point3d const& p2) + : normal(cross_product(p1, p2)) + {} + + int side_value(Point3d const& pt) const + { + return formula::sph_side_value(normal, pt); + } + + coord_t cos_angle_between(Point3d const& p1, Point3d const& p2) const + { + Point3d v1 = p1; + detail::vec_normalize(v1); + Point3d v2 = p2; + detail::vec_normalize(v2); + + return dot_product(v1, v2); + } + + coord_t cos_angle_between(Point3d const& p1, Point3d const& p2, bool & is_forward) const + { + coord_t const c0 = 0; + + Point3d v1 = p1; + detail::vec_normalize(v1); + Point3d v2 = p2; + detail::vec_normalize(v2); + + is_forward = dot_product(normal, cross_product(v1, v2)) >= c0; + return dot_product(v1, v2); + } + + Point3d normal; + }; + + template <typename Point3d> + plane<Point3d> get_plane(Point3d const& p1, Point3d const& p2) const + { + return plane<Point3d>(p1, p2); + } + + template <typename Point3d> + bool intersection_points(plane<Point3d> const& plane1, + plane<Point3d> const& plane2, + Point3d & ip1, Point3d & ip2) const + { + typedef typename coordinate_type<Point3d>::type coord_t; + + Point3d id = cross_product(plane1.normal, plane2.normal); + // NOTE: the length should be greater than 0 at this point + // NOTE: no need to normalize in this case + + ip1 = formula::projected_to_surface(id, m_spheroid); + + ip2 = ip1; + multiply_value(ip2, coord_t(-1)); + + return true; + } + +private: + Spheroid m_spheroid; +}; + +template <typename Spheroid> +struct experimental_elliptic_segments_calc_policy +{ + explicit experimental_elliptic_segments_calc_policy(Spheroid const& spheroid = Spheroid()) + : m_spheroid(spheroid) + {} + + template <typename Point, typename Point3d> + Point from_cart3d(Point3d const& point_3d) const + { + return formula::cart3d_to_geo<Point>(point_3d, m_spheroid); + } + + template <typename Point3d, typename Point> + Point3d to_cart3d(Point const& point) const + { + return formula::geo_to_cart3d<Point3d>(point, m_spheroid); + } + + // relate_xxx_calc_policy must live londer than plane because it contains + // Spheroid object and plane keeps the reference to that object. + template <typename Point3d> + struct plane + { + typedef typename coordinate_type<Point3d>::type coord_t; + + // not normalized + plane(Point3d const& p1, Point3d const& p2, Spheroid const& spheroid) + : m_spheroid(spheroid) + { + formula::experimental_elliptic_plane(p1, p2, origin, normal, m_spheroid); + } + + int side_value(Point3d const& pt) const + { + return formula::elliptic_side_value(origin, normal, pt); + } + + coord_t cos_angle_between(Point3d const& p1, Point3d const& p2) const + { + Point3d const v1 = normalized_vec(p1); + Point3d const v2 = normalized_vec(p2); + return dot_product(v1, v2); + } + + coord_t cos_angle_between(Point3d const& p1, Point3d const& p2, bool & is_forward) const + { + coord_t const c0 = 0; + + Point3d const v1 = normalized_vec(p1); + Point3d const v2 = normalized_vec(p2); + + is_forward = dot_product(normal, cross_product(v1, v2)) >= c0; + return dot_product(v1, v2); + } + + Point3d origin; + Point3d normal; + + private: + Point3d normalized_vec(Point3d const& p) const + { + Point3d v = p; + subtract_point(v, origin); + detail::vec_normalize(v); + return v; + } + + Spheroid const& m_spheroid; + }; + + template <typename Point3d> + plane<Point3d> get_plane(Point3d const& p1, Point3d const& p2) const + { + return plane<Point3d>(p1, p2, m_spheroid); + } + + template <typename Point3d> + bool intersection_points(plane<Point3d> const& plane1, + plane<Point3d> const& plane2, + Point3d & ip1, Point3d & ip2) const + { + return formula::planes_spheroid_intersection(plane1.origin, plane1.normal, + plane2.origin, plane2.normal, + ip1, ip2, m_spheroid); + } + +private: + Spheroid m_spheroid; +}; + + +template +< + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +struct great_elliptic_segments + : ecef_segments + < + great_elliptic_segments_calc_policy<Spheroid>, + CalculationType + > +{}; + +template +< + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +struct experimental_elliptic_segments + : ecef_segments + < + experimental_elliptic_segments_calc_policy<Spheroid>, + CalculationType + > +{}; + + +}} // namespace strategy::intersection + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_ELLIPTIC_HPP diff --git a/boost/geometry/strategies/geographic/mapping_ssf.hpp b/boost/geometry/strategies/geographic/mapping_ssf.hpp index 3beedc7809..20a0523616 100644 --- a/boost/geometry/strategies/geographic/mapping_ssf.hpp +++ b/boost/geometry/strategies/geographic/mapping_ssf.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2011-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014. -// Modifications copyright (c) 2014 Oracle and/or its affiliates. +// 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 @@ -136,7 +136,7 @@ public : {} template <typename P1, typename P2, typename P> - inline int apply(P1 const& p1, P2 const& p2, P const& p) + inline int apply(P1 const& p1, P2 const& p2, P const& p) const { typedef typename promote_floating_point < diff --git a/boost/geometry/strategies/geographic/parameters.hpp b/boost/geometry/strategies/geographic/parameters.hpp new file mode 100644 index 0000000000..5638db50fa --- /dev/null +++ b/boost/geometry/strategies/geographic/parameters.hpp @@ -0,0 +1,117 @@ +// Boost.Geometry + +// 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_GEOGRAPHIC_PARAMETERS_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_PARAMETERS_HPP + + +#include <boost/geometry/formulas/andoyer_inverse.hpp> +#include <boost/geometry/formulas/thomas_inverse.hpp> +#include <boost/geometry/formulas/vincenty_inverse.hpp> + +#include <boost/mpl/assert.hpp> +#include <boost/mpl/integral_c.hpp> + + +namespace boost { namespace geometry { namespace strategy +{ + +struct andoyer +{ + template + < + typename CT, + bool EnableDistance, + bool EnableAzimuth, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false + > + struct inverse + : formula::andoyer_inverse + < + CT, EnableDistance, + EnableAzimuth, EnableReverseAzimuth, + EnableReducedLength, EnableGeodesicScale + > + {}; +}; + +struct thomas +{ + template + < + typename CT, + bool EnableDistance, + bool EnableAzimuth, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false + > + struct inverse + : formula::thomas_inverse + < + CT, EnableDistance, + EnableAzimuth, EnableReverseAzimuth, + EnableReducedLength, EnableGeodesicScale + > + {}; +}; + +struct vincenty +{ + template + < + typename CT, + bool EnableDistance, + bool EnableAzimuth, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false + > + struct inverse + : formula::vincenty_inverse + < + CT, EnableDistance, + EnableAzimuth, EnableReverseAzimuth, + EnableReducedLength, EnableGeodesicScale + > + {}; +}; + + +template <typename FormulaPolicy> +struct default_order +{ + BOOST_MPL_ASSERT_MSG + ( + false, NOT_IMPLEMENTED_FOR_THIS_TYPE + , (types<FormulaPolicy>) + ); +}; + +template<> +struct default_order<andoyer> + : boost::mpl::integral_c<unsigned int, 1> +{}; + +template<> +struct default_order<thomas> + : boost::mpl::integral_c<unsigned int, 2> +{}; + +template<> +struct default_order<vincenty> + : boost::mpl::integral_c<unsigned int, 4> +{}; + +}}} // namespace boost::geometry::strategy + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_PARAMETERS_HPP diff --git a/boost/geometry/strategies/geographic/side.hpp b/boost/geometry/strategies/geographic/side.hpp new file mode 100644 index 0000000000..0d40e4da20 --- /dev/null +++ b/boost/geometry/strategies/geographic/side.hpp @@ -0,0 +1,113 @@ +// Boost.Geometry + +// 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. + +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_HPP + +#include <boost/geometry/core/cs.hpp> +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/radian_access.hpp> +#include <boost/geometry/core/radius.hpp> +#include <boost/geometry/core/srs.hpp> + +#include <boost/geometry/formulas/spherical.hpp> + +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/promote_floating_point.hpp> +#include <boost/geometry/util/select_calculation_type.hpp> + +#include <boost/geometry/strategies/geographic/parameters.hpp> +#include <boost/geometry/strategies/side.hpp> +//#include <boost/geometry/strategies/concepts/side_concept.hpp> + + +namespace boost { namespace geometry +{ + + +namespace strategy { namespace side +{ + + +/*! +\brief Check at which side of a segment a point lies + left of segment (> 0), right of segment (< 0), on segment (0) +\ingroup strategies +\tparam FormulaPolicy Geodesic solution formula policy. +\tparam Spheroid Reference model of coordinate system. +\tparam CalculationType \tparam_calculation + */ +template +< + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +class geographic +{ +public: + geographic() + {} + + explicit geographic(Spheroid const& model) + : m_model(model) + {} + + template <typename P1, typename P2, typename P> + inline int apply(P1 const& p1, P2 const& p2, P const& p) const + { + typedef typename promote_floating_point + < + typename select_calculation_type_alt + < + CalculationType, + P1, P2, P + >::type + >::type calc_t; + + typedef typename FormulaPolicy::template inverse + <calc_t, false, true, false, false, false> inverse_formula; + + calc_t a1p = azimuth<calc_t, inverse_formula>(p1, p, m_model); + calc_t a12 = azimuth<calc_t, inverse_formula>(p1, p2, m_model); + + return formula::azimuth_side_value(a1p, a12); + } + +private: + template <typename ResultType, + typename InverseFormulaType, + typename Point1, + typename Point2, + typename ModelT> + static inline ResultType azimuth(Point1 const& point1, Point2 const& point2, + ModelT const& model) + { + return InverseFormulaType::apply(get_as_radian<0>(point1), + get_as_radian<1>(point1), + get_as_radian<0>(point2), + get_as_radian<1>(point2), + model).azimuth; + } + + Spheroid m_model; +}; + + +}} // namespace strategy::side + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_HPP diff --git a/boost/geometry/strategies/geographic/side_andoyer.hpp b/boost/geometry/strategies/geographic/side_andoyer.hpp index c3e71cd1cd..204e45f6e2 100644 --- a/boost/geometry/strategies/geographic/side_andoyer.hpp +++ b/boost/geometry/strategies/geographic/side_andoyer.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014, 2015, 2016. -// Modifications copyright (c) 2014-2016 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 @@ -15,9 +15,7 @@ #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_ANDOYER_HPP -#include <boost/geometry/formulas/andoyer_inverse.hpp> - -#include <boost/geometry/strategies/geographic/side_detail.hpp> +#include <boost/geometry/strategies/geographic/side.hpp> namespace boost { namespace geometry @@ -31,17 +29,24 @@ namespace strategy { namespace side \brief Check at which side of a segment a point lies left of segment (> 0), right of segment (< 0), on segment (0) \ingroup strategies -\tparam Model Reference model of coordinate system. +\tparam Spheroid Reference model of coordinate system. \tparam CalculationType \tparam_calculation */ -template <typename Model, typename CalculationType = void> +template +< + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> class andoyer - : public detail::by_azimuth<geometry::formula::andoyer_inverse, Model, CalculationType> + : public side::geographic<strategy::andoyer, Spheroid, CalculationType> { - typedef detail::by_azimuth<geometry::formula::andoyer_inverse, Model, CalculationType> base_t; + typedef side::geographic<strategy::andoyer, Spheroid, CalculationType> base_t; public: - andoyer(Model const& model = Model()) + andoyer() + {} + + explicit andoyer(Spheroid const& model) : base_t(model) {} }; diff --git a/boost/geometry/strategies/geographic/side_detail.hpp b/boost/geometry/strategies/geographic/side_detail.hpp deleted file mode 100644 index ce1b47c88e..0000000000 --- a/boost/geometry/strategies/geographic/side_detail.hpp +++ /dev/null @@ -1,139 +0,0 @@ -// Boost.Geometry - -// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. - -// This file was modified by Oracle on 2014, 2015, 2016. -// Modifications copyright (c) 2014-2016 Oracle and/or its affiliates. - -// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle - -// Use, modification and distribution is subject to the Boost Software License, -// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) - -#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_DETAIL_HPP -#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_DETAIL_HPP - -#include <boost/geometry/core/cs.hpp> -#include <boost/geometry/core/access.hpp> -#include <boost/geometry/core/radian_access.hpp> -#include <boost/geometry/core/radius.hpp> - -#include <boost/geometry/util/math.hpp> -#include <boost/geometry/util/promote_floating_point.hpp> -#include <boost/geometry/util/select_calculation_type.hpp> - -#include <boost/geometry/strategies/side.hpp> -//#include <boost/geometry/strategies/concepts/side_concept.hpp> - - -namespace boost { namespace geometry -{ - - -namespace strategy { namespace side -{ - -#ifndef DOXYGEN_NO_DETAIL -namespace detail -{ - -/*! -\brief Check at which side of a segment a point lies - left of segment (> 0), right of segment (< 0), on segment (0) -\ingroup strategies -\tparam InverseFormula Geodesic inverse solution formula. -\tparam Model Reference model of coordinate system. -\tparam CalculationType \tparam_calculation - */ -template <template<typename, bool, bool, bool, bool, bool> class InverseFormula, - typename Model, - typename CalculationType = void> -class by_azimuth -{ -public: - by_azimuth(Model const& model = Model()) - : m_model(model) - {} - - template <typename P1, typename P2, typename P> - inline int apply(P1 const& p1, P2 const& p2, P const& p) - { - typedef typename promote_floating_point - < - typename select_calculation_type_alt - < - CalculationType, - P1, P2, P - >::type - >::type calc_t; - - typedef InverseFormula<calc_t, false, true, false, false, false> inverse_formula; - - calc_t a1p = azimuth<calc_t, inverse_formula>(p1, p, m_model); - calc_t a12 = azimuth<calc_t, inverse_formula>(p1, p2, m_model); - - calc_t const pi = math::pi<calc_t>(); - - // instead of the formula from XTD - //calc_t a_diff = asin(sin(a1p - a12)); - - calc_t a_diff = a1p - a12; - // normalize, angle in [-pi, pi] - while ( a_diff > pi ) - a_diff -= calc_t(2) * pi; - while ( a_diff < -pi ) - a_diff += calc_t(2) * pi; - - // NOTE: in general it shouldn't be required to support the pi/-pi case - // because in non-cartesian systems it makes sense to check the side - // only "between" the endpoints. - // However currently the winding strategy calls the side strategy - // for vertical segments to check if the point is "between the endpoints. - // This could be avoided since the side strategy is not required for that - // because meridian is the shortest path. So a difference of - // longitudes would be sufficient (of course normalized to [-pi, pi]). - - // NOTE: with the above said, the pi/-pi check is temporary - // however in case if this was required - // the geodesics on ellipsoid aren't "symmetrical" - // therefore instead of comparing a_diff to pi and -pi - // one should probably use inverse azimuths and compare - // the difference to 0 as well - - // positive azimuth is on the right side - return math::equals(a_diff, 0) - || math::equals(a_diff, pi) - || math::equals(a_diff, -pi) ? 0 - : a_diff > 0 ? -1 // right - : 1; // left - } - -private: - template <typename ResultType, - typename InverseFormulaType, - typename Point1, - typename Point2, - typename ModelT> - static inline ResultType azimuth(Point1 const& point1, Point2 const& point2, ModelT const& model) - { - return InverseFormulaType::apply(get_as_radian<0>(point1), - get_as_radian<1>(point1), - get_as_radian<0>(point2), - get_as_radian<1>(point2), - model).azimuth; - } - - Model m_model; -}; - -} // detail -#endif // DOXYGEN_NO_DETAIL - -}} // namespace strategy::side - - -}} // namespace boost::geometry - - -#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_DETAIL_HPP diff --git a/boost/geometry/strategies/geographic/side_thomas.hpp b/boost/geometry/strategies/geographic/side_thomas.hpp index 96b0323307..e6f8d77b58 100644 --- a/boost/geometry/strategies/geographic/side_thomas.hpp +++ b/boost/geometry/strategies/geographic/side_thomas.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014, 2015, 2016. -// Modifications copyright (c) 2014-2016 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 @@ -15,9 +15,7 @@ #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_THOMAS_HPP -#include <boost/geometry/formulas/thomas_inverse.hpp> - -#include <boost/geometry/strategies/geographic/side_detail.hpp> +#include <boost/geometry/strategies/geographic/side.hpp> namespace boost { namespace geometry @@ -31,17 +29,24 @@ namespace strategy { namespace side \brief Check at which side of a segment a point lies left of segment (> 0), right of segment (< 0), on segment (0) \ingroup strategies -\tparam Model Reference model of coordinate system. +\tparam Spheroid Reference model of coordinate system. \tparam CalculationType \tparam_calculation */ -template <typename Model, typename CalculationType = void> +template +< + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> class thomas - : public detail::by_azimuth<geometry::formula::thomas_inverse, Model, CalculationType> + : public side::geographic<strategy::thomas, Spheroid, CalculationType> { - typedef detail::by_azimuth<geometry::formula::thomas_inverse, Model, CalculationType> base_t; + typedef side::geographic<strategy::thomas, Spheroid, CalculationType> base_t; public: - thomas(Model const& model = Model()) + thomas() + {} + + explicit thomas(Spheroid const& model) : base_t(model) {} }; diff --git a/boost/geometry/strategies/geographic/side_vincenty.hpp b/boost/geometry/strategies/geographic/side_vincenty.hpp index 103277a8bd..b2f51b0901 100644 --- a/boost/geometry/strategies/geographic/side_vincenty.hpp +++ b/boost/geometry/strategies/geographic/side_vincenty.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014, 2015, 2016. -// Modifications copyright (c) 2014-2016 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 @@ -15,9 +15,7 @@ #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_VINCENTY_HPP -#include <boost/geometry/formulas/vincenty_inverse.hpp> - -#include <boost/geometry/strategies/geographic/side_detail.hpp> +#include <boost/geometry/strategies/geographic/side.hpp> namespace boost { namespace geometry @@ -31,17 +29,24 @@ namespace strategy { namespace side \brief Check at which side of a segment a point lies left of segment (> 0), right of segment (< 0), on segment (0) \ingroup strategies -\tparam Model Reference model of coordinate system. +\tparam Spheroid Reference model of coordinate system. \tparam CalculationType \tparam_calculation */ -template <typename Model, typename CalculationType = void> +template +< + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> class vincenty - : public detail::by_azimuth<geometry::formula::vincenty_inverse, Model, CalculationType> + : public side::geographic<strategy::vincenty, Spheroid, CalculationType> { - typedef detail::by_azimuth<geometry::formula::vincenty_inverse, Model, CalculationType> base_t; + typedef side::geographic<strategy::vincenty, Spheroid, CalculationType> base_t; public: - vincenty(Model const& model = Model()) + vincenty() + {} + + explicit vincenty(Spheroid const& model) : base_t(model) {} }; |