fix(s2): EdgeCrosser

s2/EdgeCrosser.ts

@@ -1,4 +1,4 @@
-import { CROSS, Crossing, DO_NOT_CROSS, MAYBE_CROSS } from './edge_crossings'
+import { CROSS, Crossing, DO_NOT_CROSS, MAYBE_CROSS, vertexCrossing } from './edge_crossings'
 import { DBL_EPSILON, expensiveSign, INDETERMINATE, robustSign, triageSign } from './predicates'
 import type { Direction } from './predicates'
 import { Point } from './Point'
@@ -22,7 +22,6 @@ import { Point } from './Point'
  *   return count;
  * }
  * ```
- * @beta incomplete
  */
 export class EdgeCrosser {
   a: Point
@@ -82,10 +81,25 @@ export class EdgeCrosser {
    * chainCrossingSign below.
    */
   crossingSign(c: Point, d: Point): Crossing {
-    if (c !== this.c) this.restartAt(c)
+    if (!c.equals(this.c)) this.restartAt(c)
     return this.chainCrossingSign(d)
   }
 
+  /**
+   * Reports whether if CrossingSign(c, d) > 0, or AB and CD share a vertex and VertexCrossing(a, b, c, d) is true.
+   *
+   * This method extends the concept of a "crossing" to the case where AB
+   * and CD have a vertex in common. The two edges may or may not cross,
+   * according to the rules defined in VertexCrossing above. The rules
+   * are designed so that point containment tests can be implemented simply
+   * by counting edge crossings. Similarly, determining whether one edge
+   * chain crosses another edge chain can be implemented by counting.
+   */
+  edgeOrVertexCrossing(c: Point, d: Point): boolean {
+    if (!c.equals(this.c)) this.restartAt(c)
+    return this.edgeOrVertexChainCrossing(d)
+  }
+
   /**
    * Sets the current point of the edge crosser to be c.
    * Call this method when your chain 'jumps' to a new place.
@@ -101,20 +115,63 @@ export class EdgeCrosser {
    * the crossing methods (or restartAt) as the first vertex of the current edge.
    */
   chainCrossingSign(d: Point): Crossing {
+    // For there to be an edge crossing, the triangles ACB, CBD, BDA, DAC must
+    // all be oriented the same way (CW or CCW). We keep the orientation of ACB
+    // as part of our state. When each new point D arrives, we compute the
+    // orientation of BDA and check whether it matches ACB. This checks whether
+    // the points C and D are on opposite sides of the great circle through AB.
+
+    // Recall that triageSign is invariant with respect to rotating its
+    // arguments, i.e. ABD has the same orientation as BDA.
     const bda = triageSign(this.a, this.b, d)
     if (this.acb === -bda && bda !== INDETERMINATE) {
       this.c = d
       this.acb = -bda
       return DO_NOT_CROSS
     }
-    return this.crossingSignHelper(d, bda)
+
+    return this._crossingSign(d, bda)
+  }
+
+  /**
+   * Like EdgeOrVertexCrossing, but uses the last vertex
+   * passed to one of the crossing methods (or RestartAt) as the first vertex of the current edge.
+   */
+  edgeOrVertexChainCrossing(d: Point): boolean {
+    // We need to copy this.c since it is clobbered by ChainCrossingSign.
+    const c = Point.fromVector(this.c.vector)
+
+    switch (this.chainCrossingSign(d)) {
+      case DO_NOT_CROSS:
+        return false
+      case CROSS:
+        return true
+    }
+
+    return vertexCrossing(this.a, this.b, c, d)
   }
 
   /**
    * Handle the slow path of crossingSign.
    */
-  private crossingSignHelper(d: Point, bda: Direction): Crossing {
+  private _crossingSign(d: Point, bda: Direction): Crossing {
     const maxError = (1.5 + 1 / Math.sqrt(3)) * DBL_EPSILON
+
+    // At this point, a very common situation is that A,B,C,D are four points on
+    // a line such that AB does not overlap CD. (For example, this happens when
+    // a line or curve is sampled finely, or when geometry is constructed by
+    // computing the union of S2CellIds.) Most of the time, we can determine
+    // that AB and CD do not intersect using the two outward-facing
+    // tangents at A and B (parallel to AB) and testing whether AB and CD are on
+    // opposite sides of the plane perpendicular to one of these tangents. This
+    // is moderately expensive but still much cheaper than expensiveSign.
+
+    // The error in RobustCrossProd is insignificant. The maximum error in
+    // the call to CrossProd (i.e., the maximum norm of the error vector) is
+    // (0.5 + 1/sqrt(3)) * dblEpsilon. The maximum error in each call to
+    // DotProd below is dblEpsilon. (There is also a small relative error
+    // term that is insignificant because we are comparing the result against a
+    // constant that is very close to zero.)
     if (
       (this.c.vector.dot(this.aTangent.vector) > maxError && d.vector.dot(this.aTangent.vector) > maxError) ||
       (this.c.vector.dot(this.bTangent.vector) > maxError && d.vector.dot(this.bTangent.vector) > maxError)
@@ -124,18 +181,25 @@ export class EdgeCrosser {
       return DO_NOT_CROSS
     }
 
-    if (this.a === this.c || this.a === d || this.b === this.c || this.b === d) {
+    // Otherwise, eliminate the cases where two vertices from different edges are
+    // equal. (These cases could be handled in the code below, but we would rather
+    // avoid calling ExpensiveSign if possible.)
+    if (this.a.equals(this.c) || this.a.equals(d) || this.b.equals(this.c) || this.b.equals(d)) {
       this.c = d
       this.acb = -bda
       return MAYBE_CROSS
     }
 
-    if (this.a === this.b || this.c === d) {
+    // Eliminate the cases where an input edge is degenerate. (Note that in
+    // most cases, if CD is degenerate then this method is not even called
+    // because acb and bda have different signs.)
+    if (this.a.equals(this.b) || this.c.equals(d)) {
       this.c = d
       this.acb = -bda
       return DO_NOT_CROSS
     }
 
+    // Otherwise it's time to break out the big guns.
     if (this.acb === INDETERMINATE) this.acb = -expensiveSign(this.a, this.b, this.c)
     if (bda === INDETERMINATE) bda = expensiveSign(this.a, this.b, d)
 

s2/EdgeCrosser_test.ts

@@ -0,0 +1,171 @@
+import { test } from 'node:test'
+import { strict as assert } from 'node:assert'
+import { nextAfter } from '../r1/math'
+import { Point } from './Point'
+import { CROSS, Crossing, DO_NOT_CROSS, MAYBE_CROSS } from './edge_crossings'
+import { EdgeCrosser } from './EdgeCrosser'
+
+test('crossings', () => {
+  const NA1 = nextAfter(1, 0)
+  const NA2 = nextAfter(1, 2)
+
+  const tests = [
+    {
+      msg: 'two regular edges that cross',
+      a: new Point(1, 2, 1),
+      b: new Point(1, -3, 0.5),
+      c: new Point(1, -0.5, -3),
+      d: new Point(0.1, 0.5, 3),
+      robust: CROSS,
+      edgeOrVertex: true
+    },
+    {
+      msg: 'two regular edges that intersect antipodal points',
+      a: new Point(1, 2, 1),
+      b: new Point(1, -3, 0.5),
+      c: new Point(-1, 0.5, 3),
+      d: new Point(-0.1, -0.5, -3),
+      robust: DO_NOT_CROSS,
+      edgeOrVertex: false
+    },
+    {
+      msg: 'two edges on the same great circle that start at antipodal points',
+      a: new Point(0, 0, -1),
+      b: new Point(0, 1, 0),
+      c: new Point(0, 0, 1),
+      d: new Point(0, 1, 1),
+      robust: DO_NOT_CROSS,
+      edgeOrVertex: false
+    },
+    {
+      msg: 'two edges that cross where one vertex is the OriginPoint',
+      a: new Point(1, 0, 0),
+      b: Point.originPoint(),
+      c: new Point(1, -0.1, 1),
+      d: new Point(1, 1, -0.1),
+      robust: CROSS,
+      edgeOrVertex: true
+    },
+    {
+      msg: 'two edges that intersect antipodal points where one vertex is the OriginPoint',
+      a: new Point(1, 0, 0),
+      b: Point.originPoint(),
+      c: new Point(1, 0.1, -1),
+      d: new Point(1, 1, -0.1),
+      robust: DO_NOT_CROSS,
+      edgeOrVertex: false
+    },
+    {
+      msg: 'two edges that cross antipodal points',
+      a: new Point(1, 0, 0),
+      b: new Point(0, 1, 0),
+      c: new Point(0, 0, -1),
+      d: new Point(-1, -1, 1),
+      robust: DO_NOT_CROSS,
+      edgeOrVertex: false
+    },
+    {
+      msg: 'two edges that share an endpoint',
+      a: new Point(2, 3, 4),
+      b: new Point(-1, 2, 5),
+      c: new Point(7, -2, 3),
+      d: new Point(2, 3, 4),
+      robust: MAYBE_CROSS,
+      edgeOrVertex: false
+    },
+    {
+      msg: 'two edges that barely cross near the middle of one edge',
+      a: new Point(1, 1, 1),
+      b: new Point(1, NA1, -1),
+      c: new Point(11, -12, -1),
+      d: new Point(10, 10, 1),
+      robust: CROSS,
+      edgeOrVertex: true
+    },
+    {
+      msg: 'two edges that barely cross near the middle separated by a distance of about 1e-15',
+      a: new Point(1, 1, 1),
+      b: new Point(1, NA2, -1),
+      c: new Point(1, -1, 0),
+      d: new Point(1, 1, 0),
+      robust: DO_NOT_CROSS,
+      edgeOrVertex: false
+    },
+    {
+      msg: 'two edges that barely cross each other near the end of both edges',
+      a: new Point(0, 0, 1),
+      b: new Point(2, -1e-323, 1),
+      c: new Point(1, -1, 1),
+      d: new Point(1e-323, 0, 1),
+      robust: CROSS,
+      edgeOrVertex: true
+    },
+    {
+      msg: 'two edges that barely cross each other near the end separated by a distance of about 1e-640',
+      a: new Point(0, 0, 1),
+      b: new Point(2, 1e-323, 1),
+      c: new Point(1, -1, 1),
+      d: new Point(1e-323, 0, 1),
+      robust: DO_NOT_CROSS,
+      edgeOrVertex: false
+    },
+    {
+      msg: 'two edges that barely cross each other near the middle of one edge',
+      a: new Point(1, -1e-323, -1e-323),
+      b: new Point(1e-323, 1, 1e-323),
+      c: new Point(1, -1, 1e-323),
+      d: new Point(1, 1, 0),
+      robust: CROSS,
+      edgeOrVertex: true
+    },
+    {
+      msg: 'two edges that barely cross each other near the middle separated by a distance of about 1e-640',
+      a: new Point(1, 1e-323, -1e-323),
+      b: new Point(-1e-323, 1, 1e-323),
+      c: new Point(1, -1, 1e-323),
+      d: new Point(1, 1, 0),
+      robust: DO_NOT_CROSS,
+      edgeOrVertex: false
+    }
+  ]
+
+  tests.forEach((test) => {
+    const a = Point.fromVector(test.a.vector.normalize())
+    const b = Point.fromVector(test.b.vector.normalize())
+    const c = Point.fromVector(test.c.vector.normalize())
+    const d = Point.fromVector(test.d.vector.normalize())
+
+    testCrossing(test.msg, a, b, c, d, test.robust, test.edgeOrVertex)
+    testCrossing(test.msg, b, a, c, d, test.robust, test.edgeOrVertex)
+    testCrossing(test.msg, a, b, d, c, test.robust, test.edgeOrVertex)
+    testCrossing(test.msg, b, a, d, c, test.robust, test.edgeOrVertex)
+
+    // test degenerate cases
+    testCrossing(test.msg, a, a, c, d, DO_NOT_CROSS, false)
+    testCrossing(test.msg, a, b, c, c, DO_NOT_CROSS, false)
+    testCrossing(test.msg, a, a, c, c, DO_NOT_CROSS, false)
+
+    testCrossing(test.msg, a, b, a, b, MAYBE_CROSS, true)
+    testCrossing(test.msg, c, d, a, b, test.robust, test.edgeOrVertex !== (test.robust === MAYBE_CROSS))
+  })
+})
+
+function testCrossing(msg: string, a: Point, b: Point, c: Point, d: Point, robust: Crossing, edgeOrVertex: boolean) {
+  if (a.equals(c) || a.equals(d) || b.equals(c) || b.equals(d)) {
+    robust = MAYBE_CROSS
+  }
+
+  const input = `${msg}: a: ${a}, b: ${b}, c: ${c}, d: ${d}`
+
+  const crosser = EdgeCrosser.newChainEdgeCrosser(a, b, c)
+  assert.equal(crosser.chainCrossingSign(d), robust, `${input}, ChainCrossingSign(d)`)
+  assert.equal(crosser.chainCrossingSign(c), robust, `${input}, ChainCrossingSign(c)`)
+  assert.equal(crosser.crossingSign(d, c), robust, `${input}, CrossingSign(d, c)`)
+  assert.equal(crosser.crossingSign(c, d), robust, `${input}, CrossingSign(c, d)`)
+
+  crosser.restartAt(c)
+  assert.equal(crosser.edgeOrVertexChainCrossing(d), edgeOrVertex, `${input}, EdgeOrVertexChainCrossing(d)`)
+  assert.equal(crosser.edgeOrVertexChainCrossing(c), edgeOrVertex, `${input}, EdgeOrVertexChainCrossing(c)`)
+  assert.equal(crosser.edgeOrVertexCrossing(d, c), edgeOrVertex, `${input}, EdgeOrVertexCrossing(d, c)`)
+  assert.equal(crosser.edgeOrVertexCrossing(c, d), edgeOrVertex, `${input}, EdgeOrVertexCrossing(c, d)`)
+}

s2/Point.ts

@@ -153,8 +153,8 @@ export class Point {
   /**
    * Reports whether this point is similar enough to be equal to another point.
    */
-  approxEqual(op: Point): boolean {
-    return this.vector.angle(op.vector) <= EPSILON
+  approxEqual(op: Point, eps: number = EPSILON): boolean {
+    return this.vector.angle(op.vector) <= eps
   }
 
   /**

s2/RectBounder.ts

@@ -88,7 +88,10 @@ export class RectBounder {
       lngAB = S1Interval.fullInterval()
     }
 
-    // Compute the latitude range spanned by the edge AB.
+    /**
+     * Next we compute the latitude range spanned by the edge AB. We start
+     * with the range spanning the two endpoints of the edge:
+     */
     let latAB = R1Interval.fromPoint(this.aLL.lat).addPoint(bLL.lat)
 
     // Check if AB crosses the plane through N and the Z-axis.

s2/edge_crossings.ts

@@ -67,16 +67,16 @@ export const crossingSign = (a: Point, b: Point, c: Point, d: Point): Crossing =
  * It is an error to call this method with 4 distinct vertices.
  */
 export const vertexCrossing = (a: Point, b: Point, c: Point, d: Point): boolean => {
-  if (a === b || c === d) return false
+  if (a.equals(b) || c.equals(d)) return false
 
   switch (true) {
-    case a === c:
-      return b === d || Point.orderedCCW(a.referenceDir(), d, b, a)
-    case b === d:
+    case a.equals(c):
+      return b.equals(d) || Point.orderedCCW(a.referenceDir(), d, b, a)
+    case b.equals(d):
       return Point.orderedCCW(b.referenceDir(), c, a, b)
-    case a === d:
-      return b === c || Point.orderedCCW(a.referenceDir(), c, b, a)
-    case b === c:
+    case a.equals(d):
+      return b.equals(c) || Point.orderedCCW(a.referenceDir(), c, b, a)
+    case b.equals(c):
       return Point.orderedCCW(b.referenceDir(), d, a, b)
   }