diff options
Diffstat (limited to 'boost/geometry')
47 files changed, 2134 insertions, 386 deletions
diff --git a/boost/geometry/algorithms/area.hpp b/boost/geometry/algorithms/area.hpp index 919ab66fe1..8193200ab9 100644 --- a/boost/geometry/algorithms/area.hpp +++ b/boost/geometry/algorithms/area.hpp @@ -92,7 +92,7 @@ struct ring_area // An open ring has at least three points, // A closed ring has at least four points, // if not, there is no (zero) area - if (boost::size(ring) + if (int(boost::size(ring)) < core_detail::closure::minimum_ring_size<Closure>::value) { return type(); diff --git a/boost/geometry/algorithms/detail/for_each_range.hpp b/boost/geometry/algorithms/detail/for_each_range.hpp index b10017c209..7cb01fa9b4 100644 --- a/boost/geometry/algorithms/detail/for_each_range.hpp +++ b/boost/geometry/algorithms/detail/for_each_range.hpp @@ -22,6 +22,7 @@ #include <boost/geometry/core/tag_cast.hpp> #include <boost/geometry/util/add_const_if_c.hpp> +#include <boost/geometry/views/box_view.hpp> namespace boost { namespace geometry diff --git a/boost/geometry/algorithms/detail/get_left_turns.hpp b/boost/geometry/algorithms/detail/get_left_turns.hpp new file mode 100644 index 0000000000..62f0f7f0f4 --- /dev/null +++ b/boost/geometry/algorithms/detail/get_left_turns.hpp @@ -0,0 +1,367 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2012 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_GET_LEFT_TURNS_HPP +#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_GET_LEFT_TURNS_HPP + +#include <boost/geometry/iterators/ever_circling_iterator.hpp> + +namespace boost { namespace geometry +{ + + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +// TODO: move this to /util/ +template <typename T> +static inline std::pair<T, T> ordered_pair(T const& first, T const& second) +{ + return first < second ? std::make_pair(first, second) : std::make_pair(second, first); +} + +template <typename AngleInfo> +inline void debug_left_turn(AngleInfo const& ai, AngleInfo const& previous) +{ +#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION + std::cout << "Angle: " << (ai.incoming ? "i" : "o") + << " " << si(ai.seg_id) + << " " << (math::r2d * (ai.angle) ) + << " turn: " << ai.turn_index << "[" << ai.operation_index << "]" + ; + + if (ai.turn_index != previous.turn_index + || ai.operation_index != previous.operation_index) + { + std::cout << " diff: " << math::r2d * math::abs(previous.angle - ai.angle); + } + std::cout << std::endl; +#endif +} + +template <typename AngleInfo> +inline void debug_left_turn(std::string const& caption, AngleInfo const& ai, AngleInfo const& previous, + int code = -99, int code2 = -99, int code3 = -99, int code4 = -99) +{ +#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION + std::cout << " " << caption + << " turn: " << ai.turn_index << "[" << ai.operation_index << "]" + << " " << si(ai.seg_id) + << " " << (ai.incoming ? "i" : "o") + << " " << (math::r2d * (ai.angle) ) + << " turn: " << previous.turn_index << "[" << previous.operation_index << "]" + << " " << si(previous.seg_id) + << " " << (previous.incoming ? "i" : "o") + << " " << (math::r2d * (previous.angle) ) + ; + + if (code != -99) + { + std::cout << " code: " << code << " , " << code2 << " , " << code3 << " , " << code4; + } + std::cout << std::endl; +#endif +} + + +template <typename Operation> +inline bool include_operation(Operation const& op, + segment_identifier const& outgoing_seg_id, + segment_identifier const& incoming_seg_id) +{ + return op.seg_id == outgoing_seg_id + && op.other_id == incoming_seg_id + && (op.operation == detail::overlay::operation_union + ||op.operation == detail::overlay::operation_continue) + ; +} + +template <typename Turn> +inline bool process_include(segment_identifier const& outgoing_seg_id, segment_identifier const& incoming_seg_id, + int turn_index, Turn& turn, + std::set<int>& keep_indices, int priority) +{ + bool result = false; + for (int i = 0; i < 2; i++) + { + if (include_operation(turn.operations[i], outgoing_seg_id, incoming_seg_id)) + { + turn.operations[i].include_in_occupation_map = true; + if (priority > turn.priority) + { + turn.priority = priority; + } + keep_indices.insert(turn_index); + result = true; + } + } + return result; +} + +template <typename AngleInfo, typename Turns, typename TurnSegmentIndices> +inline bool include_left_turn_of_all( + AngleInfo const& outgoing, AngleInfo const& incoming, + Turns& turns, TurnSegmentIndices const& turn_segment_indices, + std::set<int>& keep_indices, int priority) +{ + segment_identifier const& outgoing_seg_id = turns[outgoing.turn_index].operations[outgoing.operation_index].seg_id; + segment_identifier const& incoming_seg_id = turns[incoming.turn_index].operations[incoming.operation_index].seg_id; + + if (outgoing.turn_index == incoming.turn_index) + { + return process_include(outgoing_seg_id, incoming_seg_id, outgoing.turn_index, turns[outgoing.turn_index], keep_indices, priority); + } + + bool result = false; + std::pair<segment_identifier, segment_identifier> pair = ordered_pair(outgoing_seg_id, incoming_seg_id); + auto it = turn_segment_indices.find(pair); + if (it != turn_segment_indices.end()) + { + for (auto sit = it->second.begin(); sit != it->second.end(); ++sit) + { + if (process_include(outgoing_seg_id, incoming_seg_id, *sit, turns[*sit], keep_indices, priority)) + { + result = true; + } + } + } + return result; +} + +template <std::size_t Index, typename Turn> +inline bool corresponds(Turn const& turn, segment_identifier const& seg_id) +{ + return turn.operations[Index].operation == detail::overlay::operation_union + && turn.operations[Index].seg_id == seg_id; +} + + +template <typename Turns, typename TurnSegmentIndices> +inline bool prefer_by_other(Turns const& turns, + TurnSegmentIndices const& turn_segment_indices, + std::set<int>& indices) +{ + std::map<segment_identifier, int> map; + for (auto sit = indices.begin(); sit != indices.end(); ++sit) + { + map[turns[*sit].operations[0].seg_id]++; + map[turns[*sit].operations[1].seg_id]++; + } + + std::set<segment_identifier> segment_occuring_once; + for (auto mit = map.begin(); mit != map.end(); ++mit) + { + if (mit->second == 1) + { + segment_occuring_once.insert(mit->first); + } +#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_PREFER + std::cout << si(mit->first) << " " << mit->second << std::endl; +#endif + } + + if (segment_occuring_once.size() == 2) + { + // Try to find within all turns a turn with these two segments + + std::set<segment_identifier>::const_iterator soo_it = segment_occuring_once.begin(); + segment_identifier front = *soo_it; + soo_it++; + segment_identifier back = *soo_it; + + std::pair<segment_identifier, segment_identifier> pair = ordered_pair(front, back); + auto it = turn_segment_indices.find(pair); + if (it != turn_segment_indices.end()) + { + // debug_turns_by_indices("Found", it->second); + // Check which is the union/continue + segment_identifier good; + for (auto sit = it->second.begin(); sit != it->second.end(); ++sit) + { + if (turns[*sit].operations[0].operation == detail::overlay::operation_union) + { + good = turns[*sit].operations[0].seg_id; + } + else if (turns[*sit].operations[1].operation == detail::overlay::operation_union) + { + good = turns[*sit].operations[1].seg_id; + } + } + +#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_PREFER + std::cout << "Good: " << si(good) << std::endl; +#endif + + // Find in indexes-to-keep this segment with the union. Discard the other one + std::set<int> ok_indices; + for (auto sit = indices.begin(); sit != indices.end(); ++sit) + { + if (corresponds<0>(turns[*sit], good) || corresponds<1>(turns[*sit], good)) + { + ok_indices.insert(*sit); + } + } + if (ok_indices.size() > 0 && ok_indices.size() < indices.size()) + { + indices = ok_indices; + std::cout << "^"; + return true; + } + } + } + return false; +} + +template <typename Turns> +inline void prefer_by_priority(Turns const& turns, std::set<int>& indices) +{ + // Find max prio + int min_prio = 1024, max_prio = 0; + for (auto sit = indices.begin(); sit != indices.end(); ++sit) + { + if (turns[*sit].priority > max_prio) + { + max_prio = turns[*sit].priority; + } + if (turns[*sit].priority < min_prio) + { + min_prio = turns[*sit].priority; + } + } + + if (min_prio == max_prio) + { + return; + } + + // Only keep indices with high prio + std::set<int> ok_indices; + for (auto sit = indices.begin(); sit != indices.end(); ++sit) + { + if (turns[*sit].priority >= max_prio) + { + ok_indices.insert(*sit); + } + } + if (ok_indices.size() > 0 && ok_indices.size() < indices.size()) + { + indices = ok_indices; + std::cout << "%"; + } +} + +template <typename AngleInfo, typename Angles, typename Turns, typename TurnSegmentIndices> +inline void calculate_left_turns(Angles const& angles, + Turns& turns, TurnSegmentIndices const& turn_segment_indices, + std::set<int>& keep_indices) +{ + bool debug_indicate_size = false; + + typedef typename strategy::side::services::default_strategy + < + typename cs_tag<typename AngleInfo::point_type>::type + >::type side_strategy; + + std::size_t i = 0; + std::size_t n = boost::size(angles); + + typedef geometry::ever_circling_range_iterator<Angles const> circling_iterator; + circling_iterator cit(angles); + + debug_left_turn(*cit, *cit); + for(circling_iterator prev = cit++; i < n; prev = cit++, i++) + { + debug_left_turn(*cit, *prev); + + bool const include = ! geometry::math::equals(prev->angle, cit->angle) + && ! prev->incoming + && cit->incoming; + + if (include) + { + // Go back to possibly include more same left outgoing angles: + // Because we check on side too we can take a large "epsilon" + circling_iterator back = prev - 1; + + typename AngleInfo::angle_type eps = 0.00001; + int b = 1; + for(std::size_t d = 0; + math::abs(prev->angle - back->angle) < eps + && ! back->incoming + && d < n; + d++) + { + --back; + ++b; + } + + // Same but forward to possibly include more incoming angles + int f = 1; + circling_iterator forward = cit + 1; + for(std::size_t d = 0; + math::abs(cit->angle - forward->angle) < eps + && forward->incoming + && d < n; + d++) + { + ++forward; + ++f; + } + +#ifdef BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION + std::cout << "HANDLE " << b << "/" << f << " ANGLES" << std::endl; +#endif + for(circling_iterator bit = prev; bit != back; --bit) + { + int code = side_strategy::apply(bit->direction_point, prev->intersection_point, prev->direction_point); + int code2 = side_strategy::apply(prev->direction_point, bit->intersection_point, bit->direction_point); + for(circling_iterator fit = cit; fit != forward; ++fit) + { + int code3 = side_strategy::apply(fit->direction_point, cit->intersection_point, cit->direction_point); + int code4 = side_strategy::apply(cit->direction_point, fit->intersection_point, fit->direction_point); + + int priority = 2; + if (code == -1 && code2 == 1) + { + // This segment is lying right of the other one. + // Cannot ignore it (because of robustness, see a.o. #rt_p21 from unit test). + // But if we find more we can prefer later by priority + // (a.o. necessary for #rt_o2 from unit test) + priority = 1; + } + + bool included = include_left_turn_of_all(*bit, *fit, turns, turn_segment_indices, keep_indices, priority); + debug_left_turn(included ? "KEEP" : "SKIP", *fit, *bit, code, code2, code3, code4); + } + } + } + } + + if (debug_indicate_size) + { + std::cout << " size=" << keep_indices.size(); + } + + if (keep_indices.size() >= 2) + { + prefer_by_other(turns, turn_segment_indices, keep_indices); + } + if (keep_indices.size() >= 2) + { + prefer_by_priority(turns, keep_indices); + } +} + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_GET_LEFT_TURNS_HPP diff --git a/boost/geometry/algorithms/detail/occupation_info.hpp b/boost/geometry/algorithms/detail/occupation_info.hpp new file mode 100644 index 0000000000..e147ba12d8 --- /dev/null +++ b/boost/geometry/algorithms/detail/occupation_info.hpp @@ -0,0 +1,329 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2012 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_OCCUPATION_INFO_HPP +#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OCCUPATION_INFO_HPP + +#if ! defined(NDEBUG) + #define BOOST_GEOMETRY_DEBUG_BUFFER_OCCUPATION +#endif + +#include <algorithm> +#include <boost/range.hpp> + +#include <boost/geometry/core/coordinate_type.hpp> +#include <boost/geometry/core/point_type.hpp> + +#include <boost/geometry/algorithms/equals.hpp> +#include <boost/geometry/iterators/closing_iterator.hpp> + +#include <boost/geometry/algorithms/detail/get_left_turns.hpp> + + +namespace boost { namespace geometry +{ + + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +template <typename P> +class relaxed_less +{ + typedef typename geometry::coordinate_type<P>::type coordinate_type; + + coordinate_type epsilon; + +public : + + inline relaxed_less() + { + // TODO: adapt for ttmath, and maybe build the map in another way + // (e.g. exact constellations of segment-id's), maybe adaptive. + epsilon = std::numeric_limits<double>::epsilon() * 100.0; + } + + inline bool operator()(P const& a, P const& b) const + { + coordinate_type const dx = math::abs(geometry::get<0>(a) - geometry::get<0>(b)); + coordinate_type const dy = math::abs(geometry::get<1>(a) - geometry::get<1>(b)); + + + if (dx < epsilon && dy < epsilon) + { + return false; + } + if (dx < epsilon) + { + return geometry::get<1>(a) < geometry::get<1>(b); + } + + return geometry::get<0>(a) < geometry::get<0>(b); + } + + inline bool equals(P const& a, P const& b) const + { + typedef typename geometry::coordinate_type<P>::type coordinate_type; + + coordinate_type const dx = math::abs(geometry::get<0>(a) - geometry::get<0>(b)); + coordinate_type const dy = math::abs(geometry::get<1>(a) - geometry::get<1>(b)); + + return dx < epsilon && dy < epsilon; + }; +}; + + +template <typename T, typename P1, typename P2> +inline T calculate_angle(P1 const& from_point, P2 const& to_point) +{ + typedef P1 vector_type; + vector_type v = from_point; + geometry::subtract_point(v, to_point); + return atan2(geometry::get<1>(v), geometry::get<0>(v)); +} + +template <typename Iterator, typename Vector> +inline Iterator advance_circular(Iterator it, Vector const& vector, segment_identifier& seg_id, bool forward = true) +{ + int const increment = forward ? 1 : -1; + if (it == boost::begin(vector) && increment < 0) + { + it = boost::end(vector); + seg_id.segment_index = boost::size(vector); + } + it += increment; + seg_id.segment_index += increment; + if (it == boost::end(vector)) + { + seg_id.segment_index = 0; + it = boost::begin(vector); + } + return it; +} + +template <typename Point, typename T> +struct angle_info +{ + typedef T angle_type; + typedef Point point_type; + + segment_identifier seg_id; + int turn_index; + int operation_index; + Point intersection_point; + Point direction_point; + T angle; + bool incoming; +}; + +template <typename AngleInfo> +class occupation_info +{ + typedef std::vector<AngleInfo> collection_type; + + struct angle_sort + { + inline bool operator()(AngleInfo const& left, AngleInfo const& right) const + { + // In this case we can compare even double using equals + // return geometry::math::equals(left.angle, right.angle) + return left.angle == right.angle + ? int(left.incoming) < int(right.incoming) + : left.angle < right.angle + ; + } + }; + +public : + collection_type angles; +private : + bool m_occupied; + bool m_calculated; + + inline bool is_occupied() + { + if (boost::size(angles) <= 1) + { + return false; + } + + std::sort(angles.begin(), angles.end(), angle_sort()); + + typedef geometry::closing_iterator<collection_type const> closing_iterator; + closing_iterator vit(angles); + closing_iterator end(angles, true); + + closing_iterator prev = vit++; + for( ; vit != end; prev = vit++) + { + if (! geometry::math::equals(prev->angle, vit->angle) + && ! prev->incoming + && vit->incoming) + { + return false; + } + } + return true; + } + +public : + inline occupation_info() + : m_occupied(false) + , m_calculated(false) + {} + + template <typename PointC, typename Point1, typename Point2> + inline void add(PointC const& map_point, Point1 const& direction_point, Point2 const& intersection_point, + int turn_index, int operation_index, + segment_identifier const& seg_id, bool incoming) + { + //std::cout << "-> adding angle " << geometry::wkt(direction_point) << " .. " << geometry::wkt(intersection_point) << " " << int(incoming) << std::endl; + if (geometry::equals(direction_point, intersection_point)) + { + //std::cout << "EQUAL! Skipping" << std::endl; + return; + } + + AngleInfo info; + info.incoming = incoming; + info.angle = calculate_angle<typename AngleInfo::angle_type>(direction_point, map_point); + info.seg_id = seg_id; + info.turn_index = turn_index; + info.operation_index = operation_index; + info.intersection_point = intersection_point; + info.direction_point = direction_point; + angles.push_back(info); + + m_calculated = false; + } + + inline bool occupied() + { + if (! m_calculated) + { + m_occupied = is_occupied(); + m_calculated = true; + } + return m_occupied; + } + + template <typename Turns, typename TurnSegmentIndices> + inline void get_left_turns( + Turns& turns, TurnSegmentIndices const& turn_segment_indices, + std::set<int>& keep_indices) + { + std::sort(angles.begin(), angles.end(), angle_sort()); + calculate_left_turns<AngleInfo>(angles, turns, turn_segment_indices, keep_indices); + } +}; + + +template <typename Point, typename Ring, typename Info> +inline void add_incoming_and_outgoing_angles(Point const& map_point, Point const& intersection_point, + Ring const& ring, + int turn_index, + int operation_index, + segment_identifier seg_id, + Info& info) +{ + typedef typename boost::range_iterator + < + Ring const + >::type iterator_type; + + int const n = boost::size(ring); + if (seg_id.segment_index >= n || seg_id.segment_index < 0) + { + return; + } + + segment_identifier real_seg_id = seg_id; + iterator_type it = boost::begin(ring) + seg_id.segment_index; + + // TODO: if we use turn-info ("to", "middle"), we know if to advance without resorting to equals + relaxed_less<Point> comparator; + + if (comparator.equals(intersection_point, *it)) + { + // It should be equal only once. But otherwise we skip it (in "add") + it = advance_circular(it, ring, seg_id, false); + } + + info.add(map_point, *it, intersection_point, turn_index, operation_index, real_seg_id, true); + + if (comparator.equals(intersection_point, *it)) + { + it = advance_circular(it, ring, real_seg_id); + } + else + { + // Don't upgrade the ID + it = advance_circular(it, ring, seg_id); + } + for (int defensive_check = 0; + comparator.equals(intersection_point, *it) && defensive_check < n; + defensive_check++) + { + it = advance_circular(it, ring, real_seg_id); + } + + info.add(map_point, *it, intersection_point, turn_index, operation_index, real_seg_id, false); +} + + +// Map in two senses of the word: it is a std::map where the key is a point. +// Per point an "occupation_info" record is kept +// Used for the buffer (but will also be used for intersections/unions having complex self-tangencies) +template <typename Point, typename OccupationInfo> +class occupation_map +{ +public : + typedef std::map<Point, OccupationInfo, relaxed_less<Point> > map_type; + + map_type map; + std::set<int> turn_indices; + + inline OccupationInfo& find_or_insert(Point const& point, Point& mapped_point) + { + typename map_type::iterator it = map.find(point); + if (it == boost::end(map)) + { + std::pair<typename map_type::iterator, bool> pair + = map.insert(std::make_pair(point, OccupationInfo())); + it = pair.first; + } + mapped_point = it->first; + return it->second; + } + + inline bool contains(Point const& point) const + { + typename map_type::const_iterator it = map.find(point); + return it != boost::end(map); + } + + inline bool contains_turn_index(int index) const + { + return turn_indices.count(index) > 0; + } + + inline void insert_turn_index(int index) + { + turn_indices.insert(index); + } +}; + + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OCCUPATION_INFO_HPP diff --git a/boost/geometry/algorithms/detail/overlay/add_rings.hpp b/boost/geometry/algorithms/detail/overlay/add_rings.hpp index eb3e60e483..74595fedd0 100644 --- a/boost/geometry/algorithms/detail/overlay/add_rings.hpp +++ b/boost/geometry/algorithms/detail/overlay/add_rings.hpp @@ -9,6 +9,8 @@ #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_ADD_RINGS_HPP #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_ADD_RINGS_HPP +#include <boost/geometry/core/closure.hpp> +#include <boost/geometry/algorithms/area.hpp> #include <boost/geometry/algorithms/detail/overlay/convert_ring.hpp> #include <boost/geometry/algorithms/detail/overlay/get_ring.hpp> @@ -73,6 +75,21 @@ inline OutputIterator add_rings(SelectionMap const& map, OutputIterator out) { typedef typename SelectionMap::const_iterator iterator; + typedef typename SelectionMap::mapped_type property_type; + typedef typename property_type::area_type area_type; + + area_type const zero = 0; + std::size_t const min_num_points = core_detail::closure::minimum_ring_size + < + geometry::closure + < + typename boost::range_value + < + RingCollection const + >::type + >::value + >::value; + for (iterator it = boost::begin(map); it != boost::end(map); @@ -99,7 +116,16 @@ inline OutputIterator add_rings(SelectionMap const& map, *child_it, mit->second.reversed, true); } } - *out++ = result; + + // Only add rings if they satisfy minimal requirements. + // This cannot be done earlier (during traversal), not + // everything is figured out yet (sum of positive/negative rings) + // TODO: individual rings can still contain less than 3 points. + if (geometry::num_points(result) >= min_num_points + && math::larger(geometry::area(result), zero)) + { + *out++ = result; + } } } return out; diff --git a/boost/geometry/algorithms/detail/overlay/assign_parents.hpp b/boost/geometry/algorithms/detail/overlay/assign_parents.hpp index 133530563e..5063f49eb4 100644 --- a/boost/geometry/algorithms/detail/overlay/assign_parents.hpp +++ b/boost/geometry/algorithms/detail/overlay/assign_parents.hpp @@ -130,7 +130,7 @@ struct assign_visitor return; } - if (outer.real_area > 0) + if (math::larger(outer.real_area, 0)) { if (inner.real_area < 0 || m_check_for_orientation) { @@ -317,13 +317,14 @@ template > inline void assign_parents(Geometry const& geometry, RingCollection const& collection, - RingMap& ring_map) + RingMap& ring_map, + bool check_for_orientation) { // Call it with an empty geometry // (ring_map should be empty for source_id==1) Geometry empty; - assign_parents(geometry, empty, collection, ring_map, true); + assign_parents(geometry, empty, collection, ring_map, check_for_orientation); } diff --git a/boost/geometry/algorithms/detail/overlay/calculate_distance_policy.hpp b/boost/geometry/algorithms/detail/overlay/calculate_distance_policy.hpp index 3e6a8897f5..2003d2350d 100644 --- a/boost/geometry/algorithms/detail/overlay/calculate_distance_policy.hpp +++ b/boost/geometry/algorithms/detail/overlay/calculate_distance_policy.hpp @@ -32,9 +32,18 @@ struct calculate_distance_policy { static bool const include_no_turn = false; static bool const include_degenerate = false; - - template <typename Point1, typename Point2, typename Info> - static inline void apply(Info& info, Point1 const& p1, Point2 const& p2) + static bool const include_opposite = false; + + template + < + typename Info, + typename Point1, + typename Point2, + typename IntersectionInfo, + typename DirInfo + > + static inline void apply(Info& info, Point1 const& p1, Point2 const& p2, + IntersectionInfo const&, DirInfo const&) { info.operations[0].enriched.distance = geometry::comparable_distance(info.point, p1); diff --git a/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp b/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp index 379444428f..5e18d0453a 100644 --- a/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp +++ b/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp @@ -55,14 +55,14 @@ struct copy_segment_point_range if (second) { index++; - if (index >= boost::size(range)) + if (index >= int(boost::size(range))) { index = 0; } } // Exception? - if (index >= boost::size(range)) + if (index >= int(boost::size(range))) { return false; } diff --git a/boost/geometry/algorithms/detail/overlay/copy_segments.hpp b/boost/geometry/algorithms/detail/overlay/copy_segments.hpp index b0183a3ac2..805f3923e3 100644 --- a/boost/geometry/algorithms/detail/overlay/copy_segments.hpp +++ b/boost/geometry/algorithms/detail/overlay/copy_segments.hpp @@ -78,7 +78,7 @@ struct copy_segments_ring int const from_index = seg_id.segment_index + 1; // Sanity check - BOOST_ASSERT(from_index < boost::size(view)); + BOOST_ASSERT(from_index < int(boost::size(view))); ec_iterator it(boost::begin(view), boost::end(view), boost::begin(view) + from_index); @@ -89,7 +89,7 @@ struct copy_segments_ring typedef typename boost::range_difference<Ring>::type size_type; size_type const count = from_index <= to_index ? to_index - from_index + 1 - : boost::size(view) - from_index + to_index + 1; + : int(boost::size(view)) - from_index + to_index + 1; for (size_type i = 0; i < count; ++i, ++it) { @@ -117,7 +117,7 @@ struct copy_segments_linestring int const from_index = seg_id.segment_index + 1; // Sanity check - if (from_index > to_index || from_index < 0 || to_index >= boost::size(ls)) + if (from_index > to_index || from_index < 0 || to_index >= int(boost::size(ls))) { return; } diff --git a/boost/geometry/algorithms/detail/overlay/debug_turn_info.hpp b/boost/geometry/algorithms/detail/overlay/debug_turn_info.hpp index 0bec816946..0cc34255ca 100644 --- a/boost/geometry/algorithms/detail/overlay/debug_turn_info.hpp +++ b/boost/geometry/algorithms/detail/overlay/debug_turn_info.hpp @@ -28,6 +28,7 @@ inline char method_char(detail::overlay::method_type const& method) case method_touch_interior : return 'm'; case method_collinear : return 'c'; case method_equal : return 'e'; + case method_error : return '!'; default : return '?'; } } @@ -42,6 +43,7 @@ inline char operation_char(detail::overlay::operation_type const& operation) case operation_intersection : return 'i'; case operation_blocked : return 'x'; case operation_continue : return 'c'; + case operation_opposite : return 'o'; default : return '?'; } } diff --git a/boost/geometry/algorithms/detail/overlay/follow.hpp b/boost/geometry/algorithms/detail/overlay/follow.hpp index 087092a5f2..b110cc9602 100644 --- a/boost/geometry/algorithms/detail/overlay/follow.hpp +++ b/boost/geometry/algorithms/detail/overlay/follow.hpp @@ -279,12 +279,12 @@ class follow { // In case of turn point at the same location, we want to have continue/blocked LAST // because that should be followed (intersection) or skipped (difference). - // By chance the enumeration is ordered like that but we keep the safe way here. inline int operation_order(Turn const& turn) const { operation_type const& operation = turn.operations[0].operation; switch(operation) { + case operation_opposite : return 0; case operation_none : return 0; case operation_union : return 1; case operation_intersection : return 2; diff --git a/boost/geometry/algorithms/detail/overlay/get_ring.hpp b/boost/geometry/algorithms/detail/overlay/get_ring.hpp index ad47665500..c2c6980577 100644 --- a/boost/geometry/algorithms/detail/overlay/get_ring.hpp +++ b/boost/geometry/algorithms/detail/overlay/get_ring.hpp @@ -83,7 +83,7 @@ struct get_ring<polygon_tag> BOOST_ASSERT ( id.ring_index >= -1 - && id.ring_index < boost::size(interior_rings(polygon)) + && id.ring_index < int(boost::size(interior_rings(polygon))) ); return id.ring_index < 0 ? exterior_ring(polygon) diff --git a/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp b/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp index 663d70d9af..b8320d9b7b 100644 --- a/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp +++ b/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp @@ -250,9 +250,15 @@ struct touch : public base_turn_handler int const side_pk_q2 = SideStrategy::apply(qj, qk, pk); int const side_pk_p = SideStrategy::apply(pi, pj, pk); int const side_qk_q = SideStrategy::apply(qi, qj, qk); + + bool const both_continue = side_pk_p == 0 && side_qk_q == 0; + bool const robustness_issue_in_continue = both_continue && side_pk_q2 != 0; + bool const q_turns_left = side_qk_q == 1; bool const block_q = side_qk_p1 == 0 - && ! same(side_qi_p1, side_qk_q); + && ! same(side_qi_p1, side_qk_q) + && ! robustness_issue_in_continue + ; // If Pk at same side as Qi/Qk // (the "or" is for collinear case) @@ -278,7 +284,7 @@ struct touch : public base_turn_handler if (side_pk_q1 == 0) { ti.operations[0].operation = operation_blocked; - // Q turns right -> union (both independant), + // Q turns right -> union (both independent), // Q turns left -> intersection ti.operations[1].operation = block_q ? operation_blocked : q_turns_left ? operation_intersection @@ -466,6 +472,45 @@ struct equal : public base_turn_handler template < typename TurnInfo, + typename AssignPolicy +> +struct equal_opposite : public base_turn_handler +{ + template + < + typename Point1, + typename Point2, + typename OutputIterator, + typename IntersectionInfo, + typename DirInfo + > + static inline void apply(Point1 const& pi, Point2 const& qi, + /* by value: */ TurnInfo tp, + OutputIterator& out, + IntersectionInfo const& intersection_info, + DirInfo const& dir_info) + { + // For equal-opposite segments, normally don't do anything. + if (AssignPolicy::include_opposite) + { + tp.method = method_equal; + for (int i = 0; i < 2; i++) + { + tp.operations[i].operation = operation_opposite; + } + for (unsigned int i = 0; i < intersection_info.count; i++) + { + geometry::convert(intersection_info.intersections[i], tp.point); + AssignPolicy::apply(tp, pi, qi, intersection_info, dir_info); + *out++ = tp; + } + } + } +}; + +template +< + typename TurnInfo, typename SideStrategy > struct collinear : public base_turn_handler @@ -494,6 +539,13 @@ struct collinear : public base_turn_handler - if P arrives and P turns right: intersection for P - if Q arrives and Q turns left: union for Q (=intersection for P) - if Q arrives and Q turns right: intersection for Q (=union for P) + + ROBUSTNESS: p and q are collinear, so you would expect + that side qk//p1 == pk//q1. But that is not always the case + in near-epsilon ranges. Then decision logic is different. + If p arrives, q is further, so the angle qk//p1 is (normally) + more precise than pk//p1 + */ template < @@ -516,12 +568,18 @@ struct collinear : public base_turn_handler // Should not be 0, this is checked before BOOST_ASSERT(arrival != 0); + int const side_p = SideStrategy::apply(pi, pj, pk); + int const side_q = SideStrategy::apply(qi, qj, qk); + // If p arrives, use p, else use q int const side_p_or_q = arrival == 1 - ? SideStrategy::apply(pi, pj, pk) - : SideStrategy::apply(qi, qj, qk) + ? side_p + : side_q ; + int const side_pk = SideStrategy::apply(qi, qj, pk); + int const side_qk = SideStrategy::apply(pi, pj, qk); + // See comments above, // resulting in a strange sort of mathematic rule here: // The arrival-info multiplied by the relevant side @@ -529,7 +587,15 @@ struct collinear : public base_turn_handler int const product = arrival * side_p_or_q; - if(product == 0) + // Robustness: side_p is supposed to be equal to side_pk (because p/q are collinear) + // and side_q to side_qk + bool const robustness_issue = side_pk != side_p || side_qk != side_q; + + if (robustness_issue) + { + handle_robustness(ti, arrival, side_p, side_q, side_pk, side_qk); + } + else if(product == 0) { both(ti, operation_continue); } @@ -538,6 +604,38 @@ struct collinear : public base_turn_handler ui_else_iu(product == 1, ti); } } + + static inline void handle_robustness(TurnInfo& ti, int arrival, + int side_p, int side_q, int side_pk, int side_qk) + { + // We take the longer one, i.e. if q arrives in p (arrival == -1), + // then p exceeds q and we should take p for a union... + + bool use_p_for_union = arrival == -1; + + // ... unless one of the sides consistently directs to the other side + int const consistent_side_p = side_p == side_pk ? side_p : 0; + int const consistent_side_q = side_q == side_qk ? side_q : 0; + if (arrival == -1 && (consistent_side_p == -1 || consistent_side_q == 1)) + { + use_p_for_union = false; + } + if (arrival == 1 && (consistent_side_p == 1 || consistent_side_q == -1)) + { + use_p_for_union = true; + } + + //std::cout << "ROBUSTNESS -> Collinear " + // << " arr: " << arrival + // << " dir: " << side_p << " " << side_q + // << " rev: " << side_pk << " " << side_qk + // << " cst: " << cside_p << " " << cside_q + // << std::boolalpha << " " << use_p_for_union + // << std::endl; + + ui_else_iu(use_p_for_union, ti); + } + }; template @@ -583,6 +681,7 @@ private : TurnInfo& tp, IntersectionInfo const& intersection_info) { int const side_rk_r = SideStrategy::apply(ri, rj, rk); + operation_type blocked = operation_blocked; switch(side_rk_r) { @@ -596,15 +695,24 @@ private : break; case 0 : // No turn on opposite collinear: block, do not traverse - // But this "xx" is ignored here, it is useless to include - // two operation blocked, so the whole point does not need + // But this "xx" is usually ignored, it is useless to include + // two operations blocked, so the whole point does not need // to be generated. // So return false to indicate nothing is to be done. - return false; + if (AssignPolicy::include_opposite) + { + tp.operations[Index].operation = operation_opposite; + blocked = operation_opposite; + } + else + { + return false; + } + break; } // The other direction is always blocked when collinear opposite - tp.operations[1 - Index].operation = operation_blocked; + tp.operations[1 - Index].operation = blocked; // If P arrives within Q, set info on P (which is done above, index=0), // this turn-info belongs to the second intersection point, index=1 @@ -633,32 +741,45 @@ public: IntersectionInfo const& intersection_info, DirInfo const& dir_info) { - /* - std::cout << "arrivals: " - << dir_info.arrival[0] - << "/" << dir_info.arrival[1] - << std::endl; - */ - TurnInfo tp = tp_model; tp.method = method_collinear; - // If P arrives within Q, there is a turn dependant on P + // If P arrives within Q, there is a turn dependent on P if (dir_info.arrival[0] == 1 && set_tp<0>(pi, pj, pk, tp, intersection_info)) { - AssignPolicy::apply(tp, pi, qi); + AssignPolicy::apply(tp, pi, qi, intersection_info, dir_info); *out++ = tp; } - // If Q arrives within P, there is a turn dependant on Q + // If Q arrives within P, there is a turn dependent on Q if (dir_info.arrival[1] == 1 && set_tp<1>(qi, qj, qk, tp, intersection_info)) { - AssignPolicy::apply(tp, pi, qi); + AssignPolicy::apply(tp, pi, qi, intersection_info, dir_info); *out++ = tp; } + + if (AssignPolicy::include_opposite) + { + // Handle cases not yet handled above + if ((dir_info.arrival[1] == -1 && dir_info.arrival[0] == 0) + || (dir_info.arrival[0] == -1 && dir_info.arrival[1] == 0)) + { + for (int i = 0; i < 2; i++) + { + tp.operations[i].operation = operation_opposite; + } + for (unsigned int i = 0; i < intersection_info.count; i++) + { + geometry::convert(intersection_info.intersections[i], tp.point); + AssignPolicy::apply(tp, pi, qi, intersection_info, dir_info); + *out++ = tp; + } + } + } + } }; @@ -722,9 +843,17 @@ struct assign_null_policy { static bool const include_no_turn = false; static bool const include_degenerate = false; - - template <typename Point1, typename Point2, typename Info> - static inline void apply(Info& , Point1 const& , Point2 const& ) + static bool const include_opposite = false; + + template + < + typename Info, + typename Point1, + typename Point2, + typename IntersectionInfo, + typename DirInfo + > + static inline void apply(Info& , Point1 const& , Point2 const&, IntersectionInfo const&, DirInfo const&) {} }; @@ -763,8 +892,9 @@ struct get_turn_info typedef typename si::segment_intersection_strategy_type strategy; - - + // Intersect pi-pj with qi-qj + // The points pk and qk are only used do determine more information + // about the turn. template <typename OutputIterator> static inline OutputIterator apply( Point1 const& pi, Point1 const& pj, Point1 const& pk, @@ -796,7 +926,7 @@ struct get_turn_info { only_convert<TurnInfo>::apply(tp, result.template get<0>()); - AssignPolicy::apply(tp, pi, qi); + AssignPolicy::apply(tp, pi, qi, result.template get<0>(), result.template get<1>()); *out++ = tp; } break; @@ -824,7 +954,7 @@ struct get_turn_info policy::template apply<1>(qi, qj, qk, pi, pj, pk, tp, result.template get<0>(), result.template get<1>()); } - AssignPolicy::apply(tp, pi, qi); + AssignPolicy::apply(tp, pi, qi, result.template get<0>(), result.template get<1>()); *out++ = tp; } break; @@ -838,7 +968,7 @@ struct get_turn_info policy::apply(pi, pj, pk, qi, qj, qk, tp, result.template get<0>(), result.template get<1>()); - AssignPolicy::apply(tp, pi, qi); + AssignPolicy::apply(tp, pi, qi, result.template get<0>(), result.template get<1>()); *out++ = tp; } break; @@ -853,7 +983,7 @@ struct get_turn_info policy::apply(pi, pj, pk, qi, qj, qk, tp, result.template get<0>(), result.template get<1>()); - AssignPolicy::apply(tp, pi, qi); + AssignPolicy::apply(tp, pi, qi, result.template get<0>(), result.template get<1>()); *out++ = tp; } break; @@ -871,10 +1001,18 @@ struct get_turn_info policy::apply(pi, pj, pk, qi, qj, qk, tp, result.template get<0>(), result.template get<1>()); - AssignPolicy::apply(tp, pi, qi); + AssignPolicy::apply(tp, pi, qi, result.template get<0>(), result.template get<1>()); *out++ = tp; } - // If they ARE opposite, don't do anything. + else + { + equal_opposite + < + TurnInfo, + AssignPolicy + >::apply(pi, qi, + tp, out, result.template get<0>(), result.template get<1>()); + } } break; case 'c' : @@ -906,7 +1044,7 @@ struct get_turn_info tp, result.template get<0>(), result.template get<1>()); } - AssignPolicy::apply(tp, pi, qi); + AssignPolicy::apply(tp, pi, qi, result.template get<0>(), result.template get<1>()); *out++ = tp; } else @@ -927,7 +1065,7 @@ struct get_turn_info if (AssignPolicy::include_degenerate) { only_convert<TurnInfo>::apply(tp, result.template get<0>()); - AssignPolicy::apply(tp, pi, qi); + AssignPolicy::apply(tp, pi, qi, result.template get<0>(), result.template get<1>()); *out++ = tp; } } diff --git a/boost/geometry/algorithms/detail/overlay/get_turns.hpp b/boost/geometry/algorithms/detail/overlay/get_turns.hpp index 5740a8b7ba..26629043cb 100644 --- a/boost/geometry/algorithms/detail/overlay/get_turns.hpp +++ b/boost/geometry/algorithms/detail/overlay/get_turns.hpp @@ -495,7 +495,7 @@ struct get_turns_cs int source_id1, Range const& range, int source_id2, Box const& box, Turns& turns, - InterruptPolicy& , + InterruptPolicy& interrupt_policy, int multi_index = -1, int ring_index = -1) { if (boost::size(range) <= 1) @@ -557,7 +557,7 @@ struct get_turns_cs get_turns_with_box(seg_id, source_id2, *prev, *it, *next, bp[0], bp[1], bp[2], bp[3], - turns); + turns, interrupt_policy); // Future performance enhancement: // return if told by the interrupt policy } @@ -596,7 +596,8 @@ private: box_point_type const& bp2, box_point_type const& bp3, // Output - Turns& turns) + Turns& turns, + InterruptPolicy& interrupt_policy) { // Depending on code some relations can be left out @@ -622,6 +623,12 @@ private: ti.operations[1].seg_id = segment_identifier(source_id2, -1, -1, 3); TurnPolicy::apply(rp0, rp1, rp2, bp3, bp0, bp1, ti, std::back_inserter(turns)); + + if (InterruptPolicy::enabled) + { + interrupt_policy.apply(turns); + } + } }; diff --git a/boost/geometry/algorithms/detail/overlay/overlay.hpp b/boost/geometry/algorithms/detail/overlay/overlay.hpp index ab5b6d123d..41665e0af1 100644 --- a/boost/geometry/algorithms/detail/overlay/overlay.hpp +++ b/boost/geometry/algorithms/detail/overlay/overlay.hpp @@ -233,7 +233,7 @@ std::cout << "traverse" << std::endl; std::cout << "map_turns: " << timer.elapsed() << std::endl; #endif - typedef ring_properties<typename geometry::point_type<Geometry1>::type> properties; + typedef ring_properties<typename geometry::point_type<GeometryOut>::type> properties; std::map<ring_identifier, properties> selected; select_rings<Direction>(geometry1, geometry2, map, selected, ! turn_points.empty()); diff --git a/boost/geometry/algorithms/detail/overlay/turn_info.hpp b/boost/geometry/algorithms/detail/overlay/turn_info.hpp index aa6b428f19..89a60b21ab 100644 --- a/boost/geometry/algorithms/detail/overlay/turn_info.hpp +++ b/boost/geometry/algorithms/detail/overlay/turn_info.hpp @@ -28,7 +28,8 @@ enum operation_type operation_union, operation_intersection, operation_blocked, - operation_continue + operation_continue, + operation_opposite }; @@ -102,6 +103,12 @@ struct turn_info { return has12(type, type); } + + inline bool has(operation_type type) const + { + return this->operations[0].operation == type + || this->operations[1].operation == type; + } inline bool combination(operation_type type1, operation_type type2) const { @@ -114,10 +121,13 @@ struct turn_info { return both(operation_blocked); } + inline bool opposite() const + { + return both(operation_opposite); + } inline bool any_blocked() const { - return this->operations[0].operation == operation_blocked - || this->operations[1].operation == operation_blocked; + return has(operation_blocked); } diff --git a/boost/geometry/algorithms/detail/partition.hpp b/boost/geometry/algorithms/detail/partition.hpp index 7a7de2cdd3..45ff52ccb1 100644 --- a/boost/geometry/algorithms/detail/partition.hpp +++ b/boost/geometry/algorithms/detail/partition.hpp @@ -153,12 +153,12 @@ class partition_one_collection static inline void next_level(Box const& box, InputCollection const& collection, index_vector_type const& input, - int level, int min_elements, + int level, std::size_t min_elements, Policy& policy, VisitBoxPolicy& box_policy) { if (boost::size(input) > 0) { - if (boost::size(input) > min_elements && level < 100) + if (std::size_t(boost::size(input)) > min_elements && level < 100) { sub_divide::apply(box, collection, input, level + 1, min_elements, policy, box_policy); @@ -176,7 +176,7 @@ public : InputCollection const& collection, index_vector_type const& input, int level, - int min_elements, + std::size_t min_elements, Policy& policy, VisitBoxPolicy& box_policy) { box_policy.apply(box, level); @@ -230,13 +230,13 @@ class partition_two_collections index_vector_type const& input1, InputCollection const& collection2, index_vector_type const& input2, - int level, int min_elements, + int level, std::size_t min_elements, Policy& policy, VisitBoxPolicy& box_policy) { if (boost::size(input1) > 0 && boost::size(input2) > 0) { - if (boost::size(input1) > min_elements - && boost::size(input2) > min_elements + if (std::size_t(boost::size(input1)) > min_elements + && std::size_t(boost::size(input2)) > min_elements && level < 100) { sub_divide::apply(box, collection1, input1, collection2, @@ -257,7 +257,7 @@ public : InputCollection const& collection1, index_vector_type const& input1, InputCollection const& collection2, index_vector_type const& input2, int level, - int min_elements, + std::size_t min_elements, Policy& policy, VisitBoxPolicy& box_policy) { box_policy.apply(box, level); @@ -335,11 +335,11 @@ public : template <typename InputCollection, typename VisitPolicy> static inline void apply(InputCollection const& collection, VisitPolicy& visitor, - int min_elements = 16, + std::size_t min_elements = 16, VisitBoxPolicy box_visitor = visit_no_policy() ) { - if (boost::size(collection) > min_elements) + if (std::size_t(boost::size(collection)) > min_elements) { index_vector_type index_vector; Box total; @@ -377,12 +377,12 @@ public : static inline void apply(InputCollection const& collection1, InputCollection const& collection2, VisitPolicy& visitor, - int min_elements = 16, + std::size_t min_elements = 16, VisitBoxPolicy box_visitor = visit_no_policy() ) { - if (boost::size(collection1) > min_elements - && boost::size(collection2) > min_elements) + if (std::size_t(boost::size(collection1)) > min_elements + && std::size_t(boost::size(collection2)) > min_elements) { index_vector_type index_vector1, index_vector2; Box total; diff --git a/boost/geometry/algorithms/detail/point_on_border.hpp b/boost/geometry/algorithms/detail/point_on_border.hpp index 33177924aa..b7e15ba3f9 100644 --- a/boost/geometry/algorithms/detail/point_on_border.hpp +++ b/boost/geometry/algorithms/detail/point_on_border.hpp @@ -25,6 +25,7 @@ #include <boost/geometry/geometries/concepts/check.hpp> #include <boost/geometry/algorithms/assign.hpp> +#include <boost/geometry/algorithms/detail/convert_point_to_point.hpp> #include <boost/geometry/algorithms/detail/disjoint.hpp> @@ -50,7 +51,8 @@ struct get_point template<typename Point, std::size_t Dimension, std::size_t DimensionCount> struct midpoint_helper { - static inline bool apply(Point& p, Point const& p1, Point const& p2) + template <typename InputPoint> + static inline bool apply(Point& p, InputPoint const& p1, InputPoint const& p2) { typename coordinate_type<Point>::type const two = 2; set<Dimension>(p, @@ -63,7 +65,8 @@ struct midpoint_helper template <typename Point, std::size_t DimensionCount> struct midpoint_helper<Point, DimensionCount, DimensionCount> { - static inline bool apply(Point& , Point const& , Point const& ) + template <typename InputPoint> + static inline bool apply(Point& , InputPoint const& , InputPoint const& ) { return true; } @@ -102,7 +105,7 @@ struct point_on_range if (n > 0) { - point = *boost::begin(range); + geometry::detail::conversion::convert_point_to_point(*boost::begin(range), point); return true; } return false; diff --git a/boost/geometry/algorithms/detail/sections/sectionalize.hpp b/boost/geometry/algorithms/detail/sections/sectionalize.hpp index 36bcbdd6e7..a6e6837fe4 100644 --- a/boost/geometry/algorithms/detail/sections/sectionalize.hpp +++ b/boost/geometry/algorithms/detail/sections/sectionalize.hpp @@ -254,7 +254,7 @@ struct sectionalize_part Range const& range, ring_identifier ring_id) { - if (boost::size(range) <= index) + if (int(boost::size(range)) <= index) { return; } diff --git a/boost/geometry/algorithms/disjoint.hpp b/boost/geometry/algorithms/disjoint.hpp index f36b8dbdd0..f986cc24af 100644 --- a/boost/geometry/algorithms/disjoint.hpp +++ b/boost/geometry/algorithms/disjoint.hpp @@ -27,6 +27,7 @@ #include <boost/geometry/core/reverse_dispatch.hpp> #include <boost/geometry/algorithms/detail/disjoint.hpp> +#include <boost/geometry/algorithms/detail/for_each_range.hpp> #include <boost/geometry/algorithms/detail/point_on_border.hpp> #include <boost/geometry/algorithms/detail/overlay/get_turns.hpp> #include <boost/geometry/algorithms/within.hpp> @@ -44,15 +45,57 @@ namespace boost { namespace geometry namespace detail { namespace disjoint { +template<typename Geometry> +struct check_each_ring_for_within +{ + bool has_within; + Geometry const& m_geometry; + + inline check_each_ring_for_within(Geometry const& g) + : has_within(false) + , m_geometry(g) + {} + + template <typename Range> + inline void apply(Range const& range) + { + typename geometry::point_type<Range>::type p; + geometry::point_on_border(p, range); + if (geometry::within(p, m_geometry)) + { + has_within = true; + } + } +}; + +template <typename FirstGeometry, typename SecondGeometry> +inline bool rings_containing(FirstGeometry const& geometry1, + SecondGeometry const& geometry2) +{ + check_each_ring_for_within<FirstGeometry> checker(geometry1); + geometry::detail::for_each_range(geometry2, checker); + return checker.has_within; +} + + struct assign_disjoint_policy { // We want to include all points: static bool const include_no_turn = true; static bool const include_degenerate = true; + static bool const include_opposite = true; // We don't assign extra info: - template <typename Point1, typename Point2, typename Info> - static inline void apply(Info& , Point1 const& , Point2 const& ) + template + < + typename Info, + typename Point1, + typename Point2, + typename IntersectionInfo, + typename DirInfo + > + static inline void apply(Info& , Point1 const& , Point2 const&, + IntersectionInfo const&, DirInfo const&) {} }; @@ -107,8 +150,6 @@ struct disjoint_segment } }; - - template <typename Geometry1, typename Geometry2> struct general_areal { @@ -119,20 +160,10 @@ struct general_areal return false; } - typedef typename geometry::point_type<Geometry1>::type point_type; - // If there is no intersection of segments, they might located // inside each other - point_type p1; - geometry::point_on_border(p1, geometry1); - if (geometry::within(p1, geometry2)) - { - return false; - } - - typename geometry::point_type<Geometry1>::type p2; - geometry::point_on_border(p2, geometry2); - if (geometry::within(p2, geometry1)) + if (rings_containing(geometry1, geometry2) + || rings_containing(geometry2, geometry1)) { return false; } diff --git a/boost/geometry/algorithms/distance.hpp b/boost/geometry/algorithms/distance.hpp index eca3b03c7b..11c2bc929b 100644 --- a/boost/geometry/algorithms/distance.hpp +++ b/boost/geometry/algorithms/distance.hpp @@ -130,7 +130,7 @@ struct point_to_range // check if other segments are closer for (++prev, ++it; it != boost::end(view); ++prev, ++it) { - return_type const ds = ps_strategy.apply(point, *prev, *it); + return_type const ds = eps_strategy.apply(point, *prev, *it); if (geometry::math::equals(ds, zero)) { return ds; diff --git a/boost/geometry/algorithms/touches.hpp b/boost/geometry/algorithms/touches.hpp new file mode 100644 index 0000000000..7d424af428 --- /dev/null +++ b/boost/geometry/algorithms/touches.hpp @@ -0,0 +1,181 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2008-2012 Bruno Lalande, Paris, France. +// Copyright (c) 2009-2012 Mateusz Loskot, London, UK. + +// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library +// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_ALGORITHMS_TOUCHES_HPP +#define BOOST_GEOMETRY_ALGORITHMS_TOUCHES_HPP + + +#include <deque> + +#include <boost/geometry/geometries/concepts/check.hpp> +#include <boost/geometry/algorithms/detail/overlay/self_turn_points.hpp> +#include <boost/geometry/algorithms/detail/overlay/get_turns.hpp> +#include <boost/geometry/algorithms/disjoint.hpp> +#include <boost/geometry/algorithms/intersects.hpp> +#include <boost/geometry/algorithms/num_geometries.hpp> + + +namespace boost { namespace geometry +{ + +namespace detail { namespace touches +{ + +template <typename Turn> +inline bool ok_for_touch(Turn const& turn) +{ + return turn.both(detail::overlay::operation_union) + || turn.both(detail::overlay::operation_blocked) + || turn.combination(detail::overlay::operation_union, detail::overlay::operation_blocked) + ; +} + +template <typename Turns> +inline bool has_only_turns(Turns const& turns) +{ + bool has_touch = false; + typedef typename boost::range_iterator<Turns const>::type iterator_type; + for (iterator_type it = boost::begin(turns); it != boost::end(turns); ++it) + { + if (it->has(detail::overlay::operation_intersection)) + { + return false; + } + + switch(it->method) + { + case detail::overlay::method_crosses: + return false; + case detail::overlay::method_equal: + // Segment spatially equal means: at the right side + // the polygon internally overlaps. So return false. + return false; + case detail::overlay::method_touch: + case detail::overlay::method_touch_interior: + case detail::overlay::method_collinear: + if (ok_for_touch(*it)) + { + has_touch = true; + } + else + { + return false; + } + break; + case detail::overlay::method_none : + case detail::overlay::method_disjoint : + case detail::overlay::method_error : + break; + } + } + return has_touch; +} + +}} + +/*! +\brief \brief_check{has at least one touching point (self-tangency)} +\note This function can be called for one geometry (self-tangency) and + also for two geometries (touch) +\ingroup touches +\tparam Geometry \tparam_geometry +\param geometry \param_geometry +\return \return_check{is self-touching} + +\qbk{distinguish,one geometry} +\qbk{[def __one_parameter__]} +\qbk{[include reference/algorithms/touches.qbk]} +*/ +template <typename Geometry> +inline bool touches(Geometry const& geometry) +{ + concept::check<Geometry const>(); + + typedef detail::overlay::turn_info + < + typename geometry::point_type<Geometry>::type + > turn_info; + + typedef detail::overlay::get_turn_info + < + typename point_type<Geometry>::type, + typename point_type<Geometry>::type, + turn_info, + detail::overlay::assign_null_policy + > policy_type; + + std::deque<turn_info> turns; + detail::self_get_turn_points::no_interrupt_policy policy; + detail::self_get_turn_points::get_turns + < + Geometry, + std::deque<turn_info>, + policy_type, + detail::self_get_turn_points::no_interrupt_policy + >::apply(geometry, turns, policy); + + return detail::touches::has_only_turns(turns); +} + + +/*! +\brief \brief_check2{have at least one touching point (tangent - non overlapping)} +\ingroup touches +\tparam Geometry1 \tparam_geometry +\tparam Geometry2 \tparam_geometry +\param geometry1 \param_geometry +\param geometry2 \param_geometry +\return \return_check2{touch each other} + +\qbk{distinguish,two geometries} +\qbk{[include reference/algorithms/touches.qbk]} + */ +template <typename Geometry1, typename Geometry2> +inline bool touches(Geometry1 const& geometry1, Geometry2 const& geometry2) +{ + concept::check<Geometry1 const>(); + concept::check<Geometry2 const>(); + + + typedef detail::overlay::turn_info + < + typename geometry::point_type<Geometry1>::type + > turn_info; + + typedef detail::overlay::get_turn_info + < + typename point_type<Geometry1>::type, + typename point_type<Geometry2>::type, + turn_info, + detail::overlay::assign_null_policy + > policy_type; + + std::deque<turn_info> turns; + detail::get_turns::no_interrupt_policy policy; + boost::geometry::get_turns + < + false, false, + detail::overlay::assign_null_policy + >(geometry1, geometry2, turns, policy); + + return detail::touches::has_only_turns(turns) + && ! geometry::detail::disjoint::rings_containing(geometry1, geometry2) + && ! geometry::detail::disjoint::rings_containing(geometry2, geometry1) + ; +} + + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_ALGORITHMS_TOUCHES_HPP diff --git a/boost/geometry/algorithms/within.hpp b/boost/geometry/algorithms/within.hpp index 0e0cb68d53..f1f0993d77 100644 --- a/boost/geometry/algorithms/within.hpp +++ b/boost/geometry/algorithms/within.hpp @@ -65,7 +65,7 @@ struct point_in_ring static inline int apply(Point const& point, Ring const& ring, Strategy const& strategy) { - if (boost::size(ring) + if (int(boost::size(ring)) < core_detail::closure::minimum_ring_size<Closure>::value) { return -1; diff --git a/boost/geometry/core/access.hpp b/boost/geometry/core/access.hpp index a463b8cdc3..cdc8ff9cff 100644 --- a/boost/geometry/core/access.hpp +++ b/boost/geometry/core/access.hpp @@ -18,12 +18,13 @@ #include <cstddef> #include <boost/mpl/assert.hpp> -#include <boost/type_traits/remove_const.hpp> #include <boost/concept_check.hpp> +#include <boost/type_traits/is_pointer.hpp> #include <boost/geometry/core/coordinate_type.hpp> #include <boost/geometry/core/point_type.hpp> #include <boost/geometry/core/tag.hpp> +#include <boost/geometry/util/bare_type.hpp> namespace boost { namespace geometry @@ -79,6 +80,52 @@ struct indexed_access {}; } // namespace traits +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + +template +< + typename Geometry, + typename CoordinateType, + std::size_t Index, + std::size_t Dimension +> +struct indexed_access_non_pointer +{ + static inline CoordinateType get(Geometry const& geometry) + { + return traits::indexed_access<Geometry, Index, Dimension>::get(geometry); + } + static inline void set(Geometry& b, CoordinateType const& value) + { + traits::indexed_access<Geometry, Index, Dimension>::set(b, value); + } +}; + +template +< + typename Geometry, + typename CoordinateType, + std::size_t Index, + std::size_t Dimension +> +struct indexed_access_pointer +{ + static inline CoordinateType get(Geometry const* geometry) + { + return traits::indexed_access<typename boost::remove_pointer<Geometry>::type, Index, Dimension>::get(*geometry); + } + static inline void set(Geometry* geometry, CoordinateType const& value) + { + traits::indexed_access<typename boost::remove_pointer<Geometry>::type, Index, Dimension>::set(*geometry, value); + } +}; + + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + #ifndef DOXYGEN_NO_DISPATCH namespace core_dispatch @@ -89,7 +136,9 @@ template typename Tag, typename Geometry, typename - CoordinateType, std::size_t Dimension + CoordinateType, + std::size_t Dimension, + typename IsPointer > struct access { @@ -103,7 +152,8 @@ template typename Geometry, typename CoordinateType, std::size_t Index, - std::size_t Dimension + std::size_t Dimension, + typename IsPointer > struct indexed_access { @@ -112,7 +162,7 @@ struct indexed_access }; template <typename Point, typename CoordinateType, std::size_t Dimension> -struct access<point_tag, Point, CoordinateType, Dimension> +struct access<point_tag, Point, CoordinateType, Dimension, boost::false_type> { static inline CoordinateType get(Point const& point) { @@ -124,6 +174,20 @@ struct access<point_tag, Point, CoordinateType, Dimension> } }; +template <typename Point, typename CoordinateType, std::size_t Dimension> +struct access<point_tag, Point, CoordinateType, Dimension, boost::true_type> +{ + static inline CoordinateType get(Point const* point) + { + return traits::access<typename boost::remove_pointer<Point>::type, Dimension>::get(*point); + } + static inline void set(Point* p, CoordinateType const& value) + { + traits::access<typename boost::remove_pointer<Point>::type, Dimension>::set(*p, value); + } +}; + + template < typename Box, @@ -131,17 +195,21 @@ template std::size_t Index, std::size_t Dimension > -struct indexed_access<box_tag, Box, CoordinateType, Index, Dimension> -{ - static inline CoordinateType get(Box const& box) - { - return traits::indexed_access<Box, Index, Dimension>::get(box); - } - static inline void set(Box& b, CoordinateType const& value) - { - traits::indexed_access<Box, Index, Dimension>::set(b, value); - } -}; +struct indexed_access<box_tag, Box, CoordinateType, Index, Dimension, boost::false_type> + : detail::indexed_access_non_pointer<Box, CoordinateType, Index, Dimension> +{}; + +template +< + typename Box, + typename CoordinateType, + std::size_t Index, + std::size_t Dimension +> +struct indexed_access<box_tag, Box, CoordinateType, Index, Dimension, boost::true_type> + : detail::indexed_access_pointer<Box, CoordinateType, Index, Dimension> +{}; + template < @@ -150,17 +218,21 @@ template std::size_t Index, std::size_t Dimension > -struct indexed_access<segment_tag, Segment, CoordinateType, Index, Dimension> -{ - static inline CoordinateType get(Segment const& segment) - { - return traits::indexed_access<Segment, Index, Dimension>::get(segment); - } - static inline void set(Segment& segment, CoordinateType const& value) - { - traits::indexed_access<Segment, Index, Dimension>::set(segment, value); - } -}; +struct indexed_access<segment_tag, Segment, CoordinateType, Index, Dimension, boost::false_type> + : detail::indexed_access_non_pointer<Segment, CoordinateType, Index, Dimension> +{}; + + +template +< + typename Segment, + typename CoordinateType, + std::size_t Index, + std::size_t Dimension +> +struct indexed_access<segment_tag, Segment, CoordinateType, Index, Dimension, boost::true_type> + : detail::indexed_access_pointer<Segment, CoordinateType, Index, Dimension> +{}; } // namespace core_dispatch #endif // DOXYGEN_NO_DISPATCH @@ -199,14 +271,13 @@ inline typename coordinate_type<Geometry>::type get(Geometry const& geometry { boost::ignore_unused_variable_warning(dummy); - typedef typename boost::remove_const<Geometry>::type ncg_type; - typedef core_dispatch::access < typename tag<Geometry>::type, - ncg_type, - typename coordinate_type<ncg_type>::type, - Dimension + typename geometry::util::bare_type<Geometry>::type, + typename coordinate_type<Geometry>::type, + Dimension, + typename boost::is_pointer<Geometry>::type > coord_access_type; return coord_access_type::get(geometry); @@ -234,14 +305,13 @@ inline void set(Geometry& geometry { boost::ignore_unused_variable_warning(dummy); - typedef typename boost::remove_const<Geometry>::type ncg_type; - typedef core_dispatch::access < typename tag<Geometry>::type, - ncg_type, - typename coordinate_type<ncg_type>::type, - Dimension + typename geometry::util::bare_type<Geometry>::type, + typename coordinate_type<Geometry>::type, + Dimension, + typename boost::is_pointer<Geometry>::type > coord_access_type; coord_access_type::set(geometry, value); @@ -269,15 +339,14 @@ inline typename coordinate_type<Geometry>::type get(Geometry const& geometry { boost::ignore_unused_variable_warning(dummy); - typedef typename boost::remove_const<Geometry>::type ncg_type; - typedef core_dispatch::indexed_access < typename tag<Geometry>::type, - ncg_type, - typename coordinate_type<ncg_type>::type, + typename geometry::util::bare_type<Geometry>::type, + typename coordinate_type<Geometry>::type, Index, - Dimension + Dimension, + typename boost::is_pointer<Geometry>::type > coord_access_type; return coord_access_type::get(geometry); @@ -306,14 +375,14 @@ inline void set(Geometry& geometry { boost::ignore_unused_variable_warning(dummy); - typedef typename boost::remove_const<Geometry>::type ncg_type; - typedef core_dispatch::indexed_access < - typename tag<Geometry>::type, ncg_type, - typename coordinate_type<ncg_type>::type, + typename tag<Geometry>::type, + typename geometry::util::bare_type<Geometry>::type, + typename coordinate_type<Geometry>::type, Index, - Dimension + Dimension, + typename boost::is_pointer<Geometry>::type > coord_access_type; coord_access_type::set(geometry, value); diff --git a/boost/geometry/core/coordinate_dimension.hpp b/boost/geometry/core/coordinate_dimension.hpp index 84c11e2a55..2f3c914601 100644 --- a/boost/geometry/core/coordinate_dimension.hpp +++ b/boost/geometry/core/coordinate_dimension.hpp @@ -20,9 +20,9 @@ #include <boost/mpl/assert.hpp> #include <boost/mpl/equal_to.hpp> #include <boost/static_assert.hpp> -#include <boost/type_traits/remove_const.hpp> #include <boost/geometry/core/point_type.hpp> +#include <boost/geometry/util/bare_type.hpp> namespace boost { namespace geometry { @@ -58,7 +58,7 @@ template <typename T, typename G> struct dimension : dimension<point_tag, typename point_type<T, G>::type> {}; template <typename P> -struct dimension<point_tag, P> : traits::dimension<P> {}; +struct dimension<point_tag, P> : traits::dimension<typename geometry::util::bare_type<P>::type> {}; } // namespace core_dispatch #endif @@ -75,7 +75,7 @@ struct dimension : core_dispatch::dimension < typename tag<Geometry>::type, - typename boost::remove_const<Geometry>::type + typename geometry::util::bare_type<Geometry>::type > {}; diff --git a/boost/geometry/core/coordinate_system.hpp b/boost/geometry/core/coordinate_system.hpp index c23b8af15c..01e5ad1dd7 100644 --- a/boost/geometry/core/coordinate_system.hpp +++ b/boost/geometry/core/coordinate_system.hpp @@ -16,8 +16,9 @@ #include <boost/mpl/assert.hpp> -#include <boost/type_traits/remove_const.hpp> + #include <boost/geometry/core/point_type.hpp> +#include <boost/geometry/util/bare_type.hpp> namespace boost { namespace geometry @@ -61,10 +62,13 @@ namespace core_dispatch }; - template <typename P> - struct coordinate_system<point_tag, P> + template <typename Point> + struct coordinate_system<point_tag, Point> { - typedef typename traits::coordinate_system<P>::type type; + typedef typename traits::coordinate_system + < + typename geometry::util::bare_type<Point>::type + >::type type; }; @@ -82,11 +86,10 @@ namespace core_dispatch template <typename Geometry> struct coordinate_system { - typedef typename boost::remove_const<Geometry>::type ncg; typedef typename core_dispatch::coordinate_system < typename tag<Geometry>::type, - ncg + typename geometry::util::bare_type<Geometry>::type >::type type; }; diff --git a/boost/geometry/core/coordinate_type.hpp b/boost/geometry/core/coordinate_type.hpp index 0f901d3f19..d17f66ab73 100644 --- a/boost/geometry/core/coordinate_type.hpp +++ b/boost/geometry/core/coordinate_type.hpp @@ -16,9 +16,10 @@ #include <boost/mpl/assert.hpp> -#include <boost/type_traits/remove_const.hpp> +#include <boost/geometry/core/tag.hpp> #include <boost/geometry/core/point_type.hpp> +#include <boost/geometry/util/bare_type.hpp> #include <boost/geometry/util/promote_floating_point.hpp> @@ -63,12 +64,17 @@ struct coordinate_type template <typename Point> struct coordinate_type<point_tag, Point> { - typedef typename traits::coordinate_type<Point>::type type; + typedef typename traits::coordinate_type + < + typename geometry::util::bare_type<Point>::type + >::type type; }; + } // namespace core_dispatch #endif // DOXYGEN_NO_DISPATCH + /*! \brief \brief_meta{type, coordinate type (int\, float\, double\, etc), \meta_point_type} \tparam Geometry \tparam_geometry @@ -79,11 +85,11 @@ struct coordinate_type<point_tag, Point> template <typename Geometry> struct coordinate_type { - typedef typename core_dispatch::coordinate_type - < - typename tag<Geometry>::type, - typename boost::remove_const<Geometry>::type - >::type type; + typedef typename core_dispatch::coordinate_type + < + typename tag<Geometry>::type, + typename geometry::util::bare_type<Geometry>::type + >::type type; }; template <typename Geometry> diff --git a/boost/geometry/core/point_type.hpp b/boost/geometry/core/point_type.hpp index f49a43c56a..e148c84a57 100644 --- a/boost/geometry/core/point_type.hpp +++ b/boost/geometry/core/point_type.hpp @@ -22,6 +22,7 @@ #include <boost/geometry/core/ring_type.hpp> #include <boost/geometry/core/tag.hpp> #include <boost/geometry/core/tags.hpp> +#include <boost/geometry/util/bare_type.hpp> namespace boost { namespace geometry { @@ -115,11 +116,10 @@ struct point_type<polygon_tag, Polygon> template <typename Geometry> struct point_type { - typedef typename boost::remove_const<Geometry>::type ncg; typedef typename core_dispatch::point_type < typename tag<Geometry>::type, - ncg + typename boost::geometry::util::bare_type<Geometry>::type >::type type; }; diff --git a/boost/geometry/core/tag.hpp b/boost/geometry/core/tag.hpp index 5a99c97071..4c790fdc9e 100644 --- a/boost/geometry/core/tag.hpp +++ b/boost/geometry/core/tag.hpp @@ -16,9 +16,9 @@ #include <boost/mpl/assert.hpp> -#include <boost/type_traits/remove_const.hpp> #include <boost/geometry/core/tags.hpp> +#include <boost/geometry/util/bare_type.hpp> namespace boost { namespace geometry @@ -47,6 +47,7 @@ struct tag } // namespace traits + /*! \brief \brief_meta{type, tag, \meta_geometry_type} \details With Boost.Geometry, tags are the driving force of the tag dispatching @@ -61,7 +62,7 @@ struct tag { typedef typename traits::tag < - typename boost::remove_const<Geometry>::type + typename geometry::util::bare_type<Geometry>::type >::type type; }; diff --git a/boost/geometry/geometry.hpp b/boost/geometry/geometry.hpp index 22d08e2a45..b696b652cf 100644 --- a/boost/geometry/geometry.hpp +++ b/boost/geometry/geometry.hpp @@ -64,6 +64,7 @@ #include <boost/geometry/algorithms/reverse.hpp> #include <boost/geometry/algorithms/simplify.hpp> #include <boost/geometry/algorithms/sym_difference.hpp> +#include <boost/geometry/algorithms/touches.hpp> #include <boost/geometry/algorithms/transform.hpp> #include <boost/geometry/algorithms/union.hpp> #include <boost/geometry/algorithms/unique.hpp> diff --git a/boost/geometry/iterators/ever_circling_iterator.hpp b/boost/geometry/iterators/ever_circling_iterator.hpp index 03921096e7..566669e26d 100644 --- a/boost/geometry/iterators/ever_circling_iterator.hpp +++ b/boost/geometry/iterators/ever_circling_iterator.hpp @@ -95,67 +95,115 @@ private: bool m_skip_first; }; - - template <typename Range> -class ever_circling_range_iterator - : public boost::iterator_adaptor - < - ever_circling_range_iterator<Range>, - typename boost::range_iterator<Range>::type - > +struct ever_circling_range_iterator + : public boost::iterator_facade + < + ever_circling_range_iterator<Range>, + typename boost::range_value<Range>::type const, + boost::random_access_traversal_tag + > { -public : - typedef typename boost::range_iterator<Range>::type iterator_type; + /// Constructor including the range it is based on + explicit inline ever_circling_range_iterator(Range& range) + : m_range(&range) + , m_iterator(boost::begin(range)) + , m_size(boost::size(range)) + , m_index(0) + {} + + /// Default constructor + explicit inline ever_circling_range_iterator() + : m_range(NULL) + , m_size(0) + , m_index(0) + {} + + inline ever_circling_range_iterator<Range>& operator=(ever_circling_range_iterator<Range> const& source) + { + m_range = source.m_range; + m_iterator = source.m_iterator; + m_size = source.m_size; + m_index = source.m_index; + return *this; + } - explicit inline ever_circling_range_iterator(Range& range, - bool skip_first = false) - : m_range(range) - , m_skip_first(skip_first) + typedef std::ptrdiff_t difference_type; + +private: + friend class boost::iterator_core_access; + + inline typename boost::range_value<Range>::type const& dereference() const { - this->base_reference() = boost::begin(m_range); + return *m_iterator; } - explicit inline ever_circling_range_iterator(Range& range, iterator_type start, - bool skip_first = false) - : m_range(range) - , m_skip_first(skip_first) + inline difference_type distance_to(ever_circling_range_iterator<Range> const& other) const { - this->base_reference() = start; + return other.m_index - this->m_index; } - /// Navigate to a certain position, should be in [start .. end], if at end - /// it will circle again. - inline void moveto(iterator_type it) + inline bool equal(ever_circling_range_iterator<Range> const& other) const { - this->base_reference() = it; - check_end(); + return this->m_range == other.m_range + && this->m_index == other.m_index; } -private: + inline void increment() + { + ++m_index; + if (m_index >= 0 && m_index < m_size) + { + ++m_iterator; + } + else + { + update_iterator(); + } + } - friend class boost::iterator_core_access; + inline void decrement() + { + --m_index; + if (m_index >= 0 && m_index < m_size) + { + --m_iterator; + } + else + { + update_iterator(); + } + } - inline void increment(bool possibly_skip = true) + inline void advance(difference_type n) { - (this->base_reference())++; - check_end(possibly_skip); + if (m_index >= 0 && m_index < m_size + && m_index + n >= 0 && m_index + n < m_size) + { + m_index += n; + m_iterator += n; + } + else + { + m_index += n; + update_iterator(); + } } - inline void check_end(bool possibly_skip = true) + inline void update_iterator() { - if (this->base_reference() == boost::end(m_range)) + while (m_index < 0) { - this->base_reference() = boost::begin(m_range); - if (m_skip_first && possibly_skip) - { - increment(false); - } + m_index += m_size; } + m_index = m_index % m_size; + this->m_iterator = boost::begin(*m_range) + m_index; } - Range& m_range; - bool m_skip_first; + Range* m_range; + typename boost::range_iterator<Range>::type m_iterator; + difference_type m_size; + difference_type m_index; }; diff --git a/boost/geometry/multi/algorithms/detail/overlay/copy_segment_point.hpp b/boost/geometry/multi/algorithms/detail/overlay/copy_segment_point.hpp index 1093a84569..72be5ddbf4 100644 --- a/boost/geometry/multi/algorithms/detail/overlay/copy_segment_point.hpp +++ b/boost/geometry/multi/algorithms/detail/overlay/copy_segment_point.hpp @@ -42,7 +42,7 @@ struct copy_segment_point_multi BOOST_ASSERT ( seg_id.multi_index >= 0 - && seg_id.multi_index < boost::size(multi) + && seg_id.multi_index < int(boost::size(multi)) ); // Call the single-version diff --git a/boost/geometry/multi/algorithms/detail/overlay/copy_segments.hpp b/boost/geometry/multi/algorithms/detail/overlay/copy_segments.hpp index a234e240da..f474b12fdf 100644 --- a/boost/geometry/multi/algorithms/detail/overlay/copy_segments.hpp +++ b/boost/geometry/multi/algorithms/detail/overlay/copy_segments.hpp @@ -44,7 +44,7 @@ struct copy_segments_multi BOOST_ASSERT ( seg_id.multi_index >= 0 - && seg_id.multi_index < boost::size(multi_geometry) + && seg_id.multi_index < int(boost::size(multi_geometry)) ); // Call the single-version diff --git a/boost/geometry/multi/algorithms/detail/overlay/get_ring.hpp b/boost/geometry/multi/algorithms/detail/overlay/get_ring.hpp index 8cdd7ed970..e23acf99b3 100644 --- a/boost/geometry/multi/algorithms/detail/overlay/get_ring.hpp +++ b/boost/geometry/multi/algorithms/detail/overlay/get_ring.hpp @@ -35,7 +35,7 @@ struct get_ring<multi_polygon_tag> BOOST_ASSERT ( id.multi_index >= 0 - && id.multi_index < boost::size(multi_polygon) + && id.multi_index < int(boost::size(multi_polygon)) ); return get_ring<polygon_tag>::apply(id, multi_polygon[id.multi_index]); diff --git a/boost/geometry/multi/algorithms/detail/point_on_border.hpp b/boost/geometry/multi/algorithms/detail/point_on_border.hpp index 04d76422c6..ac462ed4c5 100644 --- a/boost/geometry/multi/algorithms/detail/point_on_border.hpp +++ b/boost/geometry/multi/algorithms/detail/point_on_border.hpp @@ -31,13 +31,13 @@ namespace detail { namespace point_on_border template < - typename MultiGeometry, typename Point, + typename MultiGeometry, typename Policy > struct point_on_multi { - static inline bool apply(MultiGeometry const& multi, Point& point) + static inline bool apply(Point& point, MultiGeometry const& multi, bool midpoint) { // Take a point on the first multi-geometry // (i.e. the first that is not empty) @@ -48,7 +48,7 @@ struct point_on_multi it != boost::end(multi); ++it) { - if (Policy::apply(*it, point)) + if (Policy::apply(point, *it, midpoint)) { return true; } @@ -69,16 +69,16 @@ namespace dispatch { -template<typename Multi, typename Point> -struct point_on_border<multi_polygon_tag, Multi, Point> +template<typename Point, typename Multi> +struct point_on_border<multi_polygon_tag, Point, Multi> : detail::point_on_border::point_on_multi < - Multi, Point, + Multi, detail::point_on_border::point_on_polygon < - typename boost::range_value<Multi>::type, - Point + Point, + typename boost::range_value<Multi>::type > > {}; diff --git a/boost/geometry/multi/algorithms/detail/sections/range_by_section.hpp b/boost/geometry/multi/algorithms/detail/sections/range_by_section.hpp index cf4bf5d83b..28a4805ed1 100644 --- a/boost/geometry/multi/algorithms/detail/sections/range_by_section.hpp +++ b/boost/geometry/multi/algorithms/detail/sections/range_by_section.hpp @@ -46,7 +46,7 @@ struct full_section_multi BOOST_ASSERT ( section.ring_id.multi_index >= 0 - && section.ring_id.multi_index < boost::size(multi) + && section.ring_id.multi_index < int(boost::size(multi)) ); return Policy::apply(multi[section.ring_id.multi_index], section); diff --git a/boost/geometry/multi/io/wkt/read.hpp b/boost/geometry/multi/io/wkt/read.hpp index 3c75fd8040..218ddf9999 100644 --- a/boost/geometry/multi/io/wkt/read.hpp +++ b/boost/geometry/multi/io/wkt/read.hpp @@ -58,6 +58,69 @@ struct multi_parser handle_close_parenthesis(it, tokens.end(), wkt); } + + check_end(it, tokens.end(), wkt); + } +}; + +template <typename P> +struct noparenthesis_point_parser +{ + static inline void apply(tokenizer::iterator& it, tokenizer::iterator end, + std::string const& wkt, P& point) + { + parsing_assigner<P, 0, dimension<P>::value>::apply(it, end, point, wkt); + } +}; + +template <typename MultiGeometry, typename PrefixPolicy> +struct multi_point_parser +{ + static inline void apply(std::string const& wkt, MultiGeometry& geometry) + { + traits::clear<MultiGeometry>::apply(geometry); + + tokenizer tokens(wkt, boost::char_separator<char>(" ", ",()")); + tokenizer::iterator it; + + if (initialize<MultiGeometry>(tokens, PrefixPolicy::apply(), wkt, it)) + { + handle_open_parenthesis(it, tokens.end(), wkt); + + // If first point definition starts with "(" then parse points as (x y) + // otherwise as "x y" + bool using_brackets = (it != tokens.end() && *it == "("); + + while(it != tokens.end() && *it != ")") + { + traits::resize<MultiGeometry>::apply(geometry, boost::size(geometry) + 1); + + if (using_brackets) + { + point_parser + < + typename boost::range_value<MultiGeometry>::type + >::apply(it, tokens.end(), wkt, geometry.back()); + } + else + { + noparenthesis_point_parser + < + typename boost::range_value<MultiGeometry>::type + >::apply(it, tokens.end(), wkt, geometry.back()); + } + + if (it != tokens.end() && *it == ",") + { + // Skip "," after point is parsed + ++it; + } + } + + handle_close_parenthesis(it, tokens.end(), wkt); + } + + check_end(it, tokens.end(), wkt); } }; @@ -69,10 +132,9 @@ namespace dispatch template <typename MultiGeometry> struct read_wkt<multi_point_tag, MultiGeometry> - : detail::wkt::multi_parser + : detail::wkt::multi_point_parser < MultiGeometry, - detail::wkt::point_parser, detail::wkt::prefix_multipoint > {}; diff --git a/boost/geometry/policies/relate/direction.hpp b/boost/geometry/policies/relate/direction.hpp index 9090d8b51e..cfbaf7dd15 100644 --- a/boost/geometry/policies/relate/direction.hpp +++ b/boost/geometry/policies/relate/direction.hpp @@ -126,7 +126,9 @@ struct segments_direction typedef typename select_most_precise<coordinate_type, double>::type rtype; + template <typename R> static inline return_type segments_intersect(side_info const& sides, + R const&, coordinate_type const& dx1, coordinate_type const& dy1, coordinate_type const& dx2, coordinate_type const& dy2, S1 const& s1, S2 const& s2) diff --git a/boost/geometry/policies/relate/intersection_points.hpp b/boost/geometry/policies/relate/intersection_points.hpp index afd1dde501..d7d5199051 100644 --- a/boost/geometry/policies/relate/intersection_points.hpp +++ b/boost/geometry/policies/relate/intersection_points.hpp @@ -36,56 +36,33 @@ struct segments_intersection_points typedef ReturnType return_type; typedef S1 segment_type1; typedef S2 segment_type2; + typedef typename select_calculation_type < S1, S2, CalculationType >::type coordinate_type; + template <typename R> static inline return_type segments_intersect(side_info const&, + R const& r, coordinate_type const& dx1, coordinate_type const& dy1, coordinate_type const& dx2, coordinate_type const& dy2, S1 const& s1, S2 const& s2) { - return_type result; typedef typename geometry::coordinate_type < typename return_type::point_type - >::type coordinate_type; + >::type return_coordinate_type; - // Get the same type, but at least a double (also used for divisions) - typedef typename select_most_precise - < - coordinate_type, double - >::type promoted_type; - - promoted_type const s1x = get<0, 0>(s1); - promoted_type const s1y = get<0, 1>(s1); - - // Calculate other determinants - Cramers rule - promoted_type const wx = get<0, 0>(s1) - get<0, 0>(s2); - promoted_type const wy = get<0, 1>(s1) - get<0, 1>(s2); - promoted_type const d = detail::determinant<promoted_type>(dx1, dy1, dx2, dy2); - promoted_type const da = detail::determinant<promoted_type>(dx2, dy2, wx, wy); - - // r: ratio 0-1 where intersection divides A/B - promoted_type r = da / d; - promoted_type const zero = promoted_type(); - promoted_type const one = 1; - // Handle robustness issues - if (r < zero) - { - r = zero; - } - else if (r > one) - { - r = one; - } + coordinate_type const s1x = get<0, 0>(s1); + coordinate_type const s1y = get<0, 1>(s1); + return_type result; result.count = 1; set<0>(result.intersections[0], - boost::numeric_cast<coordinate_type>(s1x + r * promoted_type(dx1))); + boost::numeric_cast<return_coordinate_type>(R(s1x) + r * R(dx1))); set<1>(result.intersections[0], - boost::numeric_cast<coordinate_type>(s1y + r * promoted_type(dy1))); + boost::numeric_cast<return_coordinate_type>(R(s1y) + r * R(dy1))); return result; } diff --git a/boost/geometry/policies/relate/tupled.hpp b/boost/geometry/policies/relate/tupled.hpp index 853a8a26a8..f6713c44c9 100644 --- a/boost/geometry/policies/relate/tupled.hpp +++ b/boost/geometry/policies/relate/tupled.hpp @@ -49,16 +49,18 @@ struct segments_tupled // Get the same type, but at least a double typedef typename select_most_precise<coordinate_type, double>::type rtype; + template <typename R> static inline return_type segments_intersect(side_info const& sides, + R const& r, coordinate_type const& dx1, coordinate_type const& dy1, coordinate_type const& dx2, coordinate_type const& dy2, segment_type1 const& s1, segment_type2 const& s2) { return boost::make_tuple ( - Policy1::segments_intersect(sides, + Policy1::segments_intersect(sides, r, dx1, dy1, dx2, dy2, s1, s2), - Policy2::segments_intersect(sides, + Policy2::segments_intersect(sides, r, dx1, dy1, dx2, dy2, s1, s2) ); } diff --git a/boost/geometry/strategies/cartesian/cart_intersect.hpp b/boost/geometry/strategies/cartesian/cart_intersect.hpp index 2bf7b77676..ea92cf37b0 100644 --- a/boost/geometry/strategies/cartesian/cart_intersect.hpp +++ b/boost/geometry/strategies/cartesian/cart_intersect.hpp @@ -41,21 +41,17 @@ namespace strategy { namespace intersection namespace detail { -template <typename Segment, size_t Dimension> -struct segment_arrange +template <std::size_t Dimension, typename Segment, typename T> +static inline void segment_arrange(Segment const& s, T& s_1, T& s_2, bool& swapped) { - template <typename T> - static inline void apply(Segment const& s, T& s_1, T& s_2, bool& swapped) + s_1 = get<0, Dimension>(s); + s_2 = get<1, Dimension>(s); + if (s_1 > s_2) { - s_1 = get<0, Dimension>(s); - s_2 = get<1, Dimension>(s); - if (s_1 > s_2) - { - std::swap(s_1, s_2); - swapped = true; - } + std::swap(s_1, s_2); + swapped = true; } -}; +} template <std::size_t Index, typename Segment> inline typename geometry::point_type<Segment>::type get_from_index( @@ -121,34 +117,20 @@ struct relate_cartesian_segments coordinate_type const& dx_a, coordinate_type const& dy_a, coordinate_type const& dx_b, coordinate_type const& dy_b) { - // 1) Handle "disjoint", common case. - // per dimension, 2 cases: a_1----------a_2 b_1-------b_2 or B left of A - coordinate_type ax_1, ax_2, bx_1, bx_2; - bool ax_swapped = false, bx_swapped = false; - detail::segment_arrange<segment_type1, 0>::apply(a, ax_1, ax_2, ax_swapped); - detail::segment_arrange<segment_type2, 0>::apply(b, bx_1, bx_2, bx_swapped); - if (ax_2 < bx_1 || ax_1 > bx_2) - { - return Policy::disjoint(); - } - - // 1b) In Y-dimension - coordinate_type ay_1, ay_2, by_1, by_2; - bool ay_swapped = false, by_swapped = false; - detail::segment_arrange<segment_type1, 1>::apply(a, ay_1, ay_2, ay_swapped); - detail::segment_arrange<segment_type2, 1>::apply(b, by_1, by_2, by_swapped); - if (ay_2 < by_1 || ay_1 > by_2) - { - return Policy::disjoint(); - } - typedef side::side_by_triangle<coordinate_type> side; side_info sides; - // 2) Calculate sides - // Note: Do NOT yet calculate the determinant here, but use the SIDE strategy. - // Determinant calculation is not robust; side (orient) can be made robust - // (and is much robuster even without measures) + bool collinear_use_first = math::abs(dx_a) + math::abs(dx_b) >= math::abs(dy_a) + math::abs(dy_b); + + sides.set<0> + ( + side::apply(detail::get_from_index<0>(b) + , detail::get_from_index<1>(b) + , detail::get_from_index<0>(a)), + side::apply(detail::get_from_index<0>(b) + , detail::get_from_index<1>(b) + , detail::get_from_index<1>(a)) + ); sides.set<1> ( side::apply(detail::get_from_index<0>(a) @@ -159,26 +141,18 @@ struct relate_cartesian_segments , detail::get_from_index<1>(b)) ); - if (sides.same<1>()) - { - // Both points are at same side of other segment, we can leave - return Policy::disjoint(); - } + bool collinear = sides.collinear(); - // 2b) For other segment - sides.set<0> - ( - side::apply(detail::get_from_index<0>(b) - , detail::get_from_index<1>(b) - , detail::get_from_index<0>(a)), - side::apply(detail::get_from_index<0>(b) - , detail::get_from_index<1>(b) - , detail::get_from_index<1>(a)) - ); + robustness_verify_collinear(a, b, sides, collinear); + robustness_verify_meeting(a, b, sides, collinear, collinear_use_first); - if (sides.same<0>()) + if (sides.same<0>() || sides.same<1>()) { - return Policy::disjoint(); + // Both points are at same side of other segment, we can leave + if (robustness_verify_same_side(a, b, sides)) + { + return Policy::disjoint(); + } } // Degenerate cases: segments of single point, lying on other segment, non disjoint @@ -192,82 +166,343 @@ struct relate_cartesian_segments return Policy::degenerate(b, false); } - bool collinear = sides.collinear(); - - // Get the same type, but at least a double (also used for divisions) typedef typename select_most_precise < coordinate_type, double >::type promoted_type; - // Calculate the determinant/2D cross product - // (Note, because we only check on zero, - // the order a/b does not care) - promoted_type const d = geometry::detail::determinant - < - promoted_type - >(dx_a, dy_a, dx_b, dy_b); - - // Determinant d should be nonzero. - // If it is zero, we have an robustness issue here, - // (and besides that we cannot divide by it) - promoted_type const pt_zero = promoted_type(); - if(math::equals(d, pt_zero) && ! collinear) + // r: ratio 0-1 where intersection divides A/B + // (only calculated for non-collinear segments) + promoted_type r; + if (! collinear) { -#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION - std::cout << "Determinant zero? Type : " - << typeid(coordinate_type).name() - << std::endl; - - std::cout << " dx_a : " << dx_a << std::endl; - std::cout << " dy_a : " << dy_a << std::endl; - std::cout << " dx_b : " << dx_b << std::endl; - std::cout << " dy_b : " << dy_b << std::endl; - - std::cout << " side a <-> b.first : " << sides.get<0,0>() << std::endl; - std::cout << " side a <-> b.second: " << sides.get<0,1>() << std::endl; - std::cout << " side b <-> a.first : " << sides.get<1,0>() << std::endl; - std::cout << " side b <-> a.second: " << sides.get<1,1>() << std::endl; -#endif - - if (sides.as_collinear()) + // Calculate determinants - Cramers rule + coordinate_type const wx = get<0, 0>(a) - get<0, 0>(b); + coordinate_type const wy = get<0, 1>(a) - get<0, 1>(b); + coordinate_type const d = geometry::detail::determinant<coordinate_type>(dx_a, dy_a, dx_b, dy_b); + coordinate_type const da = geometry::detail::determinant<coordinate_type>(dx_b, dy_b, wx, wy); + + coordinate_type const zero = coordinate_type(); + if (math::equals(d, zero)) { + // This is still a collinear case (because of FP imprecision this can occur here) + // sides.debug(); sides.set<0>(0,0); sides.set<1>(0,0); collinear = true; } else { - return Policy::error("Determinant zero!"); + r = promoted_type(da) / promoted_type(d); + + if (! robustness_verify_r(a, b, r)) + { + return Policy::disjoint(); + } + + robustness_handle_meeting(a, b, sides, dx_a, dy_a, wx, wy, d, r); + + if (robustness_verify_disjoint_at_one_collinear(a, b, sides)) + { + return Policy::disjoint(); + } + } } if(collinear) { - // Segments are collinear. We'll find out how. - if (math::equals(dx_b, zero)) + if (collinear_use_first) { - // Vertical -> Check y-direction - return relate_collinear(a, b, - ay_1, ay_2, by_1, by_2, - ay_swapped, by_swapped); + return relate_collinear<0>(a, b); } else { - // Check x-direction - return relate_collinear(a, b, - ax_1, ax_2, bx_1, bx_2, - ax_swapped, bx_swapped); + // Y direction contains larger segments (maybe dx is zero) + return relate_collinear<1>(a, b); } } - return Policy::segments_intersect(sides, + return Policy::segments_intersect(sides, r, dx_a, dy_a, dx_b, dy_b, a, b); } private : + + // Ratio should lie between 0 and 1 + // Also these three conditions might be of FP imprecision, the segments were actually (nearly) collinear + template <typename T> + static inline bool robustness_verify_r( + segment_type1 const& a, segment_type2 const& b, + T& r) + { + T const zero = 0; + T const one = 1; + if (r < zero || r > one) + { + if (verify_disjoint<0>(a, b) || verify_disjoint<1>(a, b)) + { + // Can still be disjoint (even if not one is left or right from another) + // This is e.g. in case #snake4 of buffer test. + return false; + } + + //std::cout << "ROBUSTNESS: correction of r " << r << std::endl; + // sides.debug(); + + // ROBUSTNESS: the r value can in epsilon-cases much larger than 1, while (with perfect arithmetic) + // it should be one. It can be 1.14 or even 1.98049 or 2 (while still intersecting) + + // If segments are crossing (we can see that with the sides) + // and one is inside the other, there must be an intersection point. + // We correct for that. + // This is (only) in case #ggl_list_20110820_christophe in unit tests + + // If segments are touching (two sides zero), of course they should intersect + // This is (only) in case #buffer_rt_i in the unit tests) + + // If one touches in the middle, they also should intersect (#buffer_rt_j) + + // Note that even for ttmath r is occasionally > 1, e.g. 1.0000000000000000000000036191231203575 + + if (r > one) + { + r = one; + } + else if (r < zero) + { + r = zero; + } + } + return true; + } + + static inline void robustness_verify_collinear( + segment_type1 const& a, segment_type2 const& b, + side_info& sides, + bool& collinear) + { + if ((sides.zero<0>() && ! sides.zero<1>()) || (sides.zero<1>() && ! sides.zero<0>())) + { + // If one of the segments is collinear, the other must be as well. + // So handle it as collinear. + // (In float/double epsilon margins it can easily occur that one or two of them are -1/1) + // sides.debug(); + sides.set<0>(0,0); + sides.set<1>(0,0); + collinear = true; + } + } + + static inline void robustness_verify_meeting( + segment_type1 const& a, segment_type2 const& b, + side_info& sides, + bool& collinear, bool& collinear_use_first) + { + if (sides.meeting()) + { + // If two segments meet each other at their segment-points, two sides are zero, + // the other two are not (unless collinear but we don't mean those here). + // However, in near-epsilon ranges it can happen that two sides are zero + // but they do not meet at their segment-points. + // In that case they are nearly collinear and handled as such. + if (! point_equals + ( + select(sides.zero_index<0>(), a), + select(sides.zero_index<1>(), b) + ) + ) + { + sides.set<0>(0,0); + sides.set<1>(0,0); + collinear = true; + + if (collinear_use_first && analyse_equal<0>(a, b)) + { + collinear_use_first = false; + } + else if (! collinear_use_first && analyse_equal<1>(a, b)) + { + collinear_use_first = true; + } + + } + } + } + + // Verifies and if necessary correct missed touch because of robustness + // This is the case at multi_polygon_buffer unittest #rt_m + static inline bool robustness_verify_same_side( + segment_type1 const& a, segment_type2 const& b, + side_info& sides) + { + int corrected = 0; + if (sides.one_touching<0>()) + { + if (point_equals( + select(sides.zero_index<0>(), a), + select(0, b) + )) + { + sides.correct_to_zero<1, 0>(); + corrected = 1; + } + if (point_equals + ( + select(sides.zero_index<0>(), a), + select(1, b) + )) + { + sides.correct_to_zero<1, 1>(); + corrected = 2; + } + } + else if (sides.one_touching<1>()) + { + if (point_equals( + select(sides.zero_index<1>(), b), + select(0, a) + )) + { + sides.correct_to_zero<0, 0>(); + corrected = 3; + } + if (point_equals + ( + select(sides.zero_index<1>(), b), + select(1, a) + )) + { + sides.correct_to_zero<0, 1>(); + corrected = 4; + } + } + + return corrected == 0; + } + + static inline bool robustness_verify_disjoint_at_one_collinear( + segment_type1 const& a, segment_type2 const& b, + side_info const& sides) + { + if (sides.one_of_all_zero()) + { + if (verify_disjoint<0>(a, b) || verify_disjoint<1>(a, b)) + { + return true; + } + } + return false; + } + + + // If r is one, or zero, segments should meet and their endpoints. + // Robustness issue: check if this is really the case. + // It turns out to be no problem, see buffer test #rt_s1 (and there are many cases generated) + // It generates an "ends in the middle" situation which is correct. + template <typename T, typename R> + static inline void robustness_handle_meeting(segment_type1 const& a, segment_type2 const& b, + side_info& sides, + T const& dx_a, T const& dy_a, T const& wx, T const& wy, + T const& d, R const& r) + { + return; + + T const db = geometry::detail::determinant<T>(dx_a, dy_a, wx, wy); + + R const zero = 0; + R const one = 1; + if (math::equals(r, zero) || math::equals(r, one)) + { + R rb = db / d; + if (rb <= 0 || rb >= 1 || math::equals(rb, 0) || math::equals(rb, 1)) + { + if (sides.one_zero<0>() && ! sides.one_zero<1>()) // or vice versa + { +#if defined(BOOST_GEOMETRY_COUNT_INTERSECTION_EQUAL) + extern int g_count_intersection_equal; + g_count_intersection_equal++; +#endif + sides.debug(); + std::cout << "E r=" << r << " r.b=" << rb << " "; + } + } + } + } + + template <std::size_t Dimension> + static inline bool verify_disjoint(segment_type1 const& a, + segment_type2 const& b) + { + coordinate_type a_1, a_2, b_1, b_2; + bool a_swapped = false, b_swapped = false; + detail::segment_arrange<Dimension>(a, a_1, a_2, a_swapped); + detail::segment_arrange<Dimension>(b, b_1, b_2, b_swapped); + return math::smaller(a_2, b_1) || math::larger(a_1, b_2); + } + + template <typename Segment> + static inline typename point_type<Segment>::type select(int index, Segment const& segment) + { + return index == 0 + ? detail::get_from_index<0>(segment) + : detail::get_from_index<1>(segment) + ; + } + + // We cannot use geometry::equals here. Besides that this will be changed + // to compare segment-coordinate-values directly (not necessary to retrieve point first) + template <typename Point1, typename Point2> + static inline bool point_equals(Point1 const& point1, Point2 const& point2) + { + return math::equals(get<0>(point1), get<0>(point2)) + && math::equals(get<1>(point1), get<1>(point2)) + ; + } + + // We cannot use geometry::equals here. Besides that this will be changed + // to compare segment-coordinate-values directly (not necessary to retrieve point first) + template <typename Point1, typename Point2> + static inline bool point_equality(Point1 const& point1, Point2 const& point2, + bool& equals_0, bool& equals_1) + { + equals_0 = math::equals(get<0>(point1), get<0>(point2)); + equals_1 = math::equals(get<1>(point1), get<1>(point2)); + return equals_0 && equals_1; + } + + template <std::size_t Dimension> + static inline bool analyse_equal(segment_type1 const& a, segment_type2 const& b) + { + coordinate_type const a_1 = geometry::get<0, Dimension>(a); + coordinate_type const a_2 = geometry::get<1, Dimension>(a); + coordinate_type const b_1 = geometry::get<0, Dimension>(b); + coordinate_type const b_2 = geometry::get<1, Dimension>(b); + return math::equals(a_1, b_1) + || math::equals(a_2, b_1) + || math::equals(a_1, b_2) + || math::equals(a_2, b_2) + ; + } + + template <std::size_t Dimension> + static inline return_type relate_collinear(segment_type1 const& a, + segment_type2 const& b) + { + coordinate_type a_1, a_2, b_1, b_2; + bool a_swapped = false, b_swapped = false; + detail::segment_arrange<Dimension>(a, a_1, a_2, a_swapped); + detail::segment_arrange<Dimension>(b, b_1, b_2, b_swapped); + if (math::smaller(a_2, b_1) || math::larger(a_1, b_2)) + //if (a_2 < b_1 || a_1 > b_2) + { + return Policy::disjoint(); + } + return relate_collinear(a, b, a_1, a_2, b_1, b_2, a_swapped, b_swapped); + } + /// Relate segments known collinear static inline return_type relate_collinear(segment_type1 const& a , segment_type2 const& b @@ -296,33 +531,36 @@ private : // Handle "equal", in polygon neighbourhood comparisons a common case - // Check if segments are equal... - bool const a1_eq_b1 = math::equals(get<0, 0>(a), get<0, 0>(b)) - && math::equals(get<0, 1>(a), get<0, 1>(b)); - bool const a2_eq_b2 = math::equals(get<1, 0>(a), get<1, 0>(b)) - && math::equals(get<1, 1>(a), get<1, 1>(b)); - if (a1_eq_b1 && a2_eq_b2) - { - return Policy::segment_equal(a, false); - } + bool const opposite = a_swapped ^ b_swapped; + bool const both_swapped = a_swapped && b_swapped; + + // Check if segments are equal or opposite equal... + bool const swapped_a1_eq_b1 = math::equals(a_1, b_1); + bool const swapped_a2_eq_b2 = math::equals(a_2, b_2); - // ... or opposite equal - bool const a1_eq_b2 = math::equals(get<0, 0>(a), get<1, 0>(b)) - && math::equals(get<0, 1>(a), get<1, 1>(b)); - bool const a2_eq_b1 = math::equals(get<1, 0>(a), get<0, 0>(b)) - && math::equals(get<1, 1>(a), get<0, 1>(b)); - if (a1_eq_b2 && a2_eq_b1) + if (swapped_a1_eq_b1 && swapped_a2_eq_b2) { - return Policy::segment_equal(a, true); + return Policy::segment_equal(a, opposite); } + bool const swapped_a2_eq_b1 = math::equals(a_2, b_1); + bool const swapped_a1_eq_b2 = math::equals(a_1, b_2); + + bool const a1_eq_b1 = both_swapped ? swapped_a2_eq_b2 : a_swapped ? swapped_a2_eq_b1 : b_swapped ? swapped_a1_eq_b2 : swapped_a1_eq_b1; + bool const a2_eq_b2 = both_swapped ? swapped_a1_eq_b1 : a_swapped ? swapped_a1_eq_b2 : b_swapped ? swapped_a2_eq_b1 : swapped_a2_eq_b2; + + bool const a1_eq_b2 = both_swapped ? swapped_a2_eq_b1 : a_swapped ? swapped_a2_eq_b2 : b_swapped ? swapped_a1_eq_b1 : swapped_a1_eq_b2; + bool const a2_eq_b1 = both_swapped ? swapped_a1_eq_b2 : a_swapped ? swapped_a1_eq_b1 : b_swapped ? swapped_a2_eq_b2 : swapped_a2_eq_b1; + + + // The rest below will return one or two intersections. // The delegated class can decide which is the intersection point, or two, build the Intersection Matrix (IM) // For IM it is important to know which relates to which. So this information is given, // without performance penalties to intersection calculation - bool const has_common_points = a1_eq_b1 || a1_eq_b2 || a2_eq_b1 || a2_eq_b2; + bool const has_common_points = swapped_a1_eq_b1 || swapped_a1_eq_b2 || swapped_a2_eq_b1 || swapped_a2_eq_b2; // "Touch" -> one intersection point -> one but not two common points @@ -342,7 +580,7 @@ private : // #4: a2<----a1 b1--->b2 (no arrival at all) // Where the arranged forms have two forms: // a_1-----a_2/b_1-------b_2 or reverse (B left of A) - if (has_common_points && (math::equals(a_2, b_1) || math::equals(b_2, a_1))) + if ((swapped_a2_eq_b1 || swapped_a1_eq_b2) && ! swapped_a1_eq_b1 && ! swapped_a2_eq_b2) { if (a2_eq_b1) return Policy::collinear_touch(get<1, 0>(a), get<1, 1>(a), 0, -1); if (a1_eq_b2) return Policy::collinear_touch(get<0, 0>(a), get<0, 1>(a), -1, 0); @@ -402,7 +640,6 @@ private : if (a1_eq_b1) return Policy::collinear_interior_boundary_intersect(a_in_b ? a : b, a_in_b, arrival_a, -arrival_a, false); } - bool const opposite = a_swapped ^ b_swapped; // "Inside", a completely within b or b completely within a @@ -464,7 +701,6 @@ private : the picture might seem wrong but it (supposed to be) is right. */ - bool const both_swapped = a_swapped && b_swapped; if (b_1 < a_2 && a_2 < b_2) { // Left column, from bottom to top @@ -486,7 +722,7 @@ private : ; } // Nothing should goes through. If any we have made an error - // Robustness: it can occur here... + // std::cout << "Robustness issue, non-logical behaviour" << std::endl; return Policy::error("Robustness issue, non-logical behaviour"); } }; diff --git a/boost/geometry/strategies/cartesian/distance_projected_point.hpp b/boost/geometry/strategies/cartesian/distance_projected_point.hpp index e704a33105..13d4168445 100644 --- a/boost/geometry/strategies/cartesian/distance_projected_point.hpp +++ b/boost/geometry/strategies/cartesian/distance_projected_point.hpp @@ -75,23 +75,27 @@ template class projected_point { public : - typedef typename strategy::distance::services::return_type<Strategy>::type calculation_type; - -private : - // The three typedefs below are necessary to calculate distances // from segments defined in integer coordinates. // Integer coordinates can still result in FP distances. // There is a division, which must be represented in FP. // So promote. - typedef typename promote_floating_point<calculation_type>::type fp_type; + typedef typename promote_floating_point + < + typename strategy::distance::services::return_type + < + Strategy + >::type + >::type calculation_type; + +private : // A projected point of points in Integer coordinates must be able to be // represented in FP. typedef model::point < - fp_type, + calculation_type, dimension<PointOfSegment>::value, typename coordinate_system<PointOfSegment>::type > fp_point_type; @@ -139,19 +143,19 @@ public : boost::ignore_unused_variable_warning(strategy); calculation_type const zero = calculation_type(); - fp_type const c1 = dot_product(w, v); + calculation_type const c1 = dot_product(w, v); if (c1 <= zero) { return strategy.apply(p, p1); } - fp_type const c2 = dot_product(v, v); + calculation_type const c2 = dot_product(v, v); if (c2 <= c1) { return strategy.apply(p, p2); } // See above, c1 > 0 AND c2 > c1 so: c2 != 0 - fp_type const b = c1 / c2; + calculation_type const b = c1 / c2; fp_strategy_type fp_strategy = strategy::distance::services::get_similar diff --git a/boost/geometry/strategies/side_info.hpp b/boost/geometry/strategies/side_info.hpp index 3fc43b8d99..f3a9da0df0 100644 --- a/boost/geometry/strategies/side_info.hpp +++ b/boost/geometry/strategies/side_info.hpp @@ -45,6 +45,19 @@ public : } template <int Which, int Index> + inline void correct_to_zero() + { + if (Index == 0) + { + sides[Which].first = 0; + } + else + { + sides[Which].second = 0; + } + } + + template <int Which, int Index> inline int get() const { return Index == 0 ? sides[Which].first : sides[Which].second; @@ -67,21 +80,81 @@ public : && sides[1].second == 0; } - // If one of the segments is collinear, the other must be as well. - // So handle it as collinear. - // (In floating point margins it can occur that one of them is 1!) - inline bool as_collinear() const + inline bool crossing() const + { + return sides[0].first * sides[0].second == -1 + && sides[1].first * sides[1].second == -1; + } + + inline bool touching() const + { + return (sides[0].first * sides[1].first == -1 + && sides[0].second == 0 && sides[1].second == 0) + || (sides[1].first * sides[0].first == -1 + && sides[1].second == 0 && sides[0].second == 0); + } + + template <int Which> + inline bool one_touching() const + { + // This is normally a situation which can't occur: + // If one is completely left or right, the other cannot touch + return one_zero<Which>() + && sides[1 - Which].first * sides[1 - Which].second == 1; + } + + inline bool meeting() const + { + // Two of them (in each segment) zero, two not + return one_zero<0>() && one_zero<1>(); + } + + template <int Which> + inline bool zero() const { - return sides[0].first * sides[0].second == 0 - || sides[1].first * sides[1].second == 0; + return sides[Which].first == 0 && sides[Which].second == 0; } + template <int Which> + inline bool one_zero() const + { + return (sides[Which].first == 0 && sides[Which].second != 0) + || (sides[Which].first != 0 && sides[Which].second == 0); + } + + inline bool one_of_all_zero() const + { + int const sum = std::abs(sides[0].first) + + std::abs(sides[0].second) + + std::abs(sides[1].first) + + std::abs(sides[1].second); + return sum == 3; + } + + + template <int Which> + inline int zero_index() const + { + return sides[Which].first == 0 ? 0 : 1; + } + + + inline void debug() const + { + std::cout << sides[0].first << " " + << sides[0].second << " " + << sides[1].first << " " + << sides[1].second + << std::endl; + } + + inline void reverse() { std::swap(sides[0], sides[1]); } -private : +//private : std::pair<int, int> sides[2]; }; diff --git a/boost/geometry/strategies/spherical/distance_cross_track.hpp b/boost/geometry/strategies/spherical/distance_cross_track.hpp index ba589223ec..7b353020eb 100644 --- a/boost/geometry/strategies/spherical/distance_cross_track.hpp +++ b/boost/geometry/strategies/spherical/distance_cross_track.hpp @@ -101,23 +101,56 @@ public : { // http://williams.best.vwh.net/avform.htm#XTE return_type d1 = m_strategy.apply(sp1, p); + return_type d3 = m_strategy.apply(sp1, sp2); + + if (geometry::math::equals(d3, 0.0)) + { + // "Degenerate" segment, return either d1 or d2 + return d1; + } - // Actually, calculation of d2 not necessary if we know that the projected point is on the great circle... return_type d2 = m_strategy.apply(sp2, p); return_type crs_AD = course(sp1, p); return_type crs_AB = course(sp1, sp2); - return_type XTD = m_radius * geometry::math::abs(asin(sin(d1 / m_radius) * sin(crs_AD - crs_AB))); + return_type crs_BA = crs_AB - geometry::math::pi<return_type>(); + return_type crs_BD = course(sp2, p); + return_type d_crs1 = crs_AD - crs_AB; + return_type d_crs2 = crs_BD - crs_BA; + + // d1, d2, d3 are in principle not needed, only the sign matters + return_type projection1 = cos( d_crs1 ) * d1 / d3; + return_type projection2 = cos( d_crs2 ) * d2 / d3; #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK -std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " " << crs_AD * geometry::math::r2d << std::endl; -std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " " << crs_AB * geometry::math::r2d << std::endl; -std::cout << "XTD: " << XTD << " d1: " << d1 << " d2: " << d2 << std::endl; + std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " " << crs_AD * geometry::math::r2d << std::endl; + std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " " << crs_AB * geometry::math::r2d << std::endl; + std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " " << crs_BD * geometry::math::r2d << std::endl; + std::cout << "Projection AD-AB " << projection1 << " : " << d_crs1 * geometry::math::r2d << std::endl; + std::cout << "Projection BD-BA " << projection2 << " : " << d_crs2 * geometry::math::r2d << std::endl; #endif + if(projection1 > 0.0 && projection2 > 0.0) + { + return_type XTD = m_radius * geometry::math::abs( asin( sin( d1 / m_radius ) * sin( d_crs1 ) )); + + #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK + std::cout << "Projection ON the segment" << std::endl; + std::cout << "XTD: " << XTD << " d1: " << d1 << " d2: " << d2 << std::endl; +#endif + + // Return shortest distance, projected point on segment sp1-sp2 + return return_type(XTD); + } + else + { +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK + std::cout << "Projection OUTSIDE the segment" << std::endl; +#endif - // Return shortest distance, either to projected point on segment sp1-sp2, or to sp1, or to sp2 - return return_type((std::min)((std::min)(d1, d2), XTD)); + // Return shortest distance, project either on point sp1 or sp2 + return return_type( (std::min)( d1 , d2 ) ); + } } inline return_type radius() const { return m_radius; } diff --git a/boost/geometry/strategies/strategy_transform.hpp b/boost/geometry/strategies/strategy_transform.hpp index 7a1f060169..61a408c617 100644 --- a/boost/geometry/strategies/strategy_transform.hpp +++ b/boost/geometry/strategies/strategy_transform.hpp @@ -23,7 +23,9 @@ #include <boost/geometry/algorithms/convert.hpp> #include <boost/geometry/arithmetic/arithmetic.hpp> #include <boost/geometry/core/access.hpp> +#include <boost/geometry/core/radian_access.hpp> #include <boost/geometry/core/coordinate_dimension.hpp> +#include <boost/geometry/strategies/transform.hpp> #include <boost/geometry/util/math.hpp> #include <boost/geometry/util/select_coordinate_type.hpp> @@ -251,6 +253,23 @@ namespace detail return false; } + template <typename P, typename T> + inline bool cartesian_to_spherical_equatorial3(T x, T y, T z, P& p) + { + assert_dimension<P, 3>(); + + // http://en.wikipedia.org/wiki/List_of_canonical_coordinate_transformations#From_Cartesian_coordinates + T const r = sqrt(x * x + y * y + z * z); + set<2>(p, r); + set_from_radian<0>(p, atan2(y, x)); + if (r > 0.0) + { + set_from_radian<1>(p, asin(z / r)); + return true; + } + return false; + } + } // namespace detail #endif // DOXYGEN_NO_DETAIL @@ -361,6 +380,16 @@ struct from_cartesian_3_to_spherical_polar_3 } }; +template <typename P1, typename P2> +struct from_cartesian_3_to_spherical_equatorial_3 +{ + inline bool apply(P1 const& p1, P2& p2) const + { + assert_dimension<P1, 3>(); + return detail::cartesian_to_spherical_equatorial3(get<0>(p1), get<1>(p1), get<2>(p1), p2); + } +}; + #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS namespace services @@ -454,6 +483,11 @@ struct default_strategy<cartesian_tag, spherical_polar_tag, CoordSys1, CoordSys2 { typedef from_cartesian_3_to_spherical_polar_3<P1, P2> type; }; +template <typename CoordSys1, typename CoordSys2, typename P1, typename P2> +struct default_strategy<cartesian_tag, spherical_equatorial_tag, CoordSys1, CoordSys2, 3, 3, P1, P2> +{ + typedef from_cartesian_3_to_spherical_equatorial_3<P1, P2> type; +}; } // namespace services diff --git a/boost/geometry/util/bare_type.hpp b/boost/geometry/util/bare_type.hpp new file mode 100644 index 0000000000..1b49de6436 --- /dev/null +++ b/boost/geometry/util/bare_type.hpp @@ -0,0 +1,38 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2012 Bruno Lalande, Paris, France. +// Copyright (c) 2012 Mateusz Loskot, London, UK. + +// 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_UTIL_BARE_TYPE_HPP +#define BOOST_GEOMETRY_UTIL_BARE_TYPE_HPP + +#include <boost/type_traits.hpp> + + +namespace boost { namespace geometry +{ + +namespace util +{ + +template <typename T> +struct bare_type +{ + typedef typename boost::remove_const + < + typename boost::remove_pointer<T>::type + >::type type; +}; + + +} // namespace util + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_UTIL_BARE_TYPE_HPP diff --git a/boost/geometry/util/math.hpp b/boost/geometry/util/math.hpp index edfa961b15..95cbdf2ce4 100644 --- a/boost/geometry/util/math.hpp +++ b/boost/geometry/util/math.hpp @@ -44,11 +44,43 @@ struct equals template <typename Type> struct equals<Type, true> { + static inline Type get_max(Type const& a, Type const& b, Type const& c) + { + return (std::max)((std::max)(a, b), c); + } + static inline bool apply(Type const& a, Type const& b) { + if (a == b) + { + return true; + } + // See http://www.parashift.com/c++-faq-lite/newbie.html#faq-29.17, // FUTURE: replace by some boost tool or boost::test::close_at_tolerance - return std::abs(a - b) <= std::numeric_limits<Type>::epsilon() * std::abs(a); + return std::abs(a - b) <= std::numeric_limits<Type>::epsilon() * get_max(std::abs(a), std::abs(b), 1.0); + } +}; + +template <typename Type, bool IsFloatingPoint> +struct smaller +{ + static inline bool apply(Type const& a, Type const& b) + { + return a < b; + } +}; + +template <typename Type> +struct smaller<Type, true> +{ + static inline bool apply(Type const& a, Type const& b) + { + if (equals<Type, true>::apply(a, b)) + { + return false; + } + return a < b; } }; @@ -116,6 +148,28 @@ inline bool equals_with_epsilon(T1 const& a, T2 const& b) >::apply(a, b); } +template <typename T1, typename T2> +inline bool smaller(T1 const& a, T2 const& b) +{ + typedef typename select_most_precise<T1, T2>::type select_type; + return detail::smaller + < + select_type, + boost::is_floating_point<select_type>::type::value + >::apply(a, b); +} + +template <typename T1, typename T2> +inline bool larger(T1 const& a, T2 const& b) +{ + typedef typename select_most_precise<T1, T2>::type select_type; + return detail::smaller + < + select_type, + boost::is_floating_point<select_type>::type::value + >::apply(b, a); +} + double const d2r = geometry::math::pi<double>() / 180.0; |