diff options
Diffstat (limited to 'boost/geometry')
36 files changed, 2363 insertions, 635 deletions
diff --git a/boost/geometry/algorithms/detail/azimuth.hpp b/boost/geometry/algorithms/detail/azimuth.hpp index 8948bf6afd..7e0d1691ed 100644 --- a/boost/geometry/algorithms/detail/azimuth.hpp +++ b/boost/geometry/algorithms/detail/azimuth.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014. -// Modifications copyright (c) 2014, Oracle and/or its affiliates. +// This file was modified by Oracle on 2014, 2016. +// Modifications copyright (c) 2014-2016, Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -22,7 +22,8 @@ #include <boost/geometry/util/math.hpp> #include <boost/geometry/algorithms/not_implemented.hpp> -#include <boost/geometry/algorithms/detail/vincenty_inverse.hpp> + +#include <boost/geometry/formulas/vincenty_inverse.hpp> namespace boost { namespace geometry { @@ -49,7 +50,7 @@ struct azimuth<ReturnType, geographic_tag> template <typename P1, typename P2, typename Spheroid> static inline ReturnType apply(P1 const& p1, P2 const& p2, Spheroid const& spheroid) { - return geometry::detail::vincenty_inverse<ReturnType, false, true>().apply + return geometry::formula::vincenty_inverse<ReturnType, false, true>().apply ( get_as_radian<0>(p1), get_as_radian<1>(p1), get_as_radian<0>(p2), get_as_radian<1>(p2), spheroid ).azimuth; diff --git a/boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp b/boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp new file mode 100644 index 0000000000..df62a1f2f6 --- /dev/null +++ b/boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp @@ -0,0 +1,100 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2016 Barend Gehrels, Amsterdam, the Netherlands. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_AGGREGATE_OPERATIONS_HPP +#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_AGGREGATE_OPERATIONS_HPP + +#include <set> + +#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp> + +namespace boost { namespace geometry +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail { namespace overlay { namespace sort_by_side +{ + +struct ring_with_direction +{ + ring_identifier ring_id; + direction_type direction; + bool only_turn_on_ring; + + + inline bool operator<(ring_with_direction const& other) const + { + return this->ring_id != other.ring_id + ? this->ring_id < other.ring_id + : this->direction < other.direction; + } + + ring_with_direction() + : direction(dir_unknown) + , only_turn_on_ring(false) + {} +}; + +struct rank_with_rings +{ + std::size_t rank; + std::set<ring_with_direction> rings; + + rank_with_rings() + : rank(0) + { + } + + inline bool all_to() const + { + for (std::set<ring_with_direction>::const_iterator it = rings.begin(); + it != rings.end(); ++it) + { + if (it->direction == sort_by_side::dir_from) + { + return false; + } + } + return true; + } +}; + +template <typename Sbs> +inline void aggregate_operations(Sbs const& sbs, std::vector<rank_with_rings>& aggregation) +{ + aggregation.clear(); + for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++) + { + typename Sbs::rp const& ranked_point = sbs.m_ranked_points[i]; + + if (aggregation.empty() || aggregation.back().rank != ranked_point.rank) + { + rank_with_rings current; + current.rank = ranked_point.rank; + aggregation.push_back(current); + } + + ring_with_direction rwd; + segment_identifier const& sid = ranked_point.seg_id; + rwd.ring_id = ring_identifier(sid.source_index, sid.multi_index, sid.ring_index); + rwd.direction = ranked_point.direction; + rwd.only_turn_on_ring = ranked_point.only_turn_on_ring; + + + aggregation.back().rings.insert(rwd); + } +} + + +}}} // namespace detail::overlay::sort_by_side +#endif //DOXYGEN_NO_DETAIL + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_AGGREGATE_OPERATIONS_HPP diff --git a/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp b/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp index cc55414870..2643415343 100644 --- a/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp +++ b/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp @@ -38,6 +38,7 @@ struct enrichment_info , count_left(0) , count_right(0) , zone(-1) + , only_turn_on_ring(false) {} // vertex to which is free travel after this IP, @@ -57,6 +58,7 @@ struct enrichment_info std::size_t count_left; std::size_t count_right; signed_size_type zone; // open zone, in cluster + bool only_turn_on_ring; // True if it is the only turn on a ring (for clusters) }; diff --git a/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp b/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp index 110e1be2f1..400ed3b881 100644 --- a/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp +++ b/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp @@ -14,6 +14,7 @@ #include <map> #include <vector> +#include <boost/core/ignore_unused.hpp> #include <boost/range.hpp> #include <boost/geometry/core/point_order.hpp> #include <boost/geometry/algorithms/detail/overlay/cluster_info.hpp> @@ -198,7 +199,9 @@ inline signed_size_type add_turn_to_cluster(Turn const& turn, // Both operations.seg_id/fraction were already part of any cluster, and // these clusters are not the same. Merge of two clusters is necessary +#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_COLOCATIONS) std::cout << " TODO: merge " << cid0 << " and " << cid1 << std::endl; +#endif return cid0; } @@ -314,11 +317,13 @@ inline void assign_cluster_to_turns(Turns& turns, typename ClusterPerSegment::const_iterator it = cluster_per_segment.find(seg_frac); if (it != cluster_per_segment.end()) { +#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_COLOCATIONS) if (turn.cluster_id != -1 && turn.cluster_id != it->second) { std::cout << " CONFLICT " << std::endl; } +#endif turn.cluster_id = it->second; clusters[turn.cluster_id].turn_indices.insert(turn_index); } @@ -392,6 +397,8 @@ inline bool is_ie_turn(segment_identifier const& ext_seg_0, bool const same_multi1 = ! Reverse1 && ext_seg_1.multi_index == other_seg_1.multi_index; + boost::ignore_unused(same_multi1); + return same_multi0 && same_multi1 && ! is_interior<Reverse0>(ext_seg_0) @@ -680,6 +687,7 @@ inline void gather_cluster_properties(Clusters& clusters, Turns& turns, sbs.apply(turn_point); sbs.find_open(); + sbs.assign_zones(for_operation); // Unset the startable flag for all 'closed' zones for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++) @@ -706,7 +714,7 @@ inline void gather_cluster_properties(Clusters& clusters, Turns& turns, } } - cinfo.open_count = sbs.open_count(turns); + cinfo.open_count = sbs.open_count(for_operation); } } diff --git a/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp b/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp index 91b0f8ae24..bbba623eee 100644 --- a/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp +++ b/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp @@ -37,10 +37,12 @@ struct ranked_point , count_left(0) , count_right(0) , operation(operation_none) + , only_turn_on_ring(false) {} - ranked_point(const Point& p, signed_size_type ti, signed_size_type oi, - direction_type d, operation_type op, segment_identifier sid) + template <typename Op> + ranked_point(const Point& p, signed_size_type ti, int oi, + direction_type d, Op op) : point(p) , rank(0) , zone(-1) @@ -49,20 +51,22 @@ struct ranked_point , direction(d) , count_left(0) , count_right(0) - , operation(op) - , seg_id(sid) + , operation(op.operation) + , seg_id(op.seg_id) + , only_turn_on_ring(op.enriched.only_turn_on_ring) {} Point point; std::size_t rank; signed_size_type zone; // index of closed zone, in uu turn there would be 2 zones signed_size_type turn_index; - signed_size_type operation_index; + int operation_index; // 0,1 direction_type direction; std::size_t count_left; std::size_t count_right; operation_type operation; segment_identifier seg_id; + bool only_turn_on_ring; }; struct less_by_turn_index @@ -181,11 +185,36 @@ private : Point m_p1, m_p2; }; +// Sorts vectors in counter clockwise order (by default) template <bool Reverse1, bool Reverse2, typename Point, typename Compare> struct side_sorter { typedef ranked_point<Point> rp; +private : + struct include_union + { + inline bool operator()(rp const& ranked_point) const + { + // New candidate if there are no polygons on left side, + // but there are on right side + return ranked_point.count_left == 0 + && ranked_point.count_right > 0; + } + }; + + struct include_intersection + { + inline bool operator()(rp const& ranked_point) const + { + // New candidate if there are two polygons on right side, + // and less on the left side + return ranked_point.count_left < 2 + && ranked_point.count_right >= 2; + } + }; + +public : inline void set_origin(Point const& origin) { m_origin = origin; @@ -202,8 +231,8 @@ struct side_sorter op.seg_id, point1, point2, point3); Point const& point_to = op.fraction.is_one() ? point3 : point2; - m_ranked_points.push_back(rp(point1, turn_index, op_index, dir_from, op.operation, op.seg_id)); - m_ranked_points.push_back(rp(point_to, turn_index, op_index, dir_to, op.operation, op.seg_id)); + m_ranked_points.push_back(rp(point1, turn_index, op_index, dir_from, op)); + m_ranked_points.push_back(rp(point_to, turn_index, op_index, dir_to, op)); if (is_origin) { @@ -292,8 +321,6 @@ struct side_sorter &segment_identifier::source_index >(handled); } - - assign_zones(); } void reverse() @@ -303,7 +330,7 @@ struct side_sorter return; } - int const last = 1 + m_ranked_points.back().rank; + std::size_t const last = 1 + m_ranked_points.back().rank; // Move iterator after rank==0 bool has_first = false; @@ -334,13 +361,18 @@ struct side_sorter } } +//private : + + typedef std::vector<rp> container_type; + container_type m_ranked_points; + Point m_origin; + +private : + //! Check how many open spaces there are - template <typename Turns> - std::size_t open_count(Turns const& turns) const + template <typename Include> + inline std::size_t open_count(Include const& include_functor) const { - typedef typename boost::range_value<Turns>::type turn_type; - typedef typename turn_type::turn_operation_type turn_operation_type; - std::size_t result = 0; std::size_t last_rank = 0; for (std::size_t i = 0; i < m_ranked_points.size(); i++) @@ -348,31 +380,16 @@ struct side_sorter rp const& ranked_point = m_ranked_points[i]; if (ranked_point.rank > last_rank - && ranked_point.direction == sort_by_side::dir_to) + && ranked_point.direction == sort_by_side::dir_to + && include_functor(ranked_point)) { - // TODO: take count-left / count_right from rank itself - turn_type const& ranked_turn = turns[ranked_point.turn_index]; - turn_operation_type const& ranked_op = ranked_turn.operations[ranked_point.operation_index]; - if (ranked_op.enriched.count_left == 0 - && ranked_op.enriched.count_right > 0) - { - result++; - last_rank = ranked_point.rank; - } + result++; + last_rank = ranked_point.rank; } } return result; } -//protected : - - typedef std::vector<rp> container_type; - container_type m_ranked_points; - Point m_origin; - -private : - - std::size_t move(std::size_t index) const { std::size_t const result = index + 1; @@ -458,7 +475,8 @@ private : } //! Find closed zones and assign it - void assign_zones() + template <typename Include> + std::size_t assign_zones(Include const& include_functor) { // Find a starting point (the first rank after an outgoing rank // with no polygons on the left side) @@ -473,8 +491,7 @@ private : max_rank = ranked_point.rank; } if (ranked_point.direction == sort_by_side::dir_to - && ranked_point.count_left == 0 - && ranked_point.count_right > 0) + && include_functor(ranked_point)) { start_rank = ranked_point.rank + 1; } @@ -510,8 +527,7 @@ private : } if (ranked_point.direction == sort_by_side::dir_to - && ranked_point.count_left == 0 - && ranked_point.count_right > 0) + && include_functor(ranked_point)) { rank_at_next_zone = ranked_point.rank + 1; if (rank_at_next_zone > max_rank) @@ -525,8 +541,25 @@ private : ranked_point.zone = zone_id; } + return zone_id; } +public : + inline std::size_t open_count(operation_type for_operation) const + { + return for_operation == operation_union + ? open_count(include_union()) + : open_count(include_intersection()) + ; + } + + inline std::size_t assign_zones(operation_type for_operation) + { + return for_operation == operation_union + ? assign_zones(include_union()) + : assign_zones(include_intersection()) + ; + } }; diff --git a/boost/geometry/algorithms/detail/overlay/traversal.hpp b/boost/geometry/algorithms/detail/overlay/traversal.hpp index 1156644bc3..5adc0fcf69 100644 --- a/boost/geometry/algorithms/detail/overlay/traversal.hpp +++ b/boost/geometry/algorithms/detail/overlay/traversal.hpp @@ -13,6 +13,7 @@ #include <boost/range.hpp> +#include <boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp> #include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp> #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp> #include <boost/geometry/core/access.hpp> @@ -152,8 +153,8 @@ struct traversal } } - inline bool is_visited(turn_type const& turn, turn_operation_type const& op, - signed_size_type turn_index, int op_index) const + inline bool is_visited(turn_type const& , turn_operation_type const& op, + signed_size_type , int) const { return op.visited.visited(); } @@ -162,39 +163,30 @@ struct traversal segment_identifier const& seg_id1, segment_identifier const& seg_id2) const { - if (target_operation == operation_intersection) + // For uu/ii, only switch sources if indicated + turn_type const& turn = m_turns[turn_index]; + + if (OverlayType == overlay_buffer) { - // For intersections always switch sources - return seg_id1.source_index != seg_id2.source_index; + // Buffer does not use source_index (always 0) + return turn.switch_source + ? seg_id1.multi_index != seg_id2.multi_index + : seg_id1.multi_index == seg_id2.multi_index; } - else if (target_operation == operation_union) - { - // For uu, only switch sources if indicated - turn_type const& turn = m_turns[turn_index]; - - if (OverlayType == overlay_buffer) - { - // Buffer does not use source_index (always 0) - return turn.switch_source - ? seg_id1.multi_index != seg_id2.multi_index - : seg_id1.multi_index == seg_id2.multi_index; - } #if defined(BOOST_GEOMETRY_DEBUG_TRAVERSAL_SWITCH_DETECTOR) - if (turn.switch_source == 1) - { - std::cout << "Switch source at " << turn_index << std::endl; - } - else - { - std::cout << "DON'T SWITCH SOURCES at " << turn_index << std::endl; - } -#endif - return turn.switch_source - ? seg_id1.source_index != seg_id2.source_index - : seg_id1.source_index == seg_id2.source_index; + if (turn.switch_source) + { + std::cout << "Switch source at " << turn_index << std::endl; } - return false; + else + { + std::cout << "DON'T SWITCH SOURCES at " << turn_index << std::endl; + } +#endif + return turn.switch_source + ? seg_id1.source_index != seg_id2.source_index + : seg_id1.source_index == seg_id2.source_index; } inline @@ -273,10 +265,6 @@ struct traversal segment_identifier const& seg_id, int& selected_op_index) const { - // For "ii", take the other one (alternate) - // UNLESS the other one is already visited - // For "uu", take the same one (see above); - bool result = false; for (int i = 0; i < 2; i++) @@ -347,13 +335,10 @@ struct traversal return true; } - inline bool select_from_cluster(signed_size_type& turn_index, + inline bool select_from_cluster_union(signed_size_type& turn_index, int& op_index, signed_size_type start_turn_index, sbs_type const& sbs, bool is_touching) const { - bool const is_union = target_operation == operation_union; - bool const is_intersection = target_operation == operation_intersection; - std::size_t selected_rank = 0; std::size_t min_rank = 0; bool result = false; @@ -386,11 +371,8 @@ struct traversal && (ranked_point.rank > min_rank || ranked_turn.both(operation_continue))) { - if ((is_union - && ranked_op.enriched.count_left == 0 + if (ranked_op.enriched.count_left == 0 && ranked_op.enriched.count_right > 0) - || (is_intersection - && ranked_op.enriched.count_right == 2)) { if (result && ranked_point.turn_index != start_turn_index) { @@ -401,16 +383,6 @@ struct traversal turn_index = ranked_point.turn_index; op_index = ranked_point.operation_index; - if (is_intersection - && ranked_turn.both(operation_intersection) - && ranked_op.visited.finalized()) - { - // Override: - // For a ii turn, even though one operation might be selected, - // it should take the other one if the first one is used in a completed ring - op_index = 1 - ranked_point.operation_index; - } - result = true; selected_rank = ranked_point.rank; } @@ -423,10 +395,94 @@ struct traversal return result; } + inline bool analyze_cluster_intersection(signed_size_type& turn_index, + int& op_index, + sbs_type const& sbs) const + { + std::vector<sort_by_side::rank_with_rings> aggregation; + sort_by_side::aggregate_operations(sbs, aggregation); + + std::size_t selected_rank = 0; + + for (std::size_t i = 0; i < aggregation.size(); i++) + { + sort_by_side::rank_with_rings const& rwr = aggregation[i]; + + if (i > 1 + && i - 1 == selected_rank + && rwr.rings.size() == 1) + { + sort_by_side::ring_with_direction const& rwd = *rwr.rings.begin(); + + if (rwd.only_turn_on_ring) + { + // Find if this arriving ring was leaving previously + sort_by_side::ring_with_direction leaving = rwd; + leaving.direction = sort_by_side::dir_to; + + sort_by_side::rank_with_rings const& previous = aggregation[i - 1]; + + if (previous.rings.size() == 1 + && previous.rings.count(leaving) == 1) + { + // It arrives back - if this is one of the selected, unselect it + selected_rank = 0; + } + } + } + + if (rwr.all_to()) + { + if (selected_rank == 0) + { + // Take the first (= right) where segments leave, + // having the polygon on the right side + selected_rank = rwr.rank; + } + } + } + + if (selected_rank > 0) + { + std::size_t selected_index = sbs.m_ranked_points.size(); + for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++) + { + typename sbs_type::rp const& ranked_point = sbs.m_ranked_points[i]; + + if (ranked_point.rank == selected_rank) + { + turn_type const& ranked_turn = m_turns[ranked_point.turn_index]; + turn_operation_type const& ranked_op = ranked_turn.operations[ranked_point.operation_index]; + + if (ranked_op.visited.finalized()) + { + // This direction is already traveled before, the same + // cannot be traveled again + return false; + } + + // Take the last turn from this rank + selected_index = i; + } + } + + if (selected_index < sbs.m_ranked_points.size()) + { + typename sbs_type::rp const& ranked_point = sbs.m_ranked_points[selected_index]; + turn_index = ranked_point.turn_index; + op_index = ranked_point.operation_index; + return true; + } + } + + return false; + } + inline bool select_turn_from_cluster(signed_size_type& turn_index, int& op_index, bool& is_touching, signed_size_type start_turn_index, - segment_identifier const& previous_seg_id) const + segment_identifier const& previous_seg_id, + bool is_start) const { bool const is_union = target_operation == operation_union; @@ -488,30 +544,76 @@ struct traversal sbs.apply(turn.point); -#if defined(BOOST_GEOMETRY_DEBUG_TRAVERSAL_SWITCH_DETECTOR) - is_touching = is_union && cinfo.open_count > 1; - if (is_touching) + bool result = false; + + if (is_union) { - if (cinfo.switch_source) + #if defined(BOOST_GEOMETRY_DEBUG_TRAVERSAL_SWITCH_DETECTOR) + is_touching = cinfo.open_count > 1; + if (is_touching) { - is_touching = false; - std::cout << "CLUSTER: SWITCH SOURCES at " << turn_index << std::endl; + if (cinfo.switch_source) + { + is_touching = false; + std::cout << "CLUSTER: SWITCH SOURCES at " << turn_index << std::endl; + } + else + { + std::cout << "CLUSTER: CONTINUE at " << turn_index << std::endl; + } } - else + #else + is_touching = cinfo.open_count > 1 && ! cinfo.switch_source; + #endif + + if (is_touching) { - std::cout << "CLUSTER: CONTINUE at " << turn_index << std::endl; + sbs.reverse(); } + + result = select_from_cluster_union(turn_index, op_index, start_turn_index, sbs, + is_touching); } -#else - is_touching = is_union && cinfo.open_count > 1 && ! cinfo.switch_source; -#endif - if (is_touching) + else { - sbs.reverse(); + if (is_start + && turn.both(operation_intersection) + && turn.operations[op_index].enriched.only_turn_on_ring) + { + // For an ii (usually interior ring), only turn on ring, + // reverse to take first exit + sbs.reverse(); + } + + result = analyze_cluster_intersection(turn_index, op_index, sbs); } + return result; + } - return select_from_cluster(turn_index, op_index, start_turn_index, sbs, - is_touching); + inline bool analyze_ii_intersection(signed_size_type& turn_index, int& op_index, + turn_type const& current_turn, + segment_identifier const& previous_seg_id) + { + sbs_type sbs; + + // Add this turn to the sort-by-side sorter + bool has_origin = false; + for (int i = 0; i < 2; i++) + { + turn_operation_type const& op = current_turn.operations[i]; + bool const is_origin = op.seg_id.source_index + == previous_seg_id.source_index; + has_origin = has_origin || is_origin; + sbs.add(op, turn_index, i, m_geometry1, m_geometry2, is_origin); + } + + if (! has_origin) + { + return false; + } + + sbs.apply(current_turn.point); + return analyze_cluster_intersection(turn_index, op_index, sbs); } inline void change_index_for_self_turn(signed_size_type& to_vertex_index, @@ -607,7 +709,7 @@ struct traversal return true; } - bool select_turn(signed_size_type start_turn_index, + bool select_turn(signed_size_type start_turn_index, int start_op_index, signed_size_type& turn_index, int& op_index, bool& is_touching, @@ -616,10 +718,37 @@ struct traversal segment_identifier const& previous_seg_id, bool is_start) { - if (m_turns[turn_index].cluster_id >= 0) + turn_type const& current_turn = m_turns[turn_index]; + + if (target_operation == operation_intersection) + { + bool const back_at_start_cluster + = current_turn.cluster_id >= 0 + && m_turns[start_turn_index].cluster_id == current_turn.cluster_id; + + if (turn_index == start_turn_index || back_at_start_cluster) + { + // Intersection can always be finished if returning + turn_index = start_turn_index; + op_index = start_op_index; + return true; + } + + if (current_turn.cluster_id < 0 + && current_turn.both(operation_intersection)) + { + if (analyze_ii_intersection(turn_index, op_index, current_turn, + previous_seg_id)) + { + return true; + } + } + } + + if (current_turn.cluster_id >= 0) { if (! select_turn_from_cluster(turn_index, op_index, is_touching, - start_turn_index, previous_seg_id)) + start_turn_index, previous_seg_id, is_start)) { return false; } @@ -631,8 +760,6 @@ struct traversal } else { - turn_type const& current_turn = m_turns[turn_index]; - op_index = starting_operation_index(current_turn); if (op_index == -1) { diff --git a/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp b/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp index 104bd6b8e7..e0dfee19a8 100644 --- a/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp +++ b/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp @@ -64,9 +64,7 @@ struct traversal_ring_creator , m_clusters(clusters) , m_robust_policy(robust_policy) , m_visitor(visitor) - , m_has_uu(false) { - } template <typename Ring> @@ -123,7 +121,8 @@ struct traversal_ring_creator } bool is_touching = false; - if (! m_trav.select_turn(start_turn_index, turn_index, op_index, + if (! m_trav.select_turn(start_turn_index, start_op_index, + turn_index, op_index, is_touching, previous_op_index, previous_turn_index, previous_seg_id, is_start)) @@ -189,6 +188,18 @@ struct traversal_ring_creator return traverse_error_none; } + if (start_turn.cluster_id >= 0) + { + turn_type const& turn = m_turns[current_turn_index]; + if (turn.cluster_id == start_turn.cluster_id) + { + turn_operation_type& op = m_turns[start_turn_index].operations[current_op_index]; + op.visited.set_finished(); + m_visitor.visit_traverse(m_turns, m_turns[current_turn_index], start_op, "Early finish (cluster)"); + return traverse_error_none; + } + } + std::size_t const max_iterations = 2 + 2 * m_turns.size(); for (std::size_t i = 0; i <= max_iterations; i++) { @@ -276,49 +287,21 @@ struct traversal_ring_creator template <typename Rings> void iterate(Rings& rings, std::size_t& finalized_ring_size, - typename Backtrack::state_type& state, - int pass) + typename Backtrack::state_type& state) { - if (pass == 1) - { - if (target_operation == operation_intersection) - { - // Second pass currently only used for uu - return; - } - if (! m_has_uu) - { - // There is no uu found in first pass - return; - } - } - - // Iterate through all unvisited points for (std::size_t turn_index = 0; turn_index < m_turns.size(); ++turn_index) { - turn_type const& start_turn = m_turns[turn_index]; + turn_type const& turn = m_turns[turn_index]; - if (start_turn.discarded || start_turn.blocked()) + if (turn.discarded || turn.blocked()) { // Skip discarded and blocked turns continue; } - if (target_operation == operation_union) - { - if (start_turn.both(operation_union)) - { - // Start with a uu-turn only in the second pass - m_has_uu = true; - if (pass == 0) - { - continue; - } - } - } for (int op_index = 0; op_index < 2; op_index++) { - traverse_with_operation(start_turn, turn_index, op_index, + traverse_with_operation(turn, turn_index, op_index, rings, finalized_ring_size, state); } } @@ -333,10 +316,6 @@ private: Clusters const& m_clusters; RobustPolicy const& m_robust_policy; Visitor& m_visitor; - - // Next member is only used for operation union - bool m_has_uu; - }; }} // namespace detail::overlay diff --git a/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp b/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp index 9381c66e0e..183131c74b 100644 --- a/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp +++ b/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp @@ -71,8 +71,8 @@ struct traversal_switch_detector { if (turn.cluster_id == -1) { - // If it is a uu-turn (non clustered), it is never same zone - return ! turn.both(operation_union); + // If it is a uu/ii-turn (non clustered), it is never same zone + return ! (turn.both(operation_union) || turn.both(operation_intersection)); } // It is a cluster, check zones of both operations @@ -110,11 +110,11 @@ struct traversal_switch_detector for (set_iterator sit = ring_turn_indices.begin(); sit != ring_turn_indices.end(); ++sit) { - int const turn_index = *sit; + signed_size_type const turn_index = *sit; turn_type const& turn = m_turns[turn_index]; if (! connects_same_zone(turn)) { - // This is a non clustered uu-turn, or a cluster connecting different 'zones' + // This is a non clustered uu/ii-turn, or a cluster connecting different 'zones' continue; } @@ -132,6 +132,56 @@ struct traversal_switch_detector } } + void check_turns_per_ring(ring_identifier const& ring_id, + std::set<signed_size_type> const& ring_turn_indices) + { + bool only_turn_on_ring = true; + if (ring_turn_indices.size() > 1) + { + // More turns on this ring. Only leave only_turn_on_ring true + // if they are all of the same cluster + int cluster_id = -1; + for (set_iterator sit = ring_turn_indices.begin(); + sit != ring_turn_indices.end(); ++sit) + { + turn_type const& turn = m_turns[*sit]; + if (turn.cluster_id == -1) + { + // Unclustered turn - and there are 2 or more turns + // so the ring has different turns + only_turn_on_ring = false; + break; + } + + // Clustered turn, check if it is the first or same as previous + if (cluster_id == -1) + { + cluster_id = turn.cluster_id; + } + else if (turn.cluster_id != cluster_id) + { + only_turn_on_ring = false; + break; + } + } + } + + // Assign result to matching operation (a turn is always on two rings) + for (set_iterator sit = ring_turn_indices.begin(); + sit != ring_turn_indices.end(); ++sit) + { + turn_type& turn = m_turns[*sit]; + for (int i = 0; i < 2; i++) + { + turn_operation_type& op = turn.operations[i]; + if (ring_id_by_seg_id(op.seg_id) == ring_id) + { + op.enriched.only_turn_on_ring = only_turn_on_ring; + } + } + } + } + void propagate_region(ring_identifier const& ring_id, int region_id) { std::map<ring_identifier, std::set<signed_size_type> >::const_iterator it = m_turns_per_ring.find(ring_id); @@ -168,6 +218,7 @@ struct traversal_switch_detector = m_turns_per_ring.begin(); it != m_turns_per_ring.end(); ++it) { create_region(it->first, it->second); + check_turns_per_ring(it->first, it->second); } // Now that all regions are filled, assign switch_source property @@ -204,7 +255,7 @@ struct traversal_switch_detector cinfo.switch_source = regions.size() == 1; } - // Iterate through all uu turns (non-clustered) + // Iterate through all uu/ii turns (non-clustered) for (std::size_t turn_index = 0; turn_index < m_turns.size(); ++turn_index) { turn_type& turn = m_turns[turn_index]; @@ -212,9 +263,9 @@ struct traversal_switch_detector if (turn.discarded || turn.blocked() || turn.cluster_id >= 0 - || ! turn.both(operation_union)) + || ! (turn.both(operation_union) || turn.both(operation_intersection))) { - // Skip discarded, blocked, non-uu and clustered turns + // Skip discarded, blocked, non-uu/ii and clustered turns continue; } @@ -242,9 +293,10 @@ struct traversal_switch_detector { turn_type const& turn = m_turns[turn_index]; - if (turn.both(operation_union) && turn.cluster_id < 0) + if ((turn.both(operation_union) || turn.both(operation_intersection)) + && turn.cluster_id < 0) { - std::cout << "UU SWITCH RESULT " + std::cout << "UU/II SWITCH RESULT " << turn_index << " -> " << turn.switch_source << std::endl; } diff --git a/boost/geometry/algorithms/detail/overlay/traverse.hpp b/boost/geometry/algorithms/detail/overlay/traverse.hpp index 2d2933ebdd..f01e50eb03 100644 --- a/boost/geometry/algorithms/detail/overlay/traverse.hpp +++ b/boost/geometry/algorithms/detail/overlay/traverse.hpp @@ -97,10 +97,7 @@ public : typename Backtrack::state_type state; - for (int pass = 0; pass < 2; pass++) - { - trav.iterate(rings, finalized_ring_size, state, pass); - } + trav.iterate(rings, finalized_ring_size, state); } }; diff --git a/boost/geometry/algorithms/detail/overlay/turn_info.hpp b/boost/geometry/algorithms/detail/overlay/turn_info.hpp index 73f266a04a..e09af126c6 100644 --- a/boost/geometry/algorithms/detail/overlay/turn_info.hpp +++ b/boost/geometry/algorithms/detail/overlay/turn_info.hpp @@ -13,6 +13,7 @@ #include <boost/array.hpp> #include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/algorithms/detail/signed_size_type.hpp> #include <boost/geometry/algorithms/detail/overlay/segment_identifier.hpp> #include <boost/geometry/algorithms/detail/overlay/overlay_type.hpp> @@ -88,7 +89,7 @@ struct turn_info Point point; method_type method; - int cluster_id; + signed_size_type cluster_id; // For multiple turns on same location, >= 0. Else -1 bool discarded; bool colocated; bool switch_source; // For u/u turns which can either switch or not diff --git a/boost/geometry/algorithms/detail/result_inverse.hpp b/boost/geometry/algorithms/detail/result_inverse.hpp deleted file mode 100644 index 01a1997e49..0000000000 --- a/boost/geometry/algorithms/detail/result_inverse.hpp +++ /dev/null @@ -1,44 +0,0 @@ -// Boost.Geometry - -// Copyright (c) 2015 Oracle and/or its affiliates. - -// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle - -// Use, modification and distribution is subject to the Boost Software License, -// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) - -#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_RESULT_INVERSE_HPP -#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RESULT_INVERSE_HPP - - -#include <boost/math/constants/constants.hpp> - -#include <boost/geometry/core/radius.hpp> -#include <boost/geometry/core/srs.hpp> - -#include <boost/geometry/util/math.hpp> - -#include <boost/geometry/algorithms/detail/flattening.hpp> - - -namespace boost { namespace geometry { namespace detail -{ - -template <typename T> -struct result_inverse -{ - void set(T const& d, T const& a) - { - distance = d; - azimuth = a; - } - - T distance; - T azimuth; -}; - -}}} // namespace boost::geometry::detail - - -#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_RESULT_INVERSE_HPP diff --git a/boost/geometry/algorithms/detail/thomas_inverse.hpp b/boost/geometry/algorithms/detail/thomas_inverse.hpp deleted file mode 100644 index 96b237e054..0000000000 --- a/boost/geometry/algorithms/detail/thomas_inverse.hpp +++ /dev/null @@ -1,191 +0,0 @@ -// Boost.Geometry - -// Copyright (c) 2015 Oracle and/or its affiliates. - -// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle - -// Use, modification and distribution is subject to the Boost Software License, -// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) - -#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_THOMAS_INVERSE_HPP -#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_THOMAS_INVERSE_HPP - - -#include <boost/math/constants/constants.hpp> - -#include <boost/geometry/core/radius.hpp> -#include <boost/geometry/core/srs.hpp> - -#include <boost/geometry/util/condition.hpp> -#include <boost/geometry/util/math.hpp> - -#include <boost/geometry/algorithms/detail/flattening.hpp> -#include <boost/geometry/algorithms/detail/result_inverse.hpp> - -namespace boost { namespace geometry { namespace detail -{ - -/*! -\brief The solution of the inverse problem of geodesics on latlong coordinates, - Forsyth-Andoyer-Lambert type approximation with second order terms. -\author See - - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965 - http://www.dtic.mil/docs/citations/AD0627893 - - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970 - http://www.dtic.mil/docs/citations/AD703541 -*/ -template <typename CT, bool EnableDistance, bool EnableAzimuth> -struct thomas_inverse -{ - typedef result_inverse<CT> result_type; - - template <typename T1, typename T2, typename Spheroid> - static inline result_type apply(T1 const& lon1, - T1 const& lat1, - T2 const& lon2, - T2 const& lat2, - Spheroid const& spheroid) - { - result_type result; - - // coordinates in radians - - if ( math::equals(lon1, lon2) - && math::equals(lat1, lat2) ) - { - result.set(CT(0), CT(0)); - return result; - } - - CT const f = detail::flattening<CT>(spheroid); - CT const one_minus_f = CT(1) - f; - -// CT const tan_theta1 = one_minus_f * tan(lat1); -// CT const tan_theta2 = one_minus_f * tan(lat2); -// CT const theta1 = atan(tan_theta1); -// CT const theta2 = atan(tan_theta2); - - CT const pi_half = math::pi<CT>() / CT(2); - CT const theta1 = math::equals(lat1, pi_half) ? lat1 : - math::equals(lat1, -pi_half) ? lat1 : - atan(one_minus_f * tan(lat1)); - CT const theta2 = math::equals(lat2, pi_half) ? lat2 : - math::equals(lat2, -pi_half) ? lat2 : - atan(one_minus_f * tan(lat2)); - - CT const theta_m = (theta1 + theta2) / CT(2); - CT const d_theta_m = (theta2 - theta1) / CT(2); - CT const d_lambda = lon2 - lon1; - CT const d_lambda_m = d_lambda / CT(2); - - CT const sin_theta_m = sin(theta_m); - CT const cos_theta_m = cos(theta_m); - CT const sin_d_theta_m = sin(d_theta_m); - CT const cos_d_theta_m = cos(d_theta_m); - CT const sin2_theta_m = math::sqr(sin_theta_m); - CT const cos2_theta_m = math::sqr(cos_theta_m); - CT const sin2_d_theta_m = math::sqr(sin_d_theta_m); - CT const cos2_d_theta_m = math::sqr(cos_d_theta_m); - CT const sin_d_lambda_m = sin(d_lambda_m); - CT const sin2_d_lambda_m = math::sqr(sin_d_lambda_m); - - CT const H = cos2_theta_m - sin2_d_theta_m; - CT const L = sin2_d_theta_m + H * sin2_d_lambda_m; - CT const cos_d = CT(1) - CT(2) * L; - CT const d = acos(cos_d); - CT const sin_d = sin(d); - - CT const one_minus_L = CT(1) - L; - - if ( math::equals(sin_d, CT(0)) - || math::equals(L, CT(0)) - || math::equals(one_minus_L, CT(0)) ) - { - result.set(CT(0), CT(0)); - return result; - } - - CT const U = CT(2) * sin2_theta_m * cos2_d_theta_m / one_minus_L; - CT const V = CT(2) * sin2_d_theta_m * cos2_theta_m / L; - CT const X = U + V; - CT const Y = U - V; - CT const T = d / sin_d; - //CT const D = CT(4) * math::sqr(T); - //CT const E = CT(2) * cos_d; - //CT const A = D * E; - //CT const B = CT(2) * D; - //CT const C = T - (A - E) / CT(2); - - if ( BOOST_GEOMETRY_CONDITION(EnableDistance) ) - { - //CT const n1 = X * (A + C*X); - //CT const n2 = Y * (B + E*Y); - //CT const n3 = D*X*Y; - - //CT const f_sqr = math::sqr(f); - //CT const f_sqr_per_64 = f_sqr / CT(64); - - CT const delta1d = f * (T*X-Y) / CT(4); - //CT const delta2d = f_sqr_per_64 * (n1 - n2 + n3); - - CT const a = get_radius<0>(spheroid); - - result.distance = a * sin_d * (T - delta1d); - //double S2 = a * sin_d * (T - delta1d + delta2d); - } - else - { - result.distance = CT(0); - } - - if ( BOOST_GEOMETRY_CONDITION(EnableAzimuth) ) - { - // NOTE: if both cos_latX == 0 then below we'd have 0 * INF - // it's a situation when the endpoints are on the poles +-90 deg - // in this case the azimuth could either be 0 or +-pi - // but above always 0 is returned - - // may also be used to calculate distance21 - //CT const D = CT(4) * math::sqr(T); - CT const E = CT(2) * cos_d; - //CT const A = D * E; - //CT const B = CT(2) * D; - // may also be used to calculate distance21 - CT const f_sqr = math::sqr(f); - CT const f_sqr_per_64 = f_sqr / CT(64); - - CT const F = CT(2)*Y-E*(CT(4)-X); - //CT const M = CT(32)*T-(CT(20)*T-A)*X-(B+CT(4))*Y; - CT const G = f*T/CT(2) + f_sqr_per_64; - CT const tan_d_lambda = tan(d_lambda); - CT const Q = -(F*G*tan_d_lambda) / CT(4); - - CT const d_lambda_p = (d_lambda + Q) / CT(2); - CT const tan_d_lambda_p = tan(d_lambda_p); - - CT const v = atan2(cos_d_theta_m, sin_theta_m * tan_d_lambda_p); - CT const u = atan2(-sin_d_theta_m, cos_theta_m * tan_d_lambda_p); - - CT const pi = math::pi<CT>(); - CT alpha1 = v + u; - if ( alpha1 > pi ) - { - alpha1 -= CT(2) * pi; - } - - result.azimuth = alpha1; - } - else - { - result.azimuth = CT(0); - } - - return result; - } -}; - -}}} // namespace boost::geometry::detail - - -#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_THOMAS_INVERSE_HPP diff --git a/boost/geometry/extensions/algorithms/inverse.hpp b/boost/geometry/extensions/algorithms/inverse.hpp deleted file mode 100644 index 44f415be42..0000000000 --- a/boost/geometry/extensions/algorithms/inverse.hpp +++ /dev/null @@ -1,59 +0,0 @@ -// Boost.Geometry (aka GGL, Generic Geometry Library) -// -// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands. - -// Use, modification and distribution is subject to the Boost Software License, -// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at -// http://www.boost.org/LICENSE_1_0.txt) - - -#ifndef BOOST_GEOMETRY_ALGORITHMS_INVERSE_HPP -#define BOOST_GEOMETRY_ALGORITHMS_INVERSE_HPP - -#include <boost/geometry.hpp> - -namespace boost { namespace geometry -{ - -// TODO: this is shared with sectionalize, move to somewhere else (assign?) -template <typename Box, typename Value> -inline void enlarge_box(Box& box, Value value) -{ - geometry::set<0, 0>(box, geometry::get<0, 0>(box) - value); - geometry::set<0, 1>(box, geometry::get<0, 1>(box) - value); - geometry::set<1, 0>(box, geometry::get<1, 0>(box) + value); - geometry::set<1, 1>(box, geometry::get<1, 1>(box) + value); -} - -// TODO: when this might be moved outside extensions it should of course -// input/output a Geometry, instead of a WKT -template <typename Geometry, typename Value> -inline std::string inverse(std::string const& wkt, Value margin) -{ - Geometry geometry; - read_wkt(wkt, geometry); - - geometry::correct(geometry); - - geometry::model::box<typename point_type<Geometry>::type> env; - geometry::envelope(geometry, env); - - // Make its envelope a bit larger - enlarge_box(env, margin); - - Geometry env_as_polygon; - geometry::convert(env, env_as_polygon); - - Geometry inversed_result; - geometry::difference(env_as_polygon, geometry, inversed_result); - - std::ostringstream out; - - out << geometry::wkt(inversed_result); - return out.str(); -} - -}} // namespace boost::geometry - - -#endif // BOOST_GEOMETRY_ALGORITHMS_INVERSE_HPP diff --git a/boost/geometry/algorithms/detail/andoyer_inverse.hpp b/boost/geometry/formulas/andoyer_inverse.hpp index 66ad6446e9..57b5ab5384 100644 --- a/boost/geometry/algorithms/detail/andoyer_inverse.hpp +++ b/boost/geometry/formulas/andoyer_inverse.hpp @@ -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_ANDOYER_INVERSE_HPP -#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_ANDOYER_INVERSE_HPP +#ifndef BOOST_GEOMETRY_FORMULAS_ANDOYER_INVERSE_HPP +#define BOOST_GEOMETRY_FORMULAS_ANDOYER_INVERSE_HPP #include <boost/math/constants/constants.hpp> @@ -21,10 +21,12 @@ #include <boost/geometry/util/math.hpp> #include <boost/geometry/algorithms/detail/flattening.hpp> -#include <boost/geometry/algorithms/detail/result_inverse.hpp> +#include <boost/geometry/formulas/differential_quantities.hpp> +#include <boost/geometry/formulas/result_inverse.hpp> -namespace boost { namespace geometry { namespace detail + +namespace boost { namespace geometry { namespace formula { /*! @@ -36,10 +38,22 @@ namespace boost { namespace geometry { namespace detail - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970 http://www.dtic.mil/docs/citations/AD703541 */ - -template <typename CT, bool EnableDistance, bool EnableAzimuth> -struct andoyer_inverse +template < + typename CT, + bool EnableDistance, + bool EnableAzimuth, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false +> +class andoyer_inverse { + static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale; + static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities; + static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities; + static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities; + +public: typedef result_inverse<CT> result_type; template <typename T1, typename T2, typename Spheroid> @@ -49,20 +63,20 @@ struct andoyer_inverse T2 const& lat2, Spheroid const& spheroid) { - CT const c0 = CT(0); - CT const c1 = CT(1); - CT const pi = math::pi<CT>(); - result_type result; // coordinates in radians if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) ) { - result.set(c0, c0); return result; } + CT const c0 = CT(0); + CT const c1 = CT(1); + CT const pi = math::pi<CT>(); + CT const f = detail::flattening<CT>(spheroid); + CT const dlon = lon2 - lon1; CT const sin_dlon = sin(dlon); CT const cos_dlon = cos(dlon); @@ -84,8 +98,6 @@ struct andoyer_inverse CT const d = acos(cos_d); // [0, pi] CT const sin_d = sin(d); // [-1, 1] - CT const f = detail::flattening<CT>(spheroid); - if ( BOOST_GEOMETRY_CONDITION(EnableDistance) ) { CT const K = math::sqr(sin_lat1-sin_lat2); @@ -109,12 +121,8 @@ struct andoyer_inverse result.distance = a * (d + dd); } - else - { - result.distance = c0; - } - if ( BOOST_GEOMETRY_CONDITION(EnableAzimuth) ) + if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) ) { // sin_d = 0 <=> antipodal points if (math::equals(sin_d, c0)) @@ -122,10 +130,7 @@ struct andoyer_inverse // T = inf // dA = inf // azimuth = -inf - if (lat1 <= lat2) - result.azimuth = c0; - else - result.azimuth = pi; + result.azimuth = lat1 <= lat2 ? c0 : pi; } else { @@ -142,64 +147,101 @@ struct andoyer_inverse U = (f/ c2)*math::sqr(cos_lat1)*sin_2A; } + CT B = c0; CT V = c0; if ( ! math::equals(cos_lat1, c0) ) { CT const tan_lat1 = sin_lat1/cos_lat1; CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon; - CT const B = atan2(sin_dlon, N); + B = atan2(sin_dlon, N); CT const sin_2B = sin(c2*B); V = (f/ c2)*math::sqr(cos_lat2)*sin_2B; } CT const T = d / sin_d; - CT const dA = V*T-U; - - result.azimuth = A - dA; // even with sin_d == 0 checked above if the second point // is somewhere in the antipodal area T may still be great - // therefore dA may be great and the resulting azimuth - // may be some more or less arbitrary angle - if (A >= c0) // A indicates Eastern hemisphere + // therefore dA and dB may be great and the resulting azimuths + // may be some more or less arbitrary angles + + if (BOOST_GEOMETRY_CONDITION(CalcFwdAzimuth)) { - if (dA >= c0) // A altered towards 0 - { - if ((result.azimuth) < c0) - result.azimuth = c0; - } - else // dA < 0, A altered towards pi - { - if (result.azimuth > pi) - result.azimuth = pi; - } + CT const dA = V*T - U; + result.azimuth = A - dA; + normalize_azimuth(result.azimuth, A, dA); } - else // A indicates Western hemisphere + + if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth)) { - if (dA <= c0) // A altered towards 0 - { - if (result.azimuth > c0) - result.azimuth = c0; - } - else // dA > 0, A altered towards -pi + CT const dB = -U*T + V; + result.reverse_azimuth = pi - B - dB; + if (result.reverse_azimuth > pi) { - CT const minus_pi = -pi; - if ((result.azimuth) < minus_pi) - result.azimuth = minus_pi; + result.reverse_azimuth -= 2 * pi; } + normalize_azimuth(result.reverse_azimuth, B, dB); } } } - else + + if (BOOST_GEOMETRY_CONDITION(CalcQuantities)) { - result.azimuth = c0; + typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 1> quantities; + quantities::apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2, + result.azimuth, result.reverse_azimuth, + get_radius<2>(spheroid), f, + result.reduced_length, result.geodesic_scale); } return result; } + +private: + static inline void normalize_azimuth(CT & azimuth, CT const& A, CT const& dA) + { + CT const c0 = 0; + + if (A >= c0) // A indicates Eastern hemisphere + { + if (dA >= c0) // A altered towards 0 + { + if (azimuth < c0) + { + azimuth = c0; + } + } + else // dA < 0, A altered towards pi + { + CT const pi = math::pi<CT>(); + if (azimuth > pi) + { + azimuth = pi; + } + } + } + else // A indicates Western hemisphere + { + if (dA <= c0) // A altered towards 0 + { + if (azimuth > c0) + { + azimuth = c0; + } + } + else // dA > 0, A altered towards -pi + { + CT const minus_pi = -math::pi<CT>(); + if (azimuth < minus_pi) + { + azimuth = minus_pi; + } + } + } + } }; -}}} // namespace boost::geometry::detail +}}} // namespace boost::geometry::formula -#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_ANDOYER_INVERSE_HPP +#endif // BOOST_GEOMETRY_FORMULAS_ANDOYER_INVERSE_HPP diff --git a/boost/geometry/formulas/differential_quantities.hpp b/boost/geometry/formulas/differential_quantities.hpp new file mode 100644 index 0000000000..9a92f14e18 --- /dev/null +++ b/boost/geometry/formulas/differential_quantities.hpp @@ -0,0 +1,300 @@ +// 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_INVERSE_DIFFERENTIAL_QUANTITIES_HPP +#define BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP + + +#include <boost/geometry/util/condition.hpp> +#include <boost/geometry/util/math.hpp> + + +namespace boost { namespace geometry { namespace formula +{ + +/*! +\brief The solution of a part of the inverse problem - differential quantities. +\author See +- Charles F.F Karney, Algorithms for geodesics, 2011 +https://arxiv.org/pdf/1109.4448.pdf +*/ +template < + typename CT, + bool EnableReducedLength, + bool EnableGeodesicScale, + unsigned int Order = 2, + bool ApproxF = true +> +class differential_quantities +{ +public: + static inline void apply(CT const& lon1, CT const& lat1, + CT const& lon2, CT const& lat2, + CT const& azimuth, CT const& reverse_azimuth, + CT const& b, CT const& f, + CT & reduced_length, CT & geodesic_scale) + { + CT const dlon = lon2 - lon1; + CT const sin_lat1 = sin(lat1); + CT const cos_lat1 = cos(lat1); + CT const sin_lat2 = sin(lat2); + CT const cos_lat2 = cos(lat2); + + apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2, + azimuth, reverse_azimuth, + b, f, + reduced_length, geodesic_scale); + } + + static inline void apply(CT const& dlon, + CT const& sin_lat1, CT const& cos_lat1, + CT const& sin_lat2, CT const& cos_lat2, + CT const& azimuth, CT const& reverse_azimuth, + CT const& b, CT const& f, + CT & reduced_length, CT & geodesic_scale) + { + CT const c0 = 0; + 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; + + // equator + if (math::equals(sin_bet1, c0) && math::equals(sin_bet2, c0)) + { + CT const sig_12 = math::abs(dlon) / one_minus_f; + if (BOOST_GEOMETRY_CONDITION(EnableReducedLength)) + { + CT m12 = sin(sig_12) * b; + reduced_length = m12; + } + + if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale)) + { + CT M12 = cos(sig_12); + geodesic_scale = M12; + } + } + else + { + CT const c2 = 2; + 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 sin_sig1 = sin_bet1; + CT cos_sig1 = cos_alp1 * cos_bet1; + CT sin_sig2 = sin_bet2; + CT cos_sig2 = cos_alp2 * cos_bet2; + + normalize(sin_sig1, cos_sig1); + normalize(sin_sig2, cos_sig2); + + CT const sin_alp0 = sin_alp1 * cos_bet1; + CT const cos_alp0_sqr = c1 - math::sqr(sin_alp0); + + CT const J12 = BOOST_GEOMETRY_CONDITION(ApproxF) ? + 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)); + + if (BOOST_GEOMETRY_CONDITION(EnableReducedLength)) + { + CT const m12_b = dn2 * (cos_sig1 * sin_sig2) + - dn1 * (sin_sig1 * cos_sig2) + - cos_sig1 * cos_sig2 * J12; + CT const m12 = m12_b * b; + + reduced_length = m12; + } + + if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale)) + { + CT const cos_sig12 = cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2; + CT const t = ep2 * (cos_bet1 - cos_bet2) * (cos_bet1 + cos_bet2) / (dn1 + dn2); + CT const M12 = cos_sig12 + (t * sin_sig2 - cos_sig2 * J12) * sin_sig1 / dn1; + + geodesic_scale = M12; + } + } + } + +private: + /*! Approximation of J12, expanded into taylor series in f + Maxima script: + ep2: f * (2-f) / ((1-f)^2); + k2: ca02 * ep2; + assume(f < 1); + assume(sig > 0); + I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig); + I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig); + J(sig):= I1(sig) - I2(sig); + S: taylor(J(sig), f, 0, 3); + S1: factor( 2*integrate(sin(s)^2,s,0,sig)*ca02*f ); + S2: factor( ((integrate(-6*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig)+integrate(-2*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig))*f^2)/4 ); + S3: factor( ((integrate(30*ca02^3*sin(s)^6-54*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig)+integrate(6*ca02^3*sin(s)^6-18*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig))*f^3)/12 ); + */ + static inline CT J12_f(CT const& sin_sig1, CT const& cos_sig1, + CT const& sin_sig2, CT const& cos_sig2, + CT const& cos_alp0_sqr, CT const& f) + { + if (Order == 0) + { + return 0; + } + + CT const c2 = 2; + + CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2, + cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2); + CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1) + CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2) + CT const sin_2sig_12 = sin_2sig2 - sin_2sig1; + CT const L1 = sig_12 - sin_2sig_12 / c2; + + if (Order == 1) + { + return cos_alp0_sqr * f * L1; + } + + CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1) + CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2) + CT const sin_4sig_12 = sin_4sig2 - sin_4sig1; + + CT const c8 = 8; + CT const c12 = 12; + CT const c16 = 16; + CT const c24 = 24; + + CT const L2 = -( cos_alp0_sqr * sin_4sig_12 + + (-c8 * cos_alp0_sqr + c12) * sin_2sig_12 + + (c12 * cos_alp0_sqr - c24) * sig_12) + / c16; + + if (Order == 2) + { + return cos_alp0_sqr * f * (L1 + f * L2); + } + + CT const c4 = 4; + CT const c9 = 9; + CT const c48 = 48; + CT const c60 = 60; + CT const c64 = 64; + CT const c96 = 96; + CT const c128 = 128; + CT const c144 = 144; + + CT const cos_alp0_quad = math::sqr(cos_alp0_sqr); + CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1; + CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2; + CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1; + + CT const A = (c9 * cos_alp0_quad - c12 * cos_alp0_sqr) * sin_4sig_12; + CT const B = c4 * cos_alp0_quad * sin3_2sig_12; + CT const C = (-c48 * cos_alp0_quad + c96 * cos_alp0_sqr - c64) * sin_2sig_12; + CT const D = (c60 * cos_alp0_quad - c144 * cos_alp0_sqr + c128) * sig_12; + + CT const L3 = (A + B + C + D) / c64; + + // Order 3 and higher + return cos_alp0_sqr * f * (L1 + f * (L2 + f * L3)); + } + + /*! Approximation of J12, expanded into taylor series in e'^2 + Maxima script: + k2: ca02 * ep2; + assume(sig > 0); + I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig); + I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig); + J(sig):= I1(sig) - I2(sig); + S: taylor(J(sig), ep2, 0, 3); + S1: factor( integrate(sin(s)^2,s,0,sig)*ca02*ep2 ); + S2: factor( (integrate(sin(s)^4,s,0,sig)*ca02^2*ep2^2)/2 ); + S3: factor( (3*integrate(sin(s)^6,s,0,sig)*ca02^3*ep2^3)/8 ); + */ + static inline CT J12_ep_sqr(CT const& sin_sig1, CT const& cos_sig1, + CT const& sin_sig2, CT const& cos_sig2, + CT const& cos_alp0_sqr, CT const& ep_sqr) + { + if (Order == 0) + { + return 0; + } + + CT const c2 = 2; + CT const c4 = 4; + + CT const c2a0ep2 = cos_alp0_sqr * ep_sqr; + + CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2, + cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2); // sig2 - sig1 + CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1) + CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2) + CT const sin_2sig_12 = sin_2sig2 - sin_2sig1; + + CT const L1 = (c2 * sig_12 - sin_2sig_12) / c4; + + if (Order == 1) + { + return c2a0ep2 * L1; + } + + CT const c8 = 8; + CT const c64 = 64; + + CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1) + CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2) + CT const sin_4sig_12 = sin_4sig2 - sin_4sig1; + + CT const L2 = (sin_4sig_12 - c8 * sin_2sig_12 + 12 * sig_12) / c64; + + if (Order == 2) + { + return c2a0ep2 * (L1 + c2a0ep2 * L2); + } + + CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1; + CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2; + CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1; + + CT const c9 = 9; + CT const c48 = 48; + CT const c60 = 60; + CT const c512 = 512; + + CT const L3 = (c9 * sin_4sig_12 + c4 * sin3_2sig_12 - c48 * sin_2sig_12 + c60 * sig_12) / c512; + + // Order 3 and higher + return c2a0ep2 * (L1 + c2a0ep2 * (L2 + c2a0ep2 * L3)); + } + + static inline void normalize(CT & x, CT & y) + { + CT const len = math::sqrt(math::sqr(x) + math::sqr(y)); + x /= len; + y /= len; + } +}; + +}}} // namespace boost::geometry::formula + + +#endif // BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP diff --git a/boost/geometry/formulas/gnomonic_intersection.hpp b/boost/geometry/formulas/gnomonic_intersection.hpp new file mode 100644 index 0000000000..7e1b7bcdab --- /dev/null +++ b/boost/geometry/formulas/gnomonic_intersection.hpp @@ -0,0 +1,148 @@ +// 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_GNOMONIC_INTERSECTION_HPP +#define BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP + +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/cs.hpp> + +#include <boost/geometry/arithmetic/cross_product.hpp> +#include <boost/geometry/formulas/gnomonic_spheroid.hpp> +#include <boost/geometry/geometries/point.hpp> +#include <boost/geometry/util/math.hpp> + + +namespace boost { namespace geometry { namespace formula +{ + +/*! +\brief The intersection of two geodesics using spheroidal gnomonic projection + as proposed by Karney. +\author See + - Charles F.F Karney, Algorithms for geodesics, 2011 + https://arxiv.org/pdf/1109.4448.pdf + - GeographicLib forum thread: Intersection between two geodesic lines + https://sourceforge.net/p/geographiclib/discussion/1026621/thread/21aaff9f/ +*/ +template +< + typename CT, + template <typename, bool, bool, bool, bool, bool> class Inverse, + template <typename, bool, bool, bool, bool> class Direct +> +class gnomonic_intersection +{ +public: + template <typename T1, typename T2, typename Spheroid> + static inline bool apply(T1 const& lona1, T1 const& lata1, + T1 const& lona2, T1 const& lata2, + T2 const& lonb1, T2 const& latb1, + T2 const& lonb2, T2 const& latb2, + CT & lon, CT & lat, + Spheroid const& spheroid) + { + CT const lon_a1 = lona1; + CT const lat_a1 = lata1; + CT const lon_a2 = lona2; + CT const lat_a2 = lata2; + CT const lon_b1 = lonb1; + CT const lat_b1 = latb1; + CT const lon_b2 = lonb2; + CT const lat_b2 = latb2; + + return apply(lon_a1, lat_a1, lon_a2, lat_a2, lon_b1, lat_b1, lon_b2, lat_b2, lon, lat, spheroid); + } + + template <typename Spheroid> + static inline bool apply(CT const& lona1, CT const& lata1, + CT const& lona2, CT const& lata2, + CT const& lonb1, CT const& latb1, + CT const& lonb2, CT const& latb2, + CT & lon, CT & lat, + Spheroid const& spheroid) + { + typedef gnomonic_spheroid<CT, Inverse, Direct> gnom_t; + + lon = (lona1 + lona2 + lonb1 + lonb2) / 4; + lat = (lata1 + lata2 + latb1 + latb2) / 4; + // TODO: consider normalizing lon + + for (int i = 0; i < 10; ++i) + { + CT xa1, ya1, xa2, ya2; + CT xb1, yb1, xb2, yb2; + CT x, y; + double lat1, lon1; + + bool ok = gnom_t::forward(lon, lat, lona1, lata1, xa1, ya1, spheroid) + && gnom_t::forward(lon, lat, lona2, lata2, xa2, ya2, spheroid) + && gnom_t::forward(lon, lat, lonb1, latb1, xb1, yb1, spheroid) + && gnom_t::forward(lon, lat, lonb2, latb2, xb2, yb2, spheroid) + && intersect(xa1, ya1, xa2, ya2, xb1, yb1, xb2, yb2, x, y) + && gnom_t::inverse(lon, lat, x, y, lon1, lat1, spheroid); + + if (! ok) + { + return false; + } + + if (math::equals(lat1, lat) && math::equals(lon1, lon)) + { + break; + } + + lat = lat1; + lon = lon1; + } + + // NOTE: true is also returned if the number of iterations is too great + // which means that the accuracy of the result is low + return true; + } + +private: + static inline bool intersect(CT const& xa1, CT const& ya1, CT const& xa2, CT const& ya2, + CT const& xb1, CT const& yb1, CT const& xb2, CT const& yb2, + CT & x, CT & y) + { + typedef model::point<CT, 3, cs::cartesian> v3d_t; + + CT const c0 = 0; + CT const c1 = 1; + + v3d_t const va1(xa1, ya1, c1); + v3d_t const va2(xa2, ya2, c1); + v3d_t const vb1(xb1, yb1, c1); + v3d_t const vb2(xb2, yb2, c1); + + v3d_t const la = cross_product(va1, va2); + v3d_t const lb = cross_product(vb1, vb2); + v3d_t const p = cross_product(la, lb); + + CT const z = get<2>(p); + + if (math::equals(z, c0)) + { + // degenerated or collinear segments + return false; + } + + x = get<0>(p) / z; + y = get<1>(p) / z; + + return true; + } +}; + +}}} // namespace boost::geometry::formula + + +#endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_INTERSECTION_HPP diff --git a/boost/geometry/formulas/gnomonic_spheroid.hpp b/boost/geometry/formulas/gnomonic_spheroid.hpp new file mode 100644 index 0000000000..3457397b0f --- /dev/null +++ b/boost/geometry/formulas/gnomonic_spheroid.hpp @@ -0,0 +1,126 @@ +// 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_GNOMONIC_SPHEROID_HPP +#define BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP + + +#include <boost/geometry/core/radius.hpp> + +#include <boost/geometry/algorithms/detail/flattening.hpp> + +#include <boost/geometry/util/condition.hpp> +#include <boost/geometry/util/math.hpp> + +#include <boost/geometry/formulas/andoyer_inverse.hpp> +#include <boost/geometry/formulas/thomas_inverse.hpp> +#include <boost/geometry/formulas/vincenty_direct.hpp> +#include <boost/geometry/formulas/vincenty_inverse.hpp> + + +namespace boost { namespace geometry { namespace formula +{ + +/*! +\brief Gnomonic projection on spheroid (ellipsoid of revolution). +\author See +- Charles F.F Karney, Algorithms for geodesics, 2011 +https://arxiv.org/pdf/1109.4448.pdf +*/ +template < + typename CT, + template <typename, bool, bool, bool, bool ,bool> class Inverse, + template <typename, bool, bool, bool, bool> class Direct +> +class gnomonic_spheroid +{ + typedef Inverse<CT, false, true, true, true, true> inverse_type; + typedef typename inverse_type::result_type inverse_result; + + typedef Direct<CT, false, false, true, true> direct_quantities_type; + typedef Direct<CT, true, false, false, false> direct_coordinates_type; + typedef typename direct_coordinates_type::result_type direct_result; + +public: + template <typename Spheroid> + static inline bool forward(CT const& lon0, CT const& lat0, + CT const& lon, CT const& lat, + CT & x, CT & y, + Spheroid const& spheroid) + { + inverse_result i_res = inverse_type::apply(lon0, lat0, lon, lat, spheroid); + CT const& m = i_res.reduced_length; + CT const& M = i_res.geodesic_scale; + + if (math::smaller_or_equals(M, CT(0))) + { + return false; + } + + CT rho = m / M; + x = sin(i_res.azimuth) * rho; + y = cos(i_res.azimuth) * rho; + + return true; + } + + template <typename Spheroid> + static inline bool inverse(CT const& lon0, CT const& lat0, + CT const& x, CT const& y, + CT & lon, CT & lat, + Spheroid const& spheroid) + { + CT const a = get_radius<0>(spheroid); + CT const ds_threshold = a * std::numeric_limits<CT>::epsilon(); // TODO: 0 for non-fundamental type + + CT const azimuth = atan2(x, y); + CT const rho = math::sqrt(math::sqr(x) + math::sqr(y)); // use hypot? + CT distance = a * atan(rho / a); + + bool found = false; + for (int i = 0 ; i < 10 ; ++i) + { + direct_result d_res = direct_quantities_type::apply(lon0, lat0, distance, azimuth, spheroid); + CT const& m = d_res.reduced_length; + CT const& M = d_res.geodesic_scale; + + if (math::smaller_or_equals(M, CT(0))) + { + // found = false; + return found; + } + + CT const drho = m / M - rho; // rho = m / M + CT const ds = drho * math::sqr(M); // drho/ds = 1/M^2 + distance -= ds; + + // ds_threshold may be 0 + if (math::abs(ds) <= ds_threshold) + { + found = true; + break; + } + } + + if (found) + { + direct_result d_res = direct_coordinates_type::apply(lon0, lat0, distance, azimuth, spheroid); + lon = d_res.lon2; + lat = d_res.lat2; + } + + return found; + } +}; + +}}} // namespace boost::geometry::formula + + +#endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP diff --git a/boost/geometry/formulas/result_direct.hpp b/boost/geometry/formulas/result_direct.hpp new file mode 100644 index 0000000000..8461d8ac9b --- /dev/null +++ b/boost/geometry/formulas/result_direct.hpp @@ -0,0 +1,39 @@ +// 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_RESULT_DIRECT_HPP +#define BOOST_GEOMETRY_FORMULAS_RESULT_DIRECT_HPP + + +namespace boost { namespace geometry { namespace formula +{ + +template <typename T> +struct result_direct +{ + result_direct() + : lon2(0) + , lat2(0) + , reverse_azimuth(0) + , reduced_length(0) + , geodesic_scale(1) + {} + + T lon2; + T lat2; + T reverse_azimuth; + T reduced_length; + T geodesic_scale; +}; + +}}} // namespace boost::geometry::formula + + +#endif // BOOST_GEOMETRY_FORMULAS_RESULT_DIRECT_HPP diff --git a/boost/geometry/formulas/result_inverse.hpp b/boost/geometry/formulas/result_inverse.hpp new file mode 100644 index 0000000000..b6faae6eaa --- /dev/null +++ b/boost/geometry/formulas/result_inverse.hpp @@ -0,0 +1,39 @@ +// Boost.Geometry + +// Copyright (c) 2015-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_RESULT_INVERSE_HPP +#define BOOST_GEOMETRY_FORMULAS_RESULT_INVERSE_HPP + + +namespace boost { namespace geometry { namespace formula +{ + +template <typename T> +struct result_inverse +{ + result_inverse() + : distance(0) + , azimuth(0) + , reverse_azimuth(0) + , reduced_length(0) + , geodesic_scale(1) + {} + + T distance; + T azimuth; + T reverse_azimuth; + T reduced_length; + T geodesic_scale; +}; + +}}} // namespace boost::geometry::formula + + +#endif // BOOST_GEOMETRY_FORMULAS_RESULT_INVERSE_HPP diff --git a/boost/geometry/formulas/sjoberg_intersection.hpp b/boost/geometry/formulas/sjoberg_intersection.hpp new file mode 100644 index 0000000000..03bd4bc97e --- /dev/null +++ b/boost/geometry/formulas/sjoberg_intersection.hpp @@ -0,0 +1,385 @@ +// 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_SJOBERG_INTERSECTION_HPP +#define BOOST_GEOMETRY_FORMULAS_SJOBERG_INTERSECTION_HPP + + +#include <boost/math/constants/constants.hpp> + +#include <boost/geometry/core/radius.hpp> +#include <boost/geometry/core/srs.hpp> + +#include <boost/geometry/util/condition.hpp> +#include <boost/geometry/util/math.hpp> + +#include <boost/geometry/algorithms/detail/flattening.hpp> + + +namespace boost { namespace geometry { namespace formula +{ + +/*! +\brief The intersection of two geodesics as proposed by Sjoberg. +\author 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 +*/ +template +< + typename CT, + template <typename, bool, bool, bool, bool, bool> class Inverse, + unsigned int Order = 4 +> +class sjoberg_intersection +{ + typedef Inverse<CT, false, true, false, false, false> inverse_type; + typedef typename inverse_type::result_type inverse_result; + +public: + template <typename T1, typename T2, typename Spheroid> + static inline bool apply(T1 const& lona1, T1 const& lata1, + T1 const& lona2, T1 const& lata2, + T2 const& lonb1, T2 const& latb1, + T2 const& lonb2, T2 const& latb2, + CT & lon, CT & lat, + Spheroid const& spheroid) + { + CT const lon_a1 = lona1; + CT const lat_a1 = lata1; + CT const lon_a2 = lona2; + CT const lat_a2 = lata2; + CT const lon_b1 = lonb1; + CT const lat_b1 = latb1; + 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; + + return apply(lon_a1, lat_a1, alpha1, lon_b1, lat_b1, alpha2, lon, lat, spheroid); + } + + template <typename Spheroid> + static inline bool apply(CT const& lon1, CT const& lat1, CT const& alpha1, + CT const& lon2, CT const& lat2, CT const& alpha2, + 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>(); + CT const pi_half = pi / c2; + CT const f = detail::flattening<CT>(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) + { + return false; + } + else if (on_equator1) + { + 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); + lat = c0; + lon = lon2 - asin_t2_t02 + dL2; + return true; + } + else if (on_equator2) + { + 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); + 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; + + 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 t = tan_beta1; + int t_id = 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; + + { + 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); + } + + // [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; + + for (int i = 0; i < 10; ++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))) + { + break; + } + + 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; + + 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; + } + + 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; + + if (Order == 1) + { + return -Cj * e_sqr * L0; + } + + CT const c1 = 1; + CT const c16 = 16; + + 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; + + 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 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; + + 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))); + } + + 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)); + } + + template <int TId> + 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) + { + CT const lon1 = lon_a1 + asin_t_t01 - asin_t1_t01; + CT const lon2 = lon_b1 + asin_t_t02 - asin_t2_t02; + + // TODO - true angle difference + CT const dlon = math::abs(lon2 - lon1); + + bool are_equal = math::equals(dlon, CT(0)); + + if ((TId == 0) || are_equal || dlon < current_dlon) + { + current_t = t; + current_lon1 = lon1; + current_lon2 = lon2; + current_dlon = dlon; + t_id = TId; + } + + return are_equal; + } +}; + +}}} // namespace boost::geometry::formula + + +#endif // BOOST_GEOMETRY_FORMULAS_SJOBERG_INTERSECTION_HPP diff --git a/boost/geometry/formulas/thomas_direct.hpp b/boost/geometry/formulas/thomas_direct.hpp new file mode 100644 index 0000000000..f8a7f83943 --- /dev/null +++ b/boost/geometry/formulas/thomas_direct.hpp @@ -0,0 +1,249 @@ +// 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_THOMAS_DIRECT_HPP +#define BOOST_GEOMETRY_FORMULAS_THOMAS_DIRECT_HPP + + +#include <boost/math/constants/constants.hpp> + +#include <boost/geometry/core/radius.hpp> +#include <boost/geometry/core/srs.hpp> + +#include <boost/geometry/util/condition.hpp> +#include <boost/geometry/util/math.hpp> + +#include <boost/geometry/algorithms/detail/flattening.hpp> + +#include <boost/geometry/formulas/differential_quantities.hpp> +#include <boost/geometry/formulas/result_direct.hpp> + + +namespace boost { namespace geometry { namespace formula +{ + + +/*! +\brief The solution of the direct problem of geodesics on latlong coordinates, + Forsyth-Andoyer-Lambert type approximation with second order terms. +\author See + - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965 + http://www.dtic.mil/docs/citations/AD0627893 + - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970 + http://www.dtic.mil/docs/citations/AD0703541 + +*/ +template < + typename CT, + bool EnableCoordinates = true, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false +> +class thomas_direct +{ + static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale; + static const bool CalcCoordinates = EnableCoordinates || CalcQuantities; + static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcCoordinates || CalcQuantities; + +public: + typedef result_direct<CT> result_type; + + template <typename T, typename Dist, typename Azi, typename Spheroid> + static inline result_type apply(T const& lo1, + T const& la1, + Dist const& distance, + Azi const& azimuth12, + Spheroid const& spheroid) + { + result_type result; + + CT const lon1 = lo1; + CT const lat1 = la1; + + if ( math::equals(distance, Dist(0)) || distance < Dist(0) ) + { + result.lon2 = lon1; + result.lat2 = lat1; + return result; + } + + CT const c0 = 0; + CT const c1 = 1; + CT const c2 = 2; + CT const c4 = 4; + + CT const a = CT(get_radius<0>(spheroid)); + CT const b = CT(get_radius<2>(spheroid)); + CT const f = detail::flattening<CT>(spheroid); + CT const one_minus_f = c1 - f; + + CT const pi = math::pi<CT>(); + CT const pi_half = pi / c2; + + // keep azimuth small - experiments show low accuracy + // if the azimuth is closer to (+-)180 deg. + CT azi12_alt = azimuth12; + CT lat1_alt = lat1; + bool alter_result = vflip_if_south(lat1, azimuth12, lat1_alt, azi12_alt); + + CT const theta1 = math::equals(lat1_alt, pi_half) ? lat1_alt : + math::equals(lat1_alt, -pi_half) ? lat1_alt : + atan(one_minus_f * tan(lat1_alt)); + CT const sin_theta1 = sin(theta1); + CT const cos_theta1 = cos(theta1); + + CT const sin_a12 = sin(azi12_alt); + CT const cos_a12 = cos(azi12_alt); + + CT const M = cos_theta1 * sin_a12; // cos_theta0 + CT const theta0 = acos(M); + CT const sin_theta0 = sin(theta0); + + CT const N = cos_theta1 * cos_a12; + CT const C1 = f * M; // lower-case c1 in the technical report + CT const C2 = f * (c1 - math::sqr(M)) / c4; // lower-case c2 in the technical report + CT const D = (c1 - C2) * (c1 - C2 - C1 * M); + CT const P = C2 * (c1 + C1 * M / c2) / D; + + // special case for equator: + // sin_theta0 = 0 <=> lat1 = 0 ^ |azimuth12| = pi/2 + // NOTE: in this case it doesn't matter what's the value of cos_sigma1 because + // theta1=0, theta0=0, M=1|-1, C2=0 so X=0 and Y=0 so d_sigma=d + // cos_a12=0 so N=0, therefore + // lat2=0, azi21=pi/2|-pi/2 + // d_eta = atan2(sin_d_sigma, cos_d_sigma) + // H = C1 * d_sigma + CT const cos_sigma1 = math::equals(sin_theta0, c0) + ? c1 + : normalized1_1(sin_theta1 / sin_theta0); + CT const sigma1 = acos(cos_sigma1); + CT const d = distance / (a * D); + CT const u = 2 * (sigma1 - d); + CT const cos_d = cos(d); + CT const sin_d = sin(d); + CT const cos_u = cos(u); + CT const sin_u = sin(u); + + CT const W = c1 - c2 * P * cos_u; + CT const V = cos_u * cos_d - sin_u * sin_d; + CT const X = math::sqr(C2) * sin_d * cos_d * (2 * math::sqr(V) - c1); + CT const Y = c2 * P * V * W * sin_d; + CT const d_sigma = d + X - Y; + CT const sin_d_sigma = sin(d_sigma); + CT const cos_d_sigma = cos(d_sigma); + + if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth)) + { + result.reverse_azimuth = atan2(M, N * cos_d_sigma - sin_theta1 * sin_d_sigma); + + if (alter_result) + { + vflip_rev_azi(result.reverse_azimuth, azimuth12); + } + } + + if (BOOST_GEOMETRY_CONDITION(CalcCoordinates)) + { + CT const S_sigma = c2 * sigma1 - d_sigma; + CT const cos_S_sigma = cos(S_sigma); + CT const d_eta = atan2(sin_d_sigma * sin_a12, cos_theta1 * cos_d_sigma - sin_theta1 * sin_d_sigma * cos_a12); + CT const H = C1 * (c1 - C2) * d_sigma - C1 * C2 * sin_d_sigma * cos_S_sigma; + CT const d_lambda = d_eta - H; + + result.lon2 = lon1 + d_lambda; + + if (! math::equals(M, c0)) + { + CT const sin_a21 = sin(result.reverse_azimuth); + CT const tan_theta2 = (sin_theta1 * cos_d_sigma + N * sin_d_sigma) * sin_a21 / M; + result.lat2 = atan(tan_theta2 / one_minus_f); + } + else + { + CT const sigma2 = S_sigma - sigma1; + //theta2 = asin(cos(sigma2)) <=> sin_theta0 = 1 + CT const tan_theta2 = cos(sigma2) / sin(sigma2); + result.lat2 = atan(tan_theta2 / one_minus_f); + } + + if (alter_result) + { + result.lat2 = -result.lat2; + } + } + + if (BOOST_GEOMETRY_CONDITION(CalcQuantities)) + { + typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities; + quantities::apply(lon1, lat1, result.lon2, result.lat2, + azimuth12, result.reverse_azimuth, + b, f, + result.reduced_length, result.geodesic_scale); + } + + return result; + } + +private: + static inline bool vflip_if_south(CT const& lat1, CT const& azi12, CT & lat1_alt, CT & azi12_alt) + { + CT const c2 = 2; + CT const pi = math::pi<CT>(); + CT const pi_half = pi / c2; + + if (azi12 > pi_half) + { + azi12_alt = pi - azi12; + lat1_alt = -lat1; + return true; + } + else if (azi12 < -pi_half) + { + azi12_alt = -pi - azi12; + lat1_alt = -lat1; + return true; + } + + return false; + } + + static inline void vflip_rev_azi(CT & rev_azi, CT const& azimuth12) + { + CT const c0 = 0; + CT const pi = math::pi<CT>(); + + if (rev_azi == c0) + { + rev_azi = azimuth12 >= 0 ? pi : -pi; + } + else if (rev_azi > c0) + { + rev_azi = pi - rev_azi; + } + else + { + rev_azi = -pi - rev_azi; + } + } + + static inline CT normalized1_1(CT const& value) + { + CT const c1 = 1; + return value > c1 ? c1 : + value < -c1 ? -c1 : + value; + } +}; + +}}} // namespace boost::geometry::formula + + +#endif // BOOST_GEOMETRY_FORMULAS_THOMAS_DIRECT_HPP diff --git a/boost/geometry/formulas/thomas_inverse.hpp b/boost/geometry/formulas/thomas_inverse.hpp new file mode 100644 index 0000000000..d68c9de054 --- /dev/null +++ b/boost/geometry/formulas/thomas_inverse.hpp @@ -0,0 +1,221 @@ +// Boost.Geometry + +// Copyright (c) 2015-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_THOMAS_INVERSE_HPP +#define BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP + + +#include <boost/math/constants/constants.hpp> + +#include <boost/geometry/core/radius.hpp> +#include <boost/geometry/core/srs.hpp> + +#include <boost/geometry/util/condition.hpp> +#include <boost/geometry/util/math.hpp> + +#include <boost/geometry/algorithms/detail/flattening.hpp> + +#include <boost/geometry/formulas/differential_quantities.hpp> +#include <boost/geometry/formulas/result_inverse.hpp> + + +namespace boost { namespace geometry { namespace formula +{ + +/*! +\brief The solution of the inverse problem of geodesics on latlong coordinates, + Forsyth-Andoyer-Lambert type approximation with second order terms. +\author See + - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965 + http://www.dtic.mil/docs/citations/AD0627893 + - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970 + http://www.dtic.mil/docs/citations/AD0703541 +*/ +template < + typename CT, + bool EnableDistance, + bool EnableAzimuth, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false +> +class thomas_inverse +{ + static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale; + static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities; + static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities; + static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities; + +public: + typedef result_inverse<CT> result_type; + + template <typename T1, typename T2, typename Spheroid> + static inline result_type apply(T1 const& lon1, + T1 const& lat1, + T2 const& lon2, + T2 const& lat2, + Spheroid const& spheroid) + { + result_type result; + + // coordinates in radians + + if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) ) + { + return result; + } + + CT const c0 = 0; + CT const c1 = 1; + CT const c2 = 2; + CT const c4 = 4; + + CT const pi_half = math::pi<CT>() / c2; + CT const f = detail::flattening<CT>(spheroid); + CT const one_minus_f = c1 - f; + +// CT const tan_theta1 = one_minus_f * tan(lat1); +// CT const tan_theta2 = one_minus_f * tan(lat2); +// CT const theta1 = atan(tan_theta1); +// CT const theta2 = atan(tan_theta2); + + CT const theta1 = math::equals(lat1, pi_half) ? lat1 : + math::equals(lat1, -pi_half) ? lat1 : + atan(one_minus_f * tan(lat1)); + CT const theta2 = math::equals(lat2, pi_half) ? lat2 : + math::equals(lat2, -pi_half) ? lat2 : + atan(one_minus_f * tan(lat2)); + + CT const theta_m = (theta1 + theta2) / c2; + CT const d_theta_m = (theta2 - theta1) / c2; + CT const d_lambda = lon2 - lon1; + CT const d_lambda_m = d_lambda / c2; + + CT const sin_theta_m = sin(theta_m); + CT const cos_theta_m = cos(theta_m); + CT const sin_d_theta_m = sin(d_theta_m); + CT const cos_d_theta_m = cos(d_theta_m); + CT const sin2_theta_m = math::sqr(sin_theta_m); + CT const cos2_theta_m = math::sqr(cos_theta_m); + CT const sin2_d_theta_m = math::sqr(sin_d_theta_m); + CT const cos2_d_theta_m = math::sqr(cos_d_theta_m); + CT const sin_d_lambda_m = sin(d_lambda_m); + CT const sin2_d_lambda_m = math::sqr(sin_d_lambda_m); + + CT const H = cos2_theta_m - sin2_d_theta_m; + CT const L = sin2_d_theta_m + H * sin2_d_lambda_m; + CT const cos_d = c1 - c2 * L; + CT const d = acos(cos_d); + CT const sin_d = sin(d); + + CT const one_minus_L = c1 - L; + + if ( math::equals(sin_d, c0) + || math::equals(L, c0) + || math::equals(one_minus_L, c0) ) + { + return result; + } + + CT const U = c2 * sin2_theta_m * cos2_d_theta_m / one_minus_L; + CT const V = c2 * sin2_d_theta_m * cos2_theta_m / L; + CT const X = U + V; + CT const Y = U - V; + CT const T = d / sin_d; + CT const D = c4 * math::sqr(T); + CT const E = c2 * cos_d; + CT const A = D * E; + CT const B = c2 * D; + CT const C = T - (A - E) / c2; + + CT const f_sqr = math::sqr(f); + CT const f_sqr_per_64 = f_sqr / CT(64); + + if ( BOOST_GEOMETRY_CONDITION(EnableDistance) ) + { + CT const n1 = X * (A + C*X); + CT const n2 = Y * (B + E*Y); + CT const n3 = D*X*Y; + + CT const delta1d = f * (T*X-Y) / c4; + CT const delta2d = f_sqr_per_64 * (n1 - n2 + n3); + + CT const a = get_radius<0>(spheroid); + + //result.distance = a * sin_d * (T - delta1d); + result.distance = a * sin_d * (T - delta1d + delta2d); + } + + if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) ) + { + // NOTE: if both cos_latX == 0 then below we'd have 0 * INF + // it's a situation when the endpoints are on the poles +-90 deg + // in this case the azimuth could either be 0 or +-pi + // but above always 0 is returned + + CT const F = c2*Y-E*(c4-X); + CT const M = CT(32)*T-(CT(20)*T-A)*X-(B+c4)*Y; + CT const G = f*T/c2 + f_sqr_per_64 * M; + + // TODO: + // If d_lambda is close to 90 or -90 deg then tan(d_lambda) is big + // and F is small. The result is not accurate. + // In the edge case the result may be 2 orders of magnitude less + // accurate than Andoyer's. + CT const tan_d_lambda = tan(d_lambda); + CT const Q = -(F*G*tan_d_lambda) / c4; + CT const d_lambda_m_p = (d_lambda + Q) / c2; + CT const tan_d_lambda_m_p = tan(d_lambda_m_p); + + CT const v = atan2(cos_d_theta_m, sin_theta_m * tan_d_lambda_m_p); + CT const u = atan2(-sin_d_theta_m, cos_theta_m * tan_d_lambda_m_p); + + CT const pi = math::pi<CT>(); + + if (BOOST_GEOMETRY_CONDITION(EnableAzimuth)) + { + CT alpha1 = v + u; + if (alpha1 > pi) + { + alpha1 -= c2 * pi; + } + + result.azimuth = alpha1; + } + + if (BOOST_GEOMETRY_CONDITION(EnableReverseAzimuth)) + { + CT alpha2 = pi - (v - u); + if (alpha2 > pi) + { + alpha2 -= c2 * pi; + } + + result.reverse_azimuth = alpha2; + } + } + + if (BOOST_GEOMETRY_CONDITION(CalcQuantities)) + { + typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities; + quantities::apply(lon1, lat1, lon2, lat2, + result.azimuth, result.reverse_azimuth, + get_radius<2>(spheroid), f, + result.reduced_length, result.geodesic_scale); + } + + return result; + } +}; + +}}} // namespace boost::geometry::formula + + +#endif // BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP diff --git a/boost/geometry/algorithms/detail/vincenty_direct.hpp b/boost/geometry/formulas/vincenty_direct.hpp index 1c47a0f68d..f3647ff4e6 100644 --- a/boost/geometry/algorithms/detail/vincenty_direct.hpp +++ b/boost/geometry/formulas/vincenty_direct.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014. -// Modifications copyright (c) 2014 Oracle and/or its affiliates. +// This file was modified by Oracle on 2014, 2016. +// Modifications copyright (c) 2014-2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -11,8 +11,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_VINCENTY_DIRECT_HPP -#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_VINCENTY_DIRECT_HPP +#ifndef BOOST_GEOMETRY_FORMULAS_VINCENTY_DIRECT_HPP +#define BOOST_GEOMETRY_FORMULAS_VINCENTY_DIRECT_HPP #include <boost/math/constants/constants.hpp> @@ -20,30 +20,22 @@ #include <boost/geometry/core/radius.hpp> #include <boost/geometry/core/srs.hpp> +#include <boost/geometry/util/condition.hpp> #include <boost/geometry/util/math.hpp> #include <boost/geometry/algorithms/detail/flattening.hpp> +#include <boost/geometry/formulas/differential_quantities.hpp> +#include <boost/geometry/formulas/result_direct.hpp> + #ifndef BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS #define BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS 1000 #endif -namespace boost { namespace geometry { namespace detail -{ - -template <typename T> -struct result_direct +namespace boost { namespace geometry { namespace formula { - void set(T const& lo2, T const& la2) - { - lon2 = lo2; - lat2 = la2; - } - T lon2; - T lat2; -}; /*! \brief The solution of the direct problem of geodesics on latlong coordinates, after Vincenty, 1975 @@ -56,12 +48,22 @@ struct result_direct - http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink */ -template <typename CT> -struct vincenty_direct +template < + typename CT, + bool EnableCoordinates = true, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false +> +class vincenty_direct { - typedef result_direct<CT> result_type; + static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale; + static const bool CalcCoordinates = EnableCoordinates || CalcQuantities; + static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities; public: + typedef result_direct<CT> result_type; + template <typename T, typename Dist, typename Azi, typename Spheroid> static inline result_type apply(T const& lo1, T const& la1, @@ -76,7 +78,8 @@ public: if ( math::equals(distance, Dist(0)) || distance < Dist(0) ) { - result.set(lon1, lat1); + result.lon2 = lon1; + result.lat2 = lat1; return result; } @@ -141,13 +144,12 @@ public: //&& geometry::math::abs(sigma) < pi && counter < BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS ); // robustness + if (BOOST_GEOMETRY_CONDITION(CalcCoordinates)) { result.lat2 = atan2( sin_U1 * cos_sigma + cos_U1 * sin_sigma * cos_azimuth12, one_min_f * math::sqrt(sin_alpha_sqr + math::sqr(sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_azimuth12))); // (8) - } - - { + CT const lambda = atan2( sin_sigma * sin_azimuth12, cos_U1 * cos_sigma - sin_U1 * sin_sigma * cos_azimuth12); // (9) CT const C = (flattening/CT(16)) * cos_alpha_sqr * ( CT(4) + flattening * ( CT(4) - CT(3) * cos_alpha_sqr ) ); // (10) @@ -157,21 +159,27 @@ public: result.lon2 = lon1 + L; } + if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth)) + { + result.reverse_azimuth + = atan2(sin_alpha, -sin_U1 * sin_sigma + cos_U1 * cos_sigma * cos_azimuth12); // (12) + } + + if (BOOST_GEOMETRY_CONDITION(CalcQuantities)) + { + typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities; + quantities::apply(lon1, lat1, result.lon2, result.lat2, + azimuth12, result.reverse_azimuth, + radius_b, flattening, + result.reduced_length, result.geodesic_scale); + } + return result; } - /* - inline CT azimuth21() const - { - // NOTE: signs of X and Y are different than in the original paper - return is_distance_zero ? - CT(0) : - atan2(-sin_alpha, sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_azimuth12); // (12) - } - */ }; -}}} // namespace boost::geometry::detail +}}} // namespace boost::geometry::formula -#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_VINCENTY_DIRECT_HPP +#endif // BOOST_GEOMETRY_FORMULAS_VINCENTY_DIRECT_HPP diff --git a/boost/geometry/algorithms/detail/vincenty_inverse.hpp b/boost/geometry/formulas/vincenty_inverse.hpp index fe05e95932..bbda00036b 100644 --- a/boost/geometry/algorithms/detail/vincenty_inverse.hpp +++ b/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. -// Modifications copyright (c) 2014 Oracle and/or its affiliates. +// This file was modified by Oracle on 2014, 2016. +// Modifications copyright (c) 2014-2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -11,8 +11,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_VINCENTY_INVERSE_HPP -#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_VINCENTY_INVERSE_HPP +#ifndef BOOST_GEOMETRY_FORMULAS_VINCENTY_INVERSE_HPP +#define BOOST_GEOMETRY_FORMULAS_VINCENTY_INVERSE_HPP #include <boost/math/constants/constants.hpp> @@ -24,7 +24,9 @@ #include <boost/geometry/util/math.hpp> #include <boost/geometry/algorithms/detail/flattening.hpp> -#include <boost/geometry/algorithms/detail/result_inverse.hpp> + +#include <boost/geometry/formulas/differential_quantities.hpp> +#include <boost/geometry/formulas/result_inverse.hpp> #ifndef BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS @@ -32,7 +34,7 @@ #endif -namespace boost { namespace geometry { namespace detail +namespace boost { namespace geometry { namespace formula { /*! @@ -46,12 +48,24 @@ namespace boost { namespace geometry { namespace detail - http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink */ -template <typename CT, bool EnableDistance, bool EnableAzimuth> +template < + typename CT, + bool EnableDistance, + bool EnableAzimuth, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false +> struct vincenty_inverse { - typedef result_inverse<CT> result_type; + static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale; + static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities; + static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities; + static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities; public: + typedef result_inverse<CT> result_type; + template <typename T1, typename T2, typename Spheroid> static inline result_type apply(T1 const& lon1, T1 const& lat1, @@ -63,7 +77,6 @@ public: if (math::equals(lat1, lat2) && math::equals(lon1, lon2)) { - result.set(CT(0), CT(0)); return result; } @@ -174,31 +187,34 @@ public: result.distance = radius_b * A * (sigma - delta_sigma); // (19) } - else - { - result.distance = CT(0); - } - if ( BOOST_GEOMETRY_CONDITION(EnableAzimuth) ) + if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) ) { - result.azimuth = atan2(cos_U2 * sin_lambda, cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda); // (20) + if (BOOST_GEOMETRY_CONDITION(CalcFwdAzimuth)) + { + result.azimuth = atan2(cos_U2 * sin_lambda, cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda); // (20) + } + + if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth)) + { + result.reverse_azimuth = atan2(cos_U1 * sin_lambda, -sin_U1 * cos_U2 + cos_U1 * sin_U2 * cos_lambda); // (21) + } } - else + + if (BOOST_GEOMETRY_CONDITION(CalcQuantities)) { - result.azimuth = CT(0); + typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities; + quantities::apply(lon1, lat1, lon2, lat2, + result.azimuth, result.reverse_azimuth, + radius_b, flattening, + result.reduced_length, result.geodesic_scale); } return result; } - -// inline CT azimuth21() const -// { -// // NOTE: signs of X and Y are different than in the original paper -// atan2(-cos_U1 * sin_lambda, sin_U1 * cos_U2 - cos_U1 * sin_U2 * cos_lambda); // (21) -// } }; -}}} // namespace boost::geometry::detail +}}} // namespace boost::geometry::formula -#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_VINCENTY_INVERSE_HPP +#endif // BOOST_GEOMETRY_FORMULAS_VINCENTY_INVERSE_HPP diff --git a/boost/geometry/geometries/adapted/std_array.hpp b/boost/geometry/geometries/adapted/std_array.hpp new file mode 100644 index 0000000000..4f5cbe0d32 --- /dev/null +++ b/boost/geometry/geometries/adapted/std_array.hpp @@ -0,0 +1,115 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2010 Alfredo Correa +// Copyright (c) 2010-2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2016 Norbert Wenzel + +// 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_GEOMETRIES_ADAPTED_STD_ARRAY_HPP +#define BOOST_GEOMETRY_GEOMETRIES_ADAPTED_STD_ARRAY_HPP + + +#define BOOST_GEOMETRY_ADAPTED_STD_ARRAY_TAG_DEFINED + + +#include <cstddef> + +#include <boost/type_traits/is_arithmetic.hpp> + +#include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/cs.hpp> +#include <boost/geometry/core/coordinate_dimension.hpp> +#include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/core/tags.hpp> + +#include <array> + +namespace boost { namespace geometry +{ + + +#ifndef DOXYGEN_NO_TRAITS_SPECIALIZATIONS +namespace traits +{ + + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + + +// Create class and specialization to indicate the tag +// for normal cases and the case that the type of the std-array is arithmetic +template <bool> +struct std_array_tag +{ + typedef geometry_not_recognized_tag type; +}; + + +template <> +struct std_array_tag<true> +{ + typedef point_tag type; +}; + + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + + +// Assign the point-tag, preventing arrays of points getting a point-tag +template <typename CoordinateType, std::size_t DimensionCount> +struct tag<std::array<CoordinateType, DimensionCount> > + : detail::std_array_tag<boost::is_arithmetic<CoordinateType>::value> {}; + + +template <typename CoordinateType, std::size_t DimensionCount> +struct coordinate_type<std::array<CoordinateType, DimensionCount> > +{ + typedef CoordinateType type; +}; + + +template <typename CoordinateType, std::size_t DimensionCount> +struct dimension<std::array<CoordinateType, DimensionCount> >: boost::mpl::int_<DimensionCount> {}; + + +template <typename CoordinateType, std::size_t DimensionCount, std::size_t Dimension> +struct access<std::array<CoordinateType, DimensionCount>, Dimension> +{ + static inline CoordinateType get(std::array<CoordinateType, DimensionCount> const& a) + { + return a[Dimension]; + } + + static inline void set(std::array<CoordinateType, DimensionCount>& a, + CoordinateType const& value) + { + a[Dimension] = value; + } +}; + + +} // namespace traits +#endif // DOXYGEN_NO_TRAITS_SPECIALIZATIONS + + +}} // namespace boost::geometry + + +#define BOOST_GEOMETRY_REGISTER_STD_ARRAY_CS(CoordinateSystem) \ + namespace boost { namespace geometry { namespace traits { \ + template <class T, std::size_t N> \ + struct coordinate_system<std::array<T, N> > \ + { \ + typedef CoordinateSystem type; \ + }; \ + }}} + + +#endif // BOOST_GEOMETRY_GEOMETRIES_ADAPTED_STD_ARRAY_HPP + diff --git a/boost/geometry/index/detail/predicates.hpp b/boost/geometry/index/detail/predicates.hpp index 227939c96d..72c6c661b7 100644 --- a/boost/geometry/index/detail/predicates.hpp +++ b/boost/geometry/index/detail/predicates.hpp @@ -11,7 +11,11 @@ #ifndef BOOST_GEOMETRY_INDEX_DETAIL_PREDICATES_HPP #define BOOST_GEOMETRY_INDEX_DETAIL_PREDICATES_HPP -#include <boost/geometry/index/predicates.hpp> +//#include <utility> + +#include <boost/mpl/assert.hpp> +#include <boost/tuple/tuple.hpp> + #include <boost/geometry/index/detail/tags.hpp> namespace boost { namespace geometry { namespace index { namespace detail { diff --git a/boost/geometry/index/predicates.hpp b/boost/geometry/index/predicates.hpp index 3bb1bf4d87..f4e22bdd55 100644 --- a/boost/geometry/index/predicates.hpp +++ b/boost/geometry/index/predicates.hpp @@ -11,10 +11,6 @@ #ifndef BOOST_GEOMETRY_INDEX_PREDICATES_HPP #define BOOST_GEOMETRY_INDEX_PREDICATES_HPP -#include <utility> -#include <boost/tuple/tuple.hpp> -#include <boost/mpl/assert.hpp> - #include <boost/geometry/index/detail/predicates.hpp> #include <boost/geometry/index/detail/tuples.hpp> diff --git a/boost/geometry/strategies/cartesian/box_in_box.hpp b/boost/geometry/strategies/cartesian/box_in_box.hpp index fe77de9e2d..28a6f29336 100644 --- a/boost/geometry/strategies/cartesian/box_in_box.hpp +++ b/boost/geometry/strategies/cartesian/box_in_box.hpp @@ -36,8 +36,7 @@ namespace within { -template <typename Geometry, std::size_t Dimension, typename CSTag> -struct box_within_range +struct box_within_coord { template <typename BoxContainedValue, typename BoxContainingValue> static inline bool apply(BoxContainedValue const& bed_min, @@ -51,8 +50,7 @@ struct box_within_range }; -template <typename Geometry, std::size_t Dimension, typename CSTag> -struct box_covered_by_range +struct box_covered_by_coord { template <typename BoxContainedValue, typename BoxContainingValue> static inline bool apply(BoxContainedValue const& bed_min, @@ -65,7 +63,19 @@ struct box_covered_by_range }; -struct box_within_longitude_check +template <typename Geometry, std::size_t Dimension, typename CSTag> +struct box_within_range + : box_within_coord +{}; + + +template <typename Geometry, std::size_t Dimension, typename CSTag> +struct box_covered_by_range + : box_covered_by_coord +{}; + + +struct box_within_longitude_diff { template <typename CalcT> static inline bool apply(CalcT const& diff_ed) @@ -74,7 +84,7 @@ struct box_within_longitude_check } }; -struct box_covered_by_longitude_check +struct box_covered_by_longitude_diff { template <typename CalcT> static inline bool apply(CalcT const&) @@ -84,6 +94,7 @@ struct box_covered_by_longitude_check }; template <typename Geometry, + typename CoordCheck, typename InteriorCheck> struct box_longitude_range { @@ -101,6 +112,11 @@ struct box_longitude_range typedef typename coordinate_system<Geometry>::type::units units_t; typedef math::detail::constants_on_spheroid<calc_t, units_t> constants; + if (CoordCheck::apply(bed_min, bed_max, bing_min, bing_max)) + { + return true; + } + // min <= max <=> diff >= 0 calc_t const diff_ed = bed_max - bed_min; calc_t const diff_ing = bing_max - bing_min; @@ -122,7 +138,8 @@ struct box_longitude_range calc_t const diff_min = math::longitude_distance_unsigned<units_t>(bing_min, bed_min); // max of contained translated into the containing origin must be lesser than max of containing - return bing_min + diff_min + diff_ed <= bing_max; + return bing_min + diff_min + diff_ed <= bing_max + /*|| bing_max - diff_min - diff_ed >= bing_min*/; } }; @@ -130,13 +147,13 @@ struct box_longitude_range // spherical_equatorial_tag, spherical_polar_tag and geographic_cat are casted to spherical_tag template <typename Geometry> struct box_within_range<Geometry, 0, spherical_tag> - : box_longitude_range<Geometry, box_within_longitude_check> + : box_longitude_range<Geometry, box_within_coord, box_within_longitude_diff> {}; template <typename Geometry> struct box_covered_by_range<Geometry, 0, spherical_tag> - : box_longitude_range<Geometry, box_covered_by_longitude_check> + : box_longitude_range<Geometry, box_covered_by_coord, box_covered_by_longitude_diff> {}; diff --git a/boost/geometry/strategies/cartesian/point_in_box.hpp b/boost/geometry/strategies/cartesian/point_in_box.hpp index 28ed8215e7..227a98f2ad 100644 --- a/boost/geometry/strategies/cartesian/point_in_box.hpp +++ b/boost/geometry/strategies/cartesian/point_in_box.hpp @@ -33,9 +33,7 @@ namespace boost { namespace geometry { namespace strategy namespace within { - -template <typename Geometry, std::size_t Dimension, typename CSTag> -struct within_range +struct within_coord { template <typename Value1, typename Value2> static inline bool apply(Value1 const& value, Value2 const& min_value, Value2 const& max_value) @@ -44,9 +42,7 @@ struct within_range } }; - -template <typename Geometry, std::size_t Dimension, typename CSTag> -struct covered_by_range +struct covered_by_coord { template <typename Value1, typename Value2> static inline bool apply(Value1 const& value, Value2 const& min_value, Value2 const& max_value) @@ -55,32 +51,47 @@ struct covered_by_range } }; +template <typename Geometry, std::size_t Dimension, typename CSTag> +struct within_range + : within_coord +{}; + + +template <typename Geometry, std::size_t Dimension, typename CSTag> +struct covered_by_range + : covered_by_coord +{}; + // NOTE: the result would be the same if instead of structs defined below // the above xxx_range were used with the following arguments: // (min_value + diff_min, min_value, max_value) -struct within_longitude_range +struct within_longitude_diff { template <typename CalcT> static inline bool apply(CalcT const& diff_min, CalcT const& min_value, CalcT const& max_value) { CalcT const c0 = 0; - return diff_min > c0 && min_value + diff_min < max_value; + return diff_min > c0 + && (min_value + diff_min < max_value + /*|| max_value - diff_min > min_value*/); } }; -struct covered_by_longitude_range +struct covered_by_longitude_diff { template <typename CalcT> static inline bool apply(CalcT const& diff_min, CalcT const& min_value, CalcT const& max_value) { - return min_value + diff_min <= max_value; + return min_value + diff_min <= max_value + /*|| max_value - diff_min >= min_value*/; } }; template <typename Geometry, - typename ResultCheck> + typename CoordCheck, + typename DiffCheck> struct longitude_range { template <typename Value1, typename Value2> @@ -93,6 +104,11 @@ struct longitude_range typedef typename coordinate_system<Geometry>::type::units units_t; typedef math::detail::constants_on_spheroid<calc_t, units_t> constants; + if (CoordCheck::apply(value, min_value, max_value)) + { + return true; + } + // min <= max <=> diff >= 0 calc_t const diff_ing = max_value - min_value; @@ -105,7 +121,7 @@ struct longitude_range // calculate positive longitude translation with min_value as origin calc_t const diff_min = math::longitude_distance_unsigned<units_t, calc_t>(min_value, value); - return ResultCheck::template apply<calc_t>(diff_min, min_value, max_value); + return DiffCheck::template apply<calc_t>(diff_min, min_value, max_value); } }; @@ -113,13 +129,13 @@ struct longitude_range // spherical_equatorial_tag, spherical_polar_tag and geographic_cat are casted to spherical_tag template <typename Geometry> struct within_range<Geometry, 0, spherical_tag> - : longitude_range<Geometry, within_longitude_range> + : longitude_range<Geometry, within_coord, within_longitude_diff> {}; template <typename Geometry> struct covered_by_range<Geometry, 0, spherical_tag> - : longitude_range<Geometry, covered_by_longitude_range> + : longitude_range<Geometry, covered_by_coord, covered_by_longitude_diff> {}; diff --git a/boost/geometry/strategies/geographic/distance_andoyer.hpp b/boost/geometry/strategies/geographic/distance_andoyer.hpp index 1646727d09..1946cd1090 100644 --- a/boost/geometry/strategies/geographic/distance_andoyer.hpp +++ b/boost/geometry/strategies/geographic/distance_andoyer.hpp @@ -20,9 +20,10 @@ #include <boost/geometry/core/radius.hpp> #include <boost/geometry/core/srs.hpp> -#include <boost/geometry/algorithms/detail/andoyer_inverse.hpp> #include <boost/geometry/algorithms/detail/flattening.hpp> +#include <boost/geometry/formulas/andoyer_inverse.hpp> + #include <boost/geometry/strategies/distance.hpp> #include <boost/geometry/util/math.hpp> @@ -89,7 +90,7 @@ public : inline typename calculation_type<Point1, Point2>::type apply(Point1 const& point1, Point2 const& point2) const { - return geometry::detail::andoyer_inverse + return geometry::formula::andoyer_inverse < typename calculation_type<Point1, Point2>::type, true, false diff --git a/boost/geometry/strategies/geographic/distance_thomas.hpp b/boost/geometry/strategies/geographic/distance_thomas.hpp index 7252b723dd..39e0ecfa6f 100644 --- a/boost/geometry/strategies/geographic/distance_thomas.hpp +++ b/boost/geometry/strategies/geographic/distance_thomas.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2015. -// Modifications copyright (c) 2015 Oracle and/or its affiliates. +// This file was modified by Oracle on 2015, 2016. +// Modifications copyright (c) 2015-2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -23,7 +23,7 @@ #include <boost/geometry/util/promote_floating_point.hpp> #include <boost/geometry/util/select_calculation_type.hpp> -#include <boost/geometry/algorithms/detail/thomas_inverse.hpp> +#include <boost/geometry/formulas/thomas_inverse.hpp> namespace boost { namespace geometry { @@ -78,7 +78,7 @@ public : inline typename calculation_type<Point1, Point2>::type apply(Point1 const& point1, Point2 const& point2) const { - return geometry::detail::thomas_inverse + return geometry::formula::thomas_inverse < typename calculation_type<Point1, Point2>::type, true, false diff --git a/boost/geometry/strategies/geographic/distance_vincenty.hpp b/boost/geometry/strategies/geographic/distance_vincenty.hpp index 65bfa19939..e79e9aeb46 100644 --- a/boost/geometry/strategies/geographic/distance_vincenty.hpp +++ b/boost/geometry/strategies/geographic/distance_vincenty.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014. -// Modifications copyright (c) 2014 Oracle and/or its affiliates. +// This file was modified by Oracle on 2014, 2016. +// Modifications copyright (c) 2014-2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -23,7 +23,7 @@ #include <boost/geometry/util/promote_floating_point.hpp> #include <boost/geometry/util/select_calculation_type.hpp> -#include <boost/geometry/algorithms/detail/vincenty_inverse.hpp> +#include <boost/geometry/formulas/vincenty_inverse.hpp> namespace boost { namespace geometry { @@ -80,7 +80,7 @@ public : inline typename calculation_type<Point1, Point2>::type apply(Point1 const& point1, Point2 const& point2) const { - return geometry::detail::vincenty_inverse + return geometry::formula::vincenty_inverse < typename calculation_type<Point1, Point2>::type, true, false diff --git a/boost/geometry/strategies/geographic/side_andoyer.hpp b/boost/geometry/strategies/geographic/side_andoyer.hpp index e0f0c04067..c3e71cd1cd 100644 --- a/boost/geometry/strategies/geographic/side_andoyer.hpp +++ b/boost/geometry/strategies/geographic/side_andoyer.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014, 2015. -// Modifications copyright (c) 2014-2015 Oracle and/or its affiliates. +// 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 @@ -15,7 +15,7 @@ #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_ANDOYER_HPP -#include <boost/geometry/algorithms/detail/andoyer_inverse.hpp> +#include <boost/geometry/formulas/andoyer_inverse.hpp> #include <boost/geometry/strategies/geographic/side_detail.hpp> @@ -36,9 +36,9 @@ namespace strategy { namespace side */ template <typename Model, typename CalculationType = void> class andoyer - : public detail::by_azimuth<geometry::detail::andoyer_inverse, Model, CalculationType> + : public detail::by_azimuth<geometry::formula::andoyer_inverse, Model, CalculationType> { - typedef detail::by_azimuth<geometry::detail::andoyer_inverse, Model, CalculationType> base_t; + typedef detail::by_azimuth<geometry::formula::andoyer_inverse, Model, CalculationType> base_t; public: andoyer(Model const& model = Model()) diff --git a/boost/geometry/strategies/geographic/side_detail.hpp b/boost/geometry/strategies/geographic/side_detail.hpp index c00213d072..ce1b47c88e 100644 --- a/boost/geometry/strategies/geographic/side_detail.hpp +++ b/boost/geometry/strategies/geographic/side_detail.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// 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, 2016. +// Modifications copyright (c) 2014-2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -46,7 +46,7 @@ namespace detail \tparam Model Reference model of coordinate system. \tparam CalculationType \tparam_calculation */ -template <template<typename, bool, bool> class InverseFormula, +template <template<typename, bool, bool, bool, bool, bool> class InverseFormula, typename Model, typename CalculationType = void> class by_azimuth @@ -68,7 +68,7 @@ public: >::type >::type calc_t; - typedef InverseFormula<calc_t, false, true> inverse_formula; + typedef InverseFormula<calc_t, false, true, false, false, false> inverse_formula; calc_t a1p = azimuth<calc_t, inverse_formula>(p1, p, m_model); calc_t a12 = azimuth<calc_t, inverse_formula>(p1, p2, m_model); diff --git a/boost/geometry/strategies/geographic/side_thomas.hpp b/boost/geometry/strategies/geographic/side_thomas.hpp index a96cabdaab..96b0323307 100644 --- a/boost/geometry/strategies/geographic/side_thomas.hpp +++ b/boost/geometry/strategies/geographic/side_thomas.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014, 2015. -// Modifications copyright (c) 2014-2015 Oracle and/or its affiliates. +// 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 @@ -15,7 +15,7 @@ #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_THOMAS_HPP -#include <boost/geometry/algorithms/detail/thomas_inverse.hpp> +#include <boost/geometry/formulas/thomas_inverse.hpp> #include <boost/geometry/strategies/geographic/side_detail.hpp> @@ -36,9 +36,9 @@ namespace strategy { namespace side */ template <typename Model, typename CalculationType = void> class thomas - : public detail::by_azimuth<geometry::detail::thomas_inverse, Model, CalculationType> + : public detail::by_azimuth<geometry::formula::thomas_inverse, Model, CalculationType> { - typedef detail::by_azimuth<geometry::detail::thomas_inverse, Model, CalculationType> base_t; + typedef detail::by_azimuth<geometry::formula::thomas_inverse, Model, CalculationType> base_t; public: thomas(Model const& model = Model()) diff --git a/boost/geometry/strategies/geographic/side_vincenty.hpp b/boost/geometry/strategies/geographic/side_vincenty.hpp index a8ef9550d6..103277a8bd 100644 --- a/boost/geometry/strategies/geographic/side_vincenty.hpp +++ b/boost/geometry/strategies/geographic/side_vincenty.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014, 2015. -// Modifications copyright (c) 2014-2015 Oracle and/or its affiliates. +// 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 @@ -15,7 +15,7 @@ #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_SIDE_VINCENTY_HPP -#include <boost/geometry/algorithms/detail/vincenty_inverse.hpp> +#include <boost/geometry/formulas/vincenty_inverse.hpp> #include <boost/geometry/strategies/geographic/side_detail.hpp> @@ -36,9 +36,9 @@ namespace strategy { namespace side */ template <typename Model, typename CalculationType = void> class vincenty - : public detail::by_azimuth<geometry::detail::vincenty_inverse, Model, CalculationType> + : public detail::by_azimuth<geometry::formula::vincenty_inverse, Model, CalculationType> { - typedef detail::by_azimuth<geometry::detail::vincenty_inverse, Model, CalculationType> base_t; + typedef detail::by_azimuth<geometry::formula::vincenty_inverse, Model, CalculationType> base_t; public: vincenty(Model const& model = Model()) |