Jason Davies’ clip-polygon, ported to d3v4 by Philippe Rivière (see d…

…3/d3-geo-projection#124)

src/cartesian.js

@@ -1,4 +1,4 @@
-import {asin, atan2, cos, sin, sqrt} from "./math";
+import {asin, atan2, cos, epsilon2, sin, sqrt} from "./math";
 
 export function spherical(cartesian) {
   return [atan2(cartesian[1], cartesian[0]), asin(cartesian[2])];
@@ -31,3 +31,10 @@ export function cartesianNormalizeInPlace(d) {
   var l = sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]);
   d[0] /= l, d[1] /= l, d[2] /= l;
 }
+
+export function cartesianEqual(a, b) {
+  var dx = b[0] - a[0],
+      dy = b[1] - a[1],
+      dz = b[2] - a[2];
+  return dx * dx + dy * dy + dz * dz < epsilon2 * epsilon2;
+}

src/clip/buffer.js

@@ -4,8 +4,11 @@ export default function() {
   var lines = [],
       line;
   return {
-    point: function(x, y) {
-      line.push([x, y]);
+    point: function(x, y, i, t) {
+      var point = [x, y];
+      // when called by clipPolygon, store index and t
+      if (arguments.length > 2) { point.index = i; point.t = t; }
+      line.push(point);
     },
     lineStart: function() {
       lines.push(line = []);

src/clip/index.js

@@ -4,7 +4,9 @@ import {epsilon, halfPi} from "../math";
 import polygonContains from "../polygonContains";
 import {merge} from "d3-array";
 
-export default function(pointVisible, clipLine, interpolate, start) {
+export default function(pointVisible, clipLine, interpolate, start, sort) {
+  if (typeof sort === "undefined") sort = compareIntersection;
+
   return function(sink) {
     var line = clipLine(sink),
         ringBuffer = clipBuffer(),
@@ -33,7 +35,7 @@ export default function(pointVisible, clipLine, interpolate, start) {
         var startInside = polygonContains(polygon, start);
         if (segments.length) {
           if (!polygonStarted) sink.polygonStart(), polygonStarted = true;
-          clipRejoin(segments, compareIntersection, startInside, interpolate, sink);
+          clipRejoin(segments, sort, startInside, interpolate, sink);
         } else if (startInside) {
           if (!polygonStarted) sink.polygonStart(), polygonStarted = true;
           sink.lineStart();
@@ -70,9 +72,9 @@ export default function(pointVisible, clipLine, interpolate, start) {
       line.lineEnd();
     }
 
-    function pointRing(lambda, phi) {
+    function pointRing(lambda, phi, close) {
       ring.push([lambda, phi]);
-      ringSink.point(lambda, phi);
+      ringSink.point(lambda, phi, close);
     }
 
     function ringStart() {
@@ -81,7 +83,7 @@ export default function(pointVisible, clipLine, interpolate, start) {
     }
 
     function ringEnd() {
-      pointRing(ring[0][0], ring[0][1]);
+      pointRing(ring[0][0], ring[0][1], true);
       ringSink.lineEnd();
 
       var clean = ringSink.clean(),

src/clip/polygon.js

@@ -1,3 +1,148 @@
-export default function() {
-  throw new Error("not yet implemented");
+import clip from "./index";
+import {atan2, cos, pi, radians, sign, sin, sqrt} from "../math";
+import {cartesian, cartesianCross, cartesianDot, cartesianEqual, spherical} from "../cartesian";
+import {intersectCoincident, intersectPointOnLine, intersectSegment, intersect} from "../intersect";
+import {default as polygonContains} from "../polygonContains";
+
+var clipNone = function(stream) { return stream; };
+
+// clipPolygon
+export default function (p) {
+  var segments = [];
+
+  if (p.type != "Polygon") return clipNone; // todo: MultiPolygon?
+
+  var polygon = p.coordinates.map(function(ring) {
+    var c, c0;
+    ring = ring.map(function(point, i) {
+      c = cartesian(point = [point[0] * radians, point[1] * radians]);
+      if (i) segments.push(new intersectSegment(c0, c));
+      c0 = c;
+      return point;
+    });
+    ring.pop();
+    return ring;
+  });
+
+
+  function visible(lambda, phi) {
+    return polygonContains(polygon, [lambda, phi]);
+  }
+
+  function clipLine(stream) {
+    var point0,
+        lambda00,
+        phi00,
+        v00,
+        v0,
+        clean;
+    return {
+      lineStart: function() {
+        point0 = null;
+        clean = 1;
+      },
+      point: function(lambda, phi, close) {
+        if (cos(lambda) == -1) lambda -= sign(sin(lambda)) * 1e-5; // move away from -180/180 https://github.com/d3/d3-geo/pull/108#issuecomment-323798937
+        if (close) lambda = lambda00, phi = phi00;
+        var point = cartesian([lambda, phi]),
+            v = v0,
+            intersection,
+            i, j, s, t;
+        if (point0) {
+          var segment = new intersectSegment(point0, point),
+              intersections = [];
+          for (i = 0, j = 100; i < segments.length && j > 0; ++i) {
+            s = segments[i];
+            intersection = intersect(segment, s);
+            if (intersection) {
+              if (intersection === intersectCoincident ||
+                  cartesianEqual(intersection, point0) || cartesianEqual(intersection, point) ||
+                  cartesianEqual(intersection, s.from) || cartesianEqual(intersection, s.to)) {
+                t = 1e-4;
+                lambda = (lambda + 3 * pi + (Math.random() < .5 ? t : -t)) % (2 * pi) - pi;
+                phi = Math.min(pi / 2 - 1e-4, Math.max(1e-4 - pi / 2, phi + (Math.random() < .5 ? t : -t)));
+                segment = new intersectSegment(point0, point = cartesian([lambda, phi]));
+                i = -1, --j;
+                intersections.length = 0;
+                continue;
+              }
+              var sph = spherical(intersection);
+              intersection.distance = clipPolygonDistance(point0, intersection);
+              intersection.index = i;
+              intersection.t = clipPolygonDistance(s.from, intersection);
+              intersection[0] = sph[0], intersection[1] = sph[1], intersection.pop();
+              intersections.push(intersection);
+            }
+          }
+          if (intersections.length) {
+            clean = 0;
+            intersections.sort(function(a, b) { return a.distance - b.distance; });
+            for (i = 0; i < intersections.length; ++i) {
+              intersection = intersections[i];
+              v = !v;
+              if (v) {
+                stream.lineStart();
+                stream.point(intersection[0], intersection[1], intersection.index, intersection.t);
+              } else {
+                stream.point(intersection[0], intersection[1], intersection.index, intersection.t);
+                stream.lineEnd();
+              }
+            }
+          }
+          if (v) stream.point(lambda, phi);
+        } else {
+          for (i = 0, j = 100; i < segments.length && j > 0; ++i) {
+            s = segments[i];
+            if (intersectPointOnLine(point, s)) {
+              t = 1e-4;
+              lambda = (lambda + 3 * pi + (Math.random() < .5 ? t : -t)) % (2 * pi) - pi;
+              phi = Math.min(pi / 2 - 1e-4, Math.max(1e-4 - pi / 2, phi + (Math.random() < .5 ? t : -t)));
+              point = cartesian([lambda, phi]);
+              i = -1, --j;
+            }
+          }
+          v00 = v = visible(lambda00 = lambda, phi00 = phi);
+          if (v) stream.lineStart(), stream.point(lambda, phi);
+        }
+        point0 = point, v0 = v;
+      },
+      lineEnd: function() {
+        if (v0) stream.lineEnd();
+      },
+      // Rejoin first and last segments if there were intersections and the first
+      // and last points were visible.
+      clean: function() {
+        return clean | ((v00 && v0) << 1);
+      }
+    };
+  }
+
+  function interpolate(from, to, direction, stream) {
+    if (from == null) {
+      var n = polygon.length;
+      polygon.forEach(function(ring, i) {
+        ring.forEach(function(point) { stream.point(point[0], point[1]); });
+        if (i < n - 1) stream.lineEnd(), stream.lineStart();
+      });
+    } else if (from.index !== to.index && from.index != null && to.index != null) {
+      for (var i = from.index; i !== to.index; i = (i + direction + segments.length) % segments.length) {
+        var segment = segments[i],
+            point = spherical(direction > 0 ? segment.to : segment.from);
+        stream.point(point[0], point[1]);
+      }
+    }
+  }
+
+  return clip(visible, clipLine, interpolate, polygon[0][0], clipPolygonSort);
 }
+
+function clipPolygonSort(a, b) {
+  a = a.x, b = b.x;
+  return a.index - b.index || a.t - b.t;
+}
+
+// Geodesic coordinates for two 3D points.
+function clipPolygonDistance(a, b) {
+  var axb = cartesianCross(a, b);
+  return atan2(sqrt(cartesianDot(axb, axb)), cartesianDot(a, b));
+}

src/intersect.js

@@ -0,0 +1,41 @@
+import {abs, epsilon2} from "./math";
+import {cartesianCross, cartesianDot, cartesianNormalizeInPlace} from "./cartesian";
+
+export function intersectSegment(from, to) {
+  this.from = from, this.to = to;
+  this.normal = cartesianCross(from, to);
+  this.fromNormal = cartesianCross(this.normal, from);
+  this.toNormal = cartesianCross(this.normal, to);
+}
+
+
+// >> here a and b are segments processed by intersectSegment
+export function intersect(a, b) {
+  var axb = cartesianCross(a.normal, b.normal);
+  cartesianNormalizeInPlace(axb);
+  var a0 = cartesianDot(axb, a.fromNormal),
+      a1 = cartesianDot(axb, a.toNormal),
+      b0 = cartesianDot(axb, b.fromNormal),
+      b1 = cartesianDot(axb, b.toNormal);
+
+  if (a0 > -epsilon2 && a1 < epsilon2 && b0 > -epsilon2 && b1 < epsilon2) return axb;
+
+  if (a0 < epsilon2 && a1 > -epsilon2 && b0 < epsilon2 && b1 > -epsilon2) {
+    axb[0] = -axb[0], axb[1] = -axb[1], axb[2] = -axb[2];
+    return axb;
+  }
+}
+
+export function intersectPointOnLine(p, a) {
+  var a0 = cartesianDot(p, a.fromNormal),
+      a1 = cartesianDot(p, a.toNormal);
+  p = cartesianDot(p, a.normal);
+
+  return abs(p) < epsilon2 && (a0 > -epsilon2 && a1 < epsilon2 || a0 < epsilon2 && a1 > -epsilon2);
+}
+
+export var intersectCoincident = {};
+
+
+// todo: publicly expose d3.geoIntersect(segment0, segment1) ??
+// cf. https://github.com/d3/d3/commit/3dbdf87974dc2588c29db0533a8500ccddb25daa#diff-65daf69cea7d039d72c1eca7c13326b0