Skip to content

Commit

Permalink
[within] Fix special case of segment pole endpoint
Browse files Browse the repository at this point in the history
  • Loading branch information
@vissarion
vissarion committed Jul 14, 2023
1 parent da3de29 commit 1a04c78
Show file tree
Hide file tree
Showing 2 changed files with 156 additions and 11 deletions.
38 changes: 34 additions & 4 deletions include/boost/geometry/strategies/spherical/point_in_poly_winding.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -273,11 +273,9 @@ class spherical_winding_base
calc_t const anti_p_lon = p_lon + (p_lon <= c0 ? pi : -pi);

eq1 = eq1_strict // lon strictly equal to s1
|| (eq1_anti = longitudes_equal<units_t>(s1_lon, anti_p_lon)) // anti-lon strictly equal to s1
|| math::equals(math::abs(s1_lat), half_pi); // s1 is pole
|| (eq1_anti = longitudes_equal<units_t>(s1_lon, anti_p_lon)); // anti-lon strictly equal to s1
eq2 = eq2_strict // lon strictly equal to s2
|| (eq2_anti = longitudes_equal<units_t>(s2_lon, anti_p_lon)) // anti-lon strictly equal to s2
|| math::equals(math::abs(s2_lat), half_pi); // s2 is pole
|| (eq2_anti = longitudes_equal<units_t>(s2_lon, anti_p_lon)); // anti-lon strictly equal to s2

// segment overlapping pole
calc_t const s_lon_diff = math::longitude_distance_signed<units_t>(s1_lon, s2_lon);
Expand All @@ -293,6 +291,28 @@ class spherical_winding_base
}
}

// check whether point is on a segment with a pole endpoint
if (math::longitude_distance_signed<units_t>(s2_lon, p_lon) == c0)
{
bool const s1_north = math::equals(get<1>(seg1), half_pi);
bool const s1_south = math::equals(get<1>(seg1), -half_pi);
if (s1_north || s1_south)
{
state.m_touches = s1_south ? s2_lat > p_lat : s2_lat < p_lat;
return state.m_touches;
}
}
if (math::longitude_distance_signed<units_t>(s1_lon, p_lon) == c0)
{
bool const s2_north = math::equals(get<1>(seg2), half_pi);
bool const s2_south = math::equals(get<1>(seg2), -half_pi);
if (s2_north || s2_south)
{
state.m_touches = s2_south ? s1_lat > p_lat : s1_lat < p_lat;
return state.m_touches;
}
}

// Both equal p -> segment vertical
// The only thing which has to be done is check if point is ON segment
if (eq1 && eq2)
Expand Down Expand Up @@ -361,7 +381,17 @@ class spherical_winding_base
// If needed (eq1 && eq2 ? 0) could be returned

calc_t const c0 = 0;
calc_t const c2 = 2;
calc_t const pi = constants::half_period();
calc_t const half_pi = pi / c2;

bool const s1_is_pole = math::equals(std::abs(get<1>(seg1)), half_pi);
bool const s2_is_pole = math::equals(std::abs(get<1>(seg2)), half_pi);

if (s1_is_pole && s2_is_pole)
{
return count_info(0, false);
}

calc_t const p = get<0>(point);
calc_t const s1 = get<0>(seg1);
Expand Down
129 changes: 122 additions & 7 deletions test/algorithms/within/within_pointlike_geometry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@ void test_p_l()

test_geometry<P, mls>("POINT(0 0)", "MULTILINESTRING((0 0,1 1,2 2),(0 0,0 1))", true);
test_geometry<P, mls>("POINT(0 0)", "MULTILINESTRING((0 0,1 1,2 2),(0 0,0 1),(0 0,1 0))", false);

test_geometry<P, mls>("POINT(1 1)", "MULTILINESTRING((0 0, 1 1),(1 1, 2 2))", true);
test_geometry<P, mls>("POINT(1 1)", "MULTILINESTRING((0 0, 1 1),(2 2, 3 3))", false);

Expand Down Expand Up @@ -289,12 +289,63 @@ void test_spherical_geographic()
// Points outside but on the same level as segment
BOOST_CHECK_EQUAL(bg::covered_by(Point(-90, 75), poly_n), false);
}*/
// Segment endpoints on pole with arbitrary longitudes
// Segment endpoints on North pole with arbitrary longitudes
{
bg::model::polygon<Point> poly_n3;
bg::read_wkt("POLYGON((45 90,45 80,0 80,45 90))", poly_n3);
BOOST_CHECK_EQUAL(bg::covered_by(Point(0, 85), poly_n3, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(45, 85), poly_n3, ws), true);
bg::model::polygon<Point> poly_n4;
bg::read_wkt("POLYGON((45 90,45 80,-10 80,45 90))", poly_n4);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-10, 85), poly_n4, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(0, 85), poly_n4, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(45, 85), poly_n4, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, 85), poly_n4, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, 85), poly_n4, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-5, 85), poly_n4, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, 70), poly_n4, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, 70), poly_n4, ws), false);

// the same polygon but with two points representing the pole
bg::model::polygon<Point> poly_n4b;
bg::read_wkt("POLYGON((45 90,45 80,-10 80,60 90,45 90))", poly_n4b);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-10, 85), poly_n4b, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(0, 85), poly_n4b, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(45, 85), poly_n4b, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, 85), poly_n4b, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, 85), poly_n4b, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-5, 85), poly_n4b, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, 70), poly_n4b, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, 70), poly_n4b, ws), false);

bg::model::polygon<Point> poly_n5;
bg::read_wkt("POLYGON((0 90,-10 80,45 80,0 90))", poly_n5);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-10, 85), poly_n5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(0, 85), poly_n5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-1, 85), poly_n5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(1, 85), poly_n5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(45, 85), poly_n5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, 85), poly_n5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, 85), poly_n5, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-5, 85), poly_n5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, 70), poly_n5, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, 70), poly_n5, ws), false);

bg::model::polygon<Point> poly_n_4edges;
bg::read_wkt("POLYGON((0 90,-10 70,5 60,20 80,0 90))", poly_n_4edges);
BOOST_CHECK_EQUAL(bg::covered_by(Point(3, 89), poly_n_4edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-1, 87), poly_n_4edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-10, 86), poly_n_4edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(15, 84), poly_n_4edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-1, 61), poly_n_4edges, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(15, 81), poly_n_4edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(7, 50), poly_n_4edges, ws), false);

bg::model::polygon<Point> poly_n_5edges;
bg::read_wkt("POLYGON((0 90,-10 70,5 60,10 85,20 80,0 90))", poly_n_5edges);
BOOST_CHECK_EQUAL(bg::covered_by(Point(3, 89), poly_n_5edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-1, 87), poly_n_5edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-10, 86), poly_n_5edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(15, 84), poly_n_5edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-1, 61), poly_n_5edges, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(15, 81), poly_n_5edges, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(7, 50), poly_n_5edges, ws), false);
}
// Segment going through pole
{
Expand All @@ -308,12 +359,76 @@ void test_spherical_geographic()
BOOST_CHECK_EQUAL(bg::covered_by(Point(0, -90), poly_s1, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(45, -90), poly_s1, ws), true);
}
// Segment endpoints on pole with arbitrary longitudes
// Segment endpoints on South pole with arbitrary longitudes
{
bg::model::polygon<Point> poly_s2;
bg::read_wkt("POLYGON((45 -90,0 -80,45 -80,45 -90))", poly_s2);
BOOST_CHECK_EQUAL(bg::covered_by(Point(0, -85), poly_s2, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(45, -85), poly_s2, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, -85), poly_s2, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, -85), poly_s2, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-5, -85), poly_s2, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, -70), poly_s2, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, -70), poly_s2, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-5, -70), poly_s2, ws), false);

bg::model::polygon<Point> poly_s3;
bg::read_wkt("POLYGON((45 -90,-10 -80,45 -80,45 -90))", poly_s3);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-10, -85), poly_s3, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(0, -85), poly_s3, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-1, -85), poly_s3, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(1, -85), poly_s3, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(45, -85), poly_s3, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, -85), poly_s3, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, -85), poly_s3, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-5, -85), poly_s3, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, -70), poly_s3, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, -70), poly_s3, ws), false);

bg::model::polygon<Point> poly_s5;
bg::read_wkt("POLYGON((0 -90,-10 -80,45 -80,0 -90))", poly_s5);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-10, -85), poly_s5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(0, -85), poly_s5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-1, -85), poly_s5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(1, -85), poly_s5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(45, -85), poly_s5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, -85), poly_s5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, -85), poly_s5, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-5, -85), poly_s5, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, -70), poly_s5, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, -70), poly_s5, ws), false);

bg::model::polygon<Point> poly_s4;
bg::read_wkt("POLYGON((0 -89,-10 -80,45 -80,0 -89))", poly_s4);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-10, -85), poly_s4, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(0, -85), poly_s4, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(45, -85), poly_s4, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, -85), poly_s4, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, -85), poly_s4, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-5, -85), poly_s4, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(30, -71), poly_s4, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(50, -70), poly_s4, ws), false);

//more complex examples
bg::model::polygon<Point> poly_s_complex_4edges;
bg::read_wkt("POLYGON((0 -90,-10 -70,5 -60,20 -80,0 -90))", poly_s_complex_4edges);
BOOST_CHECK_EQUAL(bg::covered_by(Point(3, -89), poly_s_complex_4edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-1, -87), poly_s_complex_4edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-10, -86), poly_s_complex_4edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(15, -84), poly_s_complex_4edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-1, -61), poly_s_complex_4edges, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(15, -81), poly_s_complex_4edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(7, -50), poly_s_complex_4edges, ws), false);

bg::model::polygon<Point> poly_s_complex_5edges;
bg::read_wkt("POLYGON((0 -90,-10 -70,5 -60,10 -85,20 -80,0 -90))", poly_s_complex_5edges);
BOOST_CHECK_EQUAL(bg::covered_by(Point(3, -89), poly_s_complex_5edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-1, -87), poly_s_complex_5edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-10, -86), poly_s_complex_5edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(15, -84), poly_s_complex_5edges, ws), true);
BOOST_CHECK_EQUAL(bg::covered_by(Point(-1, -61), poly_s_complex_5edges, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(15, -81), poly_s_complex_5edges, ws), false);
BOOST_CHECK_EQUAL(bg::covered_by(Point(7, -50), poly_s_complex_5edges, ws), false);
}
// Polygon covering nearly half of the globe but no poles
{
Expand Down

0 comments on commit 1a04c78

Please sign in to comment.