diff --git a/doc/doxy/Doxyfile b/doc/doxy/Doxyfile index 4c9aa41d34..aab8e217bd 100644 --- a/doc/doxy/Doxyfile +++ b/doc/doxy/Doxyfile @@ -211,6 +211,7 @@ INPUT = . .. ../../../../boost/geometry/core \ ../../../../boost/geometry/strategies/agnostic \ ../../../../boost/geometry/strategies/cartesian \ ../../../../boost/geometry/strategies/spherical \ + ../../../../boost/geometry/strategies/geographic \ ../../../../boost/geometry/strategies/transform \ ../../../../boost/geometry/util \ ../../../../boost/geometry/views \ diff --git a/doc/index/rtree/query.qbk b/doc/index/rtree/query.qbk index 0cf80acc78..75f40ec293 100644 --- a/doc/index/rtree/query.qbk +++ b/doc/index/rtree/query.qbk @@ -1,7 +1,7 @@ [/============================================================================ Boost.Geometry Index - Copyright (c) 2011-2013 Adam Wulkiewicz. + Copyright (c) 2011-2017 Adam Wulkiewicz. Use, modification and distribution is subject to the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at @@ -57,7 +57,7 @@ Query iterators returned by free functions __box__ box_region(...); std::copy(bgi::qbegin(rt, bgi::intersects(box_region)), bgi::qend(rt), std::back_inserter(returned_values)); -[h4 Spatial predicates] +[h4 Spatial queries] Queries using spatial predicates returns `__value__`s which are related somehow to some Geometry - box, polygon, etc. Names of spatial predicates correspond to names of __boost_geometry__ algorithms (boolean operations). @@ -129,6 +129,9 @@ The same way different query Geometries can be used: Segment seg(/*...*/); rt.query(bgi::nearest(seg, k), std::back_inserter(returned_values)); +[note In case of k-NN queries performed with `query()` function it's not guaranteed that the returned values will be sorted according to the distance. + It's different in case of k-NN queries performed with query iterator returned by `qbegin()` function which guarantees the iteration over the closest `__value__`s first. ] + [h4 User-defined unary predicate] The user may pass a `UnaryPredicate` - function, function object or lambda expression taking const reference to Value and returning bool. @@ -175,7 +178,7 @@ may use `index::satisfies()` function like on the example below: rt.query(index::intersects(box) && !index::satisfies(is_not_red), std::back_inserter(result)); -[h4 Passing a set of predicates] +[h4 Passing set of predicates] It's possible to use some number of predicates in one query by connecting them with `operator&&` e.g. `Pred1 && Pred2 && Pred3 && ...`. @@ -190,7 +193,7 @@ left-to-right so placing most restrictive predicates first should accelerate the rt.query(index::intersects(box1) && !index::within(box2) && index::overlaps(box3), std::back_inserter(result)); -Of course it's possible to connect different types of predicates together. +It's possible to connect different types of predicates together. index::query(rt, index::nearest(pt, k) && index::within(b), std::back_inserter(returned_values)); @@ -211,13 +214,11 @@ or may be stopped at some point, when all interesting values were gathered. The break; } -[note In the case of iterative k-NN queries it's guaranteed to iterate over the closest `__value__`s first. ] - [warning The modification of the `rtree`, e.g. insertion or removal of `__value__`s may invalidate the iterators. ] -[h4 Inserting query results into the other R-tree] +[h4 Inserting query results into another R-tree] -There are several ways of inserting Values returned by a query to the other R-tree container. +There are several ways of inserting Values returned by a query into another R-tree container. The most basic way is creating a temporary container for Values and insert them later. namespace bgi = boost::geometry::index; @@ -231,14 +232,13 @@ The most basic way is creating a temporary container for Values and insert them rt1.query(bgi::intersects(Box(/*...*/)), std::back_inserter(result)); RTree rt2(result.begin(), result.end()); -However there are better ways. One of these methods is mentioned in the "Creation and modification" section. +However there are other ways. One of these methods is mentioned in the "Creation and modification" section. The insert iterator may be passed directly into the query. RTree rt3; rt1.query(bgi::intersects(Box(/*...*/))), bgi::inserter(rt3)); -If you're a user of Boost.Range you'll appreciate the third option. You may pass the result Range directly into the -constructor or `insert()` member function. +You may also pass the resulting Range directly into the constructor or `insert()` member function using Boost.Range adaptor syntax. RTree rt4(rt1 | bgi::adaptors::queried(bgi::intersects(Box(/*...*/))))); diff --git a/doc/make_qbk.py b/doc/make_qbk.py index db6cc0ab94..184e422230 100755 --- a/doc/make_qbk.py +++ b/doc/make_qbk.py @@ -4,6 +4,7 @@ # Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. # Copyright (c) 2008-2012 Bruno Lalande, Paris, France. # Copyright (c) 2009-2012 Mateusz Loskot (mateusz@loskot.net), London, UK +# Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland # # Use, modification and distribution is subject to the Boost Software License, # Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at @@ -123,7 +124,8 @@ def cs_to_quickbook(section): , "distance::cross_track", "distance::cross_track_point_box" , "distance::projected_point" , "within::winding", "within::franklin", "within::crossings_multiply" - , "area::surveyor", "area::huiller" + , "area::surveyor", "area::spherical" + #, "area::geographic" , "buffer::point_circle", "buffer::point_square" , "buffer::join_round", "buffer::join_miter" , "buffer::end_round", "buffer::end_flat" diff --git a/doc/quickref.xml b/doc/quickref.xml index 5b095811af..755fdd255c 100644 --- a/doc/quickref.xml +++ b/doc/quickref.xml @@ -10,10 +10,10 @@ Copyright (c) 2009-2015 Bruno Lalande, Paris, France. Copyright (c) 2013-2015 Adam Wulkiewicz, Lodz, Poland. - This file was modified by Oracle on 2014, 2015. - Modifications copyright (c) 2014-2015, Oracle and/or its affiliates. - + This file was modified by Oracle on 2014, 2015, 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 Use, modification and distribution is subject to the Boost Software License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at @@ -521,7 +521,8 @@ Area strategy::area::surveyor - + strategy::area::spherical + diff --git a/doc/reference.qbk b/doc/reference.qbk index f730601816..e84115ece3 100644 --- a/doc/reference.qbk +++ b/doc/reference.qbk @@ -99,7 +99,9 @@ [include generated/crosses.qbk] [endsect] +[section:difference difference] [include generated/difference.qbk] +[endsect] [section:disjoint disjoint] [include generated/disjoint.qbk] @@ -125,9 +127,9 @@ [include generated/for_each.qbk] [endsect] -[/section:intersection intersection] +[section:intersection intersection] [include generated/intersection.qbk] -[/endsect] +[endsect] [section:intersects intersects] [include generated/intersects.qbk] @@ -135,7 +137,9 @@ [include generated/is_empty.qbk] +[section:is_simple is_simple] [include generated/is_simple.qbk] +[endsect] [section:is_valid is_valid] [include generated/is_valid.qbk] @@ -176,9 +180,9 @@ [include generated/simplify.qbk] [endsect] -[/section:sym_difference sym_difference] +[section:sym_difference sym_difference] [include generated/sym_difference.qbk] -[/endsect] +[endsect] [section:touches touches] [include generated/touches.qbk] @@ -188,9 +192,9 @@ [include generated/transform.qbk] [endsect] -[/section:union union] +[section:union_ union_] [include generated/union.qbk] -[/endsect] +[endsect] [include generated/unique.qbk] @@ -333,7 +337,8 @@ [include generated/distance_cross_track.qbk] [include generated/distance_cross_track_point_box.qbk] [include generated/area_surveyor.qbk] -[/include generated/area_huiller.qbk] +[include generated/area_spherical.qbk] +[/include generated/area_geographic.qbk] [include generated/buffer_join_round.qbk] [include generated/buffer_join_miter.qbk] [include generated/buffer_end_round.qbk] diff --git a/doc/release_notes.qbk b/doc/release_notes.qbk index 95d580be19..a640ab0e98 100644 --- a/doc/release_notes.qbk +++ b/doc/release_notes.qbk @@ -29,6 +29,8 @@ [*Breaking changes] * ublas_transformer is renamed to matrix_transformer +* explicit modifier is added to constructors of rtree index::dynamic_* parameters +* strategy::area::huiller replaced by strategy::area::spherical [*Solved issues] diff --git a/doc/src/examples/algorithms/append.cpp b/doc/src/examples/algorithms/append.cpp index f5a3085b5b..fd90a6a04d 100644 --- a/doc/src/examples/algorithms/append.cpp +++ b/doc/src/examples/algorithms/append.cpp @@ -12,12 +12,12 @@ #include -#include - #include #include #include +#include /*< At the end to avoid conflicts with Boost.QVM >*/ + BOOST_GEOMETRY_REGISTER_BOOST_TUPLE_CS(cs::cartesian) int main() diff --git a/doc/src/examples/core/tag.cpp b/doc/src/examples/core/tag.cpp index bb5b765ec8..5e4c706879 100644 --- a/doc/src/examples/core/tag.cpp +++ b/doc/src/examples/core/tag.cpp @@ -12,13 +12,13 @@ #include -#include - #include #include #include #include +#include /*< At the end to avoid conflicts with Boost.QVM >*/ + BOOST_GEOMETRY_REGISTER_BOOST_TUPLE_CS(cs::cartesian) template struct dispatch {}; diff --git a/doc/src/examples/geometries/adapted/boost_range/sliced.cpp b/doc/src/examples/geometries/adapted/boost_range/sliced.cpp index bdfd74c629..45cab8195a 100644 --- a/doc/src/examples/geometries/adapted/boost_range/sliced.cpp +++ b/doc/src/examples/geometries/adapted/boost_range/sliced.cpp @@ -12,13 +12,13 @@ #include -#include - #include #include #include #include +#include /*< At the end to avoid conflicts with Boost.QVM >*/ + int main() { diff --git a/doc/src/examples/geometries/adapted/boost_range/strided.cpp b/doc/src/examples/geometries/adapted/boost_range/strided.cpp index 8b9dc68657..1ed9b5ecc5 100644 --- a/doc/src/examples/geometries/adapted/boost_range/strided.cpp +++ b/doc/src/examples/geometries/adapted/boost_range/strided.cpp @@ -12,13 +12,13 @@ #include -#include - #include #include #include #include +#include /*< At the end to avoid conflicts with Boost.QVM >*/ + int main() { diff --git a/doc/src/examples/quick_start.cpp b/doc/src/examples/quick_start.cpp index f77c8909b7..41d3beb32a 100644 --- a/doc/src/examples/quick_start.cpp +++ b/doc/src/examples/quick_start.cpp @@ -19,6 +19,8 @@ //#pragma warning( disable : 4244 ) #endif // defined(_MSC_VER) +#include + //[quickstart_include #include diff --git a/include/boost/geometry/algorithms/area.hpp b/include/boost/geometry/algorithms/area.hpp index 4751d4e742..18aea24036 100644 --- a/include/boost/geometry/algorithms/area.hpp +++ b/include/boost/geometry/algorithms/area.hpp @@ -4,6 +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 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. @@ -303,7 +307,8 @@ inline typename default_area_result::type area(Geometry const& geometr [heading Available Strategies] \* [link geometry.reference.strategies.strategy_area_surveyor Surveyor (cartesian)] -\* [link geometry.reference.strategies.strategy_area_huiller Huiller (spherical)] +\* [link geometry.reference.strategies.strategy_area_spherical Spherical] +[/link geometry.reference.strategies.strategy_area_geographic Geographic] } */ template diff --git a/include/boost/geometry/algorithms/detail/intersection/interface.hpp b/include/boost/geometry/algorithms/detail/intersection/interface.hpp index ddf49ad031..0efc9731b5 100644 --- a/include/boost/geometry/algorithms/detail/intersection/interface.hpp +++ b/include/boost/geometry/algorithms/detail/intersection/interface.hpp @@ -15,12 +15,14 @@ #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_INTERSECTION_INTERFACE_HPP -// TODO: those headers probably may be removed -#include -#include +#include +#include +#include #include #include +#include +#include namespace boost { namespace geometry @@ -98,19 +100,74 @@ struct intersection } // namespace dispatch #endif // DOXYGEN_NO_DISPATCH - + +namespace resolve_strategy { + +struct intersection +{ + template + < + typename Geometry1, + typename Geometry2, + typename RobustPolicy, + typename GeometryOut, + typename Strategy + > + static inline bool apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + RobustPolicy const& robust_policy, + GeometryOut & geometry_out, + Strategy const& strategy) + { + return dispatch::intersection + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, robust_policy, geometry_out, + strategy); + } + + template + < + typename Geometry1, + typename Geometry2, + typename RobustPolicy, + typename GeometryOut + > + static inline bool apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + RobustPolicy const& robust_policy, + GeometryOut & geometry_out, + default_strategy) + { + typedef typename strategy::relate::services::default_strategy + < + Geometry1, Geometry2 + >::type strategy_type; + + return dispatch::intersection + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, robust_policy, geometry_out, + strategy_type()); + } +}; + +} // resolve_strategy + + namespace resolve_variant { template struct intersection { - template - static inline bool - apply( - const Geometry1& geometry1, - const Geometry2& geometry2, - GeometryOut& geometry_out) + template + static inline bool apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + GeometryOut& geometry_out, + Strategy const& strategy) { concepts::check(); concepts::check(); @@ -125,17 +182,11 @@ struct intersection = geometry::get_rescale_policy(geometry1, geometry2); - typedef typename strategy::relate::services::default_strategy - < - Geometry1, Geometry2 - >::type strategy_type; - - return dispatch::intersection - < - Geometry1, - Geometry2 - >::apply(geometry1, geometry2, robust_policy, geometry_out, - strategy_type()); + return resolve_strategy::intersection::apply(geometry1, + geometry2, + robust_policy, + geometry_out, + strategy); } }; @@ -143,40 +194,43 @@ struct intersection template struct intersection, Geometry2> { - template + template struct visitor: static_visitor { Geometry2 const& m_geometry2; GeometryOut& m_geometry_out; + Strategy const& m_strategy; visitor(Geometry2 const& geometry2, - GeometryOut& geometry_out) + GeometryOut& geometry_out, + Strategy const& strategy) : m_geometry2(geometry2) , m_geometry_out(geometry_out) + , m_strategy(strategy) {} template - result_type operator()(Geometry1 const& geometry1) const + bool operator()(Geometry1 const& geometry1) const { return intersection - < - Geometry1, - Geometry2 - >::template apply - < - GeometryOut - > - (geometry1, m_geometry2, m_geometry_out); + < + Geometry1, + Geometry2 + >::apply(geometry1, m_geometry2, m_geometry_out, m_strategy); } }; - template + template static inline bool apply(variant const& geometry1, Geometry2 const& geometry2, - GeometryOut& geometry_out) + GeometryOut& geometry_out, + Strategy const& strategy) { - return boost::apply_visitor(visitor(geometry2, geometry_out), geometry1); + return boost::apply_visitor(visitor(geometry2, + geometry_out, + strategy), + geometry1); } }; @@ -184,40 +238,43 @@ struct intersection, Geometry2> template struct intersection > { - template + template struct visitor: static_visitor { Geometry1 const& m_geometry1; GeometryOut& m_geometry_out; + Strategy const& m_strategy; visitor(Geometry1 const& geometry1, - GeometryOut& geometry_out) + GeometryOut& geometry_out, + Strategy const& strategy) : m_geometry1(geometry1) , m_geometry_out(geometry_out) + , m_strategy(strategy) {} template - result_type operator()(Geometry2 const& geometry2) const + bool operator()(Geometry2 const& geometry2) const { return intersection - < - Geometry1, - Geometry2 - >::template apply - < - GeometryOut - > - (m_geometry1, geometry2, m_geometry_out); + < + Geometry1, + Geometry2 + >::apply(m_geometry1, geometry2, m_geometry_out, m_strategy); } }; - template + template static inline bool apply(Geometry1 const& geometry1, - const variant& geometry2, - GeometryOut& geometry_out) + variant const& geometry2, + GeometryOut& geometry_out, + Strategy const& strategy) { - return boost::apply_visitor(visitor(geometry1, geometry_out), geometry2); + return boost::apply_visitor(visitor(geometry1, + geometry_out, + strategy), + geometry2); } }; @@ -225,38 +282,39 @@ struct intersection > template struct intersection, variant > { - template + template struct visitor: static_visitor { GeometryOut& m_geometry_out; + Strategy const& m_strategy; - visitor(GeometryOut& geometry_out) + visitor(GeometryOut& geometry_out, Strategy const& strategy) : m_geometry_out(geometry_out) + , m_strategy(strategy) {} template - result_type operator()(Geometry1 const& geometry1, - Geometry2 const& geometry2) const + bool operator()(Geometry1 const& geometry1, + Geometry2 const& geometry2) const { return intersection - < - Geometry1, - Geometry2 - >::template apply - < - GeometryOut - > - (geometry1, geometry2, m_geometry_out); + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, m_geometry_out, m_strategy); } }; - template + template static inline bool - apply(const variant& geometry1, - const variant& geometry2, - GeometryOut& geometry_out) + apply(variant const& geometry1, + variant const& geometry2, + GeometryOut& geometry_out, + Strategy const& strategy) { - return boost::apply_visitor(visitor(geometry_out), geometry1, geometry2); + return boost::apply_visitor(visitor(geometry_out, + strategy), + geometry1, geometry2); } }; @@ -271,32 +329,66 @@ struct intersection, variant inline bool intersection(Geometry1 const& geometry1, - Geometry2 const& geometry2, - GeometryOut& geometry_out) + Geometry2 const& geometry2, + GeometryOut& geometry_out, + Strategy const& strategy) { return resolve_variant::intersection < - Geometry1, - Geometry2 - >::template apply + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, geometry_out, strategy); +} + + +/*! +\brief \brief_calc2{intersection} +\ingroup intersection +\details \details_calc2{intersection, spatial set theoretic intersection}. +\tparam Geometry1 \tparam_geometry +\tparam Geometry2 \tparam_geometry +\tparam GeometryOut Collection of geometries (e.g. std::vector, std::deque, boost::geometry::multi*) of which + the value_type fulfills a \p_l_or_c concept, or it is the output geometry (e.g. for a box) +\param geometry1 \param_geometry +\param geometry2 \param_geometry +\param geometry_out The output geometry, either a multi_point, multi_polygon, + multi_linestring, or a box (for intersection of two boxes) + +\qbk{[include reference/algorithms/intersection.qbk]} +*/ +template +< + typename Geometry1, + typename Geometry2, + typename GeometryOut +> +inline bool intersection(Geometry1 const& geometry1, + Geometry2 const& geometry2, + GeometryOut& geometry_out) +{ + return resolve_variant::intersection < - GeometryOut - > - (geometry1, geometry2, geometry_out); + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, geometry_out, default_strategy()); } diff --git a/include/boost/geometry/algorithms/detail/is_simple/always_simple.hpp b/include/boost/geometry/algorithms/detail/is_simple/always_simple.hpp index 91e2ef76bd..5cec5e1924 100644 --- a/include/boost/geometry/algorithms/detail/is_simple/always_simple.hpp +++ b/include/boost/geometry/algorithms/detail/is_simple/always_simple.hpp @@ -1,8 +1,9 @@ // 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 Menelaos Karavelas, on behalf of Oracle +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Licensed under the Boost Software License version 1.0. // http://www.boost.org/users/license.html @@ -27,7 +28,8 @@ namespace detail { namespace is_simple template struct always_simple { - static inline bool apply(Geometry const&) + template + static inline bool apply(Geometry const&, Strategy const&) { return true; } diff --git a/include/boost/geometry/algorithms/detail/is_simple/areal.hpp b/include/boost/geometry/algorithms/detail/is_simple/areal.hpp index a2322e4831..d4d6db9bce 100644 --- a/include/boost/geometry/algorithms/detail/is_simple/areal.hpp +++ b/include/boost/geometry/algorithms/detail/is_simple/areal.hpp @@ -1,8 +1,9 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2014-2015, Oracle and/or its affiliates. +// 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 // Licensed under the Boost Software License version 1.0. // http://www.boost.org/users/license.html @@ -37,6 +38,12 @@ namespace detail { namespace is_simple template struct is_simple_ring { + template + static inline bool apply(Ring const& ring, Strategy const&) + { + return apply(ring); + } + static inline bool apply(Ring const& ring) { simplicity_failure_policy policy; @@ -69,6 +76,12 @@ class is_simple_polygon } public: + template + static inline bool apply(Polygon const& polygon, Strategy const&) + { + return apply(polygon); + } + static inline bool apply(Polygon const& polygon) { return @@ -119,7 +132,8 @@ struct is_simple template struct is_simple { - static inline bool apply(MultiPolygon const& multipolygon) + template + static inline bool apply(MultiPolygon const& multipolygon, Strategy const&) { return detail::check_iterator_range diff --git a/include/boost/geometry/algorithms/detail/is_simple/interface.hpp b/include/boost/geometry/algorithms/detail/is_simple/interface.hpp index 6d425232b0..af0127dc75 100644 --- a/include/boost/geometry/algorithms/detail/is_simple/interface.hpp +++ b/include/boost/geometry/algorithms/detail/is_simple/interface.hpp @@ -1,8 +1,9 @@ // 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 Menelaos Karavelas, on behalf of Oracle +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Licensed under the Boost Software License version 1.0. // http://www.boost.org/users/license.html @@ -17,46 +18,106 @@ #include #include +#include +#include +#include namespace boost { namespace geometry { +namespace resolve_strategy +{ + +struct is_simple +{ + template + static inline bool apply(Geometry const& geometry, + Strategy const& strategy) + { + return dispatch::is_simple::apply(geometry, strategy); + } + + template + static inline bool apply(Geometry const& geometry, + default_strategy) + { + // NOTE: Currently the strategy is only used for Linear geometries + typedef typename strategy::intersection::services::default_strategy + < + typename cs_tag::type + >::type strategy_type; + + return dispatch::is_simple::apply(geometry, strategy_type()); + } +}; -namespace resolve_variant { +} // namespace resolve_strategy + +namespace resolve_variant +{ template struct is_simple { - static inline bool apply(Geometry const& geometry) + template + static inline bool apply(Geometry const& geometry, Strategy const& strategy) { concepts::check(); - return dispatch::is_simple::apply(geometry); + + return resolve_strategy::is_simple::apply(geometry, strategy); } }; template struct is_simple > { + template struct visitor : boost::static_visitor { + Strategy const& m_strategy; + + visitor(Strategy const& strategy) + : m_strategy(strategy) + {} + template bool operator()(Geometry const& geometry) const { - return is_simple::apply(geometry); + return is_simple::apply(geometry, m_strategy); } }; + template static inline bool - apply(boost::variant const& geometry) + apply(boost::variant const& geometry, + Strategy const& strategy) { - return boost::apply_visitor(visitor(), geometry); + return boost::apply_visitor(visitor(strategy), geometry); } }; } // namespace resolve_variant +/*! +\brief \brief_check{is simple} +\ingroup is_simple +\tparam Geometry \tparam_geometry +\tparam Strategy \tparam_strategy{Is_simple} +\param geometry \param_geometry +\param strategy \param_strategy{is_simple} +\return \return_check{is simple} + +\qbk{distinguish,with strategy} +\qbk{[include reference/algorithms/is_simple.qbk]} +*/ +template +inline bool is_simple(Geometry const& geometry, Strategy const& strategy) +{ + return resolve_variant::is_simple::apply(geometry, strategy); +} + /*! \brief \brief_check{is simple} @@ -70,7 +131,7 @@ struct is_simple > template inline bool is_simple(Geometry const& geometry) { - return resolve_variant::is_simple::apply(geometry); + return resolve_variant::is_simple::apply(geometry, default_strategy()); } diff --git a/include/boost/geometry/algorithms/detail/is_simple/linear.hpp b/include/boost/geometry/algorithms/detail/is_simple/linear.hpp index 6c469f07f7..52b9d9d1c8 100644 --- a/include/boost/geometry/algorithms/detail/is_simple/linear.hpp +++ b/include/boost/geometry/algorithms/detail/is_simple/linear.hpp @@ -189,8 +189,8 @@ class is_acceptable_turn }; -template -inline bool has_self_intersections(Linear const& linear) +template +inline bool has_self_intersections(Linear const& linear, Strategy const& strategy) { typedef typename point_type::type point_type; @@ -217,16 +217,11 @@ inline bool has_self_intersections(Linear const& linear) is_acceptable_turn > interrupt_policy(predicate); - typedef typename strategy::intersection::services::default_strategy - < - typename cs_tag::type - >::type strategy_type; - detail::self_get_turn_points::get_turns < turn_policy >::apply(linear, - strategy_type(), + strategy, detail::no_rescale_policy(), turns, interrupt_policy); @@ -252,8 +247,19 @@ struct is_simple_linestring && ! detail::is_valid::has_spikes < Linestring, closed - >::apply(linestring, policy) - && ! (CheckSelfIntersections && has_self_intersections(linestring)); + >::apply(linestring, policy); + } +}; + +template +struct is_simple_linestring +{ + template + static inline bool apply(Linestring const& linestring, + Strategy const& strategy) + { + return is_simple_linestring::apply(linestring) + && ! has_self_intersections(linestring, strategy); } }; @@ -261,7 +267,9 @@ struct is_simple_linestring template struct is_simple_multilinestring { - static inline bool apply(MultiLinestring const& multilinestring) + template + static inline bool apply(MultiLinestring const& multilinestring, + Strategy const& strategy) { // check each of the linestrings for simplicity // but do not compute self-intersections yet; these will be @@ -281,7 +289,7 @@ struct is_simple_multilinestring return false; } - return ! has_self_intersections(multilinestring); + return ! has_self_intersections(multilinestring, strategy); } }; diff --git a/include/boost/geometry/algorithms/detail/is_simple/multipoint.hpp b/include/boost/geometry/algorithms/detail/is_simple/multipoint.hpp index f9f43d1cdb..61f0bc9130 100644 --- a/include/boost/geometry/algorithms/detail/is_simple/multipoint.hpp +++ b/include/boost/geometry/algorithms/detail/is_simple/multipoint.hpp @@ -1,8 +1,9 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2014-2015, Oracle and/or its affiliates. +// 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 // Licensed under the Boost Software License version 1.0. // http://www.boost.org/users/license.html @@ -38,7 +39,8 @@ namespace detail { namespace is_simple template struct is_simple_multipoint { - static inline bool apply(MultiPoint const& multipoint) + template + static inline bool apply(MultiPoint const& multipoint, Strategy const&) { if (boost::empty(multipoint)) { diff --git a/include/boost/geometry/algorithms/detail/is_valid/box.hpp b/include/boost/geometry/algorithms/detail/is_valid/box.hpp index 863ce625fe..69a4d4e78e 100644 --- a/include/boost/geometry/algorithms/detail/is_valid/box.hpp +++ b/include/boost/geometry/algorithms/detail/is_valid/box.hpp @@ -1,6 +1,6 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2014-2015, Oracle and/or its affiliates. +// 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 @@ -71,8 +71,8 @@ struct has_valid_corners template struct is_valid_box { - template - static inline bool apply(Box const& box, VisitPolicy& visitor) + template + static inline bool apply(Box const& box, VisitPolicy& visitor, Strategy const&) { return ! has_invalid_coordinate::apply(box, visitor) diff --git a/include/boost/geometry/algorithms/detail/is_valid/has_valid_self_turns.hpp b/include/boost/geometry/algorithms/detail/is_valid/has_valid_self_turns.hpp index 9f278ac041..b91dc6a697 100644 --- a/include/boost/geometry/algorithms/detail/is_valid/has_valid_self_turns.hpp +++ b/include/boost/geometry/algorithms/detail/is_valid/has_valid_self_turns.hpp @@ -48,11 +48,6 @@ class has_valid_self_turns private: typedef typename point_type::type point_type; - typedef typename strategy::intersection::services::default_strategy - < - typename cs_tag::type - >::type intersection_strategy_type; - typedef typename geometry::rescale_policy_type < point_type @@ -75,15 +70,14 @@ class has_valid_self_turns > turn_type; // returns true if all turns are valid - template + template static inline bool apply(Geometry const& geometry, Turns& turns, - VisitPolicy& visitor) + VisitPolicy& visitor, + Strategy const& strategy) { boost::ignore_unused(visitor); - intersection_strategy_type intersection_strategy; - rescale_policy_type robust_policy = geometry::get_rescale_policy(geometry); @@ -93,7 +87,7 @@ class has_valid_self_turns > interrupt_policy; geometry::self_turns(geometry, - intersection_strategy, + strategy, robust_policy, turns, interrupt_policy); @@ -110,11 +104,11 @@ class has_valid_self_turns } // returns true if all turns are valid - template - static inline bool apply(Geometry const& geometry, VisitPolicy& visitor) + template + static inline bool apply(Geometry const& geometry, VisitPolicy& visitor, Strategy const& strategy) { std::vector turns; - return apply(geometry, turns, visitor); + return apply(geometry, turns, visitor, strategy); } }; diff --git a/include/boost/geometry/algorithms/detail/is_valid/interface.hpp b/include/boost/geometry/algorithms/detail/is_valid/interface.hpp index 5a04a92824..ee013377c4 100644 --- a/include/boost/geometry/algorithms/detail/is_valid/interface.hpp +++ b/include/boost/geometry/algorithms/detail/is_valid/interface.hpp @@ -1,8 +1,9 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2014-2015, Oracle and/or its affiliates. +// 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 // Licensed under the Boost Software License version 1.0. // http://www.boost.org/users/license.html @@ -23,48 +24,88 @@ #include #include #include +#include +#include namespace boost { namespace geometry { + +namespace resolve_strategy +{ + +struct is_valid +{ + template + static inline bool apply(Geometry const& geometry, + VisitPolicy& visitor, + Strategy const& strategy) + { + return dispatch::is_valid::apply(geometry, visitor, strategy); + } + + template + static inline bool apply(Geometry const& geometry, + VisitPolicy& visitor, + default_strategy) + { + // NOTE: Currently the strategy is only used for Areal geometries + typedef typename strategy::intersection::services::default_strategy + < + typename cs_tag::type + >::type strategy_type; + + return dispatch::is_valid::apply(geometry, visitor, strategy_type()); + } +}; +} // namespace resolve_strategy -namespace resolve_variant { +namespace resolve_variant +{ template struct is_valid { - template - static inline bool apply(Geometry const& geometry, VisitPolicy& visitor) + template + static inline bool apply(Geometry const& geometry, + VisitPolicy& visitor, + Strategy const& strategy) { concepts::check(); - return dispatch::is_valid::apply(geometry, visitor); + + return resolve_strategy::is_valid::apply(geometry, visitor, strategy); } }; template struct is_valid > { - template + template struct visitor : boost::static_visitor { - visitor(VisitPolicy& policy) : m_policy(policy) {} + visitor(VisitPolicy& policy, Strategy const& strategy) + : m_policy(policy) + , m_strategy(strategy) + {} template bool operator()(Geometry const& geometry) const { - return is_valid::apply(geometry, m_policy); + return is_valid::apply(geometry, m_policy, m_strategy); } VisitPolicy& m_policy; + Strategy const& m_strategy; }; - template + template static inline bool apply(boost::variant const& geometry, - VisitPolicy& policy_visitor) + VisitPolicy& policy_visitor, + Strategy const& strategy) { - return boost::apply_visitor(visitor(policy_visitor), + return boost::apply_visitor(visitor(policy_visitor, strategy), geometry); } }; @@ -73,13 +114,38 @@ struct is_valid > // Undocumented for now -template -inline bool is_valid(Geometry const& geometry, VisitPolicy& visitor) +template +inline bool is_valid(Geometry const& geometry, + VisitPolicy& visitor, + Strategy const& strategy) { - return resolve_variant::is_valid::apply(geometry, visitor); + return resolve_variant::is_valid::apply(geometry, visitor, strategy); } +/*! +\brief \brief_check{is valid (in the OGC sense)} +\ingroup is_valid +\tparam Geometry \tparam_geometry +\tparam Strategy \tparam_strategy{Is_valid} +\param geometry \param_geometry +\param strategy \param_strategy{is_valid} +\return \return_check{is valid (in the OGC sense); +furthermore, the following geometries are considered valid: +multi-geometries with no elements, +linear geometries containing spikes, +areal geometries with duplicate (consecutive) points} + +\qbk{distinguish,with strategy} +\qbk{[include reference/algorithms/is_valid.qbk]} +*/ +template +inline bool is_valid(Geometry const& geometry, Strategy const& strategy) +{ + is_valid_default_policy<> visitor; + return resolve_variant::is_valid::apply(geometry, visitor, strategy); +} + /*! \brief \brief_check{is valid (in the OGC sense)} \ingroup is_valid @@ -96,11 +162,37 @@ inline bool is_valid(Geometry const& geometry, VisitPolicy& visitor) template inline bool is_valid(Geometry const& geometry) { - is_valid_default_policy<> policy_visitor; - return geometry::is_valid(geometry, policy_visitor); + return is_valid(geometry, default_strategy()); } +/*! +\brief \brief_check{is valid (in the OGC sense)} +\ingroup is_valid +\tparam Geometry \tparam_geometry +\tparam Strategy \tparam_strategy{Is_valid} +\param geometry \param_geometry +\param failure An enumeration value indicating that the geometry is + valid or not, and if not valid indicating the reason why +\param strategy \param_strategy{is_valid} +\return \return_check{is valid (in the OGC sense); + furthermore, the following geometries are considered valid: + multi-geometries with no elements, + linear geometries containing spikes, + areal geometries with duplicate (consecutive) points} + +\qbk{distinguish,with failure value and strategy} +\qbk{[include reference/algorithms/is_valid_with_failure.qbk]} +*/ +template +inline bool is_valid(Geometry const& geometry, validity_failure_type& failure, Strategy const& strategy) +{ + failure_type_policy<> visitor; + bool result = resolve_variant::is_valid::apply(geometry, visitor, strategy); + failure = visitor.failure(); + return result; +} + /*! \brief \brief_check{is valid (in the OGC sense)} \ingroup is_valid @@ -120,10 +212,7 @@ inline bool is_valid(Geometry const& geometry) template inline bool is_valid(Geometry const& geometry, validity_failure_type& failure) { - failure_type_policy<> policy_visitor; - bool result = geometry::is_valid(geometry, policy_visitor); - failure = policy_visitor.failure(); - return result; + return is_valid(geometry, failure, default_strategy()); } @@ -131,28 +220,52 @@ inline bool is_valid(Geometry const& geometry, validity_failure_type& failure) \brief \brief_check{is valid (in the OGC sense)} \ingroup is_valid \tparam Geometry \tparam_geometry +\tparam Strategy \tparam_strategy{Is_valid} \param geometry \param_geometry \param message A string containing a message stating if the geometry is valid or not, and if not valid a reason why +\param strategy \param_strategy{is_valid} \return \return_check{is valid (in the OGC sense); furthermore, the following geometries are considered valid: multi-geometries with no elements, linear geometries containing spikes, areal geometries with duplicate (consecutive) points} -\qbk{distinguish,with message} +\qbk{distinguish,with message and strategy} \qbk{[include reference/algorithms/is_valid_with_message.qbk]} */ -template -inline bool is_valid(Geometry const& geometry, std::string& message) +template +inline bool is_valid(Geometry const& geometry, std::string& message, Strategy const& strategy) { std::ostringstream stream; - failing_reason_policy<> policy_visitor(stream); - bool result = geometry::is_valid(geometry, policy_visitor); + failing_reason_policy<> visitor(stream); + bool result = resolve_variant::is_valid::apply(geometry, visitor, strategy); message = stream.str(); return result; } +/*! +\brief \brief_check{is valid (in the OGC sense)} +\ingroup is_valid +\tparam Geometry \tparam_geometry +\param geometry \param_geometry +\param message A string containing a message stating if the geometry + is valid or not, and if not valid a reason why +\return \return_check{is valid (in the OGC sense); + furthermore, the following geometries are considered valid: + multi-geometries with no elements, + linear geometries containing spikes, + areal geometries with duplicate (consecutive) points} + +\qbk{distinguish,with message} +\qbk{[include reference/algorithms/is_valid_with_message.qbk]} +*/ +template +inline bool is_valid(Geometry const& geometry, std::string& message) +{ + return is_valid(geometry, message, default_strategy()); +} + }} // namespace boost::geometry diff --git a/include/boost/geometry/algorithms/detail/is_valid/linear.hpp b/include/boost/geometry/algorithms/detail/is_valid/linear.hpp index a49e077237..6bc6b86cf8 100644 --- a/include/boost/geometry/algorithms/detail/is_valid/linear.hpp +++ b/include/boost/geometry/algorithms/detail/is_valid/linear.hpp @@ -1,6 +1,6 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2014-2015, Oracle and/or its affiliates. +// 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 @@ -77,6 +77,14 @@ struct is_valid_linestring } return ! has_spikes::apply(linestring, visitor); } + + template + static inline bool apply(Linestring const& linestring, + VisitPolicy& visitor, + Strategy const&) + { + return apply(linestring, visitor); + } }; @@ -142,9 +150,10 @@ class is_valid }; public: - template + template static inline bool apply(MultiLinestring const& multilinestring, - VisitPolicy& visitor) + VisitPolicy& visitor, + Strategy const&) { if (BOOST_GEOMETRY_CONDITION( AllowEmptyMultiGeometries && boost::empty(multilinestring))) diff --git a/include/boost/geometry/algorithms/detail/is_valid/multipolygon.hpp b/include/boost/geometry/algorithms/detail/is_valid/multipolygon.hpp index eafe07841a..84dacc57f1 100644 --- a/include/boost/geometry/algorithms/detail/is_valid/multipolygon.hpp +++ b/include/boost/geometry/algorithms/detail/is_valid/multipolygon.hpp @@ -237,23 +237,28 @@ class is_valid_multipolygon } - template + template struct per_polygon { - per_polygon(VisitPolicy& policy) : m_policy(policy) {} + per_polygon(VisitPolicy& policy, Strategy const& strategy) + : m_policy(policy) + , m_strategy(strategy) + {} template inline bool apply(Polygon const& polygon) const { - return base::apply(polygon, m_policy); + return base::apply(polygon, m_policy, m_strategy); } VisitPolicy& m_policy; + Strategy const& m_strategy; }; public: - template + template static inline bool apply(MultiPolygon const& multipolygon, - VisitPolicy& visitor) + VisitPolicy& visitor, + Strategy const& strategy) { typedef debug_validity_phase debug_phase; @@ -268,11 +273,11 @@ class is_valid_multipolygon if (! detail::check_iterator_range < - per_polygon, + per_polygon, false // do not check for empty multipolygon (done above) >::apply(boost::begin(multipolygon), boost::end(multipolygon), - per_polygon(visitor))) + per_polygon(visitor, strategy))) { return false; } @@ -285,7 +290,7 @@ class is_valid_multipolygon std::deque turns; bool has_invalid_turns = - ! has_valid_turns::apply(multipolygon, turns, visitor); + ! has_valid_turns::apply(multipolygon, turns, visitor, strategy); debug_print_turns(turns.begin(), turns.end()); if (has_invalid_turns) diff --git a/include/boost/geometry/algorithms/detail/is_valid/pointlike.hpp b/include/boost/geometry/algorithms/detail/is_valid/pointlike.hpp index 51035f7a73..f77f7a35eb 100644 --- a/include/boost/geometry/algorithms/detail/is_valid/pointlike.hpp +++ b/include/boost/geometry/algorithms/detail/is_valid/pointlike.hpp @@ -1,6 +1,6 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2014-2015, Oracle and/or its affiliates. +// 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 @@ -36,8 +36,8 @@ namespace dispatch template struct is_valid { - template - static inline bool apply(Point const& point, VisitPolicy& visitor) + template + static inline bool apply(Point const& point, VisitPolicy& visitor, Strategy const&) { boost::ignore_unused(visitor); return ! detail::is_valid::has_invalid_coordinate @@ -56,9 +56,10 @@ struct is_valid template struct is_valid { - template + template static inline bool apply(MultiPoint const& multipoint, - VisitPolicy& visitor) + VisitPolicy& visitor, + Strategy const&) { boost::ignore_unused(multipoint, visitor); diff --git a/include/boost/geometry/algorithms/detail/is_valid/polygon.hpp b/include/boost/geometry/algorithms/detail/is_valid/polygon.hpp index 182571c567..f7e22fb8d2 100644 --- a/include/boost/geometry/algorithms/detail/is_valid/polygon.hpp +++ b/include/boost/geometry/algorithms/detail/is_valid/polygon.hpp @@ -1,8 +1,9 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2014-2015, Oracle and/or its affiliates. +// 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 // Licensed under the Boost Software License version 1.0. // http://www.boost.org/users/license.html @@ -74,10 +75,13 @@ class is_valid_polygon protected: typedef debug_validity_phase debug_phase; - template + template struct per_ring { - per_ring(VisitPolicy& policy) : m_policy(policy) {} + per_ring(VisitPolicy& policy, Strategy const& strategy) + : m_policy(policy) + , m_strategy(strategy) + {} template inline bool apply(Ring const& ring) const @@ -85,30 +89,34 @@ class is_valid_polygon return detail::is_valid::is_valid_ring < Ring, false, true - >::apply(ring, m_policy); + >::apply(ring, m_policy, m_strategy); } VisitPolicy& m_policy; + Strategy const& m_strategy; }; - template + template static bool has_valid_interior_rings(InteriorRings const& interior_rings, - VisitPolicy& visitor) + VisitPolicy& visitor, + Strategy const& strategy) { return detail::check_iterator_range < - per_ring, + per_ring, true // allow for empty interior ring range >::apply(boost::begin(interior_rings), boost::end(interior_rings), - per_ring(visitor)); + per_ring(visitor, strategy)); } struct has_valid_rings { - template - static inline bool apply(Polygon const& polygon, VisitPolicy& visitor) + template + static inline bool apply(Polygon const& polygon, + VisitPolicy& visitor, + Strategy const& strategy) { typedef typename ring_type::type ring_type; @@ -119,7 +127,7 @@ class is_valid_polygon < ring_type, false // do not check self intersections - >::apply(exterior_ring(polygon), visitor)) + >::apply(exterior_ring(polygon), visitor, strategy)) { return false; } @@ -128,7 +136,8 @@ class is_valid_polygon debug_phase::apply(2); return has_valid_interior_rings(geometry::interior_rings(polygon), - visitor); + visitor, + strategy); } }; @@ -344,10 +353,12 @@ class is_valid_polygon }; public: - template - static inline bool apply(Polygon const& polygon, VisitPolicy& visitor) + template + static inline bool apply(Polygon const& polygon, + VisitPolicy& visitor, + Strategy const& strategy) { - if (! has_valid_rings::apply(polygon, visitor)) + if (! has_valid_rings::apply(polygon, visitor, strategy)) { return false; } @@ -364,7 +375,7 @@ class is_valid_polygon std::deque turns; bool has_invalid_turns - = ! has_valid_turns::apply(polygon, turns, visitor); + = ! has_valid_turns::apply(polygon, turns, visitor, strategy); debug_print_turns(turns.begin(), turns.end()); if (has_invalid_turns) diff --git a/include/boost/geometry/algorithms/detail/is_valid/ring.hpp b/include/boost/geometry/algorithms/detail/is_valid/ring.hpp index 925c03a472..9ab68fdc48 100644 --- a/include/boost/geometry/algorithms/detail/is_valid/ring.hpp +++ b/include/boost/geometry/algorithms/detail/is_valid/ring.hpp @@ -1,6 +1,6 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2014-2015, Oracle and/or its affiliates. +// 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 @@ -101,26 +101,21 @@ struct ring_area_predicate template struct is_properly_oriented { - typedef typename point_type::type point_type; - - typedef typename strategy::area::services::default_strategy - < - typename cs_tag::type, - point_type - >::type strategy_type; + template + static inline bool apply(Ring const& ring, VisitPolicy& visitor, + Strategy const& strategy) + { + boost::ignore_unused(visitor); - typedef detail::area::ring_area - < - order_as_direction::value>::value, - geometry::closure::value - > ring_area_type; + typedef typename point_type::type point_type; - typedef typename default_area_result::type area_result_type; + typedef detail::area::ring_area + < + order_as_direction::value>::value, + geometry::closure::value + > ring_area_type; - template - static inline bool apply(Ring const& ring, VisitPolicy& visitor) - { - boost::ignore_unused(visitor); + typedef typename default_area_result::type area_result_type; typename ring_area_predicate < @@ -128,8 +123,11 @@ struct is_properly_oriented >::type predicate; // Check area - area_result_type const zero = area_result_type(); - if (predicate(ring_area_type::apply(ring, strategy_type()), zero)) + area_result_type const zero = 0; + area_result_type const area + = ring_area_type::apply(ring, + strategy.template get_area_strategy()); + if (predicate(area, zero)) { return visitor.template apply(); } @@ -150,8 +148,9 @@ template > struct is_valid_ring { - template - static inline bool apply(Ring const& ring, VisitPolicy& visitor) + template + static inline bool apply(Ring const& ring, VisitPolicy& visitor, + Strategy const& strategy) { // return invalid if any of the following condition holds: // (a) the ring's point coordinates are not invalid (e.g., NaN) @@ -198,8 +197,8 @@ struct is_valid_ring && ! has_duplicates::apply(ring, visitor) && ! has_spikes::apply(ring, visitor) && (! CheckSelfIntersections - || has_valid_self_turns::apply(ring, visitor)) - && is_properly_oriented::apply(ring, visitor); + || has_valid_self_turns::apply(ring, visitor, strategy)) + && is_properly_oriented::apply(ring, visitor, strategy); } }; diff --git a/include/boost/geometry/algorithms/detail/is_valid/segment.hpp b/include/boost/geometry/algorithms/detail/is_valid/segment.hpp index f92f73381f..30cbf7afdb 100644 --- a/include/boost/geometry/algorithms/detail/is_valid/segment.hpp +++ b/include/boost/geometry/algorithms/detail/is_valid/segment.hpp @@ -1,6 +1,6 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2014-2015, Oracle and/or its affiliates. +// 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 @@ -44,8 +44,8 @@ namespace dispatch template struct is_valid { - template - static inline bool apply(Segment const& segment, VisitPolicy& visitor) + template + static inline bool apply(Segment const& segment, VisitPolicy& visitor, Strategy const&) { boost::ignore_unused(visitor); diff --git a/include/boost/geometry/algorithms/detail/overlay/get_turns.hpp b/include/boost/geometry/algorithms/detail/overlay/get_turns.hpp index a2e597d912..4e97a84a37 100644 --- a/include/boost/geometry/algorithms/detail/overlay/get_turns.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/get_turns.hpp @@ -45,7 +45,6 @@ #include -#include #include #include diff --git a/include/boost/geometry/algorithms/difference.hpp b/include/boost/geometry/algorithms/difference.hpp index 2d0afffe83..c11ceca243 100644 --- a/include/boost/geometry/algorithms/difference.hpp +++ b/include/boost/geometry/algorithms/difference.hpp @@ -4,6 +4,7 @@ // This file was modified by Oracle on 2017. // Modifications copyright (c) 2017, Oracle and/or its affiliates. + // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Use, modification and distribution is subject to the Boost Software License, @@ -13,10 +14,16 @@ #ifndef BOOST_GEOMETRY_ALGORITHMS_DIFFERENCE_HPP #define BOOST_GEOMETRY_ALGORITHMS_DIFFERENCE_HPP -#include + +#include +#include +#include #include #include +#include +#include + namespace boost { namespace geometry { @@ -101,18 +108,14 @@ inline OutputIterator difference_insert(Geometry1 const& geometry1, RobustPolicy const& robust_policy, OutputIterator out) { - concepts::check(); - concepts::check(); - concepts::check(); - - typename strategy::relate::services::default_strategy + typedef typename strategy::relate::services::default_strategy < Geometry1, Geometry2 - >::type strategy; + >::type strategy_type; return difference_insert(geometry1, geometry2, - robust_policy, out, strategy); + robust_policy, out, strategy_type()); } @@ -120,6 +123,215 @@ inline OutputIterator difference_insert(Geometry1 const& geometry1, #endif // DOXYGEN_NO_DETAIL +namespace resolve_strategy { + +struct difference +{ + template + < + typename Geometry1, + typename Geometry2, + typename RobustPolicy, + typename Collection, + typename Strategy + > + static inline void apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + RobustPolicy const& robust_policy, + Collection & output_collection, + Strategy const& strategy) + { + typedef typename boost::range_value::type geometry_out; + + detail::difference::difference_insert( + geometry1, geometry2, robust_policy, + range::back_inserter(output_collection), + strategy); + } + + template + < + typename Geometry1, + typename Geometry2, + typename RobustPolicy, + typename Collection + > + static inline void apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + RobustPolicy const& robust_policy, + Collection & output_collection, + default_strategy) + { + typedef typename boost::range_value::type geometry_out; + + detail::difference::difference_insert( + geometry1, geometry2, robust_policy, + range::back_inserter(output_collection)); + } +}; + +} // resolve_strategy + + +namespace resolve_variant +{ + +template +struct difference +{ + template + static inline void apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + typedef typename geometry::rescale_overlay_policy_type + < + Geometry1, + Geometry2 + >::type rescale_policy_type; + + rescale_policy_type robust_policy + = geometry::get_rescale_policy(geometry1, + geometry2); + + resolve_strategy::difference::apply(geometry1, geometry2, + robust_policy, + output_collection, + strategy); + } +}; + + +template +struct difference, Geometry2> +{ + template + struct visitor: static_visitor<> + { + Geometry2 const& m_geometry2; + Collection& m_output_collection; + Strategy const& m_strategy; + + visitor(Geometry2 const& geometry2, + Collection& output_collection, + Strategy const& strategy) + : m_geometry2(geometry2) + , m_output_collection(output_collection) + , m_strategy(strategy) + {} + + template + void operator()(Geometry1 const& geometry1) const + { + difference + < + Geometry1, + Geometry2 + >::apply(geometry1, m_geometry2, m_output_collection, m_strategy); + } + }; + + template + static inline void + apply(variant const& geometry1, + Geometry2 const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + boost::apply_visitor(visitor(geometry2, + output_collection, + strategy), + geometry1); + } +}; + + +template +struct difference > +{ + template + struct visitor: static_visitor<> + { + Geometry1 const& m_geometry1; + Collection& m_output_collection; + Strategy const& m_strategy; + + visitor(Geometry1 const& geometry1, + Collection& output_collection, + Strategy const& strategy) + : m_geometry1(geometry1) + , m_output_collection(output_collection) + , m_strategy(strategy) + {} + + template + void operator()(Geometry2 const& geometry2) const + { + difference + < + Geometry1, + Geometry2 + >::apply(m_geometry1, geometry2, m_output_collection, m_strategy); + } + }; + + template + static inline void + apply(Geometry1 const& geometry1, + variant const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + boost::apply_visitor(visitor(geometry1, + output_collection, + strategy), + geometry2); + } +}; + + +template +struct difference, variant > +{ + template + struct visitor: static_visitor<> + { + Collection& m_output_collection; + Strategy const& m_strategy; + + visitor(Collection& output_collection, Strategy const& strategy) + : m_output_collection(output_collection) + , m_strategy(strategy) + {} + + template + void operator()(Geometry1 const& geometry1, + Geometry2 const& geometry2) const + { + difference + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, m_output_collection, m_strategy); + } + }; + + template + static inline void + apply(variant const& geometry1, + variant const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + boost::apply_visitor(visitor(output_collection, + strategy), + geometry1, geometry2); + } +}; + +} // namespace resolve_variant + /*! \brief_calc2{difference} @@ -128,41 +340,63 @@ inline OutputIterator difference_insert(Geometry1 const& geometry1, \tparam Geometry1 \tparam_geometry \tparam Geometry2 \tparam_geometry \tparam Collection \tparam_output_collection +\tparam Strategy \tparam_strategy{Difference} \param geometry1 \param_geometry \param geometry2 \param_geometry \param output_collection the output collection +\param strategy \param_strategy{difference} +\qbk{distinguish,with strategy} \qbk{[include reference/algorithms/difference.qbk]} */ template < typename Geometry1, typename Geometry2, - typename Collection + typename Collection, + typename Strategy > inline void difference(Geometry1 const& geometry1, Geometry2 const& geometry2, - Collection& output_collection) + Collection& output_collection, + Strategy const& strategy) { - concepts::check(); - concepts::check(); - - typedef typename boost::range_value::type geometry_out; - concepts::check(); - - typedef typename geometry::rescale_overlay_policy_type + resolve_variant::difference < Geometry1, Geometry2 - >::type rescale_policy_type; + >::apply(geometry1, geometry2, output_collection, strategy); +} - rescale_policy_type robust_policy - = geometry::get_rescale_policy(geometry1, - geometry2); - detail::difference::difference_insert( - geometry1, geometry2, robust_policy, - range::back_inserter(output_collection)); +/*! +\brief_calc2{difference} +\ingroup difference +\details \details_calc2{difference, spatial set theoretic difference}. +\tparam Geometry1 \tparam_geometry +\tparam Geometry2 \tparam_geometry +\tparam Collection \tparam_output_collection +\param geometry1 \param_geometry +\param geometry2 \param_geometry +\param output_collection the output collection + +\qbk{[include reference/algorithms/difference.qbk]} +*/ +template +< + typename Geometry1, + typename Geometry2, + typename Collection +> +inline void difference(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Collection& output_collection) +{ + resolve_variant::difference + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, output_collection, default_strategy()); } diff --git a/include/boost/geometry/algorithms/equals.hpp b/include/boost/geometry/algorithms/equals.hpp index 550a4d66cd..1479ea66fa 100644 --- a/include/boost/geometry/algorithms/equals.hpp +++ b/include/boost/geometry/algorithms/equals.hpp @@ -138,24 +138,32 @@ struct segment_segment struct area_check { - template - static inline bool apply(Geometry1 const& geometry1, Geometry2 const& geometry2) + template + static inline bool apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Strategy const& strategy) { return geometry::math::equals( - geometry::area(geometry1), - geometry::area(geometry2)); + geometry::area(geometry1, + strategy.template get_area_strategy()), + geometry::area(geometry2, + strategy.template get_area_strategy())); } }; struct length_check { - template - static inline bool apply(Geometry1 const& geometry1, Geometry2 const& geometry2) + template + static inline bool apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Strategy const& strategy) { return geometry::math::equals( - geometry::length(geometry1), - geometry::length(geometry2)); + geometry::length(geometry1, + strategy.template get_distance_strategy()), + geometry::length(geometry2, + strategy.template get_distance_strategy())); } }; @@ -184,9 +192,11 @@ template struct equals_by_collection { template - static inline bool apply(Geometry1 const& geometry1, Geometry2 const& geometry2, Strategy const&) + static inline bool apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Strategy const& strategy) { - if (! TrivialCheck::apply(geometry1, geometry2)) + if (! TrivialCheck::apply(geometry1, geometry2, strategy)) { return false; } diff --git a/include/boost/geometry/algorithms/sym_difference.hpp b/include/boost/geometry/algorithms/sym_difference.hpp index 5d5bd904b5..725230cd5b 100644 --- a/include/boost/geometry/algorithms/sym_difference.hpp +++ b/include/boost/geometry/algorithms/sym_difference.hpp @@ -15,13 +15,21 @@ #ifndef BOOST_GEOMETRY_ALGORITHMS_SYM_DIFFERENCE_HPP #define BOOST_GEOMETRY_ALGORITHMS_SYM_DIFFERENCE_HPP + #include #include #include +#include +#include +#include + #include #include #include +#include +#include +#include namespace boost { namespace geometry @@ -289,6 +297,216 @@ inline OutputIterator sym_difference_insert(Geometry1 const& geometry1, #endif // DOXYGEN_NO_DETAIL +namespace resolve_strategy { + +struct sym_difference +{ + template + < + typename Geometry1, + typename Geometry2, + typename RobustPolicy, + typename Collection, + typename Strategy + > + static inline void apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + RobustPolicy const& robust_policy, + Collection & output_collection, + Strategy const& strategy) + { + typedef typename boost::range_value::type geometry_out; + + detail::sym_difference::sym_difference_insert( + geometry1, geometry2, robust_policy, + range::back_inserter(output_collection), + strategy); + } + + template + < + typename Geometry1, + typename Geometry2, + typename RobustPolicy, + typename Collection + > + static inline void apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + RobustPolicy const& robust_policy, + Collection & output_collection, + default_strategy) + { + typedef typename boost::range_value::type geometry_out; + + detail::sym_difference::sym_difference_insert( + geometry1, geometry2, robust_policy, + range::back_inserter(output_collection)); + } +}; + +} // resolve_strategy + + +namespace resolve_variant +{ + +template +struct sym_difference +{ + template + static inline void apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + typedef typename geometry::rescale_overlay_policy_type + < + Geometry1, + Geometry2 + >::type rescale_policy_type; + + rescale_policy_type robust_policy + = geometry::get_rescale_policy(geometry1, + geometry2); + + resolve_strategy::sym_difference::apply(geometry1, geometry2, + robust_policy, + output_collection, + strategy); + } +}; + + +template +struct sym_difference, Geometry2> +{ + template + struct visitor: static_visitor<> + { + Geometry2 const& m_geometry2; + Collection& m_output_collection; + Strategy const& m_strategy; + + visitor(Geometry2 const& geometry2, + Collection& output_collection, + Strategy const& strategy) + : m_geometry2(geometry2) + , m_output_collection(output_collection) + , m_strategy(strategy) + {} + + template + void operator()(Geometry1 const& geometry1) const + { + sym_difference + < + Geometry1, + Geometry2 + >::apply(geometry1, m_geometry2, m_output_collection, m_strategy); + } + }; + + template + static inline void + apply(variant const& geometry1, + Geometry2 const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + boost::apply_visitor(visitor(geometry2, + output_collection, + strategy), + geometry1); + } +}; + + +template +struct sym_difference > +{ + template + struct visitor: static_visitor<> + { + Geometry1 const& m_geometry1; + Collection& m_output_collection; + Strategy const& m_strategy; + + visitor(Geometry1 const& geometry1, + Collection& output_collection, + Strategy const& strategy) + : m_geometry1(geometry1) + , m_output_collection(output_collection) + , m_strategy(strategy) + {} + + template + void operator()(Geometry2 const& geometry2) const + { + sym_difference + < + Geometry1, + Geometry2 + >::apply(m_geometry1, geometry2, m_output_collection, m_strategy); + } + }; + + template + static inline void + apply(Geometry1 const& geometry1, + variant const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + boost::apply_visitor(visitor(geometry1, + output_collection, + strategy), + geometry2); + } +}; + + +template +struct sym_difference, variant > +{ + template + struct visitor: static_visitor<> + { + Collection& m_output_collection; + Strategy const& m_strategy; + + visitor(Collection& output_collection, Strategy const& strategy) + : m_output_collection(output_collection) + , m_strategy(strategy) + {} + + template + void operator()(Geometry1 const& geometry1, + Geometry2 const& geometry2) const + { + sym_difference + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, m_output_collection, m_strategy); + } + }; + + template + static inline void + apply(variant const& geometry1, + variant const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + boost::apply_visitor(visitor(output_collection, + strategy), + geometry1, geometry2); + } +}; + +} // namespace resolve_variant + + /*! \brief \brief_calc2{symmetric difference} \ingroup sym_difference @@ -297,39 +515,64 @@ inline OutputIterator sym_difference_insert(Geometry1 const& geometry1, \tparam Geometry2 \tparam_geometry \tparam Collection output collection, either a multi-geometry, or a std::vector / std::deque etc +\tparam Strategy \tparam_strategy{Sym_difference} \param geometry1 \param_geometry \param geometry2 \param_geometry \param output_collection the output collection +\param strategy \param_strategy{sym_difference} +\qbk{distinguish,with strategy} \qbk{[include reference/algorithms/sym_difference.qbk]} */ template < typename Geometry1, typename Geometry2, - typename Collection + typename Collection, + typename Strategy > inline void sym_difference(Geometry1 const& geometry1, - Geometry2 const& geometry2, Collection& output_collection) + Geometry2 const& geometry2, + Collection& output_collection, + Strategy const& strategy) { - concepts::check(); - concepts::check(); - - typedef typename boost::range_value::type geometry_out; - concepts::check(); - - typedef typename geometry::rescale_overlay_policy_type + resolve_variant::sym_difference < Geometry1, Geometry2 - >::type rescale_policy_type; + >::apply(geometry1, geometry2, output_collection, strategy); +} - rescale_policy_type robust_policy - = geometry::get_rescale_policy(geometry1, geometry2); - detail::sym_difference::sym_difference_insert( - geometry1, geometry2, robust_policy, - range::back_inserter(output_collection)); +/*! +\brief \brief_calc2{symmetric difference} +\ingroup sym_difference +\details \details_calc2{symmetric difference, spatial set theoretic symmetric difference (XOR)}. +\tparam Geometry1 \tparam_geometry +\tparam Geometry2 \tparam_geometry +\tparam Collection output collection, either a multi-geometry, + or a std::vector / std::deque etc +\param geometry1 \param_geometry +\param geometry2 \param_geometry +\param output_collection the output collection + +\qbk{[include reference/algorithms/sym_difference.qbk]} +*/ +template +< + typename Geometry1, + typename Geometry2, + typename Collection +> +inline void sym_difference(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Collection& output_collection) +{ + resolve_variant::sym_difference + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, output_collection, default_strategy()); } diff --git a/include/boost/geometry/algorithms/union.hpp b/include/boost/geometry/algorithms/union.hpp index e2a7d2b96d..d3a2daf66e 100644 --- a/include/boost/geometry/algorithms/union.hpp +++ b/include/boost/geometry/algorithms/union.hpp @@ -25,6 +25,8 @@ #include #include #include +#include +#include #include #include @@ -234,6 +236,228 @@ inline OutputIterator union_insert(Geometry1 const& geometry1, #endif // DOXYGEN_NO_DETAIL +namespace resolve_strategy { + +struct union_ +{ + template + < + typename Geometry1, + typename Geometry2, + typename RobustPolicy, + typename Collection, + typename Strategy + > + static inline void apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + RobustPolicy const& robust_policy, + Collection & output_collection, + Strategy const& strategy) + { + typedef typename boost::range_value::type geometry_out; + + dispatch::union_insert + < + Geometry1, Geometry2, geometry_out + >::apply(geometry1, geometry2, robust_policy, + range::back_inserter(output_collection), + strategy); + } + + template + < + typename Geometry1, + typename Geometry2, + typename RobustPolicy, + typename Collection + > + static inline void apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + RobustPolicy const& robust_policy, + Collection & output_collection, + default_strategy) + { + typedef typename boost::range_value::type geometry_out; + + typedef typename strategy::intersection::services::default_strategy + < + typename cs_tag::type + >::type strategy_type; + + dispatch::union_insert + < + Geometry1, Geometry2, geometry_out + >::apply(geometry1, geometry2, robust_policy, + range::back_inserter(output_collection), + strategy_type()); + } +}; + +} // resolve_strategy + + +namespace resolve_variant +{ + +template +struct union_ +{ + template + static inline void apply(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + concepts::check(); + concepts::check(); + concepts::check::type>(); + + typedef typename geometry::rescale_overlay_policy_type + < + Geometry1, + Geometry2 + >::type rescale_policy_type; + + rescale_policy_type robust_policy + = geometry::get_rescale_policy(geometry1, + geometry2); + + resolve_strategy::union_::apply(geometry1, geometry2, + robust_policy, + output_collection, + strategy); + } +}; + + +template +struct union_, Geometry2> +{ + template + struct visitor: static_visitor<> + { + Geometry2 const& m_geometry2; + Collection& m_output_collection; + Strategy const& m_strategy; + + visitor(Geometry2 const& geometry2, + Collection& output_collection, + Strategy const& strategy) + : m_geometry2(geometry2) + , m_output_collection(output_collection) + , m_strategy(strategy) + {} + + template + void operator()(Geometry1 const& geometry1) const + { + union_ + < + Geometry1, + Geometry2 + >::apply(geometry1, m_geometry2, m_output_collection, m_strategy); + } + }; + + template + static inline void + apply(variant const& geometry1, + Geometry2 const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + boost::apply_visitor(visitor(geometry2, + output_collection, + strategy), + geometry1); + } +}; + + +template +struct union_ > +{ + template + struct visitor: static_visitor<> + { + Geometry1 const& m_geometry1; + Collection& m_output_collection; + Strategy const& m_strategy; + + visitor(Geometry1 const& geometry1, + Collection& output_collection, + Strategy const& strategy) + : m_geometry1(geometry1) + , m_output_collection(output_collection) + , m_strategy(strategy) + {} + + template + void operator()(Geometry2 const& geometry2) const + { + union_ + < + Geometry1, + Geometry2 + >::apply(m_geometry1, geometry2, m_output_collection, m_strategy); + } + }; + + template + static inline void + apply(Geometry1 const& geometry1, + variant const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + boost::apply_visitor(visitor(geometry1, + output_collection, + strategy), + geometry2); + } +}; + + +template +struct union_, variant > +{ + template + struct visitor: static_visitor<> + { + Collection& m_output_collection; + Strategy const& m_strategy; + + visitor(Collection& output_collection, Strategy const& strategy) + : m_output_collection(output_collection) + , m_strategy(strategy) + {} + + template + void operator()(Geometry1 const& geometry1, + Geometry2 const& geometry2) const + { + union_ + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, m_output_collection, m_strategy); + } + }; + + template + static inline void + apply(variant const& geometry1, + variant const& geometry2, + Collection& output_collection, + Strategy const& strategy) + { + boost::apply_visitor(visitor(output_collection, + strategy), + geometry1, geometry2); + } +}; + +} // namespace resolve_variant /*! @@ -244,31 +468,66 @@ inline OutputIterator union_insert(Geometry1 const& geometry1, \tparam Geometry2 \tparam_geometry \tparam Collection output collection, either a multi-geometry, or a std::vector / std::deque etc +\tparam Strategy \tparam_strategy{Union_} \param geometry1 \param_geometry \param geometry2 \param_geometry \param output_collection the output collection +\param strategy \param_strategy{union_} \note Called union_ because union is a reserved word. +\qbk{distinguish,with strategy} \qbk{[include reference/algorithms/union.qbk]} */ template < typename Geometry1, typename Geometry2, - typename Collection + typename Collection, + typename Strategy > inline void union_(Geometry1 const& geometry1, - Geometry2 const& geometry2, - Collection& output_collection) + Geometry2 const& geometry2, + Collection& output_collection, + Strategy const& strategy) { - concepts::check(); - concepts::check(); + resolve_variant::union_ + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, output_collection, strategy); +} + - typedef typename boost::range_value::type geometry_out; - concepts::check(); +/*! +\brief Combines two geometries which each other +\ingroup union +\details \details_calc2{union, spatial set theoretic union}. +\tparam Geometry1 \tparam_geometry +\tparam Geometry2 \tparam_geometry +\tparam Collection output collection, either a multi-geometry, + or a std::vector / std::deque etc +\param geometry1 \param_geometry +\param geometry2 \param_geometry +\param output_collection the output collection +\note Called union_ because union is a reserved word. - detail::union_::union_insert(geometry1, geometry2, - range::back_inserter(output_collection)); +\qbk{[include reference/algorithms/union.qbk]} +*/ +template +< + typename Geometry1, + typename Geometry2, + typename Collection +> +inline void union_(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Collection& output_collection) +{ + resolve_variant::union_ + < + Geometry1, + Geometry2 + >::apply(geometry1, geometry2, output_collection, default_strategy()); } diff --git a/include/boost/geometry/arithmetic/normalize.hpp b/include/boost/geometry/arithmetic/normalize.hpp new file mode 100644 index 0000000000..7dfdbd2b03 --- /dev/null +++ b/include/boost/geometry/arithmetic/normalize.hpp @@ -0,0 +1,71 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 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_ARITHMETIC_NORMALIZE_HPP +#define BOOST_GEOMETRY_ARITHMETIC_NORMALIZE_HPP + + +#include + +#include +#include +#include + + +namespace boost { namespace geometry +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +template +inline typename coordinate_type::type vec_length_sqr(Point const& pt) +{ + return dot_product(pt, pt); +} + +template +inline typename coordinate_type::type vec_length(Point const& pt) +{ + // NOTE: hypot() could be used instead of sqrt() + return math::sqrt(dot_product(pt, pt)); +} + +template +inline bool vec_normalize(Point & pt, typename coordinate_type::type & len) +{ + typedef typename coordinate_type::type coord_t; + + coord_t const c0 = 0; + len = vec_length(pt); + + if (math::equals(len, c0)) + { + return false; + } + + divide_value(pt, len); + return true; +} + +template +inline bool vec_normalize(Point & pt) +{ + typedef typename coordinate_type::type coord_t; + coord_t len; + return vec_normalize(pt, len); +} + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ARITHMETIC_NORMALIZE_HPP diff --git a/include/boost/geometry/core/srs.hpp b/include/boost/geometry/core/srs.hpp index bf1b4e28a5..5c78b9f1fc 100644 --- a/include/boost/geometry/core/srs.hpp +++ b/include/boost/geometry/core/srs.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 2014. -// Modifications copyright (c) 2014 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 Adam Wulkiewicz, on behalf of Oracle @@ -51,7 +51,7 @@ class spheroid spheroid() : m_a(RadiusType(6378137.0)) - , m_b(RadiusType(6356752.314245)) + , m_b(RadiusType(6356752.3142451793)) {} template @@ -126,8 +126,9 @@ class sphere explicit sphere(RadiusType const& r) : m_r(r) {} + sphere() - : m_r(RadiusType((2.0 * 6378137.0 + 6356752.314245) / 3.0)) + : m_r(RadiusType((2.0 * 6378137.0 + 6356752.3142451793) / 3.0)) {} template diff --git a/include/boost/geometry/formulas/andoyer_inverse.hpp b/include/boost/geometry/formulas/andoyer_inverse.hpp index 04a43995b2..902fd7d8f6 100644 --- a/include/boost/geometry/formulas/andoyer_inverse.hpp +++ b/include/boost/geometry/formulas/andoyer_inverse.hpp @@ -1,6 +1,6 @@ // Boost.Geometry -// Copyright (c) 2015-2016 Oracle and/or its affiliates. +// Copyright (c) 2015-2017 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -20,9 +20,8 @@ #include #include -#include - #include +#include #include @@ -75,7 +74,7 @@ class andoyer_inverse CT const c0 = CT(0); CT const c1 = CT(1); CT const pi = math::pi(); - CT const f = detail::flattening(spheroid); + CT const f = formula::flattening(spheroid); CT const dlon = lon2 - lon1; CT const sin_dlon = sin(dlon); @@ -138,7 +137,14 @@ class andoyer_inverse CT A = c0; CT U = c0; - if ( ! math::equals(cos_lat2, c0) ) + if (math::equals(cos_lat2, c0)) + { + if (sin_lat2 < c0) + { + A = pi; + } + } + else { CT const tan_lat2 = sin_lat2/cos_lat2; CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon; @@ -149,7 +155,14 @@ class andoyer_inverse CT B = c0; CT V = c0; - if ( ! math::equals(cos_lat1, c0) ) + if (math::equals(cos_lat1, c0)) + { + if (sin_lat1 < c0) + { + B = pi; + } + } + else { CT const tan_lat1 = sin_lat1/cos_lat1; CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon; diff --git a/include/boost/geometry/formulas/area_formulas.hpp b/include/boost/geometry/formulas/area_formulas.hpp index a9c503b65b..6a0b525e25 100644 --- a/include/boost/geometry/formulas/area_formulas.hpp +++ b/include/boost/geometry/formulas/area_formulas.hpp @@ -11,7 +11,7 @@ #ifndef BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP #define BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP -#include +#include #include namespace boost { namespace geometry { namespace formula @@ -429,7 +429,7 @@ class area_formulas // Constants CT const ep = spheroid_const.m_ep; - CT const f = geometry::detail::flattening(spheroid_const.m_spheroid); + CT const f = formula::flattening(spheroid_const.m_spheroid); CT const one_minus_f = CT(1) - f; std::size_t const series_order_plus_one = SeriesOrder + 1; std::size_t const series_order_plus_two = SeriesOrder + 2; diff --git a/include/boost/geometry/formulas/differential_quantities.hpp b/include/boost/geometry/formulas/differential_quantities.hpp index 9a92f14e18..ff2ec539db 100644 --- a/include/boost/geometry/formulas/differential_quantities.hpp +++ b/include/boost/geometry/formulas/differential_quantities.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 @@ -64,8 +64,8 @@ class differential_quantities CT const c1 = 1; CT const one_minus_f = c1 - f; - CT const sin_bet1 = one_minus_f * sin_lat1; - CT const sin_bet2 = one_minus_f * sin_lat2; + CT sin_bet1 = one_minus_f * sin_lat1; + CT sin_bet2 = one_minus_f * sin_lat2; // equator if (math::equals(sin_bet1, c0) && math::equals(sin_bet2, c0)) @@ -89,14 +89,17 @@ class differential_quantities CT const e2 = f * (c2 - f); CT const ep2 = e2 / math::sqr(one_minus_f); - CT const cos_bet1 = cos_lat1; - CT const cos_bet2 = cos_lat2; - CT const sin_alp1 = sin(azimuth); CT const cos_alp1 = cos(azimuth); //CT const sin_alp2 = sin(reverse_azimuth); CT const cos_alp2 = cos(reverse_azimuth); + CT cos_bet1 = cos_lat1; + CT cos_bet2 = cos_lat2; + + normalize(sin_bet1, cos_bet1); + normalize(sin_bet2, cos_bet2); + CT sin_sig1 = sin_bet1; CT cos_sig1 = cos_alp1 * cos_bet1; CT sin_sig2 = sin_bet2; @@ -112,8 +115,8 @@ class differential_quantities J12_f(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, f) : J12_ep_sqr(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, ep2) ; - CT const dn1 = math::sqrt(c1 + e2 * math::sqr(sin_lat1)); - CT const dn2 = math::sqrt(c1 + e2 * math::sqr(sin_lat2)); + CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_bet1)); + CT const dn2 = math::sqrt(c1 + ep2 * math::sqr(sin_bet2)); if (BOOST_GEOMETRY_CONDITION(EnableReducedLength)) { diff --git a/include/boost/geometry/formulas/eccentricity_sqr.hpp b/include/boost/geometry/formulas/eccentricity_sqr.hpp new file mode 100644 index 0000000000..01a9beacb9 --- /dev/null +++ b/include/boost/geometry/formulas/eccentricity_sqr.hpp @@ -0,0 +1,70 @@ +// Boost.Geometry + +// Copyright (c) 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_FORMULAS_ECCENCRICITY_SQR_HPP +#define BOOST_GEOMETRY_FORMULAS_ECCENCRICITY_SQR_HPP + +#include +#include +#include + +#include + +namespace boost { namespace geometry +{ + +#ifndef DOXYGEN_NO_DISPATCH +namespace formula_dispatch +{ + +template ::type> +struct eccentricity_sqr + : not_implemented +{}; + +template +struct eccentricity_sqr +{ + static inline ResultType apply(Geometry const& /*geometry*/) + { + return ResultType(0); + } +}; + +template +struct eccentricity_sqr +{ + static inline ResultType apply(Geometry const& geometry) + { + // 1 - (b / a)^2 + return ResultType(1) - math::sqr(ResultType(get_radius<2>(geometry)) + / ResultType(get_radius<0>(geometry))); + } +}; + +} // namespace formula_dispatch +#endif // DOXYGEN_NO_DISPATCH + +#ifndef DOXYGEN_NO_DETAIL +namespace formula +{ + +template +ResultType eccentricity_sqr(Geometry const& geometry) +{ + return formula_dispatch::eccentricity_sqr::apply(geometry); +} + +} // namespace formula +#endif // DOXYGEN_NO_DETAIL + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_FORMULAS_ECCENCRICITY_SQR_HPP diff --git a/include/boost/geometry/algorithms/detail/flattening.hpp b/include/boost/geometry/formulas/flattening.hpp similarity index 78% rename from include/boost/geometry/algorithms/detail/flattening.hpp rename to include/boost/geometry/formulas/flattening.hpp index 8ed5fd9a89..f94ead65b0 100644 --- a/include/boost/geometry/algorithms/detail/flattening.hpp +++ b/include/boost/geometry/formulas/flattening.hpp @@ -1,6 +1,6 @@ // Boost.Geometry -// Copyright (c) 2014 Oracle and/or its affiliates. +// Copyright (c) 2014-2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -8,8 +8,8 @@ // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) -#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_FLATTENING_HPP -#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_FLATTENING_HPP +#ifndef BOOST_GEOMETRY_FORMULAS_FLATTENING_HPP +#define BOOST_GEOMETRY_FORMULAS_FLATTENING_HPP #include #include @@ -21,7 +21,7 @@ namespace boost { namespace geometry { #ifndef DOXYGEN_NO_DISPATCH -namespace detail_dispatch +namespace formula_dispatch { template ::type> @@ -43,27 +43,28 @@ struct flattening { static inline ResultType apply(Geometry const& geometry) { + // (a - b) / a return ResultType(get_radius<0>(geometry) - get_radius<2>(geometry)) / ResultType(get_radius<0>(geometry)); } }; -} // namespace detail_dispatch +} // namespace formula_dispatch #endif // DOXYGEN_NO_DISPATCH #ifndef DOXYGEN_NO_DETAIL -namespace detail +namespace formula { template ResultType flattening(Geometry const& geometry) { - return detail_dispatch::flattening::apply(geometry); + return formula_dispatch::flattening::apply(geometry); } -} // namespace detail +} // namespace formula #endif // DOXYGEN_NO_DETAIL }} // namespace boost::geometry -#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_FLATTENING_HPP +#endif // BOOST_GEOMETRY_FORMULAS_FLATTENING_HPP diff --git a/include/boost/geometry/formulas/geographic.hpp b/include/boost/geometry/formulas/geographic.hpp new file mode 100644 index 0000000000..f6feb66633 --- /dev/null +++ b/include/boost/geometry/formulas/geographic.hpp @@ -0,0 +1,457 @@ +// Boost.Geometry + +// Copyright (c) 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_FORMULAS_GEOGRAPHIC_HPP +#define BOOST_GEOMETRY_FORMULAS_GEOGRAPHIC_HPP + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +namespace boost { namespace geometry { + +namespace formula { + +template +inline Point3d geo_to_cart3d(PointGeo const& point_geo, Spheroid const& spheroid) +{ + typedef typename coordinate_type::type calc_t; + + calc_t const c1 = 1; + calc_t const e_sqr = eccentricity_sqr(spheroid); + + calc_t const lon = get_as_radian<0>(point_geo); + calc_t const lat = get_as_radian<1>(point_geo); + + Point3d res; + + calc_t const sin_lat = sin(lat); + + // "unit" spheroid, a = 1 + calc_t const N = c1 / math::sqrt(c1 - e_sqr * math::sqr(sin_lat)); + calc_t const N_cos_lat = N * cos(lat); + + set<0>(res, N_cos_lat * cos(lon)); + set<1>(res, N_cos_lat * sin(lon)); + set<2>(res, N * (c1 - e_sqr) * sin_lat); + + return res; +} + +template +inline void geo_to_cart3d(PointGeo const& point_geo, Point3d & result, Point3d & north, Point3d & east, Spheroid const& spheroid) +{ + typedef typename coordinate_type::type calc_t; + + calc_t const c1 = 1; + calc_t const e_sqr = eccentricity_sqr(spheroid); + + calc_t const lon = get_as_radian<0>(point_geo); + calc_t const lat = get_as_radian<1>(point_geo); + + calc_t const sin_lon = sin(lon); + calc_t const cos_lon = cos(lon); + calc_t const sin_lat = sin(lat); + calc_t const cos_lat = cos(lat); + + // "unit" spheroid, a = 1 + calc_t const N = c1 / math::sqrt(c1 - e_sqr * math::sqr(sin_lat)); + calc_t const N_cos_lat = N * cos_lat; + + set<0>(result, N_cos_lat * cos_lon); + set<1>(result, N_cos_lat * sin_lon); + set<2>(result, N * (c1 - e_sqr) * sin_lat); + + set<0>(east, -sin_lon); + set<1>(east, cos_lon); + set<2>(east, 0); + + set<0>(north, -sin_lat * cos_lon); + set<1>(north, -sin_lat * sin_lon); + set<2>(north, cos_lat); +} + +template +inline PointGeo cart3d_to_geo(Point3d const& point_3d, Spheroid const& spheroid) +{ + typedef typename coordinate_type::type coord_t; + typedef typename coordinate_type::type calc_t; + + calc_t const c1 = 1; + //calc_t const c2 = 2; + calc_t const e_sqr = eccentricity_sqr(spheroid); + + calc_t const x = get<0>(point_3d); + calc_t const y = get<1>(point_3d); + calc_t const z = get<2>(point_3d); + calc_t const xy_l = math::sqrt(math::sqr(x) + math::sqr(y)); + + calc_t const lonr = atan2(y, x); + + // NOTE: Alternative version + // http://www.iag-aig.org/attach/989c8e501d9c5b5e2736955baf2632f5/V60N2_5FT.pdf + // calc_t const lonr = c2 * atan2(y, x + xy_l); + + calc_t const latr = atan2(z, (c1 - e_sqr) * xy_l); + + // NOTE: If h is equal to 0 then there is no need to improve value of latitude + // because then N_i / (N_i + h_i) = 1 + // http://www.navipedia.net/index.php/Ellipsoidal_and_Cartesian_Coordinates_Conversion + + PointGeo res; + + set_from_radian<0>(res, lonr); + set_from_radian<1>(res, latr); + + coord_t lon = get<0>(res); + coord_t lat = get<1>(res); + + math::normalize_spheroidal_coordinates + < + typename coordinate_system::type::units, + coord_t + >(lon, lat); + + set<0>(res, lon); + set<1>(res, lat); + + return res; +} + +template +inline Point3d projected_to_xy(Point3d const& point_3d, Spheroid const& spheroid) +{ + typedef typename coordinate_type::type coord_t; + + // len_xy = sqrt(x^2 + y^2) + // r = len_xy - |z / tan(lat)| + // assuming h = 0 + // lat = atan2(z, (1 - e^2) * len_xy); + // |z / tan(lat)| = (1 - e^2) * len_xy + // r = e^2 * len_xy + // x_res = r * cos(lon) = e^2 * len_xy * x / len_xy = e^2 * x + // y_res = r * sin(lon) = e^2 * len_xy * y / len_xy = e^2 * y + + coord_t const c0 = 0; + coord_t const e_sqr = formula::eccentricity_sqr(spheroid); + + Point3d res; + + set<0>(res, e_sqr * get<0>(point_3d)); + set<1>(res, e_sqr * get<1>(point_3d)); + set<2>(res, c0); + + return res; +} + +template +inline Point3d projected_to_surface(Point3d const& direction, Spheroid const& spheroid) +{ + typedef typename coordinate_type::type coord_t; + + //coord_t const c0 = 0; + coord_t const c2 = 2; + coord_t const c4 = 4; + + // calculate the point of intersection of a ray and spheroid's surface + // the origin is the origin of the coordinate system + //(x*x+y*y)/(a*a) + z*z/(b*b) = 1 + // x = d.x * t + // y = d.y * t + // z = d.z * t + coord_t const dx = get<0>(direction); + coord_t const dy = get<1>(direction); + coord_t const dz = get<2>(direction); + + //coord_t const a_sqr = math::sqr(get_radius<0>(spheroid)); + //coord_t const b_sqr = math::sqr(get_radius<2>(spheroid)); + // "unit" spheroid, a = 1 + coord_t const a_sqr = 1; + coord_t const b_sqr = math::sqr(get_radius<2>(spheroid) / get_radius<0>(spheroid)); + + coord_t const param_a = (dx*dx + dy*dy) / a_sqr + dz*dz / b_sqr; + coord_t const delta = c4 * param_a; + // delta >= 0 + coord_t const t = math::sqrt(delta) / (c2 * param_a); + + // result = direction * t + Point3d result = direction; + multiply_value(result, t); + + return result; +} + +template +inline bool projected_to_surface(Point3d const& origin, Point3d const& direction, Point3d & result1, Point3d & result2, Spheroid const& spheroid) +{ + typedef typename coordinate_type::type coord_t; + + coord_t const c0 = 0; + coord_t const c1 = 1; + coord_t const c2 = 2; + coord_t const c4 = 4; + + // calculate the point of intersection of a ray and spheroid's surface + //(x*x+y*y)/(a*a) + z*z/(b*b) = 1 + // x = o.x + d.x * t + // y = o.y + d.y * t + // z = o.z + d.z * t + coord_t const ox = get<0>(origin); + coord_t const oy = get<1>(origin); + coord_t const oz = get<2>(origin); + coord_t const dx = get<0>(direction); + coord_t const dy = get<1>(direction); + coord_t const dz = get<2>(direction); + + //coord_t const a_sqr = math::sqr(get_radius<0>(spheroid)); + //coord_t const b_sqr = math::sqr(get_radius<2>(spheroid)); + // "unit" spheroid, a = 1 + coord_t const a_sqr = 1; + coord_t const b_sqr = math::sqr(get_radius<2>(spheroid) / get_radius<0>(spheroid)); + + coord_t const param_a = (dx*dx + dy*dy) / a_sqr + dz*dz / b_sqr; + coord_t const param_b = c2 * ((ox*dx + oy*dy) / a_sqr + oz*dz / b_sqr); + coord_t const param_c = (ox*ox + oy*oy) / a_sqr + oz*oz / b_sqr - c1; + + coord_t const delta = math::sqr(param_b) - c4 * param_a*param_c; + + // equals() ? + if (delta < c0 || param_a == 0) + { + return false; + } + + // result = origin + direction * t + + coord_t const sqrt_delta = math::sqrt(delta); + coord_t const two_a = c2 * param_a; + + coord_t const t1 = (-param_b + sqrt_delta) / two_a; + result1 = direction; + multiply_value(result1, t1); + add_point(result1, origin); + + coord_t const t2 = (-param_b - sqrt_delta) / two_a; + result2 = direction; + multiply_value(result2, t2); + add_point(result2, origin); + + return true; +} + +template +inline bool great_elliptic_intersection(Point3d const& a1, Point3d const& a2, + Point3d const& b1, Point3d const& b2, + Point3d & result, + Spheroid const& spheroid) +{ + typedef typename coordinate_type::type coord_t; + + coord_t c0 = 0; + coord_t c1 = 1; + + Point3d n1 = cross_product(a1, a2); + Point3d n2 = cross_product(b1, b2); + + // intersection direction + Point3d id = cross_product(n1, n2); + coord_t id_len_sqr = dot_product(id, id); + + if (math::equals(id_len_sqr, c0)) + { + return false; + } + + // no need to normalize a1 and a2 because the intersection point on + // the opposite side of the globe is at the same distance from the origin + coord_t cos_a1i = dot_product(a1, id); + coord_t cos_a2i = dot_product(a2, id); + coord_t gri = math::detail::greatest(cos_a1i, cos_a2i); + Point3d neg_id = id; + multiply_value(neg_id, -c1); + coord_t cos_a1ni = dot_product(a1, neg_id); + coord_t cos_a2ni = dot_product(a2, neg_id); + coord_t grni = math::detail::greatest(cos_a1ni, cos_a2ni); + + if (gri >= grni) + { + result = projected_to_surface(id, spheroid); + } + else + { + result = projected_to_surface(neg_id, spheroid); + } + + return true; +} + +template +static inline int elliptic_side_value(Point3d1 const& origin, Point3d1 const& norm, Point3d2 const& pt) +{ + typedef typename coordinate_type::type calc_t; + calc_t c0 = 0; + + // vector oposite to pt - origin + // only for the purpose of assigning origin + Point3d1 vec = origin; + subtract_point(vec, pt); + + calc_t d = dot_product(norm, vec); + + // since the vector is opposite the signs are opposite + return math::equals(d, c0) ? 0 + : d < c0 ? 1 + : -1; // d > 0 +} + +template +inline bool planes_spheroid_intersection(Point3d const& o1, Point3d const& n1, + Point3d const& o2, Point3d const& n2, + Point3d & ip1, Point3d & ip2, + Spheroid const& spheroid) +{ + typedef typename coordinate_type::type coord_t; + + coord_t c0 = 0; + coord_t c1 = 1; + + // Below + // n . (p - o) = 0 + // n . p - n . o = 0 + // n . p + d = 0 + // n . p = h + + // intersection direction + Point3d id = cross_product(n1, n2); + + if (math::equals(dot_product(id, id), c0)) + { + return false; + } + + coord_t dot_n1_n2 = dot_product(n1, n2); + coord_t dot_n1_n2_sqr = math::sqr(dot_n1_n2); + + coord_t h1 = dot_product(n1, o1); + coord_t h2 = dot_product(n2, o2); + + coord_t denom = c1 - dot_n1_n2_sqr; + coord_t C1 = (h1 - h2 * dot_n1_n2) / denom; + coord_t C2 = (h2 - h1 * dot_n1_n2) / denom; + + // C1 * n1 + C2 * n2 + Point3d C1_n1 = n1; + multiply_value(C1_n1, C1); + Point3d C2_n2 = n2; + multiply_value(C2_n2, C2); + Point3d io = C1_n1; + add_point(io, C2_n2); + + if (! projected_to_surface(io, id, ip1, ip2, spheroid)) + { + return false; + } + + return true; +} + + +template +inline void experimental_elliptic_plane(Point3d const& p1, Point3d const& p2, + Point3d & v1, Point3d & v2, + Point3d & origin, Point3d & normal, + Spheroid const& spheroid) +{ + typedef typename coordinate_type::type coord_t; + + Point3d xy1 = projected_to_xy(p1, spheroid); + Point3d xy2 = projected_to_xy(p2, spheroid); + + // origin = (xy1 + xy2) / 2 + origin = xy1; + add_point(origin, xy2); + multiply_value(origin, coord_t(0.5)); + + // v1 = p1 - origin + v1 = p1; + subtract_point(v1, origin); + // v2 = p2 - origin + v2 = p2; + subtract_point(v2, origin); + + normal = cross_product(v1, v2); +} + +template +inline void experimental_elliptic_plane(Point3d const& p1, Point3d const& p2, + Point3d & origin, Point3d & normal, + Spheroid const& spheroid) +{ + Point3d v1, v2; + experimental_elliptic_plane(p1, p2, v1, v2, origin, normal, spheroid); +} + +template +inline bool experimental_elliptic_intersection(Point3d const& a1, Point3d const& a2, + Point3d const& b1, Point3d const& b2, + Point3d & result, + Spheroid const& spheroid) +{ + typedef typename coordinate_type::type coord_t; + + coord_t c0 = 0; + coord_t c1 = 1; + + Point3d a1v, a2v, o1, n1; + experimental_elliptic_plane(a1, a2, a1v, a2v, o1, n1, spheroid); + Point3d b1v, b2v, o2, n2; + experimental_elliptic_plane(b1, b2, b1v, b2v, o2, n2, spheroid); + + if (! detail::vec_normalize(n1) || ! detail::vec_normalize(n2)) + { + return false; + } + + Point3d ip1_s, ip2_s; + if (! planes_spheroid_intersection(o1, n1, o2, n2, ip1_s, ip2_s, spheroid)) + { + return false; + } + + // NOTE: simplified test, may not work in all cases + coord_t dot_a1i1 = dot_product(a1, ip1_s); + coord_t dot_a2i1 = dot_product(a2, ip1_s); + coord_t gri1 = math::detail::greatest(dot_a1i1, dot_a2i1); + coord_t dot_a1i2 = dot_product(a1, ip2_s); + coord_t dot_a2i2 = dot_product(a2, ip2_s); + coord_t gri2 = math::detail::greatest(dot_a1i2, dot_a2i2); + + result = gri1 >= gri2 ? ip1_s : ip2_s; + + return true; +} + +} // namespace formula + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_FORMULAS_GEOGRAPHIC_HPP diff --git a/include/boost/geometry/formulas/gnomonic_spheroid.hpp b/include/boost/geometry/formulas/gnomonic_spheroid.hpp index 3457397b0f..4710c6c063 100644 --- a/include/boost/geometry/formulas/gnomonic_spheroid.hpp +++ b/include/boost/geometry/formulas/gnomonic_spheroid.hpp @@ -14,12 +14,11 @@ #include -#include - #include #include #include +#include #include #include #include diff --git a/include/boost/geometry/formulas/sjoberg_intersection.hpp b/include/boost/geometry/formulas/sjoberg_intersection.hpp index 03bd4bc97e..cf9e8c97b8 100644 --- a/include/boost/geometry/formulas/sjoberg_intersection.hpp +++ b/include/boost/geometry/formulas/sjoberg_intersection.hpp @@ -20,19 +20,606 @@ #include #include -#include +#include +#include namespace boost { namespace geometry { namespace formula { +/*! +\brief The intersection of two great circles as proposed by Sjoberg. +\see See + - [Sjoberg02] Lars E. Sjoberg, Intersections on the sphere and ellipsoid, 2002 + http://link.springer.com/article/10.1007/s00190-001-0230-9 +*/ +template +struct sjoberg_intersection_spherical_02 +{ + // TODO: if it will be used as standalone formula + // support segments on equator and endpoints on poles + + static inline bool apply(CT const& lon1, CT const& lat1, CT const& lon_a2, CT const& lat_a2, + CT const& lon2, CT const& lat2, CT const& lon_b2, CT const& lat_b2, + CT & lon, CT & lat) + { + CT tan_lat = 0; + bool res = apply_alt(lon1, lat1, lon_a2, lat_a2, + lon2, lat2, lon_b2, lat_b2, + lon, tan_lat); + + if (res) + { + lat = atan(tan_lat); + } + + return res; + } + + static inline bool apply_alt(CT const& lon1, CT const& lat1, CT const& lon_a2, CT const& lat_a2, + CT const& lon2, CT const& lat2, CT const& lon_b2, CT const& lat_b2, + CT & lon, CT & tan_lat) + { + CT const cos_lon1 = cos(lon1); + CT const sin_lon1 = sin(lon1); + CT const cos_lon2 = cos(lon2); + CT const sin_lon2 = sin(lon2); + CT const sin_lat1 = sin(lat1); + CT const sin_lat2 = sin(lat2); + CT const cos_lat1 = cos(lat1); + CT const cos_lat2 = cos(lat2); + + CT const tan_lat_a2 = tan(lat_a2); + CT const tan_lat_b2 = tan(lat_b2); + + return apply(lon1, lon_a2, lon2, lon_b2, + sin_lon1, cos_lon1, sin_lat1, cos_lat1, + sin_lon2, cos_lon2, sin_lat2, cos_lat2, + tan_lat_a2, tan_lat_b2, + lon, tan_lat); + } + +private: + static inline bool apply(CT const& lon1, CT const& lon_a2, CT const& lon2, CT const& lon_b2, + CT const& sin_lon1, CT const& cos_lon1, CT const& sin_lat1, CT const& cos_lat1, + CT const& sin_lon2, CT const& cos_lon2, CT const& sin_lat2, CT const& cos_lat2, + CT const& tan_lat_a2, CT const& tan_lat_b2, + CT & lon, CT & tan_lat) + { + // NOTE: + // cos_lat_ = 0 <=> segment on equator + // tan_alpha_ = 0 <=> segment vertical + + CT const tan_lat1 = sin_lat1 / cos_lat1; //tan(lat1); + CT const tan_lat2 = sin_lat2 / cos_lat2; //tan(lat2); + + CT const dlon1 = lon_a2 - lon1; + CT const sin_dlon1 = sin(dlon1); + CT const dlon2 = lon_b2 - lon2; + CT const sin_dlon2 = sin(dlon2); + + CT const c0 = 0; + bool const is_vertical1 = math::equals(sin_dlon1, c0); + bool const is_vertical2 = math::equals(sin_dlon2, c0); + + CT tan_alpha1 = 0; + CT tan_alpha2 = 0; + + if (is_vertical1 && is_vertical2) + { + // circles intersect at one of the poles or are collinear + return false; + } + else if (is_vertical1) + { + CT const cos_dlon2 = cos(dlon2); + CT const tan_alpha2_x = cos_lat2 * tan_lat_b2 - sin_lat2 * cos_dlon2; + tan_alpha2 = sin_dlon2 / tan_alpha2_x; + + lon = lon1; + } + else if (is_vertical2) + { + CT const cos_dlon1 = cos(dlon1); + CT const tan_alpha1_x = cos_lat1 * tan_lat_a2 - sin_lat1 * cos_dlon1; + tan_alpha1 = sin_dlon1 / tan_alpha1_x; + + lon = lon2; + } + else + { + CT const cos_dlon1 = cos(dlon1); + CT const cos_dlon2 = cos(dlon2); + + CT const tan_alpha1_x = cos_lat1 * tan_lat_a2 - sin_lat1 * cos_dlon1; + CT const tan_alpha2_x = cos_lat2 * tan_lat_b2 - sin_lat2 * cos_dlon2; + tan_alpha1 = sin_dlon1 / tan_alpha1_x; + tan_alpha2 = sin_dlon2 / tan_alpha2_x; + + CT const T1 = tan_alpha1 * cos_lat1; + CT const T2 = tan_alpha2 * cos_lat2; + CT const T1T2 = T1*T2; + CT const tan_lon_y = T1 * sin_lon2 - T2 * sin_lon1 + T1T2 * (tan_lat1 * cos_lon1 - tan_lat2 * cos_lon2); + CT const tan_lon_x = T1 * cos_lon2 - T2 * cos_lon1 - T1T2 * (tan_lat1 * sin_lon1 - tan_lat2 * sin_lon2); + + lon = atan2(tan_lon_y, tan_lon_x); + } + + // choose closer result + CT const pi = math::pi(); + CT const lon_2 = lon > c0 ? lon - pi : lon + pi; + CT const lon_dist1 = (std::max)((std::min)(math::longitude_difference(lon1, lon), + math::longitude_difference(lon_a2, lon)), + (std::min)(math::longitude_difference(lon2, lon), + math::longitude_difference(lon_b2, lon))); + CT const lon_dist2 = (std::max)((std::min)(math::longitude_difference(lon1, lon_2), + math::longitude_difference(lon_a2, lon_2)), + (std::min)(math::longitude_difference(lon2, lon_2), + math::longitude_difference(lon_b2, lon_2))); + if (lon_dist2 < lon_dist1) + { + lon = lon_2; + } + + CT const sin_lon = sin(lon); + CT const cos_lon = cos(lon); + + if (math::abs(tan_alpha1) >= math::abs(tan_alpha2)) // pick less vertical segment + { + CT const sin_dlon_1 = sin_lon * cos_lon1 - cos_lon * sin_lon1; + CT const cos_dlon_1 = cos_lon * cos_lon1 + sin_lon * sin_lon1; + CT const lat_y_1 = sin_dlon_1 + tan_alpha1 * sin_lat1 * cos_dlon_1; + CT const lat_x_1 = tan_alpha1 * cos_lat1; + tan_lat = lat_y_1 / lat_x_1; + } + else + { + CT const sin_dlon_2 = sin_lon * cos_lon2 - cos_lon * sin_lon2; + CT const cos_dlon_2 = cos_lon * cos_lon2 + sin_lon * sin_lon2; + CT const lat_y_2 = sin_dlon_2 + tan_alpha2 * sin_lat2 * cos_dlon_2; + CT const lat_x_2 = tan_alpha2 * cos_lat2; + tan_lat = lat_y_2 / lat_x_2; + } + + return true; + } +}; + + +/*! Approximation of dLambda_j [Sjoberg07], expanded into taylor series in e^2 + Maxima script: + dLI_j(c_j, sinB_j, sinB) := integrate(1 / (sqrt(1 - c_j ^ 2 - x ^ 2)*(1 + sqrt(1 - e2*(1 - x ^ 2)))), x, sinB_j, sinB); + dL_j(c_j, B_j, B) := -e2 * c_j * dLI_j(c_j, B_j, B); + S: taylor(dLI_j(c_j, sinB_j, sinB), e2, 0, 3); + assume(c_j < 1); + assume(c_j > 0); + L1: factor(integrate(sqrt(-x ^ 2 - c_j ^ 2 + 1) / (x ^ 2 + c_j ^ 2 - 1), x)); + L2: factor(integrate(((x ^ 2 - 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x)); + L3: factor(integrate(((x ^ 4 - 2 * x ^ 2 + 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x)); + L4: factor(integrate(((x ^ 6 - 3 * x ^ 4 + 3 * x ^ 2 - 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x)); + +\see See + - [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007 + http://link.springer.com/article/10.1007/s00190-007-0204-7 +*/ +template +inline CT sjoberg_d_lambda_e_sqr(CT const& sin_betaj, CT const& sin_beta, + CT const& Cj, CT const& sqrt_1_Cj_sqr, + CT const& e_sqr) +{ + using math::detail::bounded; + + if (Order == 0) + { + return 0; + } + + CT const c1 = 1; + CT const c2 = 2; + + CT const asin_B = asin(bounded(sin_beta / sqrt_1_Cj_sqr, -c1, c1)); + CT const asin_Bj = asin(sin_betaj / sqrt_1_Cj_sqr); + CT const L0 = (asin_B - asin_Bj) / c2; + + if (Order == 1) + { + return -Cj * e_sqr * L0; + } + + CT const c0 = 0; + CT const c16 = 16; + + CT const X = sin_beta; + CT const Xj = sin_betaj; + CT const X_sqr = math::sqr(X); + CT const Xj_sqr = math::sqr(Xj); + CT const Cj_sqr = math::sqr(Cj); + CT const Cj_sqr_plus_one = Cj_sqr + c1; + CT const one_minus_Cj_sqr = c1 - Cj_sqr; + CT const sqrt_Y = math::sqrt(bounded(-X_sqr + one_minus_Cj_sqr, c0)); + CT const sqrt_Yj = math::sqrt(-Xj_sqr + one_minus_Cj_sqr); + CT const L1 = (Cj_sqr_plus_one * (asin_B - asin_Bj) + X * sqrt_Y - Xj * sqrt_Yj) / c16; + + if (Order == 2) + { + return -Cj * e_sqr * (L0 + e_sqr * L1); + } + + CT const c3 = 3; + CT const c5 = 5; + CT const c128 = 128; + + CT const E = Cj_sqr * (c3 * Cj_sqr + c2) + c3; + CT const F = X * (-c2 * X_sqr + c3 * Cj_sqr + c5); + CT const Fj = Xj * (-c2 * Xj_sqr + c3 * Cj_sqr + c5); + CT const L2 = (E * (asin_B - asin_Bj) + F * sqrt_Y - Fj * sqrt_Yj) / c128; + + if (Order == 3) + { + return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * L2)); + } + + CT const c8 = 8; + CT const c9 = 9; + CT const c10 = 10; + CT const c15 = 15; + CT const c24 = 24; + CT const c26 = 26; + CT const c33 = 33; + CT const c6144 = 6144; + + CT const G = Cj_sqr * (Cj_sqr * (Cj_sqr * c15 + c9) + c9) + c15; + CT const H = -c10 * Cj_sqr - c26; + CT const I = Cj_sqr * (Cj_sqr * c15 + c24) + c33; + CT const J = X_sqr * (X * (c8 * X_sqr + H)) + X * I; + CT const Jj = Xj_sqr * (Xj * (c8 * Xj_sqr + H)) + Xj * I; + CT const L3 = (G * (asin_B - asin_Bj) + J * sqrt_Y - Jj * sqrt_Yj) / c6144; + + // Order 4 and higher + return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * (L2 + e_sqr * L3))); +} + +/*! +\brief The representation of geodesic as proposed by Sjoberg. +\see See + - [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007 + http://link.springer.com/article/10.1007/s00190-007-0204-7 + - [Sjoberg12] Lars E. Sjoberg, Solutions to the ellipsoidal Clairaut constant + and the inverse geodetic problem by numerical integration, 2012 + https://www.degruyter.com/view/j/jogs.2012.2.issue-3/v10156-011-0037-4/v10156-011-0037-4.xml +*/ +template +class sjoberg_geodesic +{ + sjoberg_geodesic() {} + + static int sign_C(CT const& alphaj) + { + CT const c0 = 0; + CT const c2 = 2; + CT const pi = math::pi(); + CT const pi_half = pi / c2; + + return (pi_half < alphaj && alphaj < pi) || (-pi_half < alphaj && alphaj < c0) ? -1 : 1; + } + +public: + sjoberg_geodesic(CT const& lon, CT const& lat, CT const& alpha, CT const& f) + : lonj(lon) + , latj(lat) + , alphaj(alpha) + { + CT const c0 = 0; + CT const c1 = 1; + CT const c2 = 2; + //CT const pi = math::pi(); + //CT const pi_half = pi / c2; + + one_minus_f = c1 - f; + e_sqr = f * (c2 - f); + + tan_latj = tan(lat); + tan_betaj = one_minus_f * tan_latj; + betaj = atan(tan_betaj); + sin_betaj = sin(betaj); + + cos_betaj = cos(betaj); + sin_alphaj = sin(alphaj); + // Clairaut constant (lower-case in the paper) + Cj = sign_C(alphaj) * cos_betaj * sin_alphaj; + Cj_sqr = math::sqr(Cj); + sqrt_1_Cj_sqr = math::sqrt(c1 - Cj_sqr); + + sign_lon_diff = alphaj >= 0 ? 1 : -1; // || alphaj == -pi ? + //sign_lon_diff = 1; + + is_on_equator = math::equals(sqrt_1_Cj_sqr, c0); + is_Cj_zero = math::equals(Cj, c0); + + t0j = c0; + asin_tj_t0j = c0; + + if (! is_Cj_zero) + { + t0j = sqrt_1_Cj_sqr / Cj; + } + + if (! is_on_equator) + { + //asin_tj_t0j = asin(tan_betaj / t0j); + asin_tj_t0j = asin(tan_betaj * Cj / sqrt_1_Cj_sqr); + } + } + + struct vertex_data + { + //CT beta0j; + CT sin_beta0j; + CT dL0j; + CT lon0j; + }; + + vertex_data get_vertex_data() const + { + CT const c2 = 2; + CT const pi = math::pi(); + CT const pi_half = pi / c2; + + vertex_data res; + + if (! is_Cj_zero) + { + //res.beta0j = atan(t0j); + //res.sin_beta0j = sin(res.beta0j); + res.sin_beta0j = math::sign(t0j) * sqrt_1_Cj_sqr; + res.dL0j = d_lambda(res.sin_beta0j); + res.lon0j = lonj + sign_lon_diff * (pi_half - asin_tj_t0j + res.dL0j); + } + else + { + //res.beta0j = pi_half; + //res.sin_beta0j = betaj >= 0 ? 1 : -1; + res.sin_beta0j = 1; + res.dL0j = 0; + res.lon0j = lonj; + } + + return res; + } + + bool is_sin_beta_ok(CT const& sin_beta) const + { + CT const c1 = 1; + return math::abs(sin_beta / sqrt_1_Cj_sqr) <= c1; + } + + bool k_diff(CT const& sin_beta, + CT & delta_k) const + { + if (is_Cj_zero) + { + delta_k = 0; + return true; + } + + // beta out of bounds and not close + if (! (is_sin_beta_ok(sin_beta) + || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) ) + { + return false; + } + + // NOTE: beta may be slightly out of bounds here but d_lambda handles that + CT const dLj = d_lambda(sin_beta); + delta_k = sign_lon_diff * (/*asin_t_t0j*/ - asin_tj_t0j + dLj); + + return true; + } + + bool lon_diff(CT const& sin_beta, CT const& t, + CT & delta_lon) const + { + using math::detail::bounded; + CT const c1 = 1; + + if (is_Cj_zero) + { + delta_lon = 0; + return true; + } + + CT delta_k = 0; + if (! k_diff(sin_beta, delta_k)) + { + return false; + } + + CT const t_t0j = t / t0j; + // NOTE: t may be slightly out of bounds here + CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1)); + delta_lon = sign_lon_diff * asin_t_t0j + delta_k; + + return true; + } + + bool k_diffs(CT const& sin_beta, vertex_data const& vd, + CT & delta_k_before, CT & delta_k_behind, + bool check_sin_beta = true) const + { + CT const pi = math::pi(); + + if (is_Cj_zero) + { + delta_k_before = 0; + delta_k_behind = sign_lon_diff * pi; + return true; + } + + // beta out of bounds and not close + if (check_sin_beta + && ! (is_sin_beta_ok(sin_beta) + || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) ) + { + return false; + } + + // NOTE: beta may be slightly out of bounds here but d_lambda handles that + CT const dLj = d_lambda(sin_beta); + delta_k_before = sign_lon_diff * (/*asin_t_t0j*/ - asin_tj_t0j + dLj); + + // This version require no additional dLj calculation + delta_k_behind = sign_lon_diff * (pi /*- asin_t_t0j*/ - asin_tj_t0j + vd.dL0j + (vd.dL0j - dLj)); + + // [Sjoberg12] + //CT const dL101 = d_lambda(sin_betaj, vd.sin_beta0j); + // WARNING: the following call might not work if beta was OoB because only the second argument is bounded + //CT const dL_01 = d_lambda(sin_beta, vd.sin_beta0j); + //delta_k_behind = sign_lon_diff * (pi /*- asin_t_t0j*/ - asin_tj_t0j + dL101 + dL_01); + + return true; + } + + bool lon_diffs(CT const& sin_beta, CT const& t, vertex_data const& vd, + CT & delta_lon_before, CT & delta_lon_behind) const + { + using math::detail::bounded; + CT const c1 = 1; + CT const pi = math::pi(); + + if (is_Cj_zero) + { + delta_lon_before = 0; + delta_lon_behind = sign_lon_diff * pi; + return true; + } + + CT delta_k_before = 0, delta_k_behind = 0; + if (! k_diffs(sin_beta, vd, delta_k_before, delta_k_behind)) + { + return false; + } + + CT const t_t0j = t / t0j; + // NOTE: t may be slightly out of bounds here + CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1)); + CT const sign_asin_t_t0j = sign_lon_diff * asin_t_t0j; + delta_lon_before = sign_asin_t_t0j + delta_k_before; + delta_lon_behind = -sign_asin_t_t0j + delta_k_behind; + + return true; + } + + bool lon(CT const& sin_beta, CT const& t, vertex_data const& vd, + CT & lon_before, CT & lon_behind) const + { + using math::detail::bounded; + CT const c1 = 1; + CT const pi = math::pi(); + + if (is_Cj_zero) + { + lon_before = lonj; + lon_behind = lonj + sign_lon_diff * pi; + return true; + } + + if (! (is_sin_beta_ok(sin_beta) + || math::equals(math::abs(sin_beta), sqrt_1_Cj_sqr)) ) + { + return false; + } + + CT const t_t0j = t / t0j; + CT const asin_t_t0j = asin(bounded(t_t0j, -c1, c1)); + CT const dLj = d_lambda(sin_beta); + lon_before = lonj + sign_lon_diff * (asin_t_t0j - asin_tj_t0j + dLj); + lon_behind = vd.lon0j + (vd.lon0j - lon_before); + + return true; + } + + + CT lon(CT const& delta_lon) const + { + return lonj + delta_lon; + } + + CT lat(CT const& t) const + { + // t = tan(beta) = (1-f)tan(lat) + return atan(t / one_minus_f); + } + + void vertex(CT & lon, CT & lat) const + { + lon = get_vertex_data().lon0j; + if (! is_Cj_zero) + { + lat = sjoberg_geodesic::lat(t0j); + } + else + { + CT const c2 = 2; + lat = math::pi() / c2; + } + } + + CT lon_of_equator_intersection() const + { + CT const c0 = 0; + CT const dLj = d_lambda(c0); + CT const asin_tj_t0j = asin(Cj * tan_betaj / sqrt_1_Cj_sqr); + return lonj - asin_tj_t0j + dLj; + } + + CT d_lambda(CT const& sin_beta) const + { + return sjoberg_d_lambda_e_sqr(sin_betaj, sin_beta, Cj, sqrt_1_Cj_sqr, e_sqr); + } + + // [Sjoberg12] + /*CT d_lambda(CT const& sin_beta1, CT const& sin_beta2) const + { + return sjoberg_d_lambda_e_sqr(sin_beta1, sin_beta2, Cj, sqrt_1_Cj_sqr, e_sqr); + }*/ + + CT lonj; + CT latj; + CT alphaj; + + CT one_minus_f; + CT e_sqr; + + CT tan_latj; + CT tan_betaj; + CT betaj; + CT sin_betaj; + CT cos_betaj; + CT sin_alphaj; + CT Cj; + CT Cj_sqr; + CT sqrt_1_Cj_sqr; + + int sign_lon_diff; + + bool is_on_equator; + bool is_Cj_zero; + + CT t0j; + CT asin_tj_t0j; +}; + + /*! \brief The intersection of two geodesics as proposed by Sjoberg. -\author See +\see See - [Sjoberg02] Lars E. Sjoberg, Intersections on the sphere and ellipsoid, 2002 http://link.springer.com/article/10.1007/s00190-001-0230-9 - [Sjoberg07] Lars E. Sjoberg, Geodetic intersection on the ellipsoid, 2007 http://link.springer.com/article/10.1007/s00190-007-0204-7 + - [Sjoberg12] Lars E. Sjoberg, Solutions to the ellipsoidal Clairaut constant + and the inverse geodetic problem by numerical integration, 2012 + https://www.degruyter.com/view/j/jogs.2012.2.issue-3/v10156-011-0037-4/v10156-011-0037-4.xml */ template < @@ -42,9 +629,14 @@ template > class sjoberg_intersection { + typedef sjoberg_geodesic geodesic_type; typedef Inverse inverse_type; typedef typename inverse_type::result_type inverse_result; + static bool const enable_02 = true; + static int const max_iterations_02 = 10; + static int const max_iterations_07 = 20; + public: template static inline bool apply(T1 const& lona1, T1 const& lata1, @@ -63,320 +655,536 @@ class sjoberg_intersection CT const lon_b2 = lonb2; CT const lat_b2 = latb2; - CT const alpha1 = inverse_type::apply(lon_a1, lat_a1, lon_a2, lat_a2, spheroid).azimuth; - CT const alpha2 = inverse_type::apply(lon_b1, lat_b1, lon_b2, lat_b2, spheroid).azimuth; + inverse_result const res1 = inverse_type::apply(lon_a1, lat_a1, lon_a2, lat_a2, spheroid); + inverse_result const res2 = inverse_type::apply(lon_b1, lat_b1, lon_b2, lat_b2, spheroid); - return apply(lon_a1, lat_a1, alpha1, lon_b1, lat_b1, alpha2, lon, lat, spheroid); + return apply(lon_a1, lat_a1, lon_a2, lat_a2, res1.azimuth, + lon_b1, lat_b1, lon_b2, lat_b2, res2.azimuth, + lon, lat, spheroid); } - + template - static inline bool apply(CT const& lon1, CT const& lat1, CT const& alpha1, - CT const& lon2, CT const& lat2, CT const& alpha2, + static inline bool apply(CT const& lon_a1, CT const& lat_a1, CT const& lon_a2, CT const& lat_a2, CT const& alpha_a1, + CT const& lon_b1, CT const& lat_b1, CT const& lon_b2, CT const& lat_b2, CT const& alpha_b1, CT & lon, CT & lat, Spheroid const& spheroid) { // coordinates in radians - // TODO - handle special cases like degenerated segments, equator, poles, etc. - CT const c0 = 0; CT const c1 = 1; - CT const c2 = 2; - CT const pi = math::pi(); - CT const pi_half = pi / c2; - CT const f = detail::flattening(spheroid); + CT const f = formula::flattening(spheroid); CT const one_minus_f = c1 - f; - CT const e_sqr = f * (c2 - f); - - CT const sin_alpha1 = sin(alpha1); - CT const sin_alpha2 = sin(alpha2); - - CT const tan_beta1 = one_minus_f * tan(lat1); - CT const tan_beta2 = one_minus_f * tan(lat2); - CT const beta1 = atan(tan_beta1); - CT const beta2 = atan(tan_beta2); - CT const cos_beta1 = cos(beta1); - CT const cos_beta2 = cos(beta2); - CT const sin_beta1 = sin(beta1); - CT const sin_beta2 = sin(beta2); - - // Clairaut constants (lower-case in the paper) - int const sign_C1 = math::abs(alpha1) <= pi_half ? 1 : -1; - int const sign_C2 = math::abs(alpha2) <= pi_half ? 1 : -1; - // Cj = 1 if on equator - CT const C1 = sign_C1 * cos_beta1 * sin_alpha1; - CT const C2 = sign_C2 * cos_beta2 * sin_alpha2; - - CT const sqrt_1_C1_sqr = math::sqrt(c1 - math::sqr(C1)); - CT const sqrt_1_C2_sqr = math::sqrt(c1 - math::sqr(C2)); - - // handle special case: segments on the equator - bool const on_equator1 = math::equals(sqrt_1_C1_sqr, c0); - bool const on_equator2 = math::equals(sqrt_1_C2_sqr, c0); - if (on_equator1 && on_equator2) + + geodesic_type geod1(lon_a1, lat_a1, alpha_a1, f); + geodesic_type geod2(lon_b1, lat_b1, alpha_b1, f); + + // Cj = 1 if on equator <=> sqrt_1_Cj_sqr = 0 + // Cj = 0 if vertical <=> sqrt_1_Cj_sqr = 1 + + if (geod1.is_on_equator && geod2.is_on_equator) { return false; } - else if (on_equator1) + else if (geod1.is_on_equator) { - CT const dL2 = d_lambda_e_sqr(sin_beta2, c0, C2, sqrt_1_C2_sqr, e_sqr); - CT const asin_t2_t02 = asin(C2 * tan_beta2 / sqrt_1_C2_sqr); + lon = geod2.lon_of_equator_intersection(); lat = c0; - lon = lon2 - asin_t2_t02 + dL2; return true; } - else if (on_equator2) + else if (geod2.is_on_equator) { - CT const dL1 = d_lambda_e_sqr(sin_beta1, c0, C1, sqrt_1_C1_sqr, e_sqr); - CT const asin_t1_t01 = asin(C1 * tan_beta1 / sqrt_1_C1_sqr); + lon = geod1.lon_of_equator_intersection(); lat = c0; - lon = lon1 - asin_t1_t01 + dL1; return true; } - CT const t01 = sqrt_1_C1_sqr / C1; - CT const t02 = sqrt_1_C2_sqr / C2; + // vertical segments + if (geod1.is_Cj_zero && geod2.is_Cj_zero) + { + //TODO: the geodesics may be parallel or intersect at one of the poles + return false; + } + + // (lon1 - lon2) normalized to (-180, 180] + CT const lon1_minus_lon2 = math::longitude_distance_signed(geod2.lonj, geod1.lonj); - CT const asin_t1_t01 = asin(tan_beta1 / t01); - CT const asin_t2_t02 = asin(tan_beta2 / t02); - CT const t01_t02 = t01 * t02; - CT const t01_t02_2 = c2 * t01_t02; - CT const sqr_t01_sqr_t02 = math::sqr(t01) + math::sqr(t02); + CT lon_sph = lon; - CT t = tan_beta1; - int t_id = 0; + // Starting tan(beta) + CT t = 0; - // find the initial t using simplified spherical solution - // though not entirely since the reduced latitudes and azimuths are spheroidal - // [Sjoberg07] - CT const k_base = lon1 - lon2 + asin_t2_t02 - asin_t1_t01; - + /*if (geod1.is_Cj_zero) { - CT const K = sin(k_base); - CT const d1 = sqr_t01_sqr_t02; - //CT const d2 = t01_t02_2 * math::sqrt(c1 - math::sqr(K)); - CT const d2 = t01_t02_2 * cos(k_base); - CT const D1 = math::sqrt(d1 - d2); - CT const D2 = math::sqrt(d1 + d2); - CT const K_t01_t02 = K * t01_t02; - - CT const T1 = K_t01_t02 / D1; - CT const T2 = K_t01_t02 / D2; - CT asin_T1_t01 = 0; - CT asin_T1_t02 = 0; - CT asin_T2_t01 = 0; - CT asin_T2_t02 = 0; - - // test 4 possible results - CT l1 = 0, l2 = 0, dl = 0; - bool found = check_t<0>( T1, - lon1, asin_T1_t01 = asin(T1 / t01), asin_t1_t01, - lon2, asin_T1_t02 = asin(T1 / t02), asin_t2_t02, - t, l1, l2, dl, t_id) - || check_t<1>(-T1, - lon1, -asin_T1_t01 , asin_t1_t01, - lon2, -asin_T1_t02 , asin_t2_t02, - t, l1, l2, dl, t_id) - || check_t<2>( T2, - lon1, asin_T2_t01 = asin(T2 / t01), asin_t1_t01, - lon2, asin_T2_t02 = asin(T2 / t02), asin_t2_t02, - t, l1, l2, dl, t_id) - || check_t<3>(-T2, - lon1, -asin_T2_t01 , asin_t1_t01, - lon2, -asin_T2_t02 , asin_t2_t02, - t, l1, l2, dl, t_id); - - boost::ignore_unused(found); + CT const k_base = lon1_minus_lon2 + geod2.sign_lon_diff * geod2.asin_tj_t0j; + t = sin(k_base) * geod2.t0j; + lon_sph = vertical_intersection_longitude(geod1.lonj, lon_b1, lon_b2); } + else if (geod2.is_Cj_zero) + { + CT const k_base = lon1_minus_lon2 - geod1.sign_lon_diff * geod1.asin_tj_t0j; + t = sin(-k_base) * geod1.t0j; + lon_sph = vertical_intersection_longitude(geod2.lonj, lon_a1, lon_a2); + } + else*/ + { + // TODO: Consider using betas instead of latitudes. + // Some function calls might be saved this way. + CT tan_lat_sph = 0; + sjoberg_intersection_spherical_02::apply_alt(lon_a1, lat_a1, lon_a2, lat_a2, + lon_b1, lat_b1, lon_b2, lat_b2, + lon_sph, tan_lat_sph); + t = one_minus_f * tan_lat_sph; // tan(beta) + } + + // TODO: no need to calculate atan here if reduced latitudes were used + // instead of latitudes above, in sjoberg_intersection_spherical_02 + CT const beta = atan(t); + + if (enable_02 && newton_method(geod1, geod2, beta, t, lon1_minus_lon2, lon, lat)) + { + return true; + } + + return converge_07(geod1, geod2, beta, t, lon1_minus_lon2, lon_sph, lon, lat); + } + +private: + static inline bool newton_method(geodesic_type const& geod1, geodesic_type const& geod2, // in + CT beta, CT t, CT const& lon1_minus_lon2, // in + CT & lon, CT & lat) // out + { + CT const c0 = 0; + CT const c1 = 1; + + CT const e_sqr = geod1.e_sqr; - // [Sjoberg07] - //int const d2_sign = t_id < 2 ? -1 : 1; - int const t_sign = (t_id % 2) ? -1 : 1; - // [Sjoberg02] - CT const C1_sqr = math::sqr(C1); - CT const C2_sqr = math::sqr(C2); - - CT beta = atan(t); - CT dL1 = 0, dL2 = 0; - CT asin_t_t01 = 0; - CT asin_t_t02 = 0; + CT lon1_diff = 0; + CT lon2_diff = 0; - for (int i = 0; i < 10; ++i) + CT abs_dbeta_last = 0; + + // [Sjoberg02] converges faster than solution in [Sjoberg07] + // Newton-Raphson method + for (int i = 0; i < max_iterations_02; ++i) { CT const sin_beta = sin(beta); - - // integrals approximation - dL1 = d_lambda_e_sqr(sin_beta1, sin_beta, C1, sqrt_1_C1_sqr, e_sqr); - dL2 = d_lambda_e_sqr(sin_beta2, sin_beta, C2, sqrt_1_C2_sqr, e_sqr); - - // [Sjoberg07] - /*CT const k = k_base + dL1 - dL2; - CT const K = sin(k); - CT const d1 = sqr_t01_sqr_t02; - //CT const d2 = t01_t02_2 * math::sqrt(c1 - math::sqr(K)); - CT const d2 = t01_t02_2 * cos(k); - CT const D = math::sqrt(d1 + d2_sign * d2); - CT const t_new = t_sign * K * t01_t02 / D; - CT const dt = math::abs(t_new - t); - t = t_new; - CT const new_beta = atan(t); - CT const dbeta = math::abs(new_beta - beta); - beta = new_beta;*/ - - // [Sjoberg02] - it converges faster - // Newton–Raphson method - asin_t_t01 = asin(t / t01); - asin_t_t02 = asin(t / t02); - CT const R1 = asin_t_t01 + dL1; - CT const R2 = asin_t_t02 + dL2; CT const cos_beta = cos(beta); CT const cos_beta_sqr = math::sqr(cos_beta); CT const G = c1 - e_sqr * cos_beta_sqr; - CT const f1 = C1 / cos_beta * math::sqrt(G / (cos_beta_sqr - C1_sqr)); - CT const f2 = C2 / cos_beta * math::sqrt(G / (cos_beta_sqr - C2_sqr)); - CT const abs_f1 = math::abs(f1); - CT const abs_f2 = math::abs(f2); - CT const dbeta = t_sign * (k_base - R2 + R1) / (abs_f1 + abs_f2); - - if (math::equals(dbeta, CT(0))) + + CT f1 = 0; + CT f2 = 0; + + if (!geod1.is_Cj_zero) + { + bool is_beta_ok = geod1.lon_diff(sin_beta, t, lon1_diff); + + if (is_beta_ok) + { + CT const H = cos_beta_sqr - geod1.Cj_sqr; + f1 = geod1.Cj / cos_beta * math::sqrt(G / H); + } + else + { + return false; + } + } + + if (!geod2.is_Cj_zero) + { + bool is_beta_ok = geod2.lon_diff(sin_beta, t, lon2_diff); + + if (is_beta_ok) + { + CT const H = cos_beta_sqr - geod2.Cj_sqr; + f2 = geod2.Cj / cos_beta * math::sqrt(G / H); + } + else + { + return false; + } + } + + // NOTE: Things may go wrong if the IP is near the vertex + // 1. May converge into the wrong direction (from the other way around). + // This happens when the starting point is on the other side than the vertex + // 2. During converging may "jump" into the other side of the vertex. + // In this case sin_beta/sqrt_1_Cj_sqr and t/t0j is not in [-1, 1] + // 3. f1-f2 may be 0 which means that the intermediate point is on the vertex + // In this case it's not possible to check if this is the correct result + + CT const dbeta_denom = f1 - f2; + //CT const dbeta_denom = math::abs(f1) + math::abs(f2); + + if (math::equals(dbeta_denom, c0)) + { + return false; + } + + // The sign of dbeta is changed WRT [Sjoberg02] + CT const dbeta = (lon1_minus_lon2 + lon1_diff - lon2_diff) / dbeta_denom; + + CT const abs_dbeta = math::abs(dbeta); + if (i > 0 && abs_dbeta > abs_dbeta_last) { + // The algorithm is not converging + // The intersection may be on the other side of the vertex + return false; + } + abs_dbeta_last = abs_dbeta; + + if (math::equals(dbeta, c0)) + { + // Result found break; } + // Because the sign of dbeta is changed WRT [Sjoberg02] dbeta is subtracted here beta = beta - dbeta; + t = tan(beta); } - - // t = tan(beta) = (1-f)tan(lat) - lat = atan(t / one_minus_f); - CT const l1 = lon1 + asin_t_t01 - asin_t1_t01 + dL1; - //CT const l2 = lon2 + asin_t_t02 - asin_t2_t02 + dL2; - lon = l1; + lat = geod1.lat(t); + // NOTE: if Cj is 0 then the result is lonj or lonj+180 + lon = ! geod1.is_Cj_zero + ? geod1.lon(lon1_diff) + : geod2.lon(lon2_diff); return true; } -private: - /*! Approximation of dLambda_j [Sjoberg07], expanded into taylor series in e^2 - Maxima script: - dLI_j(c_j, sinB_j, sinB) := integrate(1 / (sqrt(1 - c_j ^ 2 - x ^ 2)*(1 + sqrt(1 - e2*(1 - x ^ 2)))), x, sinB_j, sinB); - dL_j(c_j, B_j, B) := -e2 * c_j * dLI_j(c_j, B_j, B); - S: taylor(dLI_j(c_j, sinB_j, sinB), e2, 0, 3); - assume(c_j < 1); - assume(c_j > 0); - L1: factor(integrate(sqrt(-x ^ 2 - c_j ^ 2 + 1) / (x ^ 2 + c_j ^ 2 - 1), x)); - L2: factor(integrate(((x ^ 2 - 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x)); - L3: factor(integrate(((x ^ 4 - 2 * x ^ 2 + 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x)); - L4: factor(integrate(((x ^ 6 - 3 * x ^ 4 + 3 * x ^ 2 - 1)*sqrt(-x ^ 2 - c_j ^ 2 + 1)) / (x ^ 2 + c_j ^ 2 - 1), x)); - */ - static inline CT d_lambda_e_sqr(CT const& sin_betaj, CT const& sin_beta, - CT const& Cj, CT const& sqrt_1_Cj_sqr, - CT const& e_sqr) - { - if (Order == 0) - { - return 0; - } + struct geodesics_type + { + geodesics_type(geodesic_type const& g1, geodesic_type const& g2) + : geod1(g1) + , geod2(g2) + , vertex1(geod1.get_vertex_data()) + , vertex2(geod2.get_vertex_data()) + {} - CT const c2 = 2; - - CT const asin_B = asin(sin_beta / sqrt_1_Cj_sqr); - CT const asin_Bj = asin(sin_betaj / sqrt_1_Cj_sqr); - CT const L0 = (asin_B - asin_Bj) / c2; + geodesic_type const& geod1; + geodesic_type const& geod2; + typename geodesic_type::vertex_data vertex1; + typename geodesic_type::vertex_data vertex2; + }; + + struct converge_07_result + { + converge_07_result() + : lon1(0), lon2(0), k1_diff(0), k2_diff(0), t1(0), t2(0) + {} + + CT lon1, lon2; + CT k1_diff, k2_diff; + CT t1, t2; + }; - if (Order == 1) + static inline bool converge_07(geodesic_type const& geod1, geodesic_type const& geod2, + CT beta, CT t, + CT const& lon1_minus_lon2, CT const& lon_sph, + CT & lon, CT & lat) + { + //CT const c0 = 0; + //CT const c1 = 1; + //CT const c2 = 2; + //CT const pi = math::pi(); + + geodesics_type geodesics(geod1, geod2); + converge_07_result result; + + // calculate first pair of longitudes + if (!converge_07_step_one(CT(sin(beta)), t, lon1_minus_lon2, geodesics, lon_sph, result, false)) { - return -Cj * e_sqr * L0; + return false; } - CT const c1 = 1; - CT const c16 = 16; + int t_direction = 0; - CT const X = sin_beta; - CT const Xj = sin_betaj; - CT const Cj_sqr = math::sqr(Cj); - CT const Cj_sqr_plus_one = Cj_sqr + c1; - CT const one_minus_Cj_sqr = c1 - Cj_sqr; - CT const sqrt_Y = math::sqrt(-math::sqr(X) + one_minus_Cj_sqr); - CT const sqrt_Yj = math::sqrt(-math::sqr(Xj) + one_minus_Cj_sqr); - CT const L1 = (Cj_sqr_plus_one * (asin_B - asin_Bj) + X * sqrt_Y - Xj * sqrt_Yj) / c16; + CT lon_diff_prev = math::longitude_difference(result.lon1, result.lon2); - if (Order == 2) + // [Sjoberg07] + for (int i = 2; i < max_iterations_07; ++i) { - return -Cj * e_sqr * (L0 + e_sqr * L1); + // pick t candidates from previous result based on dir + CT t_cand1 = result.t1; + CT t_cand2 = result.t2; + // if direction is 0 the closer one is the first + if (t_direction < 0) + { + t_cand1 = (std::min)(result.t1, result.t2); + t_cand2 = (std::max)(result.t1, result.t2); + } + else if (t_direction > 0) + { + t_cand1 = (std::max)(result.t1, result.t2); + t_cand2 = (std::min)(result.t1, result.t2); + } + else + { + t_direction = t_cand1 < t_cand2 ? -1 : 1; + } + + CT t1 = t; + CT beta1 = beta; + // check if the further calculation is needed + if (converge_07_update(t1, beta1, t_cand1)) + { + break; + } + + bool try_t2 = false; + converge_07_result result_curr; + if (converge_07_step_one(CT(sin(beta1)), t1, lon1_minus_lon2, geodesics, lon_sph, result_curr)) + { + CT const lon_diff1 = math::longitude_difference(result_curr.lon1, result_curr.lon2); + if (lon_diff_prev > lon_diff1) + { + t = t1; + beta = beta1; + lon_diff_prev = lon_diff1; + result = result_curr; + } + else if (t_cand1 != t_cand2) + { + try_t2 = true; + } + else + { + // the result is not fully correct but it won't be more accurate + break; + } + } + // ! converge_07_step_one + else + { + if (t_cand1 != t_cand2) + { + try_t2 = true; + } + else + { + return false; + } + } + + + if (try_t2) + { + CT t2 = t; + CT beta2 = beta; + // check if the further calculation is needed + if (converge_07_update(t2, beta2, t_cand2)) + { + break; + } + + if (! converge_07_step_one(CT(sin(beta2)), t2, lon1_minus_lon2, geodesics, lon_sph, result_curr)) + { + return false; + } + + CT const lon_diff2 = math::longitude_difference(result_curr.lon1, result_curr.lon2); + if (lon_diff_prev > lon_diff2) + { + t_direction *= -1; + t = t2; + beta = beta2; + lon_diff_prev = lon_diff2; + result = result_curr; + } + else + { + // the result is not fully correct but it won't be more accurate + break; + } + } } - CT const c3 = 3; - CT const c5 = 5; - CT const c128 = 128; + lat = geod1.lat(t); + lon = ! geod1.is_Cj_zero ? result.lon1 : result.lon2; + math::normalize_longitude(lon); + + return true; + } + + static inline bool converge_07_update(CT & t, CT & beta, CT const& t_new) + { + CT const c0 = 0; + + CT const beta_new = atan(t_new); + CT const dbeta = beta_new - beta; + beta = beta_new; + t = t_new; + + return math::equals(dbeta, c0); + } + + static inline CT const& pick_t(CT const& t1, CT const& t2, int direction) + { + return direction < 0 ? (std::min)(t1, t2) : (std::max)(t1, t2); + } - CT const E = Cj_sqr * (c3 * Cj_sqr + c2) + c3; - CT const X_sqr = math::sqr(X); - CT const Xj_sqr = math::sqr(Xj); - CT const F = X * (-c2 * X_sqr + c3 * Cj_sqr + c5); - CT const Fj = Xj * (-c2 * Xj_sqr + c3 * Cj_sqr + c5); - CT const L2 = (E * (asin_B - asin_Bj) + F * sqrt_Y - Fj * sqrt_Yj) / c128; + static inline bool converge_07_step_one(CT const& sin_beta, + CT const& t, + CT const& lon1_minus_lon2, + geodesics_type const& geodesics, + CT const& lon_sph, + converge_07_result & result, + bool check_sin_beta = true) + { + bool ok = converge_07_one_geod(sin_beta, t, geodesics.geod1, geodesics.vertex1, lon_sph, + result.lon1, result.k1_diff, check_sin_beta) + && converge_07_one_geod(sin_beta, t, geodesics.geod2, geodesics.vertex2, lon_sph, + result.lon2, result.k2_diff, check_sin_beta); - if (Order == 3) + if (!ok) { - return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * L2)); + return false; } - CT const c8 = 8; - CT const c9 = 9; - CT const c10 = 10; - CT const c15 = 15; - CT const c24 = 24; - CT const c26 = 26; - CT const c33 = 33; - CT const c6144 = 6144; + CT const k = lon1_minus_lon2 + result.k1_diff - result.k2_diff; - CT const G = Cj_sqr * (Cj_sqr * (Cj_sqr * c15 + c9) + c9) + c15; - CT const H = -c10 * Cj_sqr - c26; - CT const I = Cj_sqr * (Cj_sqr * c15 + c24) + c33; - CT const J = X_sqr * (X * (c8 * X_sqr + H)) + X * I; - CT const Jj = Xj_sqr * (Xj * (c8 * Xj_sqr + H)) + Xj * I; - CT const L3 = (G * (asin_B - asin_Bj) + J * sqrt_Y - Jj * sqrt_Yj) / c6144; + // get 2 possible ts one lesser and one greater than t + // t1 is the closer one + calc_ts(t, k, geodesics.geod1, geodesics.geod2, result.t1, result.t2); - // Order 4 and higher - return -Cj * e_sqr * (L0 + e_sqr * (L1 + e_sqr * (L2 + e_sqr * L3))); + return true; } - static inline CT fj(CT const& cos_beta, CT const& cos2_beta, CT const& Cj, CT const& e_sqr) + static inline bool converge_07_one_geod(CT const& sin_beta, CT const& t, + geodesic_type const& geod, + typename geodesic_type::vertex_data const& vertex, + CT const& lon_sph, + CT & lon, CT & k_diff, + bool check_sin_beta) { + using math::detail::bounded; CT const c1 = 1; - CT const Cj_sqr = math::sqr(Cj); - return Cj / cos_beta * math::sqrt((c1 - e_sqr * cos2_beta) / (cos2_beta - Cj_sqr)); + + CT k_diff_before = 0; + CT k_diff_behind = 0; + + bool is_beta_ok = geod.k_diffs(sin_beta, vertex, k_diff_before, k_diff_behind, check_sin_beta); + + if (! is_beta_ok) + { + return false; + } + + CT const asin_t_t0j = ! geod.is_Cj_zero ? asin(bounded(t / geod.t0j, -c1, c1)) : 0; + CT const sign_asin_t_t0j = geod.sign_lon_diff * asin_t_t0j; + + CT const lon_before = geod.lonj + sign_asin_t_t0j + k_diff_before; + CT const lon_behind = geod.lonj - sign_asin_t_t0j + k_diff_behind; + + CT const lon_dist_before = math::longitude_distance_signed(lon_before, lon_sph); + CT const lon_dist_behind = math::longitude_distance_signed(lon_behind, lon_sph); + if (math::abs(lon_dist_before) <= math::abs(lon_dist_behind)) + { + k_diff = k_diff_before; + lon = lon_before; + } + else + { + k_diff = k_diff_behind; + lon = lon_behind; + } + + return true; } - template - static inline bool check_t(CT const& t, - CT const& lon_a1, CT const& asin_t_t01, CT const& asin_t1_t01, - CT const& lon_b1, CT const& asin_t_t02, CT const& asin_t2_t02, - CT & current_t, CT & current_lon1, CT & current_lon2, CT & current_dlon, - int & t_id) + static inline void calc_ts(CT const& t, CT const& k, + geodesic_type const& geod1, geodesic_type const& geod2, + CT & t1, CT& t2) { - CT const lon1 = lon_a1 + asin_t_t01 - asin_t1_t01; - CT const lon2 = lon_b1 + asin_t_t02 - asin_t2_t02; + CT const c1 = 1; + CT const c2 = 2; - // TODO - true angle difference - CT const dlon = math::abs(lon2 - lon1); + CT const K = sin(k); - bool are_equal = math::equals(dlon, CT(0)); - - if ((TId == 0) || are_equal || dlon < current_dlon) + BOOST_GEOMETRY_ASSERT(!geod1.is_Cj_zero || !geod2.is_Cj_zero); + if (geod1.is_Cj_zero) + { + t1 = K * geod2.t0j; + t2 = -t1; + } + else if (geod2.is_Cj_zero) + { + t1 = -K * geod1.t0j; + t2 = -t1; + } + else + { + CT const A = math::sqr(geod1.t0j) + math::sqr(geod2.t0j); + CT const B = c2 * geod1.t0j * geod2.t0j * math::sqrt(c1 - math::sqr(K)); + + CT const K_t01_t02 = K * geod1.t0j * geod2.t0j; + CT const D1 = math::sqrt(A + B); + CT const D2 = math::sqrt(A - B); + CT const t_new1 = K_t01_t02 / D1; + CT const t_new2 = K_t01_t02 / D2; + CT const t_new3 = -t_new1; + CT const t_new4 = -t_new2; + + // Pick 2 nearest t_new, one greater and one lesser than current t + CT const abs_t_new1 = math::abs(t_new1); + CT const abs_t_new2 = math::abs(t_new2); + CT const abs_t_max = (std::max)(abs_t_new1, abs_t_new2); + t1 = -abs_t_max; // lesser + t2 = abs_t_max; // greater + if (t1 < t) + { + if (t_new1 < t && t_new1 > t1) + t1 = t_new1; + if (t_new2 < t && t_new2 > t1) + t1 = t_new2; + if (t_new3 < t && t_new3 > t1) + t1 = t_new3; + if (t_new4 < t && t_new4 > t1) + t1 = t_new4; + } + if (t2 > t) + { + if (t_new1 > t && t_new1 < t2) + t2 = t_new1; + if (t_new2 > t && t_new2 < t2) + t2 = t_new2; + if (t_new3 > t && t_new3 < t2) + t2 = t_new3; + if (t_new4 > t && t_new4 < t2) + t2 = t_new4; + } + } + + // the first one is the closer one + if (math::abs(t - t2) < math::abs(t - t1)) { - current_t = t; - current_lon1 = lon1; - current_lon2 = lon2; - current_dlon = dlon; - t_id = TId; + std::swap(t2, t1); } + } - return are_equal; + static inline CT fj(CT const& cos_beta, CT const& cos2_beta, CT const& Cj, CT const& e_sqr) + { + CT const c1 = 1; + CT const Cj_sqr = math::sqr(Cj); + return Cj / cos_beta * math::sqrt((c1 - e_sqr * cos2_beta) / (cos2_beta - Cj_sqr)); } + + /*static inline CT vertical_intersection_longitude(CT const& ip_lon, CT const& seg_lon1, CT const& seg_lon2) + { + CT const c0 = 0; + CT const lon_2 = ip_lon > c0 ? ip_lon - pi : ip_lon + pi; + + return (std::min)(math::longitude_difference(ip_lon, seg_lon1), + math::longitude_difference(ip_lon, seg_lon2)) + <= + (std::min)(math::longitude_difference(lon_2, seg_lon1), + math::longitude_difference(lon_2, seg_lon2)) + ? ip_lon : lon_2; + }*/ }; }}} // namespace boost::geometry::formula diff --git a/include/boost/geometry/formulas/spherical.hpp b/include/boost/geometry/formulas/spherical.hpp index 23269d644a..ff24c51a88 100644 --- a/include/boost/geometry/formulas/spherical.hpp +++ b/include/boost/geometry/formulas/spherical.hpp @@ -40,38 +40,55 @@ struct result_spherical T reverse_azimuth; }; +template +static inline void sph_to_cart3d(T const& lon, T const& lat, T & x, T & y, T & z) +{ + T const cos_lat = cos(lat); + x = cos_lat * cos(lon); + y = cos_lat * sin(lon); + z = sin(lat); +} + template static inline Point3d sph_to_cart3d(PointSph const& point_sph) { typedef typename coordinate_type::type calc_t; - Point3d res; - - calc_t lon = get_as_radian<0>(point_sph); - calc_t lat = get_as_radian<1>(point_sph); + calc_t const lon = get_as_radian<0>(point_sph); + calc_t const lat = get_as_radian<1>(point_sph); + calc_t x, y, z; + sph_to_cart3d(lon, lat, x, y, z); - calc_t const cos_lat = cos(lat); - set<0>(res, cos_lat * cos(lon)); - set<1>(res, cos_lat * sin(lon)); - set<2>(res, sin(lat)); + Point3d res; + set<0>(res, x); + set<1>(res, y); + set<2>(res, z); return res; } +template +static inline void cart3d_to_sph(T const& x, T const& y, T const& z, T & lon, T & lat) +{ + lon = atan2(y, x); + lat = asin(z); +} + template static inline PointSph cart3d_to_sph(Point3d const& point_3d) { typedef typename coordinate_type::type coord_t; typedef typename coordinate_type::type calc_t; - PointSph res; - calc_t const x = get<0>(point_3d); calc_t const y = get<1>(point_3d); calc_t const z = get<2>(point_3d); + calc_t lonr, latr; + cart3d_to_sph(x, y, z, lonr, latr); - set_from_radian<0>(res, atan2(y, x)); - set_from_radian<1>(res, asin(z)); + PointSph res; + set_from_radian<0>(res, lonr); + set_from_radian<1>(res, latr); coord_t lon = get<0>(res); coord_t lat = get<1>(res); @@ -104,9 +121,9 @@ static inline int sph_side_value(Point3d1 const& norm, Point3d2 const& pt) template static inline result_spherical spherical_azimuth(T1 const& lon1, - T1 const& lat1, - T2 const& lon2, - T2 const& lat2) + T1 const& lat1, + T2 const& lon2, + T2 const& lat2) { typedef result_spherical result_type; result_type result; @@ -114,12 +131,6 @@ static inline result_spherical spherical_azimuth(T1 const& lon1, // http://williams.best.vwh.net/avform.htm#Crs // https://en.wikipedia.org/wiki/Great-circle_navigation CT dlon = lon2 - lon1; - CT cos_dlon = cos(dlon); - CT sin_dlon = sin(dlon); - CT cos_lat1 = cos(lat1); - CT cos_lat2 = cos(lat2); - CT sin_lat1 = sin(lat1); - CT sin_lat2 = sin(lat2); // An optimization which should kick in often for Boxes //if ( math::equals(dlon, ReturnType(0)) ) @@ -128,21 +139,84 @@ static inline result_spherical spherical_azimuth(T1 const& lon1, // return - sin(get_as_radian<1>(p1)) * cos_p2lat); //} - // "An alternative formula, not requiring the pre-computation of d" - // In the formula below dlon is used as "d" - result.azimuth = atan2(sin_dlon * cos_lat2, - cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon); + CT const cos_dlon = cos(dlon); + CT const sin_dlon = sin(dlon); + CT const cos_lat1 = cos(lat1); + CT const cos_lat2 = cos(lat2); + CT const sin_lat1 = sin(lat1); + CT const sin_lat2 = sin(lat2); + + { + // "An alternative formula, not requiring the pre-computation of d" + // In the formula below dlon is used as "d" + CT const y = sin_dlon * cos_lat2; + CT const x = cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon; + result.azimuth = atan2(y, x); + } if (ReverseAzimuth) { - result.reverse_azimuth = - atan2(sin_dlon * cos_lat1, - sin_lat2 * cos_lat1 * cos_dlon - cos_lat2 * sin_lat1); + CT const y = sin_dlon * cos_lat1; + CT const x = sin_lat2 * cos_lat1 * cos_dlon - cos_lat2 * sin_lat1; + result.reverse_azimuth = atan2(y, x); } return result; } +template +inline ReturnType spherical_azimuth(T1 const& lon1, T1 const& lat1, + T2 const& lon2, T2 const& lat2) +{ + return spherical_azimuth(lon1, lat1, lon2, lat2).azimuth; +} + +template +inline T spherical_azimuth(T const& lon1, T const& lat1, T const& lon2, T const& lat2) +{ + return spherical_azimuth(lon1, lat1, lon2, lat2).azimuth; +} + +template +inline int azimuth_side_value(T const& azi_a1_p, T const& azi_a1_a2) +{ + T const pi = math::pi(); + T const two_pi = math::two_pi(); + + // instead of the formula from XTD + //calc_t a_diff = asin(sin(azi_a1_p - azi_a1_a2)); + + T a_diff = azi_a1_p - azi_a1_a2; + // normalize, angle in [-pi, pi] + while (a_diff > pi) + a_diff -= two_pi; + while (a_diff < -pi) + a_diff += two_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 +} + } // namespace formula }} // namespace boost::geometry diff --git a/include/boost/geometry/formulas/thomas_direct.hpp b/include/boost/geometry/formulas/thomas_direct.hpp index f8a7f83943..f208167cf5 100644 --- a/include/boost/geometry/formulas/thomas_direct.hpp +++ b/include/boost/geometry/formulas/thomas_direct.hpp @@ -20,9 +20,8 @@ #include #include -#include - #include +#include #include @@ -82,7 +81,7 @@ class thomas_direct CT const a = CT(get_radius<0>(spheroid)); CT const b = CT(get_radius<2>(spheroid)); - CT const f = detail::flattening(spheroid); + CT const f = formula::flattening(spheroid); CT const one_minus_f = c1 - f; CT const pi = math::pi(); diff --git a/include/boost/geometry/formulas/thomas_inverse.hpp b/include/boost/geometry/formulas/thomas_inverse.hpp index d68c9de054..0853a36980 100644 --- a/include/boost/geometry/formulas/thomas_inverse.hpp +++ b/include/boost/geometry/formulas/thomas_inverse.hpp @@ -20,9 +20,8 @@ #include #include -#include - #include +#include #include @@ -78,7 +77,7 @@ class thomas_inverse CT const c4 = 4; CT const pi_half = math::pi() / c2; - CT const f = detail::flattening(spheroid); + CT const f = formula::flattening(spheroid); CT const one_minus_f = c1 - f; // CT const tan_theta1 = one_minus_f * tan(lat1); diff --git a/include/boost/geometry/formulas/vertex_latitude.hpp b/include/boost/geometry/formulas/vertex_latitude.hpp index 1d4c24a821..755067b08d 100644 --- a/include/boost/geometry/formulas/vertex_latitude.hpp +++ b/include/boost/geometry/formulas/vertex_latitude.hpp @@ -12,9 +12,9 @@ #ifndef BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP #define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP -#include -#include #include +#include +#include #include namespace boost { namespace geometry { namespace formula @@ -55,7 +55,7 @@ class vertex_latitude_on_spheroid T2 const& alp1, Spheroid const& spheroid) { - CT const f = detail::flattening(spheroid); + CT const f = formula::flattening(spheroid); CT const e2 = f * (CT(2) - f); CT const sin_alp1 = sin(alp1); @@ -76,7 +76,7 @@ class vertex_latitude_on_spheroid T2 const& alp1, Spheroid const& spheroid) { - CT const f = detail::flattening(spheroid); + CT const f = formula::flattening(spheroid); CT const one_minus_f = (CT(1) - f); diff --git a/include/boost/geometry/formulas/vincenty_direct.hpp b/include/boost/geometry/formulas/vincenty_direct.hpp index f3647ff4e6..1697e5fb63 100644 --- a/include/boost/geometry/formulas/vincenty_direct.hpp +++ b/include/boost/geometry/formulas/vincenty_direct.hpp @@ -23,9 +23,8 @@ #include #include -#include - #include +#include #include @@ -85,7 +84,7 @@ class vincenty_direct CT const radius_a = CT(get_radius<0>(spheroid)); CT const radius_b = CT(get_radius<2>(spheroid)); - CT const flattening = geometry::detail::flattening(spheroid); + CT const flattening = formula::flattening(spheroid); CT const sin_azimuth12 = sin(azimuth12); CT const cos_azimuth12 = cos(azimuth12); diff --git a/include/boost/geometry/formulas/vincenty_inverse.hpp b/include/boost/geometry/formulas/vincenty_inverse.hpp index bbda00036b..032e16e291 100644 --- a/include/boost/geometry/formulas/vincenty_inverse.hpp +++ b/include/boost/geometry/formulas/vincenty_inverse.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, 2016, 2017. +// Modifications copyright (c) 2014-2017 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -23,9 +23,8 @@ #include #include -#include - #include +#include #include @@ -41,7 +40,7 @@ namespace boost { namespace geometry { namespace formula \brief The solution of the inverse problem of geodesics on latlong coordinates, after Vincenty, 1975 \author See - http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf - - http://www.icsm.gov.au/gda/gdav2.3.pdf + - http://www.icsm.gov.au/gda/gda-v_2.4.pdf \author Adapted from various implementations to get it close to the original document - http://www.movable-type.co.uk/scripts/LatLongVincenty.html - http://exogen.case.edu/projects/geopy/source/geopy.distance.html @@ -99,10 +98,10 @@ struct vincenty_inverse CT const radius_a = CT(get_radius<0>(spheroid)); CT const radius_b = CT(get_radius<2>(spheroid)); - CT const flattening = geometry::detail::flattening(spheroid); + CT const f = formula::flattening(spheroid); // U: reduced latitude, defined by tan U = (1-f) tan phi - CT const one_min_f = c1 - flattening; + CT const one_min_f = c1 - f; CT const tan_U1 = one_min_f * tan(lat1); // above (1) CT const tan_U2 = one_min_f * tan(lat2); // above (1) @@ -113,8 +112,9 @@ struct vincenty_inverse CT const cos_U1 = c1 / temp_den_U1; CT const cos_U2 = c1 / temp_den_U2; // sin = tan / sqrt(1 + tan^2) - CT const sin_U1 = tan_U1 / temp_den_U1; - CT const sin_U2 = tan_U2 / temp_den_U2; + // sin = tan * cos + CT const sin_U1 = tan_U1 * cos_U1; + CT const sin_U2 = tan_U2 * cos_U2; // calculate sin U and cos U directly //CT const U1 = atan(tan_U1); @@ -130,7 +130,8 @@ struct vincenty_inverse CT sin_sigma; CT sin_alpha; CT cos2_alpha; - CT cos2_sigma_m; + CT cos_2sigma_m; + CT cos2_2sigma_m; CT sigma; int counter = 0; // robustness @@ -144,12 +145,13 @@ struct vincenty_inverse CT cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15) sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17) cos2_alpha = c1 - math::sqr(sin_alpha); - cos2_sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18) + cos_2sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18) + cos2_2sigma_m = math::sqr(cos_2sigma_m); - CT C = flattening/c16 * cos2_alpha * (c4 + flattening * (c4 - c3 * cos2_alpha)); // (10) + CT C = f/c16 * cos2_alpha * (c4 + f * (c4 - c3 * cos2_alpha)); // (10) sigma = atan2(sin_sigma, cos_sigma); // (16) - lambda = L + (c1 - C) * flattening * sin_alpha * - (sigma + C * sin_sigma * ( cos2_sigma_m + C * cos_sigma * (-c1 + c2 * math::sqr(cos2_sigma_m)))); // (11) + lambda = L + (c1 - C) * f * sin_alpha * + (sigma + C * sin_sigma * (cos_2sigma_m + C * cos_sigma * (-c1 + c2 * cos2_2sigma_m))); // (11) ++counter; // robustness @@ -182,8 +184,10 @@ struct vincenty_inverse CT A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3) CT B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4) - CT delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m) - - (B/c6) * cos2_sigma_m * (-c3 + c4 * math::sqr(sin_sigma)) * (-c3 + c4 * cos2_sigma_m))); // (6) + CT const cos_sigma = cos(sigma); + CT const sin2_sigma = math::sqr(sin_sigma); + CT delta_sigma = B * sin_sigma * (cos_2sigma_m + (B/c4) * (cos_sigma* (-c1 + c2 * cos2_2sigma_m) + - (B/c6) * cos_2sigma_m * (-c3 + c4 * sin2_sigma) * (-c3 + c4 * cos2_2sigma_m))); // (6) result.distance = radius_b * A * (sigma - delta_sigma); // (19) } @@ -206,7 +210,7 @@ struct vincenty_inverse typedef differential_quantities quantities; quantities::apply(lon1, lat1, lon2, lat2, result.azimuth, result.reverse_azimuth, - radius_b, flattening, + radius_b, f, result.reduced_length, result.geodesic_scale); } diff --git a/include/boost/geometry/index/parameters.hpp b/include/boost/geometry/index/parameters.hpp index 6a3a67c6cc..1a9469c103 100644 --- a/include/boost/geometry/index/parameters.hpp +++ b/include/boost/geometry/index/parameters.hpp @@ -2,7 +2,7 @@ // // R-tree algorithms parameters // -// Copyright (c) 2011-2013 Adam Wulkiewicz, Lodz, Poland. +// Copyright (c) 2011-2017 Adam Wulkiewicz, Lodz, Poland. // // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at @@ -162,8 +162,8 @@ class dynamic_linear \param max_elements Maximum number of elements in nodes. \param min_elements Minimum number of elements in nodes. Default: 0.3*Max. */ - dynamic_linear(size_t max_elements, - size_t min_elements = detail::default_min_elements_d()) + explicit dynamic_linear(size_t max_elements, + size_t min_elements = detail::default_min_elements_d()) : m_max_elements(max_elements) , m_min_elements(detail::default_min_elements_d_calc(max_elements, min_elements)) { @@ -191,8 +191,8 @@ class dynamic_quadratic \param max_elements Maximum number of elements in nodes. \param min_elements Minimum number of elements in nodes. Default: 0.3*Max. */ - dynamic_quadratic(size_t max_elements, - size_t min_elements = detail::default_min_elements_d()) + explicit dynamic_quadratic(size_t max_elements, + size_t min_elements = detail::default_min_elements_d()) : m_max_elements(max_elements) , m_min_elements(detail::default_min_elements_d_calc(max_elements, min_elements)) { @@ -228,10 +228,10 @@ class dynamic_rstar nearly minimum overlap cost, otherwise all leafs are analyzed and true minimum overlap cost is calculated. Default: 32. */ - dynamic_rstar(size_t max_elements, - size_t min_elements = detail::default_min_elements_d(), - size_t reinserted_elements = detail::default_rstar_reinserted_elements_d(), - size_t overlap_cost_threshold = 32) + explicit dynamic_rstar(size_t max_elements, + size_t min_elements = detail::default_min_elements_d(), + size_t reinserted_elements = detail::default_rstar_reinserted_elements_d(), + size_t overlap_cost_threshold = 32) : m_max_elements(max_elements) , m_min_elements(detail::default_min_elements_d_calc(max_elements, min_elements)) , m_reinserted_elements(detail::default_rstar_reinserted_elements_d_calc(max_elements, reinserted_elements)) diff --git a/include/boost/geometry/strategies/azimuth.hpp b/include/boost/geometry/strategies/azimuth.hpp index 470882d6b2..8044d3933b 100644 --- a/include/boost/geometry/strategies/azimuth.hpp +++ b/include/boost/geometry/strategies/azimuth.hpp @@ -1,7 +1,8 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2016 Oracle and/or its affiliates. +// 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 @@ -20,17 +21,18 @@ namespace strategy { namespace azimuth { namespace services { /*! - \brief Traits class binding a default azimuth strategy to a coordinate system - \ingroup azimuth - \tparam Tag tag of coordinate system +\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 +template struct default_strategy { BOOST_MPL_ASSERT_MSG ( - false, NOT_IMPLEMENTED_FOR_THIS_CALCULATION_TYPE - , (types) + false, NOT_IMPLEMENTED_FOR_THIS_TYPE + , (types) ); }; diff --git a/include/boost/geometry/strategies/cartesian/area_surveyor.hpp b/include/boost/geometry/strategies/cartesian/area_surveyor.hpp index ba76f67946..b3f19b1b1e 100644 --- a/include/boost/geometry/strategies/cartesian/area_surveyor.hpp +++ b/include/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 #include #include +#include #include diff --git a/include/boost/geometry/strategies/cartesian/azimuth_cartesian.hpp b/include/boost/geometry/strategies/cartesian/azimuth.hpp similarity index 71% rename from include/boost/geometry/strategies/cartesian/azimuth_cartesian.hpp rename to include/boost/geometry/strategies/cartesian/azimuth.hpp index 247c6cb333..62f804e8f5 100644 --- a/include/boost/geometry/strategies/cartesian/azimuth_cartesian.hpp +++ b/include/boost/geometry/strategies/cartesian/azimuth.hpp @@ -1,14 +1,15 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2016 Oracle and/or its affiliates. +// 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_CARTESIAN_HPP -#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AZIMUTH_CARTESIAN_HPP +#ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AZIMUTH_HPP +#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AZIMUTH_HPP #include @@ -20,7 +21,7 @@ namespace strategy { namespace azimuth template < - typename CalculationType + typename CalculationType = void > class cartesian {}; @@ -45,4 +46,4 @@ struct default_strategy }} // namespace boost::geometry -#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AZIMUTH_CARTESIAN_HPP +#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AZIMUTH_HPP diff --git a/include/boost/geometry/strategies/cartesian/envelope_segment.hpp b/include/boost/geometry/strategies/cartesian/envelope_segment.hpp index 63d1ff9bac..0ddbf12a41 100644 --- a/include/boost/geometry/strategies/cartesian/envelope_segment.hpp +++ b/include/boost/geometry/strategies/cartesian/envelope_segment.hpp @@ -2,6 +2,7 @@ // 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 @@ -10,8 +11,11 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_ENVELOPE_SEGMENT_HPP #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_ENVELOPE_SEGMENT_HPP -#include + #include +#include +#include + namespace boost { namespace geometry { @@ -19,7 +23,10 @@ namespace boost { namespace geometry namespace strategy { namespace envelope { -template +template +< + typename CalculationType = void +> class cartesian_segment { public : diff --git a/include/boost/geometry/strategies/cartesian/cart_intersect.hpp b/include/boost/geometry/strategies/cartesian/intersection.hpp similarity index 93% rename from include/boost/geometry/strategies/cartesian/cart_intersect.hpp rename to include/boost/geometry/strategies/cartesian/intersection.hpp index 1f6f56d66a..20b6b93367 100644 --- a/include/boost/geometry/strategies/cartesian/cart_intersect.hpp +++ b/include/boost/geometry/strategies/cartesian/intersection.hpp @@ -34,6 +34,8 @@ #include #include +#include +#include #include #include #include @@ -62,8 +64,11 @@ namespace strategy { namespace intersection /*! \see http://mathworld.wolfram.com/Line-LineIntersection.html */ -template -struct relate_cartesian_segments +template +< + typename CalculationType = void +> +struct cartesian_segments { typedef side::side_by_triangle side_strategy_type; @@ -95,6 +100,39 @@ struct relate_cartesian_segments return strategy_type(); } + template + struct area_strategy + { + typedef area::surveyor + < + typename point_type::type, + CalculationType + > type; + }; + + template + static inline typename area_strategy::type get_area_strategy() + { + typedef typename area_strategy::type strategy_type; + return strategy_type(); + } + + template + struct distance_strategy + { + typedef distance::pythagoras + < + CalculationType + > type; + }; + + template + static inline typename distance_strategy::type get_distance_strategy() + { + typedef typename distance_strategy::type strategy_type; + return strategy_type(); + } + template struct segment_intersection_info { @@ -588,7 +626,7 @@ namespace services template struct default_strategy { - typedef relate_cartesian_segments type; + typedef cartesian_segments type; }; } // namespace services @@ -606,25 +644,25 @@ namespace within { namespace services template struct default_strategy { - typedef strategy::intersection::relate_cartesian_segments<> type; + typedef strategy::intersection::cartesian_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_cartesian_segments<> type; + typedef strategy::intersection::cartesian_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_cartesian_segments<> type; + typedef strategy::intersection::cartesian_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_cartesian_segments<> type; + typedef strategy::intersection::cartesian_segments<> type; }; }} // within::services @@ -635,25 +673,25 @@ namespace covered_by { namespace services template struct default_strategy { - typedef strategy::intersection::relate_cartesian_segments<> type; + typedef strategy::intersection::cartesian_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_cartesian_segments<> type; + typedef strategy::intersection::cartesian_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_cartesian_segments<> type; + typedef strategy::intersection::cartesian_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_cartesian_segments<> type; + typedef strategy::intersection::cartesian_segments<> type; }; }} // within::services diff --git a/include/boost/geometry/strategies/envelope.hpp b/include/boost/geometry/strategies/envelope.hpp index f3ab823eab..fde9c858a6 100644 --- a/include/boost/geometry/strategies/envelope.hpp +++ b/include/boost/geometry/strategies/envelope.hpp @@ -1,7 +1,8 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2016 Oracle and/or its affiliates. +// 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 @@ -20,17 +21,18 @@ namespace strategy { namespace envelope { namespace services { /*! - \brief Traits class binding a default envelope strategy to a coordinate system - \ingroup envelope - \tparam Tag tag of coordinate system +\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 +template struct default_strategy { BOOST_MPL_ASSERT_MSG ( - false, NOT_IMPLEMENTED_FOR_THIS_CALCULATION_TYPE - , (types) + false, NOT_IMPLEMENTED_FOR_THIS_TYPE + , (types) ); }; diff --git a/include/boost/geometry/strategies/geographic/area_geographic.hpp b/include/boost/geometry/strategies/geographic/area.hpp similarity index 64% rename from include/boost/geometry/strategies/geographic/area_geographic.hpp rename to include/boost/geometry/strategies/geographic/area.hpp index efa79df467..e1d3b09b5a 100644 --- a/include/boost/geometry/strategies/geographic/area_geographic.hpp +++ b/include/boost/geometry/strategies/geographic/area.hpp @@ -1,7 +1,8 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2016 Oracle and/or its affiliates. +// 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 @@ -10,10 +11,17 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP + +#include + #include -#include +#include + +#include + #include + namespace boost { namespace geometry { @@ -21,25 +29,38 @@ namespace strategy { namespace area { /*! -\brief Geographic area calculation by trapezoidal rule plus integral -approximation that gives the ellipsoidal correction - - +\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, - template class Inverse = - geometry::formula::thomas_inverse, - std::size_t SeriesOrder = 2, - bool ExpandEpsN = true, - bool LongSegment = false, - typename Spheroid = void, + typename FormulaPolicy = strategy::andoyer, + std::size_t SeriesOrder = strategy::default_order::value, + typename Spheroid = srs::spheroid, 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 @@ -53,42 +74,22 @@ class geographic >::type, CalculationType >::type CT; -/* - typedef typename boost::mpl::if_c - < - boost::is_void::type::value, - typename geometry::formula::thomas_inverse - < - CT, - false, - true, - true - >, - Strategy - >::type AzimuthStrategy; -*/ - typedef typename boost::mpl::if_c - < - boost::is_void::type::value, - geometry::srs::spheroid, - Spheroid - >::type SpheroidType; protected : struct spheroid_constants { - SpheroidType m_spheroid; + 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(SpheroidType spheroid) + inline spheroid_constants(Spheroid const& spheroid) : m_spheroid(spheroid) , m_a2(math::sqr(get_radius<0>(spheroid))) - , m_e2(detail::flattening(spheroid) - * (CT(2.0) - CT(detail::flattening(spheroid)))) + , m_e2(formula::flattening(spheroid) + * (CT(2.0) - CT(formula::flattening(spheroid)))) , m_ep2(m_e2 / (CT(1.0) - m_e2)) , m_ep(math::sqrt(m_ep2)) , m_c2((m_a2 / CT(2.0)) + @@ -149,24 +150,26 @@ public : typedef PointOfSegment segment_point_type; typedef area_sums state_type; - inline geographic(SpheroidType spheroid = SpheroidType()) - : spheroid_const(spheroid) + explicit inline geographic(Spheroid const& spheroid = Spheroid()) + : m_spheroid_constants(spheroid) {} inline void apply(PointOfSegment const& p1, - PointOfSegment const& p2, - area_sums& state) const + PointOfSegment const& p2, + area_sums& state) const { if (! geometry::math::equals(get<0>(p1), get<0>(p2))) { - typedef geometry::formula::area_formulas area_formulas; + typedef geometry::formula::area_formulas + < + CT, SeriesOrder, ExpandEpsN + > area_formulas; typename area_formulas::return_type_ellipsoidal result = - area_formulas::template ellipsoidal - (p1, p2, spheroid_const); + area_formulas::template ellipsoidal + (p1, p2, m_spheroid_constants); state.m_excess_sum += result.spherical_term; state.m_correction_sum += result.ellipsoidal_term; @@ -179,11 +182,11 @@ public : inline return_type result(area_sums const& state) const { - return state.area(spheroid_const); + return state.area(m_spheroid_constants); } private: - spheroid_constants spheroid_const; + spheroid_constants m_spheroid_constants; }; diff --git a/include/boost/geometry/strategies/geographic/azimuth_geographic.hpp b/include/boost/geometry/strategies/geographic/azimuth.hpp similarity index 59% rename from include/boost/geometry/strategies/geographic/azimuth_geographic.hpp rename to include/boost/geometry/strategies/geographic/azimuth.hpp index 18f7fcf35c..47f59d1033 100644 --- a/include/boost/geometry/strategies/geographic/azimuth_geographic.hpp +++ b/include/boost/geometry/strategies/geographic/azimuth.hpp @@ -11,9 +11,15 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AZIMUTH_HPP #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AZIMUTH_HPP + +#include + #include +#include + +#include +#include -#include namespace boost { namespace geometry { @@ -23,10 +29,9 @@ namespace strategy { namespace azimuth template < - typename CalculationType, - typename Spheroid = geometry::srs::spheroid, - template class Inverse = - geometry::formula::thomas_inverse + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid, + typename CalculationType = void > class geographic { @@ -47,17 +52,20 @@ public : return m_spheroid; } - inline void apply(CalculationType const& lon1_rad, - CalculationType const& lat1_rad, - CalculationType const& lon2_rad, - CalculationType const& lat2_rad, - CalculationType& a1, - CalculationType& a2) const + template + 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 Inverse inverse_type; + typedef typename boost::mpl::if_ + < + boost::is_void, T, CalculationType + >::type calc_t; + + typedef typename FormulaPolicy::template inverse inverse_type; typedef typename inverse_type::result_type inverse_result; - inverse_result i_res = inverse_type::apply(lon1_rad, lat1_rad, - lon2_rad, lat2_rad, + 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; @@ -75,7 +83,12 @@ namespace services template struct default_strategy { - typedef strategy::azimuth::geographic type; + typedef strategy::azimuth::geographic + < + strategy::andoyer, + srs::spheroid, + CalculationType + > type; }; } diff --git a/include/boost/geometry/strategies/geographic/distance.hpp b/include/boost/geometry/strategies/geographic/distance.hpp new file mode 100644 index 0000000000..d3656f449c --- /dev/null +++ b/include/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 +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace distance +{ + +template +< + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid, + typename CalculationType = void +> +class geographic +{ +public : + template + 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 + inline typename calculation_type::type + apply(Point1 const& point1, Point2 const& point2) const + { + return FormulaPolicy::template inverse + < + typename calculation_type::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 > +{ + typedef strategy_tag_distance_point_point type; +}; + + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + typename P1, + typename P2 +> +struct return_type, P1, P2> + : geographic::template calculation_type +{}; + + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct comparable_type > +{ + typedef geographic type; +}; + + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct get_comparable > +{ + static inline geographic + apply(geographic const& input) + { + return input; + } +}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + typename P1, + typename P2 +> +struct result_from_distance, P1, P2> +{ + template + static inline typename return_type, P1, P2>::type + apply(geographic const& , T const& value) + { + return value; + } +}; + + +template +struct default_strategy +{ + typedef strategy::distance::geographic + < + strategy::andoyer, + srs::spheroid + < + typename select_coordinate_type::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/include/boost/geometry/strategies/geographic/distance_andoyer.hpp b/include/boost/geometry/strategies/geographic/distance_andoyer.hpp index 1946cd1090..d732951642 100644 --- a/include/boost/geometry/strategies/geographic/distance_andoyer.hpp +++ b/include/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 -#include -#include -#include - -#include - -#include - -#include - -#include -#include -#include +#include +#include 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, typename CalculationType = void > class andoyer + : public strategy::distance::geographic + < + strategy::andoyer, Spheroid, CalculationType + > { -public : - template - 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 - inline typename calculation_type::type - apply(Point1 const& point1, Point2 const& point2) const - { - return geometry::formula::andoyer_inverse - < - typename calculation_type::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, P1, P2> }; -template -struct default_strategy -{ - typedef strategy::distance::andoyer - < - srs::spheroid - < - typename select_coordinate_type::type - > - > type; -}; - - } // namespace services #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS @@ -177,4 +125,4 @@ struct default_strategy -#include +#include +#include -#include - -#include -#include - -#include namespace boost { namespace geometry { @@ -45,57 +39,28 @@ namespace strategy { namespace distance */ template < - typename Spheroid, + typename Spheroid = srs::spheroid, typename CalculationType = void > class thomas + : public strategy::distance::geographic + < + strategy::thomas, Spheroid, CalculationType + > { -public : - template - 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 - inline typename calculation_type::type - apply(Point1 const& point1, Point2 const& point2) const - { - return geometry::formula::thomas_inverse - < - typename calculation_type::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/include/boost/geometry/strategies/geographic/distance_vincenty.hpp b/include/boost/geometry/strategies/geographic/distance_vincenty.hpp index e79e9aeb46..41146db9ff 100644 --- a/include/boost/geometry/strategies/geographic/distance_vincenty.hpp +++ b/include/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 -#include +#include +#include -#include - -#include -#include - -#include namespace boost { namespace geometry { @@ -47,57 +41,28 @@ namespace strategy { namespace distance */ template < - typename Spheroid, + typename Spheroid = srs::spheroid, typename CalculationType = void > class vincenty + : public strategy::distance::geographic + < + strategy::vincenty, Spheroid, CalculationType + > { -public : - template - 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 - inline typename calculation_type::type - apply(Point1 const& point1, Point2 const& point2) const - { - return geometry::formula::vincenty_inverse - < - typename calculation_type::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/include/boost/geometry/strategies/geographic/envelope_segment.hpp b/include/boost/geometry/strategies/geographic/envelope_segment.hpp index 3d7b3cb9e1..3641b39428 100644 --- a/include/boost/geometry/strategies/geographic/envelope_segment.hpp +++ b/include/boost/geometry/strategies/geographic/envelope_segment.hpp @@ -11,10 +11,14 @@ #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ENVELOPE_SEGMENT_HPP #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_ENVELOPE_SEGMENT_HPP + #include #include +#include #include -#include +#include +#include + namespace boost { namespace geometry { @@ -24,10 +28,9 @@ namespace strategy { namespace envelope template < - typename CalculationType, - typename Spheroid = geometry::srs::spheroid, - template class Inverse = - geometry::formula::thomas_inverse + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = geometry::srs::spheroid, + typename CalculationType = void > class geographic_segment { @@ -50,9 +53,9 @@ class geographic_segment geometry::strategy::azimuth::geographic < - CalculationType, + FormulaPolicy, Spheroid, - Inverse + CalculationType > azimuth_geographic(m_spheroid); typedef typename coordinate_system::type::units units_type; @@ -81,7 +84,12 @@ namespace services template struct default_strategy { - typedef strategy::envelope::geographic_segment type; + typedef strategy::envelope::geographic_segment + < + strategy::andoyer, + srs::spheroid, + CalculationType + > type; }; } diff --git a/include/boost/geometry/strategies/geographic/intersection.hpp b/include/boost/geometry/strategies/geographic/intersection.hpp new file mode 100644 index 0000000000..5a0a8eb582 --- /dev/null +++ b/include/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 + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +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::value, + typename Spheroid = srs::spheroid, + typename CalculationType = void +> +struct geographic_segments +{ + typedef side::geographic + < + FormulaPolicy, Spheroid, CalculationType + > side_strategy_type; + + inline side_strategy_type get_side_strategy() + { + return side_strategy_type(m_spheroid); + } + + template + struct point_in_geometry_strategy + { + typedef strategy::within::winding + < + typename point_type::type, + typename point_type::type, + side_strategy_type, + CalculationType + > type; + }; + + template + inline typename point_in_geometry_strategy::type + get_point_in_geometry_strategy() + { + typedef typename point_in_geometry_strategy + < + Geometry1, Geometry2 + >::type strategy_type; + return strategy_type(get_side_strategy()); + } + + template + struct area_strategy + { + typedef area::geographic + < + typename point_type::type, + FormulaPolicy, + Order, + Spheroid, + CalculationType + > type; + }; + + template + inline typename area_strategy::type get_area_strategy() + { + typedef typename area_strategy::type strategy_type; + return strategy_type(m_spheroid); + } + + template + struct distance_strategy + { + typedef distance::geographic + < + FormulaPolicy, + Spheroid, + CalculationType + > type; + }; + + template + inline typename distance_strategy::type get_distance_strategy() + { + typedef typename distance_strategy::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 + 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 + void assign_a(Point& point, Segment1 const& a, Segment2 const& b) const + { + assign(point, a, b); + } + template + void assign_b(Point& point, Segment1 const& a, Segment2 const& b) const + { + assign(point, a, b); + } + + template + 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::type point1_t; + typedef typename point_type::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(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) ); + BOOST_CONCEPT_ASSERT( (concepts::ConstSegment) ); + + typedef typename select_calculation_type + ::type calc_t; + + // normalized spheroid + srs::spheroid spheroid = normalized_spheroid(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 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(a, true, b1, b2, a1, a2, res_b1_b2, res_b1_a1, is_b_reversed); + } + else if (b_is_point) + { + return collinear_one_degenerated(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 ra_from(dist_b1_a1, dist_b1_b2); + segment_ratio ra_to(dist_b1_a2, dist_b1_b2); + segment_ratio rb_from(dist_a1_b1, dist_a1_a2); + segment_ratio 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 + > 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 + 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(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 + 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(dist_a1_bi, dist_a1_a2).on_segment(); + } + + template + 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 + ::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 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(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(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 + 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 + static inline bool is_near(CalcT const& dist) + { + // NOTE: This strongly depends on the Inverse method + CalcT const small_number = CalcT(boost::is_same::value ? 0.0001 : 0.00000001); + return math::abs(dist) <= small_number; + } + + template + 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 + 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(azimuth1, azimuth2); + return math::abs(angle_diff) <= math::half_pi(); + } + + template + static inline void sides_reverse_segment(side_info & sides) + { + // names assuming segment A is reversed (Which == 0) + int a1_wrt_b = sides.template get(); + int a2_wrt_b = sides.template get(); + std::swap(a1_wrt_b, a2_wrt_b); + sides.template set(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 + static inline srs::spheroid normalized_spheroid(SpheroidT const& spheroid) + { + return srs::spheroid(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/include/boost/geometry/strategies/geographic/intersection_elliptic.hpp b/include/boost/geometry/strategies/geographic/intersection_elliptic.hpp new file mode 100644 index 0000000000..76e9940fe3 --- /dev/null +++ b/include/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 + +#include + +#include + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace intersection +{ + +template +struct great_elliptic_segments_calc_policy + : spherical_segments_calc_policy +{ + explicit great_elliptic_segments_calc_policy(Spheroid const& spheroid = Spheroid()) + : m_spheroid(spheroid) + {} + + template + Point from_cart3d(Point3d const& point_3d) const + { + return formula::cart3d_to_geo(point_3d, m_spheroid); + } + + template + Point3d to_cart3d(Point const& point) const + { + return formula::geo_to_cart3d(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 + struct plane + { + typedef typename coordinate_type::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 + plane get_plane(Point3d const& p1, Point3d const& p2) const + { + return plane(p1, p2); + } + + template + bool intersection_points(plane const& plane1, + plane const& plane2, + Point3d & ip1, Point3d & ip2) const + { + typedef typename coordinate_type::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 +struct experimental_elliptic_segments_calc_policy +{ + explicit experimental_elliptic_segments_calc_policy(Spheroid const& spheroid = Spheroid()) + : m_spheroid(spheroid) + {} + + template + Point from_cart3d(Point3d const& point_3d) const + { + return formula::cart3d_to_geo(point_3d, m_spheroid); + } + + template + Point3d to_cart3d(Point const& point) const + { + return formula::geo_to_cart3d(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 + struct plane + { + typedef typename coordinate_type::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 + plane get_plane(Point3d const& p1, Point3d const& p2) const + { + return plane(p1, p2, m_spheroid); + } + + template + bool intersection_points(plane const& plane1, + plane 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, + typename CalculationType = void +> +struct great_elliptic_segments + : ecef_segments + < + great_elliptic_segments_calc_policy, + CalculationType + > +{}; + +template +< + typename Spheroid = srs::spheroid, + typename CalculationType = void +> +struct experimental_elliptic_segments + : ecef_segments + < + experimental_elliptic_segments_calc_policy, + CalculationType + > +{}; + + +}} // namespace strategy::intersection + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INTERSECTION_ELLIPTIC_HPP diff --git a/include/boost/geometry/strategies/geographic/parameters.hpp b/include/boost/geometry/strategies/geographic/parameters.hpp new file mode 100644 index 0000000000..5638db50fa --- /dev/null +++ b/include/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 +#include +#include + +#include +#include + + +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 +struct default_order +{ + BOOST_MPL_ASSERT_MSG + ( + false, NOT_IMPLEMENTED_FOR_THIS_TYPE + , (types) + ); +}; + +template<> +struct default_order + : boost::mpl::integral_c +{}; + +template<> +struct default_order + : boost::mpl::integral_c +{}; + +template<> +struct default_order + : boost::mpl::integral_c +{}; + +}}} // namespace boost::geometry::strategy + + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_PARAMETERS_HPP diff --git a/include/boost/geometry/strategies/geographic/side.hpp b/include/boost/geometry/strategies/geographic/side.hpp new file mode 100644 index 0000000000..4cbfc6d314 --- /dev/null +++ b/include/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 +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +//#include + + +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, + typename CalculationType = void +> +class geographic +{ +public: + geographic() + {} + + explicit geographic(Spheroid const& model) + : m_model(model) + {} + + template + 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 typename FormulaPolicy::template inverse + inverse_formula; + + calc_t a1p = azimuth(p1, p, m_model); + calc_t a12 = azimuth(p1, p2, m_model); + + return formula::azimuth_side_value(a1p, a12); + } + +private: + template + 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/include/boost/geometry/strategies/geographic/side_andoyer.hpp b/include/boost/geometry/strategies/geographic/side_andoyer.hpp index c3e71cd1cd..204e45f6e2 100644 --- a/include/boost/geometry/strategies/geographic/side_andoyer.hpp +++ b/include/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 - -#include +#include 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 +template +< + typename Spheroid = srs::spheroid, + typename CalculationType = void +> class andoyer - : public detail::by_azimuth + : public side::geographic { - typedef detail::by_azimuth base_t; + typedef side::geographic base_t; public: - andoyer(Model const& model = Model()) + andoyer() + {} + + explicit andoyer(Spheroid const& model) : base_t(model) {} }; diff --git a/include/boost/geometry/strategies/geographic/side_detail.hpp b/include/boost/geometry/strategies/geographic/side_detail.hpp deleted file mode 100644 index ce1b47c88e..0000000000 --- a/include/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 -#include -#include -#include - -#include -#include -#include - -#include -//#include - - -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 class InverseFormula, - typename Model, - typename CalculationType = void> -class by_azimuth -{ -public: - by_azimuth(Model const& model = Model()) - : m_model(model) - {} - - template - 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 inverse_formula; - - calc_t a1p = azimuth(p1, p, m_model); - calc_t a12 = azimuth(p1, p2, m_model); - - calc_t const pi = math::pi(); - - // 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 - 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/include/boost/geometry/strategies/geographic/side_thomas.hpp b/include/boost/geometry/strategies/geographic/side_thomas.hpp index 96b0323307..e6f8d77b58 100644 --- a/include/boost/geometry/strategies/geographic/side_thomas.hpp +++ b/include/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 - -#include +#include 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 +template +< + typename Spheroid = srs::spheroid, + typename CalculationType = void +> class thomas - : public detail::by_azimuth + : public side::geographic { - typedef detail::by_azimuth base_t; + typedef side::geographic base_t; public: - thomas(Model const& model = Model()) + thomas() + {} + + explicit thomas(Spheroid const& model) : base_t(model) {} }; diff --git a/include/boost/geometry/strategies/geographic/side_vincenty.hpp b/include/boost/geometry/strategies/geographic/side_vincenty.hpp index 103277a8bd..b2f51b0901 100644 --- a/include/boost/geometry/strategies/geographic/side_vincenty.hpp +++ b/include/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 - -#include +#include 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 +template +< + typename Spheroid = srs::spheroid, + typename CalculationType = void +> class vincenty - : public detail::by_azimuth + : public side::geographic { - typedef detail::by_azimuth base_t; + typedef side::geographic base_t; public: - vincenty(Model const& model = Model()) + vincenty() + {} + + explicit vincenty(Spheroid const& model) : base_t(model) {} }; diff --git a/include/boost/geometry/strategies/intersection_strategies.hpp b/include/boost/geometry/strategies/intersection_strategies.hpp index 3b836ea512..a173505804 100644 --- a/include/boost/geometry/strategies/intersection_strategies.hpp +++ b/include/boost/geometry/strategies/intersection_strategies.hpp @@ -22,10 +22,10 @@ #include #include -#include #include +#include -#include +#include #include #include #include diff --git a/include/boost/geometry/strategies/side.hpp b/include/boost/geometry/strategies/side.hpp index 376f2fdf1b..9aaa2bdddc 100644 --- a/include/boost/geometry/strategies/side.hpp +++ b/include/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 +template struct default_strategy { BOOST_MPL_ASSERT_MSG ( false, NOT_IMPLEMENTED_FOR_THIS_TYPE - , (types) + , (types) ); }; diff --git a/include/boost/geometry/strategies/spherical/area_spherical.hpp b/include/boost/geometry/strategies/spherical/area.hpp similarity index 76% rename from include/boost/geometry/strategies/spherical/area_spherical.hpp rename to include/boost/geometry/strategies/spherical/area.hpp index faff0d9b45..206b734548 100644 --- a/include/boost/geometry/strategies/spherical/area_spherical.hpp +++ b/include/boost/geometry/strategies/spherical/area.hpp @@ -8,12 +8,15 @@ // 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_SPHERICAL_HPP -#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_SPHERICAL_HPP +#ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HPP +#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HPP + #include #include #include +#include + namespace boost { namespace geometry { @@ -22,19 +25,27 @@ namespace strategy { namespace area { /*! -\brief Spherical area calculation by trapezoidal rule - +\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, - bool LongSegment = false, 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::type::value, @@ -96,19 +107,23 @@ public : typedef excess_sum state_type; typedef geometry::srs::sphere sphere_type; - inline spherical(sphere_type sphere = sphere_type()) - : m_sphere(sphere) + // For backward compatibility reasons the radius is set to 1 + inline spherical() + : m_sphere(1.0) {} - inline spherical(CT radius) //backward compatibility - : m_sphere() - { - geometry::set_radius<0>(m_sphere, radius); - } + template + explicit inline spherical(geometry::srs::sphere 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 + PointOfSegment const& p2, + excess_sum& state) const { if (! geometry::math::equals(get<0>(p1), get<0>(p2))) { @@ -164,4 +179,4 @@ struct default_strategy }} // namespace boost::geometry -#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_SPHERICAL_HPP +#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HPP diff --git a/include/boost/geometry/strategies/spherical/azimuth_spherical.hpp b/include/boost/geometry/strategies/spherical/azimuth.hpp similarity index 60% rename from include/boost/geometry/strategies/spherical/azimuth_spherical.hpp rename to include/boost/geometry/strategies/spherical/azimuth.hpp index 52c0a4ef29..3c208fe2e2 100644 --- a/include/boost/geometry/strategies/spherical/azimuth_spherical.hpp +++ b/include/boost/geometry/strategies/spherical/azimuth.hpp @@ -8,12 +8,17 @@ // 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_SPHERICAL_HPP -#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AZIMUTH_SPHERICAL_HPP +#ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AZIMUTH_HPP +#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AZIMUTH_HPP + #include #include +#include +#include + + namespace boost { namespace geometry { @@ -22,7 +27,7 @@ namespace strategy { namespace azimuth template < - typename CalculationType + typename CalculationType = void > class spherical { @@ -31,15 +36,20 @@ public : inline spherical() {} - inline void apply(CalculationType const& lon1_rad, - CalculationType const& lat1_rad, - CalculationType const& lon2_rad, - CalculationType const& lat2_rad, - CalculationType& a1, - CalculationType& a2) const + template + static inline void apply(T const& lon1_rad, T const& lat1_rad, + T const& lon2_rad, T const& lat2_rad, + T& a1, T& a2) { - geometry::formula::result_spherical result = geometry::formula:: - spherical_azimuth(lon1_rad, lat1_rad, lon2_rad, lat2_rad); + typedef typename boost::mpl::if_ + < + boost::is_void, T, CalculationType + >::type calc_t; + + geometry::formula::result_spherical + result = geometry::formula::spherical_azimuth( + calc_t(lon1_rad), calc_t(lat1_rad), + calc_t(lon2_rad), calc_t(lat2_rad)); a1 = result.azimuth; a2 = result.reverse_azimuth; @@ -74,4 +84,4 @@ struct default_strategy }} // namespace boost::geometry -#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AZIMUTH_SPHERICAL_HPP +#endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AZIMUTH_HPP diff --git a/include/boost/geometry/strategies/spherical/envelope_segment.hpp b/include/boost/geometry/strategies/spherical/envelope_segment.hpp index 1003c18c2a..98f085fe73 100644 --- a/include/boost/geometry/strategies/spherical/envelope_segment.hpp +++ b/include/boost/geometry/strategies/spherical/envelope_segment.hpp @@ -2,6 +2,7 @@ // 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 @@ -13,7 +14,7 @@ #include #include #include -#include +#include namespace boost { namespace geometry { @@ -21,7 +22,10 @@ namespace boost { namespace geometry namespace strategy { namespace envelope { -template +template +< + typename CalculationType = void +> class spherical_segment { public : diff --git a/include/boost/geometry/strategies/spherical/intersection.hpp b/include/boost/geometry/strategies/spherical/intersection.hpp index b8b8bf93e9..5d37583333 100644 --- a/include/boost/geometry/strategies/spherical/intersection.hpp +++ b/include/boost/geometry/strategies/spherical/intersection.hpp @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -38,6 +39,8 @@ #include #include #include +#include +#include #include #include @@ -73,8 +76,12 @@ 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 -struct relate_spherical_segments +template +< + typename CalcPolicy, + typename CalculationType = void +> +struct ecef_segments { typedef side::spherical_side_formula side_strategy_type; @@ -106,8 +113,43 @@ struct relate_spherical_segments return strategy_type(); } + template + struct area_strategy + { + typedef area::spherical + < + typename point_type::type, + CalculationType + > type; + }; + + template + static inline typename area_strategy::type get_area_strategy() + { + typedef typename area_strategy::type strategy_type; + return strategy_type(); + } + + template + struct distance_strategy + { + typedef distance::haversine + < + typename coordinate_type::type, + CalculationType + > type; + }; + + template + static inline typename distance_strategy::type get_distance_strategy() + { + typedef typename distance_strategy::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 struct segment_intersection_info { @@ -116,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(); @@ -143,7 +189,7 @@ struct relate_spherical_segments if (ip_flag == ipi_inters) { // TODO: assign the rest of coordinates - point = formula::cart3d_to_sph(intersection_point); + point = calc_policy.template from_cart3d(intersection_point); } else if (ip_flag == ipi_at_a1) { @@ -167,6 +213,8 @@ struct relate_spherical_segments SegmentRatio robust_ra; SegmentRatio robust_rb; intersection_point_flag ip_flag; + + CalcPolicy const& calc_policy; }; // Relate segments a and b @@ -210,6 +258,12 @@ struct relate_spherical_segments 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) ); BOOST_CONCEPT_ASSERT( (concepts::ConstSegment) ); @@ -234,25 +288,29 @@ struct relate_spherical_segments typedef model::point vec3d_t; - using namespace formula; - vec3d_t const a1v = sph_to_cart3d(a1); - vec3d_t const a2v = sph_to_cart3d(a2); - vec3d_t const b1v = sph_to_cart3d(b1); - vec3d_t const b2v = sph_to_cart3d(b2); + vec3d_t const a1v = calc_policy.template to_cart3d(a1); + vec3d_t const a2v = calc_policy.template to_cart3d(a2); + vec3d_t const b1v = calc_policy.template to_cart3d(b1); + vec3d_t const b2v = calc_policy.template to_cart3d(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 + 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 + 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 @@ -261,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) @@ -275,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) @@ -289,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) @@ -304,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); @@ -324,12 +395,12 @@ struct relate_spherical_segments { if (a_is_point) { - return collinear_one_degenerted(a, true, b1, b2, a1, a2, b1v, b2v, norm2, a1v); + return collinear_one_degenerated(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(b, false, a1, a2, b1, b2, a1v, a2v, norm1, b1v); + return collinear_one_degenerated(b, false, a1, a2, b1, b2, a1v, a2v, plane1, b1v); } else { @@ -338,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; @@ -408,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 @@ -417,7 +489,7 @@ struct relate_spherical_segments calc_t, segment_ratio, 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); @@ -434,31 +506,32 @@ struct relate_spherical_segments } private: - template + template static inline typename Policy::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) + 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(dist_1_o, dist_1_2), degenerated_a); } - template - 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 + 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)) @@ -477,54 +550,53 @@ struct relate_spherical_segments return segment_ratio(dist_a1_i1, dist_a1_a2).on_segment(); } - template + template 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; } @@ -535,8 +607,8 @@ struct relate_spherical_segments { 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; @@ -544,8 +616,8 @@ struct relate_spherical_segments 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; @@ -556,8 +628,8 @@ struct relate_spherical_segments { 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; @@ -565,8 +637,8 @@ struct relate_spherical_segments 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; @@ -580,7 +652,7 @@ struct relate_spherical_segments { 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; @@ -588,7 +660,7 @@ struct relate_spherical_segments 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; @@ -599,7 +671,7 @@ struct relate_spherical_segments { 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; @@ -607,7 +679,7 @@ struct relate_spherical_segments 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; @@ -619,24 +691,26 @@ struct relate_spherical_segments return is_on_a && is_on_b; } - template - 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 + 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] } @@ -716,6 +790,91 @@ struct relate_spherical_segments } }; +struct spherical_segments_calc_policy +{ + template + static Point from_cart3d(Point3d const& point_3d) + { + return formula::cart3d_to_sph(point_3d); + } + + template + static Point3d to_cart3d(Point const& point) + { + return formula::sph_to_cart3d(point); + } + + template + struct plane + { + typedef typename coordinate_type::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 + static plane get_plane(Point3d const& p1, Point3d const& p2) + { + return plane(p1, p2); + } + + template + static bool intersection_points(plane const& plane1, + plane const& plane2, + Point3d & ip1, Point3d & ip2) + { + typedef typename coordinate_type::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 @@ -724,19 +883,24 @@ namespace services /*template struct default_strategy { - typedef relate_spherical_segments type; + typedef spherical_segments type; };*/ template struct default_strategy { - typedef relate_spherical_segments type; + typedef spherical_segments type; }; template struct default_strategy { - typedef relate_spherical_segments 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 type; }; } // namespace services @@ -755,25 +919,25 @@ namespace within { namespace services template struct default_strategy { - typedef strategy::intersection::relate_spherical_segments<> type; + typedef strategy::intersection::spherical_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_spherical_segments<> type; + typedef strategy::intersection::spherical_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_spherical_segments<> type; + typedef strategy::intersection::spherical_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_spherical_segments<> type; + typedef strategy::intersection::spherical_segments<> type; }; }} // within::services @@ -784,25 +948,25 @@ namespace covered_by { namespace services template struct default_strategy { - typedef strategy::intersection::relate_spherical_segments<> type; + typedef strategy::intersection::spherical_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_spherical_segments<> type; + typedef strategy::intersection::spherical_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_spherical_segments<> type; + typedef strategy::intersection::spherical_segments<> type; }; template struct default_strategy { - typedef strategy::intersection::relate_spherical_segments<> type; + typedef strategy::intersection::spherical_segments<> type; }; }} // within::services diff --git a/include/boost/geometry/strategies/strategies.hpp b/include/boost/geometry/strategies/strategies.hpp index 35763675a9..f3c6787a4b 100644 --- a/include/boost/geometry/strategies/strategies.hpp +++ b/include/boost/geometry/strategies/strategies.hpp @@ -41,7 +41,7 @@ #include #include -#include +#include #include #include #include @@ -61,13 +61,14 @@ #include #include #include +#include #include #include #include #include -#include -#include +#include +#include #include #include #include @@ -76,15 +77,19 @@ #include #include -#include -#include +#include +#include +#include #include #include #include #include -//#include -//#include -//#include +#include +//#include +#include +#include +#include +#include #include #include diff --git a/include/boost/geometry/util/math.hpp b/include/boost/geometry/util/math.hpp index 1874fe76a2..b1c3648c98 100644 --- a/include/boost/geometry/util/math.hpp +++ b/include/boost/geometry/util/math.hpp @@ -71,6 +71,19 @@ inline T const& greatest(T const& v1, T const& v2, T const& v3, T const& v4, T c } +template +inline T bounded(T const& v, T const& lower, T const& upper) +{ + return (std::min)((std::max)(v, lower), upper); +} + +template +inline T bounded(T const& v, T const& lower) +{ + return (std::max)(v, lower); +} + + template ::value> struct abs diff --git a/include/boost/geometry/util/normalize_spheroidal_coordinates.hpp b/include/boost/geometry/util/normalize_spheroidal_coordinates.hpp index 377051eef1..19d4d33d28 100644 --- a/include/boost/geometry/util/normalize_spheroidal_coordinates.hpp +++ b/include/boost/geometry/util/normalize_spheroidal_coordinates.hpp @@ -282,6 +282,40 @@ inline CoordinateType longitude_distance_unsigned(CoordinateType const& longitud return diff; } +/*! +\brief The abs difference between longitudes in range [0, 180]. +\tparam Units The units of the coordindate system in the spheroid +\tparam CoordinateType The type of the coordinates +\param longitude1 Longitude 1 +\param longitude2 Longitude 2 +\ingroup utility +*/ +template +inline CoordinateType longitude_difference(CoordinateType const& longitude1, + CoordinateType const& longitude2) +{ + return math::abs(math::longitude_distance_signed(longitude1, longitude2)); +} + +template +inline CoordinateType longitude_interval_distance_signed(CoordinateType const& longitude_a1, + CoordinateType const& longitude_a2, + CoordinateType const& longitude_b) +{ + CoordinateType const c0 = 0; + CoordinateType dist_a12 = longitude_distance_signed(longitude_a1, longitude_a2); + CoordinateType dist_a1b = longitude_distance_signed(longitude_a1, longitude_b); + if (dist_a12 < c0) + { + dist_a12 = -dist_a12; + dist_a1b = -dist_a1b; + } + + return dist_a1b < c0 ? dist_a1b + : dist_a1b > dist_a12 ? dist_a1b - dist_a12 + : c0; +} + } // namespace math diff --git a/test/algorithms/area/area.cpp b/test/algorithms/area/area.cpp index e35e8f7095..ce41150496 100644 --- a/test/algorithms/area/area.cpp +++ b/test/algorithms/area/area.cpp @@ -5,8 +5,8 @@ // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. // Copyright (c) 2009-2012 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, 2016, 2017. +// Modifications copyright (c) 2015-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 @@ -208,17 +208,20 @@ int test_main(int, char* []) typedef bg::model::point > pt_sph; typedef bg::model::point > pt_geo; + // mean Earth's radius^2 + double r2 = bg::math::sqr(bg::get_radius<0>(bg::srs::sphere())); + test_ccw(); test_ccw(); test_ccw(); test_open(2.0); - test_open(24726179921.523518); - test_open(24615760871.487991); + test_open(24726179921.523518 / r2); + test_open(24615492936.977146); test_open_ccw(2.0); - test_open_ccw(24726179921.523518); - test_open_ccw(24615760871.487991); + test_open_ccw(24726179921.523518 / r2); + test_open_ccw(24615492936.977146); test_poles_ccw(); test_poles_ccw(); diff --git a/test/algorithms/area/area_geo.cpp b/test/algorithms/area/area_geo.cpp index 9a88eccf94..284d179cdc 100644 --- a/test/algorithms/area/area_geo.cpp +++ b/test/algorithms/area/area_geo.cpp @@ -1,10 +1,10 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) // Unit Test -// 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, 2016, 2017. +// Modifications copyright (c) 2015-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 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -24,134 +24,49 @@ void test_geo_strategies() { std::string poly = "POLYGON((52 0, 41 -74, -23 -43, -26 28, 52 0))"; - typedef bg::model::point > pt_geo; - - typedef typename bg::point_type::type pt_geo_type; - - bg::strategy::area::geographic - geographic_default; - - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::andoyer_inverse, - 1, - true - > geographic_andoyer1; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::andoyer_inverse, - 2, - true - > geographic_andoyer2; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::andoyer_inverse, - 3, - true - > geographic_andoyer3; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::andoyer_inverse, - 4, - true - > geographic_andoyer4; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::andoyer_inverse, - 5, - true - > geographic_andoyer5; - - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::thomas_inverse, - 1, - true - > geographic_thomas1; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::thomas_inverse, - 2, - true - > geographic_thomas2; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::thomas_inverse, - 3, - true - > geographic_thomas3; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::thomas_inverse, - 4, - true - > geographic_thomas4; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::thomas_inverse, - 5, - true - > geographic_thomas5; - - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::vincenty_inverse, - 1, - true - > geographic_vincenty1; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::vincenty_inverse, - 2, - true - > geographic_vincenty2; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::vincenty_inverse, - 3, - true - > geographic_vincenty3; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::vincenty_inverse, - 4, - true - > geographic_vincenty4; - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::vincenty_inverse, - 5, - true - > geographic_vincenty5; - - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::vincenty_inverse, - 5 - > geographic_vincenty5_default; - - bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::vincenty_inverse - > geographic_vincenty_default; + typedef bg::model::point > pt_geo; + + bg::strategy::area::geographic geographic_default; + + bg::strategy::area::geographic + geographic_andoyer1; + bg::strategy::area::geographic + geographic_andoyer2; + bg::strategy::area::geographic + geographic_andoyer3; + bg::strategy::area::geographic + geographic_andoyer4; + bg::strategy::area::geographic + geographic_andoyer5; + + bg::strategy::area::geographic + geographic_thomas1; + bg::strategy::area::geographic + geographic_thomas2; + bg::strategy::area::geographic + geographic_thomas3; + bg::strategy::area::geographic + geographic_thomas4; + bg::strategy::area::geographic + geographic_thomas5; + + bg::strategy::area::geographic + geographic_vincenty1; + bg::strategy::area::geographic + geographic_vincenty2; + bg::strategy::area::geographic + geographic_vincenty3; + bg::strategy::area::geographic + geographic_vincenty4; + bg::strategy::area::geographic + geographic_vincenty5; + + bg::strategy::area::geographic + geographic_andoyer_default; + bg::strategy::area::geographic + geographic_thomas_default; + bg::strategy::area::geographic + geographic_vincenty_default; bg::model::polygon geometry_geo; @@ -164,12 +79,12 @@ void test_geo_strategies() CT err = 0.0000001; CT area_default = bg::area(geometry_geo); - BOOST_CHECK_CLOSE(area_default, 63316536092341.266, err); + BOOST_CHECK_CLOSE(area_default, 63316309346280.18, err); area = bg::area(geometry_geo, geographic_default); - BOOST_CHECK_CLOSE(area, 63316536092341.266, err); + BOOST_CHECK_CLOSE(area, 63316309346280.18, err); CT area_less_accurate = bg::area(geometry_geo, geographic_andoyer1); - BOOST_CHECK_CLOSE(area_less_accurate, 63316309346280.18, err); + BOOST_CHECK_CLOSE(area, 63316309346280.18, err); area = bg::area(geometry_geo, geographic_andoyer2); BOOST_CHECK_CLOSE(area, 63316309224306.5, err); area = bg::area(geometry_geo, geographic_andoyer3); @@ -201,13 +116,15 @@ void test_geo_strategies() CT area_most_accurate = bg::area(geometry_geo, geographic_vincenty5); BOOST_CHECK_CLOSE(area, 63316536351929.523, err); - area = bg::area(geometry_geo, geographic_vincenty5_default); - BOOST_CHECK_CLOSE(area, 63316536351929.523, err); + area = bg::area(geometry_geo, geographic_andoyer_default); + BOOST_CHECK_CLOSE(area, 63316309346280.18, err); + area = bg::area(geometry_geo, geographic_thomas_default); + BOOST_CHECK_CLOSE(area, 63316536092341.266, err); area = bg::area(geometry_geo, geographic_vincenty_default); - BOOST_CHECK_CLOSE(area, 63316536351824.93, err); + BOOST_CHECK_CLOSE(area, 63316536351929.523, err); BOOST_CHECK_CLOSE(area_most_accurate, area_less_accurate, .001); - BOOST_CHECK_CLOSE(area_most_accurate, area_default, .000001); + BOOST_CHECK_CLOSE(area_most_accurate, area_default, .001); /* // timings and accuracy std::cout.precision(25); diff --git a/test/algorithms/area/area_multi.cpp b/test/algorithms/area/area_multi.cpp index b351a9c017..37ac0cecf0 100644 --- a/test/algorithms/area/area_multi.cpp +++ b/test/algorithms/area/area_multi.cpp @@ -34,10 +34,13 @@ void test_all() typedef bg::model::multi_polygon > mp_sph; typedef bg::model::multi_polygon > mp_geo; + // mean Earth's radius^2 + double r2 = bg::math::sqr(bg::get_radius<0>(bg::srs::sphere())); + std::string poly = "MULTIPOLYGON(((0 0,0 7,4 2,2 0,0 0)))"; test_geometry(poly, 16.0); - test_geometry(poly, 197897454752.69489); - test_geometry(poly, 197022175077.78613); + test_geometry(poly, 197897454752.69489 / r2); + test_geometry(poly, 197018888665.8331); } int test_main( int , char* [] ) diff --git a/test/algorithms/area/area_sph_geo.cpp b/test/algorithms/area/area_sph_geo.cpp index 41e3a86a47..70185c1d8d 100644 --- a/test/algorithms/area/area_sph_geo.cpp +++ b/test/algorithms/area/area_sph_geo.cpp @@ -5,8 +5,8 @@ // Copyright (c) 2008-2012 Bruno Lalande, Paris, France. // Copyright (c) 2009-2012 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, 2016, 2017. +// Modifications copyright (c) 2015-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 @@ -31,30 +31,28 @@ void test_spherical_geo() //Geographic - typedef typename bg::model::point< - ct, - 2, - bg::cs::geographic - > pt_geo; + typedef typename bg::model::point + < + ct, 2, bg::cs::geographic + > pt_geo; typedef typename bg::point_type::type pt_geo_type; bg::strategy::area::geographic - < - pt_geo_type, - bg::formula::vincenty_inverse, - 5 - > area_geographic; + < + pt_geo_type, + bg::strategy::vincenty, + 5 + > area_geographic; bg::model::polygon geometry_geo; //Spherical - typedef typename bg::model::point< - ct, - 2, - bg::cs::spherical_equatorial - > pt; + typedef typename bg::model::point + < + ct, 2, bg::cs::spherical_equatorial + > pt; bg::model::polygon geometry; // unit-sphere has area of 4-PI. Polygon covering 1/8 of it: @@ -202,6 +200,9 @@ void test_spherical_geo() area = bg::area(geometry_geo, area_geographic); BOOST_CHECK_CLOSE(area, 133353077343.10347, 0.0001); + // mean Earth's radius^2 + double r2 = bg::math::sqr(bg::get_radius<0>(bg::srs::sphere())); + // around 0 meridian { std::string poly1 = "POLYGON((-10 0,-10 10,0 10,0 0,-10 0))"; @@ -215,7 +216,7 @@ void test_spherical_geo() ct area3 = bg::area(geometry); BOOST_CHECK_CLOSE(area1, area2, 0.001); BOOST_CHECK_CLOSE(area2, area3, 0.001); - BOOST_CHECK_CLOSE(area1, 1233204227903.1848, 0.001); + BOOST_CHECK_CLOSE(area1 * r2, 1233204227903.1848, 0.001); //geographic bg::read_wkt(poly1, geometry_geo); area1 = bg::area(geometry_geo, area_geographic); @@ -239,7 +240,7 @@ void test_spherical_geo() ct area3 = bg::area(geometry); BOOST_CHECK_CLOSE(area1, area2, 0.001); BOOST_CHECK_CLOSE(area2, area3, 0.001); - BOOST_CHECK_CLOSE(area1, 1237986107636.0261, 0.001); + BOOST_CHECK_CLOSE(area1 * r2, 1237986107636.0261, 0.001); //geographic bg::read_wkt(poly1, geometry_geo); area1 = bg::area(geometry_geo, area_geographic); @@ -272,7 +273,7 @@ void test_spherical_geo() BOOST_CHECK_CLOSE(area2, area3, 0.001); BOOST_CHECK_CLOSE(area3, area4, 0.001); BOOST_CHECK_CLOSE(area4, area5, 0.001); - BOOST_CHECK_CLOSE(area1, 1233204227903.1833, 0.001); + BOOST_CHECK_CLOSE(area1 * r2, 1233204227903.1833, 0.001); //geographic bg::read_wkt(poly1, geometry_geo); area1 = bg::area(geometry_geo, area_geographic); @@ -310,7 +311,7 @@ void test_spherical_geo() BOOST_CHECK_CLOSE(area2, area3, 0.001); BOOST_CHECK_CLOSE(area3, area4, 0.001); BOOST_CHECK_CLOSE(area4, area5, 0.001); - BOOST_CHECK_CLOSE(area1, 1237986107636.0247, 0.001); + BOOST_CHECK_CLOSE(area1 * r2, 1237986107636.0247, 0.001); //geographic bg::read_wkt(poly1, geometry_geo); area1 = bg::area(geometry_geo, area_geographic); diff --git a/test/algorithms/envelope_expand/envelope_on_spheroid.cpp b/test/algorithms/envelope_expand/envelope_on_spheroid.cpp index a6336132a4..278a44dd97 100644 --- a/test/algorithms/envelope_expand/envelope_on_spheroid.cpp +++ b/test/algorithms/envelope_expand/envelope_on_spheroid.cpp @@ -1,7 +1,7 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) // Unit Test -// Copyright (c) 2015-2016, Oracle and/or its affiliates. +// Copyright (c) 2015-2017, Oracle and/or its affiliates. // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle @@ -10,9 +10,6 @@ // Licensed under the Boost Software License version 1.0. // http://www.boost.org/users/license.html -#include -#include -#include #ifndef BOOST_TEST_MODULE #define BOOST_TEST_MODULE test_envelope_on_sphere_or_spheroid @@ -50,10 +47,8 @@ #include "test_envelope_expand_on_spheroid.hpp" -template < - template class Inverse, - typename CS_Tag - > + +template struct test_envelope { template @@ -63,8 +58,8 @@ struct test_envelope } }; -template