Fast Reversible Encoding/Decoding of DEM (16bit) into RGB

Copyright (C) 2011, Klokan Technologies GmbH. Author: Petr Pridal, Vaclav Klusak, Petr Sloup

The problem we are facing in WebGL Earth implementation of the terrain:

WebGL seems to not be able to read the 16-bit textures, which are normally used in OpenGL ES 2.0 for rendering of the terrain with elevation (in the vertex shader). To go over this limit we need to encode standard elevation model of Earth (-10,701m to 8,248m) into RGB space and prepare tile pyramid from the available data (such as SRTM3 Plus). The encoding algorithm should be:

1. very fast (only bit shifting and logical operations or basic algebra if possible, simpler means better)
2. symmetric (used colors must be decodable 1:1 to the elevation model and vice versa - without any loss)
3. optional: friendly to lossy compression (if the DEM data are saved in an image format which is using lossless compression then the differences between input elevation and decoded elevation should be minimal) - if this is not possible then lossless compression must be used (such as PNG)
4. optional: friendly to humans (people naturally expect sea with blue tones, land with green/yellow/brown, bright to dark for actual elevation values) - this can be practical for visual validation of correctness of the dataset too

Alternative A: Simple bit encoding of a float into RGB 8-bit

TERRAIN_MIN=-10000.0
TERRAIN_MAX=12000.0
r = pixel.data[0] / 256.0
g = pixel.data[1] / 256.0
elev = Math.ceil((TERRAIN_MAX-TERRAIN_MIN)*(r + g/256.0) - TERRAIN_MAX)

Encoding is possible with the utility gdaldem_web: https://github.com/webglearth/gdaldem_web

The alternative A is what we are now using in WebGL Earth.

Alternative B: Z-order curve

The suggested algorithm: The Z-order curve (known as Morton order) which goes trough the RGB cube in 2 dimensions - with the blue color set to minimal value for land and maximal for water. We have available value range (65536) which equals to 25 cm steps for elevation model. The tiles are rendered in popular mercator system, same as the map data such as OpenStreetMap.

Encoder/Decoder Source code:

function part1by1(n) {
    n &= 0x0000ffff
    n = (n | (n << 8)) & 0x00FF00FF
    n = (n | (n << 4)) & 0x0F0F0F0F
    n = (n | (n << 2)) & 0x33333333
    n = (n | (n << 1)) & 0x55555555
    return n
}

function unpart1by1(n) {
    n &= 0x55555555
    n = (n ^ (n >> 1)) & 0x33333333
    n = (n ^ (n >> 2)) & 0x0f0f0f0f
    n = (n ^ (n >> 4)) & 0x00ff00ff
    n = (n ^ (n >> 8)) & 0x0000ffff
    return n
}

function interleave2(x, y) {
    return part1by1(x) | (part1by1(y) << 1)
}

function deinterleave2(n) {
    return [ unpart1by1(n), unpart1by1(n >> 1) ]
}