feat(s2): PaddedCell

r2/Rect.ts

@@ -234,6 +234,14 @@ export class Rect {
     )
   }
 
+  /**
+   * Constructs a rectangle from another rectangle (ie. make a copy)
+   * @category Constructors
+   */
+  static fromRect(or: Rect): Rect {
+    return new Rect(new Interval(or.x.lo, or.x.hi), new Interval(or.y.lo, or.y.hi))
+  }
+
   /**
    * Constructs the canonical empty rectangle.
    * Use isEmpty() to test for empty rectangles, since they have more than one representation.

s2/PaddedCell.ts

@@ -0,0 +1,199 @@
+import type { CellID } from './cellid'
+import { Rect as R2Rect } from '../r2/Rect'
+import { Interval as R1Interval } from '../r1/Interval'
+import * as cellid from './cellid'
+import { ijToPos, INVERT_MASK, posToIJ, posToOrientation, SWAP_MASK } from './lookupIJ'
+import { Point } from './Point'
+import { faceSiTiToXYZ, siTiToST, stToUV, uvToST } from './stuv'
+import { DBL_EPSILON } from './predicates'
+import { findMSBSetNonZero64 } from '../r1/math'
+import { MAX_LEVEL } from './cellid_constants'
+
+export class PaddedCell {
+  id: CellID = 0n
+  padding: number = 0
+  bound: R2Rect = R2Rect.empty()
+  middle: R2Rect = R2Rect.empty() // A rect in (u, v)-space that belongs to all four children.
+  iLo: number = 0 // Minimum (i,j)-coordinates of this cell before padding
+  jLo: number = 0
+  orientation: number = 0 // Hilbert curve orientation of this cell.
+  level: number = 0
+
+  // Constructs a padded cell with the given padding.
+  static fromCellID(id: CellID, padding: number): PaddedCell {
+    const p = new PaddedCell()
+    p.id = id
+    p.padding = padding
+    p.middle = R2Rect.empty()
+
+    // Fast path for constructing a top-level face (the most common case).
+    if (cellid.isFace(id)) {
+      const limit = padding + 1
+      p.bound = new R2Rect(new R1Interval(-limit, limit), new R1Interval(-limit, limit))
+      p.middle = new R2Rect(new R1Interval(-padding, padding), new R1Interval(-padding, padding))
+      p.orientation = cellid.face(id) & 1
+      return p
+    }
+
+    const { i, j, orientation } = cellid.faceIJOrientation(id)
+    p.iLo = i
+    p.jLo = j
+    p.orientation = orientation
+    p.level = cellid.level(id)
+    p.bound = cellid.ijLevelToBoundUV(i, j, p.level).expandedByMargin(padding)
+    const ijSize = cellid.sizeIJ(p.level)
+    p.iLo &= -ijSize
+    p.jLo &= -ijSize
+
+    return p
+  }
+
+  // Constructs the child of parent with the given (i,j) index.
+  static fromParentIJ(parent: PaddedCell, i: number, j: number): PaddedCell {
+    const pos = ijToPos[parent.orientation][2 * i + j]
+
+    const p = new PaddedCell()
+    p.id = cellid.children(parent.id)[pos]
+    p.padding = parent.padding
+    p.bound = R2Rect.fromRect(parent.bound)
+    p.orientation = parent.orientation ^ posToOrientation[pos]
+    p.level = parent.level + 1
+    p.middle = R2Rect.empty()
+
+    const ijSize = cellid.sizeIJ(p.level)
+    p.iLo = parent.iLo + i * ijSize
+    p.jLo = parent.jLo + j * ijSize
+
+    const middle = parent.middleRect()
+    if (i === 1) {
+      p.bound.x.lo = middle.x.lo
+    } else {
+      p.bound.x.hi = middle.x.hi
+    }
+    if (j === 1) {
+      p.bound.y.lo = middle.y.lo
+    } else {
+      p.bound.y.hi = middle.y.hi
+    }
+
+    return p
+  }
+
+  // Returns the CellID this padded cell represents.
+  cellID(): CellID {
+    return this.id
+  }
+
+  // Returns the amount of padding on this cell.
+  paddingValue(): number {
+    return this.padding
+  }
+
+  // Returns the level this cell is at.
+  levelValue(): number {
+    return this.level
+  }
+
+  // Returns the center of this cell.
+  center(): Point {
+    const ijSize = cellid.sizeIJ(this.level)
+    const si = 2 * this.iLo + ijSize
+    const ti = 2 * this.jLo + ijSize
+    return Point.fromVector(faceSiTiToXYZ(cellid.face(this.id), si, ti).vector.normalize())
+  }
+
+  // Returns the rectangle in the middle of this cell that belongs to all four of its children in (u,v)-space.
+  middleRect(): R2Rect {
+    if (this.middle.isEmpty()) {
+      const ijSize = cellid.sizeIJ(this.level)
+      const u = stToUV(siTiToST(2 * this.iLo + ijSize))
+      const v = stToUV(siTiToST(2 * this.jLo + ijSize))
+      this.middle = new R2Rect(
+        new R1Interval(u - this.padding, u + this.padding),
+        new R1Interval(v - this.padding, v + this.padding)
+      )
+    }
+    return this.middle
+  }
+
+  // Returns the bounds for this cell in (u,v)-space including padding.
+  boundRect(): R2Rect {
+    return this.bound
+  }
+
+  // Returns the (i,j) coordinates for the child cell at the given traversal position.
+  childIJ(pos: number): [number, number] {
+    const ij = posToIJ[this.orientation][pos]
+    return [ij >> 1, ij & 1]
+  }
+
+  // Returns the vertex where the space-filling curve enters this cell.
+  entryVertex(): Point {
+    let i = this.iLo
+    let j = this.jLo
+    if (this.orientation & INVERT_MASK) {
+      const ijSize = cellid.sizeIJ(this.level)
+      i += ijSize
+      j += ijSize
+    }
+    return Point.fromVector(faceSiTiToXYZ(cellid.face(this.id), 2 * i, 2 * j).vector.normalize())
+  }
+
+  // Returns the vertex where the space-filling curve exits this cell.
+  exitVertex(): Point {
+    let i = this.iLo
+    let j = this.jLo
+    const ijSize = cellid.sizeIJ(this.level)
+    if (this.orientation === 0 || this.orientation === SWAP_MASK + INVERT_MASK) {
+      i += ijSize
+    } else {
+      j += ijSize
+    }
+    return Point.fromVector(faceSiTiToXYZ(cellid.face(this.id), 2 * i, 2 * j).vector.normalize())
+  }
+
+  // Returns the smallest CellID that contains all descendants of this padded cell whose bounds intersect the given rect.
+  shrinkToFit(rect: R2Rect): CellID {
+    if (this.level === 0) {
+      if (rect.x.contains(0) || rect.y.contains(0)) return this.id
+    }
+
+    const ijSize = cellid.sizeIJ(this.level)
+    if (
+      rect.x.contains(stToUV(siTiToST(2 * this.iLo + ijSize))) ||
+      rect.y.contains(stToUV(siTiToST(2 * this.jLo + ijSize)))
+    ) {
+      return this.id
+    }
+
+    const padded = rect.expandedByMargin(this.padding + 1.5 * DBL_EPSILON)
+    let iMin = this.iLo
+    let jMin = this.jLo
+    let iXor = 0
+    let jXor = 0
+
+    if (iMin < cellid.stToIJ(uvToST(padded.x.lo))) {
+      iMin = cellid.stToIJ(uvToST(padded.x.lo))
+    }
+    if (this.iLo + ijSize - 1 <= cellid.stToIJ(uvToST(padded.x.hi))) {
+      iXor = iMin ^ (this.iLo + ijSize - 1)
+    } else {
+      iXor = iMin ^ cellid.stToIJ(uvToST(padded.x.hi))
+    }
+
+    if (jMin < cellid.stToIJ(uvToST(padded.y.lo))) {
+      jMin = cellid.stToIJ(uvToST(padded.y.lo))
+    }
+    if (this.jLo + ijSize - 1 <= cellid.stToIJ(uvToST(padded.y.hi))) {
+      jXor = jMin ^ (this.jLo + ijSize - 1)
+    } else {
+      jXor = jMin ^ cellid.stToIJ(uvToST(padded.y.hi))
+    }
+
+    const levelMSB = BigInt((iXor | jXor) << 1) + 1n
+    const level = MAX_LEVEL - findMSBSetNonZero64(levelMSB)
+    if (level <= this.level) return this.id
+
+    return cellid.parent(cellid.fromFaceIJ(cellid.face(this.id), iMin, jMin), level)
+  }
+}

s2/PaddedCell_test.ts

@@ -0,0 +1,94 @@
+import { test, describe } from 'node:test'
+import { deepEqual, equal, ok } from 'node:assert/strict'
+import { PaddedCell } from './PaddedCell'
+import { Point as R2Point } from '../r2/Point'
+import { Rect as R2Rect } from '../r2/Rect'
+import { Interval as R1Interval } from '../r1/Interval'
+import * as cellid from './cellid'
+import { oneIn, randomCellID, randomFloat64, randomUniformFloat64, randomUniformInt } from './testing'
+import { Cell } from './Cell'
+
+describe('s2.PaddedCell', () => {
+  test('methods', () => {
+    for (let i = 0; i < 1000; i++) {
+      const cid = randomCellID()
+      const padding = Math.pow(1e-15, randomFloat64())
+      const cell = Cell.fromCellID(cid)
+      const pCell = PaddedCell.fromCellID(cid, padding)
+      equal(cell.id, pCell.cellID())
+      equal(cellid.level(cell.id), pCell.levelValue())
+      equal(padding, pCell.paddingValue())
+      deepEqual(pCell.boundRect(), cell.boundUV().expandedByMargin(padding))
+      const r = R2Rect.fromPoints(cellid.centerUV(cell.id)).expandedByMargin(padding)
+      deepEqual(pCell.middleRect(), r)
+      deepEqual(cellid.point(cell.id), pCell.center())
+      if (cellid.isLeaf(cid)) continue
+      const children = cell.children()
+      for (let pos = 0; pos < 4; pos++) {
+        const [i, j] = pCell.childIJ(pos)
+        const cellChild = children[pos]
+        const pCellChild = PaddedCell.fromParentIJ(pCell, i, j)
+        equal(cellChild.id, pCellChild.cellID())
+        equal(cellid.level(cellChild.id), pCellChild.levelValue())
+        equal(padding, pCellChild.paddingValue())
+        deepEqual(pCellChild.boundRect(), cellChild.boundUV().expandedByMargin(padding))
+        const r = R2Rect.fromPoints(cellid.centerUV(cellChild.id)).expandedByMargin(padding)
+        ok(r.approxEqual(pCellChild.middleRect()))
+        deepEqual(cellid.point(cellChild.id), pCellChild.center())
+      }
+    }
+  })
+
+  test('entry/exit vertices', () => {
+    for (let i = 0; i < 1000; i++) {
+      const id = randomCellID()
+      const unpadded = PaddedCell.fromCellID(id, 0)
+      const padded = PaddedCell.fromCellID(id, 0.5)
+      deepEqual(unpadded.entryVertex(), padded.entryVertex())
+      deepEqual(unpadded.exitVertex(), padded.exitVertex())
+      deepEqual(PaddedCell.fromCellID(cellid.nextWrap(id), 0).entryVertex(), unpadded.exitVertex())
+      if (!cellid.isLeaf(id)) {
+        deepEqual(PaddedCell.fromCellID(cellid.children(id)[0], 0).entryVertex(), unpadded.entryVertex())
+        deepEqual(PaddedCell.fromCellID(cellid.children(id)[3], 0).exitVertex(), unpadded.exitVertex())
+      }
+    }
+  })
+
+  test('shrinkToFit', () => {
+    for (let iter = 0; iter < 1000; iter++) {
+      const result = randomCellID()
+      const resultUV = cellid.boundUV(result)
+      const sizeUV = resultUV.size()
+      const maxPadding = 0.5 * Math.min(sizeUV.x, sizeUV.y)
+      const padding = maxPadding * randomFloat64()
+      const maxRect = resultUV.expandedByMargin(-padding)
+      const a = new R2Point(
+        randomUniformFloat64(maxRect.x.lo, maxRect.x.hi),
+        randomUniformFloat64(maxRect.y.lo, maxRect.y.hi)
+      )
+      const b = new R2Point(
+        randomUniformFloat64(maxRect.x.lo, maxRect.x.hi),
+        randomUniformFloat64(maxRect.y.lo, maxRect.y.hi)
+      )
+      if (!cellid.isLeaf(result)) {
+        const useY = oneIn(2)
+        let center = cellid.centerUV(result).x
+        if (useY) center = cellid.centerUV(result).y
+        const shared = new R1Interval(center - padding, center + padding)
+        const intersected = useY ? shared.intersection(maxRect.y) : shared.intersection(maxRect.x)
+        const mid = randomUniformFloat64(intersected.lo, intersected.hi)
+        if (useY) {
+          a.y = randomUniformFloat64(maxRect.y.lo, mid)
+          b.y = randomUniformFloat64(mid, maxRect.y.hi)
+        } else {
+          a.x = randomUniformFloat64(maxRect.x.lo, mid)
+          b.x = randomUniformFloat64(mid, maxRect.x.hi)
+        }
+      }
+      const rect = R2Rect.fromPoints(a, b)
+      const initialID = cellid.parent(result, randomUniformInt(cellid.level(result) + 1))
+      const pCell = PaddedCell.fromCellID(initialID, padding)
+      equal(pCell.shrinkToFit(rect), result)
+    }
+  })
+})

s2/_index.ts

@@ -21,3 +21,4 @@ export { Rect } from './Rect'
 export { Region } from './Region'
 export { RegionCoverer, RegionCovererOptions, Coverer } from './RegionCoverer'
 export { Polyline } from './Polyline'
+export { PaddedCell } from './PaddedCell'

s2/cellid.ts

@@ -9,9 +9,44 @@ import { Vector } from '../r3/Vector'
 import { Interval } from '../r1/Interval'
 import { Rect } from '../r2/Rect'
 import { Point as R2Point } from '../r2/Point'
+import { Rect as R2Rect } from '../r2/Rect'
 
+/**
+ * Uniquely identifies a cell in the S2 cell decomposition.
+ * The most significant 3 bits encode the face number (0-5). The
+ * remaining 61 bits encode the position of the center of this cell
+ * along the Hilbert curve on that face. The zero value and the value
+ * (1<<64)-1 are invalid cell IDs. The first compares less than any
+ * valid cell ID, the second as greater than any valid cell ID.
+ *
+ * Sequentially increasing cell IDs follow a continuous space-filling curve
+ * over the entire sphere. They have the following properties:
+ *
+ *   - The ID of a cell at level k consists of a 3-bit face number followed
+ *     by k bit pairs that recursively select one of the four children of
+ *     each cell. The next bit is always 1, and all other bits are 0.
+ *     Therefore, the level of a cell is determined by the position of its
+ *     lowest-numbered bit that is turned on (for a cell at level k, this
+ *     position is 2 * (MaxLevel - k)).
+ *
+ *   - The ID of a parent cell is at the midpoint of the range of IDs spanned
+ *     by its children (or by its descendants at any level).
+ *
+ * Leaf cells are often used to represent points on the unit sphere, and
+ * this type provides methods for converting directly between these two
+ * representations. For cells that represent 2D regions rather than
+ * discrete point, it is better to use Cells.
+ */
 export type CellID = bigint
 
+/**
+ * An invalid cell ID guaranteed to be larger than any
+ * valid cell ID. It is used primarily by ShapeIndex. The value is also used
+ * by some S2 types when encoding data.
+ * Note that the sentinel's RangeMin == RangeMax == itself.
+ */
+export const SentinelCellID = ~0n
+
 /**
  * Returns the cube face for this cell id, in the range [0,5].
  */
@@ -292,6 +327,12 @@ export const centerUV = (ci: CellID): R2Point => {
   return new R2Point(stToUV(siTiToST(si)), stToUV(siTiToST(ti)))
 }
 
+// boundUV returns the bound of this CellID in (u,v)-space.
+export const boundUV = (ci: CellID): R2Rect => {
+  const { i, j } = faceIJOrientation(ci)
+  return ijLevelToBoundUV(i, j, level(ci))
+}
+
 /**
  * Returns a leaf cell given its cube face (range 0..5) and IJ coordinates.
  * @category Constructors
@@ -570,6 +611,26 @@ export const prev = (ci: CellID): CellID => {
   return ci - (lsb(ci) << 1n)
 }
 
+/**
+ * Returns the next cell along the Hilbert curve, wrapping from last to
+ * first as necessary. This should not be used with ChildBegin and ChildEnd.
+ */
+export const nextWrap = (ci: CellID): CellID => {
+  const n = next(ci)
+  if (n < WRAP_OFFSET) return n
+  return n - WRAP_OFFSET
+}
+
+/**
+ * Returns the previous cell along the Hilbert curve, wrapping around from
+ * first to last as necessary. This should not be used with ChildBegin and ChildEnd.
+ */
+export const prevWrap = (ci: CellID): CellID => {
+  const p = prev(ci)
+  if (p < WRAP_OFFSET) return p
+  return p + WRAP_OFFSET
+}
+
 /**
  * Returns the level of the common ancestor of the two S2 CellIDs.
  * @param other The other CellID to compare with.