Skip to content

Commit e356c9c

Browse files
spepsgithub-actions[bot]
speps
authored and
github-actions[bot]
committed
Improve terrain raycast accuracy, queried X/Y position should closely match the reprojected lon/lat result now
GitOrigin-RevId: d6eb63d146d60ad0740eb40b80644aa3ab011531
1 parent 4a15e60 commit e356c9c

17 files changed

Lines changed: 324 additions & 71 deletions

File tree

‎debug/terrain-raycast.html‎

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
<!DOCTYPE html>
2+
<html>
3+
<head>
4+
<title>Terrain raycast debug</title>
5+
<meta charset="utf-8">
6+
<meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=no">
7+
<link rel="stylesheet" href="../dist/mapbox-gl.css" />
8+
<style>
9+
body { margin: 0; padding: 0; }
10+
html, body, #map { height: 100%; }
11+
#info {
12+
position: absolute;
13+
top: 10px;
14+
left: 10px;
15+
background: rgba(0,0,0,0.6);
16+
color: #fff;
17+
padding: 8px 12px;
18+
font: 12px/1.5 monospace;
19+
border-radius: 4px;
20+
pointer-events: none;
21+
}
22+
</style>
23+
</head>
24+
<body>
25+
<div id="map"></div>
26+
<div id="info">Click to sample terrain raycast</div>
27+
28+
<script src="../dist/mapbox-gl-dev.js"></script>
29+
<script type="module">
30+
import {getAccessToken} from './access_token_generated.js';
31+
32+
const map = window.map = new mapboxgl.Map({
33+
accessToken: getAccessToken(),
34+
container: 'map',
35+
zoom: 13.5,
36+
center: [-122.45814, 37.76159],
37+
style: 'mapbox://styles/mapbox/streets-v11',
38+
hash: true,
39+
projection: 'globe',
40+
});
41+
42+
map.addControl(new mapboxgl.NavigationControl());
43+
map.addControl(new mapboxgl.ScaleControl());
44+
45+
map.on('style.load', () => {
46+
map.addSource('mapbox-dem', {
47+
type: 'raster-dem',
48+
url: 'mapbox://mapbox.mapbox-terrain-dem-v1',
49+
tileSize: 514,
50+
maxzoom: 14
51+
});
52+
map.setTerrain({source: 'mapbox-dem', exaggeration: 2.0});
53+
map.setFog({});
54+
55+
map.addSource('samples', {
56+
type: 'geojson',
57+
data: {type: 'FeatureCollection', features: []}
58+
});
59+
map.addLayer({
60+
id: 'samples',
61+
type: 'circle',
62+
source: 'samples',
63+
paint: {
64+
'circle-radius': 4,
65+
'circle-color': '#ff00ff',
66+
'circle-pitch-alignment': 'viewport',
67+
'circle-pitch-scale': 'viewport'
68+
}
69+
});
70+
});
71+
72+
// Sample a 21×21 grid of points (±100px, step 10px) centred on (px, py).
73+
// Returns timing in ms.
74+
function sample(px, py) {
75+
const features = [];
76+
const t0 = performance.now();
77+
for (let dy = -100; dy <= 100; dy += 10) {
78+
for (let dx = -100; dx <= 100; dx += 10) {
79+
const ll = map.unproject(new mapboxgl.Point(px + dx, py + dy));
80+
features.push({
81+
type: 'Feature',
82+
geometry: {type: 'Point', coordinates: [ll.lng, ll.lat]}
83+
});
84+
}
85+
}
86+
const elapsed = performance.now() - t0;
87+
map.getSource('samples').setData({type: 'FeatureCollection', features});
88+
return elapsed;
89+
}
90+
91+
map.on('load', () => {
92+
const c = map.getCanvas();
93+
const ms = sample(c.width / window.devicePixelRatio / 2, c.height / window.devicePixelRatio / 2);
94+
document.getElementById('info').textContent = `Loaded — ${features_count()} samples in ${ms.toFixed(1)} ms`;
95+
});
96+
97+
map.on('click', (e) => {
98+
const ms = sample(e.point.x, e.point.y);
99+
document.getElementById('info').textContent =
100+
`Click (${e.point.x.toFixed(0)}, ${e.point.y.toFixed(0)}) — ${features_count()} samples in ${ms.toFixed(1)} ms`;
101+
console.log(`raycast: ${ms.toFixed(2)} ms for ${features_count()} points`);
102+
});
103+
104+
function features_count() { return 21 * 21; }
105+
</script>
106+
</body>
107+
</html>

‎src/data/dem_data.ts‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -99,7 +99,6 @@ export default class DEMData {
9999
values[this._idx(dim, dim)] = values[this._idx(dim - 1, dim - 1)];
100100
}
101101

102-
// Convert to float
103102
const unpack = sourceEncoding === "terrarium" ? unpackTerrarium : unpackMapbox;
104103
for (let i = 0; i < values.length; ++i) {
105104
const byteIdx = i * 4;

‎src/data/dem_tree.ts‎

Lines changed: 122 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -121,10 +121,6 @@ function triangleRayIntersect(
121121
return (acX * qvecX + acY * qvecY + acZ * qvecZ) * invDet;
122122
}
123123

124-
function frac(v: number, lo: number, hi: number) {
125-
return (v - lo) / (hi - lo);
126-
}
127-
128124
function decodeBounds(x: number, y: number, depth: number, boundsMinx: number, boundsMiny: number, boundsMaxx: number, boundsMaxy: number, outMin: Array<number>, outMax: Array<number>) {
129125
const scale = 1 << depth;
130126
const rangex = boundsMaxx - boundsMinx;
@@ -238,57 +234,108 @@ export default class DemMinMaxQuadTree {
238234
const {idx, t, nodex, nodey, depth} = stack.pop();
239235

240236
if (this.leaves[idx]) {
241-
// Create 2 triangles to approximate the surface plane for more precise tests
237+
// DDA over the leaf's DEM texels using the '/' diagonal (matching the terrain shader).
238+
// Accurate because buildDemMipmap mip-0 now stores true texel max, so tEnter = t
239+
// is the correct AABB entry point. Uses edge-sharing and a height-based skip.
242240
decodeBounds(nodex, nodey, depth, rootMinx, rootMiny, rootMaxx, rootMaxy, boundsMin, boundsMax);
243241

244242
const scale = 1 << depth;
245-
const minxUv = (nodex + 0) / scale;
243+
const minxUv = nodex / scale;
246244
const maxxUv = (nodex + 1) / scale;
247-
const minyUv = (nodey + 0) / scale;
245+
const minyUv = nodey / scale;
248246
const maxyUv = (nodey + 1) / scale;
249247

250-
// 4 corner points A, B, C and D defines the (quad) area covered by this node
251-
const az = sampleElevation(minxUv, minyUv, this.dem) * exaggeration;
252-
const bz = sampleElevation(maxxUv, minyUv, this.dem) * exaggeration;
253-
const cz = sampleElevation(maxxUv, maxyUv, this.dem) * exaggeration;
254-
const dz = sampleElevation(minxUv, maxyUv, this.dem) * exaggeration;
255-
256-
const t0 = triangleRayIntersect(
257-
258-
boundsMin[0], boundsMin[1], az, // A
259-
260-
boundsMax[0], boundsMin[1], bz, // B
261-
262-
boundsMax[0], boundsMax[1], cz, // C
263-
p, d);
264-
265-
const t1 = triangleRayIntersect(
266-
267-
boundsMax[0], boundsMax[1], cz,
268-
269-
boundsMin[0], boundsMax[1], dz,
270-
271-
boundsMin[0], boundsMin[1], az,
272-
p, d);
248+
const {dim: demSz, floatView: fv, stride: st} = this.dem;
249+
const txMin = Math.floor(minxUv * demSz);
250+
const txMax = Math.ceil(maxxUv * demSz);
251+
const tyMin = Math.floor(minyUv * demSz);
252+
const tyMax = Math.ceil(maxyUv * demSz);
253+
const txCount = txMax - txMin;
254+
const tyCount = tyMax - tyMin;
255+
if (txCount <= 0 || tyCount <= 0) continue;
256+
257+
const x0 = boundsMin[0];
258+
const y0 = boundsMin[1];
259+
const cellW = (boundsMax[0] - x0) / txCount;
260+
const cellH = (boundsMax[1] - y0) / tyCount;
261+
262+
// Starting point on the ray (clamped to tile entry)
263+
const tEnter = Math.max(0, t);
264+
const entryX = p[0] + d[0] * tEnter;
265+
const entryY = p[1] + d[1] * tEnter;
266+
267+
// Step direction
268+
const hasX = Math.abs(d[0]) > 1e-10;
269+
const hasY = Math.abs(d[1]) > 1e-10;
270+
const stepI = d[0] >= 0 ? 1 : -1;
271+
const stepJ = d[1] >= 0 ? 1 : -1;
272+
273+
let ci = Math.max(0, Math.min(txCount - 1, Math.floor((entryX - x0) / cellW)));
274+
let cj = Math.max(0, Math.min(tyCount - 1, Math.floor((entryY - y0) / cellH)));
275+
276+
// Distance along ray to cross one cell
277+
const tDeltaI = hasX ? Math.abs(cellW / d[0]) : Number.MAX_VALUE;
278+
const tDeltaJ = hasY ? Math.abs(cellH / d[1]) : Number.MAX_VALUE;
279+
// Distance to next cell boundary
280+
let tNextI = hasX ? tEnter + (x0 + (d[0] >= 0 ? ci + 1 : ci) * cellW - entryX) / d[0] : Number.MAX_VALUE;
281+
let tNextJ = hasY ? tEnter + (y0 + (d[1] >= 0 ? cj + 1 : cj) * cellH - entryY) / d[1] : Number.MAX_VALUE;
282+
283+
const elev = (tx: number, ty: number) => fv[(ty + 1) * st + tx + 1] * exaggeration;
284+
let e00 = elev(txMin + ci, tyMin + cj);
285+
let e10 = elev(txMin + ci + 1, tyMin + cj);
286+
let e01 = elev(txMin + ci, tyMin + cj + 1);
287+
let e11 = elev(txMin + ci + 1, tyMin + cj + 1);
288+
289+
let closestT: number | null = null;
290+
291+
// Traverse cells in order along the ray
292+
for (let step = 0; step < txCount + tyCount; step++) {
293+
const tExit = Math.min(tNextI, tNextJ);
294+
295+
// Height-based skip: for a downward ray, tExit is the lowest z the ray reaches in this cell
296+
if (d[2] >= 0 || p[2] + d[2] * tExit <= Math.max(e00, e10, e01, e11)) {
297+
const cx0 = x0 + ci * cellW;
298+
const cx1 = cx0 + cellW;
299+
const cy0 = y0 + cj * cellH;
300+
const cy1 = cy0 + cellH;
301+
// Triangle 1: (1,0) -> (0,0) -> (0,1)
302+
const t0 = triangleRayIntersect(cx1, cy0, e10, cx0, cy0, e00, cx0, cy1, e01, p, d);
303+
// Triangle 2: (0,1) -> (1,1) -> (1,0)
304+
const t1 = triangleRayIntersect(cx0, cy1, e01, cx1, cy1, e11, cx1, cy0, e10, p, d);
305+
if (t0 != null && t0 >= 0 && (closestT === null || t0 < closestT)) closestT = t0;
306+
if (t1 != null && t1 >= 0 && (closestT === null || t1 < closestT)) closestT = t1;
307+
}
273308

274-
const tMin = Math.min(
275-
t0 !== null ? t0 : Number.MAX_VALUE,
276-
t1 !== null ? t1 : Number.MAX_VALUE);
309+
// Move to next cell
310+
const steppedI = tNextI < tNextJ;
311+
if (steppedI) { ci += stepI; tNextI += tDeltaI; } else { cj += stepJ; tNextJ += tDeltaJ; }
277312

278-
// The ray might go below the two surface triangles but hit one of the sides.
279-
// This covers the case of skirt geometry between two dem tiles of different zoom level
280-
if (tMin === Number.MAX_VALUE) {
281-
const hitPos = vec3.scaleAndAdd([], p, d, t);
313+
// Check if we've exited the grid
314+
if (ci < 0 || ci >= txCount || cj < 0 || cj >= tyCount) break;
315+
// Early termination: found a hit and the ray has passed it in XY
316+
if (closestT !== null && tExit > closestT) break;
282317

283-
const fracx = frac(hitPos[0], boundsMin[0], boundsMax[0]);
318+
// Update shared edges (reuse 2 corners, fetch 2 new ones)
319+
if (steppedI) {
320+
if (stepI > 0) { e00 = e10; e01 = e11; e10 = elev(txMin + ci + 1, tyMin + cj); e11 = elev(txMin + ci + 1, tyMin + cj + 1); } else { e10 = e00; e11 = e01; e00 = elev(txMin + ci, tyMin + cj); e01 = elev(txMin + ci, tyMin + cj + 1); }
321+
} else {
322+
if (stepJ > 0) { e00 = e01; e10 = e11; e01 = elev(txMin + ci, tyMin + cj + 1); e11 = elev(txMin + ci + 1, tyMin + cj + 1); } else { e01 = e00; e11 = e10; e00 = elev(txMin + ci, tyMin + cj); e10 = elev(txMin + ci + 1, tyMin + cj); }
323+
}
324+
}
284325

285-
const fracy = frac(hitPos[1], boundsMin[1], boundsMax[1]);
326+
if (closestT !== null) return closestT || 0; // normalize -0 to +0
286327

287-
if (bilinearLerp(az, bz, dz, cz, fracx, fracy) >= hitPos[2])
288-
return t;
289-
} else {
290-
return tMin;
291-
}
328+
// Skirt fallback: handles rays that pass below the triangulated surface
329+
// (e.g. horizontal rays through the gap between tiles at different zoom levels).
330+
const hitPos = vec3.scaleAndAdd([], p, d, t);
331+
const fracx = (hitPos[0] - boundsMin[0]) / (boundsMax[0] - boundsMin[0]);
332+
const fracy = (hitPos[1] - boundsMin[1]) / (boundsMax[1] - boundsMin[1]);
333+
const az = sampleElevation(minxUv, minyUv, this.dem) * exaggeration;
334+
const bz = sampleElevation(maxxUv, minyUv, this.dem) * exaggeration;
335+
const cz = sampleElevation(maxxUv, maxyUv, this.dem) * exaggeration;
336+
const dz = sampleElevation(minxUv, maxyUv, this.dem) * exaggeration;
337+
if (bilinearLerp(az, bz, dz, cz, fracx, fracy) >= hitPos[2])
338+
return t;
292339

293340
continue;
294341
}
@@ -429,7 +476,6 @@ export function buildDemMipmap(dem: DEMData): Array<MipLevel> {
429476
const mips: Array<MipLevel> = [];
430477

431478
let blockCount = Math.ceil(Math.pow(2, levelCount));
432-
const blockSize = 1 / blockCount;
433479

434480
const blockSamples = (x: number, y: number, size: number, exclusive: boolean, outBounds: Array<number>) => {
435481
const padding = exclusive ? 1 : 0;
@@ -444,28 +490,48 @@ export function buildDemMipmap(dem: DEMData): Array<MipLevel> {
444490
outBounds[3] = maxy;
445491
};
446492

447-
// The first mip (0) is built by sampling 4 corner points of each 8x8 texel block
493+
// Mip 0: use bilinear corner samples for MIN (preserves existing behaviour for
494+
// getMinElevationBelowMSL and the rendering pipeline) but scan all raw DEM texels
495+
// for MAX. The corner approach underestimates peaks between samples, so the AABB
496+
// z-max was too low — causing the DDA to start at the wrong cell and miss hits.
448497
let mip = new MipLevel(blockCount);
498+
const blockSize = 1 / blockCount;
499+
const {floatView, stride} = dem;
449500
const blockBounds = [];
450501

451502
for (let idx = 0; idx < blockCount * blockCount; idx++) {
452-
const y = Math.floor(idx / blockCount);
453-
const x = idx % blockCount;
503+
const by = Math.floor(idx / blockCount);
504+
const bx = idx % blockCount;
454505

506+
// MIN: bilinear at the 4 corners (unchanged from original)
455507
// eslint-disable-next-line @typescript-eslint/no-unsafe-argument
456-
blockSamples(x, y, blockSize, false, blockBounds);
457-
508+
blockSamples(bx, by, blockSize, false, blockBounds);
458509
// eslint-disable-next-line @typescript-eslint/no-unsafe-argument
459-
const e0 = sampleElevation(blockBounds[0], blockBounds[1], dem); // minx, miny
510+
const e0 = sampleElevation(blockBounds[0], blockBounds[1], dem);
460511
// eslint-disable-next-line @typescript-eslint/no-unsafe-argument
461-
const e1 = sampleElevation(blockBounds[2], blockBounds[1], dem); // maxx, miny
512+
const e1 = sampleElevation(blockBounds[2], blockBounds[1], dem);
462513
// eslint-disable-next-line @typescript-eslint/no-unsafe-argument
463-
const e2 = sampleElevation(blockBounds[2], blockBounds[3], dem); // maxx, maxy
514+
const e2 = sampleElevation(blockBounds[2], blockBounds[3], dem);
464515
// eslint-disable-next-line @typescript-eslint/no-unsafe-argument
465-
const e3 = sampleElevation(blockBounds[0], blockBounds[3], dem); // minx, maxy
516+
const e3 = sampleElevation(blockBounds[0], blockBounds[3], dem);
517+
const blockMin = Math.min(e0, e1, e2, e3);
518+
519+
// MAX: scan all raw texels (plus one border texel per edge) for true peak detection
520+
const txStart = Math.max(0, bx * texelSizeOfMip0 - 1);
521+
const txEnd = Math.min(demSize - 1, (bx + 1) * texelSizeOfMip0);
522+
const tyStart = Math.max(0, by * texelSizeOfMip0 - 1);
523+
const tyEnd = Math.min(demSize - 1, (by + 1) * texelSizeOfMip0);
524+
let blockMax = -Infinity;
525+
for (let ty = tyStart; ty <= tyEnd; ty++) {
526+
const rowBase = (ty + 1) * stride;
527+
for (let tx = txStart; tx <= txEnd; tx++) {
528+
const v = floatView[rowBase + tx + 1];
529+
if (v > blockMax) blockMax = v;
530+
}
531+
}
466532

467-
mip.minimums.push(Math.min(e0, e1, e2, e3));
468-
mip.maximums.push(Math.max(e0, e1, e2, e3));
533+
mip.minimums.push(blockMin);
534+
mip.maximums.push(blockMax);
469535
mip.leaves.push(1);
470536
}
471537

‎test/integration/render-tests/fit-screen-coordinates/terrain/expected.png‎

69.7 KB
Loading

‎test/integration/render-tests/fog/terrain/equal-range/style.json‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,8 @@
33
"metadata": {
44
"test": {
55
"height": 256,
6-
"width": 256
6+
"width": 256,
7+
"allowed": 0.0002
78
}
89
},
910
"center": [-113.32296, 35.94662],

‎test/integration/render-tests/free-camera/terrain/expected.png‎

85.3 KB
Loading

‎test/integration/render-tests/model-layer/model-padded-terrain/expected.png‎

25.9 KB
Loading

‎test/integration/render-tests/model-layer/model-padded-terrain/style.json‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -62,7 +62,7 @@
6262
"source": "mapbox",
6363
"source-layer": "road",
6464
"paint": {
65-
"line-color": "lightyellow",
65+
"line-color": "black",
6666
"line-width": 10,
6767
"line-opacity": 0.3
6868
}

‎test/integration/render-tests/terrain/globe-high-exaggeration/expected.png‎

-7.31 KB
Loading

‎test/integration/render-tests/terrain/raycast/depth-buffer-0/expected.png‎

29 Bytes
Loading

0 commit comments

Comments
 (0)