boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-02-01 08:32:10 +00:00
Code Issues Projects Releases Wiki Activity

[strategies] Handle degenerated segments near poles in spherical and geographic intersection strategies.

Browse Source
This commit is contained in:
Adam Wulkiewicz
2017-07-02 19:04:12 +02:00
parent c194e1221e
commit 1026f9fd36
2 changed files with 236 additions and 85 deletions
Download Patch File Download Diff File Expand all files Collapse all files

139
include/boost/geometry/strategies/geographic/intersection.hpp
View File

@@ -289,6 +289,8 @@ private:
typedef typename select_calculation_type
<Segment1, Segment2, CalculationType>::type calc_t;
static const calc_t c0 = 0;
// normalized spheroid
srs::spheroid<calc_t> spheroid = normalized_spheroid<calc_t>(m_spheroid);
@@ -325,31 +327,80 @@ private:
// 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 @@ private:
{
if (a_is_point)
{
return collinear_one_degenerated<Policy, calc_t>(a, true, b1, b2, a1, a2, res_b1_b2, res_b1_a1, is_b_reversed);
return collinear_one_degenerated<Policy, calc_t>(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<Policy, calc_t>(b, false, a1, a2, b1, b2, res_a1_a2, res_a1_b1, is_a_reversed);
return collinear_one_degenerated<Policy, calc_t>(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 @@ private:
// 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 @@ private:
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 @@ private:
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 @@ private:
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<CalcT>(dist_a1_bi, dist_a1_a2).on_segment();
}
@@ -825,8 +897,9 @@ private:
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 <typename CalcT>

182
include/boost/geometry/strategies/spherical/intersection.hpp
View File

@@ -302,36 +302,64 @@ struct ecef_segments
vec3d_t const b1v = calc_policy.template to_cart3d<vec3d_t>(b1);
vec3d_t const b2v = calc_policy.template to_cart3d<vec3d_t>(b2);
bool degen_neq_coords = false;
side_info sides;
typename CalcPolicy::template plane<vec3d_t>
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<vec3d_t>
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<Policy, calc_t>(a, true, b1, b2, a1, a2, b1v, b2v, plane2, a1v);
return collinear_one_degenerated<Policy, calc_t>(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<Policy, calc_t>(b, false, a1, a2, b1, b2, a1v, a2v, plane1, b1v);
return collinear_one_degenerated<Policy, calc_t>(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 @@ private:
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<CalcT>(dist_1_o, dist_1_2), degenerated_a);
}
@@ -536,11 +570,14 @@ private:
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 @@ private:
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<CalcT>(dist_a1_i1, dist_a1_a2).on_segment();
}
@@ -568,15 +625,17 @@ private:
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 @@ private:
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 @@ private:
}
template <typename Vec3d, typename Plane, typename CalcT>
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 <typename Vec3d, typename Plane, typename CalcT>
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 @@ private:
dist_a1_i1 += c4; // += 2pi
}
}
/*
template <typename Vec3d, typename Plane, typename CalcT>
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 <typename CalcT>
static inline CalcT dist_of_i2(CalcT const& dist_a1_i1)
@@ -769,8 +846,9 @@ private:
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 <typename CalcT>