30:T8ec,md`Given [my previous work on terrestrial and martian elevation data visualization](https://beta.observablehq.com/@ehouais/heightmap-rendering-using-a-floorcasting-algorithm), I was tempted at recreating this view and make it explorable. While ESA sadly doesn't give access to the raw elevation data (unlike University of Arizona's LPL does for its own [Hirise](https://hirise.lpl.arizona.edu) data), they do provide high resolution [photography](https://www.esa.int/spaceinimages/Images/2018/12/Plan_view_of_Korolev_crater) ([608x700 jpg](https://www.esa.int/var/esa/storage/images/esa_multimedia/images/2018/12/plan_view_of_korolev_crater/18937953-1-eng-GB/Plan_view_of_Korolev_crater_node_full_image_2.jpg) / [3476x4000 jpg](https://www.esa.int/var/esa/storage/images/esa_multimedia/images/2018/12/plan_view_of_korolev_crater/18937953-1-eng-GB/Plan_view_of_Korolev_crater.jpg) / [3476x4000 tiff](https://www.esa.int/var/esa/storage/images/esa_multimedia/images/2018/12/plan_view_of_korolev_crater/18937954-1-eng-GB/Plan_view_of_Korolev_crater.tif)), and [topography](https://www.esa.int/spaceinimages/Images/2018/12/Topography_of_Korolev_crater) ([608x700 jpg](https://www.esa.int/var/esa/storage/images/esa_multimedia/images/2018/12/topography_of_korolev_crater/18938141-1-eng-GB/Topography_of_Korolev_crater_node_full_image_2.jpg) / [3475x4000 jpg](https://www.esa.int/var/esa/storage/images/esa_multimedia/images/2018/12/topography_of_korolev_crater/18938141-1-eng-GB/Topography_of_Korolev_crater.jpg) / [3475x4000 tiff](https://www.esa.int/var/esa/storage/images/esa_multimedia/images/2018/12/topography_of_korolev_crater/18938142-1-eng-GB/Topography_of_Korolev_crater.tif)) that mixes color-coded elevation and shading, and provides the color ramp that was used. plan photograph

`31:[] 2e:{"title":"Korolev crater","cells":[{"id":0,"value":"md`# Korolev crater`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1234,"value":"md`In december 2018, The European Space Agency (ESA) [showcased a beautiful view](https://www.esa.int/spaceinimages/Images/2018/12/Perspective_view_of_Korolev_crater) of [the Korolev impact crater on Mars](https://en.wikipedia.org/wiki/Korolev_(Martian_crater).\n\n![Korolev crater](https://www.esa.int/var/esa/storage/images/esa_multimedia/images/2018/12/perspective_view_of_korolev_crater/18938178-1-eng-GB/Perspective_view_of_Korolev_crater_node_full_image_2.jpg)\n\n This view is obviously *not* a photograph, [as it sometimes was wrongly advertised](https://mashable.com/article/mars-ice-korolev-crater-photo-esa), but an image generated from [Mars Express](https://en.wikipedia.org/wiki/Mars_Express)' [High Resolution Stereo Camera (HRSC)](https://en.wikipedia.org/wiki/High_Resolution_Stereo_Camera) color channels and elevation data.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1239,"value":"md`## What if...`","pinned":false,"hidden":false,"mode":"ojs"},{"id":445,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":485,"value":"md`So the goal of this notebook is to reverse-engineer the topography view and its color-ramp to recreate bare elevation data, hopefully without too much data loss and artefacts, then to visualize the result using standard heightmap floorcasting algorithms.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":128,"value":"md`## Topography image data`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1112,"value":"md`First we import and resize the high-resolution topography image data from the given ESA link (with CORS proxyfication). TIFF images are not natively decoded by browsers so I'm using [a specific library](https://beta.observablehq.com/@ehouais/decoding-tiff-image-data).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1362,"value":"resize = 1500","pinned":false,"hidden":false,"mode":"ojs"},{"id":2002,"value":"viewof topoTiffBuffer = {\n await visibility()\n return fetchProgress('https://cadasa.github.io/Topography_of_Korolev_crater.tif', {}, 23880416, 'Loading topographic data')\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":2,"value":"topographyHR = tiffArrayBufferToImageData(topoTiffBuffer)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1314,"value":"topography = resizeImageData(topographyHR, 0, 263, 3475, 3475, resize, resize)","pinned":false,"hidden":false,"mode":"ojs"},{"id":149,"value":"md`## The color ramp`","pinned":false,"hidden":false,"mode":"ojs"},{"id":277,"value":"md`Since we don't know the function that encodes altitude into RGB color in the topography view, we'll try to infer this correspondence from the color ramp displayed in the top-right corner of the image.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":416,"value":"md`First we grab the interesting part from the bitmap:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":434,"value":"rampBounds = ({x: 2272, y: 53, width: 1141, height: 40})","pinned":false,"hidden":false,"mode":"ojs"},{"id":151,"value":"{\n const { x, y, width: W, height: H } = rampBounds;\n const scontext = DOM.context2d(topographyHR.width, topographyHR.height, 1);\n scontext.putImageData(topographyHR, 0, 0);\n const context = DOM.context2d(width, H, 1);\n context.drawImage(scontext.canvas, x, y, W, H, 0, 0, width, H);\n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":422,"value":"md`Now we extract the RGB and HSV components from an horizontal line taken from this ramp. \nAll the rows don't seem to display exactly the same ramp, so we choose the first one, hoping this choice doesn't affect the result much.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":766,"value":"τ = 2 * Math.PI","pinned":false,"hidden":false,"mode":"ojs"},{"id":196,"value":"ramp = {\n const { x, y, width: W, height: H } = rampBounds;\n const data = topographyHR.data;\n let index = (y * topographyHR.width + x) * 4;\n const hsls = [];\n for (let i = 0; i < W; i++) {\n let [R, G, B] = [data[index++], data[index++], data[index++]];\n R /= 256;\n G /= 256;\n B /= 256;\n const [H, S, V] = rgb2hsv([R, G, B]);\n const [XH, YH] = [S * Math.cos(H * τ), S * Math.sin(H * τ)];\n index++;\n hsls.push({\n x: i,\n ...{ R, G, B },\n ...{ H, S, V, XH, YH },\n color: `rgb(${R},${G},${B})`\n });\n }\n return hsls;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":419,"value":"md`Then we visualize how the RGB and HSV components vary when elevation increases:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":204,"value":"vegalite({\n $schema: \"https://vega.github.io/schema/vega-lite/v3.json\",\n data: { values: ramp },\n transform: [{ fold: [\"R\", \"G\", \"B\"], as: [\"component\", \"value\"] }],\n width: width * 0.85,\n mark: \"line\",\n encoding: {\n x: {\n field: \"x\",\n type: \"ordinal\",\n axis: { ticks: false, labels: false },\n title: \"elevation\"\n },\n y: {\n field: \"value\",\n type: \"quantitative\",\n title: null,\n scale: { domain: [0, 1] }\n },\n color: {\n field: \"component\",\n type: \"nominal\",\n scale: {\n domain: [\"R\", \"G\", \"B\"],\n range: [\"red\", \"green\", \"blue\"]\n }\n }\n }\n})","pinned":true,"hidden":false,"mode":"ojs"},{"id":630,"value":"vegalite({\n $schema: \"https://vega.github.io/schema/vega-lite/v3.json\",\n data: { values: ramp },\n transform: [{ fold: [\"H\", \"S\", \"V\"], as: [\"component\", \"value\"] }],\n width: width * 0.85,\n mark: \"line\",\n encoding: {\n x: {\n field: \"x\",\n type: \"ordinal\",\n axis: { ticks: false, labels: false },\n title: \"elevation\"\n },\n y: {\n field: \"value\",\n type: \"quantitative\",\n title: null\n },\n color: {\n field: \"component\",\n type: \"nominal\",\n scale: {\n domain: [\"H\", \"S\", \"V\"],\n range: [\"cyan\", \"magenta\", \"yellow\"]\n }\n }\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1149,"value":"md`The hue (H) component from HSV looks promising since it is in bijection with elevation on roughly 3/4 of the ramp, which allows us to reverse it and get the color-to-elevation relation we need. After all, this is the point of color-coding altitude: the hue of a pixel provides 2D elevation information to human eye without using perspective.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1180,"value":"md`## Color distribution`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1156,"value":"md`The problem comes from the fact that topography *also* includes shading information in pixels to provide another hint at topographic relief to the human eye. \nWe could try to remove the shading information from topography data by supposing that shading information comes from photographic data, then guessing how elevation and photography were mixed, and finally trying to reverse the process. \nBut this doesn't give good results and it's actually not necessary, as is shown by the HSV distributions (in cylindric coordinates) in the ramp and in the whole topographic image (sampled at 1%).`","pinned":false,"hidden":false,"mode":"ojs"},{"id":376,"value":"plotHSV = (div, title, rows, [W, H]) => {\n const column = (rows, key) => rows.map(row => row[key]);\n const [X, Y, Z] = [\"XH\", \"YH\", \"V\"];\n const data = [\n {\n x: column(rows, X),\n y: column(rows, Y),\n z: column(rows, Z),\n mode: \"markers\",\n marker: {\n color: column(rows, \"color\"),\n size: 3,\n symbol: \"circle\"\n },\n type: \"scatter3d\"\n }\n ];\n\n Plotly.react(div, data, {\n width: W,\n height: H,\n margin: { l: 10, r: 10, b: 10 },\n title,\n scene: {\n xaxis: { title: X, range: [-1, 1] },\n yaxis: { title: Y, range: [-1, 1] },\n zaxis: { title: Z, range: [0, 1] },\n aspectmode: \"cube\"\n }\n });\n div.style.display = \"inline-block\";\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":367,"value":"sample = {\n const sample = [];\n let index = 0;\n decodeImage(topography, ([i, j], [R, G, B]) => {\n if (index++ % 100 == 0) {\n const [H, S, V] = rgb2hsv([R / 256, G / 256, B / 256]);\n const [XH, YH] = [S * Math.cos(H * τ), S * Math.sin(H * τ)];\n sample.push({ H, S, XH, YH, V, color: `rgb(${R},${G},${B})` });\n }\n });\n return sample;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":391,"value":"{\n let container = this;\n let rampdiv, imagediv;\n if (!container) {\n container = DOM.element(\"div\");\n rampdiv = DOM.element(\"div\");\n imagediv = DOM.element(\"div\");\n container.append(rampdiv, imagediv);\n invalidation.then(() => {\n Plotly.purge(rampdiv);\n Plotly.purge(imagediv);\n });\n }\n const children = container.children;\n const W = width / 2;\n plotHSV(children[0], \"Ramp HSV colors\", ramp, [W, W]);\n plotHSV(children[1], \"Image HSV colors (1% sample)\", sample, [W, W]);\n\n return container;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1194,"value":"md`This confirms that Hue (angle) is a good candidate for decoding elevation since shading is apparently encoded in Saturation (radius) and Value (Z position). \nThis also shows that the part of the ramp that is not reverse-coded (brownish hues) corresponds to a tiny percentage of the topographic pixels, letting us think that ignoring them won't affect the result too much.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1198,"value":"md`## Elevation look-up table`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1200,"value":"md`So we inverse the ramp into a look-up table to decode hue into elevation in the [0-860] elevation range.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":862,"value":"lut = {\n const max = 860;\n return ramp\n .slice(0, max)\n .reduce((map, { x, H }) => map.set(Math.round(H * 1000), x), new Map());\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":808,"value":"vegalite({\n $schema: \"https://vega.github.io/schema/vega-lite/v3.json\",\n data: {\n values: [...lut.entries()].map(([hue, elevation]) => ({\n hue: hue / 1000,\n elevation\n }))\n },\n width: width * 0.9,\n mark: { type: \"point\", size: 2 },\n encoding: {\n x: { field: \"hue\", type: \"quantitative\" },\n y: { field: \"elevation\", type: \"quantitative\" }\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1203,"value":"md`We must then interpolate the table to add the missing hue values`","pinned":false,"hidden":false,"mode":"ojs"},{"id":908,"value":"lut2 = {\n const table = [];\n let last = [0];\n for (let i = 0; i <= 808; i++) {\n if (lut.has(i + 9)) {\n const elev = lut.get(i + 9);\n let [lasth, laste] = last;\n if (i - lasth > 1) {\n const de = (elev - laste) / (i - lasth);\n while (++lasth < i) {\n table[lasth] = laste += de;\n }\n }\n table[i] = elev;\n last = [i, elev];\n }\n }\n return table;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":867,"value":"vegalite({\n $schema: \"https://vega.github.io/schema/vega-lite/v3.json\",\n data: {\n values: lut2.map((elevation, i) => ({\n hue: i / 1000,\n elevation: elevation || 0\n }))\n },\n width: width * 0.9,\n mark: { type: \"point\", size: 2 },\n encoding: {\n x: { field: \"hue\", type: \"quantitative\" },\n y: { field: \"elevation\", type: \"quantitative\" }\n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1232,"value":"md`## Elevation data`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1245,"value":"md`We are now ready to decode elevation from topographic colors, by extracting each pixel's hue component and converting it into elevation thanks to the look-up table.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":788,"value":"viewof elevation = {\n const { data, width: W, height: H } = topography;\n const { bar, set, dispatch } = viewofProgressBar('Converting hue to elevation');\n const E = [];\n let last;\n\n for (let j = 0, index = 0; j < H; j++) {\n E[j] = [];\n last = 0;\n for (let i = 0; i < W; i++, index += 4) {\n const [h, s, v] = rgb2hsv(data.slice(index, index + 3));\n let z = lut2[Math.floor(h * 1000)] || 0;\n // Quick fix for high luminosity pixel artefacts\n z = last > 200 && z - last > 40 ? last : z;\n last = E[j][i] = z;\n }\n set(1+100*j/H);\n yield bar;\n }\n dispatch(E);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":985,"value":"{\n const [dims, r] = [[500, 500], 255 / 860];\n return canvas(dims, (i, j) => Array(3).fill(elevation[i * 3][j * 3] * r));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1297,"value":"md`## Photographic data`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1304,"value":"md`We import photographic data at the same resolution, and quantized (5 to 1) to get a smaller palette.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2014,"value":"viewof photoTiffBuffer = {\n await visibility()\n return fetchProgress('https://cadasa.github.io/Plan_view_of_Korolev_crater.tif', {}, 22465516, 'Loading photographic data')\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":526,"value":"photography = resizeImageData(tiffArrayBufferToImageData(photoTiffBuffer), 0, 263, 3475, 3475, resize, resize)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1522,"value":"colormap = {\n const map = new Map();\n const palette = [];\n const quantize = (ratio => value => Math.floor(value / ratio) * ratio)(5);\n\n const data = [];\n for (let j = 0, index = 0; j < photography.height; j++) data[j] = [];\n\n decodeImage(photography, ([i, j], [r, g, b]) => {\n const rgb = [quantize(r), quantize(g), quantize(b)];\n const key = rgb.join();\n if (!map.has(key)) map.set(key, palette.push(rgb) - 1);\n data[j][i] = map.get(key);\n });\n\n return { data, palette };\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1052,"value":"{\n const { data, palette } = colormap;\n return canvas([500, 500], (i, j) => palette[data[i * 3][j * 3]]);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1300,"value":"md`## Perspective view`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1260,"value":"md`And finally we combine elevation and photography to recreate perspective using floorcasting. The vertical scaling is chosen so that the result resembles the original from ESA.`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1454,"value":"md `[click & move to rotate, wheel to zoom/unzoom]`","pinned":false,"hidden":false,"mode":"ojs"},{"id":942,"value":"viewof model = {\n const heightmap = { data: elevation, hunit: 1, ratio: resize / 15e3 };\n const dims = [width, (width * (3 / 5)) | 0];\n return orbitalCanvas({ heightmap, colormap }, dims);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1865,"value":"md`Let's add an arbitrary sky color and progressive distance fog to augment realism:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1461,"value":"md`[left-click/wheel to move, middle-click to rotate]`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1369,"value":"{\n const colors = { ...colormap, background: [240, 240, 235] };\n const dims = [width, (width * 0.6) | 0];\n const start = { x: 1107, y: 319, z: 65, θ: 20, φ: -0.05 };\n return explorer(elevation, colors, resize / 15e3, dims, start);\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1292,"value":"md`In retrospective, the process that produces the perspective view from topography and photography is the following:`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1275,"value":"dot`digraph imageLoading {\n rankdir = \"LR\"\n node [\n fontname = \"Arial\"\n fontsize = 10\n shape = \"record\"\n ]\n edge [\n fontname = \"Arial\"\n fontsize = 10\n ]\n\n \"Topography HR\" -> \"Color ramp\"\n \"Color ramp\" -> \"Hue component\\ndistribution\"\n \"Hue component\\ndistribution\" -> \"Look-up\\ntable\"\n \"Topography HR\" -> Topography -> \"Elevation data\"\n \"Look-up\\ntable\" -> \"Elevation data\"\n \"Elevation data\" -> \"Perspective\\nview\"\n \"Photography HR\" -> Photography -> \"Perspective\\nview\"\n}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":1471,"value":"md`## TODOs\n* Proper fix for remaining artefacts\n* Testing of other color ramp lines\n* Proper assessment of vertical scaling\n* Touch events for handheld devices\n* Use of GPU.js to accelerate computations\n* Full-screen`","pinned":false,"hidden":false,"mode":"ojs"},{"id":202,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":144,"value":"import {\n image2canvas,\n decodeImage,\n getImageData,\n canvas,\n rgb2hsv,\n fillArray,\n resizeImageData,\n viewofProgressBar,\n fetchProgress\n} from \"@ehouais/utils\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":113,"value":"import { tiffArrayBufferToImageData } from \"@ehouais/decoding-tiff-image-data\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1533,"value":"import {\n orbitalCanvas,\n explorer\n} from \"@ehouais/heightmap-rendering-using-a-floorcasting-algorithm\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":210,"value":"//vegalite = require(\"@observablehq/vega-lite@0.0\")\nvegalite = require(\"vega-embed@3.26.1\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":334,"value":"Plotly = require(\"https://cdn.plot.ly/plotly-latest.min.js\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1277,"value":"import {dot} from \"@mbostock/graphviz\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":1943,"value":"GPU = require('https://bundle.run/gpu.js@1.8.1')","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W31","stdlib":"1","titling":"auto"}