30:T9f5,md`# MARTINI: Real-Time RTIN Terrain Mesh ## Overview of Fork This is a fork of [MARTINI: Real-Time RTIN Terrain Mesh](https://observablehq.com/@mourner/martin-real-time-rtin-terrain-mesh). You may want to go read the explanations there, since I haven't adapted them to my changes. Instead, I just added the following functionality: - Mesh texture from satellite imagery - Terrain / image tiles > 256x256 pixels, using the [static-map library](https://github.com/rkaravia/static-map) - Animated GIF export - Video (.webm) export, see [this notebook](https://observablehq.com/@rzoller/canvas-recorder) for more information More ideas, maybe for someone else: - Choose location on a map - Place the mesh on a map, similar to [this example](https://docs.mapbox.com/mapbox-gl-js/example/add-3d-model/) - Export to an [OBJ file](https://en.wikipedia.org/wiki/Wavefront_.obj_file), or some other 3D mesh format ## Original text below To **render terrain in 3D**, you need a **triangular mesh** — a set of 3D triangles representing the surface. The easiest way to get one from height data is a uniform grid (two triangles per pixel), but that's _a lot_ of triangles. For fast rendering, we want to reduce this amount by adapting the mesh to the surface — using bigger triangles in relatively smooth areas, and smaller ones where there's a lot of detail. One such approach, [proposed by Will Evans et al. in 1997](https://www.cs.ubc.ca/~will/papers/rtin.pdf), is called **Right-Triangulated Irregular Networks**, or **RTIN**. An RTIN mesh consists of only right-angle triangles, which makes it less precise than Delaunay-based TIN meshes, requiring more triangles to approximate the same surface. But RTIN has two significant advantages: 1. The algorithm generates a **hierarchy of all approximations** of varying precisions — after running it once, you can quickly retrieve a mesh for any given level of detail. 2. It's very fast, making it viable for **client-side meshing** from raster terrain tiles. Surprisingly, I haven't found any prior attempts to do it in the browser. So let me present you a client-side terrain mesh generation code I called **MARTINI**: [**M**apbox](https://www.mapbox.com)'s **A**wesome **R**ight-**T**riangulated **I**rregular **N**etworks **I**mproved! It will live as an open source library at [github.com/mapbox/martini](https://github.com/mapbox/martini). Drag and zoom the terrain visualization below, and drag the slider to adjust mesh precision in real time. Then we'll dive into how this works. `31:T6fe,errors = { const errors = new Float32Array(gridSize * gridSize); const numSmallestTriangles = tileSize * tileSize; const numTriangles = numSmallestTriangles * 2 - 2; // 2 + 4 + 8 + ... 2^k = 2 * 2^k - 2 const lastLevelIndex = numTriangles - numSmallestTriangles; // iterate over all possible triangles, starting from the smallest level for (let i = numTriangles - 1; i >= 0; i--) { // get triangle coordinates from its index in an implicit binary tree let id = i + 2; let ax = 0, ay = 0, bx = 0, by = 0, cx = 0, cy = 0; if (id & 1) { bx = by = cx = tileSize; // bottom-left triangle } else { ax = ay = cy = tileSize; // top-right triangle } while ((id >>= 1) > 1) { const mx = (ax + bx) >> 1; const my = (ay + by) >> 1; if (id & 1) { // left half bx = ax; by = ay; ax = cx; ay = cy; } else { // right half ax = bx; ay = by; bx = cx; by = cy; } cx = mx; cy = my; } // calculate error in the middle of the long edge of the triangle const interpolatedHeight = (terrain[ay * gridSize + ax] + terrain[by * gridSize + bx]) / 2; const middleIndex = ((ay + by) >> 1) * gridSize + ((ax + bx) >> 1); const middleError = Math.abs(interpolatedHeight - terrain[middleIndex]); if (i >= lastLevelIndex) { // smallest triangles errors[middleIndex] = middleError; } else { // bigger triangles; accumulate error with children const leftChildError = errors[((ay + cy) >> 1) * gridSize + ((ax + cx) >> 1)]; const rightChildError = errors[((by + cy) >> 1) * gridSize + ((bx + cx) >> 1)]; errors[middleIndex] = Math.max(errors[middleIndex], middleError, leftChildError, rightChildError); } } return errors; }32:T467,{ replay; await visibility(); const size = Math.min(width, 512); const ratio = size / tileSize; const ctx = DOM.context2d(size, size, 2); ctx.fillStyle = 'white'; ctx.fillRect(0, 0, size, size); ctx.globalAlpha = 0.5; ctx.drawImage(terrainViz, 0, 0, size, size); ctx.globalAlpha = 1; function drawLine(ax, ay, bx, by, color) { ctx.beginPath(); ctx.moveTo(ax * ratio, ay * ratio); ctx.lineTo(bx * ratio, by * ratio); ctx.strokeStyle = color; ctx.stroke(); } function* split(ax, ay, bx, by, cx, cy) { const mx = (ax + bx) >> 1; const my = (ay + by) >> 1; if (Math.abs(ax - cx) + Math.abs(ay - cy) > 1 && errors[my * (tileSize + 1) + mx] > 100) { drawLine(mx, my, cx, cy, 'red'); yield ctx.canvas; yield* split(cx, cy, ax, ay, mx, my); yield* split(bx, by, cx, cy, mx, my); drawLine(mx, my, cx, cy, 'black'); } } drawLine(0, 0, tileSize, tileSize, 'red'); yield* split(0, 0, tileSize, tileSize, tileSize, 0, 0); yield* split(tileSize, tileSize, 0, 0, 0, tileSize, 0); drawLine(0, 0, tileSize, tileSize, 'black'); }33:T4a1,md`Note that **for this demo's purposes**, we allocate maximum possible vertex/index buffers at the beginning, and only update the index array and draw range on precision change; in real world apps, the approach would be different since we wouldn't need to update the mesh of a tile as often, and would instead only upload the vertices and indices used in the mesh, reducing the load on the GPU. This is reflected in the [actual library](https://github.com/mapbox/martin), where mesh retrieval routine generates an array of vertices that are only used in that mesh (avoiding duplicate vertices), and triangles referencing it. ## Thanks! Thank you for reading! Feel free to reach out to me [on Twitter](https://twitter.com/mourner) with comments & questions, check out [the repository](https://github.com/mapbox/martin) (very much a work in progress), and stay tuned for more cool stuff while we explore the 3D space! P.S. Originally I named this project MARTIN, completely forgetting there's an [existing awesome Mapbox-related open source project](https://martin.urbica.co) with the same name! After embarassing myself, I had to hastily rename this to MARTINI to save face. 😅 `34:T4fa,geometry = { const geometry = new THREE.BufferGeometry(); const vertices = new Float32Array(gridSize * gridSize * 3); const uvs = new Float32Array(gridSize * gridSize * 2); const indices = new Uint32Array(tileSize * tileSize * 6); let index = 0; let minTerrain = Infinity; for (let i = 0; i < terrain.length; i++) { minTerrain = Math.min(terrain[i], minTerrain); } for (let y = 0; y <= tileSize; y++) { for (let x = 0; x <= tileSize; x++) { const i = y * gridSize + x; vertices[3 * i + 0] = x / tileSize - 0.5; vertices[3 * i + 1] = 0.5 - y / tileSize; vertices[3 * i + 2] = (terrain[i] - minTerrain) / metersPerPixel / tileSize * terrainExaggeration; uvs[2 * i + 0] = x / tileSize; uvs[2 * i + 1] = 1 - y / tileSize; indices[index++] = i + 1; indices[index++] = i; indices[index++] = i + gridSize; indices[index++] = i + 1; indices[index++] = i + gridSize; indices[index++] = i + gridSize + 1; } } geometry.addAttribute('position', new THREE.BufferAttribute(vertices, 3)); geometry.addAttribute('uv', new THREE.BufferAttribute(uvs, 2)); geometry.setIndex(new THREE.BufferAttribute(indices, 1)); geometry.computeVertexNormals(); return geometry; }35:[] 2e:{"title":"MARTINI: Real-Time RTIN Terrain Mesh","cells":[{"id":0,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":1433,"value":"viewof location = html``;","pinned":false,"hidden":false,"mode":"ojs"},{"id":1559,"value":"viewof tilesPerAxisIdx = html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":1562,"value":"md`Size: ${tilesPerAxis}x${tilesPerAxis} tiles`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1487,"value":"viewof showTexture = {\n const div = html`Show texture`;\n const input = div.querySelector(\"input\");\n div.value = input.checked = true;\n input.oninput = () => div.value = input.checked;\n return div;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1588,"value":"DOM.download(() => {\n const gif = new GifCompiler();\n const fps = 24;\n const delay = Math.floor(1000 / fps);\n const steps = 240;\n\n renderer.setPixelRatio(1);\n for(let t = 0; t < steps; t += 1) {\n scene.rotation.z = (t / steps) * (2 * Math.PI);\n renderer.render(scene, camera);\n gif.addFrame(renderer.domElement, delay);\n }\n renderer.setPixelRatio(devicePixelRatio);\n\n return gif.render();\n}, \"mesh.gif\", \"Render animated GIF\");","pinned":false,"hidden":false,"mode":"ojs"},{"id":1715,"value":"DOM.download(async () => {\n const recorder = new CanvasRecorder(renderer.domElement);\n const steps = 360;\n \n renderer.setPixelRatio(1);\n for(let t = 0; t < steps; t += 1) {\n scene.rotation.z = (t / steps) * (2 * Math.PI);\n renderer.render(scene, camera);\n await recorder.captureFrame();\n }\n renderer.setPixelRatio(devicePixelRatio);\n\n return recorder.stop();\n}, 'martini.webm', 'Render video (.webm)');","pinned":false,"hidden":false,"mode":"ojs"},{"id":1724,"value":"md`\nNote: Adjust framerate and convert to mp4 with [ffmpeg](https://ffmpeg.org/):\n~~~\nffmpeg -i martini.webm -filter \"setpts=N/(30*TB)\" -r 30 martini.mp4\n~~~\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":864,"value":"viewof map = renderer.domElement","pinned":false,"hidden":false,"mode":"ojs"},{"id":310,"value":"viewof maxError = html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":556,"value":"{\n const numTriangles = updatedGeometry.drawRange.count / 3;\n const maxTriangles = tileSize * tileSize * 2;\n return md`Generated ${numTriangles} (${(100 * numTriangles / maxTriangles).toFixed(2)}%) triangles at ${maxError.toFixed(1)}m error from a ${tileSize}×${tileSize}px tile.`\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":220,"value":"{\n const size = Math.min(width, 640);\n const ratio = size / tileSize;\n const ctx = DOM.context2d(size, size, 2);\n \n ctx.fillStyle = 'white';\n ctx.fillRect(0, 0, size, size);\n\n ctx.globalAlpha = 0.4;\n ctx.drawImage(terrainViz, 0, 0, size, size);\n ctx.globalAlpha = 1;\n \n ctx.lineWidth = 0.5;\n ctx.beginPath();\n\n function split(ax, ay, bx, by, cx, cy) {\n const mx = (ax + bx) >> 1;\n const my = (ay + by) >> 1;\n\n if (Math.abs(ax - cx) + Math.abs(ay - cy) > 1 && errors[my * (tileSize + 1) + mx] > maxError) {\n ctx.moveTo(mx * ratio, my * ratio);\n ctx.lineTo(cx * ratio, cy * ratio);\n split(cx, cy, ax, ay, mx, my);\n split(bx, by, cx, cy, mx, my);\n }\n }\n \n ctx.moveTo(0, 0);\n ctx.lineTo(ratio * tileSize, ratio * tileSize);\n \n split(0, 0, tileSize, tileSize, tileSize, 0, 0);\n split(tileSize, tileSize, 0, 0, 0, tileSize, 0);\n \n ctx.stroke();\n \n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1073,"value":"md`For demo purposes, we render the 3D viz above in a way where the triangles would be clearly visible, but with proper lighting and shading, the terrain would look smooth. Also, it's exaggerated by a factor of ${terrainExaggeration} to make it more prominent.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":704,"value":"md`## Preparing Height Data\n\nLet's get a sample elevation data from a [Mapbox Terrain-RGB](https://docs.mapbox.com/help/troubleshooting/access-elevation-data/#mapbox-terrain-rgb) ${tileSize}×${tileSize} raster tile:\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1192,"value":"tileImg = getMapImage(`https://a.tiles.mapbox.com/v4/mapbox.terrain-rgb/{z}/{x}/{y}.png?access_token=${mapboxToken}`)","pinned":true,"hidden":false,"mode":"ojs"},{"id":1542,"value":"tilesPerAxis = 2 ** tilesPerAxisIdx","pinned":false,"hidden":false,"mode":"ojs"},{"id":1209,"value":"tileSize = 256 * tilesPerAxis","pinned":false,"hidden":false,"mode":"ojs"},{"id":1195,"value":"md`Then we'll decode it into a grid of height values in meters:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":28,"value":"terrain = {\n const ctx = DOM.context2d(tileSize, tileSize, 1);\n ctx.drawImage(tileImg, 0, 0);\n const data = ctx.getImageData(0, 0, tileSize, tileSize).data;\n\n const terrain = new Float32Array(gridSize * gridSize);\n \n // decode terrain values\n for (let y = 0; y < tileSize; y++) {\n for (let x = 0; x < tileSize; x++) {\n const k = (y * tileSize + x) * 4;\n const r = data[k + 0];\n const g = data[k + 1];\n const b = data[k + 2];\n terrain[y * gridSize + x] = (r * 256 * 256 + g * 256.0 + b) / 10.0 - 10000.0;\n }\n }\n // backfill right and bottom borders\n for (let x = 0; x < gridSize - 1; x++) {\n terrain[gridSize * (gridSize - 1) + x] = terrain[gridSize * (gridSize - 2) + x];\n }\n for (let y = 0; y < gridSize; y++) {\n terrain[gridSize * y + gridSize - 1] = terrain[gridSize * y + gridSize - 2];\n }\n return terrain;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":71,"value":"terrainViz = drawGrid(terrain, 0.6)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1124,"value":"md`The RTIN algorithm works on grids of the size _(2^k+1)_ × _(2^k+1)_ — that's why we add 1-pixel borders on the right and bottom. On interactive maps, this property would also help us stitch the tiles together for seamless terrain.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1031,"value":"gridSize = tileSize + 1","pinned":false,"hidden":false,"mode":"ojs"},{"id":709,"value":"md`## Generating RTIN Hierarchy\n\nThe RTIN algorithm generates an error map — a grid of error values that will guide subsequent mesh retrieval by indicating whether to split a certain triangle in the middle given its height error value in that position. It looks like this:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":206,"value":"{\n const canvas = drawGrid(errors, 0.5, 0.1);\n canvas.style.width = Math.min((canvas.width * 2), width, 514) + 'px';\n return canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":924,"value":"md`The idea is to start with the smallest triangles, calculating the approximation error for each, and then propagate that error to parent triangles (keeping track of maximum error), repeating the process until we reach the top two triangles and have a full error map.\n\nThe way we accumulate errors from child triangles makes the resulting meshes conformal — meaning there are no T-junctions, and consequently no gaps in the mesh.\n\nThe algorithm has an _O(numPixels)_ computational complexity and takes just **40 lines of code**, and runs in about **20ms** on my laptop:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":36,"value":"$31","pinned":true,"hidden":false,"mode":"ojs"},{"id":727,"value":"md`## Retrieving a Mesh\nTo retrieve a mesh from an RTIN hierarchy, we start with two big triangles and recursively subdivide them when they don't approximate the underlying surface well enough according the error map. This is essentially a depth-first search in an implicit binary tree of triangles, and takes _O(numTriangles)_ steps, resulting in nearly instant mesh generation. \n\nHere's a visualization of mesh retrieval for a 100 meters error:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1273,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1277,"value":"viewof replay = html``","pinned":false,"hidden":false,"mode":"ojs"},{"id":1286,"value":"md`And here's the code for updating the mesh of the actual 3D terrain visualization at the top of the notebook:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":777,"value":"updatedGeometry = {\n let i = 0;\n const indices = geometry.index.array;\n \n function processTriangle(ax, ay, bx, by, cx, cy) {\n // middle of the long edge\n const mx = (ax + bx) >> 1;\n const my = (ay + by) >> 1;\n\n if (Math.abs(ax - cx) + Math.abs(ay - cy) > 1 && errors[my * gridSize + mx] > maxError) {\n // triangle doesn't approximate the surface well enough; split it into two\n processTriangle(cx, cy, ax, ay, mx, my);\n processTriangle(bx, by, cx, cy, mx, my);\n\n } else {\n // add a triangle to the final mesh\n indices[i++] = ay * gridSize + ax;\n indices[i++] = by * gridSize + bx;\n indices[i++] = cy * gridSize + cx;\n }\n }\n\n processTriangle(0, 0, tileSize, tileSize, tileSize, 0);\n processTriangle(tileSize, tileSize, 0, 0, 0, tileSize);\n \n geometry.index.needsUpdate = true;\n geometry.setDrawRange(0, i);\n renderer.render(scene, camera);\n \n return geometry;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1052,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":883,"value":"md`## Three.js terrain setup`","pinned":false,"hidden":false,"mode":"ojs"},{"id":993,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":763,"value":"camera = {\n const camera = new THREE.PerspectiveCamera(45, Math.min(width, 640) / height, 0.0001, 1000);\n camera.position.set(0.5, 0, 0.5);\n camera.up.set(0, 0, 1);\n return camera;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":766,"value":"scene = {\n const scene = new THREE.Scene();\n scene.background = new THREE.Color(0xffffff);\n const material = new THREE.MeshNormalMaterial({flatShading: true, side: THREE.DoubleSide});\n const mesh = new THREE.Mesh(geometry, showTexture ? textureMaterial : material);\n mesh.name = 'terrain';\n scene.add(mesh);\n return scene;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":769,"value":"renderer = {\n const renderer = new THREE.WebGLRenderer({antialias: true});\n const controls = new THREE.OrbitControls(camera, renderer.domElement);\n invalidation.then(() => (controls.dispose(), renderer.dispose()));\n renderer.setSize(Math.min(width, 640), height);\n renderer.setPixelRatio(devicePixelRatio);\n controls.addEventListener(\"change\", () => renderer.render(scene, camera));\n renderer.render(scene, camera);\n return renderer;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":756,"value":"height = 420","pinned":false,"hidden":false,"mode":"ojs"},{"id":584,"value":"md`## Other stuff`","pinned":false,"hidden":false,"mode":"ojs"},{"id":957,"value":"terrainExaggeration = 1","pinned":false,"hidden":false,"mode":"ojs"},{"id":934,"value":"metersPerPixel = {\n const { lon, lat, zoom } = lonLatZoom;\n const numTiles = Math.pow(2, zoom);\n return 40075016.7 * Math.cos(lat * Math.PI / 180) / numTiles / 256;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":181,"value":"drawGrid = (data, lightness = 0.5, cutoff = 1.0) => {\n const size = Math.sqrt(data.length);\n const ctx = DOM.context2d(size, size, 1);\n ctx.canvas.style.imageRendering = 'pixelated';\n const imgData = ctx.getImageData(0, 0, size, size);\n const minHeight = data.reduce((a, b) => Math.min(a, b));\n const maxHeight = data.reduce((a, b) => Math.max(a, b));\n \n for (let y = 0; y < size; y++) {\n for (let x = 0; x < size; x++) {\n const k = y * size + x;\n const v = Math.min((data[k] - minHeight) / (maxHeight - minHeight), cutoff) / cutoff;\n const [r, g, b] = hsl2rgb(230 - v * 230, 1, lightness);\n imgData.data[4 * k + 0] = r * 255;\n imgData.data[4 * k + 1] = g * 255;\n imgData.data[4 * k + 2] = b * 255;\n imgData.data[4 * k + 3] = 255;\n }\n }\n ctx.putImageData(imgData, 0, 0);\n return ctx.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":113,"value":"hsl2rgb = {\n const f = (h, n) => { \n const k = (n + h / 30) % 12; \n return Math.max(Math.min(k - 3, 9 - k, 1), -1);\n };\n return (h, s, l) => {\n const a = s * Math.min(l, 1 - l);\n return [l - a * f(h, 0), l - a * f(h, 8), l - a * f(h, 4)];\n };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1545,"value":"lonLatZoom = {\n const [lon, lat, zoom] = location.split('/').map(parseFloat);\n return { lon, lat, zoom: zoom + (Math.log(tilesPerAxis) / Math.log(2)) };\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1491,"value":"textureImg = getMapImage(`https://api.mapbox.com/v4/mapbox.satellite/{z}/{x}/{y}.jpg?access_token=${mapboxToken}`)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1493,"value":"texture = new Promise(resolve => new THREE.TextureLoader().load(textureImg.src, resolve));","pinned":false,"hidden":false,"mode":"ojs"},{"id":1495,"value":"textureMaterial = new THREE.MeshBasicMaterial( { map: texture } );","pinned":false,"hidden":false,"mode":"ojs"},{"id":1508,"value":"function getMapImage(url) {\n return new Promise((resolve) => {\n const { lon, lat, zoom } = lonLatZoom;\n const canvas = document.createElement('canvas');\n canvas.width = canvas.height = tileSize;\n new StaticMap(url).getMap(canvas, lon, lat, zoom, () => {\n const image = new Image();\n image.onload = () => resolve(image);\n image.src = canvas.toDataURL();\n });\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1222,"value":"mapboxToken = 'pk.eyJ1IjoibW91cm5lciIsImEiOiJWWnRiWG1VIn0.j6eccFHpE3Q04XPLI7JxbA'","pinned":false,"hidden":false,"mode":"ojs"},{"id":1059,"value":"md`## Imports`","pinned":false,"hidden":false,"mode":"ojs"},{"id":753,"value":"THREE = {\n const THREE = window.THREE = await require(\"three@0.107.0/build/three.min.js\");\n await require(\"three@0.107.0/examples/js/controls/OrbitControls.js\").catch(() => {});\n return THREE;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":22,"value":"import {fetchImage} from '@mbostock/cross-origin-images'","pinned":false,"hidden":false,"mode":"ojs"},{"id":1501,"value":"StaticMap = require('@rz0/static-map@0.1.1/static-map.js')","pinned":false,"hidden":false,"mode":"ojs"},{"id":1585,"value":"import {GifCompiler} from '@rzoller/gifcompiler';","pinned":false,"hidden":false,"mode":"ojs"},{"id":1712,"value":"import { CanvasRecorder } from '@rzoller/canvas-recorder'","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}