31:T490,This notebook builds up on [Vladimir Agafonkin’s implementation](https://observablehq.com/@mourner/fast-distance-transform) of the signed distance function (SDF) algorithm by Pedro Felzenszwalb and Daniel Huttenlocher. Color map: Oléron, by [Fabio Crameri](https://observablehq.com/@fil/colormaps). The challenge in this case was to extend the algorithm to work on the sphere, using great-circle distances instead of planar/euclidian distances. To do this, I needed some additions to the original function: 1. ensure continuity at the poles and at the antimeridian. 2. take into account the varying length of a parallel with the latitude. 3. also, it had to be *signed* (with “negative” distances in land). Uncheck the box to see what happens if you compute that globe’s distances in planar coordinates. Note that in the spherical case, the result is not mathematically perfect — if you have just a tiny island and the whole sphere to cover with SDF, the result is clearly wrong. But it seems to be a good approximation. Border detection technique updated Feb. 2024, with thanks to [Harry Stevens](/@harrystevens) for noticing the flaw in the original.32:T520,// spherical version sphedt1d = (l) => { const l3 = l * 3; const v = new Uint16Array(l3); const f = new Float64Array(l3); const z = new Float64Array(l3 + 1); let i; return (grid, offset, stride, factor) => { v[0] = 0; z[0] = -INF; z[1] = INF; for (let q = 0; q < l3; q++) { if (stride === 1) { // longitudes: repeat 3 times, consider the middle, scaled by the parallel’s length i = offset + (q % l); f[q] = grid[i] / factor; } else { // latitudes // repeat with a rotation of (180, 0, 180) const q0 = q % l; const panel = Math.floor(q / l); if (panel === 1) { i = offset + q0 * stride; } else { i = l - 1 - ((offset + (l >> 1)) % l) + (l - 1 - q0) * stride; } f[q] = grid[i]; } } for (let q = 1, k = 0, s = 0; q < l3; q++) { do { const r = v[k]; s = (f[q] - f[r] + q * q - r * r) / (q - r) / 2; } while (s <= z[k] && --k > -1); k++; v[k] = q; z[k] = s; z[k + 1] = INF; } for (let q = 0, k = 0; q < l3; q++) { while (z[k + 1] < q) k++; const r = v[k]; if (q >= l && q < l * 2) grid[offset + (q - l) * stride] = factor * (f[r] + (q - r) * (q - r)); } }; }33:[] 2f:{"title":"Spherical signed distance","cells":[{"id":0,"value":"# Spherical signed distance","pinned":false,"hidden":false,"mode":"md"},{"id":686,"value":"*#30DayMapChallenge — Day 27. Heatmap.*","pinned":false,"hidden":false,"mode":"md"},{"id":321,"value":"canvas","pinned":false,"hidden":false,"mode":"ojs"},{"id":286,"value":"viewof spherical = Inputs.toggle({ label: \"spherical\", value: \"on\" })","pinned":false,"hidden":false,"mode":"ojs"},{"id":692,"value":"viewof source = Inputs.select([\"110m\", \"50m\"], {\n label: \"geojson resolution\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":661,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":319,"value":"import { canvas } with { image } from \"@fil/polyconic-glsl\" // \"@fil/satellite-glsl\";","pinned":true,"hidden":false,"mode":"ojs"},{"id":217,"value":"sdf = {\n const edt = spherical ? sphedt1d : edt1d;\n\n const grid = Float64Array.from(originalGrid.ocean);\n\n for (let x = 0; x < width; x++) {\n // transform each column\n edt(height)(grid, x, width, 1);\n }\n for (let y = 0; y < height; y++) {\n // transform each row, scaled by latitude\n const f = Math.sin(((y + 0.5) / height) * Math.PI) ** 2;\n edt(width)(grid, y * width, 1, f);\n }\n\n // apply sign\n for (let i = 0, n = width * height; i < n; i++)\n grid[i] = Math.sqrt(grid[i]) * (originalGrid.sign[i] * 2 - 1);\n\n return grid;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":21,"value":"$32","pinned":true,"hidden":false,"mode":"ojs"},{"id":25,"value":"image = {\n const M = -d3.min(sdf) / d3.max(sdf);\n const canvas = heatmap(\n sdf.map((d) => (d < 0 ? d : d * M)),\n {\n width: width,\n height: height,\n color: interpolateOleron\n }\n );\n return d3.image(canvas.toDataURL());\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":702,"value":"originalGrid","pinned":true,"hidden":false,"mode":"ojs"},{"id":707,"value":"d3.extent(originalGrid.ocean)","pinned":true,"hidden":false,"mode":"ojs"},{"id":717,"value":"d3.extent(originalGrid.sign)","pinned":true,"hidden":false,"mode":"ojs"},{"id":3,"value":"viewof originalGrid = {\n const N = width * height;\n const ocean = new Float64Array(N);\n const sign = new Uint8Array(N);\n\n const context = DOM.context2d(width, height, 1);\n const path = d3.geoPath(projection, context);\n\n context.fillStyle = \"rgb(120, 170, 0)\";\n\n context.beginPath();\n path(land);\n\n context.stroke(); // stroke gives the border\n {\n const { data } = context.getImageData(0, 0, width, height);\n for (let i = 0; i < N; i++) ocean[i] = data[(i << 2) + 3] > 100 ? 0 : INF;\n }\n\n context.fill(); // fill gives the sign (interior/exterior)\n {\n const { data } = context.getImageData(0, 0, width, height);\n for (let i = 0; i < N; i++) sign[i] = data[(i << 2) + 1] > 90;\n }\n\n const value = { ocean, sign };\n return Object.assign(context.canvas, { value, style: \"max-width: 100%\" });\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":126,"value":"import { schemes } from \"@fil/colormaps\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":130,"value":"interpolateOleron = d3.scaleSequential((t) => schemes.oleron[(t * 255) | 0])","pinned":true,"hidden":false,"mode":"ojs"},{"id":27,"value":"import { heatmap } from \"@fil/heatmap\"","pinned":true,"hidden":false,"mode":"ojs"},{"id":50,"value":"INF = 10000000000","pinned":true,"hidden":false,"mode":"ojs"},{"id":626,"value":"width = 1200","pinned":true,"hidden":false,"mode":"ojs"},{"id":7,"value":"height = width / 2","pinned":true,"hidden":false,"mode":"ojs"},{"id":8,"value":"land = d3.json(\n `https://unpkg.com/visionscarto-world-atlas@0.0.6/world/${source}_land.geojson`\n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":13,"value":"projection = d3.geoEquirectangular().fitExtent(\n [\n [0, 0],\n [width, height]\n ],\n { type: \"Sphere\" }\n)","pinned":true,"hidden":false,"mode":"ojs"},{"id":71,"value":"d3 = require(\"d3@7\", \"d3-geo-projection@4\")","pinned":true,"hidden":false,"mode":"ojs"},{"id":471,"value":"// for reference, the original algo\nedt1d = (length) => {\n const v = new Uint16Array(length);\n const z = new Float64Array(length + 1);\n const f = new Float64Array(length);\n return (grid, offset, stride) => {\n v[0] = 0;\n z[0] = -INF;\n z[1] = INF;\n for (let q = 0; q < length; q++) f[q] = grid[offset + q * stride];\n\n for (let q = 1, k = 0, s = 0; q < length; q++) {\n do {\n const r = v[k];\n s = (f[q] - f[r] + q * q - r * r) / (q - r) / 2;\n } while (s <= z[k--]);\n\n k += 2;\n v[k] = q;\n z[k] = s;\n z[k + 1] = INF;\n }\n\n for (let q = 0, k = 0; q < length; q++) {\n while (z[k + 1] < q) k++;\n const r = v[k];\n grid[offset + q * stride] = f[r] + (q - r) * (q - r);\n }\n };\n}","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}