geoIntersectArc + lots of test cases and bug fixes in the Spherical i…

…ntersection function.

Solves d3/d3-geo#160

README.md

@@ -39,6 +39,15 @@ Given a GeoJSON *polygon* or *multipolygon*, returns a clip function suitable fo
 Given a clipPolygon function, returns the GeoJSON polygon.
 
 
+<a name="geoIntersectArc" href="#geoIntersectArc">#</a> d3.<b>geoIntersectArc</b>(<i>arcs</i>) [<>](https://github.com/d3/d3-geo-polygon/blob/master/src/intersect.js "Source")
+
+Given two spherical arcs [point0, point1] and [point2, point3], returns their intersection, or undefined if there is none. See “[Spherical Intersection](https://observablehq.com/@fil/spherical-intersection)”.
+
+<a name="polygon" href="#polygon">#</a> clip.<b>polygon</b>()
+
+Given a clipPolygon function, returns the GeoJSON polygon.
+
+
 
 ## Projections
 

src/index.js

@@ -1,4 +1,5 @@
 export {default as geoClipPolygon} from "./clip/polygon";
+export {default as geoIntersectArc} from "./intersect";
 export {default as geoPolyhedral} from "./polyhedral/index";
 export {default as geoPolyhedralButterfly} from "./polyhedral/butterfly";
 export {default as geoPolyhedralCollignon} from "./polyhedral/collignon";

src/intersect.js

@@ -1,16 +1,21 @@
-import {abs, acos, cos, epsilon, epsilon2} from "./math";
-import {cartesianCross, cartesianDot, cartesianEqual, cartesianNormalizeInPlace} from "./cartesian";
+import {abs, acos, cos, degrees, epsilon, epsilon2, max, radians} from "./math";
+import {cartesian, cartesianCross, cartesianDot, cartesianEqual, cartesianNormalizeInPlace, spherical} 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);
-  this.l = acos(cartesianDot(from, to));
+  this.l = acos(max(0,cartesianDot(from, to)));
 }
 
 // >> here a and b are segments processed by intersectSegment
 export function intersect(a, b) {
+  if (cartesianEqual(a.from, b.from) || cartesianEqual(a.from, b.to))
+    return a.from;
+  if (cartesianEqual(a.to, b.from) || cartesianEqual(a.to, b.to))
+    return a.to;
+
   var lc = cos(a.l + b.l) - epsilon;
   if (cartesianDot(a.from, b.from) < lc
   || cartesianDot(a.from, b.to) < lc
@@ -30,10 +35,14 @@ export function intersect(a, b) {
   // or is almost equal to one of the four points
   if (
     (a0 > 0 && a1 < 0 && b0 > 0 && b1 < 0) ||
-    cartesianEqual(axb, a.from) ||
-    cartesianEqual(axb, a.to) ||
-    cartesianEqual(axb, b.from) ||
-    cartesianEqual(axb, b.to)
+    (a0 >= 0 &&
+      a1 <= 0 &&
+      b0 >= 0 &&
+      b1 <= 0 &&
+      (cartesianEqual(axb, a.from) ||
+        cartesianEqual(axb, a.to) ||
+        cartesianEqual(axb, b.from) ||
+        cartesianEqual(axb, b.to)))
   )
     return axb;
 
@@ -48,10 +57,14 @@ export function intersect(a, b) {
 
   if (
     (a0 > 0 && a1 < 0 && b0 > 0 && b1 < 0) ||
-    cartesianEqual(axb, a.from) ||
-    cartesianEqual(axb, a.to) ||
-    cartesianEqual(axb, b.from) ||
-    cartesianEqual(axb, b.to)
+    (a0 >= 0 &&
+      a1 <= 0 &&
+      b0 >= 0 &&
+      b1 <= 0 &&
+      (cartesianEqual(axb, a.from) ||
+        cartesianEqual(axb, a.to) ||
+        cartesianEqual(axb, b.from) ||
+        cartesianEqual(axb, b.to)))
   )
     return axb;
 }
@@ -66,6 +79,12 @@ export function intersectPointOnLine(p, a) {
 
 export var intersectCoincident = {};
 
-
-// todo: publicly expose d3.geoIntersect(segment0, segment1) ??
-// cf. https://github.com/d3/d3/commit/3dbdf87974dc2588c29db0533a8500ccddb25daa#diff-65daf69cea7d039d72c1eca7c13326b0
+export default function(a, b) {
+  var ca = a.map(p => cartesian(p.map(d => d * radians))),
+      cb = b.map(p => cartesian(p.map(d => d * radians)));
+  var i = intersect(
+    new intersectSegment(ca[0], ca[1]),
+    new intersectSegment(cb[0], cb[1])
+  );
+  return !i ? i : spherical(i).map(d => d * degrees);
+}

test/inDelta.js

@@ -0,0 +1,26 @@
+var tape = require("tape");
+
+tape.Test.prototype.inDelta = function(actual, expected, delta) {
+  delta = delta || 1e-6;
+  this._assert(inDelta(actual, expected, delta), {
+    message: "should be in delta " + delta,
+    operator: "inDelta",
+    actual: actual,
+    expected: expected
+  });
+};
+
+function inDelta(actual, expected, delta) {
+  return (Array.isArray(expected) ? inDeltaArray : inDeltaNumber)(actual, expected, delta);
+}
+
+function inDeltaArray(actual, expected, delta) {
+  var n = expected.length, i = -1;
+  if (actual.length !== n) return false;
+  while (++i < n) if (!inDelta(actual[i], expected[i], delta)) return false;
+  return true;
+}
+
+function inDeltaNumber(actual, expected, delta) {
+  return actual >= expected - delta && actual <= expected + delta;
+}

test/intersect-test.js

@@ -0,0 +1,46 @@
+var tape = require("tape"),
+    d3 = require("../"),
+    d3_geo = require("d3-geo");
+
+require("./inDelta");
+
+var e;
+
+tape("spherical intersections", function(test) {
+  e = d3.geoIntersectArc([[0,0], [0,90]], [[-10,40], [10,40]]); 
+  test.assert(e);
+  test.inDelta(e, [ 0, 40.43246108], 1e-8);
+
+  // https://observablehq.com/@fil/spherical-intersection#points
+  var p = [[0, 70], [-10, 10], [-40, 30], [10, 45]];
+
+  test.inDelta(d3.geoIntersectArc([p[0], p[1]], [p[0], p[3]]), p[0], 1e-8);
+  test.inDelta(d3.geoIntersectArc([p[0], p[1]], [p[1], p[3]]), p[1], 1e-8);
+  test.inDelta(d3.geoIntersectArc([p[0], p[1]], [p[2], p[0]]), p[0], 1e-8);
+  test.inDelta(d3.geoIntersectArc([p[0], p[1]], [p[2], p[1]]), p[1], 1e-8);
+
+  test.inDelta(
+    d3.geoIntersectArc([p[0], p[1]], [p[2], p[3]]),
+    [ -7.081398732358556, 42.94731141237317 ], 1e-8);
+
+  test.inDelta(
+    d3.geoIntersectArc([p[1], p[0]], [p[2], p[3]]),
+    [ -7.081398732358556, 42.94731141237317 ], 1e-8);
+
+  test.false(d3.geoIntersectArc([p[2], p[1]], [p[0], p[3]]));
+
+  test.inDelta(
+    d3.geoIntersectArc([[0, 89.99], [0, -89.99]], [[-89.99, 0], [89.99, 0]]),
+    [0,0], 1e-8);
+
+  test.inDelta(
+    d3.geoIntersectArc([[0, 89.99], [0, -89.99]], [[0, 0], [25, 0]]),
+    [0,0], 1e-8);
+
+  e = d3.geoIntersectArc([[0,0], [0,90]], [[0,0], [90,0]]); 
+  test.deepEqual(e, [0, 0]);
+
+  test.false(d3.geoIntersectArc([[0,0], [0,90]], [[10,0], [90,0]]));
+
+  test.end();
+});