From fd41e343b96c83dc73ca863cb016e1c4cf0f3313 Mon Sep 17 00:00:00 2001 From: Peter Johnson Date: Sat, 17 Aug 2024 20:07:56 +0200 Subject: [PATCH] fix(s2): EdgeCrosser --- s2/EdgeCrosser.ts | 78 +++++++++++++++-- s2/EdgeCrosser_test.ts | 171 ++++++++++++++++++++++++++++++++++++++ s2/Point.ts | 4 +- s2/RectBounder.ts | 5 +- s2/edge_crossings.ts | 14 ++-- s2/point_measures_test.ts | 20 ++--- s2/predicates_test.ts | 23 +++-- 7 files changed, 283 insertions(+), 32 deletions(-) create mode 100644 s2/EdgeCrosser_test.ts diff --git a/s2/EdgeCrosser.ts b/s2/EdgeCrosser.ts index deac1f9..e96a772 100644 --- a/s2/EdgeCrosser.ts +++ b/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) diff --git a/s2/EdgeCrosser_test.ts b/s2/EdgeCrosser_test.ts new file mode 100644 index 0000000..643780b --- /dev/null +++ b/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)`) +} diff --git a/s2/Point.ts b/s2/Point.ts index cc11f32..027a12b 100644 --- a/s2/Point.ts +++ b/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 } /** diff --git a/s2/RectBounder.ts b/s2/RectBounder.ts index 5374802..f14e1e0 100644 --- a/s2/RectBounder.ts +++ b/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. diff --git a/s2/edge_crossings.ts b/s2/edge_crossings.ts index 5e222f5..5f61d60 100644 --- a/s2/edge_crossings.ts +++ b/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) } diff --git a/s2/point_measures_test.ts b/s2/point_measures_test.ts index b0f1efd..a0eae9a 100644 --- a/s2/point_measures_test.ts +++ b/s2/point_measures_test.ts @@ -1,5 +1,5 @@ -import { strict as assert } from 'node:assert' import { test, describe } from 'node:test' +import { equal, ok } from 'node:assert/strict' import { Point } from './Point' import { Vector } from '../r3/Vector' import { LatLng } from './LatLng' @@ -40,7 +40,7 @@ describe('s2.point_measures', () => { for (const [d, test] of tests.entries()) { const got = pointArea(test.a, test.b, test.c) - assert.ok( + ok( Math.abs(got - test.want) <= test.nearness, `${d}, PointArea(${test.a}, ${test.b}, ${test.c}), got ${got} want ${test.want}` ) @@ -54,18 +54,18 @@ describe('s2.point_measures', () => { const p1 = Point.fromVector(p0.vector.add(d1.vector.mul(1e-15)).normalize()) const p2 = Point.fromVector(p0.vector.add(d2.vector.mul(1e-15)).normalize()) const got = pointArea(p0, p1, p2) - assert.ok(got <= 0.7e-30, `PointArea(${p1}, ${p1}, ${p2}) = ${got}, want <= ${0.7e-30}`) + ok(got <= 0.7e-30, `PointArea(${p1}, ${p1}, ${p2}) = ${got}, want <= ${0.7e-30}`) const a = girardArea(p0, p1, p2) if (a > maxGirard) maxGirard = a } - assert.ok(maxGirard <= 1e-14, `maximum GirardArea = ${maxGirard}, want <= 1e-14`) + ok(maxGirard <= 1e-14, `maximum GirardArea = ${maxGirard}, want <= 1e-14`) const a = Point.fromLatLng(LatLng.fromDegrees(-45, -170)) const b = Point.fromLatLng(LatLng.fromDegrees(45, -170)) const c = Point.fromLatLng(LatLng.fromDegrees(0, -170)) const area = pointArea(a, b, c) - assert.strictEqual(area, 0.0, `PointArea(${a}, ${b}, ${c}) = ${area}, want 0.0`) + equal(area, 0.0, `PointArea(${a}, ${b}, ${c}) = ${area}, want 0.0`) }) test('PointArea - quarter hemisphere', (t) => { @@ -94,7 +94,7 @@ describe('s2.point_measures', () => { pointArea(test.a, test.c, test.d) + pointArea(test.a, test.d, test.e) + pointArea(test.a, test.e, test.b) - assert.ok( + ok( Math.abs(area - test.want) <= epsilon, `Adding up 4 quarter hemispheres with PointArea(), got ${area} want ${test.want}` ) @@ -109,7 +109,7 @@ describe('s2.point_measures', () => { const p3 = Point.fromLatLng(new LatLng(0, lng + Math.PI).normalized()) const p4 = Point.fromLatLng(new LatLng(0, lng + 5.0).normalized()) const area = pointArea(p0, p1, p2) + pointArea(p0, p2, p3) + pointArea(p0, p3, p4) + pointArea(p0, p4, p1) - assert.ok( + ok( Math.abs(area - 2 * Math.PI) <= 2e-15, `hemisphere area of ${p1}, ${p2}, ${p3}, ${p4}, ${p1} = ${area}, want ${2 * Math.PI}` ) @@ -127,12 +127,12 @@ describe('s2.point_measures', () => { for (const test of tests) { const gotAngle = angle(test.a, test.b, test.c) - assert.ok( + ok( Math.abs(gotAngle - test.wantAngle) <= epsilon, `Angle(${test.a}, ${test.b}, ${test.c}) = ${gotAngle}, want ${test.wantAngle}` ) const gotTurnAngle = turnAngle(test.a, test.b, test.c) - assert.ok( + ok( Math.abs(gotTurnAngle - test.wantTurnAngle) <= epsilon, `TurnAngle(${test.a}, ${test.b}, ${test.c}) = ${gotTurnAngle}, want ${test.wantTurnAngle}` ) @@ -144,6 +144,6 @@ describe('s2.point_measures', () => { const b = new Point(-1.706078905422188652e-1, -8.246067119418969416e-1, 5.393669607095969987e-1) const c = new Point(-1.705800600596222294e-1, -8.244634596153025408e-1, 5.395947061167500891e-1) const area = pointArea(a, b, c) - assert.strictEqual(area, 0, `PointArea(${a}, ${b}, ${c}) should have been 0, got ${area}`) + equal(area, 0, `PointArea(${a}, ${b}, ${c}) should have been 0, got ${area}`) }) }) diff --git a/s2/predicates_test.ts b/s2/predicates_test.ts index a17bfb1..f932d57 100644 --- a/s2/predicates_test.ts +++ b/s2/predicates_test.ts @@ -7,7 +7,7 @@ import { CLOCKWISE, COUNTERCLOCKWISE, INDETERMINATE } from './predicates' import { Point } from './Point' import { Vector } from '../r3/Vector' import * as matrix from './matrix3x3.ts' -import { randomFrame } from './testing.ts' +import { EARTH_RADIUS_KM, randomFrame } from './testing.ts' describe('s2.predicates', () => { test('sign', () => { @@ -108,15 +108,28 @@ describe('s2.predicates', () => { }) test('stableSign failure rate', () => { - // The Earth's mean radius in kilometers (according to NASA). - const earthRadiusKm = 6371.01 - const iters = 1000 + + // Verify that stableSign is able to handle most cases where the three + // points are as collinear as possible. (For reference, triageSign fails + // almost 100% of the time on this test.) + // + // Note that the failure rate *decreases* as the points get closer together, + // and the decrease is approximately linear. For example, the failure rate + // is 0.4% for collinear points spaced 1km apart, but only 0.0004% for + // collinear points spaced 1 meter apart. + // + // 1km spacing: < 1% (actual is closer to 0.4%) + // 10km spacing: < 10% (actual is closer to 4%) const want = 0.01 const spacing = 1.0 + // Estimate the probability that stableSign will not be able to compute + // the determinant sign of a triangle A, B, C consisting of three points + // that are as collinear as possible and spaced the given distance apart + // by counting up the times it returns Indeterminate. let failureCount = 0 - const m = Math.tan(spacing / earthRadiusKm) + const m = Math.tan(spacing / EARTH_RADIUS_KM) for (let iter = 0; iter < iters; iter++) { const f = randomFrame() const a = matrix.col(f, 0)