From 1026f9fd366b3e67235de3f9692cdcd8e263367a Mon Sep 17 00:00:00 2001 From: Adam Wulkiewicz Date: Sun, 2 Jul 2017 19:04:12 +0200 Subject: [PATCH] [strategies] Handle degenerated segments near poles in spherical and geographic intersection strategies. --- .../strategies/geographic/intersection.hpp | 139 +++++++++---- .../strategies/spherical/intersection.hpp | 182 +++++++++++++----- 2 files changed, 236 insertions(+), 85 deletions(-) diff --git a/include/boost/geometry/strategies/geographic/intersection.hpp b/include/boost/geometry/strategies/geographic/intersection.hpp index 84acd149c0..59a40f281d 100644 --- a/include/boost/geometry/strategies/geographic/intersection.hpp +++ b/include/boost/geometry/strategies/geographic/intersection.hpp @@ -289,6 +289,8 @@ struct geographic_segments typedef typename select_calculation_type ::type calc_t; + static const calc_t c0 = 0; + // normalized spheroid srs::spheroid spheroid = normalized_spheroid(m_spheroid); @@ -325,31 +327,80 @@ struct geographic_segments // TODO: no need to call inverse formula if we know that the points are equal // distance can be set to 0 in this case and azimuth may be not calculated - bool const is_equal_a1_b1 = equals_point_point(a1, b1); - bool const is_equal_a2_b1 = equals_point_point(a2, b1); + bool is_equal_a1_b1 = equals_point_point(a1, b1); + bool is_equal_a2_b1 = equals_point_point(a2, b1); + bool degen_neq_coords = false; - inverse_result res_b1_b2 = inverse_dist_azi::apply(b1_lon, b1_lat, b2_lon, b2_lat, spheroid); - inverse_result res_b1_a1 = inverse_dist_azi::apply(b1_lon, b1_lat, a1_lon, a1_lat, spheroid); - inverse_result res_b1_a2 = inverse_dist_azi::apply(b1_lon, b1_lat, a2_lon, a2_lat, spheroid); - sides.set<0>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_b1_a1.azimuth, res_b1_b2.azimuth), - is_equal_a2_b1 ? 0 : formula::azimuth_side_value(res_b1_a2.azimuth, res_b1_b2.azimuth)); - if (sides.same<0>()) + inverse_result res_b1_b2, res_b1_a1, res_b1_a2; + if (! b_is_point) { - // Both points are at the same side of other segment, we can leave - return Policy::disjoint(); + res_b1_b2 = inverse_dist_azi::apply(b1_lon, b1_lat, b2_lon, b2_lat, spheroid); + if (math::equals(res_b1_b2.distance, c0)) + { + b_is_point = true; + degen_neq_coords = true; + } + else + { + res_b1_a1 = inverse_dist_azi::apply(b1_lon, b1_lat, a1_lon, a1_lat, spheroid); + if (math::equals(res_b1_a1.distance, c0)) + { + is_equal_a1_b1 = true; + } + res_b1_a2 = inverse_dist_azi::apply(b1_lon, b1_lat, a2_lon, a2_lat, spheroid); + if (math::equals(res_b1_a2.distance, c0)) + { + is_equal_a2_b1 = true; + } + sides.set<0>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_b1_a1.azimuth, res_b1_b2.azimuth), + is_equal_a2_b1 ? 0 : formula::azimuth_side_value(res_b1_a2.azimuth, res_b1_b2.azimuth)); + if (sides.same<0>()) + { + // Both points are at the same side of other segment, we can leave + return Policy::disjoint(); + } + } } - bool const is_equal_a1_b2 = equals_point_point(a1, b2); + bool is_equal_a1_b2 = equals_point_point(a1, b2); - inverse_result res_a1_a2 = inverse_dist_azi::apply(a1_lon, a1_lat, a2_lon, a2_lat, spheroid); - inverse_result res_a1_b1 = inverse_dist_azi::apply(a1_lon, a1_lat, b1_lon, b1_lat, spheroid); - inverse_result res_a1_b2 = inverse_dist_azi::apply(a1_lon, a1_lat, b2_lon, b2_lat, spheroid); - sides.set<1>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_a1_b1.azimuth, res_a1_a2.azimuth), - is_equal_a1_b2 ? 0 : formula::azimuth_side_value(res_a1_b2.azimuth, res_a1_a2.azimuth)); - if (sides.same<1>()) + inverse_result res_a1_a2, res_a1_b1, res_a1_b2; + if (! a_is_point) { - // Both points are at the same side of other segment, we can leave - return Policy::disjoint(); + res_a1_a2 = inverse_dist_azi::apply(a1_lon, a1_lat, a2_lon, a2_lat, spheroid); + if (math::equals(res_a1_a2.distance, c0)) + { + a_is_point = true; + degen_neq_coords = true; + } + else + { + res_a1_b1 = inverse_dist_azi::apply(a1_lon, a1_lat, b1_lon, b1_lat, spheroid); + if (math::equals(res_a1_b1.distance, c0)) + { + is_equal_a1_b1 = true; + } + res_a1_b2 = inverse_dist_azi::apply(a1_lon, a1_lat, b2_lon, b2_lat, spheroid); + if (math::equals(res_a1_b2.distance, c0)) + { + is_equal_a1_b2 = true; + } + sides.set<1>(is_equal_a1_b1 ? 0 : formula::azimuth_side_value(res_a1_b1.azimuth, res_a1_a2.azimuth), + is_equal_a1_b2 ? 0 : formula::azimuth_side_value(res_a1_b2.azimuth, res_a1_a2.azimuth)); + if (sides.same<1>()) + { + // Both points are at the same side of other segment, we can leave + return Policy::disjoint(); + } + } + } + + if(a_is_point && b_is_point) + { + return is_equal_a1_b2 + ? Policy::degenerate(a, true) + : Policy::disjoint() + ; } // NOTE: at this point the segments may still be disjoint @@ -379,11 +430,11 @@ struct geographic_segments { if (a_is_point) { - return collinear_one_degenerated(a, true, b1, b2, a1, a2, res_b1_b2, res_b1_a1, is_b_reversed); + return collinear_one_degenerated(a, true, b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, is_b_reversed, degen_neq_coords); } else if (b_is_point) { - return collinear_one_degenerated(b, false, a1, a2, b1, b2, res_a1_a2, res_a1_b1, is_a_reversed); + return collinear_one_degenerated(b, false, a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, is_a_reversed, degen_neq_coords); } else { @@ -392,16 +443,16 @@ struct geographic_segments // use shorter segment if (res_a1_a2.distance <= res_b1_b2.distance) { - calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, dist_a1_a2, dist_a1_b1); - calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b2, dist_a1_a2, dist_a1_b2); + calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_a1_a2, dist_a1_b1); + calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b2, res_a1_b1, dist_a1_a2, dist_a1_b2); dist_b1_b2 = dist_a1_b2 - dist_a1_b1; dist_b1_a1 = -dist_a1_b1; dist_b1_a2 = dist_a1_a2 - dist_a1_b1; } else { - calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a1, dist_b1_b2, dist_b1_a1); - calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a2, dist_b1_b2, dist_b1_a2); + calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a1, res_b1_a2, dist_b1_b2, dist_b1_a1); + calculate_collinear_data(b1, b2, a1, a2, res_b1_b2, res_b1_a2, res_b1_a1, dist_b1_b2, dist_b1_a2); dist_a1_a2 = dist_b1_a2 - dist_b1_a1; dist_a1_b1 = -dist_b1_a1; dist_a1_b2 = dist_b1_b2 - dist_b1_a1; @@ -549,11 +600,13 @@ struct geographic_segments Point1 const& a1, Point1 const& a2, Point2 const& b1, Point2 const& b2, ResultInverse const& res_a1_a2, - ResultInverse const& res_a1_bi, - bool is_other_reversed) + ResultInverse const& res_a1_b1, + ResultInverse const& res_a1_b2, + bool is_other_reversed, + bool degen_neq_coords) { CalcT dist_1_2, dist_1_o; - if (! calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_bi, dist_1_2, dist_1_o)) + if (! calculate_collinear_data(a1, a2, b1, b2, res_a1_a2, res_a1_b1, res_a1_b2, dist_1_2, dist_1_o, degen_neq_coords)) { return Policy::disjoint(); } @@ -574,13 +627,16 @@ struct geographic_segments static inline bool calculate_collinear_data(Point1 const& a1, Point1 const& a2, // in Point2 const& b1, Point2 const& b2, // in ResultInverse const& res_a1_a2, // in - ResultInverse const& res_a1_bi, // in - CalcT& dist_a1_a2, CalcT& dist_a1_bi) // out + ResultInverse const& res_a1_b1, // in + ResultInverse const& res_a1_b2, // in + CalcT& dist_a1_a2, // out + CalcT& dist_a1_bi, // out + bool degen_neq_coords = false) // in { dist_a1_a2 = res_a1_a2.distance; - dist_a1_bi = res_a1_bi.distance; - if (! same_direction(res_a1_bi.azimuth, res_a1_a2.azimuth)) + dist_a1_bi = res_a1_b1.distance; + if (! same_direction(res_a1_b1.azimuth, res_a1_a2.azimuth)) { dist_a1_bi = -dist_a1_bi; } @@ -598,6 +654,22 @@ struct geographic_segments return true; } + // check the other endpoint of a very short segment near the pole + if (degen_neq_coords) + { + static CalcT const c0 = 0; + if (math::equals(res_a1_b2.distance, c0)) + { + dist_a1_bi = 0; + return true; + } + else if (math::equals(dist_a1_a2 - res_a1_b2.distance, c0)) + { + dist_a1_bi = dist_a1_a2; + return true; + } + } + // or i1 is on b return segment_ratio(dist_a1_bi, dist_a1_a2).on_segment(); } @@ -825,8 +897,9 @@ struct geographic_segments static inline bool is_endpoint_equal(CalcT const& dist, P1 const& ai, P2 const& b1, P2 const& b2) { + static CalcT const c0 = 0; using geometry::detail::equals::equals_point_point; - return is_near(dist) && (equals_point_point(ai, b1) || equals_point_point(ai, b2)); + return is_near(dist) && (equals_point_point(ai, b1) || equals_point_point(ai, b2) || math::equals(dist, c0)); } template diff --git a/include/boost/geometry/strategies/spherical/intersection.hpp b/include/boost/geometry/strategies/spherical/intersection.hpp index 035961c0a8..44b1cc62bf 100644 --- a/include/boost/geometry/strategies/spherical/intersection.hpp +++ b/include/boost/geometry/strategies/spherical/intersection.hpp @@ -302,36 +302,64 @@ struct ecef_segments vec3d_t const b1v = calc_policy.template to_cart3d(b1); vec3d_t const b2v = calc_policy.template to_cart3d(b2); + bool degen_neq_coords = false; side_info sides; typename CalcPolicy::template plane plane2 = calc_policy.get_plane(b1v, b2v); - // not normalized normals, the same as in side strategy - sides.set<0>(plane2.side_value(a1v), plane2.side_value(a2v)); - if (sides.same<0>()) + calc_t dist_b1_b2 = 0; + if (! b_is_point) { - // Both points are at same side of other segment, we can leave - return Policy::disjoint(); + calculate_dist(b1v, b2v, plane2, dist_b1_b2); + if (math::equals(dist_b1_b2, c0)) + { + degen_neq_coords = true; + b_is_point = true; + dist_b1_b2 = 0; + } + else + { + // not normalized normals, the same as in side strategy + sides.set<0>(plane2.side_value(a1v), plane2.side_value(a2v)); + if (sides.same<0>()) + { + // Both points are at same side of other segment, we can leave + return Policy::disjoint(); + } + } } typename CalcPolicy::template plane plane1 = calc_policy.get_plane(a1v, a2v); - // not normalized normals, the same as in side strategy - sides.set<1>(plane1.side_value(b1v), plane1.side_value(b2v)); - if (sides.same<1>()) + calc_t dist_a1_a2 = 0; + if (! a_is_point) { - // Both points are at same side of other segment, we can leave - return Policy::disjoint(); + calculate_dist(a1v, a2v, plane1, dist_a1_a2); + if (math::equals(dist_a1_a2, c0)) + { + degen_neq_coords = true; + a_is_point = true; + dist_a1_a2 = 0; + } + else + { + // not normalized normals, the same as in side strategy + sides.set<1>(plane1.side_value(b1v), plane1.side_value(b2v)); + if (sides.same<1>()) + { + // Both points are at same side of other segment, we can leave + return Policy::disjoint(); + } + } } // NOTE: at this point the segments may still be disjoint - calc_t len1, len2; - + calc_t len1 = 0; // point or opposite sides of a sphere/spheroid, assume point - if (! detail::vec_normalize(plane1.normal, len1)) + if (! a_is_point && ! detail::vec_normalize(plane1.normal, len1)) { a_is_point = true; if (sides.get<0, 0>() == 0 || sides.get<0, 1>() == 0) @@ -340,7 +368,8 @@ struct ecef_segments } } - if (! detail::vec_normalize(plane2.normal, len2)) + calc_t len2 = 0; + if (! b_is_point && ! detail::vec_normalize(plane2.normal, len2)) { b_is_point = true; if (sides.get<1, 0>() == 0 || sides.get<1, 1>() == 0) @@ -404,30 +433,32 @@ struct ecef_segments { if (a_is_point) { - return collinear_one_degenerated(a, true, b1, b2, a1, a2, b1v, b2v, plane2, a1v); + return collinear_one_degenerated(a, true, b1, b2, a1, a2, b1v, b2v, + plane2, a1v, a2v, dist_b1_b2, degen_neq_coords); } else if (b_is_point) { // b2 used to be consistent with (degenerated) checks above (is it needed?) - return collinear_one_degenerated(b, false, a1, a2, b1, b2, a1v, a2v, plane1, b1v); + return collinear_one_degenerated(b, false, a1, a2, b1, b2, a1v, a2v, + plane1, b1v, b2v, dist_a1_a2, degen_neq_coords); } else { - calc_t dist_a1_a2, dist_a1_b1, dist_a1_b2; - calc_t dist_b1_b2, dist_b1_a1, dist_b1_a2; + calc_t dist_a1_b1, dist_a1_b2; + calc_t dist_b1_a1, dist_b1_a2; // use shorter segment if (len1 <= len2) { - calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, plane1, b1v, dist_a1_a2, dist_a1_b1); - calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, plane1, b2v, dist_a1_a2, dist_a1_b2); + calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, plane1, b1v, b2v, dist_a1_a2, dist_a1_b1); + calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, plane1, b2v, b1v, dist_a1_a2, dist_a1_b2); dist_b1_b2 = dist_a1_b2 - dist_a1_b1; dist_b1_a1 = -dist_a1_b1; dist_b1_a2 = dist_a1_a2 - dist_a1_b1; } else { - calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, plane2, a1v, dist_b1_b2, dist_b1_a1); - calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, plane2, a2v, dist_b1_b2, dist_b1_a2); + calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, plane2, a1v, a2v, dist_b1_b2, dist_b1_a1); + calculate_collinear_data(b1, b2, a1, a2, b1v, b2v, plane2, a2v, a1v, dist_b1_b2, dist_b1_a2); dist_a1_a2 = dist_b1_a2 - dist_b1_a1; dist_a1_b1 = -dist_b1_a1; dist_a1_b2 = dist_b1_b2 - dist_b1_a1; @@ -487,10 +518,11 @@ struct ecef_segments vec3d_t i1; intersection_point_flag ip_flag; - calc_t dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1; + calc_t dist_a1_i1, dist_b1_i1; if (calculate_ip_data(a1, a2, b1, b2, a1v, a2v, b1v, b2v, - plane1, plane2, calc_policy, sides, - i1, dist_a1_a2, dist_a1_i1, dist_b1_b2, dist_b1_i1, ip_flag)) + plane1, plane2, calc_policy, + sides, dist_a1_a2, dist_b1_b2, + i1, dist_a1_i1, dist_b1_i1, ip_flag)) { // intersects segment_intersection_info @@ -520,12 +552,14 @@ struct ecef_segments collinear_one_degenerated(Segment const& segment, bool degenerated_a, Point1 const& a1, Point1 const& a2, Point2 const& b1, Point2 const& b2, - Vec3d const& v1, Vec3d const& v2, + Vec3d const& a1v, Vec3d const& a2v, Plane const& plane, - Vec3d const& vother) + Vec3d const& b1v, Vec3d const& b2v, + CalcT const& dist_1_2, + bool degen_neq_coords) { - CalcT dist_1_2, dist_1_o; - return ! calculate_collinear_data(a1, a2, b1, b2, v1, v2, plane, vother, dist_1_2, dist_1_o) + CalcT dist_1_o; + return ! calculate_collinear_data(a1, a2, b1, b2, a1v, a2v, plane, b1v, b2v, dist_1_2, dist_1_o, degen_neq_coords) ? Policy::disjoint() : Policy::one_degenerate(segment, segment_ratio(dist_1_o, dist_1_2), degenerated_a); } @@ -536,11 +570,14 @@ struct ecef_segments Vec3d const& a1v, // in Vec3d const& a2v, // in Plane const& plane1, // in - Vec3d const& b1v_or_b2v, // in - CalcT& dist_a1_a2, CalcT& dist_a1_i1) // out + Vec3d const& b1v, // in + Vec3d const& b2v, // in + CalcT const& dist_a1_a2, // in + CalcT& dist_a1_i1, // out + bool degen_neq_coords = false) // in { // calculate dist_a1_a2 and dist_a1_i1 - calculate_dists(a1v, a2v, plane1, b1v_or_b2v, dist_a1_a2, dist_a1_i1); + calculate_dist(a1v, a2v, plane1, b1v, dist_a1_i1); // if i1 is close to a1 and b1 or b2 is equal to a1 if (is_endpoint_equal(dist_a1_i1, a1, b1, b2)) @@ -555,6 +592,26 @@ struct ecef_segments return true; } + // check the other endpoint of a very short segment near the pole + if (degen_neq_coords) + { + static CalcT const c0 = 0; + + CalcT dist_a1_i2 = 0; + calculate_dist(a1v, a2v, plane1, b2v, dist_a1_i2); + + if (math::equals(dist_a1_i2, c0)) + { + dist_a1_i1 = 0; + return true; + } + else if (math::equals(dist_a1_a2 - dist_a1_i2, c0)) + { + dist_a1_i1 = dist_a1_a2; + return true; + } + } + // or i1 is on b return segment_ratio(dist_a1_i1, dist_a1_a2).on_segment(); } @@ -568,15 +625,17 @@ struct ecef_segments Plane const& plane2, // in CalcPolicy const& calc_policy, // in side_info const& sides, // in - Vec3d & ip, // out - CalcT& dist_a1_a2, CalcT& dist_a1_ip, // out - CalcT& dist_b1_b2, CalcT& dist_b1_ip, // out - intersection_point_flag& ip_flag) // out + CalcT const& dist_a1_a2, // in + CalcT const& dist_b1_b2, // in + Vec3d & ip, // out + CalcT& dist_a1_ip, // out + CalcT& dist_b1_ip, // out + intersection_point_flag& ip_flag) // out { Vec3d ip1, ip2; calc_policy.intersection_points(plane1, plane2, ip1, ip2); - calculate_dists(a1v, a2v, plane1, ip1, dist_a1_a2, dist_a1_ip); + calculate_dist(a1v, a2v, plane1, ip1, dist_a1_ip); ip = ip1; // choose the opposite side of the globe if the distance is shorter @@ -602,7 +661,7 @@ struct ecef_segments return false; } - calculate_dists(b1v, b2v, plane2, ip, dist_b1_b2, dist_b1_ip); + calculate_dist(b1v, b2v, plane2, ip, dist_b1_ip); bool is_on_b = false, is_near_b1 = false, is_near_b2 = false; if (! is_potentially_crossing(dist_b1_b2, dist_b1_ip, is_on_b, is_near_b1, is_near_b2)) @@ -701,20 +760,26 @@ struct ecef_segments } template - static inline void calculate_dists(Vec3d const& a1v, // in - Vec3d const& a2v, // in - Plane const& plane1, // in - Vec3d const& i1, // in - CalcT& dist_a1_a2, // out - CalcT& dist_a1_i1) // out + static inline void calculate_dist(Vec3d const& a1v, // in + Vec3d const& a2v, // in + Plane const& plane1, // in + CalcT& dist_a1_a2) // out { - //CalcT const c0 = 0; - CalcT const c1 = 1; - CalcT const c2 = 2; - CalcT const c4 = 4; - - CalcT cos_a1_a2 = plane1.cos_angle_between(a1v, a2v); + static CalcT const c1 = 1; + CalcT const cos_a1_a2 = plane1.cos_angle_between(a1v, a2v); dist_a1_a2 = -cos_a1_a2 + c1; // [1, -1] -> [0, 2] representing [0, pi] + } + + template + static inline void calculate_dist(Vec3d const& a1v, // in + Vec3d const& /*a2v*/, // in + Plane const& plane1, // in + Vec3d const& i1, // in + CalcT& dist_a1_i1) // out + { + static CalcT const c1 = 1; + static CalcT const c2 = 2; + static CalcT const c4 = 4; bool is_forward = true; CalcT cos_a1_i1 = plane1.cos_angle_between(a1v, i1, is_forward); @@ -728,7 +793,19 @@ struct ecef_segments dist_a1_i1 += c4; // += 2pi } } - + /* + template + static inline void calculate_dists(Vec3d const& a1v, // in + Vec3d const& a2v, // in + Plane const& plane1, // in + Vec3d const& i1, // in + CalcT& dist_a1_a2, // out + CalcT& dist_a1_i1) // out + { + calculate_dist(a1v, a2v, plane1, dist_a1_a2); + calculate_dist(a1v, a2v, plane1, i1, dist_a1_i1); + } + */ // the dist of the ip on the other side of the sphere template static inline CalcT dist_of_i2(CalcT const& dist_a1_i1) @@ -769,8 +846,9 @@ struct ecef_segments static inline bool is_endpoint_equal(CalcT const& dist, P1 const& ai, P2 const& b1, P2 const& b2) { + static CalcT const c0 = 0; using geometry::detail::equals::equals_point_point; - return is_near(dist) && (equals_point_point(ai, b1) || equals_point_point(ai, b2)); + return is_near(dist) && (equals_point_point(ai, b1) || equals_point_point(ai, b2) || math::equals(dist, c0)); } template