Skip to content

Commit f86be42

Browse files
jasondaviesjasondavies
committed
Use point-in-polygon test for winding order.
This point-in-polygon algorithm uses the same area calculation as d3.geo.area, inlined for performance. Fixes d3#1300.
1 parent be2c8e7 commit f86be42

8 files changed

Lines changed: 524 additions & 111 deletions

File tree

‎d3.js‎

Lines changed: 49 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -2493,7 +2493,7 @@ d3 = function() {
24932493
a.next = b = array[0];
24942494
b.prev = a;
24952495
}
2496-
function d3_geo_clip(pointVisible, clipLine, interpolate) {
2496+
function d3_geo_clip(pointVisible, clipLine, interpolate, polygonContains) {
24972497
return function(listener) {
24982498
var line = clipLine(listener);
24992499
var clip = {
@@ -2504,9 +2504,8 @@ d3 = function() {
25042504
clip.point = pointRing;
25052505
clip.lineStart = ringStart;
25062506
clip.lineEnd = ringEnd;
2507-
invisible = false;
2508-
invisibleArea = visibleArea = 0;
25092507
segments = [];
2508+
polygon = [];
25102509
listener.polygonStart();
25112510
},
25122511
polygonEnd: function() {
@@ -2516,13 +2515,13 @@ d3 = function() {
25162515
segments = d3.merge(segments);
25172516
if (segments.length) {
25182517
d3_geo_clipPolygon(segments, d3_geo_clipSort, null, interpolate, listener);
2519-
} else if (visibleArea < -ε || invisible && invisibleArea < -ε) {
2518+
} else if (polygonContains(polygon)) {
25202519
listener.lineStart();
25212520
interpolate(null, null, 1, listener);
25222521
listener.lineEnd();
25232522
}
25242523
listener.polygonEnd();
2525-
segments = null;
2524+
segments = polygon = null;
25262525
},
25272526
sphere: function() {
25282527
listener.polygonStart();
@@ -2546,8 +2545,8 @@ d3 = function() {
25462545
clip.point = point;
25472546
line.lineEnd();
25482547
}
2549-
var segments, visibleArea, invisibleArea, invisible;
2550-
var buffer = d3_geo_clipBufferListener(), ringListener = clipLine(buffer), ring;
2548+
var segments;
2549+
var buffer = d3_geo_clipBufferListener(), ringListener = clipLine(buffer), polygon, ring;
25512550
function pointRing(λ, φ) {
25522551
ringListener.point(λ, φ);
25532552
ring.push([ λ, φ ]);
@@ -2560,16 +2559,12 @@ d3 = function() {
25602559
pointRing(ring[0][0], ring[0][1]);
25612560
ringListener.lineEnd();
25622561
var clean = ringListener.clean(), ringSegments = buffer.buffer(), segment, n = ringSegments.length;
2563-
if (!n) {
2564-
invisible = true;
2565-
invisibleArea += d3_geo_clipAreaRing(ring, -1);
2566-
ring = null;
2567-
return;
2568-
}
2562+
ring.pop();
2563+
polygon.push(ring);
25692564
ring = null;
2565+
if (!n) return;
25702566
if (clean & 1) {
25712567
segment = ringSegments[0];
2572-
visibleArea += d3_geo_clipAreaRing(segment, 1);
25732568
var n = segment.length - 1, i = -1, point;
25742569
listener.lineStart();
25752570
while (++i < n) listener.point((point = segment[i])[0], point[1]);
@@ -2606,26 +2601,40 @@ d3 = function() {
26062601
}
26072602
};
26082603
}
2609-
function d3_geo_clipAreaRing(ring, invisible) {
2610-
if (!(n = ring.length)) return 0;
2611-
var n, i = 0, area = 0, p = ring[0], λ = p[0], φ = p[1], cosφ = Math.cos(φ), x0 = Math.atan2(invisible * Math.sin(λ) * cosφ, Math.sin(φ)), y0 = 1 - invisible * Math.cos(λ) * cosφ, x1 = x0, x, y;
2612-
while (++i < n) {
2613-
p = ring[i];
2614-
cosφ = Math.cos(φ = p[1]);
2615-
x = Math.atan2(invisible * Math.sin(λ = p[0]) * cosφ, Math.sin(φ));
2616-
y = 1 - invisible * Math.cos(λ) * cosφ;
2617-
if (Math.abs(y0 - 2) < ε && Math.abs(y - 2) < ε) continue;
2618-
if (Math.abs(y) < ε || Math.abs(y0) < ε) {} else if (Math.abs(Math.abs(x - x0) - π) < ε) {
2619-
if (y + y0 > 2) area += 4 * (x - x0);
2620-
} else if (Math.abs(y0 - 2) < ε) area += 4 * (x - x1); else area += ((3 * π + x - x0) % (2 * π) - π) * (y0 + y);
2621-
x1 = x0, x0 = x, y0 = y;
2622-
}
2623-
return area;
2624-
}
26252604
function d3_geo_clipSort(a, b) {
26262605
return ((a = a.point)[0] < 0 ? a[1] - π / 2 - ε : π / 2 - a[1]) - ((b = b.point)[0] < 0 ? b[1] - π / 2 - ε : π / 2 - b[1]);
26272606
}
2628-
var d3_geo_clipAntimeridian = d3_geo_clip(d3_true, d3_geo_clipAntimeridianLine, d3_geo_clipAntimeridianInterpolate);
2607+
function d3_geo_pointInPolygon(point, polygon) {
2608+
var meridian = point[0], parallel = point[1], meridianNormal = [ Math.sin(meridian), -Math.cos(meridian), 0 ], polarAngle = 0, polar = false, southPole = false, winding = 0, area = 0;
2609+
for (var i = 0, n = polygon.length; i < n; ++i) {
2610+
var ring = polygon[i], m = ring.length;
2611+
if (!m) continue;
2612+
var point0 = ring[0], λ0 = point0[0], φ0 = point0[1] / 2 + π / 4, sinφ0 = Math.sin(φ0), cosφ0 = Math.cos(φ0), j = 1;
2613+
while (true) {
2614+
if (j === m) j = 0;
2615+
point = ring[j];
2616+
var λ = point[0], φ = point[1] / 2 + π / 4, sinφ = Math.sin(φ), cosφ = Math.cos(φ), dλ = λ - λ0, antimeridian = Math.abs(dλ) > π, k = sinφ0 * sinφ;
2617+
area += Math.atan2(k * Math.sin(dλ), cosφ0 * cosφ + k * Math.cos(dλ));
2618+
if (Math.abs(φ) < ε) southPole = true;
2619+
polarAngle += antimeridian ? dλ + (dλ >= 0 ? 2 : -2) * π : dλ;
2620+
if (antimeridian ^ λ0 >= meridian ^ λ >= meridian) {
2621+
var arc = d3_geo_cartesianCross(d3_geo_cartesian(point0), d3_geo_cartesian(point));
2622+
d3_geo_cartesianNormalize(arc);
2623+
var intersection = d3_geo_cartesianCross(meridianNormal, arc);
2624+
d3_geo_cartesianNormalize(intersection);
2625+
var φarc = (antimeridian ^ dλ >= 0 ? -1 : 1) * d3_asin(intersection[2]);
2626+
if (parallel > φarc) {
2627+
winding += antimeridian ^ dλ >= 0 ? 1 : -1;
2628+
}
2629+
}
2630+
if (!j++) break;
2631+
λ0 = λ, sinφ0 = sinφ, cosφ0 = cosφ, point0 = point;
2632+
}
2633+
if (Math.abs(polarAngle) > ε) polar = true;
2634+
}
2635+
return (!southPole && !polar && area < 0 || polarAngle < -ε) ^ winding & 1;
2636+
}
2637+
var d3_geo_clipAntimeridian = d3_geo_clip(d3_true, d3_geo_clipAntimeridianLine, d3_geo_clipAntimeridianInterpolate, d3_geo_clipAntimeridianPolygonContains);
26292638
function d3_geo_clipAntimeridianLine(listener) {
26302639
var λ0 = NaN, φ0 = NaN, sλ0 = NaN, clean;
26312640
return {
@@ -2692,9 +2701,13 @@ d3 = function() {
26922701
listener.point(to[0], to[1]);
26932702
}
26942703
}
2704+
var d3_geo_clipAntimeridianPoint = [ -π, 0 ];
2705+
function d3_geo_clipAntimeridianPolygonContains(polygon) {
2706+
return d3_geo_pointInPolygon(d3_geo_clipAntimeridianPoint, polygon);
2707+
}
26952708
function d3_geo_clipCircle(radius) {
2696-
var cr = Math.cos(radius), smallRadius = cr > 0, notHemisphere = Math.abs(cr) > ε, interpolate = d3_geo_circleInterpolate(radius, 6 * d3_radians);
2697-
return d3_geo_clip(visible, clipLine, interpolate);
2709+
var cr = Math.cos(radius), smallRadius = cr > 0, point = [ radius, 0 ], notHemisphere = Math.abs(cr) > ε, interpolate = d3_geo_circleInterpolate(radius, 6 * d3_radians);
2710+
return d3_geo_clip(visible, clipLine, interpolate, polygonContains);
26982711
function visible(λ, φ) {
26992712
return Math.cos(λ) * Math.cos(φ) > cr;
27002713
}
@@ -2787,6 +2800,9 @@ d3 = function() {
27872800
if (φ < -r) code |= 4; else if (φ > r) code |= 8;
27882801
return code;
27892802
}
2803+
function polygonContains(polygon) {
2804+
return d3_geo_pointInPolygon(point, polygon);
2805+
}
27902806
}
27912807
var d3_geo_clipViewMAX = 1e9;
27922808
function d3_geo_clipView(x0, y0, x1, y1) {

‎d3.min.js‎

Lines changed: 5 additions & 5 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

‎src/geo/clip-antimeridian.js‎

Lines changed: 12 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,13 @@
11
import "../core/true";
22
import "../math/trigonometry";
33
import "clip";
4+
import "point-in-polygon";
45

5-
var d3_geo_clipAntimeridian = d3_geo_clip(d3_true, d3_geo_clipAntimeridianLine, d3_geo_clipAntimeridianInterpolate);
6+
var d3_geo_clipAntimeridian = d3_geo_clip(
7+
d3_true,
8+
d3_geo_clipAntimeridianLine,
9+
d3_geo_clipAntimeridianInterpolate,
10+
d3_geo_clipAntimeridianPolygonContains);
611

712
// Takes a line and cuts into visible segments. Return values:
813
// 0: there were intersections or the line was empty.
@@ -88,3 +93,9 @@ function d3_geo_clipAntimeridianInterpolate(from, to, direction, listener) {
8893
listener.point(to[0], to[1]);
8994
}
9095
}
96+
97+
var d3_geo_clipAntimeridianPoint = [-π, 0];
98+
99+
function d3_geo_clipAntimeridianPolygonContains(polygon) {
100+
return d3_geo_pointInPolygon(d3_geo_clipAntimeridianPoint, polygon);
101+
}

‎src/geo/clip-circle.js‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,15 +3,17 @@ import "cartesian";
33
import "clip";
44
import "circle";
55
import "spherical";
6+
import "point-in-polygon";
67

78
// Clip features against a small circle centered at [0°, 0°].
89
function d3_geo_clipCircle(radius) {
910
var cr = Math.cos(radius),
1011
smallRadius = cr > 0,
12+
point = [radius, 0],
1113
notHemisphere = Math.abs(cr) > ε, // TODO optimise for this common case
1214
interpolate = d3_geo_circleInterpolate(radius, 6 * d3_radians);
1315

14-
return d3_geo_clip(visible, clipLine, interpolate);
16+
return d3_geo_clip(visible, clipLine, interpolate, polygonContains);
1517

1618
function visible(λ, φ) {
1719
return Math.cos(λ) * Math.cos(φ) > cr;
@@ -173,4 +175,8 @@ function d3_geo_clipCircle(radius) {
173175
else if (φ > r) code |= 8; // above
174176
return code;
175177
}
178+
179+
function polygonContains(polygon) {
180+
return d3_geo_pointInPolygon(point, polygon);
181+
}
176182
}

‎src/geo/clip.js‎

Lines changed: 10 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,7 @@ import "../core/noop";
33
import "../math/trigonometry";
44
import "clip-polygon";
55

6-
function d3_geo_clip(pointVisible, clipLine, interpolate) {
6+
function d3_geo_clip(pointVisible, clipLine, interpolate, polygonContains) {
77
return function(listener) {
88
var line = clipLine(listener);
99

@@ -15,9 +15,8 @@ function d3_geo_clip(pointVisible, clipLine, interpolate) {
1515
clip.point = pointRing;
1616
clip.lineStart = ringStart;
1717
clip.lineEnd = ringEnd;
18-
invisible = false;
19-
invisibleArea = visibleArea = 0;
2018
segments = [];
19+
polygon = [];
2120
listener.polygonStart();
2221
},
2322
polygonEnd: function() {
@@ -28,13 +27,13 @@ function d3_geo_clip(pointVisible, clipLine, interpolate) {
2827
segments = d3.merge(segments);
2928
if (segments.length) {
3029
d3_geo_clipPolygon(segments, d3_geo_clipSort, null, interpolate, listener);
31-
} else if (visibleArea < -ε || invisible && invisibleArea < -ε) {
30+
} else if (polygonContains(polygon)) {
3231
listener.lineStart();
3332
interpolate(null, null, 1, listener);
3433
listener.lineEnd();
3534
}
3635
listener.polygonEnd();
37-
segments = null;
36+
segments = polygon = null;
3837
},
3938
sphere: function() {
4039
listener.polygonStart();
@@ -50,13 +49,11 @@ function d3_geo_clip(pointVisible, clipLine, interpolate) {
5049
function lineStart() { clip.point = pointLine; line.lineStart(); }
5150
function lineEnd() { clip.point = point; line.lineEnd(); }
5251

53-
var segments,
54-
visibleArea,
55-
invisibleArea,
56-
invisible;
52+
var segments;
5753

5854
var buffer = d3_geo_clipBufferListener(),
5955
ringListener = clipLine(buffer),
56+
polygon,
6057
ring;
6158

6259
function pointRing(λ, φ) {
@@ -78,20 +75,15 @@ function d3_geo_clip(pointVisible, clipLine, interpolate) {
7875
segment,
7976
n = ringSegments.length;
8077

81-
// TODO compute on-the-fly?
82-
if (!n) {
83-
invisible = true;
84-
invisibleArea += d3_geo_clipAreaRing(ring, -1);
85-
ring = null;
86-
return;
87-
}
78+
ring.pop();
79+
polygon.push(ring);
8880
ring = null;
8981

82+
if (!n) return;
83+
9084
// No intersections.
91-
// TODO compute on-the-fly?
9285
if (clean & 1) {
9386
segment = ringSegments[0];
94-
visibleArea += d3_geo_clipAreaRing(segment, 1);
9587
var n = segment.length - 1,
9688
i = -1,
9789
point;
@@ -135,59 +127,6 @@ function d3_geo_clipBufferListener() {
135127
};
136128
}
137129

138-
// Approximate polygon ring area (×2, since we only need the sign).
139-
// For an invisible polygon ring, we rotate longitudinally by 180°.
140-
// The invisible parameter should be 1, or -1 to rotate longitudinally.
141-
// Based on Robert. G. Chamberlain and William H. Duquette,
142-
// “Some Algorithms for Polygons on a Sphere”,
143-
// http://trs-new.jpl.nasa.gov/dspace/handle/2014/40409
144-
function d3_geo_clipAreaRing(ring, invisible) {
145-
if (!(n = ring.length)) return 0;
146-
var n,
147-
i = 0,
148-
area = 0,
149-
p = ring[0],
150-
λ = p[0],
151-
φ = p[1],
152-
cosφ = Math.cos(φ),
153-
x0 = Math.atan2(invisible * Math.sin(λ) * cosφ, Math.sin(φ)),
154-
y0 = 1 - invisible * Math.cos(λ) * cosφ,
155-
x1 = x0,
156-
x, // λ'; λ rotated to south pole.
157-
y; // φ' = 1 + sin(φ); φ rotated to south pole.
158-
while (++i < n) {
159-
p = ring[i];
160-
cosφ = Math.cos(φ = p[1]);
161-
x = Math.atan2(invisible * Math.sin(λ = p[0]) * cosφ, Math.sin(φ));
162-
y = 1 - invisible * Math.cos(λ) * cosφ;
163-
164-
// If both the current point and the previous point are at the north pole,
165-
// skip this point.
166-
if (Math.abs(y0 - 2) < ε && Math.abs(y - 2) < ε) continue;
167-
168-
// If this or the previous point is at the south pole, or if this segment
169-
// goes through the south pole, the area is 0.
170-
if (Math.abs(y) < ε || Math.abs(y0) < ε) {}
171-
172-
// If this segment goes through either pole…
173-
else if (Math.abs(Math.abs(x - x0) - π) < ε) {
174-
// For the north pole, compute lune area.
175-
if (y + y0 > 2) area += 4 * (x - x0);
176-
// For the south pole, the area is zero.
177-
}
178-
179-
// If the previous point is at the north pole, then compute lune area.
180-
else if (Math.abs(y0 - 2) < ε) area += 4 * (x - x1);
181-
182-
// Otherwise, the spherical triangle area is approximately
183-
// δλ * (1 + sinφ0 + 1 + sinφ) / 2.
184-
else area += ((3 * π + x - x0) % (2 * π) - π) * (y0 + y);
185-
186-
x1 = x0, x0 = x, y0 = y;
187-
}
188-
return area;
189-
}
190-
191130
// Intersection points are sorted along the clip edge. For both antimeridian
192131
// cutting and circle clipping, the same comparison is used.
193132
function d3_geo_clipSort(a, b) {

0 commit comments

Comments
 (0)