diff --git a/README.md b/README.md index 267af46..f6918e9 100644 --- a/README.md +++ b/README.md @@ -38,17 +38,21 @@ const { s2 } = require('s2js') s2js is a full-featured JavaScript/TypeScript port of Google's S2 Geometry library. It supports: -| Feature | Supported? | Notes | -|----------------------------------|:---------:|:------| -| **CellId, Cell** | ✅ | Encoding, decoding, neighbors | -| **Point, LatLng** | ✅ | Spherical & Planar projections | -| **Loop, Polygon** | ✅ | Polygons, potentially with holes | -| **Polyline** | ✅ | Represents linear paths | -| **RegionCoverer** | ✅ | Cover regions (loops, polys) with S2 cells | -| **GeoJSON <> S2 conversions** | ✅ | Supports all GeoJSON geometry types | -| **Polygon contains/intersects** | ✅ | Point-in-polygon & shape intersection | -| **Union, Intersection, Difference** | ✅ | Basic boolean polygon operations | -| **BigInt cell/token conversion** | ✅ | Handles full 64-bit cell IDs | +| Feature | Supported? | Notes | +| ----------------------------------- | :--------: | :----------------------------------------- | +| **CellId, Cell** | ✅ | Encoding, decoding, neighbors | +| **Point, LatLng** | ✅ | Spherical & Planar projections | +| **Loop, Polygon** | ✅ | Polygons, potentially with holes | +| **Polyline** | ✅ | Represents linear paths | +| **RegionCoverer** | ✅ | Cover regions (loops, polys) with S2 cells | +| **GeoJSON <> S2 conversions** | ✅ | Supports all GeoJSON geometry types | +| **Polygon contains/intersects** | ✅ | Point-in-polygon & shape intersection | +| **Union, Intersection, Difference** | ✅ | Basic boolean polygon operations | +| **BigInt cell/token conversion** | ✅ | Handles full 64-bit cell IDs | +| **CellIndex, ClosestCellQuery** | ✅ | Query a list of cells with various shapes | + + +TODO: ADD CLOSESTEDGEQUERY HERE ### Feature Completeness diff --git a/s2/CellIndex.ts b/s2/CellIndex.ts new file mode 100644 index 0000000..8aa4748 --- /dev/null +++ b/s2/CellIndex.ts @@ -0,0 +1,502 @@ +import type { CellID } from './cellid' +import * as cellid from './cellid' +import { SentinelCellID } from './cellid' +import { MAX_LEVEL } from './cellid_constants' +import type { CellUnion } from './CellUnion' + +/** + * A special label indicating that the ContentsIterator done is true. + */ +export const CELL_INDEX_DONE_CONTENTS = -1 + +/** + * Represents a node in the CellIndex. Cells are organized in a tree such that + * the ancestors of a given node contain that node. + */ +interface CellIndexNode { + cellID: CellID + label: number + parent: number +} + +/** + * Represents a range of leaf CellIDs. The range starts at startID (a leaf cell) + * and ends at the startID field of the next RangeNode. contents points to the + * node of the CellIndex cellTree representing the cells that overlap this range. + */ +interface RangeNode { + startID: CellID + contents: number +} + +/** + * An iterator that seeks and iterates over a set of non-overlapping leaf cell + * ranges that cover the entire sphere. The indexed (CellID, label) pairs that + * intersect the current leaf cell range can be visited using + * CellIndexContentsIterator (see below). + */ +export class CellIndexRangeIterator { + rangeNodes: RangeNode[] + pos: number + + /** + * Creates an iterator for the given CellIndex. + * The iterator is initially *unpositioned*; you must call a positioning + * method such as begin() or seek() before accessing its contents. + */ + constructor( + index: CellIndex, + private nonEmpty: boolean = false + ) { + this.rangeNodes = index.rangeNodes + this.pos = 0 + } + + /** + * Reports the CellID of the start of the current range of leaf CellIDs. + * If done is true, this returns the last possible CellID. This property means + * that most loops do not need to test done explicitly. + */ + startID(): CellID { + return this.rangeNodes[this.pos].startID + } + + /** + * Reports the non-inclusive end of the current range of leaf CellIDs. + * This assumes the iterator is not done. + */ + limitID(): CellID { + return this.rangeNodes[this.pos + 1].startID + } + + /** + * Reports if no (CellID, label) pairs intersect this range. + * Also returns true if done() is true. + */ + isEmpty(): boolean { + return this.rangeNodes[this.pos].contents === CELL_INDEX_DONE_CONTENTS + } + + /** + * Positions the iterator at the first range of leaf cells (if any). + */ + begin(): void { + this.pos = 0 + while (this.nonEmpty && this.isEmpty() && !this.done()) { + this.pos++ + } + } + + /** + * Internal prev that doesn't check nonEmpty to prevent unwanted recursion. + */ + private internalPrev(): boolean { + if (this.pos === 0) { + return false + } + this.pos-- + return true + } + + /** + * Internal prev for non-empty iterator. + */ + private nonEmptyPrev(): boolean { + while (this.internalPrev()) { + if (!this.isEmpty()) { + return true + } + } + // Return the iterator to its original position. + if (this.isEmpty() && !this.done()) { + this.next() + } + return false + } + + /** + * Positions the iterator at the previous entry and reports whether it was + * not already positioned at the beginning. + */ + prev(): boolean { + if (this.nonEmpty) { + return this.nonEmptyPrev() + } + return this.internalPrev() + } + + /** + * Advances the iterator to the next range of leaf cells. + * This assumes the iterator is not done. + */ + next(): void { + this.pos++ + while (this.nonEmpty && this.isEmpty() && !this.done()) { + this.pos++ + } + } + + /** + * Reports if advancing would leave it positioned on a valid range. If the + * value would not be valid, the positioning is not changed. + */ + advance(n: number): boolean { + // Note that the last element of rangeNodes is a sentinel value. + if (n >= this.rangeNodes.length - 1 - this.pos) { + return false + } + this.pos += n + return true + } + + /** + * Positions the iterator so that done is true. + */ + finish(): void { + // Note that the last element of rangeNodes is a sentinel value. + this.pos = this.rangeNodes.length - 1 + } + + /** + * Reports if the iterator is positioned beyond the last valid range. + */ + done(): boolean { + return this.pos >= this.rangeNodes.length - 1 + } + + /** + * Positions the iterator at the first range with startID >= target. + * Such an entry always exists as long as "target" is a valid leaf cell. + * Note that it is valid to access startID even when done is true. + */ + seek(target: CellID): void { + // Binary search to find the first range with startID > target. + let lo = 0 + let hi = this.rangeNodes.length + while (lo < hi) { + const mid = (lo + hi) >>> 1 + if (this.rangeNodes[mid].startID > target) { + hi = mid + } else { + lo = mid + 1 + } + } + // Position at the previous range (which contains target), but ensure we don't go below 0. + this.pos = Math.max(lo - 1, 0) + + // Non-empty iterator needs to find the next non-empty entry. + while (this.nonEmpty && this.isEmpty() && !this.done()) { + this.pos++ + } + } +} + +/** + * An iterator that visits the (CellID, label) pairs that cover a set of leaf + * cell ranges (see CellIndexRangeIterator). Note that when multiple leaf cell + * ranges are visited, this iterator only guarantees that each result will be + * reported at least once, i.e. duplicate values may be suppressed. If you want + * duplicate values to be reported again, be sure to call clear() first. + * + * In particular, the implementation guarantees that when multiple leaf cell + * ranges are visited in monotonically increasing order, then each (CellID, label) + * pair is reported exactly once. + */ +export class CellIndexContentsIterator { + /** + * The maximum index within the cellTree slice visited during the previous + * call to startUnion. This is used to eliminate duplicate values when + * startUnion is called multiple times. + */ + private nodeCutoff: number + + /** + * The maximum index within the cellTree visited during the current call to + * startUnion. This is used to update nodeCutoff. + */ + private nextNodeCutoff: number + + /** + * The value of startID from the previous call to startUnion. This is used to + * check whether these values are monotonically increasing. + */ + private prevStartID: CellID + + /** The cell tree from CellIndex. */ + private cellTree: CellIndexNode[] + + /** A copy of the current node in the cell tree. */ + private node: CellIndexNode + + /** + * Creates a new contents iterator. + * Note that the iterator needs to be positioned using startUnion before + * it can be safely used. + */ + constructor(index: CellIndex) { + this.cellTree = index.cellTree + this.prevStartID = 0n + this.nodeCutoff = -1 + this.nextNodeCutoff = -1 + this.node = { cellID: 0n, label: CELL_INDEX_DONE_CONTENTS, parent: -1 } + } + + /** + * Clears all state with respect to which range(s) have been visited. + */ + clear(): void { + this.prevStartID = 0n + this.nodeCutoff = -1 + this.nextNodeCutoff = -1 + this.node.label = CELL_INDEX_DONE_CONTENTS + } + + /** + * Returns the current CellID. + */ + cellID(): CellID { + return this.node.cellID + } + + /** + * Returns the current Label. + */ + label(): number { + return this.node.label + } + + /** + * Advances the iterator to the next (CellID, label) pair covered by the + * current leaf cell range. This requires the iterator to not be done. + */ + next(): void { + if (this.node.parent <= this.nodeCutoff) { + // We have already processed this node and its ancestors. + this.nodeCutoff = this.nextNodeCutoff + this.node.label = CELL_INDEX_DONE_CONTENTS + } else { + this.node = { ...this.cellTree[this.node.parent] } + } + } + + /** + * Reports if all (CellID, label) pairs have been visited. + */ + done(): boolean { + return this.node.label === CELL_INDEX_DONE_CONTENTS + } + + /** + * Positions the ContentsIterator at the first (cellID, label) pair that + * covers the given leaf cell range. Note that when multiple leaf cell ranges + * are visited using the same ContentsIterator, duplicate values may be + * suppressed. If you don't want this behavior, call clear() first. + */ + startUnion(r: CellIndexRangeIterator): void { + if (r.startID() < this.prevStartID) { + this.nodeCutoff = -1 // Can't automatically eliminate duplicates. + } + this.prevStartID = r.startID() + const contents = r.rangeNodes[r.pos].contents + if (contents <= this.nodeCutoff) { + this.node.label = CELL_INDEX_DONE_CONTENTS + } else { + this.node = { ...this.cellTree[contents] } + } + // When visiting ancestors, we can stop as soon as the node index is smaller + // than any previously visited node index. Because indexes are assigned + // using a preorder traversal, such nodes are guaranteed to have already + // been reported. + this.nextNodeCutoff = contents + } +} + +/** + * Represents a delta entry used during CellIndex build. + */ +interface BuildDelta { + startID: CellID + cellID: CellID + label: number +} + +/** + * Stores a collection of (CellID, label) pairs. + * + * The CellIDs may be overlapping or contain duplicate values. For example, a + * CellIndex could store a collection of CellUnions, where each CellUnion + * gets its own non-negative int32 label. + * + * Similar to ShapeIndex and PointIndex which map each stored element to an + * identifier, CellIndex stores a label that is typically used to map the + * results of queries back to client's specific data. + * + * To build a CellIndex where each Cell has a distinct label, call add() for each + * (CellID, label) pair, and then build() the index. For example: + * + * ```typescript + * // contents is a mapping of an identifier in my system (restaurantID, + * // vehicleID, etc) to a CellID + * const contents = new Map([...]) + * + * for (const [key, val] of contents) { + * index.add(val, key) + * } + * + * index.build() + * ``` + * + * There is also a helper method that adds all elements of CellUnion with the + * same label: + * + * ```typescript + * index.addCellUnion(cellUnion, label) + * ``` + * + * Note that the index is not dynamic; the contents of the index cannot be + * changed once it has been built. Adding more after calling build() results in + * undefined behavior of the index. + * + * There are several options for retrieving data from the index. The simplest + * is to use a built-in method such as intersectingLabels (which returns + * the labels of all cells that intersect a given target CellUnion): + * + * ```typescript + * const labels = index.intersectingLabels(targetUnion) + * ``` + * + * Alternatively, you can use a ClosestCellQuery which computes the cell(s) + * that are closest to a given target geometry. + * + * Internally, the index consists of a set of non-overlapping leaf cell ranges + * that subdivide the sphere and such that each range intersects a particular + * set of (cellID, label) pairs. + * + * Most clients should use either the methods such as visitIntersectingCells + * and intersectingLabels, or a helper such as ClosestCellQuery. + */ +export class CellIndex { + /** + * A tree of (cellID, label) pairs such that if X is an ancestor of Y, then + * X.cellID contains Y.cellID. The contents of a given range of leaf cells + * can be represented by pointing to a node of this tree. + */ + cellTree: CellIndexNode[] + + /** + * The last element of rangeNodes is a sentinel value, which is necessary + * in order to represent the range covered by the previous element. + */ + rangeNodes: RangeNode[] + + constructor() { + this.cellTree = [] + this.rangeNodes = [] + } + + /** + * Adds the given CellID and Label to the index. + */ + add(id: CellID, label: number): void { + if (label < 0) { + throw new Error('labels must be non-negative') + } + this.cellTree.push({ cellID: id, label: label, parent: -1 }) + } + + /** + * Adds all of the elements of the given CellUnion to the index with the same label. + */ + addCellUnion(cu: CellUnion, label: number): void { + if (label < 0) { + throw new Error('labels must be non-negative') + } + for (const cell of cu) { + this.add(cell, label) + } + } + + /** + * Builds the index for use. This method should only be called once. + */ + build(): void { + // To build the cell tree and leaf cell ranges, we maintain a stack of + // (CellID, label) pairs that contain the current leaf cell. This struct + // represents an instruction to push or pop a (cellID, label) pair. + // + // If label >= 0, the (cellID, label) pair is pushed on the stack. + // If CellID == SentinelCellID, a pair is popped from the stack. + // Otherwise the stack is unchanged but a rangeNode is still emitted. + + const deltas: BuildDelta[] = [] + + // Create two deltas for each (cellID, label) pair: one to add the pair to + // the stack (at the start of its leaf cell range), and one to remove it from + // the stack (at the end of its leaf cell range). + for (const node of this.cellTree) { + deltas.push({ + startID: cellid.rangeMin(node.cellID), + cellID: node.cellID, + label: node.label + }) + deltas.push({ + startID: cellid.next(cellid.rangeMax(node.cellID)), + cellID: SentinelCellID, + label: -1 + }) + } + + // We also create two special deltas to ensure that a RangeNode is emitted at + // the beginning and end of the CellID range. + deltas.push({ + startID: cellid.childBeginAtLevel(cellid.fromFace(0), MAX_LEVEL), + cellID: 0n, + label: -1 + }) + deltas.push({ + startID: cellid.childEndAtLevel(cellid.fromFace(5), MAX_LEVEL), + cellID: 0n, + label: -1 + }) + + // Sort deltas: first by startID, then in reverse order by cellID, then by label. + // This is necessary to ensure that (1) larger cells are pushed on the stack + // before smaller cells, and (2) cells are popped off the stack before any + // new cells are added. + deltas.sort((a, b) => { + if (a.startID !== b.startID) { + return a.startID < b.startID ? -1 : 1 + } + if (a.cellID !== b.cellID) { + // Reverse order by cellID. + return a.cellID > b.cellID ? -1 : 1 + } + return a.label < b.label ? -1 : a.label > b.label ? 1 : 0 + }) + + // Now walk through the deltas to build the leaf cell ranges and cell tree + // (which is essentially a permanent form of the "stack" described above). + this.cellTree = [] + this.rangeNodes = [] + let contents = -1 + + for (let i = 0; i < deltas.length; ) { + const startID = deltas[i].startID + // Process all the deltas associated with the current startID. + while (i < deltas.length && deltas[i].startID === startID) { + if (deltas[i].label >= 0) { + this.cellTree.push({ + cellID: deltas[i].cellID, + label: deltas[i].label, + parent: contents + }) + contents = this.cellTree.length - 1 + } else if (deltas[i].cellID === SentinelCellID) { + contents = this.cellTree[contents].parent + } + i++ + } + this.rangeNodes.push({ startID, contents }) + } + } +} diff --git a/s2/CellIndex_test.ts b/s2/CellIndex_test.ts new file mode 100644 index 0000000..d07c947 --- /dev/null +++ b/s2/CellIndex_test.ts @@ -0,0 +1,395 @@ +import { test, describe } from 'node:test' +import { equal, ok, deepEqual } from 'node:assert/strict' +import { CellIndex, CellIndexRangeIterator, CellIndexContentsIterator, type CellIndexNode } from './CellIndex' +import type { CellID } from './cellid' +import * as cellid from './cellid' +import { MAX_LEVEL } from './cellid_constants' +import { CellUnion } from './CellUnion' +import { randomCellIDForLevel, randomUniformInt, skewedInt } from './testing' + +/** + * Represents a test input for CellIndex consisting of a cell string and label. + */ +interface CellIndexTestInput { + cellID: string + label: number +} + +/** + * Reports whether this node is less than the other for sorting purposes. + * Only compares cellID and label, not parent (which is implementation detail). + */ +const nodeLess = (a: CellIndexNode, b: CellIndexNode): number => { + if (a.cellID !== b.cellID) { + return a.cellID < b.cellID ? -1 : 1 + } + if (a.label !== b.label) { + return a.label < b.label ? -1 : 1 + } + return 0 +} + +/** + * Creates a copy of the nodes so that sorting and other tests don't alter the instance in a given CellIndex. + */ +const copyCellIndexNodes = (nodes: CellIndexNode[]): CellIndexNode[] => { + return nodes.map((n) => ({ cellID: n.cellID, label: n.label, parent: n.parent })) +} + +/** + * Reports whether two sorted arrays of nodes are equal. + * Only compares cellID and label, not parent (which is implementation detail). + */ +const cellIndexNodesEqual = (a: CellIndexNode[], b: CellIndexNode[]): boolean => { + const aSorted = copyCellIndexNodes(a).sort(nodeLess) + const bSorted = copyCellIndexNodes(b).sort(nodeLess) + + if (aSorted.length !== bSorted.length) return false + + for (let i = 0; i < aSorted.length; i++) { + if (aSorted[i].cellID !== bSorted[i].cellID) return false + if (aSorted[i].label !== bSorted[i].label) return false + } + + return true +} + +/** + * Serializes CellIndexNode arrays to a string, handling BigInt values. + */ +const stringifyNodes = (nodes: CellIndexNode[]): string => { + return JSON.stringify(nodes, (_, v) => (typeof v === 'bigint' ? v.toString() : v)) +} + +/** + * Verifies that positioning the CellIndexRangeIterator to finished results in done() being true. + */ +const verifyCellIndexCellIterator = (_desc: string, _index: CellIndex): void => { + // TODO: Once the plain iterator is implemented, add this check. + /* + const actual: CellIndexNode[] = [] + const iter = new CellIndexIterator(index) + for (iter.begin(); !iter.done(); iter.next()) { + actual.push({ cellID: iter.startID(), label: iter.label(), parent: -1 }) + } + const want = copyCellIndexNodes(index.cellTree) + ok(cellIndexNodesEqual(actual, want), `${desc}: cellIndexNodes not equal but should be`) + */ +} + +/** + * Verifies that the CellIndexRangeIterator and CellIndexNonEmptyRangeIterator work correctly. + */ +const verifyCellIndexRangeIterators = (desc: string, index: CellIndex): void => { + // Precondition: The index has been built. + // Under test: The range iterator correctly positions and navigates. + // Postcondition: All iterator operations behave correctly. + + // Tests finish(), which is not otherwise tested below. + const it = new CellIndexRangeIterator(index) + it.begin() + it.finish() + ok(it.done(), `${desc}: positioning iterator to finished should be done, but was not`) + + // And also for non-empty ranges. + const nonEmpty = new CellIndexRangeIterator(index, true) + nonEmpty.begin() + nonEmpty.finish() + ok(nonEmpty.done(), `${desc}: positioning non-empty iterator to finished should be done, but was not`) + + // Iterate through all the ranges in the index. We simultaneously iterate + // through the non-empty ranges and check that the correct ranges are found. + let prevStart: CellID = 0n + let nonEmptyPrevStart: CellID = 0n + it.begin() + nonEmpty.begin() + + while (!it.done()) { + // Check that seeking in the current range takes us to this range. + const it2 = new CellIndexRangeIterator(index) + const start = it.startID() + it2.seek(it.startID()) + equal(it2.startID(), start, `${desc}: it2.startID() = ${it2.startID()}, want ${start}`) + + it2.seek(cellid.prev(it.limitID())) + equal(it2.startID(), start, `${desc}: it2.seek(${cellid.prev(it.limitID())}) = ${it2.startID()}, want ${start}`) + + // And also for non-empty ranges. + const nonEmpty2 = new CellIndexRangeIterator(index, true) + const nonEmptyStart = nonEmpty.startID() + nonEmpty2.seek(it.startID()) + equal( + nonEmpty2.startID(), + nonEmptyStart, + `${desc}: nonEmpty2.startID() = ${nonEmpty2.startID()}, want ${nonEmptyStart}` + ) + + nonEmpty2.seek(cellid.prev(it.limitID())) + equal( + nonEmpty2.startID(), + nonEmptyStart, + `${desc}: nonEmpty2.startID() = ${nonEmpty2.startID()}, want ${nonEmptyStart}` + ) + + // Test prev() and next(). + if (it2.prev()) { + equal(it2.startID(), prevStart, `${desc}: it2.startID() = ${it2.startID()}, want ${prevStart}`) + it2.next() + equal(it2.startID(), start, `${desc}: it2.startID() = ${it2.startID()}, want ${start}`) + } else { + equal(it2.startID(), start, `${desc}: it2.startID() = ${it2.startID()}, want ${start}`) + equal(prevStart, 0n, `${desc}: prevStart = ${prevStart}, want 0n`) + } + + // And also for non-empty ranges. + if (nonEmpty2.prev()) { + equal( + nonEmpty2.startID(), + nonEmptyPrevStart, + `${desc}: nonEmpty2.startID() = ${nonEmpty2.startID()}, want ${nonEmptyPrevStart}` + ) + nonEmpty2.next() + equal( + nonEmpty2.startID(), + nonEmptyStart, + `${desc}: nonEmpty2.startID() = ${nonEmpty2.startID()}, want ${nonEmptyStart}` + ) + } else { + equal( + nonEmpty2.startID(), + nonEmptyStart, + `${desc}: nonEmpty2.startID() = ${nonEmpty2.startID()}, want ${nonEmptyStart}` + ) + equal(nonEmptyPrevStart, 0n, `${desc}: nonEmptyPrevStart = ${nonEmptyPrevStart}, want 0n`) + } + + // Keep the non-empty iterator synchronized with the regular one. + if (!it.isEmpty()) { + equal(it.startID(), nonEmpty.startID(), `${desc}: it.startID() = ${it.startID()}, want ${nonEmpty.startID()}`) + equal(it.limitID(), nonEmpty.limitID(), `${desc}: it.limitID() = ${it.limitID()}, want ${nonEmpty.limitID()}`) + ok(!nonEmpty.done(), `${desc}: nonEmpty iterator should not be done but was`) + nonEmptyPrevStart = nonEmptyStart + nonEmpty.next() + } + + prevStart = start + it.next() + } + + // Verify that the NonEmptyRangeIterator is also finished. + ok(nonEmpty.done(), `${desc}: non empty iterator should have also finished`) +} + +/** + * Verifies that RangeIterator and ContentsIterator can be used to determine + * the exact set of (CellID, label) pairs that contain any leaf cell. + */ +const verifyCellIndexContents = (desc: string, index: CellIndex): void => { + // Precondition: The index has been built. + // Under test: The contents iterator returns the correct set of (cellID, label) pairs. + // Postcondition: All ranges are verified. + + // minCellID is the first CellID that has not been validated yet. + let minCellID = cellid.childBeginAtLevel(cellid.fromFace(0), MAX_LEVEL) + const r = new CellIndexRangeIterator(index) + + for (r.begin(); !r.done(); r.next()) { + equal( + r.startID(), + minCellID, + `${desc}: minCellID should match the previous ending cellID. got ${r.startID()}, want ${minCellID}` + ) + ok( + minCellID < r.limitID(), + `${desc}: minCellID should be < the end of the current range. got ${r.limitID()}, want > ${minCellID}` + ) + ok(cellid.isLeaf(r.limitID()) || r.done(), `${desc}: ending range cell ID should be a leaf or done`) + + minCellID = r.limitID() + + // Build a list of expected (CellID, label) for this range. + const expected: CellIndexNode[] = [] + for (const x of index.cellTree) { + // The cell contains the entire range. + if (cellid.rangeMin(x.cellID) <= r.startID() && cellid.next(cellid.rangeMax(x.cellID)) >= r.limitID()) { + expected.push(x) + } else { + // Verify that the cell does not intersect the range. + if (cellid.rangeMin(x.cellID) <= cellid.prev(r.limitID()) && cellid.rangeMax(x.cellID) >= r.startID()) { + ok( + false, + `${desc}: CellID does not intersect the current range: ${cellid.rangeMin(x.cellID)} <= ${cellid.prev(r.limitID())} && ${cellid.rangeMax(x.cellID)} >= ${r.startID()}` + ) + } + } + } + + const actual: CellIndexNode[] = [] + const cIter = new CellIndexContentsIterator(index) + for (cIter.startUnion(r); !cIter.done(); cIter.next()) { + actual.push({ cellID: cIter.cellID(), label: cIter.label(), parent: -1 }) + } + + ok( + cellIndexNodesEqual(expected, actual), + `${desc}: comparing contents iterator contents to this range: got ${stringifyNodes(actual)}, want ${stringifyNodes(expected)}` + ) + } + + equal( + minCellID, + cellid.childEndAtLevel(cellid.fromFace(5), MAX_LEVEL), + `${desc}: the final cell should be the sentinel value, got ${minCellID}` + ) +} + +/** + * Verifies that the index computes the correct set of (cellID, label) pairs + * for every possible leaf cell. The running time of this function is + * quadratic in the size of the index. + */ +const cellIndexQuadraticValidate = (desc: string, index: CellIndex): void => { + // Precondition: The index has cells added. + // Under test: The index correctly builds and iterators work. + // Postcondition: All validations pass. + + index.build() + verifyCellIndexCellIterator(desc, index) + verifyCellIndexRangeIterators(desc, index) + verifyCellIndexContents(desc, index) +} + +/** + * Generates a random CellUnion with the specified number of cells. + */ +const randomCellUnion = (numCells: number): CellUnion => { + const cu = new CellUnion() + for (let i = 0; i < numCells; i++) { + cu.push(randomCellIDForLevel(randomUniformInt(MAX_LEVEL + 1))) + } + cu.normalize() + return cu +} + +describe('s2.CellIndex', () => { + test('empty index', () => { + // Precondition: An empty CellIndex is created. + // Under test: The empty index validates correctly. + // Postcondition: All iterator operations work on an empty index. + + const index = new CellIndex() + cellIndexQuadraticValidate('Empty', index) + }) + + test('one face cell', () => { + // Precondition: A CellIndex with a single face cell is created. + // Under test: The index correctly handles a single face-level cell. + // Postcondition: The index validates correctly. + + const index = new CellIndex() + index.add(cellid.fromString('0/'), 0) + cellIndexQuadraticValidate('One face cell', index) + }) + + test('one leaf cell', () => { + // Precondition: A CellIndex with a single leaf cell is created. + // Under test: The index correctly handles a single max-level cell. + // Postcondition: The index validates correctly. + + const index = new CellIndex() + index.add(cellid.fromString('1/012301230123012301230123012301'), 12) + cellIndexQuadraticValidate('One Leaf Cell', index) + }) + + test('duplicate values', () => { + // Precondition: A CellIndex with duplicate cells and different labels is created. + // Under test: The index correctly handles duplicate cell entries. + // Postcondition: The index validates correctly with all duplicates. + + const index = new CellIndex() + index.add(cellid.fromString('0/'), 0) + index.add(cellid.fromString('0/'), 0) + index.add(cellid.fromString('0/'), 1) + index.add(cellid.fromString('0/'), 17) + cellIndexQuadraticValidate('Duplicate Values', index) + }) + + test('disjoint cells', () => { + // Precondition: A CellIndex with non-overlapping cells is created. + // Under test: The index correctly handles disjoint cells. + // Postcondition: The index validates correctly. + + const index = new CellIndex() + index.add(cellid.fromString('0/'), 0) + index.add(cellid.fromString('3/'), 0) + cellIndexQuadraticValidate('Disjoint Cells', index) + }) + + test('nested cells', () => { + // Precondition: A CellIndex with nested cells is created, where some cells have the same RangeMin or RangeMax. + // Under test: The index correctly handles nested cells with randomly ordered labels. + // Postcondition: The index validates correctly. + + const index = new CellIndex() + const inputs: CellIndexTestInput[] = [ + { cellID: '1/', label: 3 }, + { cellID: '1/0', label: 15 }, + { cellID: '1/000', label: 9 }, + { cellID: '1/00000', label: 11 }, + { cellID: '1/012', label: 6 }, + { cellID: '1/01212', label: 5 }, + { cellID: '1/312', label: 17 }, + { cellID: '1/31200', label: 4 }, + { cellID: '1/3120000', label: 10 }, + { cellID: '1/333', label: 20 }, + { cellID: '1/333333', label: 18 }, + { cellID: '5/', label: 3 }, + { cellID: '5/3', label: 31 }, + { cellID: '5/3333', label: 27 } + ] + + for (const input of inputs) { + index.add(cellid.fromString(input.cellID), input.label) + } + cellIndexQuadraticValidate('Nested Cells', index) + }) + + test('contents iterator suppresses duplicates', () => { + // Precondition: A CellIndex with nested cells that would cause duplicates is created. + // Under test: The contents iterator stops reporting values once it reaches a node of the cell tree that was visited by the previous call to begin(). + // Postcondition: The index validates correctly. + + const index = new CellIndex() + const inputs: CellIndexTestInput[] = [ + { cellID: '2/1', label: 1 }, + { cellID: '2/1', label: 2 }, + { cellID: '2/10', label: 3 }, + { cellID: '2/100', label: 4 }, + { cellID: '2/102', label: 5 }, + { cellID: '2/1023', label: 6 }, + { cellID: '2/31', label: 7 }, + { cellID: '2/313', label: 8 }, + { cellID: '2/3132', label: 9 }, + { cellID: '3/1', label: 10 }, + { cellID: '3/12', label: 11 }, + { cellID: '3/13', label: 12 } + ] + + for (const input of inputs) { + index.add(cellid.fromString(input.cellID), input.label) + } + cellIndexQuadraticValidate('Contents Iterator Suppresses Duplicates', index) + }) + + test('random cell unions', () => { + // Precondition: A CellIndex with 100 random CellUnions is created. + // Under test: The index correctly handles overlapping cell unions with distinct labels. + // Postcondition: The index validates correctly. + + const index = new CellIndex() + for (let i = 0; i < 100; i++) { + index.addCellUnion(randomCellUnion(10), i) + } + cellIndexQuadraticValidate('Random Cell Unions', index) + }) +}) diff --git a/s2/ClosestCellQuery.ts b/s2/ClosestCellQuery.ts new file mode 100644 index 0000000..0486730 --- /dev/null +++ b/s2/ClosestCellQuery.ts @@ -0,0 +1,463 @@ +/** + * S2ClosestCellQuery is a helper class for finding the closest cell(s) to a + * given point, edge, S2Cell, S2CellUnion, or geometry collection. A typical + * use case would be to add a collection of S2Cell coverings to an S2CellIndex + * (representing a collection of original geometry), and then use + * S2ClosestCellQuery to find all coverings that are within a given distance + * of some target geometry (which could be represented exactly, or could also + * be a covering). The distance to the original geometry corresponding to + * each covering could then be measured more precisely if desired. + * + * For example, here is how to find all cells that are closer than + * "distanceLimit" to a given target point: + * + * ```typescript + * const query = new ClosestCellQuery(cellIndex) + * query.options.maxDistance = distanceLimit + * const target = new ClosestCellQuery.PointTarget(targetPoint) + * for (const result of query.findClosestCells(target)) { + * // result.distance is the distance to the target. + * // result.cellID is the indexed S2CellId. + * // result.label is the integer label associated with the S2CellId. + * doSomething(targetPoint, result) + * } + * ``` + * + * You can find either the k closest cells, or all cells within a given + * radius, or both (i.e., the k closest cells up to a given maximum radius). + * By default *all* cells are returned, so you should always specify either + * maxResults or maxDistance or both. You can also restrict the results + * to cells that intersect a given S2Region. + * + * There is a findClosestCell() convenience method that returns the closest + * cell. However, if you only need to test whether the distance is above or + * below a given threshold (e.g., 10 km), it is typically much faster to use + * the isDistanceLess() method instead. Unlike findClosestCell(), this method + * stops as soon as it can prove that the minimum distance is either above or + * below the threshold. + * + * @module ClosestCellQuery + */ + +import type { CellID } from './cellid' +import type { ChordAngle } from '../s1/chordangle' +import type { CellIndex } from './CellIndex' +import type { Region } from './Region' +import type { ShapeIndex } from './ShapeIndex' + +import { Cell } from './Cell' +import { CellUnion } from './CellUnion' +import { Point } from './Point' +import * as chordangle from '../s1/chordangle' +import { updateMinDistance } from './edge_distances' +import { CellIndexRangeIterator, CellIndexContentsIterator } from './CellIndex' +import { ShapeIndexTarget as ClosestEdgeShapeIndexTarget } from './ClosestEdgeQuery' + +/** + * Represents a closest (cellID, label) pair result from the query. + */ +export interface ClosestCellQueryResult { + /** The distance from the target to this cell. */ + distance: ChordAngle + + /** The indexed S2CellId. */ + cellID: CellID + + /** The integer label associated with the S2CellId. */ + label: number +} + +/** + * Returns an empty result indicating no cell was found. + */ +function emptyClosestCellQueryResult(): ClosestCellQueryResult { + return { + distance: chordangle.infChordAngle(), + cellID: 0n, + label: -1 + } +} + +/** + * Reports whether the result is empty (no cell was found). + */ +function isEmptyResult(result: ClosestCellQueryResult): boolean { + return result.label < 0 +} + +/** + * Options that control the set of cells returned. Note that by default + * *all* cells are returned, so you will always want to set either the + * maxResults option or the maxDistance option (or both). + */ +export class ClosestCellQueryOptions { + /** + * Specifies that at most "maxResults" cells should be returned. + * Default: Infinity (no limit) + */ + maxResults: number = Infinity + + /** + * Specifies that only cells whose distance to the target is less than + * "maxDistance" should be returned. + * + * Note that cells whose distance is exactly equal to "maxDistance" are + * not returned. Normally this doesn't matter, because distances are not + * computed exactly in the first place, but if such cells are needed then + * see setInclusiveMaxDistance() below. + * + * Default: Infinity + */ + maxDistance: ChordAngle = chordangle.infChordAngle() + + /** + * Specifies that cells up to maxError further away than the true + * closest cells may be substituted in the result set, as long as such + * cells satisfy all the remaining search criteria (such as maxDistance). + * This option only has an effect if maxResults is also specified; + * otherwise all cells closer than maxDistance will always be returned. + * + * Default: 0 + */ + maxError: ChordAngle = 0 + + /** + * If specified, then only cells that intersect the given region are + * returned. This can be used to restrict the results to cells that + * intersect a given S2LatLngRect, for example. + * + * Default: undefined (no region restriction) + */ + region?: Region + + /** + * Sets maxDistance to the given value such that cells whose distance + * is exactly equal to maxDistance are also returned. + * Equivalent to setting maxDistance to maxDistance.successor(). + */ + setInclusiveMaxDistance(maxDistance: ChordAngle): void { + this.maxDistance = chordangle.successor(maxDistance) + } + + /** + * Sets maxDistance such that cells whose true distance is less than + * or equal to maxDistance will be returned (along with some cells whose + * true distance is slightly greater). + * + * This ensures that all cells whose true distance is less than or equal + * to maxDistance will be returned. The maxDistance is increased by the + * maximum error in the distance calculation. + */ + setConservativeMaxDistance(maxDistance: ChordAngle): void { + this.maxDistance = chordangle.successor(chordangle.expanded(maxDistance, this.computeMaxDistanceError(maxDistance))) + } + + /** + * Computes the maximum error for the given distance. + */ + private computeMaxDistanceError(distance: ChordAngle): number { + // The maximum error depends on the distance computation method. + // For now, return a conservative estimate. + return chordangle.maxPointError(distance) + } +} + +/** + * Target represents the geometry to which the distance is measured. + * This is the base interface for all target types. + */ +export interface ClosestCellQueryTarget { + /** + * Returns the maximum number of cells in the index for which it is + * faster to use brute force search rather than the hierarchical method. + */ + maxBruteForceIndexSize(): number + + /** + * Returns the distance from the target to the given cell. + */ + distanceToCell(cell: Cell): ChordAngle + + /** + * Returns the distance from the target to the given point. + */ + distanceToPoint(point: Point): ChordAngle +} + +/** + * Target subtype that computes the closest distance to a point. + */ +export class PointTarget implements ClosestCellQueryTarget { + readonly point: Point + + constructor(point: Point) { + this.point = point + } + + maxBruteForceIndexSize(): number { + return 100 + } + + distanceToCell(cell: Cell): ChordAngle { + return cell.distance(this.point) + } + + distanceToPoint(point: Point): ChordAngle { + return Point.chordAngleBetweenPoints(this.point, point) + } +} + +/** + * Target subtype that computes the closest distance to an edge. + */ +export class EdgeTarget implements ClosestCellQueryTarget { + constructor( + readonly a: Point, + readonly b: Point + ) {} + + maxBruteForceIndexSize(): number { + return 100 + } + + distanceToCell(cell: Cell): ChordAngle { + return cell.distanceToEdge(this.a, this.b) + } + + distanceToPoint(point: Point): ChordAngle { + return updateMinDistance(point, this.a, this.b, chordangle.infChordAngle()).dist + } +} + +/** + * Target subtype that computes the closest distance to an S2Cell + * (including the interior of the cell). + */ +export class CellTarget implements ClosestCellQueryTarget { + constructor(readonly cell: Cell) {} + + maxBruteForceIndexSize(): number { + return 100 + } + + distanceToCell(cell: Cell): ChordAngle { + return this.cell.distanceToCell(cell) + } + + distanceToPoint(point: Point): ChordAngle { + return this.cell.distance(point) + } +} + +/** + * Target subtype that computes the closest distance to an S2CellUnion. + */ +export class CellUnionTarget implements ClosestCellQueryTarget { + constructor(readonly cellUnion: CellUnion) {} + + maxBruteForceIndexSize(): number { + return 100 + } + + distanceToCell(cell: Cell): ChordAngle { + let minDist = chordangle.infChordAngle() + for (const id of this.cellUnion) { + const unionCell = Cell.fromCellID(id) + const dist = unionCell.distanceToCell(cell) + if (dist < minDist) { + minDist = dist + } + if (minDist === 0) break + } + return minDist + } + + distanceToPoint(point: Point): ChordAngle { + let minDist = chordangle.infChordAngle() + for (const id of this.cellUnion) { + const cell = Cell.fromCellID(id) + const dist = cell.distance(point) + if (dist < minDist) { + minDist = dist + } + if (minDist === 0) break + } + return minDist + } +} + +/** + * Target subtype that computes the closest distance to an S2ShapeIndex + * (an arbitrary collection of points, polylines, and/or polygons). + * + * By default, distances are measured to the boundary and interior of + * polygons in the S2ShapeIndex rather than to polygon boundaries only. + * If you wish to change this behavior, you may call: + * + * target.setIncludeInteriors(false) + */ +export class ShapeIndexTarget implements ClosestCellQueryTarget { + private edgeTarget: ClosestEdgeShapeIndexTarget + + constructor(readonly index: ShapeIndex) { + this.edgeTarget = new ClosestEdgeShapeIndexTarget(index) + } + + /** + * Specifies whether distances should be measured to the boundary and + * interior of polygons (true) or only to polygon boundaries (false). + * The default is true. + */ + setIncludeInteriors(includeInteriors: boolean): void { + this.edgeTarget.setIncludeInteriors(includeInteriors) + } + + maxBruteForceIndexSize(): number { + return this.edgeTarget.maxBruteForceIndexSize() + } + + distanceToCell(cell: Cell): ChordAngle { + const result = this.edgeTarget.updateMinDistanceToCell(cell, chordangle.infChordAngle()) + return result.distance + } + + distanceToPoint(point: Point): ChordAngle { + const result = this.edgeTarget.updateMinDistanceToPoint(point, chordangle.infChordAngle()) + return result.distance + } +} + +/** + * S2ClosestCellQuery is a helper class for finding the closest cell(s) to a + * given point, edge, S2Cell, S2CellUnion, or geometry collection. + */ +export class ClosestCellQuery { + /** + * Constructs a new ClosestCellQuery for the given CellIndex. + * Options may be specified here or changed at any time using the options property. + * + */ + constructor( + private readonly index: CellIndex, + private options: ClosestCellQueryOptions = new ClosestCellQueryOptions() + ) {} + + + /** + * Find the closest cells to the target. + * @param target - The target to find the closest cells to. + * @returns An array of ClosestCellQueryResult objects. + */ + findClosestCells(target: ClosestCellQueryTarget): ClosestCellQueryResult[] { + // Iterate through all cells in the index and compute distances. + const rangeIter = new CellIndexRangeIterator(this.index, true) + const contentsIter = new CellIndexContentsIterator(this.index) + + // Collect all candidate cells. + const candidates: Array<{ cellID: CellID; label: number; distance: ChordAngle }> = [] + + for (rangeIter.begin(); !rangeIter.done(); rangeIter.next()) { + contentsIter.startUnion(rangeIter) + while (!contentsIter.done()) { + const cellID = contentsIter.cellID() + const label = contentsIter.label() + const cell = Cell.fromCellID(cellID) + + // Check if this cell passes the region filter. + if (this.options.region && !this.options.region.intersectsCell(cell)) { + contentsIter.next() + continue + } + + const distance = target.distanceToCell(cell) + + // Check distance constraint. + if (distance < this.options.maxDistance) { + candidates.push({ cellID, label, distance }) + } + + contentsIter.next() + } + } + + // Sort by distance. + candidates.sort((a, b) => a.distance - b.distance) + + // Apply maxResults limit. + const limit = Math.min(candidates.length, this.options.maxResults) + return candidates.slice(0, limit) + } + + /** + * Returns the closest cell to the target. If no cell satisfies the search + * criteria, then the Result object will have distance == Infinity and + * isEmpty() == true. + */ + findClosestCell(target: ClosestCellQueryTarget): ClosestCellQueryResult { + const savedMaxResults = this.options.maxResults + this.options.maxResults = 1 + + const results = this.findClosestCells(target) + + this.options.maxResults = savedMaxResults + + if (results.length === 0) { + return emptyClosestCellQueryResult() + } + return results[0] + } + + /** + * Returns the minimum distance to the target. If the index or target is + * empty, returns Infinity. + * + * Use isDistanceLess() if you only want to compare the distance against a + * threshold value, since it is often much faster. + */ + getDistance(target: ClosestCellQueryTarget): ChordAngle { + return this.findClosestCell(target).distance + } + + /** + * Returns true if the distance to "target" is less than "limit". + * + * This method is usually much faster than getDistance(), since it is much + * less work to determine whether the minimum distance is above or below a + * threshold than it is to calculate the actual minimum distance. + */ + isDistanceLess(target: ClosestCellQueryTarget, limit: ChordAngle): boolean { + const savedMaxResults = this.options.maxResults + const savedMaxDistance = this.options.maxDistance + + this.options.maxResults = 1 + this.options.maxDistance = limit + + const result = this.findClosestCell(target) + + this.options.maxResults = savedMaxResults + this.options.maxDistance = savedMaxDistance + + return !isEmptyResult(result) + } + + /** + * Like isDistanceLess(), but also returns true if the distance to "target" + * is exactly equal to "limit". + */ + isDistanceLessOrEqual(target: ClosestCellQueryTarget, limit: ChordAngle): boolean { + return this.isDistanceLess(target, chordangle.successor(limit)) + } + + /** + * Like isDistanceLessOrEqual(), except that "limit" is increased by the + * maximum error in the distance calculation. This ensures that this + * function returns true whenever the true, exact distance is less than + * or equal to "limit". + */ + isConservativeDistanceLessOrEqual(target: ClosestCellQueryTarget, limit: ChordAngle): boolean { + return this.isDistanceLess( + target, + chordangle.successor(chordangle.expanded(limit, chordangle.maxPointError(limit))) + ) + } +} diff --git a/s2/ClosestCellQuery_test.ts b/s2/ClosestCellQuery_test.ts new file mode 100644 index 0000000..caa8f4d --- /dev/null +++ b/s2/ClosestCellQuery_test.ts @@ -0,0 +1,507 @@ +import { test, describe } from 'node:test' +import { equal, ok, deepEqual } from 'node:assert/strict' +import { CellIndex } from './CellIndex' +import { + ClosestCellQuery, + ClosestCellQueryOptions, + PointTarget, + EdgeTarget, + CellTarget, + CellUnionTarget, + ShapeIndexTarget +} from './ClosestCellQuery' +import { Cell } from './Cell' +import { CellUnion } from './CellUnion' +import { LatLng } from './LatLng' +import { Loop } from './Loop' +import { Point } from './Point' +import { Polygon } from './Polygon' +import { ShapeIndex } from './ShapeIndex' +import * as cellid from './cellid' +import * as chordangle from '../s1/chordangle' +import { parsePoint } from './testing_textformat' +import { randomCellID, samplePointFromCap, randomPoint } from './testing' +import { Cap } from './Cap' +import { RegionCoverer } from './RegionCoverer' + +/** + * Creates a rectangular S2Polygon from lat/lng bounds in degrees. + * The loop is oriented counter-clockwise (interior on the left). + */ +const makeRectPolygon = ( + bottomLeftLat: number, + bottomLeftLng: number, + topRightLat: number, + topRightLng: number +): Polygon => { + const bottomLeft = Point.fromLatLng(LatLng.fromDegrees(bottomLeftLat, bottomLeftLng)) + const bottomRight = Point.fromLatLng(LatLng.fromDegrees(bottomLeftLat, topRightLng)) + const topRight = Point.fromLatLng(LatLng.fromDegrees(topRightLat, topRightLng)) + const topLeft = Point.fromLatLng(LatLng.fromDegrees(topRightLat, bottomLeftLng)) + // CCW order: bottomLeft -> topLeft -> topRight -> bottomRight + const loop = new Loop([bottomLeft, bottomRight, topRight, topLeft]) + return new Polygon([loop]) +} + +/** + * Helper to create a CellID from a lat:lng string. + */ +const cellIDFromLatLng = (str: string): bigint => { + const point = parsePoint(str) + return cellid.fromPoint(point) +} + +/** + * Helper to create a CellID from a cell string like "4/012". + */ +const cellIDFromString = (s: string): bigint => { + return cellid.fromString(s) +} + +describe('S2ClosestCellQuery', () => { + test('NoCells', () => { + // Precondition: An empty CellIndex is built. + const index = new CellIndex() + index.build() + + // Under test: Query for the closest cell returns an empty result. + const query = new ClosestCellQuery(index) + const target = new PointTarget(new Point(1, 0, 0)) + const result = query.findClosestCell(target) + + // Postcondition: Result indicates no cell was found. + equal(result.distance, chordangle.infChordAngle()) + equal(result.cellID, 0n) + equal(result.label, -1) + equal(query.getDistance(target), chordangle.infChordAngle()) + }) + + test('OptionsNotModified', () => { + // Precondition: An index with 3 cells and specific query options. + const options = new ClosestCellQueryOptions() + options.maxResults = 3 + options.maxDistance = chordangle.fromAngle(3 * (Math.PI / 180)) // 3 degrees + options.maxError = chordangle.fromAngle(0.001 * (Math.PI / 180)) // 0.001 degrees + + const index = new CellIndex() + index.add(cellIDFromLatLng('1:1'), 1) + index.add(cellIDFromLatLng('1:2'), 2) + index.add(cellIDFromLatLng('1:3'), 3) + index.build() + + const query = new ClosestCellQuery(index, options) + const target = new PointTarget(parsePoint('2:2')) + + // Under test: Query methods do not modify options. + const closestResult = query.findClosestCell(target) + equal(closestResult.label, 2) + + const distance = query.getDistance(target) + const distanceDegrees = chordangle.angle(distance) * (180 / Math.PI) + ok(Math.abs(distanceDegrees - 1.0) < 0.1, `Distance should be approximately 1 degree, got ${distanceDegrees}`) + + ok(query.isDistanceLess(target, chordangle.fromAngle(1.5 * (Math.PI / 180)))) + + // Postcondition: Options remain unchanged. + equal(query.options.maxResults, 3) + equal(query.options.maxDistance, options.maxDistance) + equal(query.options.maxError, options.maxError) + }) + + test('DistanceEqualToLimit', () => { + // Precondition: An index with one cell. + const id0 = cellIDFromLatLng('23:12') + const id1 = cellIDFromLatLng('47:11') + + const index = new CellIndex() + index.add(id0, 0) + index.build() + + const query = new ClosestCellQuery(index) + + // Under test: Distance comparison methods with identical cells (zero distance). + const target0 = new CellTarget(Cell.fromCellID(id0)) + const dist0 = 0 // Zero distance + + // Postcondition: isDistanceLess returns false for zero, but isDistanceLessOrEqual returns true. + ok(!query.isDistanceLess(target0, dist0), 'isDistanceLess should return false for exact distance') + ok(query.isDistanceLessOrEqual(target0, dist0), 'isDistanceLessOrEqual should return true for exact distance') + ok(query.isConservativeDistanceLessOrEqual(target0, dist0), 'isConservativeDistanceLessOrEqual should return true') + + // Under test: Distance comparison with non-zero distance between different cells. + const target1 = new CellTarget(Cell.fromCellID(id1)) + const cell0 = Cell.fromCellID(id0) + const cell1 = Cell.fromCellID(id1) + const dist1 = cell0.distanceToCell(cell1) + + // Postcondition: Same pattern for non-zero distance. + ok(!query.isDistanceLess(target1, dist1), 'isDistanceLess should return false for exact distance') + ok(query.isDistanceLessOrEqual(target1, dist1), 'isDistanceLessOrEqual should return true for exact distance') + ok(query.isConservativeDistanceLessOrEqual(target1, dist1), 'isConservativeDistanceLessOrEqual should return true') + }) + + test('TargetPointInsideIndexedCell', () => { + // Precondition: An index with one cell. + const cellId = cellIDFromString('4/012') + + const index = new CellIndex() + index.add(cellId, 1) + index.build() + + const query = new ClosestCellQuery(index) + + // Under test: Query for a point inside the indexed cell. + const target = new PointTarget(cellid.point(cellId)) + const result = query.findClosestCell(target) + + // Postcondition: Distance is zero, and the correct cell is returned. + equal(result.distance, 0) + equal(result.cellID, cellId) + equal(result.label, 1) + }) + + test('EmptyTargetOptimized', () => { + // Precondition: An index with many random cells. + const index = new CellIndex() + for (let i = 0; i < 1000; i++) { + index.add(randomCellID(), i) + } + index.build() + + const query = new ClosestCellQuery(index) + query.options.maxDistance = 1e-5 // Very small radius + + // Under test: Query with an empty ShapeIndex target. + const targetIndex = new ShapeIndex() + const target = new ShapeIndexTarget(targetIndex) + + const results = query.findClosestCells(target) + + // Postcondition: No results are returned. + equal(results.length, 0) + }) + + test('EmptyCellUnionTarget', () => { + // Precondition: An empty CellUnion target. + const target = new CellUnionTarget(new CellUnion()) + + // Under test: Query against empty index. + const emptyIndex = new CellIndex() + emptyIndex.build() + const emptyQuery = new ClosestCellQuery(emptyIndex) + + // Postcondition: Returns infinity distance. + equal(emptyQuery.getDistance(target), chordangle.infChordAngle()) + + // Under test: Query against index with one cell. + const oneCellIndex = new CellIndex() + oneCellIndex.add(cellIDFromString('1/123123'), 1) + oneCellIndex.build() + const oneCellQuery = new ClosestCellQuery(oneCellIndex) + + // Postcondition: Returns infinity distance for empty target. + equal(oneCellQuery.getDistance(target), chordangle.infChordAngle()) + }) + + test('FindClosestCellsWithMaxResults', () => { + // Precondition: An index with multiple cells. + const index = new CellIndex() + index.add(cellIDFromLatLng('0:0'), 0) + index.add(cellIDFromLatLng('0:1'), 1) + index.add(cellIDFromLatLng('0:2'), 2) + index.add(cellIDFromLatLng('0:3'), 3) + index.add(cellIDFromLatLng('0:4'), 4) + index.build() + + const query = new ClosestCellQuery(index) + query.options.maxResults = 3 + + // Under test: Find closest cells with maxResults limit. + const target = new PointTarget(parsePoint('0:1.5')) + const results = query.findClosestCells(target) + + // Postcondition: Returns at most maxResults cells. + ok(results.length <= 3, 'Should return at most 3 results') + ok(results.length > 0, 'Should return at least 1 result') + + // Results should be sorted by distance. + for (let i = 1; i < results.length; i++) { + ok(results[i].distance >= results[i - 1].distance, 'Results should be sorted by distance') + } + }) + + test('FindClosestCellsWithMaxDistance', () => { + // Precondition: An index with cells at varying distances. + const index = new CellIndex() + index.add(cellIDFromLatLng('0:0'), 0) + index.add(cellIDFromLatLng('0:10'), 10) + index.add(cellIDFromLatLng('0:20'), 20) + index.add(cellIDFromLatLng('0:30'), 30) + index.build() + + const query = new ClosestCellQuery(index) + query.options.maxDistance = chordangle.fromAngle(15 * (Math.PI / 180)) // 15 degrees + + // Under test: Find cells within maxDistance. + const target = new PointTarget(parsePoint('0:0')) + const results = query.findClosestCells(target) + + // Postcondition: Only cells within maxDistance are returned. + for (const result of results) { + ok(result.distance < query.options.maxDistance, 'All results should be within maxDistance') + } + }) + + test('EdgeTarget', () => { + // Precondition: An index with one cell. + const index = new CellIndex() + const cellId = cellIDFromLatLng('0:0') + index.add(cellId, 0) + index.build() + + const query = new ClosestCellQuery(index) + + // Under test: Query with an edge target that passes near the cell. + const a = parsePoint('1:0') + const b = parsePoint('-1:0') + const target = new EdgeTarget(a, b) + + const result = query.findClosestCell(target) + + // Postcondition: The cell is found. + equal(result.label, 0) + }) + + test('CellUnionTarget', () => { + // Precondition: An index with cells and a CellUnion target. + const index = new CellIndex() + index.add(cellIDFromLatLng('0:0'), 0) + index.add(cellIDFromLatLng('10:10'), 1) + index.add(cellIDFromLatLng('20:20'), 2) + index.build() + + // Create a CellUnion covering an area near one of the indexed cells. + const cap = Cap.fromCenterAngle(parsePoint('0:1'), 0.01) + const coverer = new RegionCoverer({ maxCells: 4 }) + const covering = coverer.covering(cap) + + const query = new ClosestCellQuery(index) + const target = new CellUnionTarget(covering) + + const result = query.findClosestCell(target) + + // Postcondition: The closest cell is found. + // The cell at 0:0 should be closest to the covering near 0:1. + equal(result.label, 0) + }) + + describe('ShapeIndexTarget', () => { + // Precondition: A CellIndex with cells covering the Bay Area region. + const index = new CellIndex() + const s2Tokens = [ + '808fb7dc', + '808fb7d4', + '808fb7c4', + '808fb7bc', + '808fb7b4', + '808fb7cc', + '808fb7ac', + '808fb7a4', + '808fb7ec', + '808fb794' + ] + for (const token of s2Tokens) { + index.add(cellid.fromToken(token), 0) + } + index.build() + + test('ShouldReturnNoTiles', () => { + // Precondition: A polygon target that does not intersect any indexed cells. + const polygon = makeRectPolygon(-78.01263, 43.7272, -78.0126, 43.72733) + const targetIndex = new ShapeIndex() + targetIndex.add(polygon) + const target = new ShapeIndexTarget(targetIndex) + + // Under test: Query with a ShapeIndex target far from the indexed cells. + const options = new ClosestCellQueryOptions() + options.setInclusiveMaxDistance(0) + const query = new ClosestCellQuery(index, options) + const results = query.findClosestCells(target) + + // Postcondition: No cells intersect the target polygon. + equal(results.length, 0) + }) + + test('ShouldReturnSomeTiles', () => { + // Precondition: A polygon target that intersects some of the indexed cells. + const polygon = makeRectPolygon(37.4074034, -122.0187345, 37.4149719, -122.0092462) + + const targetIndex = new ShapeIndex() + targetIndex.add(polygon) + const target = new ShapeIndexTarget(targetIndex) + + // Under test: Query with a ShapeIndex target that overlaps some indexed cells. + const options = new ClosestCellQueryOptions() + options.setInclusiveMaxDistance(0) + const query = new ClosestCellQuery(index, options) + const results = query.findClosestCells(target) + + // Postcondition: Should contain the expected S2 cells. + ok(results.length === 3, "Should return 3 cells") + // And they should be the expected cells. + const expectedTokens = ['808fb7bc', '808fb7c4', '808fb7cc'] + const resultTokens = results.map(r => cellid.toToken(r.cellID)) + expectedTokens.forEach(token => + ok(resultTokens.includes(token), `results should contain cell: ${token}`) + ) + }) + + test('ShouldReturnAllTiles', () => { + // Precondition: A polygon target that covers all the indexed cells. + const polygon = makeRectPolygon(28.3911156, -133.2602933, 46.417598, -110.7676874) + const targetIndex = new ShapeIndex() + targetIndex.add(polygon) + const target = new ShapeIndexTarget(targetIndex) + + // Under test: Query with a ShapeIndex target that covers all indexed cells. + const options = new ClosestCellQueryOptions() + options.setInclusiveMaxDistance(0) + const query = new ClosestCellQuery(index, options) + const results = query.findClosestCells(target) + + // Postcondition: All cells intersect the target polygon. + equal(results.length, s2Tokens.length) + }) + }) + + describe('InclusiveMaxDistance', () => { + test('RegularMaxDistanceExcludesExactDistance', () => { + // Precondition: An index with one cell and options with regular max distance. + const id0 = cellIDFromLatLng('0:0') + const id1 = cellIDFromLatLng('1:0') + + const index = new CellIndex() + index.add(id0, 0) + index.build() + + // Compute exact distance between cells. + const cell0 = Cell.fromCellID(id0) + const cell1 = Cell.fromCellID(id1) + const exactDistance = cell0.distanceToCell(cell1) + + const options = new ClosestCellQueryOptions() + options.maxDistance = exactDistance + + const query = new ClosestCellQuery(index, options) + + // Under test: With regular maxDistance, exact distance is excluded. + const target = new CellTarget(cell1) + const results = query.findClosestCells(target) + + // Postcondition: No results because exact distance is excluded. + equal(results.length, 0, 'Exact distance should be excluded with regular maxDistance') + }) + + test('InclusiveMaxDistanceIncludesExactDistance', () => { + // Precondition: An index with one cell and options with inclusive max distance. + const id0 = cellIDFromLatLng('0:0') + const id1 = cellIDFromLatLng('1:0') + + const index = new CellIndex() + index.add(id0, 0) + index.build() + + // Compute exact distance between cells. + const cell0 = Cell.fromCellID(id0) + const cell1 = Cell.fromCellID(id1) + const exactDistance = cell0.distanceToCell(cell1) + + const options = new ClosestCellQueryOptions() + options.setInclusiveMaxDistance(exactDistance) + + const query = new ClosestCellQuery(index, options) + + // Under test: With inclusive maxDistance, exact distance is included. + const target = new CellTarget(cell1) + const results = query.findClosestCells(target) + + // Postcondition: One result because exact distance is included. + equal(results.length, 1, 'Exact distance should be included with inclusive maxDistance') + }) + }) + + test('MultipleLabelsForSameCell', () => { + // Precondition: An index with the same cell added multiple times with different labels. + const cellId = cellIDFromLatLng('0:0') + const index = new CellIndex() + index.add(cellId, 1) + index.add(cellId, 2) + index.add(cellId, 3) + index.build() + + const query = new ClosestCellQuery(index) + const target = new PointTarget(cellid.point(cellId)) + + // Under test: Query returns all instances. + const results = query.findClosestCells(target) + + // Postcondition: All three labels are returned. + equal(results.length, 3) + const labels = results.map((r) => r.label).sort() + deepEqual(labels, [1, 2, 3]) + }) + + test('RandomPointTargets', () => { + // Precondition: An index with random cells in a cap. + const numCells = 50 + const indexCap = Cap.fromCenterAngle(randomPoint(), 0.1) + const index = new CellIndex() + + for (let i = 0; i < numCells; i++) { + const point = samplePointFromCap(indexCap) + index.add(cellid.fromPoint(point), i) + } + index.build() + + const query = new ClosestCellQuery(index) + query.options.maxResults = 5 + + // Under test: Query with random point targets. + for (let i = 0; i < 10; i++) { + const targetPoint = samplePointFromCap(indexCap) + const target = new PointTarget(targetPoint) + + const results = query.findClosestCells(target) + + // Postcondition: Results are valid and sorted. + ok(results.length > 0, 'Should find at least one cell') + ok(results.length <= 5, 'Should return at most maxResults cells') + + for (let j = 1; j < results.length; j++) { + ok(results[j].distance >= results[j - 1].distance, 'Results should be sorted by distance') + } + } + }) + + test('GetDistanceReturnsClosest', () => { + // Precondition: An index with cells at known distances. + const index = new CellIndex() + index.add(cellIDFromLatLng('0:0'), 0) + index.add(cellIDFromLatLng('0:5'), 5) + index.add(cellIDFromLatLng('0:10'), 10) + index.build() + + const query = new ClosestCellQuery(index) + + // Under test: GetDistance returns the minimum distance. + const target = new PointTarget(parsePoint('0:2')) + const distance = query.getDistance(target) + + // Postcondition: Distance should be to the closest cell (0:0). + const expectedCell = Cell.fromCellID(cellIDFromLatLng('0:0')) + const expectedDistance = expectedCell.distance(parsePoint('0:2')) + + // Allow some tolerance due to cell center vs actual point. + ok(Math.abs(distance - expectedDistance) < 0.1, 'Distance should be approximately to closest cell') + }) +}) diff --git a/s2/ClosestEdgeQuery.ts b/s2/ClosestEdgeQuery.ts new file mode 100644 index 0000000..6885918 --- /dev/null +++ b/s2/ClosestEdgeQuery.ts @@ -0,0 +1,1193 @@ +/** + * S2ClosestEdgeQuery is a helper class for searching within an S2ShapeIndex + * to find the closest edge(s) to a given point, edge, S2Cell, or geometry + * collection. For example, given a set of polylines, the following code + * efficiently finds the closest 5 edges to a query point: + * + * ```typescript + * const query = new ClosestEdgeQuery(index) + * query.options.maxResults = 5 + * const target = new ClosestEdgeQuery.PointTarget(point) + * for (const result of query.findClosestEdges(target)) { + * // result.distance is the distance to the edge. + * // result.shapeID identifies the Shape containing the edge. + * // result.edgeID identifies the edge within the shape. + * // result.isInterior() indicates that the result is an interior point. + * const edge = query.getEdge(result) + * const closestPoint = query.project(point, result) + * } + * ``` + * + * You can find either the k closest edges, or all edges within a given + * radius, or both (i.e., the k closest edges up to a given maximum radius). + * By default *all* edges are returned, so you should always specify either + * maxResults or maxDistance or both. + * + * Note that by default, distances are measured to the boundary and interior + * of polygons. For example, if a point is inside a polygon then its distance + * is zero. To change this behavior, set includeInteriors to false. + * + * If you only need to test whether the distance is above or below a given + * threshold (e.g., 10 km), you can use the isDistanceLess() method. This is + * much faster than actually calculating the distance with findClosestEdge(), + * since the implementation can stop as soon as it can prove that the minimum + * distance is either above or below the threshold. + * + * @module ClosestEdgeQuery + */ + +import type { CellID } from './cellid' +import type { ChordAngle } from '../s1/chordangle' +import type { Edge } from './Shape' + +import { Cell } from './Cell' +import { Point } from './Point' +import { ShapeIndex, INDEXED, SUBDIVIDED } from './ShapeIndex' +import { ShapeIndexIterator } from './ShapeIndexIterator' +import { ContainsPointQuery, VERTEX_MODEL_SEMI_OPEN } from './ContainsPointQuery' +import { CROSS, crossingSign } from './edge_crossings' +import * as cellid from './cellid' +import * as chordangle from '../s1/chordangle' +import { STRAIGHT_CHORDANGLE } from '../s1/chordangle_constants' +import { + updateMinDistance, + minUpdateDistanceMaxError, + project as projectPointToEdge +} from './edge_distances' + +/** + * A function type for filtering shapes during queries. + * Returns true if the shape should be included in the query. + */ +export type ShapeFilter = (shapeID: number) => boolean + +/** + * A function type for visiting results during queries. + * Returns true to continue visiting, false to stop. + */ +export type ResultVisitor = (result: ClosestEdgeQueryResult) => boolean + +/** + * Represents a closest edge result from the query. + * Each result object represents a closest edge. + */ +export interface ClosestEdgeQueryResult { + /** The distance from the target to this edge. */ + distance: ChordAngle + + /** Identifies the S2Shape containing the edge. */ + shapeID: number + + /** Identifies the edge within the shape. */ + edgeID: number +} + +/** + * Returns true if this result object represents the interior of a shape. + * Such results may be returned when options.includeInteriors is true. + */ +export function isInteriorResult(result: ClosestEdgeQueryResult): boolean { + return result.shapeID >= 0 && result.edgeID < 0 +} + +/** + * Returns true if this result object indicates that no edge satisfies + * the given query options. This result is only returned in one special + * case, namely when findClosestEdge() does not find any suitable edges. + * It is never returned by methods that return a vector of results. + */ +export function isEmptyResult(result: ClosestEdgeQueryResult): boolean { + return result.shapeID < 0 +} + +/** + * Returns an empty result indicating no edge was found. + */ +export function emptyClosestEdgeQueryResult(): ClosestEdgeQueryResult { + return { + distance: chordangle.infChordAngle(), + shapeID: -1, + edgeID: -1 + } +} + +/** + * Compares two results first by distance, then by (shapeID, edgeID). + */ +function compareResults(a: ClosestEdgeQueryResult, b: ClosestEdgeQueryResult): number { + if (a.distance < b.distance) return -1 + if (a.distance > b.distance) return 1 + if (a.shapeID < b.shapeID) return -1 + if (a.shapeID > b.shapeID) return 1 + if (a.edgeID < b.edgeID) return -1 + if (a.edgeID > b.edgeID) return 1 + return 0 +} + +/** + * Options that control the set of edges returned. Note that by default + * *all* edges are returned, so you will always want to set either the + * maxResults option or the maxDistance option (or both). + */ +export class ClosestEdgeQueryOptions { + /** + * Specifies that at most "maxResults" edges should be returned. + * Default: Infinity (no limit) + */ + maxResults: number = Infinity + + /** + * Specifies that only edges whose distance to the target is less than + * "maxDistance" should be returned. + * + * Note that edges whose distance is exactly equal to "maxDistance" are + * not returned. Normally this doesn't matter, because distances are not + * computed exactly in the first place, but if such edges are needed then + * see setInclusiveMaxDistance() below. + * + * Default: Infinity + */ + maxDistance: ChordAngle = chordangle.infChordAngle() + + /** + * Specifies that edges up to maxError further away than the true + * closest edges may be substituted in the result set, as long as such + * edges satisfy all the remaining search criteria (such as maxDistance). + * This option only has an effect if maxResults is also specified; + * otherwise all edges closer than maxDistance will always be returned. + * + * Default: 0 + */ + maxError: ChordAngle = 0 + + /** + * Specifies that polygon interiors should be included when measuring + * distances. If true, the distance to a point inside a polygon is zero. + * + * Default: true + */ + includeInteriors: boolean = true + + /** + * Specifies that distances should be computed by examining every edge + * rather than using the ShapeIndex. This is useful for testing and + * debugging, and also for very small indexes where the overhead of + * building the index is not worthwhile. + * + * Default: false + */ + useBruteForce: boolean = false + + /** + * Sets maxDistance to the given value such that edges whose distance + * is exactly equal to maxDistance are also returned. + * Equivalent to setting maxDistance to maxDistance.successor(). + */ + setInclusiveMaxDistance(maxDistance: ChordAngle): void { + this.maxDistance = chordangle.successor(maxDistance) + } + + /** + * Sets maxDistance such that edges whose true distance is less than + * or equal to maxDistance will be returned (along with some edges whose + * true distance is slightly greater). + * + * This ensures that all edges whose true distance is less than or equal + * to maxDistance will be returned. The maxDistance is increased by the + * maximum error in the distance calculation. + */ + setConservativeMaxDistance(maxDistance: ChordAngle): void { + this.maxDistance = chordangle.successor( + chordangle.expanded(maxDistance, minUpdateDistanceMaxError(maxDistance)) + ) + } + + /** + * Creates a copy of these options. + */ + clone(): ClosestEdgeQueryOptions { + const copy = new ClosestEdgeQueryOptions() + copy.maxResults = this.maxResults + copy.maxDistance = this.maxDistance + copy.maxError = this.maxError + copy.includeInteriors = this.includeInteriors + copy.useBruteForce = this.useBruteForce + return copy + } +} + +/** + * Target represents the geometry to which the distance is measured. + * This is the base interface for all target types. + */ +export interface ClosestEdgeQueryTarget { + /** + * Returns the maximum number of edges in the index for which it is + * faster to use brute force search rather than the hierarchical method. + */ + maxBruteForceIndexSize(): number + + /** + * Updates minDist if the distance to the edge (v0, v1) is less than minDist. + * Returns true if the distance was updated. + */ + updateMinDistanceToEdge(v0: Point, v1: Point, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } + + /** + * Updates minDist if the distance to the point is less than minDist. + * Returns true if the distance was updated. + */ + updateMinDistanceToPoint(point: Point, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } + + /** + * Updates minDist if the distance to the cell is less than minDist. + * Returns true if the distance was updated. + */ + updateMinDistanceToCell(cell: Cell, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } + + /** + * Returns true if the target includes the interior of polygons. + */ + includeInteriors(): boolean + + /** + * Sets whether to include polygon interiors when measuring distances. + */ + setIncludeInteriors(includeInteriors: boolean): void + + /** + * Sets the max error for distance calculations. + * Returns true if the target uses this error for optimization. + */ + setMaxError(maxError: ChordAngle): boolean + + /** + * Visits shape IDs in the index that contain the target. + * Used for include_interiors support. + */ + visitContainingShapeIds( + index: ShapeIndex, + visitor: (shapeID: number, point: Point) => boolean + ): boolean + + /** + * Returns a cap that bounds the target geometry. + */ + getCapBound(): { center: Point; radius: ChordAngle } | null +} + +/** + * Target subtype that computes the closest distance to a point. + */ +export class PointTarget implements ClosestEdgeQueryTarget { + readonly point: Point + private _includeInteriors: boolean = true + + constructor(point: Point) { + this.point = point + } + + maxBruteForceIndexSize(): number { + // Break-even points are approximately 80, 100, and 250 edges for point + // cloud, fractal, and regular loop geometry respectively. + return 120 + } + + updateMinDistanceToEdge(v0: Point, v1: Point, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + const result = updateMinDistance(this.point, v0, v1, minDist) + return { distance: result.dist, updated: result.less } + } + + updateMinDistanceToPoint(point: Point, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + const dist = Point.chordAngleBetweenPoints(this.point, point) + if (dist < minDist) { + return { distance: dist, updated: true } + } + return { distance: minDist, updated: false } + } + + updateMinDistanceToCell(cell: Cell, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + const dist = cell.distance(this.point) + if (dist < minDist) { + return { distance: dist, updated: true } + } + return { distance: minDist, updated: false } + } + + includeInteriors(): boolean { + return this._includeInteriors + } + + setIncludeInteriors(includeInteriors: boolean): void { + this._includeInteriors = includeInteriors + } + + setMaxError(_maxError: ChordAngle): boolean { + return false + } + + visitContainingShapeIds( + index: ShapeIndex, + visitor: (shapeID: number, point: Point) => boolean + ): boolean { + const query = new ContainsPointQuery(index, VERTEX_MODEL_SEMI_OPEN) + return query.visitContainingShapes(this.point, (shape) => { + const shapeID = index.idForShape(shape) + if (shapeID >= 0) { + return visitor(shapeID, this.point) + } + return true + }) + } + + getCapBound(): { center: Point; radius: ChordAngle } | null { + return { center: this.point, radius: 0 } + } +} + +/** + * Target subtype that computes the closest distance to an edge. + */ +export class EdgeTarget implements ClosestEdgeQueryTarget { + readonly a: Point + readonly b: Point + private _includeInteriors: boolean = true + + constructor(a: Point, b: Point) { + this.a = a + this.b = b + } + + maxBruteForceIndexSize(): number { + // Break-even points are approximately 40, 50, and 100 edges. + return 60 + } + + updateMinDistanceToEdge(v0: Point, v1: Point, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + // Compute edge-to-edge distance. + return updateEdgePairMinDistance(this.a, this.b, v0, v1, minDist) + } + + updateMinDistanceToPoint(point: Point, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + const result = updateMinDistance(point, this.a, this.b, minDist) + return { distance: result.dist, updated: result.less } + } + + updateMinDistanceToCell(cell: Cell, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + const dist = cell.distanceToEdge(this.a, this.b) + if (dist < minDist) { + return { distance: dist, updated: true } + } + return { distance: minDist, updated: false } + } + + includeInteriors(): boolean { + return this._includeInteriors + } + + setIncludeInteriors(includeInteriors: boolean): void { + this._includeInteriors = includeInteriors + } + + setMaxError(_maxError: ChordAngle): boolean { + return false + } + + visitContainingShapeIds( + index: ShapeIndex, + visitor: (shapeID: number, point: Point) => boolean + ): boolean { + // Check if either endpoint is contained. + const query = new ContainsPointQuery(index, VERTEX_MODEL_SEMI_OPEN) + const visited = new Set() + + for (const point of [this.a, this.b]) { + const result = query.visitContainingShapes(point, (shape) => { + const shapeID = index.idForShape(shape) + if (shapeID >= 0 && !visited.has(shapeID)) { + visited.add(shapeID) + return visitor(shapeID, point) + } + return true + }) + if (!result) return false + } + return true + } + + getCapBound(): { center: Point; radius: ChordAngle } | null { + // Return cap centered at midpoint with radius to cover both endpoints. + const mid = Point.fromVector(this.a.vector.add(this.b.vector).normalize()) + const radius = Math.max( + Point.chordAngleBetweenPoints(mid, this.a), + Point.chordAngleBetweenPoints(mid, this.b) + ) + return { center: mid, radius } + } +} + +/** + * Target subtype that computes the closest distance to an S2Cell + * (including the interior of the cell). + */ +export class CellTarget implements ClosestEdgeQueryTarget { + readonly cell: Cell + private _includeInteriors: boolean = true + + constructor(cell: Cell) { + this.cell = cell + } + + maxBruteForceIndexSize(): number { + // Break-even points are approximately 20, 25, and 40 edges. + return 30 + } + + updateMinDistanceToEdge(v0: Point, v1: Point, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + const dist = this.cell.distanceToEdge(v0, v1) + if (dist < minDist) { + return { distance: dist, updated: true } + } + return { distance: minDist, updated: false } + } + + updateMinDistanceToPoint(point: Point, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + const dist = this.cell.distance(point) + if (dist < minDist) { + return { distance: dist, updated: true } + } + return { distance: minDist, updated: false } + } + + updateMinDistanceToCell(cell: Cell, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + const dist = this.cell.distanceToCell(cell) + if (dist < minDist) { + return { distance: dist, updated: true } + } + return { distance: minDist, updated: false } + } + + includeInteriors(): boolean { + return this._includeInteriors + } + + setIncludeInteriors(includeInteriors: boolean): void { + this._includeInteriors = includeInteriors + } + + setMaxError(_maxError: ChordAngle): boolean { + return false + } + + visitContainingShapeIds( + index: ShapeIndex, + visitor: (shapeID: number, point: Point) => boolean + ): boolean { + // Check if the cell center is contained. + const query = new ContainsPointQuery(index, VERTEX_MODEL_SEMI_OPEN) + const center = this.cell.center() + return query.visitContainingShapes(center, (shape) => { + const shapeID = index.idForShape(shape) + if (shapeID >= 0) { + return visitor(shapeID, center) + } + return true + }) + } + + getCapBound(): { center: Point; radius: ChordAngle } | null { + const center = this.cell.center() + // Use max distance from center to any vertex as radius. + let radius: ChordAngle = 0 + for (let i = 0; i < 4; i++) { + const d = Point.chordAngleBetweenPoints(center, this.cell.vertex(i)) + if (d > radius) radius = d + } + return { center, radius } + } +} + +/** + * Target subtype that computes the closest distance to an S2ShapeIndex + * (an arbitrary collection of points, polylines, and/or polygons). + * + * By default, distances are measured to the boundary and interior of + * polygons in the S2ShapeIndex rather than to polygon boundaries only. + * If you wish to change this behavior, you may call: + * + * target.setIncludeInteriors(false) + */ +export class ShapeIndexTarget implements ClosestEdgeQueryTarget { + readonly index: ShapeIndex + private _includeInteriors: boolean = true + private query: ContainsPointQuery + + constructor(index: ShapeIndex) { + this.index = index + this.query = new ContainsPointQuery(index, VERTEX_MODEL_SEMI_OPEN) + } + + maxBruteForceIndexSize(): number { + // Break-even points are approximately 20, 30, and 40 edges. + return 25 + } + + updateMinDistanceToEdge(v0: Point, v1: Point, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + // Check if either endpoint is contained by a polygon. + if (this._includeInteriors) { + if (this.query.contains(v0) || this.query.contains(v1)) { + return { distance: 0, updated: minDist > 0 } + } + } + + // Compute distance to all edges in the target index. + let updated = false + for (const [_shapeID, shape] of this.index.shapes) { + if (!shape) continue + const numEdges = shape.numEdges() + for (let i = 0; i < numEdges; i++) { + const edge = shape.edge(i) + const result = updateEdgePairMinDistance(v0, v1, edge.v0, edge.v1, minDist) + if (result.updated) { + minDist = result.distance + updated = true + } + if (minDist === 0) return { distance: 0, updated } + } + } + return { distance: minDist, updated } + } + + updateMinDistanceToPoint(point: Point, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + // Check if the point is contained by a polygon. + if (this._includeInteriors && this.query.contains(point)) { + return { distance: 0, updated: minDist > 0 } + } + + // Compute distance to all edges in the target index. + let updated = false + for (const [_shapeID, shape] of this.index.shapes) { + if (!shape) continue + const numEdges = shape.numEdges() + for (let i = 0; i < numEdges; i++) { + const edge = shape.edge(i) + const result = updateMinDistance(point, edge.v0, edge.v1, minDist) + if (result.less) { + minDist = result.dist + updated = true + } + if (minDist === 0) return { distance: 0, updated } + } + } + return { distance: minDist, updated } + } + + updateMinDistanceToCell(cell: Cell, minDist: ChordAngle): { distance: ChordAngle; updated: boolean } { + // Check if the cell center is contained by a polygon. + if (this._includeInteriors && this.query.contains(cell.center())) { + return { distance: 0, updated: minDist > 0 } + } + + // Compute distance to all edges in the target index. + let updated = false + for (const [_shapeID, shape] of this.index.shapes) { + if (!shape) continue + const numEdges = shape.numEdges() + for (let i = 0; i < numEdges; i++) { + const edge = shape.edge(i) + const dist = cell.distanceToEdge(edge.v0, edge.v1) + if (dist < minDist) { + minDist = dist + updated = true + } + if (minDist === 0) return { distance: 0, updated } + } + } + return { distance: minDist, updated } + } + + includeInteriors(): boolean { + return this._includeInteriors + } + + setIncludeInteriors(includeInteriors: boolean): void { + this._includeInteriors = includeInteriors + } + + setMaxError(_maxError: ChordAngle): boolean { + return false + } + + visitContainingShapeIds( + index: ShapeIndex, + visitor: (shapeID: number, point: Point) => boolean + ): boolean { + // For each shape in the query index, check if it contains any point + // from the target index. + const queryContains = new ContainsPointQuery(index, VERTEX_MODEL_SEMI_OPEN) + const visited = new Set() + + // Check if target shapes contain any query shape centers. + for (const [_shapeID, shape] of this.index.shapes) { + if (!shape) continue + const numEdges = shape.numEdges() + for (let i = 0; i < numEdges; i++) { + const edge = shape.edge(i) + for (const point of [edge.v0, edge.v1]) { + const result = queryContains.visitContainingShapes(point, (queryShape) => { + const id = index.idForShape(queryShape) + if (id >= 0 && !visited.has(id)) { + visited.add(id) + return visitor(id, point) + } + return true + }) + if (!result) return false + } + } + } + return true + } + + getCapBound(): { center: Point; radius: ChordAngle } | null { + // Return null if the index is empty. + if (this.index.len() === 0) return null + + // Compute bounding cap over all edges. + let center: Point | null = null + let radius: ChordAngle = 0 + + for (const [_shapeID, shape] of this.index.shapes) { + if (!shape) continue + const numEdges = shape.numEdges() + for (let i = 0; i < numEdges; i++) { + const edge = shape.edge(i) + if (center === null) { + center = edge.v0 + } + const d0 = Point.chordAngleBetweenPoints(center, edge.v0) + const d1 = Point.chordAngleBetweenPoints(center, edge.v1) + radius = Math.max(radius, d0, d1) + } + } + + if (center === null) return null + return { center, radius } + } +} + +/** + * Computes the minimum distance between two edges. + */ +function updateEdgePairMinDistance( + a0: Point, + a1: Point, + b0: Point, + b1: Point, + minDist: ChordAngle +): { distance: ChordAngle; updated: boolean } { + if (minDist === 0) return { distance: 0, updated: false } + if (crossingSign(a0, a1, b0, b1) === CROSS) { + return { distance: 0, updated: true } + } + + // Otherwise, the minimum distance is achieved at an endpoint of at least + // one of the two edges. + let updated = false + let result = updateMinDistance(a0, b0, b1, minDist) + if (result.less) { + minDist = result.dist + updated = true + } + result = updateMinDistance(a1, b0, b1, minDist) + if (result.less) { + minDist = result.dist + updated = true + } + result = updateMinDistance(b0, a0, a1, minDist) + if (result.less) { + minDist = result.dist + updated = true + } + result = updateMinDistance(b1, a0, a1, minDist) + if (result.less) { + minDist = result.dist + updated = true + } + return { distance: minDist, updated } +} + +/** + * Priority queue entry for the optimized algorithm. + */ +interface QueueEntry { + distance: ChordAngle + id: CellID + indexCell: { shapes: { shapeID: number; edges: number[]; containsCenter: boolean; numEdges: () => number }[] } | null +} + +/** + * S2ClosestEdgeQuery is a helper class for searching within an S2ShapeIndex + * to find the closest edge(s) to a given point, edge, S2Cell, or geometry + * collection. + */ +export class ClosestEdgeQuery { + private readonly index: ShapeIndex + readonly options: ClosestEdgeQueryOptions + private iter: ShapeIndexIterator + private indexNumEdges: number = 0 + private indexNumEdgesLimit: number = 0 + + /** + * Constructs a new ClosestEdgeQuery for the given ShapeIndex. + * Options may be specified here or changed at any time using the options property. + * + * REQUIRES: "index" must persist for the lifetime of this object. + * REQUIRES: reInit() must be called if "index" is modified. + */ + constructor(index: ShapeIndex, options: ClosestEdgeQueryOptions = new ClosestEdgeQueryOptions()) { + this.index = index + this.options = options + this.iter = index.iterator() + } + + /** + * Reinitializes the query. This method must be called whenever the + * underlying S2ShapeIndex is modified. + */ + reInit(): void { + this.indexNumEdges = 0 + this.indexNumEdgesLimit = 0 + this.iter = this.index.iterator() + } + + /** + * Returns the closest edges to the given target that satisfy the current + * options. This method may be called multiple times. + * + * Note that if options.includeInteriors is true, the result vector may + * include some entries with edgeID == -1. This indicates that the target + * intersects the indexed polygon with the given shapeID. Such results may + * be identified by calling isInteriorResult(). + */ + findClosestEdges(target: ClosestEdgeQueryTarget, filter?: ShapeFilter): ClosestEdgeQueryResult[] { + return this.findClosestEdgesInternal(target, this.options, filter) + } + + /** + * Calls a callback with the closest edges to the given target that satisfy + * the given options. Edges are reported in order of increasing distance. + * + * Updating the state that the ShapeFilter accesses while visiting is allowed + * and can be used to disable reporting of results on the fly. + */ + visitClosestEdges( + target: ClosestEdgeQueryTarget, + options: ClosestEdgeQueryOptions, + visitor: ResultVisitor, + filter?: ShapeFilter + ): void { + const results = this.findClosestEdgesInternal(target, options, filter) + for (const result of results) { + if (!visitor(result)) break + } + } + + /** + * Calls a callback with the closest edge of each shape to the given target + * that satisfies the given options. Shapes are reported in order of + * increasing distance. + */ + visitClosestShapes( + target: ClosestEdgeQueryTarget, + options: ClosestEdgeQueryOptions, + visitor: ResultVisitor, + filter?: ShapeFilter + ): void { + const results = this.findClosestEdgesInternal(target, options, filter) + const seenShapes = new Set() + let lastShape = -1 + + for (const result of results) { + const shapeID = result.shapeID + if (shapeID !== lastShape && !seenShapes.has(shapeID)) { + seenShapes.add(shapeID) + lastShape = shapeID + if (!visitor(result)) break + } + } + } + + /** + * Returns the closest edge to the target. If no edge satisfies the search + * criteria, then the result object's isEmptyResult() will return true. + * + * Note that if options.includeInteriors is true, isInteriorResult() + * should be called to check whether the result represents an interior point + * (in which case edgeID == -1). + */ + findClosestEdge(target: ClosestEdgeQueryTarget, filter?: ShapeFilter): ClosestEdgeQueryResult { + const tmpOptions = this.options.clone() + tmpOptions.maxResults = 1 + const results = this.findClosestEdgesInternal(target, tmpOptions, filter) + if (results.length === 0) { + return emptyClosestEdgeQueryResult() + } + return results[0] + } + + /** + * Returns the minimum distance to the target. If the index or target is + * empty, returns Infinity. + * + * Use isDistanceLess() if you only want to compare the distance against a + * threshold value, since it is often much faster. + */ + getDistance(target: ClosestEdgeQueryTarget, filter?: ShapeFilter): ChordAngle { + return this.findClosestEdge(target, filter).distance + } + + /** + * Returns true if the distance to "target" is less than "limit". + * + * This method is usually much faster than getDistance(), since it is much + * less work to determine whether the minimum distance is above or below a + * threshold than it is to calculate the actual minimum distance. + */ + isDistanceLess(target: ClosestEdgeQueryTarget, limit: ChordAngle, filter?: ShapeFilter): boolean { + const tmpOptions = this.options.clone() + tmpOptions.maxResults = 1 + tmpOptions.maxDistance = limit + tmpOptions.maxError = STRAIGHT_CHORDANGLE + const result = this.findClosestEdgesInternal(target, tmpOptions, filter) + return result.length > 0 + } + + /** + * Like isDistanceLess(), but also returns true if the distance to "target" + * is exactly equal to "limit". + */ + isDistanceLessOrEqual(target: ClosestEdgeQueryTarget, limit: ChordAngle, filter?: ShapeFilter): boolean { + const tmpOptions = this.options.clone() + tmpOptions.maxResults = 1 + tmpOptions.setInclusiveMaxDistance(limit) + tmpOptions.maxError = STRAIGHT_CHORDANGLE + const result = this.findClosestEdgesInternal(target, tmpOptions, filter) + return result.length > 0 + } + + /** + * Like isDistanceLessOrEqual(), except that "limit" is increased by the + * maximum error in the distance calculation. This ensures that this + * function returns true whenever the true, exact distance is less than + * or equal to "limit". + * + * For example, suppose that we want to test whether two geometries might + * intersect each other after they are snapped together using S2Builder + * (using the IdentitySnapFunction with a given "snap_radius"). Since + * S2Builder uses exact distance predicates (s2predicates.h), we need to + * measure the distance between the two geometries conservatively. If the + * distance is definitely greater than "snap_radius", then the geometries + * are guaranteed to not intersect after snapping. + */ + isConservativeDistanceLessOrEqual(target: ClosestEdgeQueryTarget, limit: ChordAngle, filter?: ShapeFilter): boolean { + const tmpOptions = this.options.clone() + tmpOptions.maxResults = 1 + tmpOptions.setConservativeMaxDistance(limit) + tmpOptions.maxError = STRAIGHT_CHORDANGLE + const result = this.findClosestEdgesInternal(target, tmpOptions, filter) + return result.length > 0 + } + + /** + * Returns the endpoints of the given result edge. + * REQUIRES: !isInteriorResult(result) + */ + getEdge(result: ClosestEdgeQueryResult): Edge { + const shape = this.index.shape(result.shapeID) + return shape.edge(result.edgeID) + } + + /** + * Returns the point on given result edge that is closest to "point". + */ + project(point: Point, result: ClosestEdgeQueryResult): Point { + if (result.edgeID < 0) return point + const edge = this.getEdge(result) + return projectPointToEdge(point, edge.v0, edge.v1) + } + + /** + * Internal method that performs the actual query. + */ + private findClosestEdgesInternal( + target: ClosestEdgeQueryTarget, + options: ClosestEdgeQueryOptions, + filter?: ShapeFilter + ): ClosestEdgeQueryResult[] { + let distanceLimit = options.maxDistance + const results: ClosestEdgeQueryResult[] = [] + + if (distanceLimit === 0) return results + + // Handle include_interiors case. + if (options.includeInteriors) { + const shapeIds = new Set() + target.visitContainingShapeIds(this.index, (shapeID: number, _point: Point) => { + if (!filter || filter(shapeID)) { + shapeIds.add(shapeID) + } + return shapeIds.size < options.maxResults + }) + + for (const shapeID of shapeIds) { + results.push({ distance: 0, shapeID, edgeID: -1 }) + } + + if (distanceLimit === 0) return results + } + + // Determine whether to use brute force or optimized algorithm. + const minOptimizedEdges = target.maxBruteForceIndexSize() + 1 + if (minOptimizedEdges > this.indexNumEdgesLimit && this.indexNumEdges >= this.indexNumEdgesLimit) { + this.indexNumEdges = this.countEdgesUpTo(minOptimizedEdges) + this.indexNumEdgesLimit = minOptimizedEdges + } + + if (options.useBruteForce || this.indexNumEdges < minOptimizedEdges) { + this.findClosestEdgesBruteForce(target, options, filter, distanceLimit, results) + } else { + this.findClosestEdgesOptimized(target, options, filter, distanceLimit, results) + } + + // Sort and limit results. + results.sort(compareResults) + + // Remove duplicates. + const uniqueResults: ClosestEdgeQueryResult[] = [] + for (const result of results) { + if ( + uniqueResults.length === 0 || + uniqueResults[uniqueResults.length - 1].shapeID !== result.shapeID || + uniqueResults[uniqueResults.length - 1].edgeID !== result.edgeID + ) { + uniqueResults.push(result) + } + } + + // Apply maxResults limit. + if (uniqueResults.length > options.maxResults) { + uniqueResults.length = options.maxResults + } + + return uniqueResults + } + + /** + * Brute force algorithm that examines every edge. + */ + private findClosestEdgesBruteForce( + target: ClosestEdgeQueryTarget, + options: ClosestEdgeQueryOptions, + filter: ShapeFilter | undefined, + distanceLimit: ChordAngle, + results: ClosestEdgeQueryResult[] + ): void { + for (const [shapeID, shape] of this.index.shapes) { + if (!shape) continue + if (filter && !filter(shapeID)) continue + + const numEdges = shape.numEdges() + for (let edgeID = 0; edgeID < numEdges; edgeID++) { + const edge = shape.edge(edgeID) + const result = target.updateMinDistanceToEdge(edge.v0, edge.v1, distanceLimit) + + if (result.updated) { + const newResult: ClosestEdgeQueryResult = { + distance: result.distance, + shapeID, + edgeID + } + + if (options.maxResults === 1) { + // Keep only the single best result. + if (results.length === 0) { + results.push(newResult) + } else { + results[0] = newResult + } + distanceLimit = result.distance - options.maxError + } else if (options.maxResults === Infinity) { + // Keep all results. + results.push(newResult) + } else { + // Keep up to maxResults results. + results.push(newResult) + if (results.length > options.maxResults) { + results.sort(compareResults) + results.length = options.maxResults + distanceLimit = results[results.length - 1].distance - options.maxError + } + } + } + } + } + } + + /** + * Optimized algorithm using priority queue and cell hierarchy. + */ + private findClosestEdgesOptimized( + target: ClosestEdgeQueryTarget, + options: ClosestEdgeQueryOptions, + filter: ShapeFilter | undefined, + distanceLimit: ChordAngle, + results: ClosestEdgeQueryResult[] + ): void { + // Get the target's bounding cap. + const capBound = target.getCapBound() + if (!capBound) return // Empty target. + + // Initialize iterator. + this.iter = this.index.iterator() + + // Optimization: if looking for just the closest edge and the cap center + // happens to intersect an index cell, process that cell first. + if (options.maxResults === 1 && this.iter.locatePoint(capBound.center)) { + const cell = this.iter.indexCell() + distanceLimit = this.processIndexCell(target, options, filter, cell, distanceLimit, results) + if (distanceLimit === 0) return + } + + // Build priority queue of cells to process. + const queue: QueueEntry[] = [] + + // Start by adding all top-level cells that intersect the target. + this.iter.begin() + while (!this.iter.done()) { + const id = this.iter.cellID() + const cell = Cell.fromCellID(id) + const result = target.updateMinDistanceToCell(cell, distanceLimit) + if (result.distance < distanceLimit) { + queue.push({ + distance: result.distance, + id, + indexCell: this.iter.indexCell() as { shapes: { shapeID: number; edges: number[]; containsCenter: boolean; numEdges: () => number }[] } + }) + } + this.iter.next() + } + + // Sort queue by distance (ascending). + queue.sort((a, b) => a.distance - b.distance) + + // Process cells in order of increasing distance. + while (queue.length > 0) { + const entry = queue.shift()! + + if (entry.distance >= distanceLimit) break + + if (entry.indexCell !== null) { + // This is an index cell; process its edges. + distanceLimit = this.processIndexCell(target, options, filter, entry.indexCell, distanceLimit, results) + } else { + // This is a parent cell; add its children to the queue. + console.log("HEYYYY") + const childIds = cellid.children(entry.id) + + for (const childId of childIds) { + const relation = this.iter.locateCellID(childId) + if (relation === INDEXED) { + const cell = Cell.fromCellID(this.iter.cellID()) + const result = target.updateMinDistanceToCell(cell, distanceLimit) + if (result.distance < distanceLimit) { + queue.push({ + distance: result.distance, + id: this.iter.cellID(), + indexCell: this.iter.indexCell() as { shapes: { shapeID: number; edges: number[]; containsCenter: boolean; numEdges: () => number }[] } + }) + } + } else if (relation === SUBDIVIDED) { + const cell = Cell.fromCellID(childId) + const result = target.updateMinDistanceToCell(cell, distanceLimit) + if (result.distance < distanceLimit) { + queue.push({ + distance: result.distance, + id: childId, + indexCell: null + }) + } + } + } + + // Re-sort queue. + queue.sort((a, b) => a.distance - b.distance) + } + } + } + + /** + * Process all edges in an index cell. + */ + private processIndexCell( + target: ClosestEdgeQueryTarget, + options: ClosestEdgeQueryOptions, + filter: ShapeFilter | undefined, + indexCell: { shapes: { shapeID: number; edges: number[]; containsCenter: boolean; numEdges: () => number }[] }, + distanceLimit: ChordAngle, + results: ClosestEdgeQueryResult[] + ): ChordAngle { + for (const clipped of indexCell.shapes) { + const shapeID = clipped.shapeID + if (filter && !filter(shapeID)) continue + + const shape = this.index.shape(shapeID) + if (!shape) continue + + for (const edgeID of clipped.edges) { + const edge = shape.edge(edgeID) + const result = target.updateMinDistanceToEdge(edge.v0, edge.v1, distanceLimit) + + if (result.updated) { + const newResult: ClosestEdgeQueryResult = { + distance: result.distance, + shapeID, + edgeID + } + + if (options.maxResults === 1) { + if (results.length === 0) { + results.push(newResult) + } else { + results[0] = newResult + } + distanceLimit = result.distance - options.maxError + } else if (options.maxResults === Infinity) { + results.push(newResult) + } else { + results.push(newResult) + if (results.length > options.maxResults) { + results.sort(compareResults) + results.length = options.maxResults + distanceLimit = results[results.length - 1].distance - options.maxError + } + } + } + } + } + return distanceLimit + } + + /** + * Counts edges in the index up to a limit. + */ + private countEdgesUpTo(limit: number): number { + let count = 0 + for (const [_shapeID, shape] of this.index.shapes) { + if (!shape) continue + count += shape.numEdges() + if (count >= limit) break + } + return count + } +} diff --git a/s2/ClosestEdgeQuery_test.ts b/s2/ClosestEdgeQuery_test.ts new file mode 100644 index 0000000..e28f28d --- /dev/null +++ b/s2/ClosestEdgeQuery_test.ts @@ -0,0 +1,965 @@ +import { test, describe } from 'node:test' +import { equal, ok } from 'node:assert/strict' +import { + ClosestEdgeQuery, + ClosestEdgeQueryOptions, + ClosestEdgeQueryResult, + PointTarget, + EdgeTarget, + CellTarget, + ShapeIndexTarget, + isInteriorResult, + isEmptyResult, + emptyClosestEdgeQueryResult +} from './ClosestEdgeQuery' +import { Cell } from './Cell' +import { LatLng } from './LatLng' +import { Loop } from './Loop' +import { Point } from './Point' +import { Polygon } from './Polygon' +import { Polyline } from './Polyline' +import { PointVector } from './PointVector' +import { ShapeIndex } from './ShapeIndex' +import * as chordangle from '../s1/chordangle' +import { parsePoint, parsePoints } from './testing_textformat' +import { Vector } from '../r3/Vector' + +/** + * Creates a rectangular S2Polygon from lat/lng bounds in degrees. + * The loop is oriented counter-clockwise (interior on the left). + */ +const makeRectPolygon = ( + bottomLeftLat: number, + bottomLeftLng: number, + topRightLat: number, + topRightLng: number +): Polygon => { + const bottomLeft = Point.fromLatLng(LatLng.fromDegrees(bottomLeftLat, bottomLeftLng)) + const bottomRight = Point.fromLatLng(LatLng.fromDegrees(bottomLeftLat, topRightLng)) + const topRight = Point.fromLatLng(LatLng.fromDegrees(topRightLat, topRightLng)) + const topLeft = Point.fromLatLng(LatLng.fromDegrees(topRightLat, bottomLeftLng)) + // CCW order: bottomLeft -> topLeft -> topRight -> bottomRight + const loop = new Loop([bottomLeft, bottomRight, topRight, topLeft]) + return new Polygon([loop]) +} + +/** + * Creates a ShapeIndex containing points from a string specification. + * Format: "lat1:lng1 | lat2:lng2 | ..." + */ +const makePointIndex = (spec: string): ShapeIndex => { + const index = new ShapeIndex() + const pointStrings = spec.split('|').map((s) => s.trim()) + for (const ps of pointStrings) { + if (ps.length === 0) continue + const point = parsePoint(ps) + // Create a degenerate polyline (single point). + const polyline = new Polyline([point, point]) + index.add(polyline) + } + return index +} + +/** + * Creates a ShapeIndex containing a polyline from a string specification. + * Format: "lat1:lng1, lat2:lng2, lat3:lng3" + */ +const makePolylineIndex = (spec: string): ShapeIndex => { + const index = new ShapeIndex() + const points = parsePoints(spec) + const polyline = new Polyline(points) + index.add(polyline) + return index +} + +/** + * Creates a ShapeIndex containing a polygon from vertices. + */ +const makePolygonIndex = (spec: string): ShapeIndex => { + const index = new ShapeIndex() + const points = parsePoints(spec) + const loop = new Loop(points) + const polygon = new Polygon([loop]) + index.add(polygon) + return index +} + +describe('S2ClosestEdgeQuery', () => { + test('NoEdges', () => { + // Precondition: An empty ShapeIndex. + const index = new ShapeIndex() + + // Under test: Query for the closest edge returns an empty result. + const query = new ClosestEdgeQuery(index) + const target = new PointTarget(new Point(1, 0, 0)) + const edge = query.findClosestEdge(target) + + // Postcondition: Result indicates no edge was found. + equal(edge.distance, chordangle.infChordAngle()) + equal(edge.shapeID, -1) + equal(edge.edgeID, -1) + ok(!isInteriorResult(edge)) + ok(isEmptyResult(edge)) + equal(query.getDistance(target), chordangle.infChordAngle()) + }) + + test('OptionsNotModified', () => { + // Precondition: An index with 3 point shapes and specific query options. + const options = new ClosestEdgeQueryOptions() + options.maxResults = 3 + options.maxDistance = chordangle.fromAngle(3 * (Math.PI / 180)) // 3 degrees + options.maxError = chordangle.fromAngle(0.001 * (Math.PI / 180)) // 0.001 degrees + + const index = makePointIndex('1:1 | 1:2 | 1:3') + const query = new ClosestEdgeQuery(index, options) + const target = new PointTarget(parsePoint('2:2')) + + // Under test: Query methods should not modify options. + const closestEdge = query.findClosestEdge(target) + // The closest point is 1:2 which is shape 1 (second shape), edge 0. + equal(closestEdge.shapeID, 1) + equal(closestEdge.edgeID, 0) + + const distance = query.getDistance(target) + const distanceDegrees = chordangle.angle(distance) * (180 / Math.PI) + ok(Math.abs(distanceDegrees - 1.0) < 0.1, `Distance should be approximately 1 degree, got ${distanceDegrees}`) + + ok(query.isDistanceLess(target, chordangle.fromAngle(1.5 * (Math.PI / 180)))) + + // Postcondition: Options remain unchanged. + equal(query.options.maxResults, 3) + equal(query.options.maxDistance, options.maxDistance) + equal(query.options.maxError, options.maxError) + }) + + test('DistanceEqualToLimit', () => { + // Precondition: An index with one point. + const p0 = parsePoint('23:12') + const p1 = parsePoint('47:11') + + const index = new ShapeIndex() + const polyline = new Polyline([p0, p0]) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + + // Under test: Distance comparison methods with zero distance. + const target0 = new PointTarget(p0) + const dist0 = 0 // Zero distance + + ok(!query.isDistanceLess(target0, dist0), 'isDistanceLess should return false for exact distance') + ok(query.isDistanceLessOrEqual(target0, dist0), 'isDistanceLessOrEqual should return true for exact distance') + ok( + query.isConservativeDistanceLessOrEqual(target0, dist0), + 'isConservativeDistanceLessOrEqual should return true' + ) + + // Under test: Distance comparison with non-zero distance. + const target1 = new PointTarget(p1) + const dist1 = Point.chordAngleBetweenPoints(p0, p1) + + ok(!query.isDistanceLess(target1, dist1), 'isDistanceLess should return false for exact distance') + ok(query.isDistanceLessOrEqual(target1, dist1), 'isDistanceLessOrEqual should return true for exact distance') + ok( + query.isConservativeDistanceLessOrEqual(target1, dist1), + 'isConservativeDistanceLessOrEqual should return true' + ) + }) + + test('TargetPointInsideIndexedPolygon', () => { + // Precondition: An index with a polyline loop and a polygon. + // The polyline loop (no interior): 0:0, 0:5, 5:5, 5:0 + // The polygon: 0:10, 0:15, 5:15, 5:10 + const polylineIndex = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:5, 5:5, 5:0')) + polylineIndex.add(polyline) + + const polygonPoints = parsePoints('0:10, 0:15, 5:15, 5:10') + const loop = new Loop(polygonPoints) + const polygon = new Polygon([loop]) + polylineIndex.add(polygon) + + const options = new ClosestEdgeQueryOptions() + options.includeInteriors = true + options.maxDistance = chordangle.fromAngle(1 * (Math.PI / 180)) // 1 degree + + const query = new ClosestEdgeQuery(polylineIndex, options) + + // Under test: Query for a point inside the indexed polygon. + const target = new PointTarget(parsePoint('2:12')) + const results = query.findClosestEdges(target) + + // Postcondition: Should find the polygon interior. + equal(results.length, 1) + equal(results[0].distance, 0) + equal(results[0].shapeID, 1) + equal(results[0].edgeID, -1) + ok(isInteriorResult(results[0])) + ok(!isEmptyResult(results[0])) + }) + + test('TargetPointOutsideIndexedPolygon', () => { + // Precondition: An index with a polyline loop (no interior) and a polygon. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:5, 5:5, 5:0')) + index.add(polyline) + + const polygonPoints = parsePoints('0:10, 0:15, 5:15, 5:10') + const loop = new Loop(polygonPoints) + const polygon = new Polygon([loop]) + index.add(polygon) + + const options = new ClosestEdgeQueryOptions() + options.includeInteriors = true + options.maxDistance = chordangle.fromAngle(1 * (Math.PI / 180)) // 1 degree + + const query = new ClosestEdgeQuery(index, options) + + // Under test: Query for a point inside the polyline loop (but not inside the polygon). + const target = new PointTarget(parsePoint('2:2')) + const results = query.findClosestEdges(target) + + // Postcondition: No results since polyline has no interior. + equal(results.length, 0) + }) + + test('TargetPolygonContainingIndexedPoints', () => { + // Precondition: An index with 4 points (each as a separate degenerate polyline shape). + // Two points are inside a polyline loop (no interior): 2:2, 3:3 + // Two points are inside a polygon: 1:11, 3:13 + const index = makePointIndex('2:2 | 3:3 | 1:11 | 3:13') + + const query = new ClosestEdgeQuery(index) + query.options.maxDistance = chordangle.fromAngle(1 * (Math.PI / 180)) // 1 degree + + // Create target with polyline and polygon. + const targetIndex = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:5, 5:5, 5:0')) + targetIndex.add(polyline) + + const polygonPoints = parsePoints('0:10, 0:15, 5:15, 5:10') + const loop = new Loop(polygonPoints) + const polygon = new Polygon([loop]) + targetIndex.add(polygon) + + const target = new ShapeIndexTarget(targetIndex) + target.setIncludeInteriors(true) + + // Under test: Query returns points inside the polygon. + const results = query.findClosestEdges(target) + + // Postcondition: Only the 2 points inside the polygon are returned. + // Each point is a separate shape (shapeID 2 = 1:11, shapeID 3 = 3:13). + equal(results.length, 2) + equal(results[0].distance, 0) + equal(results[0].shapeID, 2) // Shape containing point 1:11 + equal(results[0].edgeID, 0) + ok(!isInteriorResult(results[0])) + + equal(results[1].distance, 0) + equal(results[1].shapeID, 3) // Shape containing point 3:13 + equal(results[1].edgeID, 0) + ok(!isInteriorResult(results[1])) + }) + + test('EmptyTargetOptimized', () => { + // Precondition: An index with a large polygon. + const index = new ShapeIndex() + const loop = Loop.regularLoop(new Point(1, 0, 0), 0.1, 1000) + const polygon = new Polygon([loop]) + index.add(polygon) + + const query = new ClosestEdgeQuery(index) + query.options.maxDistance = 1e-5 // Very small radius + + // Under test: Query with an empty ShapeIndex target. + const targetIndex = new ShapeIndex() + const target = new ShapeIndexTarget(targetIndex) + + const results = query.findClosestEdges(target) + + // Postcondition: No results are returned. + equal(results.length, 0) + }) + + test('EmptyResult', () => { + // Precondition: Create an empty result. + const result = emptyClosestEdgeQueryResult() + + // Under test: Check that the result is empty. + ok(isEmptyResult(result)) + equal(result.shapeID, -1) + equal(result.edgeID, -1) + equal(result.distance, chordangle.infChordAngle()) + }) + + test('GetEdge', () => { + // Precondition: An index with a polyline. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 1:1, 2:2')) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + const target = new PointTarget(parsePoint('0.5:0.5')) + + // Under test: Find closest edge and retrieve its endpoints. + const result = query.findClosestEdge(target) + const edge = query.getEdge(result) + + // Postcondition: Edge endpoints are valid. + ok(edge.v0 instanceof Point) + ok(edge.v1 instanceof Point) + }) + + describe('ClosestEdgeQueryOptions', () => { + test('DefaultOptions', () => { + // Precondition: Create default options. + const options = new ClosestEdgeQueryOptions() + + // Under test: Check default values. + equal(options.maxResults, Infinity) + equal(options.maxDistance, chordangle.infChordAngle()) + equal(options.maxError, 0) + equal(options.includeInteriors, true) + equal(options.useBruteForce, false) + }) + + test('SetInclusiveMaxDistance', () => { + // Precondition: Create options and set inclusive max distance. + const options = new ClosestEdgeQueryOptions() + const distance = chordangle.fromAngle(1 * (Math.PI / 180)) // 1 degree + + // Under test: setInclusiveMaxDistance sets the successor. + options.setInclusiveMaxDistance(distance) + + // Postcondition: maxDistance is the successor of the input. + equal(options.maxDistance, chordangle.successor(distance)) + }) + + test('SetConservativeMaxDistance', () => { + // Precondition: Create options and set conservative max distance. + const options = new ClosestEdgeQueryOptions() + const distance = chordangle.fromAngle(1 * (Math.PI / 180)) // 1 degree + + // Under test: setConservativeMaxDistance sets distance with error. + options.setConservativeMaxDistance(distance) + + // Postcondition: maxDistance is greater than the original distance. + ok(options.maxDistance > distance) + }) + }) + + describe('PointTarget', () => { + test('Construction', () => { + // Precondition: Create a point. + const point = parsePoint('45:90') + + // Under test: Create a PointTarget. + const target = new PointTarget(point) + + // Postcondition: Target contains the point. + ok(target.point.equals(point)) + }) + + test('IncludeInteriors', () => { + // Precondition: Create a PointTarget. + const target = new PointTarget(new Point(1, 0, 0)) + + // Under test: Default includeInteriors is true. + ok(target.includeInteriors()) + + // Under test: Can set includeInteriors to false. + target.setIncludeInteriors(false) + ok(!target.includeInteriors()) + }) + }) + + describe('EdgeTarget', () => { + test('Construction', () => { + // Precondition: Create two points. + const a = parsePoint('0:0') + const b = parsePoint('1:1') + + // Under test: Create an EdgeTarget. + const target = new EdgeTarget(a, b) + + // Postcondition: Target contains the edge endpoints. + ok(target.a.equals(a)) + ok(target.b.equals(b)) + }) + + test('IncludeInteriors', () => { + // Precondition: Create an EdgeTarget. + const target = new EdgeTarget(new Point(1, 0, 0), new Point(0, 1, 0)) + + // Under test: Default includeInteriors is true. + ok(target.includeInteriors()) + + // Under test: Can set includeInteriors to false. + target.setIncludeInteriors(false) + ok(!target.includeInteriors()) + }) + }) + + describe('CellTarget', () => { + test('Construction', () => { + // Precondition: Create a cell. + const point = parsePoint('45:90') + const cell = Cell.fromPoint(point) + + // Under test: Create a CellTarget. + const target = new CellTarget(cell) + + // Postcondition: Target contains the cell. + equal(target.cell.id, cell.id) + }) + + test('IncludeInteriors', () => { + // Precondition: Create a CellTarget. + const cell = Cell.fromPoint(new Point(1, 0, 0)) + const target = new CellTarget(cell) + + // Under test: Default includeInteriors is true. + ok(target.includeInteriors()) + + // Under test: Can set includeInteriors to false. + target.setIncludeInteriors(false) + ok(!target.includeInteriors()) + }) + }) + + describe('ShapeIndexTarget', () => { + test('Construction', () => { + // Precondition: Create a ShapeIndex. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 1:1')) + index.add(polyline) + + // Under test: Create a ShapeIndexTarget. + const target = new ShapeIndexTarget(index) + + // Postcondition: Target contains the index. + equal(target.index, index) + }) + + test('IncludeInteriors', () => { + // Precondition: Create a ShapeIndexTarget. + const index = new ShapeIndex() + const target = new ShapeIndexTarget(index) + + // Under test: Default includeInteriors is true. + ok(target.includeInteriors()) + + // Under test: Can set includeInteriors to false. + target.setIncludeInteriors(false) + ok(!target.includeInteriors()) + }) + }) + + describe('FindClosestEdges', () => { + test('WithMaxResults', () => { + // Precondition: An index with multiple edges. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:1, 0:2, 0:3, 0:4, 0:5')) + index.add(polyline) + + const options = new ClosestEdgeQueryOptions() + options.maxResults = 3 + + const query = new ClosestEdgeQuery(index, options) + + // Under test: Find closest edges with maxResults limit. + const target = new PointTarget(parsePoint('0:2.5')) + const results = query.findClosestEdges(target) + + // Postcondition: Returns at most maxResults edges. + ok(results.length <= 3, 'Should return at most 3 results') + ok(results.length > 0, 'Should return at least 1 result') + + // Results should be sorted by distance. + for (let i = 1; i < results.length; i++) { + ok(results[i].distance >= results[i - 1].distance, 'Results should be sorted by distance') + } + }) + + test('WithMaxDistance', () => { + // Precondition: An index with edges at varying distances. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:10, 0:20, 0:30')) + index.add(polyline) + + const options = new ClosestEdgeQueryOptions() + options.maxDistance = chordangle.fromAngle(15 * (Math.PI / 180)) // 15 degrees + + const query = new ClosestEdgeQuery(index, options) + + // Under test: Find edges within maxDistance. + const target = new PointTarget(parsePoint('0:0')) + const results = query.findClosestEdges(target) + + // Postcondition: Only edges within maxDistance are returned. + for (const result of results) { + ok(result.distance < query.options.maxDistance, 'All results should be within maxDistance') + } + }) + }) + + test('ReuseOfQuery', () => { + // Precondition: An index with one point. + const index = new ShapeIndex() + const polyline = new Polyline([parsePoint('2:2'), parsePoint('2:2')]) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + query.options.maxError = chordangle.fromAngle(1 * (Math.PI / 180)) // 1 degree + + // Create target. + const targetIndex = makePolygonIndex('0:0, 0:5, 5:5, 5:0') + const target = new ShapeIndexTarget(targetIndex) + + // Under test: Query can be reused between calls. + const results1 = query.findClosestEdges(target) + const results2 = query.findClosestEdges(target) + + // Postcondition: Both queries return the same results. + equal(results1.length, results2.length) + }) + + test('IsDistanceLess', () => { + // Precondition: An index with one edge. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 1:1')) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + const target = new PointTarget(parsePoint('0:0')) + + // Under test: isDistanceLess returns correct values. + ok(query.isDistanceLess(target, chordangle.fromAngle(1 * (Math.PI / 180)))) + ok(!query.isDistanceLess(target, 0)) + }) + + test('IsDistanceLessOrEqual', () => { + // Precondition: An index with one point. + const p0 = parsePoint('0:0') + const index = new ShapeIndex() + const polyline = new Polyline([p0, p0]) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + const target = new PointTarget(p0) + + // Under test: isDistanceLessOrEqual returns true for zero distance. + ok(query.isDistanceLessOrEqual(target, 0)) + }) + + test('IsConservativeDistanceLessOrEqual', () => { + // Precondition: An index with one point. + const p0 = parsePoint('0:0') + const index = new ShapeIndex() + const polyline = new Polyline([p0, p0]) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + const target = new PointTarget(p0) + + // Under test: isConservativeDistanceLessOrEqual returns true for zero distance. + ok(query.isConservativeDistanceLessOrEqual(target, 0)) + }) + + test('BasicTestWithPointVector', () => { + // Precondition: An index with 3 points in a single PointVector shape. + const index = new ShapeIndex() + const points = [parsePoint('1:1'), parsePoint('1:2'), parsePoint('1:3')] + const pv = new PointVector(points) + index.add(pv) + + const options = new ClosestEdgeQueryOptions() + options.maxResults = 1 + options.maxDistance = chordangle.fromAngle(3 * (Math.PI / 180)) // 3 degrees + options.maxError = chordangle.fromAngle(0.001 * (Math.PI / 180)) // 0.001 degrees + + const query = new ClosestEdgeQuery(index, options) + const target = new PointTarget(parsePoint('2:2')) + + // Under test: Query returns the closest edge (point 1:2 is edge 1). + const result = query.findClosestEdge(target) + + // Postcondition: Edge 1 (point 1:2) is the closest. + equal(result.edgeID, 1) + + const distance = query.getDistance(target) + const distanceDegrees = chordangle.angle(distance) * (180 / Math.PI) + ok(Math.abs(distanceDegrees - 1.0) < 0.1, `Distance should be approximately 1 degree, got ${distanceDegrees}`) + + ok(query.isDistanceLess(target, chordangle.fromAngle(1.5 * (Math.PI / 180)))) + }) + + test('TrueDistanceLessThanChordAngleDistance', () => { + // Precondition: Two points with worst-case ChordAngle error. + // These points were chosen because ChordAngle distance is ~4 ulps greater than true distance. + const p0 = Point.fromVector( + new Vector(0.78516762584829192, -0.5020040069084597, -0.36263449417782678) + ) + const p1 = Point.fromVector( + new Vector(0.78563011732429433, -0.50187655940493503, -0.36180828883938054) + ) + + const index = new ShapeIndex() + const pv = new PointVector([p0]) + index.add(pv) + + const query = new ClosestEdgeQuery(index) + + // Under test: The ChordAngle distance has error compared to true distance. + const dist = Point.chordAngleBetweenPoints(p0, p1) + // Go 4 ulps back from the computed distance. + const limit = chordangle.predecessor( + chordangle.predecessor(chordangle.predecessor(chordangle.predecessor(dist))) + ) + + const target = new PointTarget(p1) + + // Postcondition: isDistanceLess should return false for the predecessor limit. + ok(!query.isDistanceLess(target, limit), 'isDistanceLess should return false') + + // Postcondition: isConservativeDistanceLessOrEqual should still return true. + ok( + query.isConservativeDistanceLessOrEqual(target, limit), + 'isConservativeDistanceLessOrEqual should return true' + ) + }) + + test('DistanceEqualToLimitWithSuccessor', () => { + // Precondition: An index with one point. + const p0 = parsePoint('23:12') + const p1 = parsePoint('47:11') + + const index = new ShapeIndex() + const pv = new PointVector([p0]) + index.add(pv) + + const query = new ClosestEdgeQuery(index) + + // Under test: isDistanceLess with successor returns true. + const target0 = new PointTarget(p0) + const dist0 = 0 + + ok(!query.isDistanceLess(target0, dist0), 'isDistanceLess should return false for exact distance') + ok( + query.isDistanceLess(target0, chordangle.successor(dist0)), + 'isDistanceLess should return true for successor' + ) + ok(query.isConservativeDistanceLessOrEqual(target0, dist0), 'isConservativeDistanceLessOrEqual should return true') + + // Under test: With non-zero distance. + const target1 = new PointTarget(p1) + const dist1 = Point.chordAngleBetweenPoints(p0, p1) + + ok(!query.isDistanceLess(target1, dist1), 'isDistanceLess should return false for exact distance') + ok( + query.isDistanceLess(target1, chordangle.successor(dist1)), + 'isDistanceLess should return true for successor' + ) + ok(query.isConservativeDistanceLessOrEqual(target1, dist1), 'isConservativeDistanceLessOrEqual should return true') + }) + + test('TargetPolygonContainingIndexedPointsWithPointVector', () => { + // Precondition: An index with 4 points in a single PointVector. + // Two points are inside a polyline loop (no interior): 2:2, 3:3 + // Two points are inside a polygon: 1:11, 3:13 + const index = new ShapeIndex() + const points = [parsePoint('2:2'), parsePoint('3:3'), parsePoint('1:11'), parsePoint('3:13')] + const pv = new PointVector(points) + index.add(pv) + + const options = new ClosestEdgeQueryOptions() + options.useBruteForce = false + options.maxDistance = chordangle.fromAngle(1 * (Math.PI / 180)) // 1 degree + + const query = new ClosestEdgeQuery(index, options) + + // Create target with polyline and polygon. + const targetIndex = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:5, 5:5, 5:0')) + targetIndex.add(polyline) + + const polygonPoints = parsePoints('0:10, 0:15, 5:15, 5:10') + const loop = new Loop(polygonPoints) + const polygon = new Polygon([loop]) + targetIndex.add(polygon) + + const target = new ShapeIndexTarget(targetIndex) + target.setIncludeInteriors(true) + + // Under test: Query returns points inside the polygon. + const results = query.findClosestEdges(target) + + // Postcondition: Only the 2 points inside the polygon are returned (edgeIDs 2 and 3). + equal(results.length, 2) + equal(results[0].distance, 0) + equal(results[0].shapeID, 0) + equal(results[0].edgeID, 2) // Point 1:11 + ok(!isInteriorResult(results[0])) + + equal(results[1].distance, 0) + equal(results[1].shapeID, 0) + equal(results[1].edgeID, 3) // Point 3:13 + ok(!isInteriorResult(results[1])) + }) + + test('BruteForceVsOptimizedConsistency', () => { + // Precondition: An index with multiple edges. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 1:1, 2:2, 3:3, 4:4')) + index.add(polyline) + + const target = new PointTarget(parsePoint('1.5:1.5')) + + // Under test: Both brute force and optimized should give same results. + const optsBruteForce = new ClosestEdgeQueryOptions() + optsBruteForce.useBruteForce = true + optsBruteForce.maxResults = 2 + + const optsOptimized = new ClosestEdgeQueryOptions() + optsOptimized.useBruteForce = false + optsOptimized.maxResults = 2 + + const queryBruteForce = new ClosestEdgeQuery(index, optsBruteForce) + const queryOptimized = new ClosestEdgeQuery(index, optsOptimized) + + const resultsBruteForce = queryBruteForce.findClosestEdges(target) + const resultsOptimized = queryOptimized.findClosestEdges(target) + + // Postcondition: Both return same number of results. + equal(resultsBruteForce.length, resultsOptimized.length) + + // Postcondition: Both return same closest edge. + if (resultsBruteForce.length > 0 && resultsOptimized.length > 0) { + equal(resultsBruteForce[0].shapeID, resultsOptimized[0].shapeID) + equal(resultsBruteForce[0].edgeID, resultsOptimized[0].edgeID) + } + }) + + test('EdgeTargetQuery', () => { + // Precondition: An index with a polyline. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:10, 10:10, 10:0')) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + + // Under test: Query with an edge target. + const target = new EdgeTarget(parsePoint('5:5'), parsePoint('5:15')) + const result = query.findClosestEdge(target) + + // Postcondition: Found a closest edge. + ok(!isEmptyResult(result)) + ok(result.distance >= 0) + }) + + test('CellTargetQuery', () => { + // Precondition: An index with a polyline. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:10, 10:10, 10:0')) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + + // Under test: Query with a cell target. + const cell = Cell.fromPoint(parsePoint('5:5')) + const target = new CellTarget(cell) + const result = query.findClosestEdge(target) + + // Postcondition: Found a closest edge. + ok(!isEmptyResult(result)) + ok(result.distance >= 0) + }) + + test('ShapeIndexTargetQuery', () => { + // Precondition: An index with a polyline. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:10, 10:10, 10:0')) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + + // Under test: Query with a ShapeIndex target containing another polyline. + const targetIndex = new ShapeIndex() + const targetPolyline = new Polyline(parsePoints('5:5, 5:15, 15:15')) + targetIndex.add(targetPolyline) + + const target = new ShapeIndexTarget(targetIndex) + const result = query.findClosestEdge(target) + + // Postcondition: Found a closest edge. + ok(!isEmptyResult(result)) + ok(result.distance >= 0) + }) + + test('ProjectPointToEdge', () => { + // Precondition: An index with a vertical polyline (lng = 0, lat from 0 to 10). + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 10:0')) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + const queryPoint = parsePoint('5:5') + const target = new PointTarget(queryPoint) + + // Under test: Project finds closest point on edge. + const result = query.findClosestEdge(target) + const projected = query.project(queryPoint, result) + + // Postcondition: Projected point is on the edge. + const projectedLatLng = LatLng.fromPoint(projected) + const DEGREE = Math.PI / 180 + // The edge goes from lat 0 to lat 10 degrees at lng 0. + // The projected point should have lng near 0. + const lngDegrees = projectedLatLng.lng / DEGREE + ok(Math.abs(lngDegrees) < 1, `Projected point should have lng near 0, got ${lngDegrees}`) + // The projected point should have lat between 0 and 10 degrees. + const latDegrees = projectedLatLng.lat / DEGREE + ok(latDegrees >= 0 && latDegrees <= 10, `Projected point lat should be on edge, got ${latDegrees}`) + }) + + test('ShapeFilterExcludesShapes', () => { + // Precondition: An index with multiple shapes. + const index = new ShapeIndex() + const polyline1 = new Polyline(parsePoints('0:0, 0:10')) + const polyline2 = new Polyline(parsePoints('1:0, 1:10')) + const polyline3 = new Polyline(parsePoints('2:0, 2:10')) + index.add(polyline1) // shapeID 0 + index.add(polyline2) // shapeID 1 + index.add(polyline3) // shapeID 2 + + const query = new ClosestEdgeQuery(index) + const target = new PointTarget(parsePoint('0.5:5')) + + // Under test: Filter to only include shapeID 1. + const filter = (shapeID: number) => shapeID === 1 + const result = query.findClosestEdge(target, filter) + + // Postcondition: Only shapeID 1 is returned. + equal(result.shapeID, 1) + }) + + test('ShapeFilterExcludesAllShapes', () => { + // Precondition: An index with shapes. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:10')) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + const target = new PointTarget(parsePoint('0:5')) + + // Under test: Filter that excludes all shapes. + const filter = (_shapeID: number) => false + const result = query.findClosestEdge(target, filter) + + // Postcondition: No result found. + ok(isEmptyResult(result)) + }) + + test('MultipleCallsReturnConsistentResults', () => { + // Precondition: An index with edges. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 1:1, 2:2')) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + const target = new PointTarget(parsePoint('0.5:0.5')) + + // Under test: Multiple calls return same results. + const result1 = query.findClosestEdge(target) + const result2 = query.findClosestEdge(target) + const result3 = query.findClosestEdge(target) + + // Postcondition: All results are identical. + equal(result1.shapeID, result2.shapeID) + equal(result1.edgeID, result2.edgeID) + equal(result1.distance, result2.distance) + equal(result2.shapeID, result3.shapeID) + equal(result2.edgeID, result3.edgeID) + equal(result2.distance, result3.distance) + }) + + test('VisitClosestEdges', () => { + // Precondition: An index with multiple edges. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:1, 0:2, 0:3')) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + const options = new ClosestEdgeQueryOptions() + options.maxResults = 10 + + const target = new PointTarget(parsePoint('0:1.5')) + + // Under test: Visit all closest edges. + const visited: ClosestEdgeQueryResult[] = [] + query.visitClosestEdges(target, options, (result) => { + visited.push({ ...result }) + return true + }) + + // Postcondition: Visited some edges. + ok(visited.length > 0) + + // Postcondition: Results are sorted by distance. + for (let i = 1; i < visited.length; i++) { + ok(visited[i].distance >= visited[i - 1].distance) + } + }) + + test('VisitClosestEdgesCanStopEarly', () => { + // Precondition: An index with multiple edges. + const index = new ShapeIndex() + const polyline = new Polyline(parsePoints('0:0, 0:1, 0:2, 0:3, 0:4, 0:5')) + index.add(polyline) + + const query = new ClosestEdgeQuery(index) + const options = new ClosestEdgeQueryOptions() + options.maxResults = 10 + + const target = new PointTarget(parsePoint('0:2.5')) + + // Under test: Stop after visiting 2 edges. + const visited: ClosestEdgeQueryResult[] = [] + query.visitClosestEdges(target, options, (result) => { + visited.push({ ...result }) + return visited.length < 2 + }) + + // Postcondition: Stopped after 2 edges. + equal(visited.length, 2) + }) + + test('VisitClosestShapes', () => { + // Precondition: An index with multiple shapes. + const index = new ShapeIndex() + const polyline1 = new Polyline(parsePoints('0:0, 0:5')) + const polyline2 = new Polyline(parsePoints('1:0, 1:5')) + index.add(polyline1) + index.add(polyline2) + + const query = new ClosestEdgeQuery(index) + const options = new ClosestEdgeQueryOptions() + options.maxResults = 10 + + const target = new PointTarget(parsePoint('0.5:2.5')) + + // Under test: Visit closest edge per shape. + const visited: ClosestEdgeQueryResult[] = [] + query.visitClosestShapes(target, options, (result) => { + visited.push({ ...result }) + return true + }) + + // Postcondition: Visited 2 shapes (one result per shape). + equal(visited.length, 2) + // Postcondition: Different shape IDs. + ok(visited[0].shapeID !== visited[1].shapeID) + }) +}) + diff --git a/s2/_index.ts b/s2/_index.ts index ef2ae11..91ee804 100644 --- a/s2/_index.ts +++ b/s2/_index.ts @@ -13,7 +13,29 @@ export * as cellid from './cellid' export { Cap } from './Cap' export { Cell } from './Cell' +export { CellIndex, CellIndexRangeIterator, CellIndexContentsIterator } from './CellIndex' export { CellUnion } from './CellUnion' +export { + ClosestCellQuery, + ClosestCellQueryOptions, + ClosestCellQueryResult, + ClosestCellQueryTarget, + PointTarget, + EdgeTarget, + CellTarget, + CellUnionTarget, + ShapeIndexTarget +} from './ClosestCellQuery' +export { + ClosestEdgeQuery, + ClosestEdgeQueryOptions, + ClosestEdgeQueryResult, + ClosestEdgeQueryTarget, + PointTarget as EdgeQueryPointTarget, + EdgeTarget as EdgeQueryEdgeTarget, + CellTarget as EdgeQueryCellTarget, + ShapeIndexTarget as EdgeQueryShapeIndexTarget, +} from './ClosestEdgeQuery' // export { ContainsPointQuery } from './ContainsPointQuery' // export { ContainsVertexQuery } from './ContainsVertexQuery' // export { CrossingEdgeQuery } from './CrossingEdgeQuery'