update volume viewer

_posts/2025-01-25-volume-rendering.md

@@ -161,6 +161,8 @@ f.savefig("hist.svg")
 
 We could probably get away with clamping all the data from -1000 to -500 to one air value, which would free up a lot of our limited 0-225 for the more interesting stuff happening between -100 and 400. But I didn't really notice an issues with the quantisation so I didn't pursue this. 
 
+EDIT: After I implemented the iso-surface rendering mode and found that I could see interesting regions like my windpipe and inside my sinuses I wondered if having more density precision would help see them. So I using float16 or float32 textures but didn't see much improvement at the expense of doubling or quadrupling the file size, so I switched back to 8 bit values. 
+
 ## Viewing the Data
 
 For the viewer I mostly copied the code from [this excellent tutorial](https://observablehq.com/@mroehlig/3d-volume-rendering-with-webgl-three-js) and integrated it into my existing three.js helper methods.

assets/js/outline-model-viewer/VolumeShaders.js

@@ -29,10 +29,12 @@ precision highp float; // Precision for floating point numbers.
 
 uniform sampler3D dataTexture; // Sampler for the volume data texture.
 // uniform sampler2D colorTexture; // Sampler for the color palette texture.
+uniform int renderMode; // Rendering mode.
 uniform float samplingRate; // The sampling rate.
 uniform float clampMin; // Clamp values below this value to 0.
 uniform float clampMax; // Clamp values above this value to 1.
-uniform float threshold; // Threshold to use for isosurface-style rendering.
+uniform float iso_threshold; // Threshold to use for isosurface-style rendering.
+uniform float iso_width; // Threshold to use for isosurface-style rendering.
 uniform float alphaScale; // Scaling of the color alpha value.
 uniform bool invertColor; // Option to invert the color palette.
 
@@ -70,9 +72,15 @@ vec2 intersectAABB(vec3 rayOrigin, vec3 rayDir, vec3 boxMin, vec3 boxMax) {
 // Volume sampling and composition.
 // Note that the code is inserted based on the selected algorithm in the user interface.
 vec4 compose(vec4 color, vec3 entryPoint, vec3 rayDir, float samples, float tStart, float tEnd, float tIncr) {
-  // Composition of samples using maximum intensity projection.
   // Loop through all samples along the ray.
-  float density = 0.0;
+  float max_density = 0.0;
+  float min_density = 1.0;
+
+  float mean_density = 0.0;
+  int mean_samples = 0;
+
+  float iso_depth = 0.0;
+
   for (float i = 0.0; i < samples; i += 1.0) {
     // Determine the sampling position.
     float t = tStart + tIncr * i; // Current distance along ray.
@@ -83,19 +91,39 @@ vec4 compose(vec4 color, vec3 entryPoint, vec3 rayDir, float samples, float tSta
     value = value < clampMin ? 0. : value;   
     value = value > clampMax ? 0. : value;
     
-
-    // Keep track of the maximum value.
-    if (value > density) {
-      // Store the value if it is greater than the previous values.
-      density = value;
+    if (value > max_density) {
+      max_density = value;
+    }
+    if (value < min_density && value > 0.0) {
+        min_density = value;
+    }
+    if (value > 0.0) {
+        mean_density += value;
+        mean_samples += 1;
+    }
+    if (abs(value - iso_threshold) < iso_width && iso_depth == 0.0) {
+        iso_depth = 1.;
     }
 
-    // Early exit the loop when the maximum possible value is found or the exit point is reached. 
-    if (density >= 1.0 || t > tEnd) {
+    // Early exit if the exit point is reached. 
+    if (t > tEnd) {
       break;
     }
   }
 
+  // Compute the final density value based on the selected rendering mode.
+  mean_density = mean_samples > 0 ? mean_density / float(mean_samples) : 0.0;
+  float density = 0.0;
+  if (renderMode == 0) {
+    density = max_density;
+  } else if (renderMode == 1) {
+    density = mean_density;
+  } else if (renderMode == 2) {
+    density = min_density;
+  } else if (renderMode == 3) {
+    density = iso_depth;
+  }
+
   // Convert the found value to a color by sampling the color palette texture.
   color.rgb = sampleColor(density).rgb;
   // Modify the alpha value of the color to make lower values more transparent.

assets/js/outline-model-viewer/VolumeViewer.js

@@ -11,56 +11,81 @@ import {
   deserialiseControls,
 } from "./helpers.js";
 
+// See https://stackoverflow.com/questions/62003464/what-is-relation-between-type-and-format-of-texture
+// https://webgl2fundamentals.org/webgl/lessons/webgl-data-textures.html
+const dtypes = {
+  uint8: {
+    internalFormat: "R8",
+    format: THREE.RedFormat,
+    type: THREE.UnsignedByteType,
+    array_type: Uint8Array,
+  },
+  float16: {
+    internalFormat: "R16F",
+    format: THREE.RedFormat,
+    type: THREE.HalfFloatType,
+    array_type: Uint16Array,
+  },
+  float32: {
+    internalFormat: "R32F",
+    format: THREE.RedFormat,
+    type: THREE.FloatType,
+    array_type: Float32Array,
+  },
+};
+
 async function load_metadata(metadata_path) {
   console.log("Loading metadata from", metadata_path);
   const metadata_res = await fetch(metadata_path);
   return await metadata_res.json();
 }
 
-async function load_model_bytes(model_path) {
-  console.log("Loading model from", model_path);
-  const res = await fetch(model_path);
-  const buffer = await res.arrayBuffer();
-  return new Uint8Array(buffer); // Create an uint8-array-view from the file buffer.
-}
-
-async function load_model_bytes_gzip(model_path, metadata_path, scene) {
+async function load_model_compressed_bytes(model_path) {
+  const model_response = await fetch(model_path);
   const ds = new DecompressionStream("gzip");
-  const response = await fetch(model_path);
-  const blob_in = await response.blob();
+  const blob_in = await model_response.blob();
   console.log("Compressed Model size", blob_in.size);
   const stream_in = blob_in.stream().pipeThrough(ds);
   const buffer = await new Response(stream_in).arrayBuffer();
   console.log("Decompressed Model size", buffer.byteLength);
-  return new Uint8Array(buffer);
+  return buffer;
+}
+
+async function load_model_bytes_gzip(model_path, metadata_path) {
+  const [metadata, model_buffer] = await Promise.all([
+    load_metadata(metadata_path),
+    load_model_compressed_bytes(model_path),
+  ]);
+
+  const array_type = dtypes[metadata.dtype].array_type;
+  return [metadata, new array_type(model_buffer)];
 }
 
 async function load_model(model_path, metadata_path, scene) {
   // If the model path ends in ".gz", we assume that the model is compressed.
-  const model_promise = model_path.endsWith(".gz")
-    ? load_model_bytes_gzip(model_path, metadata_path, scene)
-    : load_model_bytes(model_path);
-
-  const [byteArray, metadata] = await Promise.all([
-    model_promise,
-    load_metadata(metadata_path),
-  ]);
+  const [metadata, model_data] = await load_model_bytes_gzip(
+    model_path,
+    metadata_path
+  );
 
   console.log("Loaded model with metadata", metadata);
   console.log("Model shape", metadata.shape);
   console.log("Model dtype", metadata.dtype);
 
   const texture = new THREE.Data3DTexture(
-    byteArray, // The data values stored in the pixels of the texture.
+    model_data, // The data values stored in the pixels of the texture.
     metadata.shape[2], // Width of texture.
     metadata.shape[1], // Height of texture.
     metadata.shape[0] // Depth of texture.
   );
+  texture.internalFormat = dtypes[metadata.dtype].internalFormat;
+  texture.format = dtypes[metadata.dtype].format;
+  texture.type = dtypes[metadata.dtype].type;
 
-  texture.format = THREE.RedFormat; // Our texture has only one channel (red).
-  texture.type = THREE.UnsignedByteType; // The data type is 8 bit unsighed integer.
   texture.minFilter = THREE.LinearFilter; // Linear filter for minification.
   texture.magFilter = THREE.LinearFilter; // Linear filter for maximization.
+  //   texture.minFilter = THREE.NearestFilter; // Nearest filter for minification.
+  //   texture.magFilter = THREE.NearestFilter; // Nearest filter for maximization.
 
   // Repeat edge values when sampling outside of texture boundaries.
   texture.wrapS = THREE.ClampToEdgeWrapping;
@@ -93,13 +118,16 @@ function volumeMaterial(texture, renderProps) {
     uniforms: {
       dataTexture: { value: texture }, // Volume data texture.
       //   colorTexture: { value: colorTexture }, // Color palette texture.
+      renderMode: { value: renderProps.renderMode }, // Rendering mode.
       cameraPosition: { value: new THREE.Vector3() }, // Current camera position.
       samplingRate: { value: renderProps.samplingRate }, // Sampling rate of the volume.
 
       clampMin: { value: renderProps.clampMin }, // Clamp values below this value to 0.
       clampMax: { value: renderProps.clampMax }, // Clamp values above this value to 1.
 
-      threshold: { value: renderProps.threshold }, // Threshold for adjusting volume rendering.
+      iso_threshold: { value: renderProps.iso_threshold }, // Threshold for adjusting volume rendering.
+      iso_width: { value: renderProps.iso_width }, // Threshold for adjusting volume rendering.
+
       alphaScale: { value: renderProps.alphaScale }, // Alpha scale of volume rendering.
       invertColor: { value: renderProps.invertColor }, // Invert color palette.
     },
@@ -127,11 +155,21 @@ export class VolumeViewer extends HTMLElement {
     const box = make_box();
     scene.add(box);
 
+    const renderModes = {
+      "Max Intensity": 0,
+      "Mean Intensity": 1,
+      "Min Intensity": 2,
+      Isosurface: 3,
+    };
+
     let material = null;
     load_model(model, model_metadata, scene).then(({ texture, metadata }) => {
       // Create the custom material with attached shaders.
-      material = volumeMaterial(texture, renderProps);
+      material = volumeMaterial(texture, presets.Default);
       box.material = material;
+      gui
+        .add(material.uniforms.renderMode, "value", renderModes)
+        .name("Render Mode");
       gui
         .add(material.uniforms.samplingRate, "value", 0.1, 2.0, 0.1)
         .name("Sampling Rate");
@@ -142,23 +180,76 @@ export class VolumeViewer extends HTMLElement {
         .add(material.uniforms.clampMax, "value", 0.0, 1.0, 0.01)
         .name("Clamp Max");
       gui
-        .add(material.uniforms.threshold, "value", 0.0, 1.0, 0.01)
-        .name("Threshold");
+        .add(material.uniforms.iso_threshold, "value", 0.0, 1.0, 0.01)
+        .name("Isosurface Threshold");
+      gui
+        .add(material.uniforms.iso_width, "value", 0.0, 0.05, 0.001)
+        .name("Isosurface Width");
       gui
         .add(material.uniforms.alphaScale, "value", 0.1, 2.0, 0.1)
         .name("Alpha Scale");
       gui.add(material.uniforms.invertColor, "value").name("Invert Color");
     });
 
-    const renderProps = {
-      samplingRate: 1.0,
-      clampMin: 0.0,
-      clampMax: 1.0,
-      threshold: 0.1,
-      alphaScale: 1.0,
-      invertColor: false,
+    const presets = {
+      Default: {
+        renderMode: 0,
+        samplingRate: 1.0,
+        clampMin: 0.0,
+        clampMax: 1.0,
+        iso_threshold: 0.1,
+        iso_width: 0.01,
+        alphaScale: 1.0,
+        invertColor: false,
+      },
+      "Air Pockets": {
+        alphaScale: 2,
+        clampMax: 1,
+        clampMin: 0,
+        invertColor: false,
+        iso_threshold: 0.06,
+        iso_width: 0.002,
+        renderMode: 3,
+        samplingRate: 1,
+      },
     };
 
+    // Add a button to print the current settings to the console
+    gui
+      .add(
+        {
+          printSettings: () =>
+            console.log(
+              Object.fromEntries(
+                Object.keys(presets.Default).map((key) => [
+                  key,
+                  material?.uniforms[key]?.value,
+                ])
+              )
+            ),
+        },
+        "printSettings"
+      )
+      .name("Print Current Settings");
+
+    // Add a dropdown to select a preset
+    let renderProps = {
+      presets: "Default",
+    };
+    gui
+      .add(renderProps, "preset", presets)
+      .onChange((preset) => {
+        Object.keys(preset).forEach((key) => {
+          if (material.uniforms[key]) {
+            material.uniforms[key].value = preset[key];
+          } else {
+            console.warn(`No uniform found for ${key}`);
+          }
+        });
+        gui.controllers.forEach((control) => control.updateDisplay());
+      })
+      .name("Presets");
+
     const render = () => renderer.render(scene, this.camera);
     this.render = render;
 
@@ -188,6 +279,8 @@ export class VolumeViewer extends HTMLElement {
       this.camera.updateProjectionMatrix();
       renderer.setSize(canvas.clientWidth, canvas.clientHeight);
     };
+    this.onWindowResize();
+
     const timer = new Timer();
 
     const update = () => {

assets/js/outline-model-viewer/helpers.js

@@ -169,7 +169,7 @@ function componentHTML(component_rect) {
 
         #container.fullscreen .lil-gui.root {
             margin-top: 0;
-            width: 200px;
+            width: 50%;
         }
       </style>
     `;