Support float and half-float data type

- direct-volume rendering
- GPU isosurface extraction

CHANGELOG.md

@@ -5,14 +5,15 @@ Note that since we don't clearly distinguish between a public and private interf
 
 ## [Unreleased]
 
-- Volume UI improvements 
+- Volume UI improvements
   - Render all volume entries instead of selecting them one-by-one
   - Toggle visibility of all volumes
   - More accessible iso value control
 - Support wheel event on sliders
 - MolViewSpec extension:
   - Add validation for discriminated union params
   - Primitives: remove triangle_colors, line_colors, have implicit grouping instead; rename many parameters
+- Support float and half-float data type for direct-volume rendering and GPU isosurface extraction
 
 ## [v4.10.0] - 2024-12-15
 

src/mol-geo/geometry/direct-volume/direct-volume.ts

@@ -49,6 +49,7 @@ export interface DirectVolume {
     readonly cartnToUnit: ValueCell<Mat4>
     readonly packedGroup: ValueCell<boolean>
     readonly axisOrder: ValueCell<Vec3>
+    readonly dataType: ValueCell<'byte' | 'float' | 'halfFloat'>
 
     /** Bounding sphere of the volume */
     readonly boundingSphere: Sphere3D
@@ -57,10 +58,10 @@ export interface DirectVolume {
 }
 
 export namespace DirectVolume {
-    export function create(bbox: Box3D, gridDimension: Vec3, transform: Mat4, unitToCartn: Mat4, cellDim: Vec3, texture: Texture, stats: Grid['stats'], packedGroup: boolean, axisOrder: Vec3, directVolume?: DirectVolume): DirectVolume {
+    export function create(bbox: Box3D, gridDimension: Vec3, transform: Mat4, unitToCartn: Mat4, cellDim: Vec3, texture: Texture, stats: Grid['stats'], packedGroup: boolean, axisOrder: Vec3, dataType: 'byte' | 'float' | 'halfFloat', directVolume?: DirectVolume): DirectVolume {
         return directVolume ?
-            update(bbox, gridDimension, transform, unitToCartn, cellDim, texture, stats, packedGroup, axisOrder, directVolume) :
-            fromData(bbox, gridDimension, transform, unitToCartn, cellDim, texture, stats, packedGroup, axisOrder);
+            update(bbox, gridDimension, transform, unitToCartn, cellDim, texture, stats, packedGroup, axisOrder, dataType, directVolume) :
+            fromData(bbox, gridDimension, transform, unitToCartn, cellDim, texture, stats, packedGroup, axisOrder, dataType);
     }
 
     function hashCode(directVolume: DirectVolume) {
@@ -71,7 +72,7 @@ export namespace DirectVolume {
         ]);
     }
 
-    function fromData(bbox: Box3D, gridDimension: Vec3, transform: Mat4, unitToCartn: Mat4, cellDim: Vec3, texture: Texture, stats: Grid['stats'], packedGroup: boolean, axisOrder: Vec3): DirectVolume {
+    function fromData(bbox: Box3D, gridDimension: Vec3, transform: Mat4, unitToCartn: Mat4, cellDim: Vec3, texture: Texture, stats: Grid['stats'], packedGroup: boolean, axisOrder: Vec3, dataType: 'byte' | 'float' | 'halfFloat'): DirectVolume {
         const boundingSphere = Sphere3D();
         let currentHash = -1;
 
@@ -103,6 +104,7 @@ export namespace DirectVolume {
             },
             packedGroup: ValueCell.create(packedGroup),
             axisOrder: ValueCell.create(axisOrder),
+            dataType: ValueCell.create(dataType),
             setBoundingSphere(sphere: Sphere3D) {
                 Sphere3D.copy(boundingSphere, sphere);
                 currentHash = hashCode(directVolume);
@@ -111,7 +113,7 @@ export namespace DirectVolume {
         return directVolume;
     }
 
-    function update(bbox: Box3D, gridDimension: Vec3, transform: Mat4, unitToCartn: Mat4, cellDim: Vec3, texture: Texture, stats: Grid['stats'], packedGroup: boolean, axisOrder: Vec3, directVolume: DirectVolume): DirectVolume {
+    function update(bbox: Box3D, gridDimension: Vec3, transform: Mat4, unitToCartn: Mat4, cellDim: Vec3, texture: Texture, stats: Grid['stats'], packedGroup: boolean, axisOrder: Vec3, dataType: 'byte' | 'float' | 'halfFloat', directVolume: DirectVolume): DirectVolume {
         const width = texture.getWidth();
         const height = texture.getHeight();
         const depth = texture.getDepth();
@@ -129,6 +131,7 @@ export namespace DirectVolume {
         ValueCell.update(directVolume.cartnToUnit, Mat4.invert(Mat4(), unitToCartn));
         ValueCell.updateIfChanged(directVolume.packedGroup, packedGroup);
         ValueCell.updateIfChanged(directVolume.axisOrder, Vec3.fromArray(directVolume.axisOrder.ref.value, axisOrder, 0));
+        ValueCell.updateIfChanged(directVolume.dataType, dataType);
         return directVolume;
     }
 
@@ -142,7 +145,8 @@ export namespace DirectVolume {
         const stats = Grid.One.stats;
         const packedGroup = false;
         const axisOrder = Vec3.create(0, 1, 2);
-        return create(bbox, gridDimension, transform, unitToCartn, cellDim, texture, stats, packedGroup, axisOrder, directVolume);
+        const dataType = 'byte';
+        return create(bbox, gridDimension, transform, unitToCartn, cellDim, texture, stats, packedGroup, axisOrder, dataType, directVolume);
     }
 
     export const Params = {

src/mol-gl/compute/marching-cubes/active-voxels.ts

@@ -1,13 +1,13 @@
 /**
- * Copyright (c) 2019-2022 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2019-2024 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
 
 import { ComputeRenderable, createComputeRenderable } from '../../renderable';
 import { WebGLContext } from '../../webgl/context';
 import { createComputeRenderItem } from '../../webgl/render-item';
-import { Values, TextureSpec, UniformSpec } from '../../renderable/schema';
+import { Values, TextureSpec, UniformSpec, DefineSpec } from '../../renderable/schema';
 import { Texture } from '../../../mol-gl/webgl/texture';
 import { ShaderCode } from '../../../mol-gl/shader-code';
 import { ValueCell } from '../../../mol-util';
@@ -17,12 +17,14 @@ import { getTriCount } from './tables';
 import { quad_vert } from '../../../mol-gl/shader/quad.vert';
 import { activeVoxels_frag } from '../../../mol-gl/shader/marching-cubes/active-voxels.frag';
 import { isTimingMode } from '../../../mol-util/debug';
+import { isWebGL2 } from '../../webgl/compat';
 
 const ActiveVoxelsSchema = {
     ...QuadSchema,
 
     tTriCount: TextureSpec('image-uint8', 'alpha', 'ubyte', 'nearest'),
     tVolumeData: TextureSpec('texture', 'rgba', 'ubyte', 'nearest'),
+    dValueChannel: DefineSpec('string', ['red', 'alpha']),
     uIsoValue: UniformSpec('f'),
 
     uGridDim: UniformSpec('v3'),
@@ -34,12 +36,17 @@ type ActiveVoxelsValues = Values<typeof ActiveVoxelsSchema>
 
 const ActiveVoxelsName = 'active-voxels';
 
+function valueChannel(ctx: WebGLContext, volumeData: Texture) {
+    return isWebGL2(ctx.gl) && volumeData.format === ctx.gl.RED ? 'red' : 'alpha';
+}
+
 function getActiveVoxelsRenderable(ctx: WebGLContext, volumeData: Texture, gridDim: Vec3, gridTexDim: Vec3, isoValue: number, scale: Vec2): ComputeRenderable<ActiveVoxelsValues> {
     if (ctx.namedComputeRenderables[ActiveVoxelsName]) {
         const v = ctx.namedComputeRenderables[ActiveVoxelsName].values as ActiveVoxelsValues;
 
         ValueCell.update(v.uQuadScale, scale);
         ValueCell.update(v.tVolumeData, volumeData);
+        ValueCell.update(v.dValueChannel, valueChannel(ctx, volumeData));
         ValueCell.updateIfChanged(v.uIsoValue, isoValue);
         ValueCell.update(v.uGridDim, gridDim);
         ValueCell.update(v.uGridTexDim, gridTexDim);
@@ -59,6 +66,7 @@ function createActiveVoxelsRenderable(ctx: WebGLContext, volumeData: Texture, gr
 
         uQuadScale: ValueCell.create(scale),
         tVolumeData: ValueCell.create(volumeData),
+        dValueChannel: ValueCell.create(valueChannel(ctx, volumeData)),
         uIsoValue: ValueCell.create(isoValue),
         uGridDim: ValueCell.create(gridDim),
         uGridTexDim: ValueCell.create(gridTexDim),

src/mol-gl/compute/marching-cubes/isosurface.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2019-2022 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2019-2024 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
@@ -28,6 +28,7 @@ const IsosurfaceSchema = {
     tActiveVoxelsPyramid: TextureSpec('texture', 'rgba', 'float', 'nearest'),
     tActiveVoxelsBase: TextureSpec('texture', 'rgba', 'float', 'nearest'),
     tVolumeData: TextureSpec('texture', 'rgba', 'ubyte', 'nearest'),
+    dValueChannel: DefineSpec('string', ['red', 'alpha']),
     uIsoValue: UniformSpec('f'),
 
     uSize: UniformSpec('f'),
@@ -48,13 +49,18 @@ type IsosurfaceValues = Values<typeof IsosurfaceSchema>
 
 const IsosurfaceName = 'isosurface';
 
+function valueChannel(ctx: WebGLContext, volumeData: Texture) {
+    return isWebGL2(ctx.gl) && volumeData.format === ctx.gl.RED ? 'red' : 'alpha';
+}
+
 function getIsosurfaceRenderable(ctx: WebGLContext, activeVoxelsPyramid: Texture, activeVoxelsBase: Texture, volumeData: Texture, gridDim: Vec3, gridTexDim: Vec3, transform: Mat4, isoValue: number, levels: number, scale: Vec2, count: number, invert: boolean, packedGroup: boolean, axisOrder: Vec3, constantGroup: boolean): ComputeRenderable<IsosurfaceValues> {
     if (ctx.namedComputeRenderables[IsosurfaceName]) {
         const v = ctx.namedComputeRenderables[IsosurfaceName].values as IsosurfaceValues;
 
         ValueCell.update(v.tActiveVoxelsPyramid, activeVoxelsPyramid);
         ValueCell.update(v.tActiveVoxelsBase, activeVoxelsBase);
         ValueCell.update(v.tVolumeData, volumeData);
+        ValueCell.update(v.dValueChannel, valueChannel(ctx, volumeData));
 
         ValueCell.updateIfChanged(v.uIsoValue, isoValue);
         ValueCell.updateIfChanged(v.uSize, Math.pow(2, levels));
@@ -87,6 +93,7 @@ function createIsosurfaceRenderable(ctx: WebGLContext, activeVoxelsPyramid: Text
         tActiveVoxelsPyramid: ValueCell.create(activeVoxelsPyramid),
         tActiveVoxelsBase: ValueCell.create(activeVoxelsBase),
         tVolumeData: ValueCell.create(volumeData),
+        dValueChannel: ValueCell.create(valueChannel(ctx, volumeData)),
 
         uIsoValue: ValueCell.create(isoValue),
         uSize: ValueCell.create(Math.pow(2, levels)),

src/mol-gl/shader/marching-cubes/active-voxels.frag.ts

@@ -38,24 +38,29 @@ vec4 texture3dFrom2dNearest(sampler2D tex, vec3 pos, vec3 gridDim, vec2 texDim)
     return texture2D(tex, coord);
 }
 
-vec4 voxel(vec3 pos) {
+float voxelValue(vec3 pos) {
     pos = min(max(vec3(0.0), pos), uGridDim - vec3(1.0));
-    return texture3dFrom2dNearest(tVolumeData, pos / uGridDim, uGridDim, uGridTexDim.xy);
+    vec4 v = texture3dFrom2dNearest(tVolumeData, pos / uGridDim, uGridDim, uGridTexDim.xy);
+    #ifdef dValueChannel_red
+        return v.r;
+    #else
+        return v.a;
+    #endif
 }
 
 void main(void) {
     vec2 uv = gl_FragCoord.xy / uGridTexDim.xy;
     vec3 posXYZ = index3dFrom2d(uv);
 
     // get MC case as the sum of corners that are below the given iso level
-    float c = step(voxel(posXYZ).a, uIsoValue)
-        + 2. * step(voxel(posXYZ + c1).a, uIsoValue)
-        + 4. * step(voxel(posXYZ + c2).a, uIsoValue)
-        + 8. * step(voxel(posXYZ + c3).a, uIsoValue)
-        + 16. * step(voxel(posXYZ + c4).a, uIsoValue)
-        + 32. * step(voxel(posXYZ + c5).a, uIsoValue)
-        + 64. * step(voxel(posXYZ + c6).a, uIsoValue)
-        + 128. * step(voxel(posXYZ + c7).a, uIsoValue);
+    float c = step(voxelValue(posXYZ), uIsoValue)
+        + 2. * step(voxelValue(posXYZ + c1), uIsoValue)
+        + 4. * step(voxelValue(posXYZ + c2), uIsoValue)
+        + 8. * step(voxelValue(posXYZ + c3), uIsoValue)
+        + 16. * step(voxelValue(posXYZ + c4), uIsoValue)
+        + 32. * step(voxelValue(posXYZ + c5), uIsoValue)
+        + 64. * step(voxelValue(posXYZ + c6), uIsoValue)
+        + 128. * step(voxelValue(posXYZ + c7), uIsoValue);
     c *= step(c, 254.);
 
     // handle out of bounds positions

src/mol-gl/shader/marching-cubes/isosurface.frag.ts

@@ -59,9 +59,14 @@ vec4 voxel(vec3 pos) {
     return texture3dFrom2dNearest(tVolumeData, pos / uGridDim, uGridDim, uGridTexDim.xy);
 }
 
-vec4 voxelPadded(vec3 pos) {
+float voxelValuePadded(vec3 pos) {
     pos = min(max(vec3(0.0), pos), uGridDim - vec3(vec2(2.0), 1.0)); // remove xy padding
-    return texture3dFrom2dNearest(tVolumeData, pos / uGridDim, uGridDim, uGridTexDim.xy);
+    vec4 v = texture3dFrom2dNearest(tVolumeData, pos / uGridDim, uGridDim, uGridTexDim.xy);
+    #ifdef dValueChannel_red
+        return v.r;
+    #else
+        return v.a;
+    #endif
 }
 
 int idot2(const in ivec2 a, const in ivec2 b) {
@@ -261,8 +266,13 @@ void main(void) {
     vec4 d0 = voxel(b0);
     vec4 d1 = voxel(b1);
 
-    float v0 = d0.a;
-    float v1 = d1.a;
+    #ifdef dValueChannel_red
+        float v0 = d0.r;
+        float v1 = d1.r;
+    #else
+        float v0 = d0.a;
+        float v1 = d1.a;
+    #endif
 
     float t = (uIsoValue - v0) / (v0 - v1);
     gl_FragData[0].xyz = (uGridTransform * vec4(b0 + t * (b0 - b1), 1.0)).xyz;
@@ -286,14 +296,14 @@ void main(void) {
 
     // normals from gradients
     vec3 n0 = -normalize(vec3(
-        voxelPadded(b0 - c1).a - voxelPadded(b0 + c1).a,
-        voxelPadded(b0 - c3).a - voxelPadded(b0 + c3).a,
-        voxelPadded(b0 - c4).a - voxelPadded(b0 + c4).a
+        voxelValuePadded(b0 - c1) - voxelValuePadded(b0 + c1),
+        voxelValuePadded(b0 - c3) - voxelValuePadded(b0 + c3),
+        voxelValuePadded(b0 - c4) - voxelValuePadded(b0 + c4)
     ));
     vec3 n1 = -normalize(vec3(
-        voxelPadded(b1 - c1).a - voxelPadded(b1 + c1).a,
-        voxelPadded(b1 - c3).a - voxelPadded(b1 + c3).a,
-        voxelPadded(b1 - c4).a - voxelPadded(b1 + c4).a
+        voxelValuePadded(b1 - c1) - voxelValuePadded(b1 + c1),
+        voxelValuePadded(b1 - c3) - voxelValuePadded(b1 + c3),
+        voxelValuePadded(b1 - c4) - voxelValuePadded(b1 + c4)
     ));
     gl_FragData[2].xyz = -vec3(
         n0.x + t * (n0.x - n1.x),

src/mol-gl/webgl/texture.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2018-2023 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2024 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  * @author Gianluca Tomasello <giagitom@gmail.com>
@@ -59,14 +59,14 @@ export function getTarget(gl: GLRenderingContext, kind: TextureKind): number {
 export function getFormat(gl: GLRenderingContext, format: TextureFormat, type: TextureType): number {
     switch (format) {
         case 'alpha':
-            if (isWebGL2(gl) && type === 'float') return gl.RED;
+            if (isWebGL2(gl) && (type === 'float' || type === 'fp16')) return gl.RED;
             else if (isWebGL2(gl) && type === 'int') return gl.RED_INTEGER;
             else return gl.ALPHA;
         case 'rgb':
             if (isWebGL2(gl) && type === 'int') return gl.RGB_INTEGER;
             return gl.RGB;
         case 'rg':
-            if (isWebGL2(gl) && type === 'float') return gl.RG;
+            if (isWebGL2(gl) && (type === 'float' || type === 'fp16')) return gl.RG;
             else if (isWebGL2(gl) && type === 'int') return gl.RG_INTEGER;
             else throw new Error('texture format "rg" requires webgl2 and type "float" or int"');
         case 'rgba':

src/mol-model/custom-property.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2018-2021 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2024 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
@@ -9,6 +9,7 @@ import { CifWriter } from '../mol-io/writer/cif';
 import { CifExportContext } from './structure/export/mmcif';
 import { QuerySymbolRuntime } from '../mol-script/runtime/query/compiler';
 import { UUID } from '../mol-util';
+import { arrayRemoveInPlace } from '../mol-util/array';
 
 export { CustomPropertyDescriptor, CustomProperties };
 
@@ -61,6 +62,14 @@ class CustomProperties {
         this._set.add(desc);
     }
 
+    remove(desc: CustomPropertyDescriptor<any>) {
+        if (!this._set.has(desc)) return;
+
+        arrayRemoveInPlace(this._list, desc);
+        this._set.delete(desc);
+        this.assets(desc);
+    }
+
     reference(desc: CustomPropertyDescriptor<any>, add: boolean) {
         let refs = this._refs.get(desc) || 0;
         refs += add ? 1 : -1;

src/mol-repr/structure/visual/gaussian-density-volume.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2018-2021 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2024 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
@@ -31,7 +31,7 @@ async function createGaussianDensityVolume(ctx: VisualContext, structure: Struct
     const unitToCartn = Mat4.mul(Mat4(), transform, Mat4.fromScaling(Mat4(), gridDim));
     const cellDim = Mat4.getScaling(Vec3(), transform);
     const axisOrder = Vec3.create(0, 1, 2);
-    const vol = DirectVolume.create(bbox, gridDim, transform, unitToCartn, cellDim, texture, stats, true, axisOrder, directVolume);
+    const vol = DirectVolume.create(bbox, gridDim, transform, unitToCartn, cellDim, texture, stats, true, axisOrder, 'byte', directVolume);
 
     const sphere = Sphere3D.expand(Sphere3D(), structure.boundary.sphere, densityTextureData.maxRadius);
     vol.setBoundingSphere(sphere);
@@ -89,7 +89,7 @@ async function createUnitsGaussianDensityVolume(ctx: VisualContext, unit: Unit,
     const unitToCartn = Mat4.mul(Mat4(), transform, Mat4.fromScaling(Mat4(), gridDim));
     const cellDim = Mat4.getScaling(Vec3(), transform);
     const axisOrder = Vec3.create(0, 1, 2);
-    const vol = DirectVolume.create(bbox, gridDim, transform, unitToCartn, cellDim, texture, stats, true, axisOrder, directVolume);
+    const vol = DirectVolume.create(bbox, gridDim, transform, unitToCartn, cellDim, texture, stats, true, axisOrder, 'byte', directVolume);
 
     const sphere = Sphere3D.expand(Sphere3D(), unit.boundary.sphere, densityTextureData.maxRadius);
     vol.setBoundingSphere(sphere);