diff --git a/dev-docs/specs/2026-05-27-globe-reprojector.md b/dev-docs/specs/2026-05-27-globe-reprojector.md new file mode 100644 index 00000000..1ea9548e --- /dev/null +++ b/dev-docs/specs/2026-05-27-globe-reprojector.md @@ -0,0 +1,377 @@ +# GlobeReprojector: sphere-aware adaptive meshing + +- **Date:** 2026-05-27 +- **Issues:** Follow-up to [#563](https://github.com/developmentseed/deck.gl-raster/pull/563) (Initial Globe view support); replaces the throwaway uniform-grid scaffold from [#563](https://github.com/developmentseed/deck.gl-raster/pull/563). +- **Status:** Draft (pending review) + +## Summary + +Replace the throwaway uniform-grid globe mesh (`buildUniformGridMesh`) with a +`GlobeReprojector` that adds a **curvature (sag)** term to the mesh-refinement +metric, so the adaptive Delatin mesh stays accurate on the sphere instead of +faceting at low zoom. + +`GlobeReprojector` is a subclass of `RasterReprojector` that lives in +**`deck.gl-raster`**, not the generic `raster-reproject` package — because it +legitimately needs deck.gl [`GlobeView`](https://deck.gl/docs/api-reference/core/globe-view)'s projection to measure sag. The base +`RasterReprojector` gains a handful of small, **behavior-preserving** +extensibility seams and learns nothing about spheres. + +## Background / Problem + +`RasterReprojector` ([Delatin](https://github.com/mapbox/delatin)) refines triangles by **reprojection error**, +measured in **input raster pixels** — the deviation between the linearly +interpolated sample and the exact reprojection ([delatin.ts](../../packages/raster-reproject/src/delatin.ts) `_findReprojectionCandidate`). This +metric deliberately mirrors GDAL's [approximate transformer](https://gdal.org/en/stable/programs/gdalwarp.html#approximate-transformation) (`gdalwarp -et`), and +**must be preserved** for the planar case. + +That metric is **blind to sphere curvature**. For an EPSG:4326 source the +reprojection is near-linear, so pixel error is ~0 and Delatin emits just two +triangles. On a globe those two triangles chord straight through the sphere and +visibly facet at low zoom. + +The current stopgap is [`buildUniformGridMesh`](../../packages/deck.gl-raster/src/globe-grid-mesh.ts) — a 32×32 uniform grid per +tile, explicitly marked THROWAWAY (see the [globe-view design](./2026-05-21-globe-view-design.md)). It is dense everywhere regardless of where +curvature actually matters. This spec replaces it. + +### Key facts that shape the design + +- **The rendered mesh position is lng/lat, not sphere xyz.** deck.gl `GlobeView` + projects lng/lat → sphere in the vertex shader. So the mesh `POSITION` + attribute stays lng/lat (fp64 hi/lo), exactly as today. The sphere positions + we compute are used **only** by the sag metric — they are never emitted. + `GlobeReprojector` therefore produces the **same mesh format** as + `RasterReprojector`; `reprojectorToMesh` and the debug layer work on it + unchanged. +- **The reprojector is viewport-independent, and stays that way.** Each tile + gets its own reprojector → its own mesh, built once. Zoom changes which + tiles/overviews are resident; it never regenerates a tile's mesh. We do **not** + make the reprojector zoom-aware. + +## Goals / Non-goals + +**Goals** + +- Sphere-aware adaptive refinement: dense where the surface is curved, sparse + where it is flat. +- Zoom-independent: one mesh per tile, never rebuilt on zoom. +- Preserve the planar/GDAL pixel metric byte-for-byte. +- Keep `raster-reproject` free of any globe/sphere/deck.gl concept. + +**Non-goals (explicitly deferred, but not foreclosed)** + +- **GCP-based seeding.** The `InitialTriangulation` seed mechanism added in #574 + (constructor `options.initialTriangulation` + the `triangulateRectangle` + helper) is the natural future home for triangulating from Ground Control + Points, but we do not build that here. The globe uses the default full-image + seed. +- Antimeridian handling (tracked separately). +- Pole singularity. Near the poles the reprojection never fully converges; the + existing `maxIterations` safeguard still applies. Sag refinement does not make + this worse. + +## Design overview + +Three pieces: + +1. **`RasterReprojector` (raster-reproject)** — add behavior-preserving + extensibility seams. No globe concepts. +2. **`GlobeReprojector` (deck.gl-raster)** — subclass adding the sphere cache, + the sag metric, and a two-tolerance `run`. +3. **`RasterLayer` (deck.gl-raster)** — the globe branch builds a + `GlobeReprojector` through the existing `reprojectorToMesh` path, and the + uniform-grid hack is deleted. + +## Component 1 — `RasterReprojector` extensibility seams + +All edits are **behavior-preserving**: the planar output is identical, existing +tests pass unchanged. The seams are projection-agnostic (equally useful for the +future GCP case or any custom metric). + +### 1a. Make `_addPoint` overridable + +`_addPoint` is `private` today ([delatin.ts:465](../../packages/raster-reproject/src/delatin.ts#L465)). Change to +`protected`. No logic change. + +### 1b. Extract the scoring seam + +Today the per-sample pixel error is computed inline in the sampling loop of +`_findReprojectionCandidate` ([delatin.ts:387](../../packages/raster-reproject/src/delatin.ts#L387)). Extract the per-sample scalar +into a `protected` method, default = today's pixel error: + +```ts +/** + * Per-sample refinement error. The triangle's queue priority is the max of + * this over its sample points, and the split candidate is the argmax sample. + * + * Default: the reprojection error in input pixels (GDAL-like). Subclasses may + * return a different scalar; the value's units define what `run`'s tolerance + * means. + */ +protected _sampleError(ctx: SampleErrorContext): number { + return ctx.pixelError; // hypot(dx, dy), unchanged behavior +} +``` + +`SampleErrorContext` carries everything the base already computes for the +sample, so the seam never recomputes: the barycentric weights, the three +triangle vertex indices, the interpolated output position, and the base's +`pixelError`. The base loop calls `_sampleError`, tracks its max + the +corresponding uv as the split candidate, and pushes the max to the queue — +structurally identical to today, just routed through the seam. + +### 1c. Keep eager seeding; the init-order hazard moves to the subclass + +#574 already extracted seeding into `private _seed(seed: InitialTriangulation)` +([delatin.ts:221](../../packages/raster-reproject/src/delatin.ts#L221)), still +called **eagerly** from the constructor: `_seed(...)` → `_addPoint` ×N → +`_flush` ([delatin.ts:211-212](../../packages/raster-reproject/src/delatin.ts#L211-L212)). +We keep that as-is — **no constructor lifecycle change.** The globe uses the +default full-image seed (the Web-Mercator latitude clamp from #574 is wired only +into the Mercator branch; the globe shows the poles), so `_seed` stays `private`. + +The catch: because seeding runs inside `super()`, a subclass's `_addPoint` +override fires before the subclass's own fields exist, and the construction-time +`_flush` runs `_sampleError` before any run-time tolerances are set. Rather than +make seeding lazy, **`GlobeReprojector` absorbs both** (see 2a/2b/2d/2e): +`renderPositions` is lazy-created inside the `_addPoint` override, and `run()` +re-scores via `_reevaluate()` after setting tolerances, discarding the +provisional construction-time scores. So the base needs only **1a/1b/1d** — +three small `protected` edits, none touching the constructor. + +### 1d. `_reevaluate()` for tolerance changes + +```ts +/** Re-score every existing triangle through `_sampleError` and rebuild the + * queue. Used when run-time scoring inputs (e.g. tolerances) change between + * runs so a stale, differently-scored queue isn't reused. */ +protected _reevaluate(): void { /* mark all triangles pending, _flush */ } +``` + +Projection-agnostic. The base never calls it (its queue holds raw, +tolerance-free pixel error); it exists for subclasses whose scoring depends on +run-time inputs. + +### Public API + +No new public exports from `raster-reproject`. The seams are `protected`; +`RasterReprojector` is already exported and subclassable. The package boundary +**enforces** the clean seam: a subclass in `deck.gl-raster` physically cannot +touch the base's `private` queue internals, so the extension surface stays the +narrow set of `protected` hooks above. + +## Component 2 — `GlobeReprojector` (deck.gl-raster) + +New file `packages/deck.gl-raster/src/globe-reprojector.ts`. Internal to the +package (not exported from `index.ts`) until there's an external use case. + +### 2a. The render-position cache + +```ts +/** Sphere positions (deck.gl GlobeView common space), stride 3, indexed by + * vertex — parallel to the base's `exactOutputPositions` (stride 2). + * No initializer: the base constructor's eager `_seed` calls `_addPoint` + * before this subclass's field initializers would run, so it is created + * lazily in the override (2b). */ +renderPositions!: number[]; +``` + +A definite-assignment declaration (no initializer to clobber); the array is +lazy-created via `??=` in `_addPoint` (2b) on the first vertex — which always +precedes any read of it. + +### 2b. `_addPoint` override + +```ts +protected override _addPoint(u: number, v: number): number { + this.renderPositions ??= []; // created on the first (construction-time) vertex + const i = super._addPoint(u, v); // pushes uv + exact lng/lat as today + const lng = this.exactOutputPositions[2 * i]!; + const lat = this.exactOutputPositions[2 * i + 1]!; + const [x, y, z] = this._projectToSphere(lng, lat); + this.renderPositions.push(x, y, z); // slot i, stride 3 + return i; +} +``` + +`_projectToSphere` is deck.gl `GlobeView`'s lng/lat → common-space sphere +mapping (view-independent, constant radius `R`). Because we're in +`deck.gl-raster`, we use deck.gl's projection directly — no injection. + +### 2c. The sag metric — definition + +**Sag is the gap between the flat triangle the GPU draws and the round sphere +surface deck.gl wants.** deck.gl fills a triangle as a flat facet between its +three projected vertices; it does not re-project interior pixels. So the facet +dips inside the sphere. The deeper the dip, the more faceted it looks. + +Decompose an interior point's error into two orthogonal components: + +- **Radial** (toward/away from the globe center): the facet sits below the + sphere. Pure geometry — exists even with a perfect texture. **This is sag.** +- **Tangential** (along the surface): the texture lands slightly off because the + CRS→lng/lat map is nonlinear. **This is the existing pixel error.** + +They don't overlap, which is exactly why two independent tolerances are clean. + +**Formula.** Every flat-raster vertex sits at the same radius `R` from the globe +center (z=0, no elevation). The rendered interior point is the barycentric mix +of the three cached corner positions, which lands at radius `< R`. The closest +point on a sphere is always radial, so the perpendicular distance from the facet +point to the sphere is exact: + +``` +sag(sample) = R − | barycentricMix(renderPos_a, renderPos_b, renderPos_c) | +``` + +No extra projection, no lng/lat interpolation — just the cached corner positions +and one vector length. Normalizing by `R` gives a dimensionless **relative dip** +(`1 − |mix|/R`), making the sag tolerance resolution- and zoom-independent. + +> **Rejected alternative:** the full 3-D distance to `sphere(exactReproject(sample))`. +> It costs an extra projection per sample *and* folds the tangential/texture +> error back into sag, double-counting the pixel metric. The radial measure is +> cheaper and properly orthogonal. + +### 2d. The scoring seam override — two tolerances, normalized + +The priority queue needs **one** scalar per triangle to decide what to split +next. Pixel error (input pixels) and sag (relative dip) are different units, so +rank by how badly each blows its own budget: + +```ts +protected override _sampleError(ctx: SampleErrorContext): number { + const sag = this._relativeSag(ctx); // (R − |mix|)/R from cached corners + weights + return Math.max( + ctx.pixelError / this._pixelTolerance, + sag / this._sagTolerance, + ); +} +``` + +`max(px/pxTol, sag/sagTol) > 1` is identically `px > pxTol OR sag > sagTol` — so +this **is** two independent tolerances in their native units; the division only +makes them rankable in one queue. The split candidate is the sample maximizing +this normalized value (so the most-over-budget point is added first). + +The base planar metric is this with `sagTolerance = ∞` (sag term vanishes, +score = `px/pxTol`). `GlobeReprojector` strictly generalizes the base. + +Because the base seeds eagerly, this seam also runs **once during construction**, +before `run()` sets the tolerances. `_pixelTolerance`/`_sagTolerance` therefore +read as defaults on that pass, and `run()`'s `_reevaluate()` recomputes those +provisional scores with the real tolerances (2e). Those construction-time values +are never consumed by refinement before the recompute. + +### 2e. `run` — tolerances are per-run inputs + +Tolerance is a property of *a run*, not of the reprojector, and may change +between runs. #574 kept base `run(maxError, { maxIterations })` +([delatin.ts:246](../../packages/raster-reproject/src/delatin.ts#L246)), so the +globe overrides it **compatibly** — no base API change. `maxError` stays the +pixel tolerance; `sagTolerance` joins the options bag: + +```ts +// base (raster-reproject), unchanged +run(maxError = DEFAULT_MAX_ERROR, { maxIterations = 10_000 } = {}): void + +// GlobeReprojector — signature-compatible override (options bag extended) +run( + maxError = DEFAULT_MAX_ERROR, + { sagTolerance = DEFAULT_SAG_TOLERANCE, maxIterations = 10_000 } = {}, +): void { + this._pixelTolerance = maxError; + this._sagTolerance = sagTolerance; + this._reevaluate(); // re-score the seed with the real tolerances + while (this.getMaxError() > 1) { // seam is normalized → fixed threshold 1 + this.refine(); + // ...same maxIterations safeguard as base + } +} +``` + +`_pixelTolerance`/`_sagTolerance` are **transient run state** (set at `run` +entry, read by the seam during that run), not persistent configuration — the +public contract is still "supply tolerances per `run` call." `_reevaluate()` +re-scores the existing triangles (just the seed on the first run) so the queue +reflects the current tolerances rather than the construction-time defaults (2d). + +## Component 3 — `RasterLayer` wiring + +In `_generateMesh` ([raster-layer.ts:205-228](../../packages/deck.gl-raster/src/raster-layer.ts#L205-L228)), the globe branch currently calls +`buildUniformGridMesh`. Replace with a `GlobeReprojector`, then reuse the +**existing** `reprojectorToMesh` path (positions = lng/lat, unchanged): + +```ts +const reprojector = isGlobe + ? new GlobeReprojector(reprojectionFns, width + 1, height + 1) + : new RasterReprojector(reprojectionFns, width + 1, height + 1); +// maxError is the pixel tolerance for both; the globe adds sagTolerance +if (reprojector instanceof GlobeReprojector) { + reprojector.run(maxError, { sagTolerance }); +} else { + reprojector.run(maxError); +} +const mesh = reprojectorToMesh(reprojector); +``` + +- `isGlobe` detection stays as-is (`viewport.resolution !== undefined`). +- `reprojector` is stored in state for both branches → the debug layer now works + on the globe (today it's `undefined` for globe). +- **Delete** `globe-grid-mesh.ts` and its `buildUniformGridMesh` import. +- Add a `sagTolerance` prop to `RasterLayerProps` (default tuned visually), + threaded down from the COG/globe example. + +## Data flow + +``` +tile load → RasterLayer._generateMesh + → new GlobeReprojector(reprojectionFns, w+1, h+1) + constructor → _seed(default) → _addPoint ×4 (lazy-creates + caches sphere xyz) + → _flush → _sampleError (provisional: tolerances still default) + → run(maxError /* = pixelTolerance */, { sagTolerance }) + set tolerances → _reevaluate (re-score the seed triangles) + refine loop: _step → _addPoint (caches) ; _flush → _findReprojectionCandidate + per sample → _sampleError = max(px/pxTol, sag/sagTol) [sag from renderPositions] + until getMaxError() ≤ 1 + → reprojectorToMesh (lng/lat positions, fp64 hi/lo) → MeshTextureLayer + → GlobeView projects lng/lat → sphere in the vertex shader +``` + +## Error handling / edge cases + +- **Poles:** unchanged non-convergence risk; `maxIterations` warns and bails + ([delatin.ts:257-264](../../packages/raster-reproject/src/delatin.ts#L257-L264)). Sag pushes more refinement near curved regions but does + not introduce a new failure mode. +- **Degenerate `R`:** all vertices share `R`; deriving it from the projection is + safe. Guard against `sagTolerance ≤ 0` (throw, like `maxError ≤ 0` today). +- **`renderPositions` / `exactOutputPositions` alignment:** guaranteed by + pushing in lockstep inside `_addPoint`; both indexed by vertex ordinal. + +## Testing strategy + +- **Base unchanged:** existing `raster-reproject` tests pass after the seam + extraction (behavior-preserving; the constructor and seeding are untouched). +- **Init order:** constructing a `GlobeReprojector` (eager seed) populates + `renderPositions` via the `??=` guard, and the construction-time scoring pass + does not throw with tolerances unset; `run()`'s `_reevaluate()` then yields the + correct mesh. +- **Sag formula:** given three known sphere corner positions, `sag` equals the + analytic radial gap (e.g. a chord subtending angle θ → `R(1−cos(θ/2))` at the + midpoint). +- **Sag-driven refinement:** a large EPSG:4326 tile (pixel error ≈ 0) refines + into many triangles under `GlobeReprojector`, vs. 2 under the base — proving + the curvature term drives subdivision the pixel metric misses. +- **Generalization:** `GlobeReprojector` with `sagTolerance = ∞` reproduces the + base mesh (pixel-only). +- **Tolerance change:** `_reevaluate()` re-scores existing triangles; a second + `run` with a tighter `sagTolerance` yields a denser mesh. + +## Open questions / to confirm during implementation + +1. **`_projectToSphere` source.** Confirm the exact deck.gl API for a + view-independent lng/lat → common-space sphere position (and `R`). If deck.gl + only exposes it via a viewport, use a canonical globe viewport or replicate + the small spherical formula. Design is unaffected. +2. **Default `sagTolerance`.** Needs visual tuning on the globe example; ship a + sensible default (relative dip) and expose the prop. diff --git a/packages/raster-reproject/src/delatin.ts b/packages/raster-reproject/src/delatin.ts index 6183e20e..541e7bd9 100644 --- a/packages/raster-reproject/src/delatin.ts +++ b/packages/raster-reproject/src/delatin.ts @@ -131,6 +131,20 @@ export interface ReprojectionFns { inverseReproject(x: number, y: number): [number, number]; } +/** + * Per-sample context handed to {@link RasterReprojector._sampleError}. Carries + * everything the base already computed for one barycentric sample of a + * triangle, so a scoring override never recomputes it. + */ +export interface SampleErrorContext { + /** Reprojection error at this sample, in input pixels — `hypot(dx, dy)`. */ + pixelError: number; + /** Barycentric weights of the sample within the triangle (sum to 1). */ + weights: readonly [number, number, number]; + /** The triangle's three vertex indices (index into `uvs`/2, `renderPositions`/3). */ + vertices: readonly [number, number, number]; +} + /** * RasterReprojector performs a Delaunay triangulation-based reprojection of a * raster image. @@ -271,6 +285,24 @@ export class RasterReprojector { this._flush(); } + /** + * Re-score every triangle through {@link _sampleError} and rebuild the + * priority queue. Use when a run-time scoring input changed since the last + * pass so the queue isn't a mix of stale and current scores. The queue is + * cleared first because `_queuePush` appends unconditionally. + */ + protected _reevaluate(): void { + this._queue.length = 0; + this._errors.length = 0; + const numTriangles = this.triangles.length / 3; + for (let t = 0; t < numTriangles; t++) { + this._queueIndices[t] = -1; + this._pending[t] = t; + } + this._pendingLen = numTriangles; + this._flush(); + } + // max error of the current mesh getMaxError(): number { return this._errors[0]!; @@ -285,6 +317,18 @@ export class RasterReprojector { this._pendingLen = 0; } + /** + * Per-sample refinement error. The triangle's queue priority is the max of + * this over its sample points, and the split candidate is the argmax sample. + * + * Default: the reprojection error in input pixels (matches GDAL's approximate + * transformer). Subclasses may return a different scalar; the value's units + * define what `run`'s tolerance threshold means. + */ + protected _sampleError(ctx: SampleErrorContext): number { + return ctx.pixelError; + } + /** * Conversion of upstream's `_findCandidate` for reprojection error handling. * @@ -294,9 +338,13 @@ export class RasterReprojector { */ private _findReprojectionCandidate(t: number): void { // Find the three vertices of this triangle - const a = 2 * this.triangles[t * 3 + 0]!; - const b = 2 * this.triangles[t * 3 + 1]!; - const c = 2 * this.triangles[t * 3 + 2]!; + const v0 = this.triangles[t * 3 + 0]!; + const v1 = this.triangles[t * 3 + 1]!; + const v2 = this.triangles[t * 3 + 2]!; + const a = 2 * v0; + const b = 2 * v1; + const c = 2 * v2; + const vertices: [number, number, number] = [v0, v1, v2]; // Get the UV coordinates of each vertex const p0u = this.uvs[a]!; @@ -384,10 +432,17 @@ export class RasterReprojector { // 4. error in pixel space const dx = pixelExactX - pixelSampled[0]; const dy = pixelExactY - pixelSampled[1]; - const err = Math.hypot(dx, dy); + const pixelError = Math.hypot(dx, dy); + + // Route through the scoring seam (default: the pixel error itself). + const score = this._sampleError({ + pixelError, + weights: samplePoint, + vertices, + }); - if (err > maxError) { - maxError = err; + if (score > maxError) { + maxError = score; maxErrorU = uvSampleU; maxErrorV = uvSampleV; } @@ -462,7 +517,7 @@ export class RasterReprojector { } // add coordinates for a new vertex - private _addPoint(u: number, v: number): number { + protected _addPoint(u: number, v: number): number { const i = this.uvs.length >> 1; this.uvs.push(u, v); diff --git a/packages/raster-reproject/src/index.ts b/packages/raster-reproject/src/index.ts index 78df1855..817806bc 100644 --- a/packages/raster-reproject/src/index.ts +++ b/packages/raster-reproject/src/index.ts @@ -1,2 +1,6 @@ -export type { InitialTriangulation, ReprojectionFns } from "./delatin.js"; +export type { + InitialTriangulation, + ReprojectionFns, + SampleErrorContext, +} from "./delatin.js"; export { RasterReprojector, triangulateRectangle } from "./delatin.js"; diff --git a/packages/raster-reproject/tests/reevaluate.test.ts b/packages/raster-reproject/tests/reevaluate.test.ts new file mode 100644 index 00000000..51ec08ec --- /dev/null +++ b/packages/raster-reproject/tests/reevaluate.test.ts @@ -0,0 +1,33 @@ +import { describe, expect, it } from "vitest"; +import type { ReprojectionFns, SampleErrorContext } from "../src/delatin.js"; +import { RasterReprojector } from "../src/delatin.js"; + +const identity: ReprojectionFns = { + forwardTransform: (x, y) => [x, y], + inverseTransform: (x, y) => [x, y], + forwardReproject: (x, y) => [x, y], + inverseReproject: (x, y) => [x, y], +}; + +// A reprojector whose per-sample error is a tunable constant, so we can change +// a scoring input after construction and verify _reevaluate re-scores cleanly. +class TunableReprojector extends RasterReprojector { + error = 0; + protected override _sampleError(_ctx: SampleErrorContext): number { + return this.error ?? 0; // ?? 0: this field is undefined during super()'s seed flush + } + reevaluate(): void { + this._reevaluate(); + } +} + +describe("RasterReprojector._reevaluate", () => { + it("re-scores existing triangles when a scoring input changes", () => { + const r = new TunableReprojector(identity, 64, 64); + expect(r.getMaxError()).toBe(0); // seeded with error = 0 + + r.error = 5; + r.reevaluate(); + expect(r.getMaxError()).toBe(5); // re-scored; no stale/duplicate entries + }); +});