Commit eff0a116 authored by Jason Davies's avatar Jason Davies
Browse files

Fix point-in-polygon for multiple polar rings.

Instead of detecting if any single polygon ring winds around a pole, we
consider the cumulative winding of all polygon rings together.  This is
consistent with the area calculation, which considers the cumulative
area total of all rings.

This fixes #1521: an issue with the Hammer Retroazimuthal projection,
which uses such a polygon with two rings, covering most of the globe.

In addition, drop the special handling of points at the south pole,
which might have been there to pass an incorrect test: a CCW triangle
touching the south pole, which was probably incorrectly thought to be
clockwise.  This fixes an issue with a “stripe” polygon rotated so that
a point is at the south pole, mentioned in #1453.
parent 5f8056c1
......@@ -2821,7 +2821,7 @@ d3 = function() {
return ((a = a.point)[0] < 0 ? a[1] - π / 2 - ε : π / 2 - a[1]) - ((b = b.point)[0] < 0 ? b[1] - π / 2 - ε : π / 2 - b[1]);
}
function d3_geo_pointInPolygon(point, polygon) {
var meridian = point[0], parallel = point[1], meridianNormal = [ Math.sin(meridian), -Math.cos(meridian), 0 ], polarAngle = 0, polar = false, southPole = false, winding = 0;
var meridian = point[0], parallel = point[1], meridianNormal = [ Math.sin(meridian), -Math.cos(meridian), 0 ], polarAngle = 0, winding = 0;
d3_geo_areaRingSum.reset();
for (var i = 0, n = polygon.length; i < n; ++i) {
var ring = polygon[i], m = ring.length;
......@@ -2832,7 +2832,6 @@ d3 = function() {
point = ring[j];
var λ = point[0], φ = point[1] / 2 + π / 4, sinφ = Math.sin(φ), cosφ = Math.cos(φ), = λ - λ0, antimeridian = Math.abs() > π, k = sinφ0 * sinφ;
d3_geo_areaRingSum.add(Math.atan2(k * Math.sin(), cosφ0 * cosφ + k * Math.cos()));
if (Math.abs(φ) < ε) southPole = true;
polarAngle += antimeridian ? + ( >= 0 ? 2 : -2) * π : ;
if (antimeridian ^ λ0 >= meridian ^ λ >= meridian) {
var arc = d3_geo_cartesianCross(d3_geo_cartesian(point0), d3_geo_cartesian(point));
......@@ -2847,9 +2846,8 @@ d3 = function() {
if (!j++) break;
λ0 = λ, sinφ0 = sinφ, cosφ0 = cosφ, point0 = point;
}
if (Math.abs(polarAngle) > ε) polar = true;
}
return (!southPole && !polar && d3_geo_areaRingSum < 0 || polarAngle < -ε) ^ winding & 1;
return (polarAngle < -ε || polarAngle < ε && d3_geo_areaRingSum < 0) ^ winding & 1;
}
var d3_geo_clipAntimeridian = d3_geo_clip(d3_true, d3_geo_clipAntimeridianLine, d3_geo_clipAntimeridianInterpolate, d3_geo_clipAntimeridianPolygonContains);
function d3_geo_clipAntimeridianLine(listener) {
......
This diff is collapsed.
......@@ -8,8 +8,6 @@ function d3_geo_pointInPolygon(point, polygon) {
parallel = point[1],
meridianNormal = [Math.sin(meridian), -Math.cos(meridian), 0],
polarAngle = 0,
polar = false,
southPole = false,
winding = 0;
d3_geo_areaRingSum.reset();
......@@ -36,7 +34,6 @@ function d3_geo_pointInPolygon(point, polygon) {
k = sinφ0 * sinφ;
d3_geo_areaRingSum.add(Math.atan2(k * Math.sin(), cosφ0 * cosφ + k * Math.cos()));
if (Math.abs(φ) < ε) southPole = true;
polarAngle += antimeridian ? + ( >= 0 ? 2 : -2) * π : ;
// Are the longitudes either side of the point's meridian, and are the
......@@ -54,18 +51,18 @@ function d3_geo_pointInPolygon(point, polygon) {
if (!j++) break;
λ0 = λ, sinφ0 = sinφ, cosφ0 = cosφ, point0 = point;
}
if (Math.abs(polarAngle) > ε) polar = true;
}
// First, determine whether the South pole is inside or outside:
//
// It is inside if:
// * the polygon doesn't wind around it, and its area is negative (counter-clockwise).
// * otherwise, if the polygon winds around it in a clockwise direction.
// * the polygon winds around it in a clockwise direction.
// * the polygon does not (cumulatively) wind around it, but has a negative
// (counter-clockwise) area.
//
// Second, count the (signed) number of times a segment crosses a meridian
// from the point to the South pole. If it is zero, then the point is the
// same side as the South pole.
return (!southPole && !polar && d3_geo_areaRingSum < 0 || polarAngle < -ε) ^ (winding & 1);
return (polarAngle < -ε || polarAngle < ε && d3_geo_areaRingSum < 0) ^ (winding & 1);
}
......@@ -122,6 +122,36 @@ suite.addBatch({
"inside": function(pointInPolygon) {
assert.ok(pointInPolygon([0, 20]));
}
},
"narrow equatorial hole": {
topic: function(pointInPolygon) {
var circle = _.geo.circle().origin([0, -90]);
return pointInPolygon([
circle.angle(90 - .01)().coordinates[0],
circle.angle(90 + .01)().coordinates[0].reverse()
]);
},
"outside": function(pointInPolygon) {
assert.ok(!pointInPolygon([0, 0]));
},
"inside": function(pointInPolygon) {
assert.ok(pointInPolygon([0, -90]));
}
},
"narrow equatorial strip": {
topic: function(pointInPolygon) {
var circle = _.geo.circle().origin([0, -90]);
return pointInPolygon([
circle.angle(90 + .01)().coordinates[0],
circle.angle(90 - .01)().coordinates[0].reverse()
]);
},
"outside": function(pointInPolygon) {
assert.ok(!pointInPolygon([0, -90]));
},
"inside": function(pointInPolygon) {
assert.ok(pointInPolygon([0, 0]));
}
}
},
"ring": {
......@@ -317,15 +347,15 @@ suite.addBatch({
return pointInPolygon([[[180, -90], [-135, 0], [135, 0], [180, -90]]]);
},
"inside": function(pointInPolygon) {
assert.ok(pointInPolygon([180, 0]));
assert.ok(pointInPolygon([150, 0]));
assert.ok(pointInPolygon([180, -30]));
assert.ok(pointInPolygon([150, -80]));
assert.ok(pointInPolygon([0, 0]));
assert.ok(pointInPolygon([180, 1]));
assert.ok(pointInPolygon([-90, -80]));
},
"outside": function(pointInPolygon) {
assert.ok(!pointInPolygon([0, 0]));
assert.ok(!pointInPolygon([180, 1]));
assert.ok(!pointInPolygon([-90, -80]));
assert.ok(!pointInPolygon([180, 0]));
assert.ok(!pointInPolygon([150, 0]));
assert.ok(!pointInPolygon([180, -30]));
assert.ok(!pointInPolygon([150, -80]));
}
},
"triangle touching the North pole": {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment