diff --git a/packages/deck.gl-raster/src/raster-layer.ts b/packages/deck.gl-raster/src/raster-layer.ts index 6534c7c1..bd3181d6 100644 --- a/packages/deck.gl-raster/src/raster-layer.ts +++ b/packages/deck.gl-raster/src/raster-layer.ts @@ -247,6 +247,68 @@ export class RasterLayer extends CompositeLayer { // To account for this, we add 1 to both width and height when generating // the mesh. This also solves obvious gaps in between neighboring tiles in // the COGLayer. + // + // For tiles that straddle the CRS domain boundary (e.g. a Mollweide tile + // with one corner outside the ellipse), the adaptive reprojector wastes its + // entire iteration budget refining OOD-vertex triangles, leaving the valid + // area under-refined (reprojection error 100s of pixels). Detect these + // "border tiles" up front and use a dense uniform grid instead: OOD-vertex + // triangles are simply filtered out, producing a clean domain-edge cutoff + // at predictable 1/(BORDER_GRID_SIZE) tile-fraction resolution. + const { forwardTransform, forwardReproject, inverseReproject } = + reprojectionFns; + const borderCorners: [number, number][] = [ + [0, 0], + [1, 0], + [0, 1], + [1, 1], + ]; + // Tile width in projected CRS units for the round-trip error threshold. + const tileWidthCRS = Math.abs( + forwardTransform(width, 0)[0] - forwardTransform(0, 0)[0], + ); + const roundTripThreshold = tileWidthCRS * 0.01; + const isBorderTile = borderCorners.some(([u, v]) => { + const [ix, iy] = forwardTransform(u * width, v * height); + const [ox, oy] = forwardReproject(ix, iy); + if (!Number.isFinite(ox) || !Number.isFinite(oy)) { + return true; + } + // Round-trip check: OOD corners where forwardReproject clamps to a + // finite domain-boundary position instead of returning NaN. The clamped + // value re-projects back to a CRS position far from the original, + // revealing the OOD via large round-trip error. + if (inverseReproject && roundTripThreshold > 0) { + const [ix2, iy2] = inverseReproject(ox, oy); + if (!Number.isFinite(ix2) || !Number.isFinite(iy2)) { + return true; + } + const err = Math.sqrt( + (ix2 - ix) * (ix2 - ix) + (iy2 - iy) * (iy2 - iy), + ); + if (err > roundTripThreshold) { + return true; + } + } + return false; + }); + if (isBorderTile) { + const { indices, positions64High, positions64Low, texCoords } = + buildClippedGridMesh(reprojectionFns, width + 1, height + 1); + this.setState({ + reprojector: undefined, + mesh: { + indices: { value: indices, size: 1 }, + attributes: { + POSITION: { value: positions64High, size: 3 }, + TEXCOORD_0: { value: texCoords, size: 2 }, + }, + }, + positions64Low, + }); + return; + } + const reprojector = new RasterReprojector( reprojectionFns, width + 1, @@ -379,19 +441,38 @@ function reprojectorToMesh(reprojector: RasterReprojector): { const texCoords = new Float32Array(reprojector.uvs); const positions = new Float64Array(numVertices * 3); + // Track which vertices are outside the CRS domain (NaN output position). + const isOOD = new Uint8Array(numVertices); for (let i = 0; i < numVertices; i++) { - positions[i * 3] = reprojector.exactOutputPositions[i * 2]!; - positions[i * 3 + 1] = reprojector.exactOutputPositions[i * 2 + 1]!; + const x = reprojector.exactOutputPositions[i * 2]!; + const y = reprojector.exactOutputPositions[i * 2 + 1]!; + positions[i * 3] = x; + positions[i * 3 + 1] = y; // z (flat on the ground) positions[i * 3 + 2] = 0; + if (!Number.isFinite(x) || !Number.isFinite(y)) { + isOOD[i] = 1; + } + } + + // Filter out any triangle that contains an out-of-domain vertex. This + // clips the rendered mesh cleanly at the CRS boundary without relying on + // undefined GPU NaN behaviour. + const allTriangles = reprojector.triangles; + const filteredTriangles: number[] = []; + for (let t = 0; t * 3 < allTriangles.length; t++) { + const a = allTriangles[t * 3]!; + const b = allTriangles[t * 3 + 1]!; + const c = allTriangles[t * 3 + 2]!; + if (!isOOD[a] && !isOOD[b] && !isOOD[c]) { + filteredTriangles.push(a, b, c); + } } // Split the float64 positions into high and low parts for fp64 emulation in // the shader. const [positions64Low, positions64High] = splitFloat64Array(positions); - - // TODO: Consider using 16-bit indices if the mesh is small enough - const indices = new Uint32Array(reprojector.triangles); + const indices = new Uint32Array(filteredTriangles); return { indices, @@ -400,3 +481,96 @@ function reprojectorToMesh(reprojector: RasterReprojector): { texCoords, }; } + +/** + * Build a dense uniform grid mesh over a tile, filtering out any triangle + * that contains a vertex whose output position is outside the CRS domain + * (i.e. forwardReproject returned NaN or fails the round-trip check). Used + * for "border tiles" where the CRS boundary passes through the tile; the + * adaptive reprojector handles these poorly because OOD-vertex triangles + * consume its entire iteration budget and leave the valid area under-refined. + */ +const BORDER_GRID_SIZE = 64; + +function buildClippedGridMesh( + reprojectionFns: ReprojectionFns, + width: number, + height: number, + gridSize = BORDER_GRID_SIZE, +): { + indices: Uint32Array; + positions64High: Float32Array; + positions64Low: Float32Array; + texCoords: Float32Array; +} { + const { forwardTransform, forwardReproject, inverseReproject } = + reprojectionFns; + const cols = gridSize; + const rows = gridSize; + const numVerts = (cols + 1) * (rows + 1); + + const positions = new Float64Array(numVerts * 3); + const texCoords = new Float32Array(numVerts * 2); + const isOOD = new Uint8Array(numVerts); + + // Tile CRS width for round-trip error threshold (same logic as isBorderTile). + const tileWidthCRS = Math.abs( + forwardTransform(width - 1, 0)[0] - forwardTransform(0, 0)[0], + ); + const roundTripThreshold = tileWidthCRS * 0.01; + + let vi = 0; + for (let r = 0; r <= rows; r++) { + for (let c = 0; c <= cols; c++) { + const u = c / cols; + const v = r / rows; + const pixelX = u * (width - 1); + const pixelY = v * (height - 1); + const [ix, iy] = forwardTransform(pixelX, pixelY); + const [ox, oy] = forwardReproject(ix, iy); + positions[vi * 3] = ox; + positions[vi * 3 + 1] = oy; + positions[vi * 3 + 2] = 0; + texCoords[vi * 2] = u; + texCoords[vi * 2 + 1] = v; + let ood = !Number.isFinite(ox) || !Number.isFinite(oy); + if (!ood && inverseReproject && roundTripThreshold > 0) { + const [ix2, iy2] = inverseReproject(ox, oy); + if (!Number.isFinite(ix2) || !Number.isFinite(iy2)) { + ood = true; + } else { + const err = Math.sqrt( + (ix2 - ix) * (ix2 - ix) + (iy2 - iy) * (iy2 - iy), + ); + if (err > roundTripThreshold) { + ood = true; + } + } + } + if (ood) { + isOOD[vi] = 1; + } + vi++; + } + } + + const filteredIndices: number[] = []; + for (let r = 0; r < rows; r++) { + for (let c = 0; c < cols; c++) { + const i0 = r * (cols + 1) + c; + const i1 = i0 + 1; + const i2 = i0 + (cols + 1); + const i3 = i2 + 1; + if (!isOOD[i0] && !isOOD[i2] && !isOOD[i1]) { + filteredIndices.push(i0, i2, i1); + } + if (!isOOD[i1] && !isOOD[i2] && !isOOD[i3]) { + filteredIndices.push(i1, i2, i3); + } + } + } + + const [positions64Low, positions64High] = splitFloat64Array(positions); + const indices = new Uint32Array(filteredIndices); + return { indices, positions64High, positions64Low, texCoords }; +} diff --git a/packages/deck.gl-raster/src/raster-tileset/affine-tileset-level.ts b/packages/deck.gl-raster/src/raster-tileset/affine-tileset-level.ts index 671a1f0d..9f60486d 100644 --- a/packages/deck.gl-raster/src/raster-tileset/affine-tileset-level.ts +++ b/packages/deck.gl-raster/src/raster-tileset/affine-tileset-level.ts @@ -44,12 +44,16 @@ export class AffineTilesetLevel implements RasterTilesetLevel { private readonly _affine: Affine; private readonly _invAffine: Affine; + private readonly _arrayWidth: number; + private readonly _arrayHeight: number; constructor(options: AffineTilesetLevelOptions) { this._affine = options.affine; this._invAffine = affine.invert(options.affine); this.tileWidth = options.tileWidth; this.tileHeight = options.tileHeight; + this._arrayWidth = options.arrayWidth; + this._arrayHeight = options.arrayHeight; this.matrixWidth = Math.ceil(options.arrayWidth / options.tileWidth); this.matrixHeight = Math.ceil(options.arrayHeight / options.tileHeight); @@ -78,12 +82,19 @@ export class AffineTilesetLevel implements RasterTilesetLevel { const tw = this.tileWidth; const th = this.tileHeight; const af = this._affine; - + // Clip to actual array extent so corners of the last tile row/column don't + // extrapolate past the data boundary. For projections whose valid domain + // isn't axis-aligned with the affine (Mollweide, Sinusoidal, Equal Earth), + // extrapolated corners fall outside the CRS domain — proj4 maps them to the + // pole, collapsing every such corner onto a single ±85.05° Mercator line and + // producing a zero-height bounding volume that causes every tile to be culled. + const right = Math.min((col + 1) * tw, this._arrayWidth); + const bottom = Math.min((row + 1) * th, this._arrayHeight); return { topLeft: affine.apply(af, col * tw, row * th), - topRight: affine.apply(af, (col + 1) * tw, row * th), - bottomLeft: affine.apply(af, col * tw, (row + 1) * th), - bottomRight: affine.apply(af, (col + 1) * tw, (row + 1) * th), + topRight: affine.apply(af, right, row * th), + bottomLeft: affine.apply(af, col * tw, bottom), + bottomRight: affine.apply(af, right, bottom), }; } diff --git a/packages/deck.gl-raster/src/raster-tileset/raster-tile-traversal.ts b/packages/deck.gl-raster/src/raster-tileset/raster-tile-traversal.ts index fac63439..230080a3 100644 --- a/packages/deck.gl-raster/src/raster-tileset/raster-tile-traversal.ts +++ b/packages/deck.gl-raster/src/raster-tileset/raster-tile-traversal.ts @@ -524,7 +524,7 @@ export class RasterTileNode { boundingVolume: OrientedBoundingBox; commonSpaceBounds: Bounds; } { - const [minZ, maxZ] = zRange; + const [minZ] = zRange; const tileCorners = this.level.projectedTileCorners(this.x, this.y); @@ -539,24 +539,38 @@ export class RasterTileNode { rescaleEPSG3857ToCommonSpace(xy), ); - const refPointPositions: [number, number, number][] = []; - for (const p of commonSpacePositions) { - refPointPositions.push([p[0], p[1], minZ]); + // Filter out NaN positions that arise when tile corners lie outside the + // CRS domain (e.g. corners of a bounding rectangle that extend beyond + // a Mollweide ellipse). Use only the valid subset for the OBB and AABB. + const validPositions = commonSpacePositions.filter( + ([x, y]) => Number.isFinite(x) && Number.isFinite(y), + ); - if (minZ !== maxZ) { - // Also sample at maximum elevation to capture the full 3D volume - refPointPositions.push([p[0], p[1], maxZ]); - } + // If no reference point projects successfully the tile lies entirely + // outside the CRS valid domain (e.g. a corner tile whose entire extent + // is outside the Mollweide ellipse). Return a degenerate bounding volume + // placed far off-screen so both the AABB bounds check and frustum + // culling reject it — never selecting it for rendering. + if (validPositions.length === 0) { + const OFF = -1e8; + return { + boundingVolume: makeOrientedBoundingBoxFromPoints([ + [OFF, OFF, minZ], + [OFF + 100, OFF, minZ], + [OFF, OFF + 100, minZ], + [OFF + 100, OFF + 100, minZ], + ]), + commonSpaceBounds: [OFF, OFF, OFF + 100, OFF + 100], + }; } - // Compute [minx, miny, maxx, maxy] in common space for quick bounds check - // TODO: this doesn't densify edges + // Compute [minx, miny, maxx, maxy] in common space for the bounds check let minX = Number.POSITIVE_INFINITY; let minY = Number.POSITIVE_INFINITY; let maxX = Number.NEGATIVE_INFINITY; let maxY = Number.NEGATIVE_INFINITY; - for (const [x, y] of commonSpacePositions) { + for (const [x, y] of validPositions) { if (x < minX) { minX = x; } @@ -572,8 +586,23 @@ export class RasterTileNode { } const commonSpaceBounds: Bounds = [minX, minY, maxX, maxY]; + + // Build the OBB from AABB corners with a minimum 1-unit extent in each + // axis. makeOrientedBoundingBoxFromPoints produces a degenerate OBB when + // the input points are collinear — which happens for projections like + // Mollweide where the only valid reference points all fall on the same + // latitude (the equatorial band). A degenerate OBB causes + // computeVisibility to return -1 (outside frustum), incorrectly culling + // the tile. The AABB insideBounds check handles precise culling. + const safeMaxX = Math.max(maxX, minX + 1); + const safeMaxY = Math.max(maxY, minY + 1); return { - boundingVolume: makeOrientedBoundingBoxFromPoints(refPointPositions), + boundingVolume: makeOrientedBoundingBoxFromPoints([ + [minX, minY, minZ], + [safeMaxX, minY, minZ], + [minX, safeMaxY, minZ], + [safeMaxX, safeMaxY, minZ], + ]), commonSpaceBounds, }; } diff --git a/packages/deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts b/packages/deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts index f99cb7a8..a5ca904f 100644 --- a/packages/deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts +++ b/packages/deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts @@ -180,8 +180,18 @@ export class RasterTileset2D extends Tileset2D { // the tile to keep `RasterLayer`'s reprojection-equality check stable // across renders (deck.gl recreates the layer instance every render, so // per-render-derived closures would regenerate the mesh every frame). - this.projectPosition = (x, y) => - rescaleEPSG3857ToCommonSpace(descriptor.projectTo3857(x, y)); + this.projectPosition = (x, y) => { + const result = rescaleEPSG3857ToCommonSpace( + descriptor.projectTo3857(x, y), + ); + if (Number.isFinite(result[0]) && Number.isFinite(result[1])) { + return result; + } + // Point is outside the CRS domain (e.g. beyond the Mollweide ellipse). + // Return NaN so the GPU rasterizer discards any triangle that touches + // this vertex, producing a clean domain-boundary cutoff with no artifacts. + return [NaN, NaN]; + }; this.unprojectPosition = (cx, cy) => { const [mx, my] = rescaleCommonSpaceToEPSG3857([cx, cy]); return descriptor.projectFrom3857(mx, my); diff --git a/packages/deck.gl-raster/tests/raster-tileset/affine-tileset-level.test.ts b/packages/deck.gl-raster/tests/raster-tileset/affine-tileset-level.test.ts index 8125e5b7..ad0dccdd 100644 --- a/packages/deck.gl-raster/tests/raster-tileset/affine-tileset-level.test.ts +++ b/packages/deck.gl-raster/tests/raster-tileset/affine-tileset-level.test.ts @@ -86,6 +86,37 @@ describe("AffineTilesetLevel", () => { expect(corners.bottomRight).toEqual([140, 160]); }); + it("clips the right edge of the last-column tile to arrayWidth", () => { + // arrayWidth=6 means tile 1 (pixels 4..7) is clipped to pixel 6. + const level = new AffineTilesetLevel({ + affine: SQUARE_AFFINE, + arrayWidth: 6, + arrayHeight: 8, + tileWidth: 4, + tileHeight: 4, + mpu: 1, + }); + // Interior tile (col=0): right edge is at pixel 4, unclipped. + expect(level.projectedTileCorners(0, 0).topRight).toEqual([140, 200]); + // Last-column tile (col=1): right edge clipped from pixel 8 to pixel 6. + expect(level.projectedTileCorners(1, 0).topRight).toEqual([160, 200]); + }); + + it("clips the bottom edge of the last-row tile to arrayHeight", () => { + const level = new AffineTilesetLevel({ + affine: SQUARE_AFFINE, + arrayWidth: 8, + arrayHeight: 6, + tileWidth: 4, + tileHeight: 4, + mpu: 1, + }); + // Interior tile (row=0): bottom at pixel 4, unclipped. + expect(level.projectedTileCorners(0, 0).bottomLeft).toEqual([100, 160]); + // Last-row tile (row=1): bottom clipped from pixel 8 to pixel 6. + expect(level.projectedTileCorners(0, 1).bottomLeft).toEqual([100, 140]); + }); + it("returns rotated quadrilateral corners for a rotated affine", () => { const level = new AffineTilesetLevel({ affine: ROTATED_AFFINE,