diff options
Diffstat (limited to 'boost/geometry/strategies')
44 files changed, 4106 insertions, 1162 deletions
diff --git a/boost/geometry/strategies/agnostic/point_in_point.hpp b/boost/geometry/strategies/agnostic/point_in_point.hpp index e4f9bec4c6..c6e2b63623 100644 --- a/boost/geometry/strategies/agnostic/point_in_point.hpp +++ b/boost/geometry/strategies/agnostic/point_in_point.hpp @@ -1,12 +1,13 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2014 Oracle and/or its affiliates. +// 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) -// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle #ifndef BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POINT_HPP #define BOOST_GEOMETRY_STRATEGY_AGNOSTIC_POINT_IN_POINT_HPP @@ -41,41 +42,12 @@ struct point_in_point namespace services { -template <typename Point1, typename Point2> -struct default_strategy<point_tag, point_tag, point_tag, point_tag, cartesian_tag, cartesian_tag, Point1, Point2> -{ - typedef strategy::within::point_in_point<Point1, Point2> type; -}; - -template <typename Point1, typename Point2> -struct default_strategy<point_tag, point_tag, point_tag, point_tag, spherical_tag, spherical_tag, Point1, Point2> -{ - typedef strategy::within::point_in_point<Point1, Point2> type; -}; - -template <typename Point1, typename Point2, typename AnyCS1, typename AnyCS2> -struct default_strategy<point_tag, point_tag, point_tag, point_tag, AnyCS1, AnyCS2, Point1, Point2> +template <typename Point, typename PointLike, typename Tag2, typename AnyCS1, typename AnyCS2> +struct default_strategy<Point, PointLike, point_tag, Tag2, pointlike_tag, pointlike_tag, AnyCS1, AnyCS2> { - typedef strategy::within::point_in_point<Point1, Point2> type; + typedef strategy::within::point_in_point<Point, typename point_type<PointLike>::type> type; }; -template <typename Point, typename MultiPoint> -struct default_strategy<point_tag, multi_point_tag, point_tag, multi_point_tag, cartesian_tag, cartesian_tag, Point, MultiPoint> -{ - typedef strategy::within::point_in_point<Point, typename point_type<MultiPoint>::type> type; -}; - -template <typename Point, typename MultiPoint> -struct default_strategy<point_tag, multi_point_tag, point_tag, multi_point_tag, spherical_tag, spherical_tag, Point, MultiPoint> -{ - typedef strategy::within::point_in_point<Point, typename point_type<MultiPoint>::type> type; -}; - -template <typename Point, typename MultiPoint, typename AnyCS1, typename AnyCS2> -struct default_strategy<point_tag, multi_point_tag, point_tag, multi_point_tag, AnyCS1, AnyCS2, Point, MultiPoint> -{ - typedef strategy::within::point_in_point<Point, typename point_type<MultiPoint>::type> type; -}; } // namespace services @@ -85,45 +57,14 @@ struct default_strategy<point_tag, multi_point_tag, point_tag, multi_point_tag, }} // namespace strategy::within - #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS namespace strategy { namespace covered_by { namespace services { -template <typename Point1, typename Point2> -struct default_strategy<point_tag, point_tag, point_tag, point_tag, cartesian_tag, cartesian_tag, Point1, Point2> -{ - typedef strategy::within::point_in_point<Point1, Point2> type; -}; - -template <typename Point1, typename Point2> -struct default_strategy<point_tag, point_tag, point_tag, point_tag, spherical_tag, spherical_tag, Point1, Point2> -{ - typedef strategy::within::point_in_point<Point1, Point2> type; -}; - -template <typename Point1, typename Point2, typename AnyCS1, typename AnyCS2> -struct default_strategy<point_tag, point_tag, point_tag, point_tag, AnyCS1, AnyCS2, Point1, Point2> -{ - typedef strategy::within::point_in_point<Point1, Point2> type; -}; - -template <typename Point, typename MultiPoint> -struct default_strategy<point_tag, multi_point_tag, point_tag, multi_point_tag, cartesian_tag, cartesian_tag, Point, MultiPoint> -{ - typedef strategy::within::point_in_point<Point, typename point_type<MultiPoint>::type> type; -}; - -template <typename Point, typename MultiPoint> -struct default_strategy<point_tag, multi_point_tag, point_tag, multi_point_tag, spherical_tag, spherical_tag, Point, MultiPoint> -{ - typedef strategy::within::point_in_point<Point, typename point_type<MultiPoint>::type> type; -}; - -template <typename Point, typename MultiPoint, typename AnyCS1, typename AnyCS2> -struct default_strategy<point_tag, multi_point_tag, point_tag, multi_point_tag, AnyCS1, AnyCS2, Point, MultiPoint> +template <typename Point, typename PointLike, typename Tag2, typename AnyCS1, typename AnyCS2> +struct default_strategy<Point, PointLike, point_tag, Tag2, pointlike_tag, pointlike_tag, AnyCS1, AnyCS2> { - typedef strategy::within::point_in_point<Point, typename point_type<MultiPoint>::type> type; + typedef strategy::within::point_in_point<Point, typename point_type<PointLike>::type> type; }; }}} // namespace strategy::covered_by::services diff --git a/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp b/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp index 9e2ec2c4ff..3e7e258788 100644 --- a/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp +++ b/boost/geometry/strategies/agnostic/point_in_poly_winding.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. // Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland. -// This file was modified by Oracle on 2013, 2014, 2016. -// Modifications copyright (c) 2013-2016 Oracle and/or its affiliates. +// This file was modified by Oracle on 2013, 2014, 2016, 2017. +// Modifications copyright (c) 2013-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 @@ -314,6 +314,7 @@ struct winding_calculate_count<Point, CalculationType, cartesian_tag> \ingroup strategies \tparam Point \tparam_point \tparam PointOfSegment \tparam_segment_point +\tparam SideStrategy Side strategy \tparam CalculationType \tparam_calculation \author Barend Gehrels \note The implementation is inspired by terralib http://www.terralib.org (LGPL) @@ -329,6 +330,10 @@ template < typename Point, typename PointOfSegment = Point, + typename SideStrategy = typename strategy::side::services::default_strategy + < + typename cs_tag<Point>::type + >::type, typename CalculationType = void > class winding @@ -339,14 +344,7 @@ class winding PointOfSegment, CalculationType >::type calculation_type; - - - typedef typename strategy::side::services::default_strategy - < - typename cs_tag<Point>::type - >::type strategy_side_type; - - + /*! subclass to keep state */ class counter { @@ -371,7 +369,6 @@ class winding }; - static inline int check_segment(Point const& point, PointOfSegment const& seg1, PointOfSegment const& seg2, counter& state, bool& eq1, bool& eq2) @@ -390,16 +387,22 @@ class winding } -public : +public: + winding() + {} + + explicit winding(SideStrategy const& side_strategy) + : m_side_strategy(side_strategy) + {} // Typedefs and static methods to fulfill the concept typedef Point point_type; typedef PointOfSegment segment_point_type; typedef counter state_type; - static inline bool apply(Point const& point, - PointOfSegment const& s1, PointOfSegment const& s2, - counter& state) + inline bool apply(Point const& point, + PointOfSegment const& s1, PointOfSegment const& s2, + counter& state) const { typedef typename cs_tag<Point>::type cs_t; @@ -418,7 +421,7 @@ public : else // count == 2 || count == -2 { // 1 left, -1 right - side = strategy_side_type::apply(s1, s2, point); + side = m_side_strategy.apply(s1, s2, point); } if (side == 0) @@ -445,6 +448,9 @@ public : { return state.code(); } + +private: + SideStrategy m_side_strategy; }; @@ -453,29 +459,26 @@ public : namespace services { -// Register using "areal_tag" for ring, polygon, multi-polygon -template <typename AnyTag, typename Point, typename Geometry> -struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, cartesian_tag, cartesian_tag, Point, Geometry> +template <typename Point, typename Geometry, typename AnyTag> +struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, polygonal_tag, cartesian_tag, cartesian_tag> { typedef winding<Point, typename geometry::point_type<Geometry>::type> type; }; -template <typename AnyTag, typename Point, typename Geometry> -struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, spherical_tag, spherical_tag, Point, Geometry> +template <typename Point, typename Geometry, typename AnyTag> +struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag> { typedef winding<Point, typename geometry::point_type<Geometry>::type> type; }; -// TODO: use linear_tag and pointlike_tag the same way how areal_tag is used - template <typename Point, typename Geometry, typename AnyTag> -struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, cartesian_tag, cartesian_tag, Point, Geometry> +struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, linear_tag, cartesian_tag, cartesian_tag> { typedef winding<Point, typename geometry::point_type<Geometry>::type> type; }; template <typename Point, typename Geometry, typename AnyTag> -struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, spherical_tag, spherical_tag, Point, Geometry> +struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, linear_tag, spherical_tag, spherical_tag> { typedef winding<Point, typename geometry::point_type<Geometry>::type> type; }; @@ -488,34 +491,30 @@ struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, spherical_tag, sph }} // namespace strategy::within - #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS namespace strategy { namespace covered_by { namespace services { -// Register using "areal_tag" for ring, polygon, multi-polygon -template <typename AnyTag, typename Point, typename Geometry> -struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, cartesian_tag, cartesian_tag, Point, Geometry> +template <typename Point, typename Geometry, typename AnyTag> +struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, polygonal_tag, cartesian_tag, cartesian_tag> { typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type; }; -template <typename AnyTag, typename Point, typename Geometry> -struct default_strategy<point_tag, AnyTag, point_tag, areal_tag, spherical_tag, spherical_tag, Point, Geometry> +template <typename Point, typename Geometry, typename AnyTag> +struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, polygonal_tag, spherical_tag, spherical_tag> { typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type; }; -// TODO: use linear_tag and pointlike_tag the same way how areal_tag is used - template <typename Point, typename Geometry, typename AnyTag> -struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, cartesian_tag, cartesian_tag, Point, Geometry> +struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, linear_tag, cartesian_tag, cartesian_tag> { typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type; }; template <typename Point, typename Geometry, typename AnyTag> -struct default_strategy<point_tag, AnyTag, point_tag, AnyTag, spherical_tag, spherical_tag, Point, Geometry> +struct default_strategy<Point, Geometry, point_tag, AnyTag, pointlike_tag, linear_tag, spherical_tag, spherical_tag> { typedef strategy::within::winding<Point, typename geometry::point_type<Geometry>::type> type; }; diff --git a/boost/geometry/strategies/agnostic/relate.hpp b/boost/geometry/strategies/agnostic/relate.hpp deleted file mode 100644 index 676207694f..0000000000 --- a/boost/geometry/strategies/agnostic/relate.hpp +++ /dev/null @@ -1,122 +0,0 @@ -// Boost.Geometry (aka GGL, Generic Geometry Library) - -// Copyright (c) 2014-2015 Oracle and/or its affiliates. - -// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle - -// Use, modification and distribution is subject to the Boost Software License, -// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) - -#ifndef BOOST_GEOMETRY_STRATEGY_AGNOSTIC_RELATE_HPP -#define BOOST_GEOMETRY_STRATEGY_AGNOSTIC_RELATE_HPP - -#include <boost/geometry/algorithms/relate.hpp> - - -namespace boost { namespace geometry -{ - -namespace strategy { namespace relate -{ - -template <typename Geometry1, typename Geometry2, typename StaticMask> -struct relate -{ - static inline bool apply(Geometry1 const& geometry1, Geometry2 const& geometry2) - { - return geometry::relate(geometry1, geometry2, StaticMask()); - } -}; - -} // namespace relate - -namespace within -{ - -#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - -namespace services -{ - - -template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2, typename AnyCS> -struct default_strategy<AnyTag1, AnyTag2, AnyTag1, AnyTag2, AnyCS, AnyCS, Geometry1, Geometry2> -{ - typedef strategy::relate::relate - < - Geometry1, - Geometry2, - typename detail::de9im::static_mask_within_type - < - Geometry1, Geometry2 - >::type - > type; -}; - -template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2, typename AnyCS> -struct default_strategy<AnyTag1, AnyTag2, AnyTag1, areal_tag, AnyCS, AnyCS, Geometry1, Geometry2> -{ - typedef strategy::relate::relate - < - Geometry1, - Geometry2, - typename detail::de9im::static_mask_within_type - < - Geometry1, Geometry2 - >::type - > type; -}; - - -} // namespace services - -#endif - - -}} // namespace strategy::within - - - -#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS -namespace strategy { namespace covered_by { namespace services -{ - - -template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2, typename AnyCS> -struct default_strategy<AnyTag1, AnyTag2, AnyTag1, AnyTag2, AnyCS, AnyCS, Geometry1, Geometry2> -{ - typedef strategy::relate::relate - < - Geometry1, - Geometry2, - typename detail::de9im::static_mask_covered_by_type - < - Geometry1, Geometry2 - >::type - > type; -}; - -template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2, typename AnyCS> -struct default_strategy<AnyTag1, AnyTag2, AnyTag1, areal_tag, AnyCS, AnyCS, Geometry1, Geometry2> -{ - typedef strategy::relate::relate - < - Geometry1, - Geometry2, - typename detail::de9im::static_mask_covered_by_type - < - Geometry1, Geometry2 - >::type - > type; -}; - - -}}} // namespace strategy::covered_by::services -#endif - - -}} // namespace boost::geometry - - -#endif // BOOST_GEOMETRY_STRATEGY_AGNOSTIC_RELATE_HPP diff --git a/boost/geometry/strategies/azimuth.hpp b/boost/geometry/strategies/azimuth.hpp new file mode 100644 index 0000000000..8044d3933b --- /dev/null +++ b/boost/geometry/strategies/azimuth.hpp @@ -0,0 +1,44 @@ +// 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_AZIMUTH_HPP +#define BOOST_GEOMETRY_STRATEGIES_AZIMUTH_HPP + +#include <boost/mpl/assert.hpp> + +namespace boost { namespace geometry +{ + + +namespace strategy { namespace azimuth { namespace services +{ + +/*! +\brief Traits class binding a default azimuth strategy to a coordinate system +\ingroup util +\tparam CSTag tag of coordinate system +\tparam CalculationType \tparam_calculation +*/ +template <typename CSTag, typename CalculationType = void> +struct default_strategy +{ + BOOST_MPL_ASSERT_MSG + ( + false, NOT_IMPLEMENTED_FOR_THIS_TYPE + , (types<CSTag>) + ); +}; + +}}} // namespace strategy::azimuth::services + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_AZIMUTH_HPP diff --git a/boost/geometry/strategies/cartesian/area_surveyor.hpp b/boost/geometry/strategies/cartesian/area_surveyor.hpp index ba76f67946..b3f19b1b1e 100644 --- a/boost/geometry/strategies/cartesian/area_surveyor.hpp +++ b/boost/geometry/strategies/cartesian/area_surveyor.hpp @@ -4,8 +4,8 @@ // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. // Copyright (c) 2009-2012 Mateusz Loskot, London, UK. -// This file was modified by Oracle on 2016. -// Modifications copyright (c) 2016, Oracle and/or its affiliates. +// 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 @@ -25,6 +25,7 @@ //#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> diff --git a/boost/geometry/strategies/cartesian/azimuth.hpp b/boost/geometry/strategies/cartesian/azimuth.hpp new file mode 100644 index 0000000000..62f804e8f5 --- /dev/null +++ b/boost/geometry/strategies/cartesian/azimuth.hpp @@ -0,0 +1,49 @@ +// 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_CARTESIAN_AZIMUTH_HPP +#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AZIMUTH_HPP + +#include <boost/geometry/core/tags.hpp> + +namespace boost { namespace geometry +{ + +namespace strategy { namespace azimuth +{ + +template +< + typename CalculationType = void +> +class cartesian +{}; + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +namespace services +{ + +template <typename CalculationType> +struct default_strategy<cartesian_tag, CalculationType> +{ + typedef strategy::azimuth::cartesian<CalculationType> type; +}; + +} + +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +}} // namespace strategy::azimuth + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AZIMUTH_HPP diff --git a/boost/geometry/strategies/cartesian/box_in_box.hpp b/boost/geometry/strategies/cartesian/box_in_box.hpp index 28a6f29336..4641aafcb1 100644 --- a/boost/geometry/strategies/cartesian/box_in_box.hpp +++ b/boost/geometry/strategies/cartesian/box_in_box.hpp @@ -5,8 +5,8 @@ // Copyright (c) 2009-2015 Mateusz Loskot, London, UK. // Copyright (c) 2013-2015 Adam Wulkiewicz, Lodz, Poland. -// This file was modified by Oracle on 2015, 2016. -// Modifications copyright (c) 2016, Oracle and/or its affiliates. +// This file was modified by Oracle on 2015, 2016, 2017. +// Modifications copyright (c) 2016-2017, Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -238,10 +238,10 @@ namespace within { namespace services template <typename BoxContained, typename BoxContaining> struct default_strategy < + BoxContained, BoxContaining, box_tag, box_tag, - box_tag, areal_tag, - cartesian_tag, cartesian_tag, - BoxContained, BoxContaining + areal_tag, areal_tag, + cartesian_tag, cartesian_tag > { typedef within::box_in_box<BoxContained, BoxContaining> type; @@ -251,10 +251,10 @@ struct default_strategy template <typename BoxContained, typename BoxContaining> struct default_strategy < + BoxContained, BoxContaining, box_tag, box_tag, - box_tag, areal_tag, - spherical_tag, spherical_tag, - BoxContained, BoxContaining + areal_tag, areal_tag, + spherical_tag, spherical_tag > { typedef within::box_in_box<BoxContained, BoxContaining> type; @@ -269,10 +269,10 @@ namespace covered_by { namespace services template <typename BoxContained, typename BoxContaining> struct default_strategy < + BoxContained, BoxContaining, box_tag, box_tag, - box_tag, areal_tag, - cartesian_tag, cartesian_tag, - BoxContained, BoxContaining + areal_tag, areal_tag, + cartesian_tag, cartesian_tag > { typedef within::box_in_box @@ -286,10 +286,10 @@ struct default_strategy template <typename BoxContained, typename BoxContaining> struct default_strategy < + BoxContained, BoxContaining, box_tag, box_tag, - box_tag, areal_tag, - spherical_tag, spherical_tag, - BoxContained, BoxContaining + areal_tag, areal_tag, + spherical_tag, spherical_tag > { typedef within::box_in_box diff --git a/boost/geometry/strategies/cartesian/disjoint_segment_box.hpp b/boost/geometry/strategies/cartesian/disjoint_segment_box.hpp new file mode 100644 index 0000000000..4270803322 --- /dev/null +++ b/boost/geometry/strategies/cartesian/disjoint_segment_box.hpp @@ -0,0 +1,320 @@ +// Boost.Geometry + +// Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2008-2014 Bruno Lalande, Paris, France. +// Copyright (c) 2009-2014 Mateusz Loskot, London, UK. +// Copyright (c) 2013-2014 Adam Wulkiewicz, Lodz, Poland. + +// This file was modified by Oracle on 2013-2017. +// Modifications copyright (c) 2013-2017, Oracle and/or its affiliates. + +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle +// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_DISJOINT_SEGMENT_BOX_HPP +#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_DISJOINT_SEGMENT_BOX_HPP + + +#include <cstddef> +#include <utility> + +#include <boost/numeric/conversion/cast.hpp> + +#include <boost/geometry/util/math.hpp> +#include <boost/geometry/util/calculation_type.hpp> + +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/tags.hpp> +#include <boost/geometry/core/coordinate_dimension.hpp> +#include <boost/geometry/core/point_type.hpp> + +#include <boost/geometry/algorithms/detail/assign_indexed_point.hpp> + +#include <boost/geometry/strategies/disjoint.hpp> + + +namespace boost { namespace geometry { namespace strategy { namespace disjoint +{ + +namespace detail +{ + +template <std::size_t I> +struct compute_tmin_tmax_per_dim +{ + template <typename SegmentPoint, typename Box, typename RelativeDistance> + static inline void apply(SegmentPoint const& p0, + SegmentPoint const& p1, + Box const& box, + RelativeDistance& ti_min, + RelativeDistance& ti_max, + RelativeDistance& diff) + { + typedef typename coordinate_type<Box>::type box_coordinate_type; + typedef typename coordinate_type + < + SegmentPoint + >::type point_coordinate_type; + + RelativeDistance c_p0 = boost::numeric_cast + < + point_coordinate_type + >( geometry::get<I>(p0) ); + + RelativeDistance c_p1 = boost::numeric_cast + < + point_coordinate_type + >( geometry::get<I>(p1) ); + + RelativeDistance c_b_min = boost::numeric_cast + < + box_coordinate_type + >( geometry::get<geometry::min_corner, I>(box) ); + + RelativeDistance c_b_max = boost::numeric_cast + < + box_coordinate_type + >( geometry::get<geometry::max_corner, I>(box) ); + + if ( geometry::get<I>(p1) >= geometry::get<I>(p0) ) + { + diff = c_p1 - c_p0; + ti_min = c_b_min - c_p0; + ti_max = c_b_max - c_p0; + } + else + { + diff = c_p0 - c_p1; + ti_min = c_p0 - c_b_max; + ti_max = c_p0 - c_b_min; + } + } +}; + + +template +< + typename RelativeDistance, + typename SegmentPoint, + typename Box, + std::size_t I, + std::size_t Dimension +> +struct disjoint_segment_box_impl +{ + template <typename RelativeDistancePair> + static inline bool apply(SegmentPoint const& p0, + SegmentPoint const& p1, + Box const& box, + RelativeDistancePair& t_min, + RelativeDistancePair& t_max) + { + RelativeDistance ti_min, ti_max, diff; + + compute_tmin_tmax_per_dim<I>::apply(p0, p1, box, ti_min, ti_max, diff); + + if ( geometry::math::equals(diff, 0) ) + { + if ( (geometry::math::equals(t_min.second, 0) + && t_min.first > ti_max) + || + (geometry::math::equals(t_max.second, 0) + && t_max.first < ti_min) + || + (math::sign(ti_min) * math::sign(ti_max) > 0) ) + { + return true; + } + } + + RelativeDistance t_min_x_diff = t_min.first * diff; + RelativeDistance t_max_x_diff = t_max.first * diff; + + if ( t_min_x_diff > ti_max * t_min.second + || t_max_x_diff < ti_min * t_max.second ) + { + return true; + } + + if ( ti_min * t_min.second > t_min_x_diff ) + { + t_min.first = ti_min; + t_min.second = diff; + } + if ( ti_max * t_max.second < t_max_x_diff ) + { + t_max.first = ti_max; + t_max.second = diff; + } + + if ( t_min.first > t_min.second || t_max.first < 0 ) + { + return true; + } + + return disjoint_segment_box_impl + < + RelativeDistance, + SegmentPoint, + Box, + I + 1, + Dimension + >::apply(p0, p1, box, t_min, t_max); + } +}; + + +template +< + typename RelativeDistance, + typename SegmentPoint, + typename Box, + std::size_t Dimension +> +struct disjoint_segment_box_impl + < + RelativeDistance, SegmentPoint, Box, 0, Dimension + > +{ + static inline bool apply(SegmentPoint const& p0, + SegmentPoint const& p1, + Box const& box) + { + std::pair<RelativeDistance, RelativeDistance> t_min, t_max; + RelativeDistance diff; + + compute_tmin_tmax_per_dim<0>::apply(p0, p1, box, + t_min.first, t_max.first, diff); + + if ( geometry::math::equals(diff, 0) ) + { + if ( geometry::math::equals(t_min.first, 0) ) { t_min.first = -1; } + if ( geometry::math::equals(t_max.first, 0) ) { t_max.first = 1; } + + if (math::sign(t_min.first) * math::sign(t_max.first) > 0) + { + return true; + } + } + + if ( t_min.first > diff || t_max.first < 0 ) + { + return true; + } + + t_min.second = t_max.second = diff; + + return disjoint_segment_box_impl + < + RelativeDistance, SegmentPoint, Box, 1, Dimension + >::apply(p0, p1, box, t_min, t_max); + } +}; + + +template +< + typename RelativeDistance, + typename SegmentPoint, + typename Box, + std::size_t Dimension +> +struct disjoint_segment_box_impl + < + RelativeDistance, SegmentPoint, Box, Dimension, Dimension + > +{ + template <typename RelativeDistancePair> + static inline bool apply(SegmentPoint const&, SegmentPoint const&, + Box const&, + RelativeDistancePair&, RelativeDistancePair&) + { + return false; + } +}; + +} // namespace detail + +// NOTE: This may be temporary place for this or corresponding strategy +// It seems to be more appropriate to implement the opposite of it +// e.g. intersection::segment_box because in disjoint() algorithm +// other strategies that are used are intersection and covered_by strategies. +struct segment_box +{ + template <typename Segment, typename Box> + struct point_in_geometry_strategy + : services::default_strategy + < + typename point_type<Segment>::type, + Box + > + {}; + + template <typename Segment, typename Box> + static inline typename point_in_geometry_strategy<Segment, Box>::type + get_point_in_geometry_strategy() + { + typedef typename point_in_geometry_strategy<Segment, Box>::type strategy_type; + + return strategy_type(); + } + + template <typename Segment, typename Box> + static inline bool apply(Segment const& segment, Box const& box) + { + assert_dimension_equal<Segment, Box>(); + + typedef typename util::calculation_type::geometric::binary + < + Segment, Box, void + >::type relative_distance_type; + + typedef typename point_type<Segment>::type segment_point_type; + segment_point_type p0, p1; + geometry::detail::assign_point_from_index<0>(segment, p0); + geometry::detail::assign_point_from_index<1>(segment, p1); + + return detail::disjoint_segment_box_impl + < + relative_distance_type, segment_point_type, Box, + 0, dimension<Box>::value + >::apply(p0, p1, box); + } +}; + + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +namespace services +{ + +// Currently used in all coordinate systems + +template <typename Linear, typename Box, typename LinearTag> +struct default_strategy<Linear, Box, LinearTag, box_tag, 1, 2> +{ + typedef disjoint::segment_box type; +}; + +template <typename Box, typename Linear, typename LinearTag> +struct default_strategy<Box, Linear, box_tag, LinearTag, 2, 1> +{ + typedef disjoint::segment_box type; +}; + + +} // namespace services + + +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}}}} // namespace boost::geometry::strategy::disjoint + + +#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_DISJOINT_SEGMENT_BOX_HPP diff --git a/boost/geometry/strategies/cartesian/envelope_segment.hpp b/boost/geometry/strategies/cartesian/envelope_segment.hpp new file mode 100644 index 0000000000..0ddbf12a41 --- /dev/null +++ b/boost/geometry/strategies/cartesian/envelope_segment.hpp @@ -0,0 +1,71 @@ +// 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_CARTESIAN_ENVELOPE_SEGMENT_HPP +#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_ENVELOPE_SEGMENT_HPP + + +#include <boost/geometry/algorithms/detail/envelope/segment.hpp> +#include <boost/geometry/core/tags.hpp> +#include <boost/geometry/util/select_calculation_type.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace envelope +{ + +template +< + typename CalculationType = void +> +class cartesian_segment +{ +public : + + template <typename Point1, typename Point2, typename Box> + inline void + apply(Point1 const& point1, Point2 const& point2, Box& box) const + { + geometry::detail::envelope::envelope_one_segment + < + 0, + dimension<Point1>::value + > + ::apply(point1, + point2, + box, + strategy::envelope::cartesian_segment<CalculationType>()); + } + +}; + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +namespace services +{ + +template <typename CalculationType> +struct default_strategy<cartesian_tag, CalculationType> +{ + typedef strategy::envelope::cartesian_segment<CalculationType> type; +}; + +} + +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::envelope + +}} //namepsace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_ENVELOPE_SEGMENT_HPP diff --git a/boost/geometry/strategies/cartesian/cart_intersect.hpp b/boost/geometry/strategies/cartesian/intersection.hpp index 0cb5d75457..20b6b93367 100644 --- a/boost/geometry/strategies/cartesian/cart_intersect.hpp +++ b/boost/geometry/strategies/cartesian/intersection.hpp @@ -3,8 +3,8 @@ // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. // Copyright (c) 2013-2014 Adam Wulkiewicz, Lodz, Poland. -// This file was modified by Oracle on 2014, 2016. -// Modifications copyright (c) 2014-2016, Oracle and/or its affiliates. +// This file was modified by Oracle on 2014, 2016, 2017. +// Modifications copyright (c) 2014-2017, Oracle and/or its affiliates. // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -33,13 +33,16 @@ #include <boost/geometry/util/promote_integral.hpp> #include <boost/geometry/util/select_calculation_type.hpp> -// Temporary / will be Strategy as template parameter -#include <boost/geometry/strategies/side.hpp> +#include <boost/geometry/strategies/agnostic/point_in_poly_winding.hpp> +#include <boost/geometry/strategies/cartesian/area_surveyor.hpp> +#include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp> #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp> - -#include <boost/geometry/strategies/side_info.hpp> +#include <boost/geometry/strategies/covered_by.hpp> #include <boost/geometry/strategies/intersection.hpp> #include <boost/geometry/strategies/intersection_result.hpp> +#include <boost/geometry/strategies/side.hpp> +#include <boost/geometry/strategies/side_info.hpp> +#include <boost/geometry/strategies/within.hpp> #include <boost/geometry/policies/robustness/robust_point_type.hpp> #include <boost/geometry/policies/robustness/segment_ratio_type.hpp> @@ -61,10 +64,74 @@ namespace strategy { namespace intersection /*! \see http://mathworld.wolfram.com/Line-LineIntersection.html */ -template <typename Policy, typename CalculationType = void> -struct relate_cartesian_segments +template +< + typename CalculationType = void +> +struct cartesian_segments { - typedef typename Policy::return_type return_type; + typedef side::side_by_triangle<CalculationType> side_strategy_type; + + static inline side_strategy_type get_side_strategy() + { + return side_strategy_type(); + } + + 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> + static inline typename point_in_geometry_strategy<Geometry1, Geometry2>::type + get_point_in_geometry_strategy() + { + typedef typename point_in_geometry_strategy + < + Geometry1, Geometry2 + >::type strategy_type; + return strategy_type(); + } + + template <typename Geometry> + struct area_strategy + { + typedef area::surveyor + < + typename point_type<Geometry>::type, + CalculationType + > type; + }; + + template <typename Geometry> + static inline typename area_strategy<Geometry>::type get_area_strategy() + { + typedef typename area_strategy<Geometry>::type strategy_type; + return strategy_type(); + } + + template <typename Geometry> + struct distance_strategy + { + typedef distance::pythagoras + < + CalculationType + > type; + }; + + template <typename Geometry> + static inline typename distance_strategy<Geometry>::type get_distance_strategy() + { + typedef typename distance_strategy<Geometry>::type strategy_type; + return strategy_type(); + } template <typename CoordinateType, typename SegmentRatio> struct segment_intersection_info @@ -144,9 +211,16 @@ struct relate_cartesian_segments // Relate segments a and b - template <typename Segment1, typename Segment2, typename RobustPolicy> - static inline return_type apply(Segment1 const& a, Segment2 const& b, - RobustPolicy const& robust_policy) + template + < + typename Segment1, + typename Segment2, + typename Policy, + typename RobustPolicy + > + static inline typename Policy::return_type + apply(Segment1 const& a, Segment2 const& b, + Policy const& policy, RobustPolicy const& robust_policy) { // type them all as in Segment1 - TODO reconsider this, most precise? typedef typename geometry::point_type<Segment1>::type point_type; @@ -169,18 +243,25 @@ struct relate_cartesian_segments geometry::recalculate(b0_rob, b0, robust_policy); geometry::recalculate(b1_rob, b1, robust_policy); - return apply(a, b, robust_policy, a0_rob, a1_rob, b0_rob, b1_rob); + return apply(a, b, policy, robust_policy, a0_rob, a1_rob, b0_rob, b1_rob); } // The main entry-routine, calculating intersections of segments a / b // NOTE: Robust* types may be the same as Segments' point types - template <typename Segment1, typename Segment2, - typename RobustPolicy, - typename RobustPoint1, typename RobustPoint2> - static inline return_type apply(Segment1 const& a, Segment2 const& b, - RobustPolicy const& /*robust_policy*/, - RobustPoint1 const& robust_a1, RobustPoint1 const& robust_a2, - RobustPoint2 const& robust_b1, RobustPoint2 const& robust_b2) + template + < + typename Segment1, + typename Segment2, + typename Policy, + typename RobustPolicy, + typename RobustPoint1, + typename RobustPoint2 + > + static inline typename Policy::return_type + apply(Segment1 const& a, Segment2 const& b, + Policy const&, RobustPolicy const& /*robust_policy*/, + RobustPoint1 const& robust_a1, RobustPoint1 const& robust_a2, + RobustPoint2 const& robust_b1, RobustPoint2 const& robust_b2) { BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) ); BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) ); @@ -197,14 +278,9 @@ struct relate_cartesian_segments ; } - typedef typename select_calculation_type - <Segment1, Segment2, CalculationType>::type coordinate_type; - - typedef side::side_by_triangle<coordinate_type> side; - side_info sides; - sides.set<0>(side::apply(robust_b1, robust_b2, robust_a1), - side::apply(robust_b1, robust_b2, robust_a2)); + sides.set<0>(side_strategy_type::apply(robust_b1, robust_b2, robust_a1), + side_strategy_type::apply(robust_b1, robust_b2, robust_a2)); if (sides.same<0>()) { @@ -212,8 +288,8 @@ struct relate_cartesian_segments return Policy::disjoint(); } - sides.set<1>(side::apply(robust_a1, robust_a2, robust_b1), - side::apply(robust_a1, robust_a2, robust_b2)); + sides.set<1>(side_strategy_type::apply(robust_a1, robust_a2, robust_b1), + side_strategy_type::apply(robust_a1, robust_a2, robust_b2)); if (sides.same<1>()) { @@ -230,16 +306,16 @@ struct relate_cartesian_segments >::type robust_coordinate_type; typedef typename segment_ratio_type - < - typename geometry::point_type<Segment1>::type, // TODO: most precise point? - RobustPolicy - >::type ratio_type; + < + typename geometry::point_type<Segment1>::type, // TODO: most precise point? + RobustPolicy + >::type ratio_type; segment_intersection_info - < - coordinate_type, - ratio_type - > sinfo; + < + typename select_calculation_type<Segment1, Segment2, CalculationType>::type, + ratio_type + > sinfo; sinfo.dx_a = get<1, 0>(a) - get<0, 0>(a); // distance in x-dir sinfo.dx_b = get<1, 0>(b) - get<0, 0>(b); @@ -303,14 +379,14 @@ struct relate_cartesian_segments if (collinear_use_first.first) { - return relate_collinear<0, ratio_type>(a, b, + return relate_collinear<0, Policy, ratio_type>(a, b, robust_a1, robust_a2, robust_b1, robust_b2, a_is_point, b_is_point); } else { // Y direction contains larger segments (maybe dx is zero) - return relate_collinear<1, ratio_type>(a, b, + return relate_collinear<1, Policy, ratio_type>(a, b, robust_a1, robust_a2, robust_b1, robust_b2, a_is_point, b_is_point); } @@ -358,33 +434,35 @@ private: template < std::size_t Dimension, + typename Policy, typename RatioType, typename Segment1, typename Segment2, typename RobustPoint1, typename RobustPoint2 > - static inline return_type relate_collinear(Segment1 const& a, - Segment2 const& b, - RobustPoint1 const& robust_a1, RobustPoint1 const& robust_a2, - RobustPoint2 const& robust_b1, RobustPoint2 const& robust_b2, - bool a_is_point, bool b_is_point) + static inline typename Policy::return_type + relate_collinear(Segment1 const& a, + Segment2 const& b, + RobustPoint1 const& robust_a1, RobustPoint1 const& robust_a2, + RobustPoint2 const& robust_b1, RobustPoint2 const& robust_b2, + bool a_is_point, bool b_is_point) { if (a_is_point) { - return relate_one_degenerate<RatioType>(a, + return relate_one_degenerate<Policy, RatioType>(a, get<Dimension>(robust_a1), get<Dimension>(robust_b1), get<Dimension>(robust_b2), true); } if (b_is_point) { - return relate_one_degenerate<RatioType>(b, + return relate_one_degenerate<Policy, RatioType>(b, get<Dimension>(robust_b1), get<Dimension>(robust_a1), get<Dimension>(robust_a2), false); } - return relate_collinear<RatioType>(a, b, + return relate_collinear<Policy, RatioType>(a, b, get<Dimension>(robust_a1), get<Dimension>(robust_a2), get<Dimension>(robust_b1), @@ -394,17 +472,17 @@ private: /// Relate segments known collinear template < + typename Policy, typename RatioType, typename Segment1, typename Segment2, typename RobustType1, typename RobustType2 > - static inline return_type relate_collinear(Segment1 const& a - , Segment2 const& b - , RobustType1 oa_1, RobustType1 oa_2 - , RobustType2 ob_1, RobustType2 ob_2 - ) + static inline typename Policy::return_type + relate_collinear(Segment1 const& a, Segment2 const& b, + RobustType1 oa_1, RobustType1 oa_2, + RobustType2 ob_1, RobustType2 ob_2) { // Calculate the ratios where a starts in b, b starts in a // a1--------->a2 (2..7) @@ -496,17 +574,16 @@ private: /// Relate segments where one is degenerate template < + typename Policy, typename RatioType, typename DegenerateSegment, typename RobustType1, typename RobustType2 > - static inline return_type relate_one_degenerate( - DegenerateSegment const& degenerate_segment - , RobustType1 d - , RobustType2 s1, RobustType2 s2 - , bool a_degenerate - ) + static inline typename Policy::return_type + relate_one_degenerate(DegenerateSegment const& degenerate_segment, + RobustType1 d, RobustType2 s1, RobustType2 s2, + bool a_degenerate) { // Calculate the ratios where ds starts in s // a1--------->a2 (2..6) @@ -546,10 +623,10 @@ private: namespace services { -template <typename Policy, typename CalculationType> -struct default_strategy<cartesian_tag, Policy, CalculationType> +template <typename CalculationType> +struct default_strategy<cartesian_tag, CalculationType> { - typedef relate_cartesian_segments<Policy, CalculationType> type; + typedef cartesian_segments<CalculationType> type; }; } // namespace services @@ -558,6 +635,69 @@ struct default_strategy<cartesian_tag, Policy, CalculationType> }} // namespace strategy::intersection +namespace strategy +{ + +namespace within { namespace services +{ + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, cartesian_tag, cartesian_tag> +{ + typedef strategy::intersection::cartesian_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, cartesian_tag, cartesian_tag> +{ + typedef strategy::intersection::cartesian_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, cartesian_tag, cartesian_tag> +{ + typedef strategy::intersection::cartesian_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, cartesian_tag, cartesian_tag> +{ + typedef strategy::intersection::cartesian_segments<> type; +}; + +}} // within::services + +namespace covered_by { namespace services +{ + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, cartesian_tag, cartesian_tag> +{ + typedef strategy::intersection::cartesian_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, cartesian_tag, cartesian_tag> +{ + typedef strategy::intersection::cartesian_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, cartesian_tag, cartesian_tag> +{ + typedef strategy::intersection::cartesian_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, cartesian_tag, cartesian_tag> +{ + typedef strategy::intersection::cartesian_segments<> type; +}; + +}} // within::services + +} // strategy + }} // namespace boost::geometry diff --git a/boost/geometry/strategies/cartesian/point_in_box.hpp b/boost/geometry/strategies/cartesian/point_in_box.hpp index 227a98f2ad..1c14125a6d 100644 --- a/boost/geometry/strategies/cartesian/point_in_box.hpp +++ b/boost/geometry/strategies/cartesian/point_in_box.hpp @@ -4,8 +4,8 @@ // Copyright (c) 2008-2015 Bruno Lalande, Paris, France. // Copyright (c) 2009-2015 Mateusz Loskot, London, UK. -// 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 @@ -219,10 +219,10 @@ namespace within { namespace services template <typename Point, typename Box> struct default_strategy < + Point, Box, point_tag, box_tag, - point_tag, areal_tag, - cartesian_tag, cartesian_tag, - Point, Box + pointlike_tag, areal_tag, + cartesian_tag, cartesian_tag > { typedef within::point_in_box<Point, Box> type; @@ -232,10 +232,10 @@ struct default_strategy template <typename Point, typename Box> struct default_strategy < + Point, Box, point_tag, box_tag, - point_tag, areal_tag, - spherical_tag, spherical_tag, - Point, Box + pointlike_tag, areal_tag, + spherical_tag, spherical_tag > { typedef within::point_in_box<Point, Box> type; @@ -252,10 +252,10 @@ namespace covered_by { namespace services template <typename Point, typename Box> struct default_strategy < + Point, Box, point_tag, box_tag, - point_tag, areal_tag, - cartesian_tag, cartesian_tag, - Point, Box + pointlike_tag, areal_tag, + cartesian_tag, cartesian_tag > { typedef within::point_in_box @@ -269,10 +269,10 @@ struct default_strategy template <typename Point, typename Box> struct default_strategy < + Point, Box, point_tag, box_tag, - point_tag, areal_tag, - spherical_tag, spherical_tag, - Point, Box + pointlike_tag, areal_tag, + spherical_tag, spherical_tag > { typedef within::point_in_box diff --git a/boost/geometry/strategies/covered_by.hpp b/boost/geometry/strategies/covered_by.hpp index a5aae7703b..363e34b68f 100644 --- a/boost/geometry/strategies/covered_by.hpp +++ b/boost/geometry/strategies/covered_by.hpp @@ -4,6 +4,11 @@ // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. // Copyright (c) 2009-2012 Mateusz Loskot, London, UK. +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 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. @@ -16,6 +21,12 @@ #include <boost/mpl/assert.hpp> +#include <boost/geometry/core/cs.hpp> +#include <boost/geometry/core/point_type.hpp> +#include <boost/geometry/core/tag.hpp> +#include <boost/geometry/core/tags.hpp> +#include <boost/geometry/core/tag_cast.hpp> + namespace boost { namespace geometry { @@ -31,23 +42,39 @@ namespace services /*! \brief Traits class binding a covered_by determination strategy to a coordinate system \ingroup covered_by -\tparam TagContained tag (possibly casted) of point-type -\tparam TagContained tag (possibly casted) of (possibly) containing type -\tparam CsTagContained tag of coordinate system of point-type -\tparam CsTagContaining tag of coordinate system of (possibly) containing type -\tparam Geometry geometry-type of input (often point, or box) +\tparam GeometryContained geometry-type of input (possibly) contained type \tparam GeometryContaining geometry-type of input (possibly) containing type +\tparam TagContained casted tag of (possibly) contained type +\tparam TagContaining casted tag of (possibly) containing type +\tparam CsTagContained tag of coordinate system of (possibly) contained type +\tparam CsTagContaining tag of coordinate system of (possibly) containing type */ template < - typename TagContained, - typename TagContaining, - typename CastedTagContained, - typename CastedTagContaining, - typename CsTagContained, - typename CsTagContaining, typename GeometryContained, - typename GeometryContaining + typename GeometryContaining, + typename TagContained = typename tag<GeometryContained>::type, + typename TagContaining = typename tag<GeometryContaining>::type, + typename CastedTagContained = typename tag_cast + < + typename tag<GeometryContained>::type, + pointlike_tag, linear_tag, polygonal_tag, areal_tag + >::type, + typename CastedTagContaining = typename tag_cast + < + typename tag<GeometryContaining>::type, + pointlike_tag, linear_tag, polygonal_tag, areal_tag + >::type, + typename CsTagContained = typename tag_cast + < + typename cs_tag<typename point_type<GeometryContained>::type>::type, + spherical_tag + >::type, + typename CsTagContaining = typename tag_cast + < + typename cs_tag<typename point_type<GeometryContaining>::type>::type, + spherical_tag + >::type > struct default_strategy { diff --git a/boost/geometry/strategies/disjoint.hpp b/boost/geometry/strategies/disjoint.hpp new file mode 100644 index 0000000000..f844cd86ad --- /dev/null +++ b/boost/geometry/strategies/disjoint.hpp @@ -0,0 +1,90 @@ +// 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_DISJOINT_HPP +#define BOOST_GEOMETRY_STRATEGIES_DISJOINT_HPP + + +#include <boost/mpl/assert.hpp> +#include <boost/type_traits/is_same.hpp> + +#include <boost/geometry/core/cs.hpp> +#include <boost/geometry/core/point_type.hpp> +#include <boost/geometry/core/topological_dimension.hpp> + +#include <boost/geometry/strategies/covered_by.hpp> +#include <boost/geometry/strategies/default_strategy.hpp> +#include <boost/geometry/strategies/relate.hpp> + + +namespace boost { namespace geometry { namespace strategy { namespace disjoint +{ + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +template +< + typename Geometry1, + typename Geometry2, + typename Tag1 = typename geometry::tag<Geometry1>::type, + typename Tag2 = typename geometry::tag<Geometry2>::type, + int TopDim1 = geometry::topological_dimension<Geometry1>::value, + int TopDim2 = geometry::topological_dimension<Geometry2>::value +> +struct default_strategy + : relate::services::default_strategy + < + Geometry1, Geometry2 + > +{}; + +template <typename Point, typename Box> +struct default_strategy<Point, Box, point_tag, box_tag, 0, 2> + : strategy::covered_by::services::default_strategy<Point, Box> +{}; + +template <typename Box, typename Point> +struct default_strategy<Box, Point, box_tag, point_tag, 2, 0> + : strategy::covered_by::services::default_strategy<Point, Box> +{}; + +template <typename MultiPoint, typename Box> +struct default_strategy<MultiPoint, Box, multi_point_tag, box_tag, 0, 2> + : strategy::covered_by::services::default_strategy + < + typename point_type<MultiPoint>::type, + Box + > +{}; + +template <typename Box, typename MultiPoint> +struct default_strategy<Box, MultiPoint, box_tag, multi_point_tag, 2, 0> + : strategy::covered_by::services::default_strategy + < + typename point_type<MultiPoint>::type, + Box + > +{}; + +template <typename Box1, typename Box2> +struct default_strategy<Box1, Box2, box_tag, box_tag, 2, 2> +{ + // dummy strategy which will be ignored + typedef geometry::default_strategy type; +}; + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +}}}} // namespace boost::geometry::strategy::disjoint + + +#endif // BOOST_GEOMETRY_STRATEGIES_DISJOINT_HPP diff --git a/boost/geometry/strategies/envelope.hpp b/boost/geometry/strategies/envelope.hpp new file mode 100644 index 0000000000..fde9c858a6 --- /dev/null +++ b/boost/geometry/strategies/envelope.hpp @@ -0,0 +1,45 @@ +// 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_ENVELOPE_HPP +#define BOOST_GEOMETRY_STRATEGIES_ENVELOPE_HPP + +#include <boost/mpl/assert.hpp> + +namespace boost { namespace geometry +{ + + +namespace strategy { namespace envelope { namespace services +{ + +/*! +\brief Traits class binding a default envelope strategy to a coordinate system +\ingroup util +\tparam CSTag tag of coordinate system +\tparam CalculationType \tparam_calculation +*/ +template <typename CSTag, typename CalculationType = void> +struct default_strategy +{ + BOOST_MPL_ASSERT_MSG + ( + false, NOT_IMPLEMENTED_FOR_THIS_TYPE + , (types<CSTag>) + ); +}; + +}}} // namespace strategy::envelope::services + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_ENVELOPE_HPP + 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) {} }; diff --git a/boost/geometry/strategies/intersection.hpp b/boost/geometry/strategies/intersection.hpp index f51c5cb206..e5662c9e46 100644 --- a/boost/geometry/strategies/intersection.hpp +++ b/boost/geometry/strategies/intersection.hpp @@ -1,6 +1,6 @@ // Boost.Geometry -// Copyright (c) 2016, Oracle and/or its affiliates. +// 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, @@ -30,10 +30,9 @@ namespace services \brief Traits class binding a segments intersection strategy to a coordinate system \ingroup util \tparam CSTag tag of coordinate system of point-type -\tparam Policy intersection policy \tparam CalculationType \tparam_calculation */ -template <typename CSTag, typename Policy, typename CalculationType = void> +template <typename CSTag, typename CalculationType = void> struct default_strategy { BOOST_MPL_ASSERT_MSG diff --git a/boost/geometry/strategies/intersection_strategies.hpp b/boost/geometry/strategies/intersection_strategies.hpp index 0452c4692c..a173505804 100644 --- a/boost/geometry/strategies/intersection_strategies.hpp +++ b/boost/geometry/strategies/intersection_strategies.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2016. -// Modifications copyright (c) 2016, Oracle and/or its affiliates. +// 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 // Use, modification and distribution is subject to the Boost Software License, @@ -22,10 +22,10 @@ #include <boost/geometry/policies/relate/tupled.hpp> #include <boost/geometry/strategies/intersection.hpp> -#include <boost/geometry/strategies/side.hpp> #include <boost/geometry/strategies/intersection_result.hpp> +#include <boost/geometry/strategies/side.hpp> -#include <boost/geometry/strategies/cartesian/cart_intersect.hpp> +#include <boost/geometry/strategies/cartesian/intersection.hpp> #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp> #include <boost/geometry/strategies/spherical/intersection.hpp> #include <boost/geometry/strategies/spherical/ssf.hpp> @@ -70,17 +70,18 @@ private : > ip_type; public: + typedef policies::relate::segments_tupled + < + policies::relate::segments_intersection_points + < + ip_type + > , + policies::relate::segments_direction + > intersection_policy_type; + typedef typename strategy::intersection::services::default_strategy < Tag, - policies::relate::segments_tupled - < - policies::relate::segments_intersection_points - < - ip_type - > , - policies::relate::segments_direction - >, CalculationType >::type segment_intersection_strategy_type; diff --git a/boost/geometry/strategies/relate.hpp b/boost/geometry/strategies/relate.hpp new file mode 100644 index 0000000000..ffeed78959 --- /dev/null +++ b/boost/geometry/strategies/relate.hpp @@ -0,0 +1,177 @@ +// 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_RELATE_HPP +#define BOOST_GEOMETRY_STRATEGIES_RELATE_HPP + + +#include <boost/mpl/assert.hpp> +#include <boost/type_traits/is_same.hpp> + +#include <boost/geometry/core/cs.hpp> +#include <boost/geometry/core/point_type.hpp> +#include <boost/geometry/core/topological_dimension.hpp> + +#include <boost/geometry/strategies/covered_by.hpp> +#include <boost/geometry/strategies/intersection.hpp> +#include <boost/geometry/strategies/within.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy +{ + +namespace point_in_geometry +{ + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +template +< + typename Point, + typename Geometry, + typename Tag1 = typename tag<Point>::type, + typename Tag2 = typename tag<Geometry>::type +> +struct default_strategy + : strategy::within::services::default_strategy + < + Point, + Geometry + > +{ + typedef typename default_strategy::type within_strategy_type; + + typedef typename strategy::covered_by::services::default_strategy + < + Point, + Geometry + >::type covered_by_strategy_type; + + static const bool same_strategies = boost::is_same<within_strategy_type, covered_by_strategy_type>::value; + BOOST_MPL_ASSERT_MSG((same_strategies), + DEFAULT_WITHIN_AND_COVERED_BY_STRATEGIES_NOT_COMPATIBLE, + (within_strategy_type, covered_by_strategy_type)); +}; + +template<typename Point, typename Geometry> +struct default_strategy<Point, Geometry, point_tag, point_tag> + : strategy::within::services::default_strategy<Point, Geometry> +{}; + +template<typename Point, typename Geometry> +struct default_strategy<Point, Geometry, point_tag, multi_point_tag> + : strategy::within::services::default_strategy<Point, Geometry> +{}; + + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +} // namespace point_in_geometry + +namespace relate +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +template <typename Geometry> +struct default_intersection_strategy + : strategy::intersection::services::default_strategy + < + typename cs_tag<Geometry>::type + > +{}; + +template <typename PointLike, typename Geometry> +struct default_point_in_geometry_strategy + : point_in_geometry::services::default_strategy + < + typename point_type<PointLike>::type, + Geometry + > +{}; + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +template +< + typename Geometry1, + typename Geometry2, + int TopDim1 = geometry::topological_dimension<Geometry1>::value, + int TopDim2 = geometry::topological_dimension<Geometry2>::value +> +struct default_strategy +{ + BOOST_MPL_ASSERT_MSG + ( + false, NOT_IMPLEMENTED_FOR_THESE_TYPES + , (types<Geometry1, Geometry2>) + ); +}; + +template <typename PointLike1, typename PointLike2> +struct default_strategy<PointLike1, PointLike2, 0, 0> + : detail::default_point_in_geometry_strategy<PointLike1, PointLike2> +{}; + +template <typename PointLike, typename Geometry, int TopDim2> +struct default_strategy<PointLike, Geometry, 0, TopDim2> + : detail::default_point_in_geometry_strategy<PointLike, Geometry> +{}; + +template <typename Geometry, typename PointLike, int TopDim1> +struct default_strategy<Geometry, PointLike, TopDim1, 0> + : detail::default_point_in_geometry_strategy<PointLike, Geometry> +{}; + +template <typename Geometry1, typename Geometry2> +struct default_strategy<Geometry1, Geometry2, 1, 1> + : detail::default_intersection_strategy<Geometry1> +{}; + +template <typename Geometry1, typename Geometry2> +struct default_strategy<Geometry1, Geometry2, 1, 2> + : detail::default_intersection_strategy<Geometry1> +{}; + +template <typename Geometry1, typename Geometry2> +struct default_strategy<Geometry1, Geometry2, 2, 1> + : detail::default_intersection_strategy<Geometry1> +{}; + +template <typename Geometry1, typename Geometry2> +struct default_strategy<Geometry1, Geometry2, 2, 2> + : detail::default_intersection_strategy<Geometry1> +{}; + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +} // namespace relate + +} // namespace strategy + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_RELATE_HPP diff --git a/boost/geometry/strategies/side.hpp b/boost/geometry/strategies/side.hpp index 376f2fdf1b..9aaa2bdddc 100644 --- a/boost/geometry/strategies/side.hpp +++ b/boost/geometry/strategies/side.hpp @@ -30,16 +30,16 @@ namespace services /*! \brief Traits class binding a side determination strategy to a coordinate system \ingroup util -\tparam Tag tag of coordinate system of point-type +\tparam CSTag tag of coordinate system of point-type \tparam CalculationType \tparam_calculation */ -template <typename Tag, typename CalculationType = void> +template <typename CSTag, typename CalculationType = void> struct default_strategy { BOOST_MPL_ASSERT_MSG ( false, NOT_IMPLEMENTED_FOR_THIS_TYPE - , (types<Tag>) + , (types<CSTag>) ); }; diff --git a/boost/geometry/strategies/spherical/area.hpp b/boost/geometry/strategies/spherical/area.hpp new file mode 100644 index 0000000000..206b734548 --- /dev/null +++ b/boost/geometry/strategies/spherical/area.hpp @@ -0,0 +1,182 @@ +// 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_SPHERICAL_AREA_HPP +#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HPP + + +#include <boost/geometry/formulas/area_formulas.hpp> +#include <boost/geometry/core/radius.hpp> +#include <boost/geometry/core/srs.hpp> +#include <boost/geometry/strategies/area.hpp> + + +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 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 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 + < + typename coordinate_type<PointOfSegment>::type, + double + >::type, + CalculationType + >::type CT; + +protected : + struct excess_sum + { + CT m_sum; + + // Keep track if encircles some pole + size_t m_crosses_prime_meridian; + + inline excess_sum() + : m_sum(0) + , m_crosses_prime_meridian(0) + {} + template <typename SphereType> + inline CT area(SphereType sphere) const + { + CT result; + CT radius = geometry::get_radius<0>(sphere); + + // Encircles pole + if(m_crosses_prime_meridian % 2 == 1) + { + size_t times_crosses_prime_meridian + = 1 + (m_crosses_prime_meridian / 2); + + result = CT(2) + * geometry::math::pi<CT>() + * times_crosses_prime_meridian + - geometry::math::abs(m_sum); + + if(geometry::math::sign<CT>(m_sum) == 1) + { + result = - result; + } + + } else { + result = m_sum; + } + + result *= radius * radius; + + return result; + } + }; + +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)) + {} + + explicit inline spherical(CT const& radius) + : m_sphere(radius) + {} + + inline void apply(PointOfSegment const& p1, + PointOfSegment const& p2, + excess_sum& state) const + { + if (! geometry::math::equals(get<0>(p1), get<0>(p2))) + { + + state.m_sum += geometry::formula::area_formulas + <CT>::template spherical<LongSegment>(p1, p2); + + // Keep track whenever a segment crosses the prime meridian + geometry::formula::area_formulas + <CT>::crosses_prime_meridian(p1, p2, state); + + } + } + + inline return_type result(excess_sum const& state) const + { + return state.area(m_sphere); + } + +private : + /// srs Sphere + sphere_type m_sphere; +}; + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +namespace services +{ + + +template <typename Point> +struct default_strategy<spherical_equatorial_tag, Point> +{ + typedef strategy::area::spherical<Point> type; +}; + +// Note: spherical polar coordinate system requires "get_as_radian_equatorial" +template <typename Point> +struct default_strategy<spherical_polar_tag, Point> +{ + typedef strategy::area::spherical<Point> type; +}; + +} // namespace services + +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::area + + + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HPP diff --git a/boost/geometry/strategies/spherical/area_huiller.hpp b/boost/geometry/strategies/spherical/area_huiller.hpp deleted file mode 100644 index 37d8d20124..0000000000 --- a/boost/geometry/strategies/spherical/area_huiller.hpp +++ /dev/null @@ -1,214 +0,0 @@ -// Boost.Geometry (aka GGL, Generic Geometry Library) - -// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. - -// This file was modified by Oracle on 2015. -// Modifications copyright (c) 2015, Oracle and/or its affiliates. - -// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle -// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle - -// Use, modification and distribution is subject to the Boost Software License, -// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) - -#ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP -#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP - - - -#include <boost/geometry/strategies/spherical/distance_haversine.hpp> - -#include <boost/geometry/core/radian_access.hpp> -#include <boost/geometry/util/math.hpp> - - -namespace boost { namespace geometry -{ - -namespace strategy { namespace area -{ - - - -/*! -\brief Area calculation by spherical excess / Huiller's formula -\ingroup strategies -\tparam PointOfSegment point type of segments of rings/polygons -\tparam CalculationType \tparam_calculation -\author Barend Gehrels. Adapted from: -- http://webdocs.cs.ualberta.ca/~graphics/books/GraphicsGems/gemsiv/sph_poly.c -- http://tog.acm.org/resources/GraphicsGems/gemsiv/sph_poly.c -- http://williams.best.vwh.net/avform.htm -\note The version in Graphics Gems IV (page 132-137) didn't account for -polygons crossing the 0 and 180 meridians. The fix for this algorithm -can be found in Graphics Gems V (pages 45-46). See: -- http://kysmykseka.net/koti/wizardry/Game%20Development/Programming/Graphics%20Gems%204.pdf -- http://kysmykseka.net/koti/wizardry/Game%20Development/Programming/Graphics%20Gems%205.pdf -\note This version works for convex and non-convex polygons, for 180 meridian -crossing polygons and for polygons with holes. However, some cases (especially -180 meridian cases) must still be checked. -\note The version which sums angles, which is often seen, doesn't handle non-convex -polygons correctly. -\note The version which sums longitudes, see http://hdl.handle.net/2014/40409, -is simple and works well in most cases but not in 180 meridian crossing cases. -This probably could be solved. - -\note This version is made for spherical equatorial coordinate systems - -\qbk{ - -[heading Example] -[area_with_strategy] -[area_with_strategy_output] - - -[heading See also] -[link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)] -} - -*/ -template -< - typename PointOfSegment, - typename CalculationType = void -> -class huiller -{ -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 calculation_type; - -protected : - struct excess_sum - { - calculation_type sum; - - // Distances are calculated on unit sphere here - strategy::distance::haversine<calculation_type> distance_over_unit_sphere; - - - inline excess_sum() - : sum(0) - , distance_over_unit_sphere(1) - {} - inline calculation_type area(calculation_type radius) const - { - return - sum * radius * radius; - } - }; - -public : - typedef calculation_type return_type; - typedef PointOfSegment segment_point_type; - typedef excess_sum state_type; - - inline huiller(calculation_type radius = 1.0) - : m_radius(radius) - {} - - inline void apply(PointOfSegment const& p1, - PointOfSegment const& p2, - excess_sum& state) const - { - if (! geometry::math::equals(get<0>(p1), get<0>(p2))) - { - calculation_type const half = 0.5; - calculation_type const two = 2.0; - calculation_type const four = 4.0; - calculation_type const pi - = geometry::math::pi<calculation_type>(); - calculation_type const two_pi - = geometry::math::two_pi<calculation_type>(); - calculation_type const half_pi - = geometry::math::half_pi<calculation_type>(); - - // Distance p1 p2 - calculation_type a = state.distance_over_unit_sphere.apply(p1, p2); - - // Sides on unit sphere to south pole - calculation_type b = half_pi - geometry::get_as_radian<1>(p2); - calculation_type c = half_pi - geometry::get_as_radian<1>(p1); - - // Semi parameter - calculation_type s = half * (a + b + c); - - // E: spherical excess, using l'Huiller's formula - // [tg(e / 4)]2 = tg[s / 2] tg[(s-a) / 2] tg[(s-b) / 2] tg[(s-c) / 2] - calculation_type excess = four - * atan(geometry::math::sqrt(geometry::math::abs(tan(s / two) - * tan((s - a) / two) - * tan((s - b) / two) - * tan((s - c) / two)))); - - excess = geometry::math::abs(excess); - - // In right direction: positive, add area. In left direction: negative, subtract area. - // Longitude comparisons are not so obvious. If one is negative and other is positive, - // we have to take the dateline into account. - - calculation_type lon_diff = geometry::get_as_radian<0>(p2) - - geometry::get_as_radian<0>(p1); - if (lon_diff <= 0) - { - lon_diff += two_pi; - } - - if (lon_diff > pi) - { - excess = -excess; - } - - state.sum += excess; - } - } - - inline return_type result(excess_sum const& state) const - { - return state.area(m_radius); - } - -private : - /// Radius of the sphere - calculation_type m_radius; -}; - -#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - -namespace services -{ - - -template <typename Point> -struct default_strategy<spherical_equatorial_tag, Point> -{ - typedef strategy::area::huiller<Point> type; -}; - -// Note: spherical polar coordinate system requires "get_as_radian_equatorial" -/***template <typename Point> -struct default_strategy<spherical_polar_tag, Point> -{ - typedef strategy::area::huiller<Point> type; -};***/ - -} // namespace services - -#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS - - -}} // namespace strategy::area - - - - -}} // namespace boost::geometry - -#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HUILLER_HPP diff --git a/boost/geometry/strategies/spherical/azimuth.hpp b/boost/geometry/strategies/spherical/azimuth.hpp new file mode 100644 index 0000000000..3c208fe2e2 --- /dev/null +++ b/boost/geometry/strategies/spherical/azimuth.hpp @@ -0,0 +1,87 @@ +// 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_SPHERICAL_AZIMUTH_HPP +#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AZIMUTH_HPP + + +#include <boost/geometry/strategies/azimuth.hpp> +#include <boost/geometry/formulas/spherical.hpp> + +#include <boost/mpl/if.hpp> +#include <boost/type_traits/is_void.hpp> + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace azimuth +{ + +template +< + typename CalculationType = void +> +class spherical +{ +public : + + inline spherical() + {} + + template <typename T> + static inline void apply(T const& lon1_rad, T const& lat1_rad, + T const& lon2_rad, T const& lat2_rad, + T& a1, T& a2) + { + typedef typename boost::mpl::if_ + < + boost::is_void<CalculationType>, T, CalculationType + >::type calc_t; + + geometry::formula::result_spherical<calc_t> + result = geometry::formula::spherical_azimuth<calc_t, true>( + calc_t(lon1_rad), calc_t(lat1_rad), + calc_t(lon2_rad), calc_t(lat2_rad)); + + a1 = result.azimuth; + a2 = result.reverse_azimuth; + } + +}; + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +namespace services +{ + +template <typename CalculationType> +struct default_strategy<spherical_equatorial_tag, CalculationType> +{ + typedef strategy::azimuth::spherical<CalculationType> type; +}; + +/* +template <typename CalculationType> +struct default_strategy<spherical_polar_tag, CalculationType> +{ + typedef strategy::azimuth::spherical<CalculationType> type; +}; +*/ +} + +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +}} // namespace strategy::azimuth + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AZIMUTH_HPP diff --git a/boost/geometry/strategies/spherical/envelope_segment.hpp b/boost/geometry/strategies/spherical/envelope_segment.hpp new file mode 100644 index 0000000000..98f085fe73 --- /dev/null +++ b/boost/geometry/strategies/spherical/envelope_segment.hpp @@ -0,0 +1,86 @@ +// 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_SPHERICAL_ENVELOPE_SEGMENT_HPP +#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_ENVELOPE_SEGMENT_HPP + +#include <boost/geometry/algorithms/detail/envelope/segment.hpp> +#include <boost/geometry/algorithms/detail/normalize.hpp> +#include <boost/geometry/strategies/envelope.hpp> +#include <boost/geometry/strategies/spherical/azimuth.hpp> + +namespace boost { namespace geometry +{ + +namespace strategy { namespace envelope +{ + +template +< + typename CalculationType = void +> +class spherical_segment +{ +public : + + inline spherical_segment() + {} + + 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::spherical<CalculationType> azimuth_spherical; + + typedef typename coordinate_system<Point1>::type::units units_type; + + geometry::detail::envelope::envelope_segment_impl<spherical_equatorial_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_spherical); + + } +}; + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +namespace services +{ + +template <typename CalculationType> +struct default_strategy<spherical_equatorial_tag, CalculationType> +{ + typedef strategy::envelope::spherical_segment<CalculationType> type; +}; + + +template <typename CalculationType> +struct default_strategy<spherical_polar_tag, CalculationType> +{ + typedef strategy::envelope::spherical_segment<CalculationType> type; +}; + +} + +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + + +}} // namespace strategy::envelope + +}} //namepsace boost::geometry + +#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_ENVELOPE_SEGMENT_HPP + diff --git a/boost/geometry/strategies/spherical/intersection.hpp b/boost/geometry/strategies/spherical/intersection.hpp index 4ffc853aad..5d37583333 100644 --- a/boost/geometry/strategies/spherical/intersection.hpp +++ b/boost/geometry/strategies/spherical/intersection.hpp @@ -1,6 +1,6 @@ // Boost.Geometry -// Copyright (c) 2016, Oracle and/or its affiliates. +// 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, @@ -25,6 +25,7 @@ #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/formulas/spherical.hpp> #include <boost/geometry/geometries/concepts/point_concept.hpp> @@ -32,9 +33,16 @@ #include <boost/geometry/policies/robustness/segment_ratio.hpp> -#include <boost/geometry/strategies/side_info.hpp> +#include <boost/geometry/strategies/agnostic/point_in_poly_winding.hpp> +#include <boost/geometry/strategies/covered_by.hpp> #include <boost/geometry/strategies/intersection.hpp> #include <boost/geometry/strategies/intersection_result.hpp> +#include <boost/geometry/strategies/side.hpp> +#include <boost/geometry/strategies/side_info.hpp> +#include <boost/geometry/strategies/spherical/area.hpp> +#include <boost/geometry/strategies/spherical/distance_haversine.hpp> +#include <boost/geometry/strategies/spherical/ssf.hpp> +#include <boost/geometry/strategies/within.hpp> #include <boost/geometry/util/math.hpp> #include <boost/geometry/util/select_calculation_type.hpp> @@ -68,13 +76,80 @@ namespace strategy { namespace intersection // For now, intersection points near the endpoints are checked explicitly if needed (if the IP is near the endpoint) // to generate precise result for them. Only the crossing (i) case may suffer from lower precision. -template <typename Policy, typename CalculationType = void> -struct relate_spherical_segments +template +< + typename CalcPolicy, + typename CalculationType = void +> +struct ecef_segments { - typedef typename Policy::return_type return_type; + typedef side::spherical_side_formula<CalculationType> side_strategy_type; + + static inline side_strategy_type get_side_strategy() + { + return side_strategy_type(); + } + + 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> + static inline typename point_in_geometry_strategy<Geometry1, Geometry2>::type + get_point_in_geometry_strategy() + { + typedef typename point_in_geometry_strategy + < + Geometry1, Geometry2 + >::type strategy_type; + return strategy_type(); + } + + template <typename Geometry> + struct area_strategy + { + typedef area::spherical + < + typename point_type<Geometry>::type, + CalculationType + > type; + }; + + template <typename Geometry> + static inline typename area_strategy<Geometry>::type get_area_strategy() + { + typedef typename area_strategy<Geometry>::type strategy_type; + return strategy_type(); + } + + template <typename Geometry> + struct distance_strategy + { + typedef distance::haversine + < + typename coordinate_type<Geometry>::type, + CalculationType + > type; + }; + + template <typename Geometry> + static inline typename distance_strategy<Geometry>::type get_distance_strategy() + { + typedef typename distance_strategy<Geometry>::type strategy_type; + return strategy_type(); + } enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 }; + // segment_intersection_info cannot outlive relate_ecef_segments template <typename CoordinateType, typename SegmentRatio, typename Vector3d> struct segment_intersection_info { @@ -83,6 +158,10 @@ struct relate_spherical_segments CoordinateType, double >::type promoted_type; + segment_intersection_info(CalcPolicy const& calc) + : calc_policy(calc) + {} + promoted_type comparable_length_a() const { return robust_ra.denominator(); @@ -110,7 +189,7 @@ struct relate_spherical_segments if (ip_flag == ipi_inters) { // TODO: assign the rest of coordinates - point = formula::cart3d_to_sph<Point>(intersection_point); + point = calc_policy.template from_cart3d<Point>(intersection_point); } else if (ip_flag == ipi_at_a1) { @@ -134,12 +213,21 @@ struct relate_spherical_segments SegmentRatio robust_ra; SegmentRatio robust_rb; intersection_point_flag ip_flag; + + CalcPolicy const& calc_policy; }; // Relate segments a and b - template <typename Segment1, typename Segment2, typename RobustPolicy> - static inline return_type apply(Segment1 const& a, Segment2 const& b, - RobustPolicy const& robust_policy) + template + < + typename Segment1, + typename Segment2, + typename Policy, + typename RobustPolicy + > + static inline typename Policy::return_type + apply(Segment1 const& a, Segment2 const& b, + Policy const& policy, RobustPolicy const& robust_policy) { typedef typename point_type<Segment1>::type point1_t; typedef typename point_type<Segment2>::type point2_t; @@ -152,15 +240,30 @@ struct relate_spherical_segments detail::assign_point_from_index<0>(b, b1); detail::assign_point_from_index<1>(b, b2); - return apply(a, b, robust_policy, a1, a2, b1, b2); + return apply(a, b, policy, robust_policy, a1, a2, b1, b2); } // Relate segments a and b - template <typename Segment1, typename Segment2, typename RobustPolicy, typename Point1, typename Point2> - static inline return_type apply(Segment1 const& a, Segment2 const& b, - RobustPolicy const&, - Point1 const& a1, Point1 const& a2, Point2 const& b1, Point2 const& b2) + template + < + typename Segment1, + typename Segment2, + typename Policy, + typename RobustPolicy, + typename Point1, + typename Point2 + > + static inline typename Policy::return_type + apply(Segment1 const& a, Segment2 const& b, + Policy const&, RobustPolicy const&, + Point1 const& a1, Point1 const& a2, Point2 const& b1, Point2 const& b2) { + // For now create it using default constructor. In the future it could + // be stored in strategy. However then apply() wouldn't be static and + // all relops and setops would have to take the strategy or model. + // Initialize explicitly to prevent compiler errors in case of PoD type + CalcPolicy const calc_policy = CalcPolicy(); + BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment1>) ); BOOST_CONCEPT_ASSERT( (concepts::ConstSegment<Segment2>) ); @@ -185,25 +288,29 @@ struct relate_spherical_segments typedef model::point<calc_t, 3, cs::cartesian> vec3d_t; - using namespace formula; - vec3d_t const a1v = sph_to_cart3d<vec3d_t>(a1); - vec3d_t const a2v = sph_to_cart3d<vec3d_t>(a2); - vec3d_t const b1v = sph_to_cart3d<vec3d_t>(b1); - vec3d_t const b2v = sph_to_cart3d<vec3d_t>(b2); + vec3d_t const a1v = calc_policy.template to_cart3d<vec3d_t>(a1); + vec3d_t const a2v = calc_policy.template to_cart3d<vec3d_t>(a2); + vec3d_t const b1v = calc_policy.template to_cart3d<vec3d_t>(b1); + vec3d_t const b2v = calc_policy.template to_cart3d<vec3d_t>(b2); - vec3d_t norm1 = cross_product(a1v, a2v); - vec3d_t norm2 = cross_product(b1v, b2v); - side_info sides; - // not normalized normals, the same as in SSF - sides.set<0>(sph_side_value(norm2, a1v), sph_side_value(norm2, a2v)); + + typename CalcPolicy::template plane<vec3d_t> + plane2 = calc_policy.get_plane(b1v, b2v); + + // not normalized normals, the same as in side strategy + sides.set<0>(plane2.side_value(a1v), plane2.side_value(a2v)); if (sides.same<0>()) { // Both points are at same side of other segment, we can leave return Policy::disjoint(); } - // not normalized normals, the same as in SSF - sides.set<1>(sph_side_value(norm1, b1v), sph_side_value(norm1, b2v)); + + typename CalcPolicy::template plane<vec3d_t> + plane1 = calc_policy.get_plane(a1v, a2v); + + // not normalized normals, the same as in side strategy + sides.set<1>(plane1.side_value(b1v), plane1.side_value(b2v)); if (sides.same<1>()) { // Both points are at same side of other segment, we can leave @@ -212,13 +319,10 @@ struct relate_spherical_segments // NOTE: at this point the segments may still be disjoint - bool collinear = sides.collinear(); - - calc_t const len1 = math::sqrt(dot_product(norm1, norm1)); - calc_t const len2 = math::sqrt(dot_product(norm2, norm2)); + calc_t len1, len2; - // point or opposite sides of a sphere, assume point - if (math::equals(len1, c0)) + // point or opposite sides of a sphere/spheroid, assume point + if (! detail::vec_normalize(plane1.normal, len1)) { a_is_point = true; if (sides.get<0, 0>() == 0 || sides.get<0, 1>() == 0) @@ -226,13 +330,8 @@ struct relate_spherical_segments sides.set<0>(0, 0); } } - else - { - // normalize - divide_value(norm1, len1); - } - if (math::equals(len2, c0)) + if (! detail::vec_normalize(plane2.normal, len2)) { b_is_point = true; if (sides.get<1, 0>() == 0 || sides.get<1, 1>() == 0) @@ -240,11 +339,6 @@ struct relate_spherical_segments sides.set<1>(0, 0); } } - else - { - // normalize - divide_value(norm2, len2); - } // check both degenerated once more if (a_is_point && b_is_point) @@ -255,16 +349,42 @@ struct relate_spherical_segments ; } + // NOTE: at this point the segments may still be disjoint // NOTE: at this point one of the segments may be degenerated - // and the segments may still be disjoint - calc_t dot_n1n2 = dot_product(norm1, norm2); + bool collinear = sides.collinear(); + + if (! collinear) + { + // NOTE: for some approximations it's possible that both points may lie + // on the same geodesic but still some of the sides may be != 0. + // This is e.g. true for long segments represented as elliptic arcs + // with origin different than the center of the coordinate system. + // So make the sides consistent + + // 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); + } + } + + calc_t dot_n1n2 = dot_product(plane1.normal, plane2.normal); // NOTE: this is technically not needed since theoretically above sides // are calculated, but just in case check the normals. // Have in mind that SSF side strategy doesn't check this. // collinear if normals are equal or opposite: cos(a) in {-1, 1} - if (!collinear && math::equals(math::abs(dot_n1n2), c1)) + if (! collinear && math::equals(math::abs(dot_n1n2), c1)) { collinear = true; sides.set<0>(0, 0); @@ -275,12 +395,12 @@ struct relate_spherical_segments { if (a_is_point) { - return collinear_one_degenerted<calc_t>(a, true, b1, b2, a1, a2, b1v, b2v, norm2, a1v); + return collinear_one_degenerated<Policy, calc_t>(a, true, b1, b2, a1, a2, b1v, b2v, plane2, a1v); } else if (b_is_point) { // b2 used to be consistent with (degenerated) checks above (is it needed?) - return collinear_one_degenerted<calc_t>(b, false, a1, a2, b1, b2, a1v, a2v, norm1, b1v); + return collinear_one_degenerated<Policy, calc_t>(b, false, a1, a2, b1, b2, a1v, a2v, plane1, b1v); } else { @@ -289,16 +409,16 @@ struct relate_spherical_segments // use shorter segment if (len1 <= len2) { - calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, norm1, b1v, dist_a1_a2, dist_a1_b1); - calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, norm1, b2v, dist_a1_a2, dist_a1_b2); + calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, plane1, b1v, dist_a1_a2, dist_a1_b1); + calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, plane1, b2v, 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, b1v, b2v, norm2, a1v, dist_b1_b2, dist_b1_a1); - calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, norm2, a2v, dist_b1_b2, dist_b1_a2); + calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, plane2, a1v, dist_b1_b2, dist_b1_a1); + calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, plane2, a2v, 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; @@ -359,7 +479,8 @@ struct relate_spherical_segments vec3d_t i1; 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, a1v, a2v, b1v, b2v, norm1, norm2, sides, + if (calculate_ip_data(a1, a2, b1, b2, a1v, a2v, b1v, b2v, + plane1, plane2, calc_policy, sides, i1, dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1, ip_flag)) { // intersects @@ -368,7 +489,7 @@ struct relate_spherical_segments calc_t, segment_ratio<calc_t>, vec3d_t - > sinfo; + > sinfo(calc_policy); sinfo.robust_ra.assign(dist_a1_i1, dist_a1_a2); sinfo.robust_rb.assign(dist_b1_i1, dist_b1_b2); @@ -385,30 +506,32 @@ struct relate_spherical_segments } private: - template <typename CalcT, typename Segment, typename Point1, typename Point2, typename Vec3d> - static inline return_type collinear_one_degenerted(Segment const& segment, bool degenerated_a, - Point1 const& a1, Point1 const& a2, - Point2 const& b1, Point2 const& b2, - Vec3d const& v1, Vec3d const& v2, Vec3d const& norm, - Vec3d const& vother) + template <typename Policy, typename CalcT, typename Segment, typename Point1, typename Point2, typename Vec3d, typename Plane> + 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, + Vec3d const& v1, Vec3d const& v2, + Plane const& plane, + Vec3d const& vother) { CalcT dist_1_2, dist_1_o; - return ! calculate_collinear_data(a1, a2, b1, b2, v1, v2, norm, vother, dist_1_2, dist_1_o) + return ! calculate_collinear_data(a1, a2, b1, b2, v1, v2, plane, vother, dist_1_2, dist_1_o) ? Policy::disjoint() : Policy::one_degenerate(segment, segment_ratio<CalcT>(dist_1_o, dist_1_2), degenerated_a); } - template <typename Point1, typename Point2, typename Vec3d, typename CalcT> - static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, - Point2 const& b1, Point2 const& b2, - Vec3d const& a1v, // in - Vec3d const& a2v, // in - Vec3d const& norm1, // in - Vec3d const& b1v_or_b2v, // in + template <typename Point1, typename Point2, typename Vec3d, typename Plane, typename CalcT> + static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in + Point2 const& b1, Point2 const& b2, // in + Vec3d const& a1v, // in + Vec3d const& a2v, // in + Plane const& plane1, // in + Vec3d const& b1v_or_b2v, // in CalcT& dist_a1_a2, CalcT& dist_a1_i1) // out { // calculate dist_a1_a2 and dist_a1_i1 - calculate_dists(a1v, a2v, norm1, b1v_or_b2v, dist_a1_a2, dist_a1_i1); + calculate_dists(a1v, a2v, plane1, b1v_or_b2v, dist_a1_a2, dist_a1_i1); // if i1 is close to a1 and b1 or b2 is equal to a1 if (is_endpoint_equal(dist_a1_i1, a1, b1, b2)) @@ -427,54 +550,53 @@ private: return segment_ratio<CalcT>(dist_a1_i1, dist_a1_a2).on_segment(); } - template <typename Point1, typename Point2, typename Vec3d, typename CalcT> + template <typename Point1, typename Point2, typename Vec3d, typename Plane, typename CalcT> static inline bool calculate_ip_data(Point1 const& a1, Point1 const& a2, // in Point2 const& b1, Point2 const& b2, // in Vec3d const& a1v, Vec3d const& a2v, // in Vec3d const& b1v, Vec3d const& b2v, // in - Vec3d const& norm1, Vec3d const& norm2, // in - side_info const& sides, // in - Vec3d & i1, // in-out - CalcT& dist_a1_a2, CalcT& dist_a1_i1, // out - CalcT& dist_b1_b2, CalcT& dist_b1_i1, // out + Plane const& plane1, // in + Plane const& plane2, // in + CalcPolicy const& calc_policy, // in + side_info const& sides, // in + Vec3d & ip, // 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 { - // great circles intersections - i1 = cross_product(norm1, norm2); - // NOTE: the length should be greater than 0 at this point - // if the normals were not normalized and their dot product - // not checked before this function is called the length - // should be checked here (math::equals(len, c0)) - CalcT const len = math::sqrt(dot_product(i1, i1)); - divide_value(i1, len); // normalize i1 - - calculate_dists(a1v, a2v, norm1, i1, dist_a1_a2, dist_a1_i1); + Vec3d ip1, ip2; + calc_policy.intersection_points(plane1, plane2, ip1, ip2); + calculate_dists(a1v, a2v, plane1, ip1, dist_a1_a2, dist_a1_ip); + ip = ip1; + // choose the opposite side of the globe if the distance is shorter { - CalcT const d = abs_distance(dist_a1_a2, dist_a1_i1); + CalcT const d = abs_distance(dist_a1_a2, dist_a1_ip); if (d > CalcT(0)) { - CalcT const dist_a1_i2 = dist_of_i2(dist_a1_i1); + // TODO: this should be ok not only for sphere + // but requires more investigation + CalcT const dist_a1_i2 = dist_of_i2(dist_a1_ip); CalcT const d2 = abs_distance(dist_a1_a2, dist_a1_i2); if (d2 < d) { - dist_a1_i1 = dist_a1_i2; - multiply_value(i1, CalcT(-1)); // the opposite intersection + dist_a1_ip = dist_a1_i2; + ip = ip2; } } } bool is_on_a = false, is_near_a1 = false, is_near_a2 = false; - if (! is_potentially_crossing(dist_a1_a2, dist_a1_i1, is_on_a, is_near_a1, is_near_a2)) + if (! is_potentially_crossing(dist_a1_a2, dist_a1_ip, is_on_a, is_near_a1, is_near_a2)) { return false; } - calculate_dists(b1v, b2v, norm2, i1, dist_b1_b2, dist_b1_i1); + calculate_dists(b1v, b2v, plane2, ip, dist_b1_b2, dist_b1_ip); bool is_on_b = false, is_near_b1 = false, is_near_b2 = false; - if (! is_potentially_crossing(dist_b1_b2, dist_b1_i1, is_on_b, is_near_b1, is_near_b2)) + if (! is_potentially_crossing(dist_b1_b2, dist_b1_ip, is_on_b, is_near_b1, is_near_b2)) { return false; } @@ -485,8 +607,8 @@ private: { if (is_near_b1 && equals_point_point(a1, b1)) { - dist_a1_i1 = 0; - dist_b1_i1 = 0; + dist_a1_ip = 0; + dist_b1_ip = 0; //i1 = a1v; ip_flag = ipi_at_a1; return true; @@ -494,8 +616,8 @@ private: if (is_near_b2 && equals_point_point(a1, b2)) { - dist_a1_i1 = 0; - dist_b1_i1 = dist_b1_b2; + dist_a1_ip = 0; + dist_b1_ip = dist_b1_b2; //i1 = a1v; ip_flag = ipi_at_a1; return true; @@ -506,8 +628,8 @@ private: { if (is_near_b1 && equals_point_point(a2, b1)) { - dist_a1_i1 = dist_a1_a2; - dist_b1_i1 = 0; + dist_a1_ip = dist_a1_a2; + dist_b1_ip = 0; //i1 = a2v; ip_flag = ipi_at_a2; return true; @@ -515,8 +637,8 @@ private: if (is_near_b2 && equals_point_point(a2, b2)) { - dist_a1_i1 = dist_a1_a2; - dist_b1_i1 = dist_b1_b2; + dist_a1_ip = dist_a1_a2; + dist_b1_ip = dist_b1_b2; //i1 = a2v; ip_flag = ipi_at_a2; return true; @@ -530,7 +652,7 @@ private: { if (is_near_b1 && sides.template get<1, 0>() == 0) // b1 wrt a { - dist_b1_i1 = 0; + dist_b1_ip = 0; //i1 = b1v; ip_flag = ipi_at_b1; return true; @@ -538,7 +660,7 @@ private: if (is_near_b2 && sides.template get<1, 1>() == 0) // b2 wrt a { - dist_b1_i1 = dist_b1_b2; + dist_b1_ip = dist_b1_b2; //i1 = b2v; ip_flag = ipi_at_b2; return true; @@ -549,7 +671,7 @@ private: { if (is_near_a1 && sides.template get<0, 0>() == 0) // a1 wrt b { - dist_a1_i1 = 0; + dist_a1_ip = 0; //i1 = a1v; ip_flag = ipi_at_a1; return true; @@ -557,7 +679,7 @@ private: if (is_near_a2 && sides.template get<0, 1>() == 0) // a2 wrt b { - dist_a1_i1 = dist_a1_a2; + dist_a1_ip = dist_a1_a2; //i1 = a2v; ip_flag = ipi_at_a2; return true; @@ -569,24 +691,26 @@ private: return is_on_a && is_on_b; } - template <typename Vec3d, typename CalcT> - static inline void calculate_dists(Vec3d const& a1v, // in - Vec3d const& a2v, // in - Vec3d const& norm1, // in - Vec3d const& i1, // in - CalcT& dist_a1_a2, CalcT& dist_a1_i1) // out + template <typename Vec3d, typename Plane, typename CalcT> + static inline void calculate_dists(Vec3d const& a1v, // in + Vec3d const& a2v, // in + Plane const& plane1, // in + Vec3d const& i1, // in + CalcT& dist_a1_a2, // out + CalcT& dist_a1_i1) // out { - CalcT const c0 = 0; + //CalcT const c0 = 0; CalcT const c1 = 1; CalcT const c2 = 2; CalcT const c4 = 4; - CalcT cos_a1_a2 = dot_product(a1v, a2v); + CalcT cos_a1_a2 = plane1.cos_angle_between(a1v, a2v); dist_a1_a2 = -cos_a1_a2 + c1; // [1, -1] -> [0, 2] representing [0, pi] - CalcT cos_a1_i1 = dot_product(a1v, i1); + bool is_forward = true; + CalcT cos_a1_i1 = plane1.cos_angle_between(a1v, i1, is_forward); dist_a1_i1 = -cos_a1_i1 + c1; // [0, 2] representing [0, pi] - if (dot_product(norm1, cross_product(a1v, i1)) < c0) // left or right of a1 on a + if (! is_forward) // left or right of a1 on a { dist_a1_i1 = -dist_a1_i1; // [0, 2] -> [0, -2] representing [0, -pi] } @@ -666,27 +790,117 @@ private: } }; +struct spherical_segments_calc_policy +{ + template <typename Point, typename Point3d> + static Point from_cart3d(Point3d const& point_3d) + { + return formula::cart3d_to_sph<Point>(point_3d); + } + + template <typename Point3d, typename Point> + static Point3d to_cart3d(Point const& point) + { + return formula::sph_to_cart3d<Point3d>(point); + } + + 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); + } + + static coord_t cos_angle_between(Point3d const& p1, Point3d const& p2) + { + return dot_product(p1, p2); + } + + coord_t cos_angle_between(Point3d const& p1, Point3d const& p2, bool & is_forward) const + { + coord_t const c0 = 0; + is_forward = dot_product(normal, cross_product(p1, p2)) >= c0; + return dot_product(p1, p2); + } + + Point3d normal; + }; + + template <typename Point3d> + static plane<Point3d> get_plane(Point3d const& p1, Point3d const& p2) + { + return plane<Point3d>(p1, p2); + } + + template <typename Point3d> + static bool intersection_points(plane<Point3d> const& plane1, + plane<Point3d> const& plane2, + Point3d & ip1, Point3d & ip2) + { + typedef typename coordinate_type<Point3d>::type coord_t; + + ip1 = cross_product(plane1.normal, plane2.normal); + // NOTE: the length should be greater than 0 at this point + // if the normals were not normalized and their dot product + // not checked before this function is called the length + // should be checked here (math::equals(len, c0)) + coord_t const len = math::sqrt(dot_product(ip1, ip1)); + divide_value(ip1, len); // normalize i1 + + ip2 = ip1; + multiply_value(ip2, coord_t(-1)); + + return true; + } +}; + + +template +< + typename CalculationType = void +> +struct spherical_segments + : ecef_segments + < + spherical_segments_calc_policy, + CalculationType + > +{}; + #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS namespace services { -/*template <typename Policy, typename CalculationType> -struct default_strategy<spherical_polar_tag, Policy, CalculationType> +/*template <typename CalculationType> +struct default_strategy<spherical_polar_tag, CalculationType> { - typedef relate_spherical_segments<Policy, CalculationType> type; + typedef spherical_segments<CalculationType> type; };*/ -template <typename Policy, typename CalculationType> -struct default_strategy<spherical_equatorial_tag, Policy, CalculationType> +template <typename CalculationType> +struct default_strategy<spherical_equatorial_tag, CalculationType> { - typedef relate_spherical_segments<Policy, CalculationType> type; + typedef spherical_segments<CalculationType> type; }; -template <typename Policy, typename CalculationType> -struct default_strategy<geographic_tag, Policy, CalculationType> +template <typename CalculationType> +struct default_strategy<geographic_tag, CalculationType> { - typedef relate_spherical_segments<Policy, CalculationType> type; + // NOTE: Spherical strategy returns the same result as the geographic one + // representing segments as great elliptic arcs. If the elliptic arcs are + // not great elliptic arcs (the origin not in the center of the coordinate + // system) then there may be problems with consistency of the side and + // intersection strategies. + typedef spherical_segments<CalculationType> type; }; } // namespace services @@ -695,6 +909,71 @@ struct default_strategy<geographic_tag, Policy, CalculationType> }} // namespace strategy::intersection + +namespace strategy +{ + +namespace within { namespace services +{ + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +}} // within::services + +namespace covered_by { namespace services +{ + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, linear_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, linear_tag, polygonal_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, linear_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +template <typename Geometry1, typename Geometry2, typename AnyTag1, typename AnyTag2> +struct default_strategy<Geometry1, Geometry2, AnyTag1, AnyTag2, polygonal_tag, polygonal_tag, spherical_tag, spherical_tag> +{ + typedef strategy::intersection::spherical_segments<> type; +}; + +}} // within::services + +} // strategy + + }} // namespace boost::geometry diff --git a/boost/geometry/strategies/strategies.hpp b/boost/geometry/strategies/strategies.hpp index 342485cc4c..f3c6787a4b 100644 --- a/boost/geometry/strategies/strategies.hpp +++ b/boost/geometry/strategies/strategies.hpp @@ -4,9 +4,10 @@ // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. // Copyright (c) 2009-2012 Mateusz Loskot, London, UK. -// This file was modified by Oracle on 2014-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 Vissarion Fysikopoulos, on behalf of Oracle // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library @@ -23,18 +24,24 @@ #include <boost/geometry/strategies/tags.hpp> #include <boost/geometry/strategies/area.hpp> +#include <boost/geometry/strategies/azimuth.hpp> #include <boost/geometry/strategies/buffer.hpp> #include <boost/geometry/strategies/centroid.hpp> #include <boost/geometry/strategies/compare.hpp> #include <boost/geometry/strategies/convex_hull.hpp> +#include <boost/geometry/strategies/covered_by.hpp> +#include <boost/geometry/strategies/disjoint.hpp> #include <boost/geometry/strategies/distance.hpp> +#include <boost/geometry/strategies/envelope.hpp> #include <boost/geometry/strategies/intersection.hpp> #include <boost/geometry/strategies/intersection_strategies.hpp> // for backward compatibility +#include <boost/geometry/strategies/relate.hpp> #include <boost/geometry/strategies/side.hpp> #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/azimuth.hpp> #include <boost/geometry/strategies/cartesian/box_in_box.hpp> #include <boost/geometry/strategies/cartesian/buffer_end_flat.hpp> #include <boost/geometry/strategies/cartesian/buffer_end_round.hpp> @@ -47,30 +54,42 @@ #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/disjoint_segment_box.hpp> #include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp> #include <boost/geometry/strategies/cartesian/distance_pythagoras_point_box.hpp> #include <boost/geometry/strategies/cartesian/distance_pythagoras_box_box.hpp> #include <boost/geometry/strategies/cartesian/distance_projected_point.hpp> #include <boost/geometry/strategies/cartesian/distance_projected_point_ax.hpp> +#include <boost/geometry/strategies/cartesian/envelope_segment.hpp> +#include <boost/geometry/strategies/cartesian/intersection.hpp> #include <boost/geometry/strategies/cartesian/point_in_box.hpp> #include <boost/geometry/strategies/cartesian/point_in_poly_franklin.hpp> #include <boost/geometry/strategies/cartesian/point_in_poly_crossings_multiply.hpp> #include <boost/geometry/strategies/cartesian/side_by_triangle.hpp> -#include <boost/geometry/strategies/spherical/area_huiller.hpp> +#include <boost/geometry/strategies/spherical/area.hpp> +#include <boost/geometry/strategies/spherical/azimuth.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_point_box.hpp> #include <boost/geometry/strategies/spherical/compare_circular.hpp> +#include <boost/geometry/strategies/spherical/envelope_segment.hpp> #include <boost/geometry/strategies/spherical/intersection.hpp> #include <boost/geometry/strategies/spherical/ssf.hpp> +#include <boost/geometry/strategies/geographic/area.hpp> +#include <boost/geometry/strategies/geographic/azimuth.hpp> +#include <boost/geometry/strategies/geographic/distance.hpp> #include <boost/geometry/strategies/geographic/distance_andoyer.hpp> #include <boost/geometry/strategies/geographic/distance_thomas.hpp> #include <boost/geometry/strategies/geographic/distance_vincenty.hpp> -//#include <boost/geometry/strategies/geographic/side_andoyer.hpp> -//#include <boost/geometry/strategies/geographic/side_thomas.hpp> -//#include <boost/geometry/strategies/geographic/side_vincenty.hpp> +#include <boost/geometry/strategies/geographic/envelope_segment.hpp> +#include <boost/geometry/strategies/geographic/intersection.hpp> +//#include <boost/geometry/strategies/geographic/intersection_elliptic.hpp> +#include <boost/geometry/strategies/geographic/side.hpp> +#include <boost/geometry/strategies/geographic/side_andoyer.hpp> +#include <boost/geometry/strategies/geographic/side_thomas.hpp> +#include <boost/geometry/strategies/geographic/side_vincenty.hpp> #include <boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp> #include <boost/geometry/strategies/agnostic/buffer_distance_asymmetric.hpp> @@ -80,8 +99,6 @@ #include <boost/geometry/strategies/agnostic/point_in_poly_winding.hpp> #include <boost/geometry/strategies/agnostic/simplify_douglas_peucker.hpp> -#include <boost/geometry/strategies/agnostic/relate.hpp> - #include <boost/geometry/strategies/strategy_transform.hpp> #include <boost/geometry/strategies/transform/matrix_transformers.hpp> diff --git a/boost/geometry/strategies/transform/inverse_transformer.hpp b/boost/geometry/strategies/transform/inverse_transformer.hpp index e64a46e4a8..5213bce47a 100644 --- a/boost/geometry/strategies/transform/inverse_transformer.hpp +++ b/boost/geometry/strategies/transform/inverse_transformer.hpp @@ -14,15 +14,8 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_TRANSFORM_INVERSE_TRANSFORMER_HPP #define BOOST_GEOMETRY_STRATEGIES_TRANSFORM_INVERSE_TRANSFORMER_HPP -// Remove the ublas checking, otherwise the inverse might fail -// (while nothing seems to be wrong) -#ifdef BOOST_UBLAS_TYPE_CHECK -#undef BOOST_UBLAS_TYPE_CHECK -#endif -#define BOOST_UBLAS_TYPE_CHECK 0 - -#include <boost/numeric/ublas/lu.hpp> -#include <boost/numeric/ublas/io.hpp> +#include <boost/qvm/mat.hpp> +#include <boost/qvm/mat_operations.hpp> #include <boost/geometry/strategies/transform/matrix_transformers.hpp> @@ -44,31 +37,13 @@ template std::size_t Dimension2 > class inverse_transformer - : public ublas_transformer<CalculationType, Dimension1, Dimension2> + : public matrix_transformer<CalculationType, Dimension1, Dimension2> { public : template <typename Transformer> inline inverse_transformer(Transformer const& input) { - typedef boost::numeric::ublas::matrix<CalculationType> matrix_type; - - // create a working copy of the input - matrix_type copy(input.matrix()); - - // create a permutation matrix for the LU-factorization - typedef boost::numeric::ublas::permutation_matrix<> permutation_matrix; - permutation_matrix pm(copy.size1()); - - // perform LU-factorization - int res = boost::numeric::ublas::lu_factorize<matrix_type>(copy, pm); - if( res == 0 ) - { - // create identity matrix - this->m_matrix.assign(boost::numeric::ublas::identity_matrix<CalculationType>(copy.size1())); - - // backsubstitute to get the inverse - boost::numeric::ublas::lu_substitute(copy, pm, this->m_matrix); - } + this->m_matrix = boost::qvm::inverse(input.matrix()); } }; diff --git a/boost/geometry/strategies/transform/map_transformer.hpp b/boost/geometry/strategies/transform/map_transformer.hpp index 1109e814b9..01ea6168c5 100644 --- a/boost/geometry/strategies/transform/map_transformer.hpp +++ b/boost/geometry/strategies/transform/map_transformer.hpp @@ -46,9 +46,10 @@ template bool SameScale = true > class map_transformer - : public ublas_transformer<CalculationType, Dimension1, Dimension2> + : public matrix_transformer<CalculationType, Dimension1, Dimension2> { - typedef boost::numeric::ublas::matrix<CalculationType> M; + typedef boost::qvm::mat<CalculationType, Dimension1 + 1, Dimension2 + 1> M; + typedef boost::qvm::mat<CalculationType, 3, 3> matrix33; public : template <typename B, typename D> @@ -76,26 +77,26 @@ private : { // Translate to a coordinate system centered on world coordinates (-wx, -wy) - M t1(3,3); - t1(0,0) = 1; t1(0,1) = 0; t1(0,2) = -wx; - t1(1,0) = 0; t1(1,1) = 1; t1(1,2) = -wy; - t1(2,0) = 0; t1(2,1) = 0; t1(2,2) = 1; + matrix33 t1; + qvm::A<0,0>(t1) = 1; qvm::A<0,1>(t1) = 0; qvm::A<0,2>(t1) = -wx; + qvm::A<1,0>(t1) = 0; qvm::A<1,1>(t1) = 1; qvm::A<1,2>(t1) = -wy; + qvm::A<2,0>(t1) = 0; qvm::A<2,1>(t1) = 0; qvm::A<2,2>(t1) = 1; // Scale the map - M s(3,3); - s(0,0) = scalex; s(0,1) = 0; s(0,2) = 0; - s(1,0) = 0; s(1,1) = scaley; s(1,2) = 0; - s(2,0) = 0; s(2,1) = 0; s(2,2) = 1; + matrix33 s; + qvm::A<0,0>(s) = scalex; qvm::A<0,1>(s) = 0; qvm::A<0,2>(s) = 0; + qvm::A<1,0>(s) = 0; qvm::A<1,1>(s) = scaley; qvm::A<1,2>(s) = 0; + qvm::A<2,0>(s) = 0; qvm::A<2,1>(s) = 0; qvm::A<2,2>(s) = 1; // Translate to a coordinate system centered on the specified pixels (+px, +py) - M t2(3, 3); - t2(0,0) = 1; t2(0,1) = 0; t2(0,2) = px; - t2(1,0) = 0; t2(1,1) = 1; t2(1,2) = py; - t2(2,0) = 0; t2(2,1) = 0; t2(2,2) = 1; + matrix33 t2; + qvm::A<0,0>(t2) = 1; qvm::A<0,1>(t2) = 0; qvm::A<0,2>(t2) = px; + qvm::A<1,0>(t2) = 0; qvm::A<1,1>(t2) = 1; qvm::A<1,2>(t2) = py; + qvm::A<2,0>(t2) = 0; qvm::A<2,1>(t2) = 0; qvm::A<2,2>(t2) = 1; // Calculate combination matrix in two steps - this->m_matrix = boost::numeric::ublas::prod(s, t1); - this->m_matrix = boost::numeric::ublas::prod(t2, this->m_matrix); + this->m_matrix = s * t1; + this->m_matrix = t2 * this->m_matrix; } @@ -140,20 +141,20 @@ private : if (Mirror) { // Mirror in y-direction - M m(3,3); - m(0,0) = 1; m(0,1) = 0; m(0,2) = 0; - m(1,0) = 0; m(1,1) = -1; m(1,2) = 0; - m(2,0) = 0; m(2,1) = 0; m(2,2) = 1; + matrix33 m; + qvm::A<0,0>(m) = 1; qvm::A<0,1>(m) = 0; qvm::A<0,2>(m) = 0; + qvm::A<1,0>(m) = 0; qvm::A<1,1>(m) = -1; qvm::A<1,2>(m) = 0; + qvm::A<2,0>(m) = 0; qvm::A<2,1>(m) = 0; qvm::A<2,2>(m) = 1; // Translate in y-direction such that it fits again - M y(3, 3); - y(0,0) = 1; y(0,1) = 0; y(0,2) = 0; - y(1,0) = 0; y(1,1) = 1; y(1,2) = height; - y(2,0) = 0; y(2,1) = 0; y(2,2) = 1; + matrix33 y; + qvm::A<0,0>(y) = 1; qvm::A<0,1>(y) = 0; qvm::A<0,2>(y) = 0; + qvm::A<1,0>(y) = 0; qvm::A<1,1>(y) = 1; qvm::A<1,2>(y) = height; + qvm::A<2,0>(y) = 0; qvm::A<2,1>(y) = 0; qvm::A<2,2>(y) = 1; // Calculate combination matrix in two steps - this->m_matrix = boost::numeric::ublas::prod(m, this->m_matrix); - this->m_matrix = boost::numeric::ublas::prod(y, this->m_matrix); + this->m_matrix = m * this->m_matrix; + this->m_matrix = y * this->m_matrix; } } }; diff --git a/boost/geometry/strategies/transform/matrix_transformers.hpp b/boost/geometry/strategies/transform/matrix_transformers.hpp index d891263a7d..e0ac6496f3 100644 --- a/boost/geometry/strategies/transform/matrix_transformers.hpp +++ b/boost/geometry/strategies/transform/matrix_transformers.hpp @@ -22,27 +22,9 @@ #include <cstddef> -// Remove the ublas checking, otherwise the inverse might fail -// (while nothing seems to be wrong) -#ifdef BOOST_UBLAS_TYPE_CHECK -#undef BOOST_UBLAS_TYPE_CHECK -#endif -#define BOOST_UBLAS_TYPE_CHECK 0 - -#include <boost/numeric/conversion/cast.hpp> - -#if defined(__clang__) -// Avoid warning about unused UBLAS function: boost_numeric_ublas_abs -#pragma clang diagnostic push -#pragma clang diagnostic ignored "-Wunused-function" -#endif - -#include <boost/numeric/ublas/vector.hpp> -#include <boost/numeric/ublas/matrix.hpp> - -#if defined(__clang__) -#pragma clang diagnostic pop -#endif +#include <boost/qvm/mat.hpp> +#include <boost/qvm/mat_access.hpp> +#include <boost/qvm/mat_operations.hpp> #include <boost/geometry/core/access.hpp> #include <boost/geometry/core/coordinate_dimension.hpp> @@ -75,38 +57,37 @@ template std::size_t Dimension1, std::size_t Dimension2 > -class ublas_transformer +class matrix_transformer { }; template <typename CalculationType> -class ublas_transformer<CalculationType, 2, 2> +class matrix_transformer<CalculationType, 2, 2> { protected : typedef CalculationType ct; - typedef boost::numeric::ublas::matrix<ct> matrix_type; + typedef boost::qvm::mat<ct, 3, 3> matrix_type; matrix_type m_matrix; public : - inline ublas_transformer( + inline matrix_transformer( ct const& m_0_0, ct const& m_0_1, ct const& m_0_2, ct const& m_1_0, ct const& m_1_1, ct const& m_1_2, ct const& m_2_0, ct const& m_2_1, ct const& m_2_2) - : m_matrix(3, 3) { - m_matrix(0,0) = m_0_0; m_matrix(0,1) = m_0_1; m_matrix(0,2) = m_0_2; - m_matrix(1,0) = m_1_0; m_matrix(1,1) = m_1_1; m_matrix(1,2) = m_1_2; - m_matrix(2,0) = m_2_0; m_matrix(2,1) = m_2_1; m_matrix(2,2) = m_2_2; + qvm::A<0,0>(m_matrix) = m_0_0; qvm::A<0,1>(m_matrix) = m_0_1; qvm::A<0,2>(m_matrix) = m_0_2; + qvm::A<1,0>(m_matrix) = m_1_0; qvm::A<1,1>(m_matrix) = m_1_1; qvm::A<1,2>(m_matrix) = m_1_2; + qvm::A<2,0>(m_matrix) = m_2_0; qvm::A<2,1>(m_matrix) = m_2_1; qvm::A<2,2>(m_matrix) = m_2_2; } - inline ublas_transformer(matrix_type const& matrix) + inline matrix_transformer(matrix_type const& matrix) : m_matrix(matrix) {} - inline ublas_transformer() : m_matrix(3, 3) {} + inline matrix_transformer() {} template <typename P1, typename P2> inline bool apply(P1 const& p1, P2& p2) const @@ -117,8 +98,8 @@ public : ct const& c1 = get<0>(p1); ct const& c2 = get<1>(p1); - ct p2x = c1 * m_matrix(0,0) + c2 * m_matrix(0,1) + m_matrix(0,2); - ct p2y = c1 * m_matrix(1,0) + c2 * m_matrix(1,1) + m_matrix(1,2); + ct p2x = c1 * qvm::A<0,0>(m_matrix) + c2 * qvm::A<0,1>(m_matrix) + qvm::A<0,2>(m_matrix); + ct p2y = c1 * qvm::A<1,0>(m_matrix) + c2 * qvm::A<1,1>(m_matrix) + qvm::A<1,2>(m_matrix); typedef typename geometry::coordinate_type<P2>::type ct2; set<0>(p2, boost::numeric_cast<ct2>(p2x)); @@ -133,51 +114,50 @@ public : // It IS possible to go from 3 to 2 coordinates template <typename CalculationType> -class ublas_transformer<CalculationType, 3, 2> : public ublas_transformer<CalculationType, 2, 2> +class matrix_transformer<CalculationType, 3, 2> : public matrix_transformer<CalculationType, 2, 2> { typedef CalculationType ct; public : - inline ublas_transformer( + inline matrix_transformer( ct const& m_0_0, ct const& m_0_1, ct const& m_0_2, ct const& m_1_0, ct const& m_1_1, ct const& m_1_2, ct const& m_2_0, ct const& m_2_1, ct const& m_2_2) - : ublas_transformer<CalculationType, 2, 2>( + : matrix_transformer<CalculationType, 2, 2>( m_0_0, m_0_1, m_0_2, m_1_0, m_1_1, m_1_2, m_2_0, m_2_1, m_2_2) {} - inline ublas_transformer() - : ublas_transformer<CalculationType, 2, 2>() + inline matrix_transformer() + : matrix_transformer<CalculationType, 2, 2>() {} }; template <typename CalculationType> -class ublas_transformer<CalculationType, 3, 3> +class matrix_transformer<CalculationType, 3, 3> { protected : typedef CalculationType ct; - typedef boost::numeric::ublas::matrix<ct> matrix_type; + typedef boost::qvm::mat<ct, 4, 4> matrix_type; matrix_type m_matrix; public : - inline ublas_transformer( + inline matrix_transformer( ct const& m_0_0, ct const& m_0_1, ct const& m_0_2, ct const& m_0_3, ct const& m_1_0, ct const& m_1_1, ct const& m_1_2, ct const& m_1_3, ct const& m_2_0, ct const& m_2_1, ct const& m_2_2, ct const& m_2_3, ct const& m_3_0, ct const& m_3_1, ct const& m_3_2, ct const& m_3_3 ) - : m_matrix(4, 4) { - m_matrix(0,0) = m_0_0; m_matrix(0,1) = m_0_1; m_matrix(0,2) = m_0_2; m_matrix(0,3) = m_0_3; - m_matrix(1,0) = m_1_0; m_matrix(1,1) = m_1_1; m_matrix(1,2) = m_1_2; m_matrix(1,3) = m_1_3; - m_matrix(2,0) = m_2_0; m_matrix(2,1) = m_2_1; m_matrix(2,2) = m_2_2; m_matrix(2,3) = m_2_3; - m_matrix(3,0) = m_3_0; m_matrix(3,1) = m_3_1; m_matrix(3,2) = m_3_2; m_matrix(3,3) = m_3_3; + qvm::A<0,0>(m_matrix) = m_0_0; qvm::A<0,1>(m_matrix) = m_0_1; qvm::A<0,2>(m_matrix) = m_0_2; qvm::A<0,3>(m_matrix) = m_0_3; + qvm::A<1,0>(m_matrix) = m_1_0; qvm::A<1,1>(m_matrix) = m_1_1; qvm::A<1,2>(m_matrix) = m_1_2; qvm::A<1,3>(m_matrix) = m_1_3; + qvm::A<2,0>(m_matrix) = m_2_0; qvm::A<2,1>(m_matrix) = m_2_1; qvm::A<2,2>(m_matrix) = m_2_2; qvm::A<2,3>(m_matrix) = m_2_3; + qvm::A<3,0>(m_matrix) = m_3_0; qvm::A<3,1>(m_matrix) = m_3_1; qvm::A<3,2>(m_matrix) = m_3_2; qvm::A<3,3>(m_matrix) = m_3_3; } - inline ublas_transformer() : m_matrix(4, 4) {} + inline matrix_transformer() {} template <typename P1, typename P2> inline bool apply(P1 const& p1, P2& p2) const @@ -222,7 +202,7 @@ class translate_transformer template<typename CalculationType> -class translate_transformer<CalculationType, 2, 2> : public ublas_transformer<CalculationType, 2, 2> +class translate_transformer<CalculationType, 2, 2> : public matrix_transformer<CalculationType, 2, 2> { public : // To have translate transformers compatible for 2/3 dimensions, the @@ -230,7 +210,7 @@ public : inline translate_transformer(CalculationType const& translate_x, CalculationType const& translate_y, CalculationType const& = 0) - : ublas_transformer<CalculationType, 2, 2>( + : matrix_transformer<CalculationType, 2, 2>( 1, 0, translate_x, 0, 1, translate_y, 0, 0, 1) @@ -239,13 +219,13 @@ public : template <typename CalculationType> -class translate_transformer<CalculationType, 3, 3> : public ublas_transformer<CalculationType, 3, 3> +class translate_transformer<CalculationType, 3, 3> : public matrix_transformer<CalculationType, 3, 3> { public : inline translate_transformer(CalculationType const& translate_x, CalculationType const& translate_y, CalculationType const& translate_z) - : ublas_transformer<CalculationType, 3, 3>( + : matrix_transformer<CalculationType, 3, 3>( 1, 0, 0, translate_x, 0, 1, 0, translate_y, 0, 0, 1, translate_z, @@ -275,14 +255,14 @@ class scale_transformer template <typename CalculationType> -class scale_transformer<CalculationType, 2, 2> : public ublas_transformer<CalculationType, 2, 2> +class scale_transformer<CalculationType, 2, 2> : public matrix_transformer<CalculationType, 2, 2> { public : inline scale_transformer(CalculationType const& scale_x, CalculationType const& scale_y, CalculationType const& = 0) - : ublas_transformer<CalculationType, 2, 2>( + : matrix_transformer<CalculationType, 2, 2>( scale_x, 0, 0, 0, scale_y, 0, 0, 0, 1) @@ -290,7 +270,7 @@ public : inline scale_transformer(CalculationType const& scale) - : ublas_transformer<CalculationType, 2, 2>( + : matrix_transformer<CalculationType, 2, 2>( scale, 0, 0, 0, scale, 0, 0, 0, 1) @@ -299,13 +279,13 @@ public : template <typename CalculationType> -class scale_transformer<CalculationType, 3, 3> : public ublas_transformer<CalculationType, 3, 3> +class scale_transformer<CalculationType, 3, 3> : public matrix_transformer<CalculationType, 3, 3> { public : inline scale_transformer(CalculationType const& scale_x, CalculationType const& scale_y, CalculationType const& scale_z) - : ublas_transformer<CalculationType, 3, 3>( + : matrix_transformer<CalculationType, 3, 3>( scale_x, 0, 0, 0, 0, scale_y, 0, 0, 0, 0, scale_z, 0, @@ -314,7 +294,7 @@ public : inline scale_transformer(CalculationType const& scale) - : ublas_transformer<CalculationType, 3, 3>( + : matrix_transformer<CalculationType, 3, 3>( scale, 0, 0, 0, 0, scale, 0, 0, 0, 0, scale, 0, @@ -363,11 +343,11 @@ template std::size_t Dimension2 > class rad_rotate_transformer - : public ublas_transformer<CalculationType, Dimension1, Dimension2> + : public matrix_transformer<CalculationType, Dimension1, Dimension2> { public : inline rad_rotate_transformer(CalculationType const& angle) - : ublas_transformer<CalculationType, Dimension1, Dimension2>( + : matrix_transformer<CalculationType, Dimension1, Dimension2>( cos(angle), sin(angle), 0, -sin(angle), cos(angle), 0, 0, 0, 1) diff --git a/boost/geometry/strategies/within.hpp b/boost/geometry/strategies/within.hpp index d625bc40e6..28a3283db0 100644 --- a/boost/geometry/strategies/within.hpp +++ b/boost/geometry/strategies/within.hpp @@ -4,6 +4,11 @@ // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. // Copyright (c) 2009-2012 Mateusz Loskot, London, UK. +// This file was modified by Oracle on 2017. +// Modifications copyright (c) 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. @@ -16,6 +21,12 @@ #include <boost/mpl/assert.hpp> +#include <boost/geometry/core/cs.hpp> +#include <boost/geometry/core/point_type.hpp> +#include <boost/geometry/core/tag.hpp> +#include <boost/geometry/core/tags.hpp> +#include <boost/geometry/core/tag_cast.hpp> + namespace boost { namespace geometry { @@ -30,23 +41,39 @@ namespace services /*! \brief Traits class binding a within determination strategy to a coordinate system \ingroup within -\tparam TagContained tag (possibly casted) of point-type -\tparam TagContained tag (possibly casted) of (possibly) containing type -\tparam CsTagContained tag of coordinate system of point-type -\tparam CsTagContaining tag of coordinate system of (possibly) containing type -\tparam Geometry geometry-type of input (often point, or box) +\tparam GeometryContained geometry-type of input (possibly) contained type \tparam GeometryContaining geometry-type of input (possibly) containing type +\tparam TagContained casted tag of (possibly) contained type +\tparam TagContaining casted tag of (possibly) containing type +\tparam CsTagContained tag of coordinate system of (possibly) contained type +\tparam CsTagContaining tag of coordinate system of (possibly) containing type */ template < - typename TagContained, - typename TagContaining, - typename CastedTagContained, - typename CastedTagContaining, - typename CsTagContained, - typename CsTagContaining, typename GeometryContained, - typename GeometryContaining + typename GeometryContaining, + typename TagContained = typename tag<GeometryContained>::type, + typename TagContaining = typename tag<GeometryContaining>::type, + typename CastedTagContained = typename tag_cast + < + typename tag<GeometryContained>::type, + pointlike_tag, linear_tag, polygonal_tag, areal_tag + >::type, + typename CastedTagContaining = typename tag_cast + < + typename tag<GeometryContaining>::type, + pointlike_tag, linear_tag, polygonal_tag, areal_tag + >::type, + typename CsTagContained = typename tag_cast + < + typename cs_tag<typename point_type<GeometryContained>::type>::type, + spherical_tag + >::type, + typename CsTagContaining = typename tag_cast + < + typename cs_tag<typename point_type<GeometryContaining>::type>::type, + spherical_tag + >::type > struct default_strategy { |