wip, gpu mc & vol isosurface

src/mol-geo/geometry/texture-mesh/texture-mesh.ts

@@ -104,7 +104,8 @@ export namespace TextureMesh {
 
         const counts = { drawCount: textureMesh.vertexCount, vertexCount: textureMesh.vertexCount / 3, groupCount, instanceCount };
 
-        const transformBoundingSphere = calculateTransformBoundingSphere(textureMesh.boundingSphere, transform.aTransform.ref.value, transform.instanceCount.ref.value);
+        const invariantBoundingSphere = Sphere3D.clone(textureMesh.boundingSphere);
+        const boundingSphere = calculateTransformBoundingSphere(invariantBoundingSphere, transform.aTransform.ref.value, instanceCount);
 
         return {
             uGeoTexDim: textureMesh.geoTextureDim,
@@ -113,9 +114,9 @@ export namespace TextureMesh {
 
             // aGroup is used as a vertex index here and the group id is retirieved from tPositionGroup
             aGroup: ValueCell.create(fillSerial(new Float32Array(textureMesh.vertexCount))),
-            boundingSphere: ValueCell.create(transformBoundingSphere),
-            invariantBoundingSphere: ValueCell.create(Sphere3D.clone(textureMesh.boundingSphere)),
-            uInvariantBoundingSphere: ValueCell.create(Vec4.ofSphere(textureMesh.boundingSphere)),
+            boundingSphere: ValueCell.create(boundingSphere),
+            invariantBoundingSphere: ValueCell.create(invariantBoundingSphere),
+            uInvariantBoundingSphere: ValueCell.create(Vec4.ofSphere(invariantBoundingSphere)),
 
             ...color,
             ...marker,
@@ -141,8 +142,7 @@ export namespace TextureMesh {
     }
 
     function updateValues(values: TextureMeshValues, props: PD.Values<Params>) {
-        ValueCell.updateIfChanged(values.alpha, props.alpha); // `uAlpha` is set in renderable.render
-
+        BaseGeometry.updateValues(values, props);
         ValueCell.updateIfChanged(values.dDoubleSided, props.doubleSided);
         ValueCell.updateIfChanged(values.dFlatShaded, props.flatShaded);
         ValueCell.updateIfChanged(values.dFlipSided, props.flipSided);
@@ -156,8 +156,9 @@ export namespace TextureMesh {
     }
 
     function updateBoundingSphere(values: TextureMeshValues, textureMesh: TextureMesh) {
-        const invariantBoundingSphere = textureMesh.boundingSphere;
+        const invariantBoundingSphere = Sphere3D.clone(textureMesh.boundingSphere);
         const boundingSphere = calculateTransformBoundingSphere(invariantBoundingSphere, values.aTransform.ref.value, values.instanceCount.ref.value);
+
         if (!Sphere3D.equals(boundingSphere, values.boundingSphere.ref.value)) {
             ValueCell.update(values.boundingSphere, boundingSphere);
         }

src/mol-gl/compute/histogram-pyramid/reduction.ts

@@ -124,11 +124,11 @@ export function createHistogramPyramid(ctx: WebGLContext, inputTexture: Texture,
     const maxSize = Math.pow(2, levels);
     // console.log('levels', levels, 'maxSize', maxSize, 'input', w);
 
-    const pyramidTexture = getTexture('pyramid', ctx, 'image-float32', 'rgba', 'float', 'nearest');
-    pyramidTexture.define(maxSize, maxSize);
+    const pyramidTex = getTexture('pyramid', ctx, 'image-float32', 'rgba', 'float', 'nearest');
+    pyramidTex.define(maxSize, maxSize);
 
     const framebuffer = getFramebuffer('pyramid', ctx);
-    pyramidTexture.attachFramebuffer(framebuffer, 0);
+    pyramidTex.attachFramebuffer(framebuffer, 0);
     gl.viewport(0, 0, maxSize, maxSize);
     gl.clear(gl.COLOR_BUFFER_BIT);
 
@@ -162,29 +162,24 @@ export function createHistogramPyramid(ctx: WebGLContext, inputTexture: Texture,
         gl.scissor(0, 0, gridTexDim[0], gridTexDim[1]);
         renderable.render();
 
-        pyramidTexture.bind(0);
+        pyramidTex.bind(0);
         gl.copyTexSubImage2D(gl.TEXTURE_2D, 0, offset, 0, 0, 0, size, size);
-        pyramidTexture.unbind(0);
+        pyramidTex.unbind(0);
 
         offset += size;
     }
 
     gl.finish();
 
-    // printTexture(ctx, pyramidTexture, 2)
+    // printTexture(ctx, pyramidTex, 2)
 
     //
 
-    const finalCount = getHistopyramidSum(ctx, levelTexturesFramebuffers[0].texture);
-    const height = Math.ceil(finalCount / Math.pow(2, levels));
+    // return at least a count of one to avoid issues downstram
+    const count = Math.max(1, getHistopyramidSum(ctx, levelTexturesFramebuffers[0].texture));
+    const height = Math.ceil(count / Math.pow(2, levels));
     // const scale = Vec2.create(maxSize / inputTexture.width, maxSize / inputTexture.height);
     // console.log('height', height, 'finalCount', finalCount, 'scale', scale);
 
-    return {
-        pyramidTex: pyramidTexture,
-        count: finalCount,
-        height,
-        levels,
-        scale
-    };
+    return { pyramidTex, count, height, levels, scale };
 }

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

@@ -109,9 +109,9 @@ export function calcActiveVoxels(ctx: WebGLContext, volumeData: Texture, gridDim
     gl.scissor(0, 0, gridTexDim[0], gridTexDim[1]);
     renderable.render();
 
-    // console.log('gridScale', gridScale, 'gridTexDim', gridTexDim, 'gridDim', gridDim)
-    // console.log('volumeData', volumeData)
-    // console.log('at', readTexture(ctx, activeVoxelsTex))
+    // console.log('gridScale', gridScale, 'gridTexDim', gridTexDim, 'gridDim', gridDim);
+    // console.log('volumeData', volumeData);
+    // console.log('at', readTexture(ctx, activeVoxelsTex));
 
     gl.finish();
 

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

@@ -11,7 +11,7 @@ uniform vec3 uGridTexDim;
 uniform vec2 uScale;
 
 // cube corners
-const vec3 c0 = vec3(0., 0., 0.);
+// const vec3 c0 = vec3(0., 0., 0.);
 const vec3 c1 = vec3(1., 0., 0.);
 const vec3 c2 = vec3(1., 1., 0.);
 const vec3 c3 = vec3(0., 1., 0.);
@@ -25,8 +25,7 @@ vec3 index3dFrom2d(vec2 coord) {
     vec2 columnRow = floor(gridTexPos / uGridDim.xy);
     vec2 posXY = gridTexPos - columnRow * uGridDim.xy;
     float posZ = columnRow.y * floor(uGridTexDim.x / uGridDim.x) + columnRow.x;
-    vec3 posXYZ = vec3(posXY, posZ) / uGridDim;
-    return posXYZ;
+    return vec3(posXY, posZ);
 }
 
 float intDiv(float a, float b) { return float(int(a) / int(b)); }
@@ -41,7 +40,7 @@ vec4 texture3dFrom2dNearest(sampler2D tex, vec3 pos, vec3 gridDim, vec2 texDim)
 }
 
 vec4 voxel(vec3 pos) {
-    return texture3dFrom2dNearest(tVolumeData, pos, uGridDim, uGridTexDim.xy);
+    return texture3dFrom2dNearest(tVolumeData, pos / uGridDim, uGridDim, uGridTexDim.xy);
 }
 
 void main(void) {
@@ -50,20 +49,26 @@ void main(void) {
 
     // 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 / uGridDim).a, uIsoValue)
-        + 4. * step(voxel(posXYZ + c2 / uGridDim).a, uIsoValue)
-        + 8. * step(voxel(posXYZ + c3 / uGridDim).a, uIsoValue)
-        + 16. * step(voxel(posXYZ + c4 / uGridDim).a, uIsoValue)
-        + 32. * step(voxel(posXYZ + c5 / uGridDim).a, uIsoValue)
-        + 64. * step(voxel(posXYZ + c6 / uGridDim).a, uIsoValue)
-        + 128. * step(voxel(posXYZ + c7 / uGridDim).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);
     c *= step(c, 254.);
 
+    // handle out of bounds positions
+    posXYZ += 1.0;
+    posXYZ.xy += 1.0; // one pixel padding (usually ok even if the texture has no padding)
+    if (posXYZ.x >= uGridDim.x || posXYZ.y >= uGridDim.y || posXYZ.z >= uGridDim.z)
+        c = 0.0;
+
     // get total triangles to generate for calculated MC case from triCount texture
     float totalTrianglesToGenerate = texture2D(tTriCount, vec2(intMod(c, 16.), floor(c / 16.)) / 16.).a;
     gl_FragColor = vec4(vec3(totalTrianglesToGenerate * 3.0), c / 255.0);
 
-    // gl_FragColor = vec4(255.0, 0.0, 0.0, voxel(posXYZ + c4 / uGridDim).a * 255.0);
+    // gl_FragColor = vec4(255.0, 0.0, 0.0, voxel(posXYZ + c4).a * 255.0);
     // gl_FragColor = vec4(255.0, 0.0, 0.0, voxel(posXYZ).a * 255.0);
 
     // vec2 uv = vCoordinate;

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

@@ -31,15 +31,12 @@ const vec3 c5 = vec3(1., 0., 1.);
 const vec3 c6 = vec3(1., 1., 1.);
 const vec3 c7 = vec3(0., 1., 1.);
 
-const float EPS = 0.00001;
-
 vec3 index3dFrom2d(vec2 coord) {
     vec2 gridTexPos = coord * uGridTexDim.xy;
     vec2 columnRow = floor(gridTexPos / uGridDim.xy);
     vec2 posXY = gridTexPos - columnRow * uGridDim.xy;
     float posZ = columnRow.y * floor(uGridTexDim.x / uGridDim.x) + columnRow.x;
-    vec3 posXYZ = vec3(posXY, posZ);
-    return posXYZ;
+    return vec3(posXY, posZ);
 }
 
 vec4 texture3dFrom2dNearest(sampler2D tex, vec3 pos, vec3 gridDim, vec2 texDim) {
@@ -51,7 +48,7 @@ vec4 texture3dFrom2dNearest(sampler2D tex, vec3 pos, vec3 gridDim, vec2 texDim)
 }
 
 vec4 voxel(vec3 pos) {
-    return texture3dFrom2dNearest(tVolumeData, pos, uGridDim, uGridTexDim.xy);
+    return texture3dFrom2dNearest(tVolumeData, pos / uGridDim, uGridDim, uGridTexDim.xy);
 }
 
 void main(void) {
@@ -173,27 +170,34 @@ void main(void) {
     // b0 = floor(b0 + 0.5);
     // b1 = floor(b1 + 0.5);
 
-    vec4 d0 = voxel(b0 / uGridDim);
-    vec4 d1 = voxel(b1 / uGridDim);
+    vec4 d0 = voxel(b0);
+    vec4 d1 = voxel(b1);
 
     float v0 = d0.a;
     float v1 = d1.a;
 
     float t = (uIsoValue - v0) / (v0 - v1);
-    // t = -0.5;
     gl_FragData[0].xyz = (uGridTransform * vec4(b0 + t * (b0 - b1), 1.0)).xyz;
-    gl_FragData[0].w = decodeFloatRGB(d0.rgb); // group id
+
+    // group id
+    #if __VERSION__ == 100
+        // webgl1 does not support 'flat' interpolation (i.e. no interpolation)
+        // so we provide a constant that stands for 'unknown-group'
+        gl_FragData[0].w = 16777216.0;
+    #else
+        gl_FragData[0].w = t < 0.5 ? decodeFloatRGB(d0.rgb) : decodeFloatRGB(d1.rgb);
+    #endif
 
     // normals from gradients
     vec3 n0 = -normalize(vec3(
-        voxel((b0 - c1) / uGridDim).a - voxel((b0 + c1) / uGridDim).a,
-        voxel((b0 - c3) / uGridDim).a - voxel((b0 + c3) / uGridDim).a,
-        voxel((b0 - c4) / uGridDim).a - voxel((b0 + c4) / uGridDim).a
+        voxel(b0 - c1).a - voxel(b0 + c1).a,
+        voxel(b0 - c3).a - voxel(b0 + c3).a,
+        voxel(b0 - c4).a - voxel(b0 + c4).a
     ));
     vec3 n1 = -normalize(vec3(
-        voxel((b1 - c1) / uGridDim).a - voxel((b1 + c1) / uGridDim).a,
-        voxel((b1 - c3) / uGridDim).a - voxel((b1 + c3) / uGridDim).a,
-        voxel((b1 - c4) / uGridDim).a - voxel((b1 + c4) / uGridDim).a
+        voxel(b1 - c1).a - voxel(b1 + c1).a,
+        voxel(b1 - c3).a - voxel(b1 + c3).a,
+        voxel(b1 - c4).a - voxel(b1 + c4).a
     ));
     gl_FragData[1].xyz = -vec3(
         n0.x + t * (n0.x - n1.x),

src/mol-gl/webgl/program.ts

@@ -24,7 +24,7 @@ export interface Program {
     use: () => void
     setUniforms: (uniformValues: UniformsList) => void
     bindAttributes: (attribueBuffers: AttributeBuffers) => void
-    bindTextures: (textures: Textures, startingTargetUnit?: number) => void
+    bindTextures: (textures: Textures, startingTargetUnit: number) => void
 
     reset: () => void
     destroy: () => void
@@ -198,9 +198,7 @@ export function createProgram(gl: GLRenderingContext, state: WebGLState, extensi
                 if (l !== -1) buffer.bind(l);
             }
         },
-        bindTextures: (textures: Textures, startingTargetUnit?: number) => {
-            startingTargetUnit = startingTargetUnit ?? 0;
-
+        bindTextures: (textures: Textures, startingTargetUnit: number) => {
             for (let i = 0, il = textures.length; i < il; ++i) {
                 const [k, texture] = textures[i];
                 const l = locations[k];

src/mol-gl/webgl/render-item.ts

@@ -173,7 +173,7 @@ export function createRenderItem<T extends string>(ctx: WebGLContext, drawMode:
                     program.bindTextures(sharedTexturesList, 0);
                     program.bindTextures(textures, sharedTexturesList.length);
                 } else {
-                    program.bindTextures(textures);
+                    program.bindTextures(textures, 0);
                 }
             } else {
                 const vertexArray = vertexArrays[variant];
@@ -191,7 +191,7 @@ export function createRenderItem<T extends string>(ctx: WebGLContext, drawMode:
                     program.bindTextures(sharedTexturesList, 0);
                     program.bindTextures(textures, sharedTexturesList.length);
                 } else {
-                    program.bindTextures(textures);
+                    program.bindTextures(textures, 0);
                 }
                 if (vertexArray) {
                     vertexArray.bind();

src/mol-gl/webgl/texture.ts

@@ -173,7 +173,11 @@ export interface Texture {
     getByteCount: () => number
 
     define: (width: number, height: number, depth?: number) => void
-    load: (image: TextureImage<any> | TextureVolume<any>) => void
+    /**
+     * The `sub` option requires an existing allocation on the GPU, that is, either
+     * `define` or `load` without `sub` must have been called before.
+     */
+    load: (image: TextureImage<any> | TextureVolume<any>, sub?: boolean) => void
     bind: (id: TextureId) => void
     unbind: (id: TextureId) => void
     /** Use `layer` to attach a z-slice of a 3D texture */
@@ -246,7 +250,7 @@ export function createTexture(gl: GLRenderingContext, extensions: WebGLExtension
         }
     }
 
-    function load(data: TextureImage<any> | TextureVolume<any>) {
+    function load(data: TextureImage<any> | TextureVolume<any>, sub = false) {
         gl.bindTexture(target, texture);
         // unpack alignment of 1 since we use textures only for data
         gl.pixelStorei(gl.UNPACK_ALIGNMENT, 1);
@@ -254,12 +258,20 @@ export function createTexture(gl: GLRenderingContext, extensions: WebGLExtension
         gl.pixelStorei(gl.UNPACK_PREMULTIPLY_ALPHA_WEBGL, 0);
         if (isTexture2d(data, target, gl)) {
             gl.pixelStorei(gl.UNPACK_FLIP_Y_WEBGL, !!data.flipY);
-            width = data.width, height = data.height;
-            gl.texImage2D(target, 0, internalFormat, width, height, 0, format, type, data.array);
+            if (sub) {
+                gl.texSubImage2D(target, 0, 0, 0, data.width, data.height, format, type, data.array);
+            } else {
+                width = data.width, height = data.height;
+                gl.texImage2D(target, 0, internalFormat, width, height, 0, format, type, data.array);
+            }
         } else if (isWebGL2(gl) && isTexture3d(data, target, gl)) {
             gl.pixelStorei(gl.UNPACK_FLIP_Y_WEBGL, false);
-            width = data.width, height = data.height, depth = data.depth;
-            gl.texImage3D(target, 0, internalFormat, width, height, depth, 0, format, type, data.array);
+            if (sub) {
+                gl.texSubImage3D(target, 0, 0, 0, 0, data.width, data.height, data.depth, format, type, data.array);
+            } else {
+                width = data.width, height = data.height, depth = data.depth;
+                gl.texImage3D(target, 0, internalFormat, width, height, depth, 0, format, type, data.array);
+            }
         } else {
             throw new Error('unknown texture target');
         }
@@ -317,11 +329,8 @@ export function createTexture(gl: GLRenderingContext, extensions: WebGLExtension
             texture = getTexture(gl);
             init();
 
-            if (loadedData) {
-                load(loadedData);
-            } else {
-                define(width, height, depth);
-            }
+            define(width, height, depth);
+            if (loadedData) load(loadedData);
         },
         destroy: () => {
             if (destroyed) return;

src/mol-repr/structure/visual/gaussian-surface-mesh.ts

@@ -157,7 +157,7 @@ async function createGaussianSurfaceTextureMesh(ctx: VisualContext, unit: Unit,
     // ctx.webgl.waitForGpuCommandsCompleteSync();
     // console.timeEnd('createIsosurfaceBuffers');
 
-    const boundingSphere = Sphere3D.expand(Sphere3D(), structure.boundary.sphere, props.radiusOffset + getStructureExtraRadius(structure));
+    const boundingSphere = Sphere3D.expand(Sphere3D(), unit.boundary.sphere, props.radiusOffset + getStructureExtraRadius(structure));
     const surface = TextureMesh.create(gv.vertexCount, 1, gv.vertexGroupTexture, gv.normalTexture, boundingSphere, textureMesh);
     // console.log({
     //     renderables: ctx.webgl.namedComputeRenderables,

src/mol-repr/volume/direct-volume.ts

@@ -21,15 +21,7 @@ import { RepresentationContext, RepresentationParamsGetter } from '../representa
 import { Interval } from '../../mol-data/int';
 import { Loci, EmptyLoci } from '../../mol-model/loci';
 import { PickingId } from '../../mol-geo/geometry/picking';
-import { eachVolumeLoci } from './util';
-
-// avoiding namespace lookup improved performance in Chrome (Aug 2020)
-const v3set = Vec3.set;
-const v3normalize = Vec3.normalize;
-const v3sub = Vec3.sub;
-const v3addScalar = Vec3.addScalar;
-const v3scale = Vec3.scale;
-const v3toArray = Vec3.toArray;
+import { createVolumeTexture2d, createVolumeTexture3d, eachVolumeLoci } from './util';
 
 function getBoundingBox(gridDimension: Vec3, transform: Mat4) {
     const bbox = Box3D();
@@ -40,75 +32,9 @@ function getBoundingBox(gridDimension: Vec3, transform: Mat4) {
 
 // 2d volume texture
 
-function getVolumeTexture2dLayout(dim: Vec3, maxTextureSize: number) {
-    let width = 0;
-    let height = dim[1];
-    let rows = 1;
-    let columns = dim[0];
-    if (maxTextureSize < dim[0] * dim[2]) {
-        columns =  Math.floor(maxTextureSize / dim[0]);
-        rows = Math.ceil(dim[2] / columns);
-        width = columns * dim[0];
-        height *= rows;
-    } else {
-        width = dim[0] * dim[2];
-    }
-    return { width, height, columns, rows };
-}
-
-function createVolumeTexture2d(volume: Volume, maxTextureSize: number) {
-    const { cells: { space, data }, stats: { max, min } } = volume.grid;
-    const dim = space.dimensions as Vec3;
-    const { dataOffset: o } = space;
-    const { width, height } = getVolumeTexture2dLayout(dim, maxTextureSize);
-
-    const array = new Uint8Array(width * height * 4);
-    const textureImage = { array, width, height };
-
-    const diff = max - min;
-    const [ xn, yn, zn ] = dim;
-
-    const n0 = Vec3();
-    const n1 = Vec3();
-
-    const xn1 = xn - 1;
-    const yn1 = yn - 1;
-    const zn1 = zn - 1;
-
-    for (let z = 0; z < zn; ++z) {
-        for (let y = 0; y < yn; ++y) {
-            for (let x = 0; x < xn; ++x) {
-                const column = Math.floor(((z * xn) % width) / xn);
-                const row = Math.floor((z * xn) / width);
-                const px = column * xn + x;
-                const index = 4 * ((row * yn * width) + (y * width) + px);
-                const offset = o(x, y, z);
-
-                v3set(n0,
-                    data[o(Math.max(0, x - 1), y, z)],
-                    data[o(x, Math.max(0, y - 1), z)],
-                    data[o(x, y, Math.max(0, z - 1))]
-                );
-                v3set(n1,
-                    data[o(Math.min(xn1, x + 1), y, z)],
-                    data[o(x, Math.min(yn1, y + 1), z)],
-                    data[o(x, y, Math.min(zn1, z + 1))]
-                );
-                v3normalize(n0, v3sub(n0, n0, n1));
-                v3addScalar(n0, v3scale(n0, n0, 0.5), 0.5);
-                v3toArray(v3scale(n0, n0, 255), array, index);
-
-                array[index + 3] = ((data[offset] - min) / diff) * 255;
-            }
-        }
-    }
-
-    return textureImage;
-}
-
 export function createDirectVolume2d(ctx: RuntimeContext, webgl: WebGLContext, volume: Volume, directVolume?: DirectVolume) {
     const gridDimension = volume.grid.cells.space.dimensions as Vec3;
-    const textureImage = createVolumeTexture2d(volume, webgl.maxTextureSize);
+    const textureImage = createVolumeTexture2d(volume, 'normals');
     // debugTexture(createImageData(textureImage.array, textureImage.width, textureImage.height), 1/3)
     const transform = Grid.getGridToCartesianTransform(volume.grid);
     const bbox = getBoundingBox(gridDimension, transform);
@@ -123,51 +49,6 @@ export function createDirectVolume2d(ctx: RuntimeContext, webgl: WebGLContext, v
 
 // 3d volume texture
 
-function createVolumeTexture3d(volume: Volume) {
-    const { cells: { space, data }, stats: { max, min } } = volume.grid;
-    const [ width, height, depth ] = space.dimensions as Vec3;
-    const { dataOffset: o } = space;
-
-    const array = new Uint8Array(width * height * depth * 4);
-    const textureVolume = { array, width, height, depth };
-    const diff = max - min;
-
-    const n0 = Vec3();
-    const n1 = Vec3();
-
-    const width1 = width - 1;
-    const height1 = height - 1;
-    const depth1 = depth - 1;
-
-    let i = 0;
-    for (let z = 0; z < depth; ++z) {
-        for (let y = 0; y < height; ++y) {
-            for (let x = 0; x < width; ++x) {
-                const offset = o(x, y, z);
-
-                v3set(n0,
-                    data[o(Math.max(0, x - 1), y, z)],
-                    data[o(x, Math.max(0, y - 1), z)],
-                    data[o(x, y, Math.max(0, z - 1))]
-                );
-                v3set(n1,
-                    data[o(Math.min(width1, x + 1), y, z)],
-                    data[o(x, Math.min(height1, y + 1), z)],
-                    data[o(x, y, Math.min(depth1, z + 1))]
-                );
-                v3normalize(n0, v3sub(n0, n0, n1));
-                v3addScalar(n0, v3scale(n0, n0, 0.5), 0.5);
-                v3toArray(v3scale(n0, n0, 255), array, i);
-
-                array[i + 3] = ((data[offset] - min) / diff) * 255;
-                i += 4;
-            }
-        }
-    }
-
-    return textureVolume;
-}
-
 function getUnitToCartn(grid: Grid) {
     if (grid.transform.kind === 'matrix') {
         return {

src/mol-repr/volume/isosurface.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2018-2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2021 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>
@@ -20,9 +20,13 @@ import { RepresentationContext, RepresentationParamsGetter, Representation } fro
 import { PickingId } from '../../mol-geo/geometry/picking';
 import { EmptyLoci, Loci } from '../../mol-model/loci';
 import { Interval } from '../../mol-data/int';
-import { Tensor } from '../../mol-math/linear-algebra';
+import { Tensor, Vec2, Vec3 } from '../../mol-math/linear-algebra';
 import { fillSerial } from '../../mol-util/array';
-import { eachVolumeLoci } from './util';
+import { createVolumeTexture2d, eachVolumeLoci, getVolumeTexture2dLayout } from './util';
+import { TextureMesh } from '../../mol-geo/geometry/texture-mesh/texture-mesh';
+import { calcActiveVoxels } from '../../mol-gl/compute/marching-cubes/active-voxels';
+import { createHistogramPyramid } from '../../mol-gl/compute/histogram-pyramid/reduction';
+import { createIsosurfaceBuffers } from '../../mol-gl/compute/marching-cubes/isosurface';
 
 export const VolumeIsosurfaceParams = {
     isoValue: Volume.IsoValueParam
@@ -91,6 +95,80 @@ export function IsosurfaceMeshVisual(materialId: number): VolumeVisual<Isosurfac
 
 //
 
+async function createVolumeIsosurfaceTextureMesh(ctx: VisualContext, volume: Volume, theme: Theme, props: VolumeIsosurfaceProps, textureMesh?: TextureMesh) {
+    if (!ctx.webgl) throw new Error('webgl context required to create volume isosurface texture-mesh');
+
+    const { resources } = ctx.webgl;
+    if (!volume._propertyData['texture2d']) {
+        // TODO: handle disposal
+        volume._propertyData['texture2d'] = resources.texture('image-uint8', 'rgba', 'ubyte', 'linear');
+    }
+    const texture = volume._propertyData['texture2d'];
+
+    const padding = 1;
+    const transform = Grid.getGridToCartesianTransform(volume.grid);
+    const gridDimension = Vec3.clone(volume.grid.cells.space.dimensions as Vec3);
+    const { width, height, powerOfTwoSize: texDim } = getVolumeTexture2dLayout(gridDimension, padding);
+    const gridTexDim = Vec3.create(width, height, 0);
+    const gridTexScale = Vec2.create(width / texDim, height / texDim);
+    // console.log({ texDim, width, height, gridDimension });
+
+    if (!textureMesh) {
+        // set to power-of-two size required for histopyramid calculation
+        texture.define(texDim, texDim);
+        // load volume into sub-section of texture
+        texture.load(createVolumeTexture2d(volume, 'groups', padding), true);
+    }
+
+    const { max, min } = volume.grid.stats;
+    const diff = max - min;
+    const value = Volume.IsoValue.toAbsolute(props.isoValue, volume.grid.stats).absoluteValue;
+    const isoLevel = ((value - min) / diff);
+
+    gridDimension[0] += padding;
+    gridDimension[1] += padding;
+
+    // console.time('calcActiveVoxels');
+    const activeVoxelsTex = calcActiveVoxels(ctx.webgl, texture, gridDimension, gridTexDim, isoLevel, gridTexScale);
+    // ctx.webgl.waitForGpuCommandsCompleteSync();
+    // console.timeEnd('calcActiveVoxels');
+
+    // console.time('createHistogramPyramid');
+    const compacted = createHistogramPyramid(ctx.webgl, activeVoxelsTex, gridTexScale, gridTexDim);
+    // ctx.webgl.waitForGpuCommandsCompleteSync();
+    // console.timeEnd('createHistogramPyramid');
+
+    // console.time('createIsosurfaceBuffers');
+    const gv = createIsosurfaceBuffers(ctx.webgl, activeVoxelsTex, texture, compacted, gridDimension, gridTexDim, transform, isoLevel, textureMesh ? textureMesh.vertexGroupTexture.ref.value : undefined, textureMesh ? textureMesh.normalTexture.ref.value : undefined);
+    // ctx.webgl.waitForGpuCommandsCompleteSync();
+    // console.timeEnd('createIsosurfaceBuffers');
+
+    const surface = TextureMesh.create(gv.vertexCount, 1, gv.vertexGroupTexture, gv.normalTexture, Volume.getBoundingSphere(volume), textureMesh);
+    // console.log({
+    //     renderables: ctx.webgl.namedComputeRenderables,
+    //     framebuffers: ctx.webgl.namedFramebuffers,
+    //     textures: ctx.webgl.namedTextures,
+    // });
+    // ctx.webgl.waitForGpuCommandsCompleteSync();
+    return surface;
+}
+
+export function IsosurfaceTextureMeshVisual(materialId: number): VolumeVisual<IsosurfaceMeshParams> {
+    return VolumeVisual<TextureMesh, IsosurfaceMeshParams>({
+        defaultProps: PD.getDefaultValues(IsosurfaceMeshParams),
+        createGeometry: createVolumeIsosurfaceTextureMesh,
+        createLocationIterator: (volume: Volume) => LocationIterator(volume.grid.cells.data.length, 1, 1, () => NullLocation),
+        getLoci: getIsosurfaceLoci,
+        eachLocation: eachIsosurface,
+        setUpdateState: (state: VisualUpdateState, volume: Volume, newProps: PD.Values<IsosurfaceMeshParams>, currentProps: PD.Values<IsosurfaceMeshParams>) => {
+            if (!Volume.IsoValue.areSame(newProps.isoValue, currentProps.isoValue, volume.grid.stats)) state.createGeometry = true;
+        },
+        geometryUtils: TextureMesh.Utils
+    }, materialId);
+}
+
+//
+
 export async function createVolumeIsosurfaceWireframe(ctx: VisualContext, volume: Volume, theme: Theme, props: VolumeIsosurfaceProps, lines?: Lines) {
     ctx.runtime.update({ message: 'Marching cubes...' });
 
@@ -136,6 +214,8 @@ export function IsosurfaceWireframeVisual(materialId: number): VolumeVisual<Isos
 
 const IsosurfaceVisuals = {
     'solid': (ctx: RepresentationContext, getParams: RepresentationParamsGetter<Volume, IsosurfaceMeshParams>) => VolumeRepresentation('Isosurface mesh', ctx, getParams, IsosurfaceMeshVisual, getLoci),
+    // TODO: don't enable yet as it breaks state sessions
+    // 'solid-gpu': (ctx: RepresentationContext, getParams: RepresentationParamsGetter<Volume, IsosurfaceMeshParams>) => VolumeRepresentation('Isosurface texture-mesh', ctx, getParams, IsosurfaceTextureMeshVisual, getLoci),
     'wireframe': (ctx: RepresentationContext, getParams: RepresentationParamsGetter<Volume, IsosurfaceWireframeParams>) => VolumeRepresentation('Isosurface wireframe', ctx, getParams, IsosurfaceWireframeVisual, getLoci),
 };
 

src/mol-repr/volume/util.ts

@@ -8,6 +8,16 @@ import { Volume } from '../../mol-model/volume';
 import { Loci } from '../../mol-model/loci';
 import { Interval, OrderedSet } from '../../mol-data/int';
 import { equalEps } from '../../mol-math/linear-algebra/3d/common';
+import Vec3 from '../../mol-math/linear-algebra/3d/vec3';
+import { encodeFloatRGBtoArray } from '../../mol-util/float-packing';
+
+// avoiding namespace lookup improved performance in Chrome (Aug 2020)
+const v3set = Vec3.set;
+const v3normalize = Vec3.normalize;
+const v3sub = Vec3.sub;
+const v3addScalar = Vec3.addScalar;
+const v3scale = Vec3.scale;
+const v3toArray = Vec3.toArray;
 
 export function eachVolumeLoci(loci: Loci, volume: Volume, isoValue: Volume.IsoValue | undefined, apply: (interval: Interval) => boolean) {
     let changed = false;
@@ -41,4 +51,127 @@ export function eachVolumeLoci(loci: Loci, volume: Volume, isoValue: Volume.IsoV
         }
     }
     return changed;
+}
+
+//
+
+export function getVolumeTexture2dLayout(dim: Vec3, padding = 0) {
+    const area = dim[0] * dim[1] * dim[2];
+    const squareDim = Math.sqrt(area);
+    const powerOfTwoSize = Math.pow(2, Math.ceil(Math.log(squareDim) / Math.log(2)));
+
+    let width = dim[0] + padding;
+    let height = dim[1] + padding;
+    let rows = 1;
+    let columns = width;
+    if (powerOfTwoSize < width * dim[2]) {
+        columns =  Math.floor(powerOfTwoSize / width);
+        rows = Math.ceil(dim[2] / columns);
+        width *= columns;
+        height *= rows;
+    } else {
+        width *= dim[2];
+    }
+    return { width, height, columns, rows, powerOfTwoSize: height < powerOfTwoSize ? powerOfTwoSize : powerOfTwoSize * 2 };
+}
+
+export function createVolumeTexture2d(volume: Volume, variant: 'normals' | 'groups', padding = 0) {
+    const { cells: { space, data }, stats: { max, min } } = volume.grid;
+    const dim = space.dimensions as Vec3;
+    const { dataOffset: o } = space;
+    const { width, height } = getVolumeTexture2dLayout(dim, padding);
+
+    const array = new Uint8Array(width * height * 4);
+    const textureImage = { array, width, height };
+
+    const diff = max - min;
+    const [ xn, yn, zn ] = dim;
+    const xnp = xn + padding;
+    const ynp = yn + padding;
+
+    const n0 = Vec3();
+    const n1 = Vec3();
+
+    const xn1 = xn - 1;
+    const yn1 = yn - 1;
+    const zn1 = zn - 1;
+
+    for (let z = 0; z < zn; ++z) {
+        for (let y = 0; y < yn; ++y) {
+            for (let x = 0; x < xn; ++x) {
+                const column = Math.floor(((z * xnp) % width) / xnp);
+                const row = Math.floor((z * xnp) / width);
+                const px = column * xnp + x;
+                const index = 4 * ((row * ynp * width) + (y * width) + px);
+                const offset = o(x, y, z);
+
+                if (variant === 'groups') {
+                    encodeFloatRGBtoArray(offset, array, index);
+                } else {
+                    v3set(n0,
+                        data[o(Math.max(0, x - 1), y, z)],
+                        data[o(x, Math.max(0, y - 1), z)],
+                        data[o(x, y, Math.max(0, z - 1))]
+                    );
+                    v3set(n1,
+                        data[o(Math.min(xn1, x + 1), y, z)],
+                        data[o(x, Math.min(yn1, y + 1), z)],
+                        data[o(x, y, Math.min(zn1, z + 1))]
+                    );
+                    v3normalize(n0, v3sub(n0, n0, n1));
+                    v3addScalar(n0, v3scale(n0, n0, 0.5), 0.5);
+                    v3toArray(v3scale(n0, n0, 255), array, index);
+                }
+
+                array[index + 3] = ((data[offset] - min) / diff) * 255;
+            }
+        }
+    }
+
+    return textureImage;
+}
+
+export function createVolumeTexture3d(volume: Volume) {
+    const { cells: { space, data }, stats: { max, min } } = volume.grid;
+    const [ width, height, depth ] = space.dimensions as Vec3;
+    const { dataOffset: o } = space;
+
+    const array = new Uint8Array(width * height * depth * 4);
+    const textureVolume = { array, width, height, depth };
+    const diff = max - min;
+
+    const n0 = Vec3();
+    const n1 = Vec3();
+
+    const width1 = width - 1;
+    const height1 = height - 1;
+    const depth1 = depth - 1;
+
+    let i = 0;
+    for (let z = 0; z < depth; ++z) {
+        for (let y = 0; y < height; ++y) {
+            for (let x = 0; x < width; ++x) {
+                const offset = o(x, y, z);
+
+                v3set(n0,
+                    data[o(Math.max(0, x - 1), y, z)],
+                    data[o(x, Math.max(0, y - 1), z)],
+                    data[o(x, y, Math.max(0, z - 1))]
+                );
+                v3set(n1,
+                    data[o(Math.min(width1, x + 1), y, z)],
+                    data[o(x, Math.min(height1, y + 1), z)],
+                    data[o(x, y, Math.min(depth1, z + 1))]
+                );
+                v3normalize(n0, v3sub(n0, n0, n1));
+                v3addScalar(n0, v3scale(n0, n0, 0.5), 0.5);
+                v3toArray(v3scale(n0, n0, 255), array, i);
+
+                array[i + 3] = ((data[offset] - min) / diff) * 255;
+                i += 4;
+            }
+        }
+    }
+
+    return textureVolume;
 }