diff --git a/doc/release_notes.qbk b/doc/release_notes.qbk index c93dd7265..a144628c8 100644 --- a/doc/release_notes.qbk +++ b/doc/release_notes.qbk @@ -24,6 +24,8 @@ [*Improvements] +* Add distance for geographic PointLike/AnyGeometry. + [*Solved issues] * [@https://svn.boost.org/trac10/ticket/12503 12503] Validity of complex polygon @@ -33,6 +35,7 @@ * Fixes in validity of union/intersection/difference which were sometimes invalid. In most cases, results are valid now. * Fixes in validity of buffer which were sometimes invalid. In most cases, results are valid now. * Fixes in results of union/intersection/difference which could be incorrect in very complex cases +* Fixes in set and relational operations for non-cartesian coordinate systems. [/=================] [heading Boost 1.65] diff --git a/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp b/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp index 06570cf87..990081a86 100644 --- a/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp +++ b/include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp @@ -412,7 +412,7 @@ inline void buffer_point(Point const& point, Collection& collection, DistanceStrategy const& distance_strategy, PointStrategy const& point_strategy) { - collection.start_new_ring(); + collection.start_new_ring(false); std::vector range_out; point_strategy.apply(point, distance_strategy, range_out); collection.add_piece(geometry::strategy::buffer::buffered_point, range_out, false); @@ -627,7 +627,7 @@ struct buffer_inserter RobustPolicy const& robust_policy, Strategy const& strategy) // side strategy { - collection.start_new_ring(); + collection.start_new_ring(distance.negative()); geometry::strategy::buffer::result_code const code = buffer_inserter_ring::apply(ring, collection, distance, @@ -745,7 +745,7 @@ struct buffer_inserter std::size_t n = boost::size(simplified); if (n > 1) { - collection.start_new_ring(); + collection.start_new_ring(false); output_point_type first_p1; code = iterate(collection, boost::begin(simplified), boost::end(simplified), @@ -819,7 +819,13 @@ private: { for (Iterator it = begin; it != end; ++it) { - collection.start_new_ring(); + // For exterior rings, it deflates if distance is negative. + // For interior rings, it is vice versa + bool const deflate = is_interior + ? ! distance.negative() + : distance.negative(); + + collection.start_new_ring(deflate); geometry::strategy::buffer::result_code const code = policy::apply(*it, collection, distance, side_strategy, join_strategy, end_strategy, point_strategy, @@ -881,7 +887,7 @@ public: Strategy const& strategy) // side strategy { { - collection.start_new_ring(); + collection.start_new_ring(distance.negative()); geometry::strategy::buffer::result_code const code = policy::apply(exterior_ring(polygon), collection, 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 a9160b683..bc9c1ca7f 100644 --- a/include/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp +++ b/include/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp @@ -220,6 +220,7 @@ struct buffered_piece_collection bool is_flat_start; bool is_flat_end; + bool is_deflated; bool is_convex; bool is_monotonic_increasing[2]; // 0=x, 1=y bool is_monotonic_decreasing[2]; // 0=x, 1=y @@ -249,6 +250,7 @@ struct buffered_piece_collection , offsetted_count(-1) , is_flat_start(false) , is_flat_end(false) + , is_deflated(false) , is_convex(false) , robust_min_comparable_radius(0) , robust_max_comparable_radius(0) @@ -300,6 +302,8 @@ struct buffered_piece_collection piece_vector_type m_pieces; turn_vector_type m_turns; signed_size_type m_first_piece_index; + bool m_deflate; + bool m_has_deflated; buffered_ring_collection > offsetted_rings; // indexed by multi_index std::vector robust_originals; // robust representation of the original(s) @@ -338,6 +342,8 @@ struct buffered_piece_collection buffered_piece_collection(IntersectionStrategy const& intersection_strategy, RobustPolicy const& robust_policy) : m_first_piece_index(-1) + , m_deflate(false) + , m_has_deflated(false) , m_intersection_strategy(intersection_strategy) , m_side_strategy(intersection_strategy.get_side_strategy()) , m_area_strategy(intersection_strategy.template get_area_strategy()) @@ -529,6 +535,90 @@ struct buffered_piece_collection } } + struct deflate_properties + { + bool has_inflated; + std::size_t count; + + inline deflate_properties() + : has_inflated(false) + , count(0u) + {} + }; + + inline void discard_turns_for_deflate() + { + // Deflate cases should have at least 3 points PER deflated original + // to form a correct triangle + + // But if there are intersections between a deflated ring and another + // ring, it is all accepted + + // In deflate most turns are i/u by nature, but u/u is also possible + + std::map properties; + + for (typename boost::range_iterator::type it = + boost::begin(m_turns); it != boost::end(m_turns); ++it) + { + const buffer_turn_info_type& turn = *it; + if (turn.location == location_ok) + { + const buffer_turn_operation_type& op0 = turn.operations[0]; + const buffer_turn_operation_type& op1 = turn.operations[1]; + + if (! m_pieces[op0.seg_id.piece_index].is_deflated + || ! m_pieces[op1.seg_id.piece_index].is_deflated) + { + properties[op0.seg_id.multi_index].has_inflated = true; + properties[op1.seg_id.multi_index].has_inflated = true; + continue; + } + + // It is deflated, update counts + for (int i = 0; i < 2; i++) + { + const buffer_turn_operation_type& op = turn.operations[i]; + if (op.operation == detail::overlay::operation_union + || op.operation == detail::overlay::operation_continue) + { + properties[op.seg_id.multi_index].count++; + } + } + } + } + + for (typename boost::range_iterator::type it = + boost::begin(m_turns); it != boost::end(m_turns); ++it) + { + buffer_turn_info_type& turn = *it; + + if (turn.location == location_ok) + { + const buffer_turn_operation_type& op0 = turn.operations[0]; + const buffer_turn_operation_type& op1 = turn.operations[1]; + signed_size_type const multi0 = op0.seg_id.multi_index; + signed_size_type const multi1 = op1.seg_id.multi_index; + + if (multi0 == multi1) + { + const deflate_properties& prop = properties[multi0]; + if (! prop.has_inflated && prop.count < 3) + { + // Property is not inflated + // Not enough points, this might be caused by where + // detection turn-in-original failed because of numeric errors + turn.location = location_discard; + } + } + else + { + // Two different (possibly deflated) rings + } + } + } + } + template inline void check_remaining_points(DistanceStrategy const& distance_strategy) { @@ -554,7 +644,7 @@ struct buffered_piece_collection { if (deflate && turn.count_in_original <= 0) { - // For deflate: it is not in original, discard + // For deflate/negative buffers: it is not in original, discard turn.location = location_discard; } else if (! deflate && turn.count_in_original > 0) @@ -564,6 +654,12 @@ struct buffered_piece_collection } } } + + if (m_has_deflated) + { + // Either strategy was negative, or there were interior rings + discard_turns_for_deflate(); + } } inline bool assert_indices_in_robust_rings() const @@ -832,7 +928,7 @@ struct buffered_piece_collection } } - inline void start_new_ring() + inline void start_new_ring(bool deflate) { signed_size_type const n = static_cast(offsetted_rings.size()); current_segment_id.source_index = 0; @@ -844,6 +940,13 @@ struct buffered_piece_collection current_robust_ring.clear(); m_first_piece_index = static_cast(boost::size(m_pieces)); + m_deflate = deflate; + if (deflate) + { + // Pieces contain either deflated exterior rings, or inflated + // interior rings which are effectively deflated too + m_has_deflated = true; + } } inline void abort_ring() @@ -978,6 +1081,7 @@ struct buffered_piece_collection piece pc; pc.type = type; pc.index = static_cast(boost::size(m_pieces)); + pc.is_deflated = m_deflate; current_segment_id.piece_index = pc.index; diff --git a/include/boost/geometry/formulas/thomas_direct.hpp b/include/boost/geometry/formulas/thomas_direct.hpp index ef8f119b2..748960b01 100644 --- a/include/boost/geometry/formulas/thomas_direct.hpp +++ b/include/boost/geometry/formulas/thomas_direct.hpp @@ -169,9 +169,8 @@ public: { CT const sigma2 = S_sigma - sigma1; //theta2 = asin(cos(sigma2)) <=> sin_theta0 = 1 - // calculate absolute value because the geodesic is reversed - // if azimuth is pi and the result is altered below - CT const tan_theta2 = math::abs(cos(sigma2) / sin(sigma2)); + // NOTE: cos(sigma2) defines the sign of tan_theta2 + CT const tan_theta2 = cos(sigma2) / math::abs(sin(sigma2)); result.lat2 = atan(tan_theta2 / one_minus_f); } @@ -179,7 +178,7 @@ public: { result.lat2 = -result.lat2; } - } + } if (BOOST_GEOMETRY_CONDITION(CalcQuantities)) { diff --git a/include/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp b/include/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp index 73bd21ac7..fe973e3f3 100644 --- a/include/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp +++ b/include/boost/geometry/strategies/agnostic/buffer_distance_symmetric.hpp @@ -68,7 +68,7 @@ public : return negative() ? -1 : 1; } - //! Returns true if distance is negative + //! Returns true if distance is negative (aka deflate) inline bool negative() const { return m_distance < 0; diff --git a/test/algorithms/buffer/buffer_multi_polygon.cpp b/test/algorithms/buffer/buffer_multi_polygon.cpp index 961415941..ec9b276d1 100644 --- a/test/algorithms/buffer/buffer_multi_polygon.cpp +++ b/test/algorithms/buffer/buffer_multi_polygon.cpp @@ -464,8 +464,8 @@ void test_all() test_one("rt_u2", rt_u2, join_round, end_flat, 138.8001, 1.0); test_one("rt_u3", rt_u3, join_round, end_flat, 133.4526, 1.0); test_one("rt_u4", rt_u4, join_round, end_flat, 126.9268, 1.0); - test_one("rt_u5", rt_u5, join_round, end_flat, 78.4906, 1.0, ut_settings::ignore_validity()); - test_one("rt_u6", rt_u6, join_round, end_flat, 115.4461, 1.0, ut_settings::ignore_validity()); + test_one("rt_u5", rt_u5, join_round, end_flat, 78.4906, 1.0); + test_one("rt_u6", rt_u6, join_round, end_flat, 115.4461, 1.0); test_one("rt_u7", rt_u7, join_miter, end_flat, 42.6421, 1.0); test_one("rt_u7", rt_u7, join_round, end_flat, 35.6233, 1.0); diff --git a/test/algorithms/buffer/buffer_polygon.cpp b/test/algorithms/buffer/buffer_polygon.cpp index 544159d60..0a4e0ecb8 100644 --- a/test/algorithms/buffer/buffer_polygon.cpp +++ b/test/algorithms/buffer/buffer_polygon.cpp @@ -114,6 +114,10 @@ static std::string const italy_part1 static std::string const nl_part1 = "POLYGON ((471871.966884626832325 6683515.683521211147308,470695.876464393688366 6681756.129975977353752,469211.542374156008009 6679924.978601523675025,464542.357652322971262 6674631.078279769048095,463243.59315323777264 6673494.345109832473099,459502.033748187706806 6670774.304660517722368,452204.484529234410729 6666030.372161027044058,439990.287360413349234 6659313.515823394991457,434547.988775020290632 6657783.034025833010674,433867.715366783668287 6657685.311832630075514,433063.65468478546245 6657783.034025833010674,423725.285241116478574 6659758.338574352674186,422581.143514745577704 6660350.880751021206379,422333.791606202081311 6660844.461689073592424,421993.599242338153999 6662622.453031159006059,421993.599242338153999 6663363.130126810632646,422612.090333183819894 6667314.244697011075914,422241.062470370030496 6667759.324870104901493,421096.80942450783914 6668303.522556710988283,408449.802075482031796 6673245.655735924839973,401646.845354124438018 6675273.665065255947411,400842.895991614612285 6675372.844085373915732,400255.351719210040756 6675224.699206517077982,392370.258227849320974 6672455.311684883199632,391968.283546594379004 6672009.975794731639326,391875.554410762328189 6671219.573235900141299,391937.336728153284639 6670527.121663697995245,392122.906319305824582 6669933.841785561293364,392339.311409408226609 6667957.501854715868831,392308.475910458364524 6665733.178529324010015,391937.336728153284639 6665289.631341196596622,387793.802641845482867 6664449.731981083750725,385814.7647345236619 6664202.740090588107705,384268.648326894966885 6664300.5401631873101,382846.319193030707538 6664646.406163938343525,382289.610419573145919 6664943.560305980034173,377186.502322628628463 6668957.888622742146254,376352.608017096004914 6669834.728738470934331,376197.985244384559337 6670428.001423852518201,375548.658654586179182 6676313.056885857135057,375243.086652359867003 6687692.866469252854586,379228.324422756850254 6689778.692146780900657,391782.713955441897269 6694388.125634397380054,393885.427817036863416 6694487.537080916576087,395215.139134563156404 6694091.147844122722745,405171.999669074953999 6689084.665672835893929,414263.128523688821588 6684478.40333267301321,415778.298112876422238 6683439.398042757064104,416396.677884233708028 6683192.009851422160864,419025.042381353967357 6682597.809754087589681,429909.63955213711597 6681508.792763692326844,430497.183824544539675 6681656.873437026515603,440979.806314075656701 6686955.330480101518333,467325.344922157470137 6687797.546342735178769,468129.294284664443694 6687698.216345063410699,468747.785375512961764 6687450.699095219373703,469211.542374156008009 6687054.651439971290529,469489.952420631831046 6686707.835750493220985,471871.966884626832325 6683515.683521211147308))"; +static std::string const erode_triangle ="POLYGON((-262.6446228027343700 -261.9999389648437500, -261.0000000000000000 -262.0000000000000000, -261.0000915527343700 -262.5285644531250000, -262.6446228027343700 -261.9999389648437500))"; + +static std::string const forming_uu_a ="POLYGON((0 0,0 10,5 6,10 10,10 0,5 4,0 0))"; +static std::string const forming_uu_b ="POLYGON((0 0,0 10,5 6,10 10,10 0,5 4,0 0),(1.25 2.5, 4.25 5,1.25 7.5,1.25 5.5,2.25 5.0,1.25 4.5,1.25 2.5))"; // Ticket 10398, fails at next distances ( /10.0 ): // #1: 5,25,84 @@ -724,7 +728,42 @@ void test_all() join_round(49), end_flat, distance(1.0), custom_side_strategy, point_strategy, 31.1087); } +} +template +void test_deflate_special_cases() +{ + typedef bg::model::polygon polygon_type; + + bg::strategy::buffer::join_miter join_miter(5); + bg::strategy::buffer::join_round join_round(8); + bg::strategy::buffer::end_flat end_flat; + + // This case fails for because there is an IP formed at the border of the original + test_one("erode_50", erode_triangle, join_miter, end_flat, 0, 0, 0.0, -0.50); + test_one("erode_40", erode_triangle, join_miter, end_flat, 0, 0, 0.0, -0.40); + test_one("erode_60", erode_triangle, join_miter, end_flat, 0, 0, 0.0, -0.60); + + // This case generates a uu-turn in deflate, at 1.0 + test_one("forming_uu_a_95", forming_uu_a, join_round, end_flat, 1, 0, 23.0516, -0.95); + test_one("forming_uu_a_10", forming_uu_a, join_round, end_flat, 2, 0, 21.4606, -1.0); + test_one("forming_uu_a_15", forming_uu_a, join_round, end_flat, 2, 0, 8.8272, -1.5); + test_one("forming_uu_a_20", forming_uu_a, join_round, end_flat, 2, 0, 1.7588, -2.0); + test_one("forming_uu_a_21", forming_uu_a, join_round, end_flat, 2, 0, 0.9944, -2.1); + + test_one("forming_uu_b_25", forming_uu_b, join_round, end_flat, 1, 1, 38.4064, -0.25); + test_one("forming_uu_b_50", forming_uu_b, join_round, end_flat, 1, 1, 24.4551, -0.50); + test_one("forming_uu_b_95", forming_uu_b, join_round, end_flat, 1, 0, 11.5009, -0.95); + test_one("forming_uu_b_10", forming_uu_b, join_round, end_flat, 1, 0, 10.7152, -1.0); + test_one("forming_uu_b_15", forming_uu_b, join_round, end_flat, 1, 0, 4.4136, -1.5); + + test_one("forming_uu_b_25", forming_uu_b, join_round, end_flat, 1, 1, 67.7139, 0.25); + test_one("forming_uu_b_50", forming_uu_b, join_round, end_flat, 1, 1, 82.0260, 0.50); + test_one("forming_uu_b_70", forming_uu_b, join_round, end_flat, 1, 3, 93.0760, 0.70); + + // From here on a/b have the same result + test_one("forming_uu_a_10", forming_uu_a, join_round, end_flat, 1, 0, 108.0959, 1.0); + test_one("forming_uu_b_10", forming_uu_b, join_round, end_flat, 1, 0, 108.0959, 1.0); } template @@ -764,6 +803,10 @@ int test_main(int, char* []) test_all(); test_all(); + typedef bg::model::point fpoint; + test_deflate_special_cases(); + test_deflate_special_cases(); + #if ! defined(BOOST_GEOMETRY_TEST_ONLY_ONE_TYPE) test_mixed(); diff --git a/test/robustness/overlay/buffer/recursive_polygons_buffer.cpp b/test/robustness/overlay/buffer/recursive_polygons_buffer.cpp index 8816e2b59..6a6d431a9 100644 --- a/test/robustness/overlay/buffer/recursive_polygons_buffer.cpp +++ b/test/robustness/overlay/buffer/recursive_polygons_buffer.cpp @@ -87,8 +87,8 @@ bool verify(std::string const& caseid, MultiPolygon const& mp, MultiPolygon cons bool result = true; // Area of buffer must be larger than of original polygon - BOOST_AUTO(area_mp, bg::area(mp)); - BOOST_AUTO(area_buf, bg::area(buffer)); + double area_mp = bg::area(mp); + double area_buf = bg::area(buffer); if (area_buf < area_mp) { @@ -109,6 +109,16 @@ bool verify(std::string const& caseid, MultiPolygon const& mp, MultiPolygon cons } } + if (result) + { + std::string message; + if (! bg::is_valid(buffer, message)) + { + std::cout << "Buffer is not valid: " << message << std::endl; + result = false; + } + } + bool svg = settings.svg; bool wkt = settings.wkt; if (! result) @@ -204,7 +214,7 @@ bool test_buffer(MultiPolygon& result, int& index, bg::strategy::buffer::end_round end_strategy; bg::strategy::buffer::point_circle point_strategy; bg::strategy::buffer::side_straight side_strategy; - bg::strategy::buffer::join_round join_round_strategy(100); // Compatible with unit tests + bg::strategy::buffer::join_round join_round_strategy(32); // Compatible with MySQL bg::strategy::buffer::join_miter join_miter_strategy; try