DSM creation from 3D Scene

CesiumJS
1

Hi,

Cesium version: 1.131

I try to create a DSM (Digital Surface Model) as raster out of the existing 3D Scene in Cesium. The scene has 3D Tiles, Terrain, and GLTF objects in it and i would like to create a DSM for a ceratin extent out of it. The raster should have different resolution like 5x5m, 1x1m or 10x10m. The extent could be easily 1kmx1km. I started with scene.sampleHeightMostDetailed(), but for bigger extents it takes really long and that is why i am looking for a faster approach.

I tried so far:
Direct depth buffer sampling with world position reconstruction, which reconstructs world positions from the depth buffer and extracts heights.It works like that

  1. Setting up an orthographic camera positioned above the target area looking straight down

  2. Createe a post-process shader that:

    • Reads depth from Cesium’s depth buffer using czm_readDepth(czm_globeDepthTexture, ...)

    • Reconstructs world-space positions from depth using czm_windowToEyeCoordinates() + czm_inverseView

    • Extracts geodetic height from those world positions (ECEF coordinates)

    • Encodes the height values into RGB channels for readback

  3. Render the scene with the post-process stage applied to get the encoded height data

  4. Reads back pixels from the framebuffer using the WebGL context

  5. Decoding the RGB pixels back into height values and calculates statistics

Essentially, it’s leveraging Cesium’s depth buffer and built-in matrix uniforms (czm_inverseView, etc.) to convert rendered pixels into elevation data for the geographic extent.
here is a code snippet:

const worldPositionShader = `
      uniform sampler2D colorTexture;
      uniform float u_minHeight;
      uniform float u_heightRange;

      in vec2 v_textureCoordinates;

      void main() {

        // Read depth using Cesium's built-in function

        float depth = czm_readDepth(czm_globeDepthTexture, v_textureCoordinates);

        if (depth >= 0.9999) {
          // No geometry
          out_FragColor = vec4(1.0, 0.0, 1.0, 0.0);
          return;
        }

        // Reconstruct window coordinates

        vec2 windowCoord = v_textureCoordinates * czm_viewport.zw + czm_viewport.xy;

        // Get raw depth for position reconstruction

        float rawDepth = texture(czm_globeDepthTexture, v_textureCoordinates).r;

        // Convert window coordinates to eye space

        vec4 eyePosition = czm_windowToEyeCoordinates(windowCoord, rawDepth);

        // Transform from eye space to world space

        vec4 worldPosition = czm_inverseView * eyePosition;

        worldPosition /= worldPosition.w;


        // Convert Cartesian world position to Cartographic (lon, lat, height)

        // We need the height above ellipsoid

        vec3 posECEF = worldPosition.xyz;


        // Approximate height calculation using ellipsoid radius

        // For WGS84: semi-major axis = 6378137.0, semi-minor = 6356752.3142

        float a = 6378137.0;

        float b = 6356752.3142;

        // Distance from Earth center

        float r = length(posECEF);


        // Approximate geodetic height (simplified, assumes roughly spherical)

        // For more accuracy, we'd need full geodetic conversion

        float lat = asin(posECEF.z / r);

        float N = a / sqrt(1.0 - (1.0 - (b*b)/(a*a)) * sin(lat) * sin(lat));

        float height = r - N;

        

        // Adjust for latitude (rough approximation)

        float cosLat = cos(lat);

        if (cosLat > 0.001) {

          height = length(posECEF.xy) / cosLat - N;

        }

        

        // Normalize and encode height

        float heightNormalized = (height - u_minHeight) / u_heightRange;

        heightNormalized = clamp(heightNormalized, 0.0, 1.0);


        float encoded = heightNormalized * 16777215.0;

        float red = floor(encoded / 65536.0);

        float green = floor(mod(encoded, 65536.0) / 256.0);

        float blue = mod(encoded, 256.0);

        out_FragColor = vec4(red / 255.0, green / 255.0, blue / 255.0, 1.0);

      }

    `;


    const { minHeightRange = -1000, maxHeightRange = 10000 } = config;

    const heightRange = maxHeightRange - minHeightRange;

    postProcessStage = new PostProcessStage({

      fragmentShader: worldPositionShader,

      uniforms: {

        // eslint-disable-next-line @typescript-eslint/naming-convention

        u_minHeight: minHeightRange,

        // eslint-disable-next-line @typescript-eslint/naming-convention

        u_heightRange: heightRange,

      },

    });

    scene.postProcessStages.add(postProcessStage);

The second approach i tried is creating a DSM (Digital Surface Model) raster from the current Cesium scene by using an off-screen framebuffer with proper logarithmic depth decoding. It works like that:

  1. Saves the camera state and scene settings

  2. Enables depth features - Turns on depthTestAgainstTerrain and useDepthPicking so the depth buffer gets populated

  3. Positions camera above the target area looking straight down (keeps perspective frustum, not orthographic)

  4. Waits for terrain tiles to load at the new camera position, then renders several frames to ensure the depth buffer is populated

  5. Creates a post-process shader that:

    • Reads depth directly from gl_FragCoord.z

    • Uses czm_windowToEyeCoordinates() to convert window-space depth to eye-space position

    • Transforms to world space with czm_inverseView

    • Calculates ellipsoid height from the world position

    • Encodes the height into RGB channels

  6. Adds the shader as a post-process stage to the scene and renders

  7. Reads back the framebuffer pixels using the WebGL context

  8. Decodes RGB pixels back to height values and calculates statistics

  9. Restores the original camera and scene state

/**
* Post-process shader that reconstructs world position from GL depth and extracts ellipsoid height.
*
* PostProcessStage doesn’t provide depth texture access by default.
* Instead, we use gl_FragCoord.z to read the depth value directly from the depth buffer.
*
* For perspective projection:
* 1. Read depth directly from gl_FragCoord (window-space depth)
* 2. Use czm_windowToEyeCoordinates to get eye-space position
* 3. Transform to world coordinates using czm_inverseView
* 4. Calculate the ellipsoid height from the world position
*
* => encode the height into RGB channels for readback (24-bit precision).
*/
const depthToHeightShaderV2 = `
uniform sampler2D colorTexture;
uniform float u_minHeight;
uniform float u_heightRange;
uniform vec3 u_radiiSquared;

 in vec2 v_textureCoordinates;
  
  // Calculate height above ellipsoid from world position
  float getEllipsoidHeight(vec3 worldPos, vec3 radiiSquared) {
    // Normalize to ellipsoid surface
    vec3 oneOverRadiiSquared = 1.0 / radiiSquared;
    
    // Find the closest point on ellipsoid
    vec3 pos = worldPos;
    float beta = 1.0 / sqrt(
      pos.x * pos.x * oneOverRadiiSquared.x +
      pos.y * pos.y * oneOverRadiiSquared.y +
      pos.z * pos.z * oneOverRadiiSquared.z
    );
    
    vec3 surfacePoint = pos * beta;
    float surfaceHeight = length(surfacePoint);
    float actualHeight = length(worldPos);
    
    return actualHeight - surfaceHeight;
  }
  
  void main() {
    // WORKAROUND: Since PostProcessStage doesn't expose depth texture,
    // use the fragment's own depth value (gl_FragCoord.z)
    // This only works if the post-process stage writes to the same framebuffer
    float depth = gl_FragCoord.z;
    
    // Check for sky/no geometry
    if (depth >= 0.9999 || depth <= 0.0001) {
      out_FragColor = vec4(1.0, 0.0, 1.0, 0.0); // Magenta = no data
      return;
    }
    
    // Get window coordinates from texture coordinates
    vec2 windowCoord = v_textureCoordinates * czm_viewport.zw + czm_viewport.xy;
    
    // Convert window position + depth to eye coordinates
    // czm_windowToEyeCoordinates expects the depth value from the depth buffer
    vec4 eyeCoord = czm_windowToEyeCoordinates(windowCoord, depth);
    
    // Transform from eye space to world space
    vec4 worldCoord = czm_inverseView * eyeCoord;
    vec3 worldPos = worldCoord.xyz / worldCoord.w;
    
    // Calculate ellipsoid height
    float elevationMeters = getEllipsoidHeight(worldPos, u_radiiSquared);
    
    // Normalize and encode to RGB
    float heightNormalized = (elevationMeters - u_minHeight) / u_heightRange;
    heightNormalized = clamp(heightNormalized, 0.0, 1.0);
    
    // Encode to 24-bit RGB
    float encoded = heightNormalized * 16777215.0;
    float r = floor(encoded / 65536.0);
    float g = floor(mod(encoded, 65536.0) / 256.0);
    float b = mod(encoded, 256.0);
    
    out_FragColor = vec4(r / 255.0, g / 255.0, b / 255.0, 1.0);
  }
`;

// Get ellipsoid radii squared for height calculation
const { ellipsoid } = scene.globe;
const { radii } = ellipsoid;
const radiiSquared = new Cartesian3(
  radii.x * radii.x,
  radii.y * radii.y,
  radii.z * radii.z,
);

// Wait for terrain tiles to load at the new camera position
// eslint-disable-next-line no-console
console.log('Waiting for terrain tiles to load...');

// Helper to wait for tiles
const waitForTiles = async (maxWaitMs: number): Promise<void> => {
  const startTime = Date.now();
  let continueLoop = true;
  while (continueLoop) {
    scene.render();
    const { tilesLoaded } = scene.globe;

    const tilesToRender =
      (scene.globe as any)._surface?._tilesToRender?.length ?? 0;

    if (tilesLoaded && tilesToRender > 0) {
      // eslint-disable-next-line no-console
      console.log(`Tiles loaded. Tiles to render: ${tilesToRender}`);
      continueLoop = false;
    } else if (Date.now() - startTime > maxWaitMs) {
      // eslint-disable-next-line no-console
      console.log('Timeout waiting for tiles, proceeding anyway...');
      continueLoop = false;
    } else {
      // eslint-disable-next-line no-await-in-loop
      await new Promise<void>((resolve) => {
        setTimeout(() => {
          resolve();
        }, 100);
      });
    }
  }
};

await waitForTiles(3000);

// Render several more frames to ensure depth buffer is populated
// eslint-disable-next-line no-console
console.log('Rendering frames to populate depth buffer...');
for (let i = 0; i < 5; i += 1) {
  scene.render();
  // eslint-disable-next-line no-await-in-loop
  await new Promise<void>((resolve) => {
    setTimeout(() => {
      resolve();
    }, 50);
  });
}

// Debug: Check if globe depth texture exists
// eslint-disable-next-line @typescript-eslint/no-explicit-any, no-underscore-dangle
const frameState = (scene as any)._frameState;
// eslint-disable-next-line no-console
console.log('Frame state useLogDepth:', frameState?.useLogDepth);

// Check multiple places where depth might be stored
// eslint-disable-next-line @typescript-eslint/no-explicit-any, no-underscore-dangle
const globeDepth = (scene as any)._globeDepth;
// eslint-disable-next-line @typescript-eslint/no-explicit-any, no-underscore-dangle
const view = (scene as any)._view;
// eslint-disable-next-line no-underscore-dangle
const viewGlobeDepth = view?._globeDepth;

// eslint-disable-next-line no-console
console.log('Scene._globeDepth exists:', !!globeDepth);
// eslint-disable-next-line no-console
console.log('Scene._view exists:', !!view);
// eslint-disable-next-line no-console
console.log('View._globeDepth exists:', !!viewGlobeDepth);

if (viewGlobeDepth) {
  // eslint-disable-next-line no-console, no-underscore-dangle
  console.log('View globe depth framebuffer:', viewGlobeDepth._framebuffer);
  // eslint-disable-next-line no-console, no-underscore-dangle
  console.log('View globe depth texture:', viewGlobeDepth._colorTexture);
}

// Also check picking depth
// eslint-disable-next-line @typescript-eslint/no-explicit-any, no-underscore-dangle
const pickingDepth = (scene as any)._picking?.pickDepth;
// eslint-disable-next-line no-console
console.log('Picking depth exists:', !!pickingDepth);

// Create post-process stage with updated uniforms
postProcessStage = new PostProcessStage({
  fragmentShader: depthToHeightShaderV2,
  uniforms: {
    // eslint-disable-next-line @typescript-eslint/naming-convention
    u_minHeight: minHeightRange,
    // eslint-disable-next-line @typescript-eslint/naming-convention
    u_heightRange: heightRange,
    // eslint-disable-next-line @typescript-eslint/naming-convention
    u_radiiSquared: radiiSquared,
  },
});

// Add post-process stage
scene.postProcessStages.add(postProcessStage);

// Render with the post-process stage
scene.render();

// Wait a frame for async operations and render again
await new Promise<void>((resolve) => {
  setTimeout(() => {
    resolve();
  }, 100);
});
scene.render();

But for both approaches i am getting DSM statistics like this:

DSM created with statistics:
Min Height: -1000.00 m
Max Height: -1000.00 m
Center Height: -1000.00 m
Average Height: -1000.00 m
Valid Pixels: 1048576 / 1048576

So it does not work out. It is more or less R&D and just a prototype figuring out if that would be possible. Has anyone an idea if that would be kind of possible or is it just a stupid idea to create a DSM export like that not using sampleHeightMostDetailed()?

Many thanks in advance. Every comment is highly welcome :slight_smile:

BR Thomas