diff --git a/doc/release_notes.qbk b/doc/release_notes.qbk index 66fef803e..5dd271517 100644 --- a/doc/release_notes.qbk +++ b/doc/release_notes.qbk @@ -27,6 +27,8 @@ * Added rtree const_iterator, begin(), end() and the support for Boost.Range. * The support for C++11 `std::initializer_list` in geometries models. * Disjoint and intersects support the following geometry combinations: multipoint/linestring, multipoint/multilinestring +* Intersection has been implemented for combinations of pointlike and linear geometries +* Added implementation for difference(pointlike, linear) [*Improvements] diff --git a/extensions/test/gis/projections/projections_static.cpp b/extensions/test/gis/projections/projections_static.cpp index 2bfdc4158..1f7432eb0 100644 --- a/extensions/test/gis/projections/projections_static.cpp +++ b/extensions/test/gis/projections/projections_static.cpp @@ -225,7 +225,6 @@ void test_all() test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m", 2); test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m"); - test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m", 18); test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m", 5); test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m +lat_1=60 +lat_2=0"); @@ -241,7 +240,6 @@ void test_all() test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m", 7); test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m +m=0.5 +n=1.785"); - test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m +m=0.5 +n=1.785", 18); test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m", 50); test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m"); @@ -267,7 +265,6 @@ void test_all() test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m"); test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m", 2); - test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m", 18); test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m"); test_forward(amsterdam, utrecht, "+ellps=WGS84 +units=m", 22); diff --git a/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp b/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp index 127e4c3fb..67396123c 100644 --- a/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp +++ b/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp @@ -383,6 +383,7 @@ inline void buffer_point(Point const& point, Collection& collection, std::vector range_out; point_strategy.apply(point, distance_strategy, range_out); collection.add_piece(strategy::buffer::buffered_point, range_out, false); + collection.set_piece_center(point); collection.finish_ring(); } diff --git a/include/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp b/include/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp index c5bb8acc0..8def4d5f4 100644 --- a/include/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp +++ b/include/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp @@ -9,6 +9,7 @@ #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_BUFFER_POLICIES_HPP #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_BUFFER_BUFFER_POLICIES_HPP +#define BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION #include @@ -122,6 +123,10 @@ struct buffer_turn_info intersection_location_type location; +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + robust_point_type rob_pi, rob_pj, rob_qi, rob_qj; +#endif + int count_within; bool within_original; @@ -130,7 +135,9 @@ struct buffer_turn_info int count_on_offsetted; int count_on_helper; +#if ! defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) int count_within_near_offsetted; +#endif bool remove_on_multi; @@ -147,7 +154,9 @@ struct buffer_turn_info , count_in_original(0) , count_on_offsetted(0) , count_on_helper(0) +#if ! defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) , count_within_near_offsetted(0) +#endif , remove_on_multi(false) , count_on_occupied(0) , count_on_multi(0) diff --git a/include/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp b/include/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp index dea50b7dc..bff581392 100644 --- a/include/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp +++ b/include/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp @@ -19,8 +19,10 @@ #include #include +#include #include #include +#include #include @@ -126,6 +128,11 @@ struct buffered_piece_collection typedef geometry::model::ring robust_ring_type; typedef geometry::model::box robust_box_type; + typedef typename default_comparable_distance_result + < + robust_point_type + >::type robust_comparable_radius_type; + typedef typename strategy::side::services::default_strategy < typename cs_tag::type @@ -191,13 +198,13 @@ struct buffered_piece_collection std::vector helper_points; // 4 points for side, 3 points for join - 0 points for flat-end #endif + bool is_convex; bool is_monotonic_increasing[2]; // 0=x, 1=y bool is_monotonic_decreasing[2]; // 0=x, 1=y // Monotonic sections of pieces around points std::vector sections; - // Robust representations // 3: complete ring robust_ring_type robust_ring; @@ -206,6 +213,27 @@ struct buffered_piece_collection robust_box_type robust_offsetted_envelope; std::vector robust_turns; // Used only in insert_rescaled_piece_turns - we might use a map instead + + robust_point_type robust_center; + robust_comparable_radius_type robust_min_comparable_radius; + robust_comparable_radius_type robust_max_comparable_radius; + + piece() + : type(strategy::buffer::piece_type_unknown) + , index(-1) + , left_index(-1) + , right_index(-1) + , last_segment_index(-1) + , offsetted_count(-1) + , is_convex(false) + , robust_min_comparable_radius(0) + , robust_max_comparable_radius(0) + { + is_monotonic_increasing[0] = false; + is_monotonic_increasing[1] = false; + is_monotonic_decreasing[0] = false; + is_monotonic_decreasing[1] = false; + } }; struct robust_original @@ -444,6 +472,7 @@ struct buffered_piece_collection { it->location = inside_buffer; } +#if ! defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) if (it->count_within_near_offsetted > 0) { // Within can have in rare cases a rounding issue. We don't discard this @@ -451,6 +480,7 @@ struct buffered_piece_collection // will never start a new ring from this type of points. it->selectable_start = false; } +#endif } } @@ -544,6 +574,7 @@ struct buffered_piece_collection } } +#if ! defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) // Insert all rescaled turn-points into these rings, to form a // reliable integer-based ring. All turns can be compared (inside) to this // rings to see if they are inside. @@ -582,6 +613,7 @@ struct buffered_piece_collection } BOOST_ASSERT(assert_indices_in_robust_rings()); +#endif } template @@ -606,6 +638,8 @@ struct buffered_piece_collection pc.is_monotonic_decreasing[0] = true; pc.is_monotonic_decreasing[1] = true; + pc.is_convex = geometry::is_convex(pc.robust_ring); + if (pc.offsetted_count < 2) { return; @@ -621,7 +655,6 @@ struct buffered_piece_collection current = next; ++next; } - } void determine_properties() @@ -658,7 +691,31 @@ struct buffered_piece_collection geometry::sectionalize(pc.robust_ring, detail::no_rescale_policy(), pc.sections); - // TODO (next phase) determine min/max radius + // Determine min/max radius + typedef geometry::model::referring_segment + robust_segment_type; + + typename robust_ring_type::const_iterator current = pc.robust_ring.begin(); + typename robust_ring_type::const_iterator next = current + 1; + + for (int i = 1; i < pc.offsetted_count; i++) + { + robust_segment_type s(*current, *next); + robust_comparable_radius_type const d + = geometry::comparable_distance(pc.robust_center, s); + + if (i == 1 || d < pc.robust_min_comparable_radius) + { + pc.robust_min_comparable_radius = d; + } + if (i == 1 || d > pc.robust_max_comparable_radius) + { + pc.robust_max_comparable_radius = d; + } + + current = next; + ++next; + } } inline void prepare_buffered_point_pieces() @@ -772,6 +829,13 @@ struct buffered_piece_collection } } + inline void set_piece_center(point_type const& center) + { + BOOST_ASSERT(! m_pieces.empty()); + geometry::recalculate(m_pieces.back().robust_center, center, + m_robust_policy); + } + inline void finish_ring(bool is_interior = false, bool has_interiors = false) { if (m_first_piece_index == -1) diff --git a/include/boost/geometry/algorithms/detail/buffer/get_piece_turns.hpp b/include/boost/geometry/algorithms/detail/buffer/get_piece_turns.hpp index 6a3daa282..76980074c 100644 --- a/include/boost/geometry/algorithms/detail/buffer/get_piece_turns.hpp +++ b/include/boost/geometry/algorithms/detail/buffer/get_piece_turns.hpp @@ -17,6 +17,7 @@ #include #include #include +#include namespace boost { namespace geometry @@ -28,6 +29,31 @@ namespace detail { namespace buffer { +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) +struct buffer_assign_turn +{ + static bool const include_no_turn = false; + static bool const include_degenerate = false; + static bool const include_opposite = false; + + template + < + typename Info, + typename Point1, + typename Point2, + typename IntersectionInfo + > + static inline void apply(Info& info, Point1 const& p1, Point2 const& p2, IntersectionInfo const& iinfo) + { + info.rob_pi = iinfo.rpi(); + info.rob_pj = iinfo.rpj(); + info.rob_qi = iinfo.rqi(); + info.rob_qj = iinfo.rqj(); + } + +}; +#endif + template < typename Pieces, @@ -204,7 +230,11 @@ class piece_turn_visitor // and iterating in sync with them... typedef detail::overlay::get_turn_info < +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + buffer_assign_turn +#else detail::overlay::assign_null_policy +#endif > turn_policy; turn_policy::apply(*prev1, *it1, *next1, diff --git a/include/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp b/include/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp index 8803efdec..4e386ca96 100644 --- a/include/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp +++ b/include/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp @@ -25,6 +25,11 @@ #include #include #include +#include + +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) +#include +#endif namespace boost { namespace geometry @@ -89,8 +94,10 @@ enum analyse_result analyse_disjoint, analyse_within, analyse_on_original_boundary, - analyse_on_offsetted, - analyse_near_offsetted + analyse_on_offsetted +#if ! defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + , analyse_near_offsetted +#endif }; template @@ -112,6 +119,31 @@ inline analyse_result check_segment(Point const& previous, Point const& current, Turn const& turn, bool from_monotonic) { + +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + typedef geometry::model::referring_segment segment_type; + segment_type const p(turn.rob_pi, turn.rob_pj); + segment_type const q(turn.rob_qi, turn.rob_qj); + segment_type const r(previous, current); + int const side = strategy::side::side_of_intersection::apply(p, q, r, + turn.robust_point); + + if (side == 0) + { + return analyse_on_offsetted; + } + if (side == -1 && from_monotonic) + { + return analyse_within; + } + if (side == 1 && from_monotonic) + { + return analyse_disjoint; + } + return analyse_continue; + +#else + typedef typename strategy::side::services::default_strategy < typename cs_tag::type @@ -156,6 +188,7 @@ inline analyse_result check_segment(Point const& previous, // Not monotonic, on left or right side: continue analysing return analyse_continue; +#endif } @@ -169,14 +202,22 @@ public : typedef typename Turn::robust_point_type point_type; typedef typename geometry::coordinate_type::type coordinate_type; - coordinate_type const point_y = geometry::get<1>(turn.robust_point); - +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + typedef geometry::model::referring_segment segment_type; + segment_type const p(turn.rob_pi, turn.rob_pj); + segment_type const q(turn.rob_qi, turn.rob_qj); +#else typedef strategy::within::winding strategy_type; typename strategy_type::state_type state; strategy_type strategy; boost::ignore_unused(strategy); - +#endif + + BOOST_ASSERT(! piece.sections.empty()); + + coordinate_type const point_y = geometry::get<1>(turn.robust_point); + for (std::size_t s = 0; s < piece.sections.size(); s++) { section_type const& section = piece.sections[s]; @@ -191,6 +232,38 @@ public : point_type const& previous = piece.robust_ring[i - 1]; point_type const& current = piece.robust_ring[i]; +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + + // First check if it is in range - if it is not, the + // expensive side_of_intersection does not need to be + // applied + coordinate_type y1 = geometry::get<1>(previous); + coordinate_type y2 = geometry::get<1>(current); + + if (y1 > y2) + { + std::swap(y1, y2); + } + + if (point_y >= y1 - 1 && point_y <= y2 + 1) + { + segment_type const r(previous, current); + int const side = strategy::side::side_of_intersection::apply(p, q, r, + turn.robust_point); + + // Sections are monotonic in y-dimension + if (side == 1) + { + // Left on segment + return analyse_disjoint; + } + else if (side == 0) + { + // Collinear - TODO: check if really on segment + return analyse_on_offsetted; + } + } +#else analyse_result code = check_segment(previous, current, turn, false); if (code != analyse_continue) { @@ -200,10 +273,15 @@ public : // Get the state (to determine it is within), we don't have // to cover the on-segment case (covered above) strategy.apply(turn.robust_point, previous, current, state); +#endif } } } +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + // It is nowhere outside, and not on segment, so it is within + return analyse_within; +#else int const code = strategy.result(state); if (code == 1) { @@ -216,6 +294,7 @@ public : // Should normally not occur - on-segment is covered return analyse_unknown; +#endif } }; @@ -228,6 +307,48 @@ class analyse_turn_wrt_piece bool is_original, Point const& offsetted) { +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + typedef geometry::model::referring_segment segment_type; + segment_type const p(turn.rob_pi, turn.rob_pj); + segment_type const q(turn.rob_qi, turn.rob_qj); + segment_type const r(s1, s2); + int const side = strategy::side::side_of_intersection::apply(p, q, r, + turn.robust_point); + + if (side == 1) + { + // left of segment + return analyse_disjoint; + } + else if (side == 0) + { + // If is collinear, either on segment or before/after + typedef geometry::model::box box_type; + + box_type box; + geometry::assign_inverse(box); + geometry::expand(box, s1); + geometry::expand(box, s2); + + if (geometry::covered_by(turn.robust_point, box)) + { + // Points on helper-segments are considered as within + // Points on original boundary are processed differently + return is_original + ? analyse_on_original_boundary + : analyse_within; + } + + // It is collinear but not on the segment. Because these + // segments are convex, it is outside + // Unless the offsetted ring is collinear or concave w.r.t. + // helper-segment but that scenario is not yet supported + return analyse_disjoint; + } + + // right of segment + return analyse_continue; +#else typedef typename strategy::side::services::default_strategy < typename cs_tag::type @@ -276,6 +397,7 @@ class analyse_turn_wrt_piece // right of segment return analyse_continue; +#endif } template @@ -492,6 +614,60 @@ class turn_in_piece_visitor return false; } +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + template + static inline int turn_in_convex_piece(Turn const& turn, + Piece const& piece) + { + typedef typename Turn::robust_point_type point_type; + typedef typename Piece::piece_robust_ring_type ring_type; + typedef geometry::model::referring_segment segment; + + segment const p(turn.rob_pi, turn.rob_pj); + segment const q(turn.rob_qi, turn.rob_qj); + + typedef typename boost::range_iterator::type iterator_type; + iterator_type it = boost::begin(piece.robust_ring); + iterator_type end = boost::end(piece.robust_ring); + + // A robust ring is always closed, and always clockwise + for (iterator_type previous = it++; it != end; ++previous, ++it) + { + geometry::equal_to comparator; + if (comparator(*previous, *it)) + { + // Points are the same + continue; + } + + segment r(*previous, *it); + + int const side = strategy::side::side_of_intersection::apply(p, q, r, + turn.robust_point); + + if (side == 1) + { + // IP is left of segment, so it is outside + return -1; // outside + } + else if (side == 0) + { + // IP is collinear with segment. TODO: we should analyze this further + // For now we use the fallback point + if (in_box(*previous, *it, turn.robust_point)) + { + return 0; + } + else + { + return -1; // outside + } + } + } + return 1; // inside + } +#endif + public: @@ -530,12 +706,37 @@ public: } // TODO: mutable_piece to make some on-demand preparations in analyse + Turn& mutable_turn = m_turns[turn.turn_index]; + + if (piece.type == geometry::strategy::buffer::buffered_point) + { + // Optimization for buffer around points: if distance from center + // is not between min/max radius, the result is clear + typedef typename default_comparable_distance_result + < + typename Turn::robust_point_type + >::type distance_type; + + distance_type const cd + = geometry::comparable_distance(piece.robust_center, + turn.robust_point); + + if (cd < piece.robust_min_comparable_radius) + { + mutable_turn.count_within++; + return; + } + if (cd > piece.robust_max_comparable_radius) + { + return; + } + } + analyse_result analyse_code = piece.type == geometry::strategy::buffer::buffered_point ? analyse_turn_wrt_point_piece::apply(turn, piece) : analyse_turn_wrt_piece::apply(turn, piece); - Turn& mutable_turn = m_turns[turn.turn_index]; switch(analyse_code) { case analyse_disjoint : @@ -549,16 +750,32 @@ public: case analyse_within : mutable_turn.count_within++; return; +#if ! defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) case analyse_near_offsetted : mutable_turn.count_within_near_offsetted++; return; +#endif default : break; } - // TODO: this point_in_geometry is a performance-bottleneck here and - // will be replaced completely by extending analyse_piece functionality +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + // We don't know (yet) + int geometry_code = 0; + if (piece.is_convex) + { + geometry_code = turn_in_convex_piece(turn, piece); + } + else + { + + // TODO: this point_in_geometry is a performance-bottleneck here and + // will be replaced completely by extending analyse_piece functionality + geometry_code = detail::within::point_in_geometry(turn.robust_point, piece.robust_ring); + } +#else int geometry_code = detail::within::point_in_geometry(turn.robust_point, piece.robust_ring); +#endif if (geometry_code == 1) { diff --git a/include/boost/geometry/algorithms/detail/disjoint/point_point.hpp b/include/boost/geometry/algorithms/detail/disjoint/point_point.hpp index b1d32bf95..2261efd3a 100644 --- a/include/boost/geometry/algorithms/detail/disjoint/point_point.hpp +++ b/include/boost/geometry/algorithms/detail/disjoint/point_point.hpp @@ -1,12 +1,12 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. -// Copyright (c) 2008-2014 Bruno Lalande, Paris, France. -// Copyright (c) 2009-2014 Mateusz Loskot, London, UK. -// Copyright (c) 2013-2014 Adam Wulkiewicz, Lodz, Poland +// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2008-2015 Bruno Lalande, Paris, France. +// Copyright (c) 2009-2015 Mateusz Loskot, London, UK. +// Copyright (c) 2013-2015 Adam Wulkiewicz, Lodz, Poland -// This file was modified by Oracle on 2013-2014. -// Modifications copyright (c) 2013-2014, Oracle and/or its affiliates. +// This file was modified by Oracle on 2013, 2014, 2015. +// Modifications copyright (c) 2013-2015, Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle @@ -23,11 +23,19 @@ #include +#include + #include -#include +#include +#include +#include +#include #include +#include #include +#include +#include #include @@ -40,12 +48,103 @@ namespace boost { namespace geometry namespace detail { namespace disjoint { + +class point_point_on_spheroid +{ +private: + template + < + typename Units, + typename CoordinateType1, + typename CoordinateType2 + > + static inline bool apply_same_units(CoordinateType1& lon1, + CoordinateType1& lat1, + CoordinateType2& lon2, + CoordinateType2& lat2) + { + math::normalize_spheroidal_coordinates(lon1, lat1); + math::normalize_spheroidal_coordinates(lon2, lat2); + + return ! math::equals(lat1, lat2) || ! math::equals(lon1, lon2); + } + +public: + template + static inline bool apply(Point1 const& p1, Point2 const& p2) + { + typedef typename coordinate_type::type coordinate_type1; + typedef typename coordinate_type::type coordinate_type2; + + typedef typename coordinate_system::type::units units1; + typedef typename coordinate_system::type::units units2; + + if (BOOST_GEOMETRY_CONDITION((boost::is_same::value))) + { + coordinate_type1 lon1 = geometry::get<0>(p1); + coordinate_type1 lat1 = geometry::get<1>(p1); + coordinate_type2 lon2 = geometry::get<0>(p2); + coordinate_type2 lat2 = geometry::get<1>(p2); + + return apply_same_units(lon1, lat1, lon2, lat2); + } + + typedef typename geometry::select_most_precise + < + typename fp_coordinate_type::type, + typename fp_coordinate_type::type + >::type calculation_type; + + calculation_type lon1 = get_as_radian<0>(p1); + calculation_type lat1 = get_as_radian<1>(p1); + + calculation_type lon2 = get_as_radian<0>(p2); + calculation_type lat2 = get_as_radian<1>(p2); + + return apply_same_units(lon1, lat1, lon2, lat2); + } +}; + + +template +< + typename Point1, typename Point2, + std::size_t Dimension, std::size_t DimensionCount, + typename CSTag1 = typename cs_tag::type, + typename CSTag2 = CSTag1 +> +struct point_point + : point_point +{}; + template < typename Point1, typename Point2, std::size_t Dimension, std::size_t DimensionCount > struct point_point + < + Point1, Point2, Dimension, DimensionCount, spherical_equatorial_tag + > : point_point_on_spheroid +{}; + +template +< + typename Point1, typename Point2, + std::size_t Dimension, std::size_t DimensionCount +> +struct point_point + < + Point1, Point2, Dimension, DimensionCount, geographic_tag + > : point_point_on_spheroid +{}; + +template +< + typename Point1, typename Point2, + std::size_t Dimension, std::size_t DimensionCount +> +struct point_point { static inline bool apply(Point1 const& p1, Point2 const& p2) { @@ -61,9 +160,11 @@ struct point_point } }; - template -struct point_point +struct point_point + < + Point1, Point2, DimensionCount, DimensionCount, cartesian_tag + > { static inline bool apply(Point1 const& , Point2 const& ) { diff --git a/include/boost/geometry/algorithms/detail/overlay/clip_linestring.hpp b/include/boost/geometry/algorithms/detail/overlay/clip_linestring.hpp index fc4f65732..b1a25c9f5 100644 --- a/include/boost/geometry/algorithms/detail/overlay/clip_linestring.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/clip_linestring.hpp @@ -1,6 +1,11 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. + +// This file was modified by Oracle on 2015. +// Modifications copyright (c) 2015 Oracle and/or its affiliates. + +// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at @@ -160,10 +165,12 @@ template typename OutputLinestring, typename OutputIterator, typename Range, + typename RobustPolicy, typename Box, typename Strategy > OutputIterator clip_range_with_box(Box const& b, Range const& range, + RobustPolicy const&, OutputIterator out, Strategy const& strategy) { if (boost::begin(range) == boost::end(range)) diff --git a/include/boost/geometry/algorithms/detail/overlay/intersection_insert.hpp b/include/boost/geometry/algorithms/detail/overlay/intersection_insert.hpp index 491a540c1..5b9562674 100644 --- a/include/boost/geometry/algorithms/detail/overlay/intersection_insert.hpp +++ b/include/boost/geometry/algorithms/detail/overlay/intersection_insert.hpp @@ -410,13 +410,13 @@ struct intersection_insert template static inline OutputIterator apply(Linestring const& linestring, Box const& box, - RobustPolicy const& , + RobustPolicy const& robust_policy, OutputIterator out, Strategy const& ) { typedef typename point_type::type point_type; strategy::intersection::liang_barsky lb_strategy; return detail::intersection::clip_range_with_box - (box, linestring, out, lb_strategy); + (box, linestring, robust_policy, out, lb_strategy); } }; @@ -488,7 +488,7 @@ struct intersection_insert template static inline OutputIterator apply(Segment const& segment, Box const& box, - RobustPolicy const& ,// TODO: propagate to clip_range_with_box + RobustPolicy const& robust_policy, OutputIterator out, Strategy const& ) { geometry::segment_view range(segment); @@ -496,7 +496,7 @@ struct intersection_insert typedef typename point_type::type point_type; strategy::intersection::liang_barsky lb_strategy; return detail::intersection::clip_range_with_box - (box, range, out, lb_strategy); + (box, range, robust_policy, out, lb_strategy); } }; diff --git a/include/boost/geometry/algorithms/detail/relate/less.hpp b/include/boost/geometry/algorithms/detail/relate/less.hpp index 3f11d4e87..462fcc35f 100644 --- a/include/boost/geometry/algorithms/detail/relate/less.hpp +++ b/include/boost/geometry/algorithms/detail/relate/less.hpp @@ -14,6 +14,10 @@ #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_LESS_HPP #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RELATE_LESS_HPP +#include +#include +#include + namespace boost { namespace geometry { diff --git a/include/boost/geometry/algorithms/is_convex.hpp b/include/boost/geometry/algorithms/is_convex.hpp new file mode 100644 index 000000000..8feb48db6 --- /dev/null +++ b/include/boost/geometry/algorithms/is_convex.hpp @@ -0,0 +1,160 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_ALGORITHMS_IS_CONVEX_HPP +#define BOOST_GEOMETRY_ALGORITHMS_IS_CONVEX_HPP + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { namespace geometry +{ + + +#ifndef DOXYGEN_NO_DETAIL +namespace detail { namespace is_convex +{ + +struct ring_is_convex +{ + template + static inline bool apply(Ring const& ring) + { + typedef typename geometry::point_type::type point_type; + typedef typename strategy::side::services::default_strategy + < + typename cs_tag::type + >::type side_strategy_type; + + std::size_t n = boost::size(ring); + if (boost::size(ring) < core_detail::closure::minimum_ring_size + < + geometry::closure::value + >::value) + { + // (Too) small rings are considered as non-concave, is convex + return true; + } + + // Walk in clockwise direction, consider ring as closed + // (though closure is not important in this algorithm - any dupped + // point is skipped) + typedef detail::normalized_view view_type; + view_type view(ring); + + typedef geometry::ever_circling_range_iterator it_type; + it_type previous(view); + it_type current(view); + current++; + + std::size_t index = 1; + while (equals::equals_point_point(*current, *previous) && index < n) + { + current++; + index++; + } + + if (index == n) + { + // All points are apparently equal + return true; + } + + it_type next = current; + next++; + while (equals::equals_point_point(*current, *next)) + { + next++; + } + + // We have now three different points on the ring + // Walk through all points, use a counter because of the ever-circling + // iterator + for (std::size_t i = 0; i < n; i++) + { + int const side = side_strategy_type::apply(*previous, *current, *next); + if (side == 1) + { + // Next is on the left side of clockwise ring: + // the piece is not convex + return false; + } + + previous = current; + current = next; + + // Advance next to next different point + // (because there are non-equal points, this loop is not infinite) + next++; + while (equals::equals_point_point(*current, *next)) + { + next++; + } + } + return true; + } +}; + + +}} // namespace detail::is_convex +#endif // DOXYGEN_NO_DETAIL + + +#ifndef DOXYGEN_NO_DISPATCH +namespace dispatch +{ + +template +< + typename Geometry, + typename Tag = typename tag::type +> +struct is_convex : not_implemented +{}; + +template +struct is_convex +{ + static inline bool apply(Box const& ) + { + // Any box is convex (TODO: consider spherical boxes) + return true; + } +}; + +template +struct is_convex : detail::is_convex::ring_is_convex +{}; + + +} // namespace dispatch +#endif // DOXYGEN_NO_DISPATCH + +// TODO: variants + +// TODO: documentation / qbk +template +inline bool is_convex(Geometry const& geometry) +{ + return dispatch::is_convex::apply(geometry); +} + + +}} // namespace boost::geometry + + +#endif // BOOST_GEOMETRY_ALGORITHMS_IS_CONVEX_HPP diff --git a/include/boost/geometry/core/radian_access.hpp b/include/boost/geometry/core/radian_access.hpp index bac77d7bc..374c2954a 100644 --- a/include/boost/geometry/core/radian_access.hpp +++ b/include/boost/geometry/core/radian_access.hpp @@ -1,8 +1,13 @@ // 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. +// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2008-2015 Bruno Lalande, Paris, France. +// Copyright (c) 2009-2015 Mateusz Loskot, London, UK. + +// This file was modified by Oracle on 2015. +// Modifications copyright (c) 2015, Oracle and/or its affiliates. + +// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -47,7 +52,8 @@ struct degree_radian_converter return boost::numeric_cast < coordinate_type - >(geometry::get(geometry) * geometry::math::d2r); + >(geometry::get(geometry) + * math::d2r()); } static inline void set(Geometry& geometry, coordinate_type const& radians) @@ -55,7 +61,7 @@ struct degree_radian_converter geometry::set(geometry, boost::numeric_cast < coordinate_type - >(radians * geometry::math::r2d)); + >(radians * math::r2d())); } }; diff --git a/include/boost/geometry/extensions/gis/projections/proj/aeqd.hpp b/include/boost/geometry/extensions/gis/projections/proj/aeqd.hpp index 35c41676a..125294235 100644 --- a/include/boost/geometry/extensions/gis/projections/proj/aeqd.hpp +++ b/include/boost/geometry/extensions/gis/projections/proj/aeqd.hpp @@ -419,12 +419,14 @@ namespace boost { namespace geometry { namespace projections public : virtual projection* create_new(const Parameters& par) const { - if (! par.es) - return new base_v_fi, Geographic, Cartesian, Parameters>(par); - else if (pj_param(par.params, "bguam").i) + bool const guam = pj_param(par.params, "bguam").i; + + if (par.es && ! guam) + return new base_v_fi, Geographic, Cartesian, Parameters>(par); + else if (par.es && guam) return new base_v_fi, Geographic, Cartesian, Parameters>(par); else - return new base_v_fi, Geographic, Cartesian, Parameters>(par); + return new base_v_fi, Geographic, Cartesian, Parameters>(par); } }; diff --git a/include/boost/geometry/extensions/gis/projections/proj/gn_sinu.hpp b/include/boost/geometry/extensions/gis/projections/proj/gn_sinu.hpp index ac0d42d4e..ee7fe43ca 100644 --- a/include/boost/geometry/extensions/gis/projections/proj/gn_sinu.hpp +++ b/include/boost/geometry/extensions/gis/projections/proj/gn_sinu.hpp @@ -213,9 +213,9 @@ namespace boost { namespace geometry { namespace projections \image html ex_gn_sinu.gif */ template - struct gn_sinu_ellipsoid : public detail::gn_sinu::base_gn_sinu_ellipsoid + struct gn_sinu_spheroid : public detail::gn_sinu::base_gn_sinu_spheroid { - inline gn_sinu_ellipsoid(const Parameters& par) : detail::gn_sinu::base_gn_sinu_ellipsoid(par) + inline gn_sinu_spheroid(const Parameters& par) : detail::gn_sinu::base_gn_sinu_spheroid(par) { detail::gn_sinu::setup_gn_sinu(this->m_par, this->m_proj_parm); } @@ -243,70 +243,6 @@ namespace boost { namespace geometry { namespace projections } }; - /*! - \brief Eckert VI projection - \ingroup projections - \tparam Geographic latlong point type - \tparam Cartesian xy point type - \tparam Parameters parameter type - \par Projection characteristics - - Pseudocylindrical - - Spheroid - \par Example - \image html ex_eck6.gif - */ - template - struct eck6_ellipsoid : public detail::gn_sinu::base_gn_sinu_ellipsoid - { - inline eck6_ellipsoid(const Parameters& par) : detail::gn_sinu::base_gn_sinu_ellipsoid(par) - { - detail::gn_sinu::setup_eck6(this->m_par, this->m_proj_parm); - } - }; - - /*! - \brief McBryde-Thomas Flat-Polar Sinusoidal projection - \ingroup projections - \tparam Geographic latlong point type - \tparam Cartesian xy point type - \tparam Parameters parameter type - \par Projection characteristics - - Pseudocylindrical - - Spheroid - \par Example - \image html ex_mbtfps.gif - */ - template - struct mbtfps_ellipsoid : public detail::gn_sinu::base_gn_sinu_ellipsoid - { - inline mbtfps_ellipsoid(const Parameters& par) : detail::gn_sinu::base_gn_sinu_ellipsoid(par) - { - detail::gn_sinu::setup_mbtfps(this->m_par, this->m_proj_parm); - } - }; - - /*! - \brief General Sinusoidal Series projection - \ingroup projections - \tparam Geographic latlong point type - \tparam Cartesian xy point type - \tparam Parameters parameter type - \par Projection characteristics - - Pseudocylindrical - - Spheroid - - m= n= - \par Example - \image html ex_gn_sinu.gif - */ - template - struct gn_sinu_spheroid : public detail::gn_sinu::base_gn_sinu_spheroid - { - inline gn_sinu_spheroid(const Parameters& par) : detail::gn_sinu::base_gn_sinu_spheroid(par) - { - detail::gn_sinu::setup_gn_sinu(this->m_par, this->m_proj_parm); - } - }; - /*! \brief Sinusoidal (Sanson-Flamsteed) projection \ingroup projections diff --git a/include/boost/geometry/extensions/gis/projections/proj/isea.hpp b/include/boost/geometry/extensions/gis/projections/proj/isea.hpp index a2a3707a1..a38e3028f 100644 --- a/include/boost/geometry/extensions/gis/projections/proj/isea.hpp +++ b/include/boost/geometry/extensions/gis/projections/proj/isea.hpp @@ -56,11 +56,11 @@ namespace boost { namespace geometry { namespace projections static const double DEG144 = 2.51327412287183459075; static const double DEG36 = 0.62831853071795864768; static const double DEG108 = 1.88495559215387594306; - static const double DEG180 = M_PI; + static const double DEG180 = boost::math::constants::pi(); static const double ISEA_SCALE = 0.8301572857837594396028083; static const double V_LAT = 0.46364760899944494524; - static const double RAD2DEG = (180.0/M_PI); - static const double DEG2RAD = (M_PI/180.0); + static const double RAD2DEG = (180.0/boost::math::constants::pi()); + static const double DEG2RAD = (boost::math::constants::pi()/180.0); static const double E_RAD = 0.91843818702186776133; static const double F_RAD = 0.18871053072122403508; static const double TABLE_G = 0.6615845383; @@ -116,7 +116,7 @@ namespace boost { namespace geometry { namespace projections int ix, iy, iz, s; struct hex h; - x = x / cos(30 * M_PI / 180.0); /* rotated X coord */ + x = x / cos(30 * boost::math::constants::pi() / 180.0); /* rotated X coord */ y = y - x / 2.0; /* adjustment for rotated X */ /* adjust for actual hexwidth */ @@ -215,7 +215,7 @@ namespace boost { namespace geometry { namespace projections }; - /* sqrt(5)/M_PI */ + /* sqrt(5)/boost::math::constants::pi() */ /* 26.565051177 degrees */ @@ -409,7 +409,7 @@ namespace boost { namespace geometry { namespace projections /* TODO I don't know why we do this. It's not in snyder */ /* maybe because we should have picked a better vertex */ if (Az < 0.0) { - Az += 2.0 * M_PI; + Az += 2.0 * boost::math::constants::pi(); } Az_adjust_multiples = 0; @@ -447,7 +447,7 @@ namespace boost { namespace geometry { namespace projections H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G)); /* eq 7 */ - /* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */ + /* Ag = (Az + G + H - DEG180) * boost::math::constants::pi() * R * R / DEG180; */ Ag = Az + G + H - DEG180; /* eq 8 */ @@ -523,11 +523,11 @@ namespace boost { namespace geometry { namespace projections /* normalize longitude */ /* TODO can we just do a modulus ? */ - lambdap = fmod(lambdap, 2 * M_PI); - while (lambdap > M_PI) - lambdap -= 2 * M_PI; - while (lambdap < -M_PI) - lambdap += 2 * M_PI; + lambdap = fmod(lambdap, 2 * boost::math::constants::pi()); + while (lambdap > boost::math::constants::pi()) + lambdap -= 2 * boost::math::constants::pi(); + while (lambdap < -boost::math::constants::pi()) + lambdap += 2 * boost::math::constants::pi(); phip = asin(sin_phip); @@ -543,19 +543,19 @@ namespace boost { namespace geometry { namespace projections { struct isea_geo npt; - np->lon += M_PI; + np->lon += boost::math::constants::pi(); npt = snyder_ctran(np, pt); - np->lon -= M_PI; + np->lon -= boost::math::constants::pi(); - npt.lon -= (M_PI - lon0 + np->lon); + npt.lon -= (boost::math::constants::pi() - lon0 + np->lon); - npt.lon += M_PI; + npt.lon += boost::math::constants::pi(); /* normalize longitude */ - npt.lon = fmod(npt.lon, 2 * M_PI); - while (npt.lon > M_PI) - npt.lon -= 2 * M_PI; - while (npt.lon < -M_PI) - npt.lon += 2 * M_PI; + npt.lon = fmod(npt.lon, 2 * boost::math::constants::pi()); + while (npt.lon > boost::math::constants::pi()) + npt.lon -= 2 * boost::math::constants::pi(); + while (npt.lon < -boost::math::constants::pi()) + npt.lon += 2 * boost::math::constants::pi(); return npt; } @@ -601,7 +601,7 @@ namespace boost { namespace geometry { namespace projections { if (!g) return 0; - g->o_lat = M_PI / 2.0; + g->o_lat = boost::math::constants::pi() / 2.0; g->o_lon = 0.0; g->o_az = 0; return 1; @@ -637,9 +637,9 @@ namespace boost { namespace geometry { namespace projections double x, y; - rad = -degrees * M_PI / 180.0; - while (rad >= 2.0 * M_PI) rad -= 2.0 * M_PI; - while (rad <= -2.0 * M_PI) rad += 2.0 * M_PI; + rad = -degrees * boost::math::constants::pi() / 180.0; + while (rad >= 2.0 * boost::math::constants::pi()) rad -= 2.0 * boost::math::constants::pi(); + while (rad <= -2.0 * boost::math::constants::pi()) rad += 2.0 * boost::math::constants::pi(); x = pt->x * cos(rad) + pt->y * sin(rad); y = -pt->x * sin(rad) + pt->y * cos(rad); @@ -676,7 +676,7 @@ namespace boost { namespace geometry { namespace projections isea_rotate(pt, downtri ? 240.0 : 60.0); if (downtri) { pt->x += 0.5; - /* pt->y += cos(30.0 * M_PI / 180.0); */ + /* pt->y += cos(30.0 * boost::math::constants::pi() / 180.0); */ pt->y += .86602540378443864672; } return quad; @@ -697,7 +697,7 @@ namespace boost { namespace geometry { namespace projections sidelength = (pow(2.0, g->resolution) + 1.0) / 2.0; /* apex to base is cos(30deg) */ - hexwidth = cos(M_PI / 6.0) / sidelength; + hexwidth = cos(boost::math::constants::pi() / 6.0) / sidelength; /* TODO I think sidelength is always x.5, so * (int)sidelength * 2 + 1 might be just as good diff --git a/include/boost/geometry/extensions/gis/projections/proj/ob_tran.hpp b/include/boost/geometry/extensions/gis/projections/proj/ob_tran.hpp index cd4928893..fd99c475f 100644 --- a/include/boost/geometry/extensions/gis/projections/proj/ob_tran.hpp +++ b/include/boost/geometry/extensions/gis/projections/proj/ob_tran.hpp @@ -292,6 +292,7 @@ namespace boost { namespace geometry { namespace projections detail::ob_tran::par_ob_tran proj_parm; Parameters p = par; double phip = setup_ob_tran(p, proj_parm, false); + if (fabs(phip) > detail::ob_tran::TOL) return new base_v_fi, Geographic, Cartesian, Parameters>(par); else diff --git a/include/boost/geometry/extensions/gis/projections/proj/stere.hpp b/include/boost/geometry/extensions/gis/projections/proj/stere.hpp index 8f258a699..ee23389b9 100644 --- a/include/boost/geometry/extensions/gis/projections/proj/stere.hpp +++ b/include/boost/geometry/extensions/gis/projections/proj/stere.hpp @@ -360,29 +360,6 @@ namespace boost { namespace geometry { namespace projections } }; - /*! - \brief Universal Polar Stereographic projection - \ingroup projections - \tparam Geographic latlong point type - \tparam Cartesian xy point type - \tparam Parameters parameter type - \par Projection characteristics - - Azimuthal - - Spheroid - - Ellipsoid - - south - \par Example - \image html ex_ups.gif - */ - template - struct ups_ellipsoid : public detail::stere::base_stere_ellipsoid - { - inline ups_ellipsoid(const Parameters& par) : detail::stere::base_stere_ellipsoid(par) - { - detail::stere::setup_ups(this->m_par, this->m_proj_parm); - } - }; - /*! \brief Stereographic projection \ingroup projections @@ -406,6 +383,29 @@ namespace boost { namespace geometry { namespace projections } }; + /*! + \brief Universal Polar Stereographic projection + \ingroup projections + \tparam Geographic latlong point type + \tparam Cartesian xy point type + \tparam Parameters parameter type + \par Projection characteristics + - Azimuthal + - Spheroid + - Ellipsoid + - south + \par Example + \image html ex_ups.gif + */ + template + struct ups_ellipsoid : public detail::stere::base_stere_ellipsoid + { + inline ups_ellipsoid(const Parameters& par) : detail::stere::base_stere_ellipsoid(par) + { + detail::stere::setup_ups(this->m_par, this->m_proj_parm); + } + }; + /*! \brief Universal Polar Stereographic projection \ingroup projections @@ -453,7 +453,10 @@ namespace boost { namespace geometry { namespace projections public : virtual projection* create_new(const Parameters& par) const { - return new base_v_fi, Geographic, Cartesian, Parameters>(par); + if (par.es) + return new base_v_fi, Geographic, Cartesian, Parameters>(par); + else + return new base_v_fi, Geographic, Cartesian, Parameters>(par); } }; diff --git a/include/boost/geometry/extensions/gis/projections/proj/tmerc.hpp b/include/boost/geometry/extensions/gis/projections/proj/tmerc.hpp index daff56690..c68209f54 100644 --- a/include/boost/geometry/extensions/gis/projections/proj/tmerc.hpp +++ b/include/boost/geometry/extensions/gis/projections/proj/tmerc.hpp @@ -284,28 +284,6 @@ namespace boost { namespace geometry { namespace projections } }; - /*! - \brief Universal Transverse Mercator (UTM) projection - \ingroup projections - \tparam Geographic latlong point type - \tparam Cartesian xy point type - \tparam Parameters parameter type - \par Projection characteristics - - Cylindrical - - Spheroid - - zone= south - \par Example - \image html ex_utm.gif - */ - template - struct utm_ellipsoid : public detail::tmerc::base_tmerc_ellipsoid - { - inline utm_ellipsoid(const Parameters& par) : detail::tmerc::base_tmerc_ellipsoid(par) - { - detail::tmerc::setup_utm(this->m_par, this->m_proj_parm); - } - }; - /*! \brief Transverse Mercator projection \ingroup projections @@ -328,6 +306,28 @@ namespace boost { namespace geometry { namespace projections } }; + /*! + \brief Universal Transverse Mercator (UTM) projection + \ingroup projections + \tparam Geographic latlong point type + \tparam Cartesian xy point type + \tparam Parameters parameter type + \par Projection characteristics + - Cylindrical + - Spheroid + - zone= south + \par Example + \image html ex_utm.gif + */ + template + struct utm_ellipsoid : public detail::tmerc::base_tmerc_ellipsoid + { + inline utm_ellipsoid(const Parameters& par) : detail::tmerc::base_tmerc_ellipsoid(par) + { + detail::tmerc::setup_utm(this->m_par, this->m_proj_parm); + } + }; + /*! \brief Universal Transverse Mercator (UTM) projection \ingroup projections @@ -374,7 +374,10 @@ namespace boost { namespace geometry { namespace projections public : virtual projection* create_new(const Parameters& par) const { - return new base_v_fi, Geographic, Cartesian, Parameters>(par); + if (par.es) + return new base_v_fi, Geographic, Cartesian, Parameters>(par); + else + return new base_v_fi, Geographic, Cartesian, Parameters>(par); } }; diff --git a/include/boost/geometry/extensions/gis/projections/project_inverse_transformer.hpp b/include/boost/geometry/extensions/gis/projections/project_inverse_transformer.hpp index 4386c8eea..136d51133 100644 --- a/include/boost/geometry/extensions/gis/projections/project_inverse_transformer.hpp +++ b/include/boost/geometry/extensions/gis/projections/project_inverse_transformer.hpp @@ -60,8 +60,8 @@ struct project_inverse_transformer { // Latlong (LL -> XY) will be projected, rest will be copied. // So first copy third or higher dimensions - geometry::detail::convert::point_to_point::value> ::copy(p1, p2); + geometry::detail::conversion::point_to_point::value> ::apply(p1, p2); return m_prj->inverse(p1, p2); } diff --git a/include/boost/geometry/strategies/buffer.hpp b/include/boost/geometry/strategies/buffer.hpp index 7dbe03b4a..be3daefec 100644 --- a/include/boost/geometry/strategies/buffer.hpp +++ b/include/boost/geometry/strategies/buffer.hpp @@ -66,7 +66,8 @@ enum piece_type buffered_round_end, buffered_flat_end, buffered_point, - buffered_concave // always on the inside + buffered_concave, // always on the inside + piece_type_unknown }; diff --git a/include/boost/geometry/strategies/cartesian/side_of_intersection.hpp b/include/boost/geometry/strategies/cartesian/side_of_intersection.hpp index 39487676c..d95510223 100644 --- a/include/boost/geometry/strategies/cartesian/side_of_intersection.hpp +++ b/include/boost/geometry/strategies/cartesian/side_of_intersection.hpp @@ -13,8 +13,13 @@ #include #include #include +#include +#include #include +#include +#include +#include namespace boost { namespace geometry { @@ -29,15 +34,74 @@ namespace strategy { namespace side // It can be used for either integer (rescaled) points, and also for FP class side_of_intersection { +private : + template + static inline + int sign_of_product(T const& a, U const& b) + { + return a == 0 || b == 0 ? 0 + : a > 0 && b > 0 ? 1 + : a < 0 && b < 0 ? 1 + : -1; + } + + template + static inline + int sign_of_compare(T const& a, T const& b, T const& c, T const& d) + { + // Both a*b and c*d are positive + // We have to judge if a*b > c*d + + using namespace boost::multiprecision; + cpp_int const lab = cpp_int(a) * cpp_int(b); + cpp_int const lcd = cpp_int(c) * cpp_int(d); + + return lab > lcd ? 1 + : lab < lcd ? -1 + : 0 + ; + } + + template + static inline + int sign_of_addition_of_two_products(T const& a, T const& b, T const& c, T const& d) + { + // sign of a*b+c*d, 1 if positive, -1 if negative, else 0 + int const ab = sign_of_product(a, b); + int const cd = sign_of_product(c, d); + if (ab == 0) + { + return cd; + } + if (cd == 0) + { + return ab; + } + + if (ab == cd) + { + // Both positive or both negative + return ab; + } + + // One is positive, one is negative, both are non zero + // If ab is positive, we have to judge if a*b > -c*d (then 1 because sum is positive) + // If ab is negative, we have to judge if c*d > -a*b (idem) + return ab == 1 + ? sign_of_compare(a, b, -c, d) + : sign_of_compare(c, d, -a, b); + } + + public : // Calculates the side of the intersection-point (if any) of // of segment a//b w.r.t. segment c // This is calculated without (re)calculating the IP itself again and fully // based on integer mathematics - template + template static inline T side_value(Segment const& a, Segment const& b, - Segment const& c) + Segment const& c, Point const& fallback_point) { // The first point of the three segments is reused several times T const ax = get<0, 0>(a); @@ -67,9 +131,15 @@ public : if (d == zero) { // There is no IP of a//b, they are collinear or parallel - // We don't have to divide but we can already conclude the side-value - // is meaningless and the resulting determinant will be 0 - return zero; + // Assuming they intersect (this method should be called for + // segments known to intersect), they are collinear and overlap. + // They have one or two intersection points - we don't know and + // have to rely on the fallback intersection point + + Point c1, c2; + geometry::detail::assign_point_from_index<0>(c, c1); + geometry::detail::assign_point_from_index<1>(c, c2); + return side_by_triangle<>::apply(c1, c2, fallback_point); } // Cramer's rule: da (see cart_intersect.hpp) @@ -82,7 +152,9 @@ public : // IP is at (ax + (da/d) * dx_a, ay + (da/d) * dy_a) // Side of IP is w.r.t. c is: determinant(dx_c, dy_c, ipx-cx, ipy-cy) // We replace ipx by expression above and multiply each term by d - T const result = geometry::detail::determinant + +#ifdef BOOST_GEOMETRY_SIDE_OF_INTERSECTION_DEBUG + T const result1 = geometry::detail::determinant ( dx_c * d, dy_c * d, d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da @@ -93,15 +165,51 @@ public : // Therefore, the sign is always the same as that result, and the // resulting side (left,right,collinear) is the same + // The first row we divide again by d because of determinant multiply rule + T const result2 = d * geometry::detail::determinant + ( + dx_c, dy_c, + d * (ax - cx) + dx_a * da, d * (ay - cy) + dy_a * da + ); + // Write out: + T const result3 = d * (dx_c * (d * (ay - cy) + dy_a * da) + - dy_c * (d * (ax - cx) + dx_a * da)); + // Write out in braces: + T const result4 = d * (dx_c * d * (ay - cy) + dx_c * dy_a * da + - dy_c * d * (ax - cx) - dy_c * dx_a * da); + // Write in terms of d * XX + da * YY + T const result5 = d * (d * (dx_c * (ay - cy) - dy_c * (ax - cx)) + + da * (dx_c * dy_a - dy_c * dx_a)); + + //return result; +#endif + + // We consider the results separately + // (in the end we only have to return the side-value 1,0 or -1) + + // To avoid multiplications we judge the product (easy, avoids *d) + // and the sign of p*q+r*s (more elaborate) + T const result = sign_of_product + ( + d, + sign_of_addition_of_two_products + ( + d, dx_c * (ay - cy) - dy_c * (ax - cx), + da, dx_c * dy_a - dy_c * dx_a + ) + ); return result; + } - template - static inline int apply(Segment const& a, Segment const& b, Segment const& c) + template + static inline int apply(Segment const& a, Segment const& b, + Segment const& c, + Point const& fallback_point) { typedef typename geometry::coordinate_type::type coordinate_type; - coordinate_type const s = side_value(a, b, c); + coordinate_type const s = side_value(a, b, c, fallback_point); coordinate_type const zero = coordinate_type(); return math::equals(s, zero) ? 0 : s > zero ? 1 diff --git a/include/boost/geometry/strategies/spherical/distance_cross_track.hpp b/include/boost/geometry/strategies/spherical/distance_cross_track.hpp index 486c7effb..31993e76c 100644 --- a/include/boost/geometry/strategies/spherical/distance_cross_track.hpp +++ b/include/boost/geometry/strategies/spherical/distance_cross_track.hpp @@ -372,11 +372,16 @@ public : typedef typename return_type::type return_type; #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 << "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; + 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 // http://williams.best.vwh.net/avform.htm#XTE diff --git a/include/boost/geometry/strategies/strategy_transform.hpp b/include/boost/geometry/strategies/strategy_transform.hpp index 9d3b1d910..99fcb720c 100644 --- a/include/boost/geometry/strategies/strategy_transform.hpp +++ b/include/boost/geometry/strategies/strategy_transform.hpp @@ -1,8 +1,13 @@ // 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. +// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2008-2015 Bruno Lalande, Paris, France. +// Copyright (c) 2009-2015 Mateusz Loskot, London, UK. + +// This file was modified by Oracle on 2015. +// Modifications copyright (c) 2015 Oracle and/or its affiliates. + +// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -28,6 +33,7 @@ #include #include +#include #include namespace boost { namespace geometry @@ -130,7 +136,15 @@ struct degree_radian_vv assert_dimension(); assert_dimension(); - detail::transform_coordinates::transform(p1, p2, math::d2r); + typedef typename promote_floating_point + < + typename select_coordinate_type::type + >::type calculation_type; + + detail::transform_coordinates + < + P1, P2, 0, 2, F + >::transform(p1, p2, math::d2r()); return true; } }; @@ -143,7 +157,16 @@ struct degree_radian_vv_3 assert_dimension(); assert_dimension(); - detail::transform_coordinates::transform(p1, p2, math::d2r); + typedef typename promote_floating_point + < + typename select_coordinate_type::type + >::type calculation_type; + + detail::transform_coordinates + < + P1, P2, 0, 2, F + >::transform(p1, p2, math::d2r()); + // Copy height or other third dimension set<2>(p2, get<2>(p1)); return true; diff --git a/include/boost/geometry/strategies/transform/matrix_transformers.hpp b/include/boost/geometry/strategies/transform/matrix_transformers.hpp index 27a3a2ae8..699b91b3a 100644 --- a/include/boost/geometry/strategies/transform/matrix_transformers.hpp +++ b/include/boost/geometry/strategies/transform/matrix_transformers.hpp @@ -1,8 +1,13 @@ // 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. +// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2008-2015 Bruno Lalande, Paris, France. +// Copyright (c) 2009-2015 Mateusz Loskot, London, UK. + +// This file was modified by Oracle on 2015. +// Modifications copyright (c) 2015 Oracle and/or its affiliates. + +// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands. @@ -40,6 +45,7 @@ #include #include #include +#include #include #include @@ -340,7 +346,8 @@ struct as_radian template static inline T get(T const& value) { - return value * math::d2r; + typedef typename promote_floating_point::type promoted_type; + return value * math::d2r(); } }; diff --git a/include/boost/geometry/util/math.hpp b/include/boost/geometry/util/math.hpp index 8361ee76f..7b18b879c 100644 --- a/include/boost/geometry/util/math.hpp +++ b/include/boost/geometry/util/math.hpp @@ -493,9 +493,20 @@ inline bool larger(T1 const& a, T2 const& b) } +template +inline T d2r() +{ + static T const conversion_coefficient = geometry::math::pi() / T(180.0); + return conversion_coefficient; +} + +template +inline T r2d() +{ + static T const conversion_coefficient = T(180.0) / geometry::math::pi(); + return conversion_coefficient; +} -double const d2r = geometry::math::pi() / 180.0; -double const r2d = 1.0 / d2r; /*! \brief Calculates the haversine of an angle diff --git a/include/boost/geometry/util/normalize_spheroidal_coordinates.hpp b/include/boost/geometry/util/normalize_spheroidal_coordinates.hpp new file mode 100644 index 000000000..2379b2d22 --- /dev/null +++ b/include/boost/geometry/util/normalize_spheroidal_coordinates.hpp @@ -0,0 +1,254 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2015, Oracle and/or its affiliates. + +// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle + +// Licensed under the Boost Software License version 1.0. +// http://www.boost.org/users/license.html + +#ifndef BOOST_GEOMETRY_UTIL_NORMALIZE_SPHEROIDAL_COORDINATES_HPP +#define BOOST_GEOMETRY_UTIL_NORMALIZE_SPHEROIDAL_COORDINATES_HPP + +#include + +#include +#include + + +namespace boost { namespace geometry +{ + +namespace math +{ + +#ifndef DOXYGEN_NO_DETAIL +namespace detail +{ + + +template +struct constants_on_spheroid +{ + static inline CoordinateType period() + { + return math::two_pi(); + } + + static inline CoordinateType half_period() + { + return math::pi(); + } + + static inline CoordinateType min_longitude() + { + static CoordinateType const minus_pi = -math::pi(); + return minus_pi; + } + + static inline CoordinateType max_longitude() + { + return math::pi(); + } + + static inline CoordinateType min_latitude() + { + static CoordinateType const minus_half_pi + = -math::half_pi(); + return minus_half_pi; + } + + static inline CoordinateType max_latitude() + { + return math::half_pi(); + } +}; + +template +struct constants_on_spheroid +{ + static inline CoordinateType period() + { + return CoordinateType(360.0); + } + + static inline CoordinateType half_period() + { + return CoordinateType(180.0); + } + + static inline CoordinateType min_longitude() + { + return CoordinateType(-180.0); + } + + static inline CoordinateType max_longitude() + { + return CoordinateType(180.0); + } + + static inline CoordinateType min_latitude() + { + return CoordinateType(-90.0); + } + + static inline CoordinateType max_latitude() + { + return CoordinateType(90.0); + } +}; + + +template +class normalize_spheroidal_coordinates +{ + typedef constants_on_spheroid constants; + +protected: + static inline CoordinateType normalize_up(CoordinateType const& value) + { + return + math::mod(value + constants::half_period(), constants::period()) + - constants::half_period(); + } + + static inline CoordinateType normalize_down(CoordinateType const& value) + { + return + math::mod(value - constants::half_period(), constants::period()) + + constants::half_period(); + } + +public: + static inline void apply(CoordinateType& longitude, + CoordinateType& latitude, + bool normalize_poles = true) + { + + // first normalize latitude + if (math::larger(latitude, constants::half_period())) + { + latitude = normalize_up(latitude); + } + else if (math::smaller(latitude, -constants::half_period())) + { + latitude = normalize_down(latitude); + } + + // fix latitude range + if (latitude < constants::min_latitude()) + { + latitude = -constants::half_period() - latitude; + longitude -= constants::half_period(); + } + else if (latitude > constants::max_latitude()) + { + latitude = constants::half_period() - latitude; + longitude -= constants::half_period(); + } + + // now normalize longitude + if (math::equals(longitude, -constants::half_period()) + || math::equals(longitude, constants::half_period())) + { + longitude = constants::half_period(); + } + else if (longitude > constants::half_period()) + { + longitude = normalize_up(longitude); + if (math::equals(longitude, -constants::half_period())) + { + longitude = constants::half_period(); + } + } + else if (longitude < -constants::half_period()) + { + longitude = normalize_down(longitude); + } + + // finally normalize poles + if (normalize_poles) + { + if (math::equals(math::abs(latitude), constants::max_latitude())) + { + // for the north and south pole we set the longitude to 0 + // (works for both radians and degrees) + longitude = CoordinateType(0); + } + } + + BOOST_ASSERT(! math::larger(constants::min_latitude(), latitude)); + BOOST_ASSERT(! math::larger(latitude, constants::max_latitude())); + + BOOST_ASSERT(math::smaller(constants::min_longitude(), longitude)); + BOOST_ASSERT(! math::larger(longitude, constants::max_longitude())); + } +}; + + +} // namespace detail +#endif // DOXYGEN_NO_DETAIL + + +template +inline CoordinateType min_longitude() +{ + return detail::constants_on_spheroid + < + CoordinateType, Units + >::min_longitude(); +} + +template +inline CoordinateType max_longitude() +{ + return detail::constants_on_spheroid + < + CoordinateType, Units + >::max_longitude(); +} + +template +inline CoordinateType min_latitude() +{ + return detail::constants_on_spheroid + < + CoordinateType, Units + >::min_latitude(); +} + +template +inline CoordinateType max_latitude() +{ + return detail::constants_on_spheroid + < + CoordinateType, Units + >::max_latitude(); +} + + +/*! +\brief Short utility to normalize the coordinates on a spheroid +\tparam Units The units of the coordindate system in the spheroid +\tparam CoordinateType The type of the coordinates +\param longitude Longitude +\param latitude Latitude +\ingroup utility +*/ +template +inline void normalize_spheroidal_coordinates(CoordinateType& longitude, + CoordinateType& latitude) +{ + detail::normalize_spheroidal_coordinates + < + Units, CoordinateType + >::apply(longitude, latitude); +} + + +} // namespace math + + +}} // namespace boost::geometry + +#endif // BOOST_GEOMETRY_UTIL_NORMALIZE_SPHEROIDAL_COORDINATES_HPP diff --git a/test/algorithms/Jamfile.v2 b/test/algorithms/Jamfile.v2 index 89814a628..18d0a56a2 100644 --- a/test/algorithms/Jamfile.v2 +++ b/test/algorithms/Jamfile.v2 @@ -30,6 +30,7 @@ test-suite boost-geometry-algorithms [ run expand.cpp ] [ run expand_on_spheroid.cpp ] [ run for_each.cpp ] + [ run is_convex.cpp ] [ run is_simple.cpp ] [ run is_valid.cpp ] [ run is_valid_failure.cpp ] diff --git a/test/algorithms/buffer/linestring_buffer.cpp b/test/algorithms/buffer/linestring_buffer.cpp index 5179eff9b..1ce527420 100644 --- a/test/algorithms/buffer/linestring_buffer.cpp +++ b/test/algorithms/buffer/linestring_buffer.cpp @@ -203,9 +203,10 @@ void test_all() double const d15 = 1.5; test_one("mysql_report_2015_03_02c_asym1", mysql_report_2015_03_02c, join_round(7), end_round(7), 39.714, d10, d15); test_one("mysql_report_2015_03_02c_asym2", mysql_report_2015_03_02c, join_round(7), end_round(7), 46.116, d15, d10); -#if defined(BOOST_GEOMETRY_BUFFER_INCLUDE_FAILING_TESTS) + +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) double const d100 = 10; - test_one("mysql_report_2015_04_01", mysql_report_2015_04_01, join_round(32), end_round(32), 1.0/*NON ZERO*/, d100); + test_one("mysql_report_2015_04_01", mysql_report_2015_04_01, join_round(32), end_round(32), 632.234, d100); #endif } diff --git a/test/algorithms/buffer/multi_linestring_buffer.cpp b/test/algorithms/buffer/multi_linestring_buffer.cpp index 1f09ceeb8..f4ff52be0 100644 --- a/test/algorithms/buffer/multi_linestring_buffer.cpp +++ b/test/algorithms/buffer/multi_linestring_buffer.cpp @@ -100,7 +100,9 @@ void test_all() test_one("mikado4_small", mikado4, join_round32, end_flat, 1930.785, 10.0); #ifdef BOOST_GEOMETRY_BUFFER_INCLUDE_FAILING_TESTS - test_one("mysql_15_04_10", mysql_15_04_10, join_round32, end_round32, 29151950703.779/*something big*/, 0x98); + // Cause of failure: coordinates vary so much that + // length = geometry::math::sqrt(dx * dx + dy * dy); returns a value of inf for length + test_one("mysql_15_04_10", mysql_15_04_10, join_round32, end_round32, 99999.999, 0.98); #endif } diff --git a/test/algorithms/buffer/multi_point_buffer.cpp b/test/algorithms/buffer/multi_point_buffer.cpp index d5fcf0583..0e7e686c1 100644 --- a/test/algorithms/buffer/multi_point_buffer.cpp +++ b/test/algorithms/buffer/multi_point_buffer.cpp @@ -30,7 +30,7 @@ void test_all() typedef bg::model::polygon polygon; typedef bg::model::multi_point

multi_point_type; - bg::strategy::buffer::join_miter join_miter; + bg::strategy::buffer::join_round join; bg::strategy::buffer::end_flat end_flat; typedef bg::strategy::buffer::distance_symmetric < @@ -40,41 +40,39 @@ void test_all() double const pi = boost::geometry::math::pi(); - test_one("simplex1", simplex, join_miter, end_flat, 2.0 * pi, 1.0, 1.0); - test_one("simplex2", simplex, join_miter, end_flat, 22.8372, 2.0, 2.0); - test_one("simplex3", simplex, join_miter, end_flat, 44.5692, 3.0, 3.0); + test_one("simplex1", simplex, join, end_flat, 2.0 * pi, 1.0, 1.0); + test_one("simplex2", simplex, join, end_flat, 22.8372, 2.0, 2.0); + test_one("simplex3", simplex, join, end_flat, 44.5692, 3.0, 3.0); - test_one("three1", three, join_miter, end_flat, 3.0 * pi, 1.0, 1.0); - test_one("three2", three, join_miter, end_flat, 36.7592, 2.0, 2.0); - test_one("three19", three, join_miter, end_flat, 33.6914, 1.9, 1.9); - test_one("three21", three, join_miter, end_flat, 39.6394, 2.1, 2.1); - test_one("three3", three, join_miter, end_flat, 65.533, 3.0, 3.0); - - test_one("multipoint_a", multipoint_a, join_miter, end_flat, 2049.98, 14.0, 14.0); - test_one("multipoint_b", multipoint_b, join_miter, end_flat, 7109.88, 15.0, 15.0); - test_one("multipoint_b1", multipoint_b, join_miter, end_flat, 6911.89, 14.7, 14.7); - test_one("multipoint_b2", multipoint_b, join_miter, end_flat, 7174.79, 15.1, 15.1); + test_one("three1", three, join, end_flat, 3.0 * pi, 1.0, 1.0); + test_one("three2", three, join, end_flat, 36.7592, 2.0, 2.0); + test_one("three19", three, join, end_flat, 33.6914, 1.9, 1.9); + test_one("three21", three, join, end_flat, 39.6394, 2.1, 2.1); + test_one("three3", three, join, end_flat, 65.533, 3.0, 3.0); + test_one("multipoint_a", multipoint_a, join, end_flat, 2049.98, 14.0, 14.0); + test_one("multipoint_b", multipoint_b, join, end_flat, 7109.88, 15.0, 15.0); + test_one("multipoint_b1", multipoint_b, join, end_flat, 6911.89, 14.7, 14.7); + test_one("multipoint_b2", multipoint_b, join, end_flat, 7174.79, 15.1, 15.1); // Grid tests { bg::strategy::buffer::point_square point_strategy; - test_with_custom_strategies("grid_a50", - grid_a, join_miter, end_flat, + grid_a, join, end_flat, distance_strategy(0.5), side_strategy, point_strategy, 7.0); -#if defined(BOOST_GEOMETRY_BUFFER_INCLUDE_FAILING_TESTS) +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) test_with_custom_strategies("grid_a54", - grid_a, join_miter, end_flat, - distance_strategy(0.54), side_strategy, point_strategy, 99); + grid_a, join, end_flat, + distance_strategy(0.54), side_strategy, point_strategy, 7.819); #endif } test_with_custom_strategies("mysql_report_2015_02_25_1_800", - mysql_report_2015_02_25_1, join_miter, end_flat, + mysql_report_2015_02_25_1, join, end_flat, distance_strategy(6051788), side_strategy, bg::strategy::buffer::point_circle(800), 115057490003226.125, 1.0); } @@ -92,7 +90,7 @@ void test_many_points_per_circle() typedef bg::model::polygon polygon; typedef bg::model::multi_point

multi_point_type; - bg::strategy::buffer::join_miter join_miter; + bg::strategy::buffer::join_round join; bg::strategy::buffer::end_flat end_flat; typedef bg::strategy::buffer::distance_symmetric < @@ -104,41 +102,62 @@ void test_many_points_per_circle() double const tolerance = 1.0; + // Area should be somewhat larger (~>) than pi*distance^2 + // 6051788: area ~> 115058122875258 + // Strategies with many points, which are (very) slow in debug mode test_with_custom_strategies( "mysql_report_2015_02_25_1_8000", - mysql_report_2015_02_25_1, join_miter, end_flat, + mysql_report_2015_02_25_1, join, end_flat, distance_strategy(6051788), side_strategy, point_circle(8000), 115058661065242.812, tolerance); test_with_custom_strategies( "mysql_report_2015_02_25_1", - mysql_report_2015_02_25_1, join_miter, end_flat, + mysql_report_2015_02_25_1, join, end_flat, distance_strategy(6051788), side_strategy, point_circle(83585), 115058672785611.219, tolerance); - // Takes about 20 seconds in release mode + // Takes about 7 seconds in release mode test_with_custom_strategies( "mysql_report_2015_02_25_1_250k", - mysql_report_2015_02_25_1, join_miter, end_flat, + mysql_report_2015_02_25_1, join, end_flat, distance_strategy(6051788), side_strategy, point_circle(250000), 115058672880671.531, tolerance); -#if defined(BOOST_GEOMETRY_BUFFER_INCLUDE_FAILING_TESTS) - // Takes too long, TODO improve turn_in_piece_visitor +#if defined(BOOST_GEOMETRY_BUFFER_INCLUDE_SLOW_TESTS) + // Takes about 110 seconds in release mode test_with_custom_strategies( - "mysql_report_2015_02_25_1", - mysql_report_2015_02_25_1, join_miter, end_flat, + "mysql_report_2015_02_25_1_800k", + mysql_report_2015_02_25_1, join, end_flat, distance_strategy(6051788), side_strategy, point_circle(800000), - 115058672799999.999, tolerance); // area to be determined + 115058672871849.219, tolerance); #endif + // 5666962: area ~> 100890546298964 test_with_custom_strategies( "mysql_report_2015_02_25_2", - mysql_report_2015_02_25_2, join_miter, end_flat, + mysql_report_2015_02_25_2, join, end_flat, distance_strategy(5666962), side_strategy, point_circle(46641), 100891031341757.344, tolerance); + // Multipoint b with large distances/many points + // Area ~> pi * 10x + test_with_custom_strategies( + "multipoint_b_50k", + multipoint_b, join, end_flat, + distance_strategy(1000000), side_strategy, point_circle(50000), + 3141871558215.26465, tolerance); + +#if defined(BOOST_GEOMETRY_BUFFER_INCLUDE_SLOW_TESTS) + // Tests optimization min/max radius + // Takes about 55 seconds in release mode + test_with_custom_strategies( + "multipoint_b_500k", + multipoint_b, join, end_flat, + distance_strategy(10000000), side_strategy, point_circle(500000), + 314162054419515.562, tolerance); +#endif } int test_main(int, char* []) diff --git a/test/algorithms/buffer/multi_polygon_buffer.cpp b/test/algorithms/buffer/multi_polygon_buffer.cpp index eef6d5155..62a47f8da 100644 --- a/test/algorithms/buffer/multi_polygon_buffer.cpp +++ b/test/algorithms/buffer/multi_polygon_buffer.cpp @@ -413,7 +413,11 @@ void test_all() test_one("rt_p22", rt_p22, join_miter, end_flat, 26.5711, 1.0); test_one("rt_q1", rt_q1, join_miter, end_flat, 27, 1.0); +#if ! defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + // Fails because all IP's on external borders are now selectable, two are on top of each other + // and this should be fixed differently test_one("rt_q2", rt_q2, join_miter, end_flat, 26.4853, 1.0); +#endif test_one("rt_q2", rt_q2, join_miter, end_flat, 0.9697, -0.25); test_one("rt_r", rt_r, join_miter, end_flat, 21.0761, 1.0); @@ -449,8 +453,8 @@ void test_all() test_one("rt_u11_50", rt_u11, join_miter, end_flat, 0.04289, -0.50); test_one("rt_u11_25", rt_u11, join_miter, end_flat, 10.1449, -0.25); -#if defined(BOOST_GEOMETRY_BUFFER_INCLUDE_FAILING_TESTS) - test_one("rt_u12", rt_u12, join_miter, end_flat, 999, 1.0); +#if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + test_one("rt_u12", rt_u12, join_miter, end_flat, 142.1348, 1.0); test_one("rt_u13", rt_u13, join_miter, end_flat, 115.4853, 1.0); #endif diff --git a/test/algorithms/buffer/test_buffer.hpp b/test/algorithms/buffer/test_buffer.hpp index 6b0d91105..29f35ce9d 100644 --- a/test/algorithms/buffer/test_buffer.hpp +++ b/test/algorithms/buffer/test_buffer.hpp @@ -1,7 +1,7 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Unit Test +// Unit Test Helper -// Copyright (c) 2010-2014 Barend Gehrels, Amsterdam, the Netherlands. +// Copyright (c) 2010-2015 Barend Gehrels, Amsterdam, the Netherlands. // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) @@ -10,29 +10,18 @@ #ifndef BOOST_GEOMETRY_TEST_BUFFER_HPP #define BOOST_GEOMETRY_TEST_BUFFER_HPP +#if defined(TEST_WITH_SVG) + // Define before including any buffer headerfile + #define BOOST_GEOMETRY_BUFFER_USE_HELPER_POINTS +#endif + #include #include #include -#if defined(TEST_WITH_SVG) -#define BOOST_GEOMETRY_BUFFER_USE_HELPER_POINTS - -// Uncomment next lines if you want to have a zoomed view -// #define BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX - -// If possible define box before including this unit with the right view -#ifdef BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX -# ifndef BOOST_GEOMETRY_BUFFER_TEST_SVG_ALTERNATE_BOX -# define BOOST_GEOMETRY_BUFFER_TEST_SVG_ALTERNATE_BOX "BOX(0 0,100 100)" -# endif -#endif - -#endif // TEST_WITH_SVG - #include #include - #include #include #include @@ -57,296 +46,8 @@ #if defined(TEST_WITH_SVG) - -#include - - -inline char piece_type_char(bg::strategy::buffer::piece_type const& type) -{ - using namespace bg::strategy::buffer; - switch(type) - { - case buffered_segment : return 's'; - case buffered_join : return 'j'; - case buffered_round_end : return 'r'; - case buffered_flat_end : return 'f'; - case buffered_point : return 'p'; - case buffered_concave : return 'c'; - default : return '?'; - } -} - -template -void post_map(Geometry const& geometry, Mapper& mapper, RescalePolicy const& rescale_policy) -{ - typedef typename bg::point_type::type point_type; - typedef bg::detail::overlay::turn_info - < - point_type, - typename bg::segment_ratio_type::type - > turn_info; - - std::vector turns; - - bg::detail::self_get_turn_points::no_interrupt_policy policy; - bg::self_turns - < - bg::detail::overlay::assign_null_policy - >(geometry, rescale_policy, turns, policy); - - BOOST_FOREACH(turn_info const& turn, turns) - { - mapper.map(turn.point, "fill:rgb(255,128,0);stroke:rgb(0,0,100);stroke-width:1", 3); - } -} - -template -struct svg_visitor -{ -#ifdef BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX - Box m_alternate_box; -#endif - - class si - { - private : - bg::segment_identifier m_id; - - public : - inline si(bg::segment_identifier const& id) - : m_id(id) - {} - - template - inline friend std::basic_ostream& operator<<( - std::basic_ostream& os, - si const& s) - { - os << s.m_id.multi_index << "." << s.m_id.segment_index; - return os; - } - }; - - - SvgMapper& m_mapper; - - svg_visitor(SvgMapper& mapper) - : m_mapper(mapper) - {} - - template - inline void map_turns(Turns const& turns, bool label_good_turns, bool label_wrong_turns) - { - namespace bgdb = boost::geometry::detail::buffer; - typedef typename boost::range_value::type turn_type; - typedef typename turn_type::point_type point_type; - typedef typename turn_type::robust_point_type robust_point_type; - - std::map > offsets; - - for (typename boost::range_iterator::type it = - boost::begin(turns); it != boost::end(turns); ++it) - { -#ifdef BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX - if (bg::disjoint(it->point, m_alternate_box)) - { - continue; - } -#endif - - bool is_good = true; - char color = 'g'; - std::string fill = "fill:rgb(0,255,0);"; - switch(it->location) - { - case bgdb::inside_buffer : - fill = "fill:rgb(255,0,0);"; - color = 'r'; - is_good = false; - break; - case bgdb::location_discard : - fill = "fill:rgb(0,0,255);"; - color = 'b'; - is_good = false; - break; - } - if (!it->selectable_start) - { - fill = "fill:rgb(255,192,0);"; - color = 'o'; // orange - } - if (it->blocked()) - { - fill = "fill:rgb(128,128,128);"; - color = '-'; - is_good = false; - } - - fill += "fill-opacity:0.7;"; - - m_mapper.map(it->point, fill, 4); - - if ((label_good_turns && is_good) || (label_wrong_turns && ! is_good)) - { - std::ostringstream out; - out << it->turn_index - << " " << it->operations[0].piece_index << "/" << it->operations[1].piece_index - << " " << si(it->operations[0].seg_id) << "/" << si(it->operations[1].seg_id) - - // If you want to see travel information - << std::endl - << " nxt " << it->operations[0].enriched.travels_to_ip_index - << "/" << it->operations[1].enriched.travels_to_ip_index - << " or " << it->operations[0].enriched.next_ip_index - << "/" << it->operations[1].enriched.next_ip_index - //<< " frac " << it->operations[0].fraction - - // If you want to see robust-point coordinates (e.g. to find duplicates) - // << std::endl - // << " " << bg::get<0>(it->robust_point) << " , " << bg::get<1>(it->robust_point) - - << std::endl; - out << " " << bg::method_char(it->method) - << ":" << bg::operation_char(it->operations[0].operation) - << "/" << bg::operation_char(it->operations[1].operation); - out << " " - << (it->count_on_offsetted > 0 ? "b" : "") // b: offsetted border - << (it->count_within_near_offsetted > 0 ? "n" : "") - << (it->count_within > 0 ? "w" : "") - << (it->count_on_helper > 0 ? "h" : "") - << (it->count_on_multi > 0 ? "m" : "") - ; - - offsets[it->get_robust_point()] += 10; - int offset = offsets[it->get_robust_point()]; - - m_mapper.text(it->point, out.str(), "fill:rgb(0,0,0);font-family='Arial';font-size:9px;", 5, offset); - - offsets[it->get_robust_point()] += 25; - } - } - } - - template - inline void map_pieces(Pieces const& pieces, - OffsettedRings const& offsetted_rings, - bool do_pieces, bool do_indices) - { - typedef typename boost::range_value::type piece_type; - typedef typename boost::range_value::type ring_type; - - for(typename boost::range_iterator::type it = boost::begin(pieces); - it != boost::end(pieces); - ++it) - { - const piece_type& piece = *it; - bg::segment_identifier seg_id = piece.first_seg_id; - if (seg_id.segment_index < 0) - { - continue; - } - - ring_type corner; - - - ring_type const& ring = offsetted_rings[seg_id.multi_index]; - - std::copy(boost::begin(ring) + seg_id.segment_index, - boost::begin(ring) + piece.last_segment_index, - std::back_inserter(corner)); - std::copy(boost::begin(piece.helper_points), - boost::end(piece.helper_points), - std::back_inserter(corner)); - - if (corner.empty()) - { - continue; - } -#ifdef BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX - if (bg::disjoint(corner, m_alternate_box)) - { - continue; - } -#endif - - if (do_pieces) - { - std::string style = "opacity:0.3;stroke:rgb(0,0,0);stroke-width:1;"; - m_mapper.map(corner, - piece.type == bg::strategy::buffer::buffered_segment - ? style + "fill:rgb(255,128,0);" - : style + "fill:rgb(255,0,0);"); - } - - if (do_indices) - { - // Label starting piece_index / segment_index - typedef typename bg::point_type::type point_type; - - std::ostringstream out; - out << piece.index << " (" << piece_type_char(piece.type) << ") " << piece.first_seg_id.segment_index << ".." << piece.last_segment_index - 1; - point_type label_point = bg::return_centroid(corner); - - if ((piece.type == bg::strategy::buffer::buffered_concave - || piece.type == bg::strategy::buffer::buffered_flat_end) - && corner.size() >= 2u) - { - bg::set<0>(label_point, (bg::get<0>(corner[0]) + bg::get<0>(corner[1])) / 2.0); - bg::set<1>(label_point, (bg::get<1>(corner[0]) + bg::get<1>(corner[1])) / 2.0); - } - m_mapper.text(label_point, out.str(), "fill:rgb(255,0,0);font-family='Arial';font-size:10px;", 5, 5); - } - } - } - - template - inline void map_traversed_rings(TraversedRings const& traversed_rings) - { - for(typename boost::range_iterator::type it - = boost::begin(traversed_rings); it != boost::end(traversed_rings); ++it) - { - m_mapper.map(*it, "opacity:0.4;fill:none;stroke:rgb(0,255,0);stroke-width:2"); - } - } - - template - inline void map_offsetted_rings(OffsettedRings const& offsetted_rings) - { - for(typename boost::range_iterator::type it - = boost::begin(offsetted_rings); it != boost::end(offsetted_rings); ++it) - { - if (it->discarded()) - { - m_mapper.map(*it, "opacity:0.4;fill:none;stroke:rgb(255,0,0);stroke-width:2"); - } - else - { - m_mapper.map(*it, "opacity:0.4;fill:none;stroke:rgb(0,0,255);stroke-width:2"); - } - } - } - - template - inline void apply(PieceCollection const& collection, int phase) - { - // Comment next return if you want to see pieces, turns, etc. - return; - - if(phase == 0) - { - map_pieces(collection.m_pieces, collection.offsetted_rings, true, true); - map_turns(collection.m_turns, true, false); - } -#if !defined(BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX) - if (phase == 1) - { -// map_traversed_rings(collection.traversed_rings); -// map_offsetted_rings(collection.offsetted_rings); - } -#endif - } -}; - +# include +# include #endif //----------------------------------------------------------------------------- @@ -472,40 +173,25 @@ void test_buffer(std::string const& caseid, Geometry const& geometry, //std::cout << complete.str() << std::endl; +#if defined(TEST_WITH_SVG_PER_TURN) + save_turns_visitor visitor; +#elif defined(TEST_WITH_SVG) + + buffer_svg_mapper buffer_mapper(complete.str()); + std::ostringstream filename; filename << "buffer_" << complete.str() << ".svg"; - -#if defined(TEST_WITH_SVG) std::ofstream svg(filename.str().c_str()); typedef bg::svg_mapper mapper_type; - mapper_type mapper(svg, 1000, 1000); + mapper_type mapper(svg, 1000, 800); - svg_visitor, tag> visitor(mapper); - -#ifdef BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX - // Create a zoomed-in view - bg::model::box alternate_box; - bg::read_wkt(BOOST_GEOMETRY_BUFFER_TEST_SVG_ALTERNATE_BOX, alternate_box); - mapper.add(alternate_box); - - // Take care non-visible elements are skipped - visitor.m_alternate_box = alternate_box; -#else - - { - bg::model::box box = envelope; - if (distance_strategy.negative()) - { - bg::buffer(box, box, 1.0); - } - else - { - bg::buffer(box, box, 1.1 * distance_strategy.max_distance(join_strategy, end_strategy)); - } - mapper.add(box); - } -#endif + svg_visitor > visitor(mapper); + buffer_mapper.prepare(mapper, visitor, envelope, + distance_strategy.negative() + ? 1.0 + : 1.1 * distance_strategy.max_distance(join_strategy, end_strategy) + ); #else bg::detail::buffer::visit_pieces_default_policy visitor; #endif @@ -557,14 +243,22 @@ void test_buffer(std::string const& caseid, Geometry const& geometry, if (check_self_intersections) { - // Be sure resulting polygon does not contain - // self-intersections - BOOST_CHECK_MESSAGE - ( - ! bg::detail::overlay::has_self_intersections(buffered, - rescale_policy, false), - complete.str() << " output is self-intersecting. " - ); + try + { + bool has_self_ips = bg::detail::overlay::has_self_intersections(buffered, + rescale_policy, false); + // Be sure resulting polygon does not contain + // self-intersections + BOOST_CHECK_MESSAGE + ( + ! has_self_ips, + complete.str() << " output is self-intersecting. " + ); + } + catch(...) + { + BOOST_MESSAGE("Exception in checking self-intersections"); + } } } @@ -579,48 +273,28 @@ void test_buffer(std::string const& caseid, Geometry const& geometry, // BOOST_CHECK_MESSAGE(bg::intersects(buffered) == false, complete.str() << " intersects"); #endif -#if defined(TEST_WITH_SVG) - bool const areal = boost::is_same - < - typename bg::tag_cast::type, bg::areal_tag - >::type::value; - - // Map input geometry in green - if (areal) +#if defined(TEST_WITH_SVG_PER_TURN) { -#ifdef BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX_FOR_INPUT - // Assuming input is areal - bg::model::multi_polygon clipped; - bg::intersection(geometry, alternate_box, clipped); - mapper.map(clipped, "opacity:0.5;fill:rgb(0,128,0);stroke:rgb(0,64,0);stroke-width:2"); -#else - mapper.map(geometry, "opacity:0.5;fill:rgb(0,128,0);stroke:rgb(0,64,0);stroke-width:2"); + // Create a per turn visitor to map per turn, and buffer again with it + per_turn_visitor ptv(complete.str(), visitor.get_points()); + bg::detail::buffer::buffer_inserter(geometry, + std::back_inserter(buffered), + distance_strategy, + side_strategy, + join_strategy, + end_strategy, + point_strategy, + rescale_policy, + ptv); + ptv.map_input_output(geometry, buffered, distance_strategy.negative()); + // self_ips NYI here + } +#elif defined(TEST_WITH_SVG) + buffer_mapper.map_input_output(mapper, geometry, buffered, distance_strategy.negative()); + buffer_mapper.map_self_ips(mapper, buffered, rescale_policy); #endif - } - else - { - // TODO: clip input points/linestring - mapper.map(geometry, "opacity:0.5;stroke:rgb(0,128,0);stroke-width:10"); - } - - { - // Map buffer in yellow (inflate) and with orange-dots (deflate) - std::string style = distance_strategy.negative() - ? "opacity:0.4;fill:rgb(255,255,192);stroke:rgb(255,128,0);stroke-width:3" - : "opacity:0.4;fill:rgb(255,255,128);stroke:rgb(0,0,0);stroke-width:3"; - -#ifdef BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX - // Clip output multi-polygon with box - bg::model::multi_polygon clipped; - bg::intersection(buffered, alternate_box, clipped); - mapper.map(clipped, style); -#else - mapper.map(buffered, style); -#endif - } - post_map(buffered, mapper, rescale_policy); -#endif // TEST_WITH_SVG + // Check for self-intersections if (self_ip_count != NULL) { std::size_t count = 0; diff --git a/test/algorithms/buffer/test_buffer_svg.hpp b/test/algorithms/buffer/test_buffer_svg.hpp new file mode 100644 index 000000000..2acf3376e --- /dev/null +++ b/test/algorithms/buffer/test_buffer_svg.hpp @@ -0,0 +1,468 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit Test Helper + +// Copyright (c) 2010-2015 Barend Gehrels, Amsterdam, the Netherlands. +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + + +#ifndef BOOST_GEOMETRY_TEST_BUFFER_SVG_HPP +#define BOOST_GEOMETRY_TEST_BUFFER_SVG_HPP + +#include +#include + +// Uncomment next lines if you want to have a zoomed view +//#define BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX + +// If possible define box before including this unit with the right view +#ifdef BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX +# ifndef BOOST_GEOMETRY_BUFFER_TEST_SVG_ALTERNATE_BOX +# define BOOST_GEOMETRY_BUFFER_TEST_SVG_ALTERNATE_BOX "BOX(0 0,100 100)" +# endif +#endif + +#include +#include + + +inline char piece_type_char(bg::strategy::buffer::piece_type const& type) +{ + using namespace bg::strategy::buffer; + switch(type) + { + case buffered_segment : return 's'; + case buffered_join : return 'j'; + case buffered_round_end : return 'r'; + case buffered_flat_end : return 'f'; + case buffered_point : return 'p'; + case buffered_concave : return 'c'; + default : return '?'; + } +} + +template +class svg_visitor +{ +public : + svg_visitor(SvgMapper& mapper) + : m_mapper(mapper) + , m_zoom(false) + {} + + void set_alternate_box(Box const& box) + { + m_alternate_box = box; + m_zoom = true; + } + + template + inline void apply(PieceCollection const& collection, int phase) + { + // Comment next return if you want to see pieces, turns, etc. + return; + + if(phase == 0) + { + map_pieces(collection.m_pieces, collection.offsetted_rings, true, true); + map_turns(collection.m_turns, true, false); + } + if (phase == 1 && ! m_zoom) + { +// map_traversed_rings(collection.traversed_rings); +// map_offsetted_rings(collection.offsetted_rings); + } + } + +private : + class si + { + private : + bg::segment_identifier m_id; + + public : + inline si(bg::segment_identifier const& id) + : m_id(id) + {} + + template + inline friend std::basic_ostream& operator<<( + std::basic_ostream& os, + si const& s) + { + os << s.m_id.multi_index << "." << s.m_id.segment_index; + return os; + } + }; + + template + inline void map_turns(Turns const& turns, bool label_good_turns, bool label_wrong_turns) + { + namespace bgdb = boost::geometry::detail::buffer; + typedef typename boost::range_value::type turn_type; + typedef typename turn_type::point_type point_type; + typedef typename turn_type::robust_point_type robust_point_type; + + std::map > offsets; + + for (typename boost::range_iterator::type it = + boost::begin(turns); it != boost::end(turns); ++it) + { + if (m_zoom && bg::disjoint(it->point, m_alternate_box)) + { + continue; + } + + bool is_good = true; + char color = 'g'; + std::string fill = "fill:rgb(0,255,0);"; + switch(it->location) + { + case bgdb::inside_buffer : + fill = "fill:rgb(255,0,0);"; + color = 'r'; + is_good = false; + break; + case bgdb::location_discard : + fill = "fill:rgb(0,0,255);"; + color = 'b'; + is_good = false; + break; + } + if (!it->selectable_start) + { + fill = "fill:rgb(255,192,0);"; + color = 'o'; // orange + } + if (it->blocked()) + { + fill = "fill:rgb(128,128,128);"; + color = '-'; + is_good = false; + } + + fill += "fill-opacity:0.7;"; + + m_mapper.map(it->point, fill, 4); + + if ((label_good_turns && is_good) || (label_wrong_turns && ! is_good)) + { + std::ostringstream out; + out << it->turn_index + << " " << it->operations[0].piece_index << "/" << it->operations[1].piece_index + << " " << si(it->operations[0].seg_id) << "/" << si(it->operations[1].seg_id) + + // If you want to see travel information + << std::endl + << " nxt " << it->operations[0].enriched.travels_to_ip_index + << "/" << it->operations[1].enriched.travels_to_ip_index + << " or " << it->operations[0].enriched.next_ip_index + << "/" << it->operations[1].enriched.next_ip_index + //<< " frac " << it->operations[0].fraction + + // If you want to see robust-point coordinates (e.g. to find duplicates) + // << std::endl + // << " " << bg::get<0>(it->robust_point) << " , " << bg::get<1>(it->robust_point) + + << std::endl; + out << " " << bg::method_char(it->method) + << ":" << bg::operation_char(it->operations[0].operation) + << "/" << bg::operation_char(it->operations[1].operation); + out << " " + << (it->count_on_offsetted > 0 ? "b" : "") // b: offsetted border +#if ! defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION) + << (it->count_within_near_offsetted > 0 ? "n" : "") +#endif + << (it->count_within > 0 ? "w" : "") + << (it->count_on_helper > 0 ? "h" : "") + << (it->count_on_multi > 0 ? "m" : "") + ; + + offsets[it->get_robust_point()] += 10; + int offset = offsets[it->get_robust_point()]; + + m_mapper.text(it->point, out.str(), "fill:rgb(0,0,0);font-family='Arial';font-size:9px;", 5, offset); + + offsets[it->get_robust_point()] += 25; + } + } + } + + template + inline void map_pieces(Pieces const& pieces, + OffsettedRings const& offsetted_rings, + bool do_pieces, bool do_indices) + { + typedef typename boost::range_value::type piece_type; + typedef typename boost::range_value::type ring_type; + + for(typename boost::range_iterator::type it = boost::begin(pieces); + it != boost::end(pieces); + ++it) + { + const piece_type& piece = *it; + bg::segment_identifier seg_id = piece.first_seg_id; + if (seg_id.segment_index < 0) + { + continue; + } + + ring_type corner; + + + ring_type const& ring = offsetted_rings[seg_id.multi_index]; + + std::copy(boost::begin(ring) + seg_id.segment_index, + boost::begin(ring) + piece.last_segment_index, + std::back_inserter(corner)); + std::copy(boost::begin(piece.helper_points), + boost::end(piece.helper_points), + std::back_inserter(corner)); + + if (corner.empty()) + { + continue; + } + if (m_zoom && bg::disjoint(corner, m_alternate_box)) + { + continue; + } + + if (m_zoom && do_pieces) + { + try + { + std::string style = "opacity:0.3;stroke:rgb(0,0,0);stroke-width:1;"; + typedef typename bg::point_type::type point_type; + bg::model::multi_polygon > clipped; + bg::intersection(ring, m_alternate_box, clipped); + m_mapper.map(clipped, + piece.type == bg::strategy::buffer::buffered_segment + ? style + "fill:rgb(255,128,0);" + : style + "fill:rgb(255,0,0);"); + } + catch (...) + { + std::cout << "Error for piece " << piece.index << std::endl; + } + } + else if (do_pieces) + { + std::string style = "opacity:0.3;stroke:rgb(0,0,0);stroke-width:1;"; + m_mapper.map(corner, + piece.type == bg::strategy::buffer::buffered_segment + ? style + "fill:rgb(255,128,0);" + : style + "fill:rgb(255,0,0);"); + } + + if (do_indices) + { + // Label starting piece_index / segment_index + typedef typename bg::point_type::type point_type; + + std::ostringstream out; + out << piece.index << " (" << piece_type_char(piece.type) << ") " << piece.first_seg_id.segment_index << ".." << piece.last_segment_index - 1; + point_type label_point = bg::return_centroid(corner); + + if ((piece.type == bg::strategy::buffer::buffered_concave + || piece.type == bg::strategy::buffer::buffered_flat_end) + && corner.size() >= 2u) + { + bg::set<0>(label_point, (bg::get<0>(corner[0]) + bg::get<0>(corner[1])) / 2.0); + bg::set<1>(label_point, (bg::get<1>(corner[0]) + bg::get<1>(corner[1])) / 2.0); + } + m_mapper.text(label_point, out.str(), "fill:rgb(255,0,0);font-family='Arial';font-size:10px;", 5, 5); + } + } + } + + template + inline void map_traversed_rings(TraversedRings const& traversed_rings) + { + for(typename boost::range_iterator::type it + = boost::begin(traversed_rings); it != boost::end(traversed_rings); ++it) + { + m_mapper.map(*it, "opacity:0.4;fill:none;stroke:rgb(0,255,0);stroke-width:2"); + } + } + + template + inline void map_offsetted_rings(OffsettedRings const& offsetted_rings) + { + for(typename boost::range_iterator::type it + = boost::begin(offsetted_rings); it != boost::end(offsetted_rings); ++it) + { + if (it->discarded()) + { + m_mapper.map(*it, "opacity:0.4;fill:none;stroke:rgb(255,0,0);stroke-width:2"); + } + else + { + m_mapper.map(*it, "opacity:0.4;fill:none;stroke:rgb(0,0,255);stroke-width:2"); + } + } + } + + + SvgMapper& m_mapper; + Box m_alternate_box; + bool m_zoom; + +}; + +template +class buffer_svg_mapper +{ +public : + + buffer_svg_mapper(std::string const& casename) + : m_casename(casename) + , m_zoom(false) + {} + + template + void prepare(Mapper& mapper, Visitor& visitor, Envelope const& envelope, double box_buffer_distance) + { +#ifdef BOOST_GEOMETRY_BUFFER_TEST_SVG_USE_ALTERNATE_BOX + // Create a zoomed-in view + bg::read_wkt(BOOST_GEOMETRY_BUFFER_TEST_SVG_ALTERNATE_BOX, alternate_box); + mapper.add(alternate_box); + + // Take care non-visible elements are skipped + visitor.m_alternate_box = alternate_box; + m_zoom = true; +#else + bg::model::box box = envelope; + bg::buffer(box, box, box_buffer_distance); + mapper.add(box); +#endif + + boost::ignore_unused(visitor); + } + + void set_alternate_box(bg::model::box const& box) + { + m_alternate_box = box; + m_zoom = true; + } + + template + void map_input_output(Mapper& mapper, Geometry const& geometry, + GeometryBuffer const& buffered, bool negative) + { + bool const areal = boost::is_same + < + typename bg::tag_cast + < + typename bg::tag::type, + bg::areal_tag + >::type, bg::areal_tag + >::type::value; + + if (m_zoom) + { + map_io_zoomed(mapper, geometry, buffered, negative, areal); + } + else + { + map_io(mapper, geometry, buffered, negative, areal); + } + } + + template + void map_self_ips(Mapper& mapper, Geometry const& geometry, RescalePolicy const& rescale_policy) + { + typedef bg::detail::overlay::turn_info + < + Point, + typename bg::segment_ratio_type::type + > turn_info; + + std::vector turns; + + bg::detail::self_get_turn_points::no_interrupt_policy policy; + bg::self_turns + < + bg::detail::overlay::assign_null_policy + >(geometry, rescale_policy, turns, policy); + + BOOST_FOREACH(turn_info const& turn, turns) + { + mapper.map(turn.point, "fill:rgb(255,128,0);stroke:rgb(0,0,100);stroke-width:1", 3); + } + } + +private : + + template + void map_io(Mapper& mapper, Geometry const& geometry, + GeometryBuffer const& buffered, bool negative, bool areal) + { + // Map input geometry in green + if (areal) + { + mapper.map(geometry, "opacity:0.5;fill:rgb(0,128,0);stroke:rgb(0,64,0);stroke-width:2"); + } + else + { + // TODO: clip input points/linestring + mapper.map(geometry, "opacity:0.5;stroke:rgb(0,128,0);stroke-width:10"); + } + + { + // Map buffer in yellow (inflate) and with orange-dots (deflate) + std::string style = negative + ? "opacity:0.4;fill:rgb(255,255,192);stroke:rgb(255,128,0);stroke-width:3" + : "opacity:0.4;fill:rgb(255,255,128);stroke:rgb(0,0,0);stroke-width:3"; + + mapper.map(buffered, style); + } + } + + template + void map_io_zoomed(Mapper& mapper, Geometry const& geometry, + GeometryBuffer const& buffered, bool negative, bool areal) + { + // Map input geometry in green + if (areal) + { + // Assuming input is areal + GeometryBuffer clipped; +// TODO: the next line does NOT compile for multi-point, TODO: implement that line +// bg::intersection(geometry, m_alternate_box, clipped); + mapper.map(clipped, "opacity:0.5;fill:rgb(0,128,0);stroke:rgb(0,64,0);stroke-width:2"); + } + else + { + // TODO: clip input (multi)point/linestring + mapper.map(geometry, "opacity:0.5;stroke:rgb(0,128,0);stroke-width:10"); + } + + { + // Map buffer in yellow (inflate) and with orange-dots (deflate) + std::string style = negative + ? "opacity:0.4;fill:rgb(255,255,192);stroke:rgb(255,128,0);stroke-width:3" + : "opacity:0.4;fill:rgb(255,255,128);stroke:rgb(0,0,0);stroke-width:3"; + + try + { + // Clip output multi-polygon with box + GeometryBuffer clipped; + bg::intersection(buffered, m_alternate_box, clipped); + mapper.map(clipped, style); + } + catch (...) + { + std::cout << "Error for buffered output " << m_casename << std::endl; + } + } + } + + bg::model::box m_alternate_box; + bool m_zoom; + std::string m_casename; +}; + + +#endif diff --git a/test/algorithms/buffer/test_buffer_svg_per_turn.hpp b/test/algorithms/buffer/test_buffer_svg_per_turn.hpp new file mode 100644 index 000000000..1a06d6fbf --- /dev/null +++ b/test/algorithms/buffer/test_buffer_svg_per_turn.hpp @@ -0,0 +1,164 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit Test Helper + +// Copyright (c) 2010-2015 Barend Gehrels, Amsterdam, the Netherlands. +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + + +#ifndef BOOST_GEOMETRY_TEST_BUFFER_SVG_PER_TURN_HPP +#define BOOST_GEOMETRY_TEST_BUFFER_SVG_PER_TURN_HPP + +#include +#include + +#include +#include + +template +class save_turns_visitor +{ +public : + typedef std::vector > vector_type; + + template + inline void get_turns(Turns const& turns) + { + for (typename boost::range_iterator::type it = + boost::begin(turns); it != boost::end(turns); ++it) + { + m_points.push_back(std::make_pair(it->point, it->turn_index)); + } + } + + template + inline void apply(PieceCollection const& collection, int phase) + { + if (phase == 0) + { + get_turns(collection.m_turns); + } + } + + vector_type const& get_points() { return m_points; } + +private : + vector_type m_points; +}; + + + +template +class mapper_visitor +{ +public : + mapper_visitor(std::string const& complete, int index, Point const& point) + : m_filename(get_filename(complete, index)) + , m_svg(m_filename.c_str()) + , m_mapper(m_svg, 1000, 800) + , m_visitor(m_mapper) + , m_buffer_mapper(m_filename) + { + box_type box; + double half_size = 75.0; // specific for multi_point buffer + bg::set(box, bg::get<0>(point) - half_size); + bg::set(box, bg::get<1>(point) - half_size); + bg::set(box, bg::get<0>(point) + half_size); + bg::set(box, bg::get<1>(point) + half_size); + + m_mapper.add(box); + m_visitor.set_alternate_box(box); + m_buffer_mapper.set_alternate_box(box); + } + + // It is used in a ptr vector + virtual ~mapper_visitor() + { + } + + template + inline void apply(PieceCollection const& collection, int phase) + { + m_visitor.apply(collection, phase); + } + + template + void map_input_output(Geometry const& geometry, + GeometryBuffer const& buffered, bool negative) + { + m_buffer_mapper.map_input_output(m_mapper, geometry, buffered, negative); + } + +private : + std::string get_filename(std::string const& complete, int index) + { + std::ostringstream filename; + filename << "buffer_per_turn_" << complete << "_" << index << ".svg"; + return filename.str(); + } + + typedef bg::svg_mapper mapper_type; + typedef bg::model::box box_type; + + std::string m_filename; + std::ofstream m_svg; + mapper_type m_mapper; + svg_visitor m_visitor; + buffer_svg_mapper m_buffer_mapper; +}; + +template +class per_turn_visitor +{ + // Both fstreams and svg mappers are non-copyable, + // therefore we need to use dynamic memory + typedef boost::ptr_vector > container_type; + container_type mappers; + +public : + + typedef std::pair pair_type; + typedef std::vector vector_type; + + per_turn_visitor(std::string const& complete_caseid, + vector_type const& points) + { + namespace bg = boost::geometry; + + if (points.size() > 100u) + { + // Defensive check. Too much intersections. Don't create anything + return; + } + + BOOST_FOREACH(pair_type const& p, points) + { + mappers.push_back(new mapper_visitor(complete_caseid, p.second, p.first)); + } + } + + template + inline void apply(PieceCollection const& collection, int phase) + { + for(typename container_type::iterator it = mappers.begin(); + it != mappers.end(); ++it) + { + it->apply(collection, phase); + } + } + + template + void map_input_output(Geometry const& geometry, + GeometryBuffer const& buffered, bool negative) + { + for(typename container_type::iterator it = mappers.begin(); + it != mappers.end(); ++it) + { + it->map_input_output(geometry, buffered, negative); + } + } +}; + + +#endif // BOOST_GEOMETRY_TEST_BUFFER_SVG_PER_TURN_HPP diff --git a/test/algorithms/distance/distance_se_pl_l.cpp b/test/algorithms/distance/distance_se_pl_l.cpp index 65c80f8d2..4fcfaee3f 100644 --- a/test/algorithms/distance/distance_se_pl_l.cpp +++ b/test/algorithms/distance/distance_se_pl_l.cpp @@ -122,19 +122,19 @@ void test_distance_point_segment(Strategy const& strategy) #endif typedef test_distance_of_geometries tester; + double const d2r = bg::math::d2r(); + tester::apply("p-s-01", "POINT(0 0)", "SEGMENT(2 0,3 0)", - 2.0 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 2.0 * bg::math::d2r * strategy.radius()), + 2.0 * d2r * strategy.radius(), + to_comparable(strategy, 2.0 * d2r * strategy.radius()), strategy); tester::apply("p-s-02", "POINT(2.5 3)", "SEGMENT(2 0,3 0)", - 3.0 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 3.0 * bg::math::d2r * strategy.radius()), + 3.0 * d2r * strategy.radius(), + to_comparable(strategy, 3.0 * d2r * strategy.radius()), strategy); tester::apply("p-s-03", "POINT(2 0)", @@ -163,44 +163,38 @@ void test_distance_point_segment(Strategy const& strategy) tester::apply("p-s-07", "POINT(90 1e-3)", "SEGMENT(0.5 0,175.5 0)", - 1e-3 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 1e-3 * bg::math::d2r * strategy.radius()), + 1e-3 * d2r * strategy.radius(), + to_comparable(strategy, 1e-3 * d2r * strategy.radius()), strategy); tester::apply("p-s-08", "POINT(90 1e-4)", "SEGMENT(0.5 0,175.5 0)", - 1e-4 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 1e-4 * bg::math::d2r * strategy.radius()), + 1e-4 * d2r * strategy.radius(), + to_comparable(strategy, 1e-4 * d2r * strategy.radius()), strategy); tester::apply("p-s-09", "POINT(90 1e-5)", "SEGMENT(0.5 0,175.5 0)", - 1e-5 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 1e-5 * bg::math::d2r * strategy.radius()), + 1e-5 * d2r * strategy.radius(), + to_comparable(strategy, 1e-5 * d2r * strategy.radius()), strategy); tester::apply("p-s-10", "POINT(90 1e-6)", "SEGMENT(0.5 0,175.5 0)", - 1e-6 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 1e-6 * bg::math::d2r * strategy.radius()), + 1e-6 * d2r * strategy.radius(), + to_comparable(strategy, 1e-6 * d2r * strategy.radius()), strategy); tester::apply("p-s-11", "POINT(90 1e-7)", "SEGMENT(0.5 0,175.5 0)", - 1e-7 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 1e-7 * bg::math::d2r * strategy.radius()), + 1e-7 * d2r * strategy.radius(), + to_comparable(strategy, 1e-7 * d2r * strategy.radius()), strategy); tester::apply("p-s-12", "POINT(90 1e-8)", "SEGMENT(0.5 0,175.5 0)", - 1e-8 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 1e-8 * bg::math::d2r * strategy.radius()), + 1e-8 * d2r * strategy.radius(), + to_comparable(strategy, 1e-8 * d2r * strategy.radius()), strategy); } @@ -215,26 +209,25 @@ void test_distance_point_linestring(Strategy const& strategy) #endif typedef test_distance_of_geometries tester; + double const d2r = bg::math::d2r(); + tester::apply("p-l-01", "POINT(0 0)", "LINESTRING(2 0,2 0)", - 2.0 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 2.0 * bg::math::d2r * strategy.radius()), + 2.0 * d2r * strategy.radius(), + to_comparable(strategy, 2.0 * d2r * strategy.radius()), strategy); tester::apply("p-l-02", "POINT(0 0)", "LINESTRING(2 0,3 0)", - 2.0 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 2.0 * bg::math::d2r * strategy.radius()), + 2.0 * d2r * strategy.radius(), + to_comparable(strategy, 2.0 * d2r * strategy.radius()), strategy); tester::apply("p-l-03", "POINT(2.5 3)", "LINESTRING(2 0,3 0)", - 3.0 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 3.0 * bg::math::d2r * strategy.radius()), + 3.0 * d2r * strategy.radius(), + to_comparable(strategy, 3.0 * d2r * strategy.radius()), strategy); tester::apply("p-l-04", "POINT(2 0)", @@ -254,9 +247,8 @@ void test_distance_point_linestring(Strategy const& strategy) tester::apply("p-l-07", "POINT(7.5 10)", "LINESTRING(1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 0)", - 10.0 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 10.0 * bg::math::d2r * strategy.radius()), + 10.0 * d2r * strategy.radius(), + to_comparable(strategy, 10.0 * d2r * strategy.radius()), strategy); tester::apply("p-l-08", "POINT(7.5 10)", @@ -282,19 +274,19 @@ void test_distance_point_multilinestring(Strategy const& strategy) point_type, multi_linestring_type > tester; + double const d2r = bg::math::d2r(); + tester::apply("p-ml-01", "POINT(0 0)", "MULTILINESTRING((-5 0,-3 0),(2 0,3 0))", - 2.0 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 2.0 * bg::math::d2r * strategy.radius()), + 2.0 * d2r * strategy.radius(), + to_comparable(strategy, 2.0 * d2r * strategy.radius()), strategy); tester::apply("p-ml-02", "POINT(2.5 3)", "MULTILINESTRING((-5 0,-3 0),(2 0,3 0))", - 3.0 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 3.0 * bg::math::d2r * strategy.radius()), + 3.0 * d2r * strategy.radius(), + to_comparable(strategy, 3.0 * d2r * strategy.radius()), strategy); tester::apply("p-ml-03", "POINT(2 0)", @@ -431,6 +423,8 @@ void test_distance_multipoint_segment(Strategy const& strategy) #endif typedef test_distance_of_geometries tester; + double d2r = bg::math::d2r(); + tester::apply("mp-s-01", "MULTIPOINT(0 0,1 0,0 1,1 1)", "SEGMENT(2 0,0 2)", @@ -442,9 +436,8 @@ void test_distance_multipoint_segment(Strategy const& strategy) tester::apply("mp-s-02", "MULTIPOINT(0 0,1 0,0 1,1 1)", "SEGMENT(0 -3,1 -10)", - 3.0 * bg::math::d2r * strategy.radius(), - to_comparable(strategy, - 3.0 * bg::math::d2r * strategy.radius()), + 3.0 * d2r * strategy.radius(), + to_comparable(strategy, 3.0 * d2r * strategy.radius()), strategy); tester::apply("mp-s-03", "MULTIPOINT(0 0,1 0,0 1,1 1)", diff --git a/test/algorithms/distance/distance_se_pl_pl.cpp b/test/algorithms/distance/distance_se_pl_pl.cpp index 2e14e8892..c7dfda2a9 100644 --- a/test/algorithms/distance/distance_se_pl_pl.cpp +++ b/test/algorithms/distance/distance_se_pl_pl.cpp @@ -89,7 +89,7 @@ void test_distance_point_point(Strategy const& strategy, "POINT(180 -10)", (is_comparable_strategy ? 1.0 - : (180.0 * bg::math::d2r * strategy.radius())), + : (bg::math::pi() * strategy.radius())), 1.0, strategy); tester::apply("p-p-04", @@ -97,7 +97,7 @@ void test_distance_point_point(Strategy const& strategy, "POINT(180 0)", (is_comparable_strategy ? 1.0 - : (180.0 * bg::math::d2r * strategy.radius())), + : (bg::math::pi() * strategy.radius())), 1.0, strategy); } diff --git a/test/algorithms/is_convex.cpp b/test/algorithms/is_convex.cpp new file mode 100644 index 000000000..7fcc8a913 --- /dev/null +++ b/test/algorithms/is_convex.cpp @@ -0,0 +1,61 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit Test + +// Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands. + +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#include +#include + +#include + +#include + +#include +#include +#include +#include +#include + + +template +void test_one(std::string const& case_id, std::string const& wkt, bool expected) +{ + Geometry geometry; + bg::read_wkt(wkt, geometry); + bg::correct(geometry); + + bool detected = bg::is_convex(geometry); + BOOST_CHECK_MESSAGE(detected == expected, + "Not as expected, case: " << case_id + << " / expected: " << expected + << " / detected: " << detected); +} + + +template +void test_all() +{ + // rectangular, with concavity + std::string const concave1 = "polygon((1 1, 1 4, 3 4, 3 3, 4 3, 4 4, 5 4, 5 1, 1 1))"; + std::string const triangle = "polygon((1 1, 1 4, 5 1, 1 1))"; + + test_one >("triangle", triangle, true); + test_one >("concave1", concave1, false); + test_one >("triangle", triangle, true); + test_one >("concave1", concave1, false); + + test_one >("box", "box(0 0,2 2)", true); +} + + +int test_main(int, char* []) +{ + test_all >(); + test_all >(); + + return 0; +} diff --git a/test/algorithms/relational_operations/Jamfile.v2 b/test/algorithms/relational_operations/Jamfile.v2 index 2c0c34c9e..39cf1381e 100644 --- a/test/algorithms/relational_operations/Jamfile.v2 +++ b/test/algorithms/relational_operations/Jamfile.v2 @@ -1,11 +1,11 @@ # Boost.Geometry (aka GGL, Generic Geometry Library) # -# Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. -# Copyright (c) 2008-2014 Bruno Lalande, Paris, France. -# Copyright (c) 2009-2014 Mateusz Loskot, London, UK. +# Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands. +# Copyright (c) 2008-2015 Bruno Lalande, Paris, France. +# Copyright (c) 2009-2015 Mateusz Loskot, London, UK. # -# This file was modified by Oracle on 2014. -# Modifications copyright (c) 2014, Oracle and/or its affiliates. +# This file was modified by Oracle on 2014, 2015. +# Modifications copyright (c) 2014-2015, Oracle and/or its affiliates. # # Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle # Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -19,6 +19,7 @@ test-suite boost-geometry-algorithms-relational [ run covered_by.cpp ] [ run crosses.cpp ] [ run equals.cpp ] + [ run equals_on_spheroid.cpp ] [ run intersects.cpp ] [ run multi_covered_by.cpp ] [ run multi_equals.cpp ] diff --git a/test/algorithms/relational_operations/equals_on_spheroid.cpp b/test/algorithms/relational_operations/equals_on_spheroid.cpp new file mode 100644 index 000000000..80c8593e8 --- /dev/null +++ b/test/algorithms/relational_operations/equals_on_spheroid.cpp @@ -0,0 +1,183 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit test + +// Copyright (c) 2015, Oracle and/or its affiliates. + +// Licensed under the Boost Software License version 1.0. +// http://www.boost.org/users/license.html + +// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle + +#ifndef BOOST_TEST_MODULE +#define BOOST_TEST_MODULE test_equals_on_spheroid +#endif + +#include + +#include + +#include "test_equals.hpp" + +#include + +#include + +namespace bgm = bg::model; + +template +struct test_point_point +{ + static inline void apply(std::string const& header) + { + std::string const str = header + "-"; + + test_geometry(str + "pp_01", "POINT(0 0)", "POINT(0 0)", true); + test_geometry(str + "pp_02", "POINT(0 0)", "POINT(10 0)", false); + + // points whose longitudes differ by 360 degrees + test_geometry(str + "pp_03", "POINT(0 0)", "POINT(360 0)", true); + test_geometry(str + "pp_04", "POINT(10 0)", "POINT(370 0)", true); + test_geometry(str + "pp_05", "POINT(10 0)", "POINT(-350 0)", true); + test_geometry(str + "pp_06", "POINT(180 10)", "POINT(-180 10)", true); + test_geometry(str + "pp_06a", "POINT(540 10)", "POINT(-540 10)", true); + test_geometry(str + "pp_06b", "POINT(540 370)", "POINT(-540 -350)", true); + test_geometry(str + "pp_06c", "POINT(1260 370)", "POINT(-1260 -350)", true); + test_geometry(str + "pp_06d", "POINT(2340 370)", "POINT(-2340 -350)", true); + test_geometry(str + "pp_06e", "POINT(-180 10)", "POINT(-540 10)", true); + test_geometry(str + "pp_06f", "POINT(180 10)", "POINT(-540 10)", true); + + // north & south pole + test_geometry(str + "pp_07", "POINT(0 90)", "POINT(0 90)", true); + test_geometry(str + "pp_07a", "POINT(0 450)", "POINT(10 -270)", true); + test_geometry(str + "pp_07b", "POINT(0 270)", "POINT(10 90)", false); + test_geometry(str + "pp_07c", "POINT(0 -450)", "POINT(10 90)", false); + test_geometry(str + "pp_08", "POINT(0 90)", "POINT(10 90)", true); + test_geometry(str + "pp_09", "POINT(0 90)", "POINT(0 -90)", false); + test_geometry(str + "pp_10", "POINT(0 -90)", "POINT(0 -90)", true); + test_geometry(str + "pp_11", "POINT(0 -90)", "POINT(10 -90)", true); + test_geometry(str + "pp_11a", "POINT(0 -90)", "POINT(10 90)", false); + test_geometry(str + "pp_12", "POINT(0 -90)", "POINT(0 -85)", false); + test_geometry(str + "pp_13", "POINT(0 90)", "POINT(0 85)", false); + test_geometry(str + "pp_14", "POINT(0 90)", "POINT(10 85)", false); + + // symmetric wrt prime meridian + test_geometry(str + "pp_15", "POINT(-10 45)", "POINT(10 45)", false); + test_geometry(str + "pp_16", "POINT(-170 45)", "POINT(170 45)", false); + + // other points + test_geometry(str + "pp_17", "POINT(-10 45)", "POINT(10 -45)", false); + test_geometry(str + "pp_18", "POINT(-10 -45)", "POINT(10 45)", false); + test_geometry(str + "pp_19", "POINT(10 -135)", "POINT(10 45)", false); + test_geometry(str + "pp_20", "POINT(190 135)", "POINT(10 45)", true); + test_geometry(str + "pp_21", "POINT(190 150)", "POINT(10 30)", true); + test_geometry(str + "pp_21a", "POINT(-170 150)", "POINT(10 30)", true); + test_geometry(str + "pp_22", "POINT(190 -135)", "POINT(10 -45)", true); + test_geometry(str + "pp_23", "POINT(190 -150)", "POINT(10 -30)", true); + test_geometry(str + "pp_23a", "POINT(-170 -150)", "POINT(10 -30)", true); + } +}; + + +template +void test_segment_segment(std::string const& header) +{ + typedef bgm::segment

seg; + + std::string const str = header + "-"; + + test_geometry(str + "ss_01", + "SEGMENT(10 0,180 0)", + "SEGMENT(10 0,-180 0)", + true); + test_geometry(str + "ss_02", + "SEGMENT(0 90,180 0)", + "SEGMENT(10 90,-180 0)", + true); + test_geometry(str + "ss_03", + "SEGMENT(0 90,0 -90)", + "SEGMENT(10 90,20 -90)", + true); + test_geometry(str + "ss_04", + "SEGMENT(10 80,10 -80)", + "SEGMENT(10 80,20 -80)", + false); + test_geometry(str + "ss_05", + "SEGMENT(170 10,-170 10)", + "SEGMENT(170 10,350 10)", + false); +} + + +BOOST_AUTO_TEST_CASE( equals_point_point_se ) +{ + typedef bg::cs::spherical_equatorial cs_type; + + test_point_point >::apply("se"); + test_point_point >::apply("se"); + test_point_point >::apply("se"); + + // mixed point types + test_point_point + < + bgm::point, bgm::point + >::apply("se"); + + test_point_point + < + bgm::point, bgm::point + >::apply("se"); + +#if defined(HAVE_TTMATH) + test_point_point >::apply("se"); +#endif +} + +BOOST_AUTO_TEST_CASE( equals_point_point_geo ) +{ + typedef bg::cs::geographic cs_type; + + test_point_point >::apply("geo"); + test_point_point >::apply("geo"); + test_point_point >::apply("geo"); + + // mixed point types + test_point_point + < + bgm::point, bgm::point + >::apply("se"); + + test_point_point + < + bgm::point, bgm::point + >::apply("se"); + +#if defined(HAVE_TTMATH) + test_point_point >::apply("geo"); +#endif +} + +BOOST_AUTO_TEST_CASE( equals_segment_segment_se ) +{ + typedef bg::cs::spherical_equatorial cs_type; + + test_segment_segment >("se"); + test_segment_segment >("se"); + test_segment_segment >("se"); + +#if defined(HAVE_TTMATH) + test_segment_segment >("se"); +#endif +} + +BOOST_AUTO_TEST_CASE( equals_segment_segment_geo ) +{ + typedef bg::cs::geographic cs_type; + + test_segment_segment >("geo"); + test_segment_segment >("geo"); + test_segment_segment >("geo"); + +#if defined(HAVE_TTMATH) + test_segment_segment >("geo"); +#endif +} diff --git a/test/algorithms/set_operations/intersection/intersection.cpp b/test/algorithms/set_operations/intersection/intersection.cpp index f89e70fab..2380fe582 100644 --- a/test/algorithms/set_operations/intersection/intersection.cpp +++ b/test/algorithms/set_operations/intersection/intersection.cpp @@ -434,6 +434,38 @@ void test_areal_linear() } + +template +void test_linear_box() +{ + typedef bg::model::multi_linestring multi_linestring_type; + + test_one_lp + ("case-l-b-01", + "BOX(-10 -10,10 10)", + "LINESTRING(-20 -20, 0 0,20 20)", + 1, 3, 20 * sqrt(2.0)); + + test_one_lp + ("case-l-b-02", + "BOX(-10 -10,10 10)", + "LINESTRING(-20 -20, 20 20)", + 1, 2, 20.0 * sqrt(2.0)); + + test_one_lp + ("case-l-b-02", + "BOX(-10 -10,10 10)", + "LINESTRING(-20 -20, 20 20,15 0,0 -15)", + 2, 4, 25.0 * sqrt(2.0)); + + test_one_lp + ("case-ml-b-01", + "BOX(-10 -10,10 10)", + "MULTILINESTRING((-20 -20, 20 20),(0 -15,15 0))", + 2, 4, 25.0 * sqrt(2.0)); +} + + template void test_all() { @@ -456,6 +488,8 @@ void test_all() test_areal_linear(); #endif + test_linear_box(); + // Test polygons clockwise and counter clockwise test_areal(); diff --git a/test/algorithms/transform.cpp b/test/algorithms/transform.cpp index 1746c7d23..f0414ff51 100644 --- a/test/algorithms/transform.cpp +++ b/test/algorithms/transform.cpp @@ -130,24 +130,34 @@ int test_main(int, char* []) test_all, bg::model::d2::point_xy >(1.0); test_all >, - bg::model::point > >(bg::math::d2r); + bg::model::point > >(bg::math::d2r()); test_all >, - bg::model::point > >(bg::math::r2d); + bg::model::point > >(bg::math::r2d()); test_all >, - bg::model::point > >(bg::math::d2r); + bg::model::point > >(bg::math::d2r()); test_transformations(4, 52, 1); test_transformations(4, 52, 1); - test_transformations(3 * bg::math::d2r, 51 * bg::math::d2r, 1); - test_transformations(3 * bg::math::d2r, 51 * bg::math::d2r, 1); + test_transformations + < + float, bg::radian + >(3 * bg::math::d2r(), 51 * bg::math::d2r(), 1); + + test_transformations + < + double, bg::radian + >(3 * bg::math::d2r(), 51 * bg::math::d2r(), 1); #if defined(HAVE_TTMATH) typedef bg::model::d2::point_xy PT; test_all(); test_transformations(4, 52, 1); - test_transformations(3 * bg::math::d2r, 51 * bg::math::d2r, 1); + test_transformations + < + ttmath_big, bg::radian + >(3 * bg::math::d2r(), 51 * bg::math::d2r(), 1); #endif return 0; diff --git a/test/strategies/side_of_intersection.cpp b/test/strategies/side_of_intersection.cpp index c5de7e0ea..b42270e3d 100644 --- a/test/strategies/side_of_intersection.cpp +++ b/test/strategies/side_of_intersection.cpp @@ -19,9 +19,13 @@ namespace bg = boost::geometry; int test_main(int, char* []) { - typedef bg::model::d2::point_xy point; + typedef bg::model::d2::point_xy point; typedef bg::model::segment segment; + typedef bg::strategy::side::side_of_intersection side; + + point no_fb(-99, -99); + segment a(point(20, 10), point(10, 20)); segment b1(point(11, 16), point(20, 14)); // IP with a: (14.857, 15.143) @@ -31,27 +35,70 @@ int test_main(int, char* []) segment c2(point(15, 16), point(14, 8)); segment c3(point(15, 16), point(15, 8)); - typedef bg::strategy::side::side_of_intersection side; - BOOST_CHECK_EQUAL( 1, side::apply(a, b1, c1)); - BOOST_CHECK_EQUAL(-1, side::apply(a, b1, c2)); - BOOST_CHECK_EQUAL(-1, side::apply(a, b1, c3)); + BOOST_CHECK_EQUAL( 1, side::apply(a, b1, c1, no_fb)); + BOOST_CHECK_EQUAL(-1, side::apply(a, b1, c2, no_fb)); + BOOST_CHECK_EQUAL(-1, side::apply(a, b1, c3, no_fb)); - BOOST_CHECK_EQUAL( 1, side::apply(a, b2, c1)); - BOOST_CHECK_EQUAL( 1, side::apply(a, b2, c2)); - BOOST_CHECK_EQUAL( 0, side::apply(a, b2, c3)); + BOOST_CHECK_EQUAL( 1, side::apply(a, b2, c1, no_fb)); + BOOST_CHECK_EQUAL( 1, side::apply(a, b2, c2, no_fb)); + BOOST_CHECK_EQUAL( 0, side::apply(a, b2, c3, no_fb)); - // Check internal calculation-method: - BOOST_CHECK_EQUAL(-1400, side::side_value(a, b1, c2)); - BOOST_CHECK_EQUAL( 2800, side::side_value(a, b1, c1)); + // Collinear cases + // 1: segments intersecting are collinear (with a fallback point): + { + point fb(5, 5); + segment col1(point(0, 5), point(5, 5)); + segment col2(point(5, 5), point(10, 5)); // One IP with col1 at (5,5) + segment col3(point(5, 0), point(5, 5)); + BOOST_CHECK_EQUAL( 0, side::apply(col1, col2, col3, fb)); + } + // 2: segment of side calculation collinear with one of the segments + { + point fb(5, 5); + segment col1(point(0, 5), point(10, 5)); + segment col2(point(5, 5), point(5, 12)); // IP with col1 at (5,5) + segment col3(point(5, 1), point(5, 5)); // collinear with col2 + BOOST_CHECK_EQUAL( 0, side::apply(col1, col2, col3, fb)); + } + { + point fb(5, 5); + segment col1(point(10, 5), point(0, 5)); + segment col2(point(5, 5), point(5, 12)); // IP with col1 at (5,5) + segment col3(point(5, 1), point(5, 5)); // collinear with col2 + BOOST_CHECK_EQUAL( 0, side::apply(col1, col2, col3, fb)); + } + { + point fb(5, 5); + segment col1(point(0, 5), point(10, 5)); + segment col2(point(5, 12), point(5, 5)); // IP with col1 at (5,5) + segment col3(point(5, 1), point(5, 5)); // collinear with col2 + BOOST_CHECK_EQUAL( 0, side::apply(col1, col2, col3, fb)); + } - BOOST_CHECK_EQUAL (2800, side::side_value(a, b1, c1)); - BOOST_CHECK_EQUAL(-1400, side::side_value(a, b1, c2)); - BOOST_CHECK_EQUAL(-5600, side::side_value(a, b1, c3)); + { + point fb(517248, -517236); + segment col1(point(-172408, -517236), point(862076, -517236)); + segment col2(point(517248, -862064), point(517248, -172408)); + segment col3(point(517248, -172408), point(517248, -517236)); + BOOST_CHECK_EQUAL( 0, side::apply(col1, col2, col3, fb)); + } + { + point fb(-221647, -830336); + segment col1(point(-153817, -837972), point(-222457, -830244)); + segment col2(point(-221139, -833615), point(-290654, -384388)); + segment col3(point(-255266, -814663), point(-264389, -811197)); + BOOST_CHECK_EQUAL(1, side::apply(col1, col2, col3, fb)); // Left of segment... + } - BOOST_CHECK_EQUAL(12800, side::side_value(a, b2, c1)); - BOOST_CHECK_EQUAL( 6400, side::side_value(a, b2, c2)); - BOOST_CHECK_EQUAL( 0, side::side_value(a, b2, c3)); + + { + point fb(27671131, 30809240); + segment col1(point(27671116, 30809247), point(27675474, 30807351)); + segment col2(point(27666779, 30811130), point(27671139, 30809237)); + segment col3(point(27671122, 30809244), point(27675480, 30807348)); + BOOST_CHECK_EQUAL(1, side::apply(col1, col2, col3, fb)); // Left of segment... + } // TODO: we might add a check calculating the IP, determining the side // with the normal side strategy, and verify the results are equal diff --git a/test/strategies/vincenty.cpp b/test/strategies/vincenty.cpp index b80da4bb5..e963a5549 100644 --- a/test/strategies/vincenty.cpp +++ b/test/strategies/vincenty.cpp @@ -121,24 +121,27 @@ void test_vincenty(double lon1, double lat1, double lon2, double lat2, // formula { - bg::detail::vincenty_inverse vi(lon1 * bg::math::d2r, - lat1 * bg::math::d2r, - lon2 * bg::math::d2r, - lat2 * bg::math::d2r, + double const d2r = bg::math::d2r(); + double const r2d = bg::math::r2d(); + + bg::detail::vincenty_inverse vi(lon1 * d2r, + lat1 * d2r, + lon2 * d2r, + lat2 * d2r, spheroid); calc_t dist = vi.distance(); calc_t az12 = vi.azimuth12(); calc_t az21 = vi.azimuth21(); - calc_t az12_deg = az12 * bg::math::r2d; - calc_t az21_deg = az21 * bg::math::r2d; + calc_t az12_deg = az12 * r2d; + calc_t az21_deg = az21 * r2d; BOOST_CHECK_CLOSE(dist, calc_t(expected_distance), tolerance); check_deg("az12_deg", az12_deg, calc_t(expected_azimuth_12), tolerance, error); check_deg("az21_deg", az21_deg, calc_t(expected_azimuth_21), tolerance, error); - bg::detail::vincenty_direct vd(lon1 * bg::math::d2r, - lat1 * bg::math::d2r, + bg::detail::vincenty_direct vd(lon1 * d2r, + lat1 * d2r, dist, az12, spheroid); @@ -146,9 +149,9 @@ void test_vincenty(double lon1, double lat1, double lon2, double lat2, calc_t direct_lat2 = vd.lat2(); calc_t direct_az21 = vd.azimuth21(); - calc_t direct_lon2_deg = direct_lon2 * bg::math::r2d; - calc_t direct_lat2_deg = direct_lat2 * bg::math::r2d; - calc_t direct_az21_deg = direct_az21 * bg::math::r2d; + calc_t direct_lon2_deg = direct_lon2 * r2d; + calc_t direct_lat2_deg = direct_lat2 * r2d; + calc_t direct_az21_deg = direct_az21 * r2d; check_deg("direct_lon2_deg", direct_lon2_deg, calc_t(lon2), tolerance, error); check_deg("direct_lat2_deg", direct_lat2_deg, calc_t(lat2), tolerance, error);