marching-cubes and histogram-pyramid compute shader

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

+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { createComputeRenderable } from '../../renderable'
+import { WebGLContext } from '../../webgl/context';
+import { createComputeRenderItem } from '../../webgl/render-item';
+import { AttributeSpec, Values, TextureSpec, ValueSpec, UniformSpec } from '../../renderable/schema';
+import { Texture, createTexture } from 'mol-gl/webgl/texture';
+import { ShaderCode } from 'mol-gl/shader-code';
+import { ValueCell } from 'mol-util';
+import { GLRenderingContext } from 'mol-gl/webgl/compat';
+import { QuadPositions } from '../util';
+import { Vec2 } from 'mol-math/linear-algebra';
+import { getHistopyramidSum } from './sum';
+
+const HistopyramidReductionSchema = {
+    drawCount: ValueSpec('number'),
+    instanceCount: ValueSpec('number'),
+    aPosition: AttributeSpec('float32', 2, 0),
+
+    tPreviousLevel: TextureSpec('texture', 'rgba', 'float', 'nearest'),
+    uSize: UniformSpec('f'),
+}
+
+function getHistopyramidReductionRenderable(ctx: WebGLContext, initialTexture: Texture) {
+    const values: Values<typeof HistopyramidReductionSchema> = {
+        drawCount: ValueCell.create(6),
+        instanceCount: ValueCell.create(1),
+        aPosition: ValueCell.create(QuadPositions),
+
+        tPreviousLevel: ValueCell.create(initialTexture),
+        uSize: ValueCell.create(0),
+    }
+
+    const schema = { ...HistopyramidReductionSchema }
+    const shaderCode = ShaderCode(
+        require('mol-gl/shader/quad.vert').default,
+        require('mol-gl/shader/histogram-pyramid/reduction.frag').default,
+        { standardDerivatives: false, fragDepth: false }
+    )
+    const renderItem = createComputeRenderItem(ctx, 'triangles', shaderCode, schema, values)
+
+    return createComputeRenderable(renderItem, values);
+}
+
+/** name for shared framebuffer used for histogram-pyramid operations */
+const FramebufferName = 'histogram-pyramid'
+
+function setRenderingDefaults(gl: GLRenderingContext) {
+    gl.disable(gl.CULL_FACE)
+    gl.disable(gl.BLEND)
+    gl.disable(gl.DEPTH_TEST)
+    gl.depthMask(false)
+}
+
+export interface HistogramPyramid {
+    pyramidTex: Texture
+    totalTex: Texture
+    initialTex: Texture
+    count: number
+    height: number
+    levels: number
+    scale: Vec2
+}
+
+export function createHistogramPyramid(ctx: WebGLContext, inputTexture: Texture): HistogramPyramid {
+    const { gl, framebufferCache } = ctx
+
+    const inputTextureMaxDim = Math.max(inputTexture.width, inputTexture.height)
+
+    // This part set the levels
+    const levels = Math.ceil(Math.log(inputTextureMaxDim) / Math.log(2))
+
+    const initialTexture = createTexture(ctx, 'image-float32', 'rgba', 'float', 'nearest')
+    initialTexture.load({ array: new Float32Array(4), width: 1, height: 1 })
+    initialTexture.define(Math.pow(2, levels), Math.pow(2, levels))
+
+    const framebuffer = framebufferCache.get(FramebufferName).value
+    inputTexture.attachFramebuffer(framebuffer, 0)
+    initialTexture.define(Math.pow(2, levels), Math.pow(2, levels))
+    // TODO need to initialize texSubImage2D to make Firefox happy
+    gl.copyTexSubImage2D(gl.TEXTURE_2D, 0, 0, 0, 0, 0, inputTexture.width, inputTexture.height);
+
+    const initialTextureMaxDim = Math.max(initialTexture.width, initialTexture.height)
+
+    const pyramidTexture = createTexture(ctx, 'image-float32', 'rgba', 'float', 'nearest')
+    pyramidTexture.define(Math.pow(2, levels), Math.pow(2, levels))
+
+    const levelTextures: Texture[] = []
+    for (let i = 0; i < levels; ++i) {
+        const tex = createTexture(ctx, 'image-float32', 'rgba', 'float', 'nearest')
+        tex.define(Math.pow(2, i), Math.pow(2, i))
+        levelTextures.push(tex)
+    }
+
+    const renderable = getHistopyramidReductionRenderable(ctx, initialTexture)
+    renderable.update()
+    renderable.use()
+
+    let offset = 0;
+    for (let i = 0; i < levels; i++) {
+        const currLevel = levels - 1 - i
+        levelTextures[currLevel].attachFramebuffer(framebuffer, 0)
+
+        const size = Math.pow(2, currLevel)
+        // console.log('size', size, 'draw-level', currLevel, 'read-level', levels - i)
+        gl.clear(gl.COLOR_BUFFER_BIT)
+
+        ValueCell.update(renderable.values.uSize, Math.pow(2, i + 1) / initialTextureMaxDim)
+        const readTex = i === 0 ? initialTexture : levelTextures[levels - i]
+        // console.log(readTex.width, readTex.height)
+        ValueCell.update(renderable.values.tPreviousLevel, readTex)
+
+        renderable.update()
+        renderable.use()
+        setRenderingDefaults(gl)
+        gl.viewport(0, 0, size, size)
+        renderable.render()
+
+        pyramidTexture.bind(0)
+        // TODO need to initialize texSubImage2D to make Firefox happy
+        gl.copyTexSubImage2D(gl.TEXTURE_2D, 0, offset, 0, 0, 0, size, size);
+        pyramidTexture.unbind(0)
+
+        // if (i >= levels - 4) {
+        //     console.log('==============', i)
+        //     const rt = readTexture(ctx, levelTextures[currLevel])
+        //     console.log('array', rt.array)
+        //     for (let i = 0, il = rt.width * rt.height; i < il; ++i) {
+        //         // const v = decodeFloatRGB(rt.array[i * 4], rt.array[i * 4 + 1], rt.array[i * 4 + 2])
+        //         // console.log(i, 'v', v, 'rgb', rt.array[i * 4], rt.array[i * 4 + 1], rt.array[i * 4 + 2])
+        //         console.log(i, 'f', rt.array[i * 4])
+        //     }
+        // }
+
+        offset += size;
+    }
+
+    // printTexture(ctx, pyramidTexture, 3)
+
+    //
+
+    const finalCount = getHistopyramidSum(ctx, levelTextures[0])
+    const height = Math.ceil(finalCount / Math.pow(2, levels))
+    // console.log('height', height)
+
+    //
+
+    const scale = Vec2.create(
+        initialTexture.width / inputTexture.width,
+        initialTexture.height / inputTexture.height
+    )
+    return {
+        pyramidTex: pyramidTexture,
+        totalTex: levelTextures[0],
+        initialTex: initialTexture,
+        count: finalCount,
+        height,
+        levels,
+        scale
+    }
+}
\ No newline at end of file

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

+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { createComputeRenderable } from '../../renderable'
+import { WebGLContext } from '../../webgl/context';
+import { createComputeRenderItem } from '../../webgl/render-item';
+import { AttributeSpec, Values, TextureSpec, ValueSpec } from '../../renderable/schema';
+import { Texture, createTexture } from 'mol-gl/webgl/texture';
+import { ShaderCode } from 'mol-gl/shader-code';
+import { ValueCell } from 'mol-util';
+import { decodeFloatRGB } from 'mol-util/float-packing';
+import { QuadPositions, readTexture } from '../util';
+
+const HistopyramidSumSchema = {
+    drawCount: ValueSpec('number'),
+    instanceCount: ValueSpec('number'),
+    aPosition: AttributeSpec('float32', 2, 0),
+
+    tTexture: TextureSpec('texture', 'rgba', 'float', 'nearest'),
+}
+
+function getHistopyramidSumRenderable(ctx: WebGLContext, texture: Texture) {
+    const values: Values<typeof HistopyramidSumSchema> = {
+        drawCount: ValueCell.create(6),
+        instanceCount: ValueCell.create(1),
+        aPosition: ValueCell.create(QuadPositions),
+
+        tTexture: ValueCell.create(texture),
+    }
+
+    const schema = { ...HistopyramidSumSchema }
+    const shaderCode = ShaderCode(
+        require('mol-gl/shader/quad.vert').default,
+        require('mol-gl/shader/histogram-pyramid/sum.frag').default,
+        { standardDerivatives: false, fragDepth: false }
+    )
+    const renderItem = createComputeRenderItem(ctx, 'triangles', shaderCode, schema, values)
+
+    return createComputeRenderable(renderItem, values);
+}
+
+/** name for shared framebuffer used for histogram-pyramid operations */
+const FramebufferName = 'histogram-pyramid'
+
+export function getHistopyramidSum(ctx: WebGLContext, pyramidTopTexture: Texture) {
+    const { gl, framebufferCache } = ctx
+
+    const framebuffer = framebufferCache.get(FramebufferName).value
+
+    const encodeFloatRenderable = getHistopyramidSumRenderable(ctx, pyramidTopTexture)
+    encodeFloatRenderable.update()
+    encodeFloatRenderable.use()
+
+    const encodedFloatTexture = createTexture(ctx, 'image-uint8', 'rgba', 'ubyte', 'nearest')
+    encodedFloatTexture.define(1, 1)
+    encodedFloatTexture.attachFramebuffer(framebuffer, 0)
+
+    gl.viewport(0, 0, 1, 1)
+    encodeFloatRenderable.render()
+
+    const sumImage = readTexture(ctx, encodedFloatTexture)
+    return decodeFloatRGB(sumImage.array[0], sumImage.array[1], sumImage.array[2])
+}
\ No newline at end of file

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

+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { createComputeRenderable } from '../../renderable'
+import { WebGLContext } from '../../webgl/context';
+import { createComputeRenderItem } from '../../webgl/render-item';
+import { AttributeSpec, Values, TextureSpec, ValueSpec, UniformSpec } from '../../renderable/schema';
+import { Texture, createTexture } from 'mol-gl/webgl/texture';
+import { ShaderCode } from 'mol-gl/shader-code';
+import { ValueCell } from 'mol-util';
+import { GLRenderingContext } from 'mol-gl/webgl/compat';
+import { Vec3 } from 'mol-math/linear-algebra';
+import { QuadPositions } from '../util';
+import { getTriCount } from './tables';
+
+/** name for shared framebuffer used for gpu marching cubes operations */
+const FramebufferName = 'marching-cubes'
+
+const ActiveVoxelsSchema = {
+    drawCount: ValueSpec('number'),
+    instanceCount: ValueSpec('number'),
+    aPosition: AttributeSpec('float32', 2, 0),
+
+    tTriCount: TextureSpec('image-uint8', 'alpha', 'ubyte', 'nearest'),
+    tVolumeData: TextureSpec('texture', 'rgba', 'ubyte', 'nearest'),
+    uIsoValue: UniformSpec('f'),
+
+    uGridDim: UniformSpec('v3'),
+    uGridTexDim: UniformSpec('v3'),
+}
+
+function getActiveVoxelsRenderable(ctx: WebGLContext, volumeData: Texture, gridDimensions: Vec3, isoValue: number) {
+    const values: Values<typeof ActiveVoxelsSchema> = {
+        drawCount: ValueCell.create(6),
+        instanceCount: ValueCell.create(1),
+        aPosition: ValueCell.create(QuadPositions),
+
+        tTriCount: ValueCell.create(getTriCount()),
+        tVolumeData: ValueCell.create(volumeData),
+        uIsoValue: ValueCell.create(isoValue),
+
+        uGridDim: ValueCell.create(gridDimensions),
+        uGridTexDim: ValueCell.create(Vec3.create(volumeData.width, volumeData.height, 0)),
+    }
+
+    const schema = { ...ActiveVoxelsSchema }
+    const shaderCode = ShaderCode(
+        require('mol-gl/shader/quad.vert').default,
+        require('mol-gl/shader/marching-cubes/active-voxels.frag').default,
+        { standardDerivatives: false, fragDepth: false }
+    )
+    const renderItem = createComputeRenderItem(ctx, 'triangles', shaderCode, schema, values)
+
+    return createComputeRenderable(renderItem, values);
+}
+
+function setRenderingDefaults(gl: GLRenderingContext) {
+    gl.disable(gl.CULL_FACE)
+    gl.disable(gl.BLEND)
+    gl.disable(gl.DEPTH_TEST)
+    gl.depthMask(false)
+}
+
+export function calcActiveVoxels(ctx: WebGLContext, cornerTex: Texture, gridDimensions: Vec3, isoValue: number) {
+    const { gl, framebufferCache } = ctx
+    const { width, height } = cornerTex
+
+    const framebuffer = framebufferCache.get(FramebufferName).value
+    framebuffer.bind()
+
+    const activeVoxelsTex = createTexture(ctx, 'image-float32', 'rgba', 'float', 'nearest')
+    activeVoxelsTex.define(width, height)
+
+    const renderable = getActiveVoxelsRenderable(ctx, cornerTex, gridDimensions, isoValue)
+    renderable.update()
+    renderable.use()
+
+    activeVoxelsTex.attachFramebuffer(framebuffer, 0)
+    setRenderingDefaults(gl)
+    gl.viewport(0, 0, width, height)
+    renderable.render()
+
+    // const at = readTexture(ctx, activeVoxelsTex)
+    // console.log('at', at)
+
+    return activeVoxelsTex
+}
\ No newline at end of file

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

+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { createComputeRenderable } from '../../renderable'
+import { WebGLContext } from '../../webgl/context';
+import { createComputeRenderItem } from '../../webgl/render-item';
+import { AttributeSpec, Values, TextureSpec, ValueSpec, UniformSpec } from '../../renderable/schema';
+import { Texture, createTexture } from 'mol-gl/webgl/texture';
+import { ShaderCode } from 'mol-gl/shader-code';
+import { ValueCell } from 'mol-util';
+import { GLRenderingContext } from 'mol-gl/webgl/compat';
+import { Vec3, Vec2, Mat4 } from 'mol-math/linear-algebra';
+import { QuadPositions } from '../util';
+import { HistogramPyramid } from '../histogram-pyramid/reduction';
+import { getTriIndices } from './tables';
+
+/** name for shared framebuffer used for gpu marching cubes operations */
+const FramebufferName = 'marching-cubes'
+
+const IsosurfaceSchema = {
+    drawCount: ValueSpec('number'),
+    instanceCount: ValueSpec('number'),
+    aPosition: AttributeSpec('float32', 2, 0),
+
+    tTriIndices: TextureSpec('image-uint8', 'alpha', 'ubyte', 'nearest'),
+    tActiveVoxelsPyramid: TextureSpec('texture', 'rgba', 'float', 'nearest'),
+    tActiveVoxelsBase: TextureSpec('texture', 'rgba', 'float', 'nearest'),
+    tVolumeData: TextureSpec('texture', 'rgba', 'ubyte', 'nearest'),
+    tActiveVoxelsTotal: TextureSpec('texture', 'rgba', 'float', 'nearest'),
+    uIsoValue: UniformSpec('f'),
+
+    uSize: UniformSpec('f'),
+    uLevels: UniformSpec('f'),
+
+    uGridDim: UniformSpec('v3'),
+    uGridTexDim: UniformSpec('v3'),
+    uGridTransform: UniformSpec('m4'),
+
+    uScale: UniformSpec('v2'),
+}
+
+function getIsosurfaceRenderable(ctx: WebGLContext, activeVoxelsPyramid: Texture, activeVoxelsBase: Texture, volumeData: Texture, activeVoxelsTotal: Texture, gridDimensions: Vec3, transform: Mat4, isoValue: number, levels: number, scale: Vec2) {
+    // console.log('uSize', Math.pow(2, levels))
+    const values: Values<typeof IsosurfaceSchema> = {
+        drawCount: ValueCell.create(6),
+        instanceCount: ValueCell.create(1),
+        aPosition: ValueCell.create(QuadPositions),
+
+        tTriIndices: ValueCell.create(getTriIndices()),
+        tActiveVoxelsPyramid: ValueCell.create(activeVoxelsPyramid),
+        tActiveVoxelsBase: ValueCell.create(activeVoxelsBase),
+        tVolumeData: ValueCell.create(volumeData),
+        tActiveVoxelsTotal: ValueCell.create(activeVoxelsTotal),
+        uIsoValue: ValueCell.create(isoValue),
+
+        uSize: ValueCell.create(Math.pow(2, levels)),
+        uLevels: ValueCell.create(levels),
+
+        uGridDim: ValueCell.create(gridDimensions),
+        uGridTexDim: ValueCell.create(Vec3.create(volumeData.width, volumeData.height, 0)),
+        uGridTransform: ValueCell.create(transform),
+
+        uScale: ValueCell.create(scale),
+    }
+
+    const schema = { ...IsosurfaceSchema }
+    const shaderCode = ShaderCode(
+        require('mol-gl/shader/quad.vert').default,
+        require('mol-gl/shader/marching-cubes/isosurface.frag').default,
+        { drawBuffers: true }
+    )
+    const renderItem = createComputeRenderItem(ctx, 'triangles', shaderCode, schema, values)
+
+    return createComputeRenderable(renderItem, values);
+}
+
+function setRenderingDefaults(gl: GLRenderingContext) {
+    gl.disable(gl.CULL_FACE)
+    gl.disable(gl.BLEND)
+    gl.disable(gl.DEPTH_TEST)
+    gl.depthMask(false)
+}
+
+export function createIsosurfaceBuffers(ctx: WebGLContext, activeVoxelsBase: Texture, volumeData: Texture, histogramPyramid: HistogramPyramid, gridDimensions: Vec3, transform: Mat4, isoValue: number) {
+    const { gl, framebufferCache } = ctx
+    const { pyramidTex, totalTex, height, levels, scale, count } = histogramPyramid
+
+    const framebuffer = framebufferCache.get(FramebufferName).value
+    framebuffer.bind()
+
+    const verticesTex = createTexture(ctx, 'image-float32', 'rgba', 'float', 'nearest')
+    verticesTex.define(pyramidTex.width, pyramidTex.height)
+
+    // const infoTex = createTexture(ctx, 'image-float32', 'rgba', 'float', 'nearest')
+    // infoTex.define(pyramidTex.width, pyramidTex.height)
+
+    // const pointTexA = createTexture(ctx, 'image-float32', 'rgba', 'float', 'nearest')
+    // pointTexA.define(pyramidTex.width, pyramidTex.height)
+
+    // const pointTexB = createTexture(ctx, 'image-float32', 'rgba', 'float', 'nearest')
+    // pointTexB.define(pyramidTex.width, pyramidTex.height)
+
+    // const coordTex = createTexture(ctx, 'image-float32', 'rgba', 'float', 'nearest')
+    // coordTex.define(pyramidTex.width, pyramidTex.height)
+
+    // const indexTex = createTexture(ctx, 'image-float32', 'rgba', 'float', 'nearest')
+    // indexTex.define(pyramidTex.width, pyramidTex.height)
+
+    const pr = getIsosurfaceRenderable(ctx, pyramidTex, activeVoxelsBase, volumeData, totalTex, gridDimensions, transform, isoValue, levels, scale)
+    pr.update()
+    pr.use()
+
+    verticesTex.attachFramebuffer(framebuffer, 0)
+    // infoTex.attachFramebuffer(framebuffer, 1)
+    // pointTexA.attachFramebuffer(framebuffer, 2)
+    // pointTexB.attachFramebuffer(framebuffer, 3)
+    // coordTex.attachFramebuffer(framebuffer, 4)
+    // indexTex.attachFramebuffer(framebuffer, 5)
+
+    // const { drawBuffers } = ctx.extensions
+    // if (!drawBuffers) throw new Error('need draw buffers')
+
+    // drawBuffers.drawBuffers([
+    //     drawBuffers.COLOR_ATTACHMENT0,
+    //     drawBuffers.COLOR_ATTACHMENT1,
+    //     drawBuffers.COLOR_ATTACHMENT2,
+    //     drawBuffers.COLOR_ATTACHMENT3,
+    //     drawBuffers.COLOR_ATTACHMENT4,
+    //     drawBuffers.COLOR_ATTACHMENT5
+    // ])
+
+    setRenderingDefaults(gl)
+    gl.viewport(0, 0, pyramidTex.width, pyramidTex.height)
+    gl.scissor(0, 0, pyramidTex.width, height)
+    pr.render()
+    gl.disable(gl.SCISSOR_TEST)
+
+    // const vt = readTexture(ctx, verticesTex, pyramidTex.width, height)
+    // console.log('vt', vt)
+    // const vertices = new Float32Array(3 * compacted.count)
+    // for (let i = 0; i < compacted.count; ++i) {
+    //     vertices[i * 3] = vt.array[i * 4]
+    //     vertices[i * 3 + 1] = vt.array[i * 4 + 1]
+    //     vertices[i * 3 + 2] = vt.array[i * 4 + 2]
+    // }
+    // console.log('vertices', vertices)
+
+    // const it = readTexture(ctx, infoTex, pyramidTex.width, height)
+    // console.log('info', it.array.subarray(0, 4 * compacted.count))
+
+    // const pat = readTexture(ctx, pointTexA, pyramidTex.width, height)
+    // console.log('point a', pat.array.subarray(0, 4 * compacted.count))
+
+    // const pbt = readTexture(ctx, pointTexB, pyramidTex.width, height)
+    // console.log('point b', pbt.array.subarray(0, 4 * compacted.count))
+
+    // const ct = readTexture(ctx, coordTex, pyramidTex.width, height)
+    // console.log('coord', ct.array.subarray(0, 4 * compacted.count))
+
+    // const idxt = readTexture(ctx, indexTex, pyramidTex.width, height)
+    // console.log('index', idxt.array.subarray(0, 4 * compacted.count))
+
+    // const { field, idField } = await fieldFromTexture2d(ctx, volumeData, gridDimensions)
+    // console.log({ field, idField })
+
+    // const valuesA = new Float32Array(compacted.count)
+    // const valuesB = new Float32Array(compacted.count)
+    // for (let i = 0; i < compacted.count; ++i) {
+    //     valuesA[i] = field.space.get(field.data, pat.array[i * 4], pat.array[i * 4 + 1], pat.array[i * 4 + 2])
+    //     valuesB[i] = field.space.get(field.data, pbt.array[i * 4], pbt.array[i * 4 + 1], pbt.array[i * 4 + 2])
+    // }
+    // console.log('valuesA', valuesA)
+    // console.log('valuesB', valuesB)
+
+    const vertexTexture = verticesTex
+    const normalBuffer = new Float32Array(count * 3)
+    const groupBuffer = new Float32Array(count)
+
+    return { vertexTexture, normalBuffer, groupBuffer, vertexCount: count }
+}
\ No newline at end of file

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

+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { TriTable,  } from 'mol-geo/util/marching-cubes/tables';
+import { TextureImage, createTextureImage } from 'mol-gl/renderable/util';
+
+let TriCount: TextureImage<Uint8Array> | undefined
+export function getTriCount(): TextureImage<Uint8Array> {
+    if (TriCount !== undefined) return TriCount
+    TriCount = createTextureImage(16 * 16, 1, Uint8Array)
+    const { array } = TriCount
+    for (let i = 0, il = TriTable.length; i < il; ++i) {
+        array[i] = TriTable[i].length / 3
+    }
+    return TriCount
+}
+
+let TriIndices: TextureImage<Uint8Array> | undefined
+export function getTriIndices(): TextureImage<Uint8Array> {
+    if (TriIndices !== undefined) return TriIndices
+    TriIndices = createTextureImage(64 * 64, 1, Uint8Array)
+    const { array } = TriIndices
+    for (let i = 0, il = TriTable.length; i < il; ++i) {
+        for (let j = 0; j < 16; ++j) {
+            if (j < TriTable[i].length) {
+                array[i * 16 + j] = TriTable[i][j]
+            } else {
+                array[i * 16 + j] = 255
+            }
+        }
+    }
+    return TriIndices
+}
\ No newline at end of file

src/mol-gl/compute/util.ts

+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { WebGLContext } from 'mol-gl/webgl/context';
+import { Texture } from 'mol-gl/webgl/texture';
+import { printTextureImage } from 'mol-gl/renderable/util';
+import { defaults } from 'mol-util';
+
+export const QuadPositions = new Float32Array([
+     1.0,  1.0,  -1.0,  1.0,  -1.0, -1.0, // First triangle:
+    -1.0, -1.0,   1.0, -1.0,   1.0,  1.0  // Second triangle:
+])
+
+export function readTexture(ctx: WebGLContext, texture: Texture, width?: number, height?: number) {
+    const { gl, framebufferCache } = ctx
+    width = defaults(width, texture.width)
+    height = defaults(height, texture.height)
+    const framebuffer = framebufferCache.get('read-texture').value
+    const array = texture.type === gl.UNSIGNED_BYTE ? new Uint8Array(width * height * 4) : new Float32Array(width * height * 4)
+    framebuffer.bind()
+    texture.attachFramebuffer(framebuffer, 0)
+    ctx.readPixels(0, 0, width, height, array)
+
+    return { array, width, height }
+}
+
+export function printTexture(ctx: WebGLContext, texture: Texture, scale: number) {
+    printTextureImage(readTexture(ctx, texture), scale)
+}
\ No newline at end of file

src/mol-gl/shader/histogram-pyramid/reduction.frag

+precision mediump float;
+precision mediump sampler2D;
+
+// input texture (previous level used to evaluate the new level)
+uniform sampler2D tPreviousLevel;
+
+// 1/size of the previous level texture.
+uniform float uSize;
+
+varying vec2 vCoordinate;
+
+void main(void) {
+    float k = 0.5 * uSize;
+    vec2 position = floor(vCoordinate / uSize) * uSize;
+    float a = texture2D(tPreviousLevel, position + vec2(0., 0.)).r;
+    float b = texture2D(tPreviousLevel, position + vec2(k, 0.)).r;
+    float c = texture2D(tPreviousLevel, position + vec2(0., k)).r;
+    float d = texture2D(tPreviousLevel, position + vec2(k, k)).r;
+    gl_FragColor.a = a;
+    gl_FragColor.b = a + b;
+    gl_FragColor.g = gl_FragColor.b + c;
+    gl_FragColor.r = gl_FragColor.g + d;
+}
\ No newline at end of file

src/mol-gl/shader/histogram-pyramid/sum.frag

+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+precision mediump float;
+precision mediump sampler2D;
+
+uniform sampler2D tTexture;
+
+// uv coordinate of each fragment obtained from vertex shader.
+varying vec2 vCoordinate;
+
+#pragma glslify: encodeFloatRGB = require(../utils/encode-float-rgb.glsl)
+
+void main(void) {
+    gl_FragColor = vec4(encodeFloatRGB(texture2D(tTexture, vCoordinate).r), 1.0);
+}
\ No newline at end of file

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

+precision mediump float;
+precision mediump sampler2D;
+
+uniform sampler2D tTriCount;
+uniform sampler2D tVolumeData;
+
+uniform float uIsoValue;
+uniform vec3 uGridDim;
+uniform vec3 uGridTexDim;
+
+varying vec2 vCoordinate;
+
+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) / uGridDim;
+    return posXYZ;
+}
+
+float intDiv(float a, float b) { return float(int(a) / int(b)); }
+float intMod(float a, float b) { return a - b * float(int(a) / int(b)); }
+
+vec4 texture3dFrom2dNearest(sampler2D tex, vec3 pos, vec3 gridDim, vec2 texDim) {
+    float zSlice = floor(pos.z * gridDim.z + 0.5); // round to nearest z-slice
+    float column = intMod(zSlice * gridDim.x, texDim.x) / gridDim.x;
+    float row = floor(intDiv(zSlice * gridDim.x, texDim.x));
+    vec2 coord = (vec2(column * gridDim.x, row * gridDim.y) + (pos.xy * gridDim.xy)) / texDim;
+    return texture2D(tex, coord);
+}
+
+vec4 voxel(vec3 pos) {
+    return texture3dFrom2dNearest(tVolumeData, pos, uGridDim, uGridTexDim.xy);
+}
+
+void main(void) {
+    vec3 posXYZ = index3dFrom2d(vCoordinate);
+
+    // 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 + vec3(1., 0., 0.) / uGridDim).a, uIsoValue)
+        + 4. * step(voxel(posXYZ + vec3(1., 1., 0.) / uGridDim).a, uIsoValue)
+        + 8. * step(voxel(posXYZ + vec3(0., 1., 0.) / uGridDim).a, uIsoValue)
+        + 16. * step(voxel(posXYZ + vec3(0., 0., 1.) / uGridDim).a, uIsoValue)
+        + 32. * step(voxel(posXYZ + vec3(1., 0., 1.) / uGridDim).a, uIsoValue)
+        + 64. * step(voxel(posXYZ + vec3(1., 1., 1.) / uGridDim).a, uIsoValue)
+        + 128. * step(voxel(posXYZ + vec3(0., 1., 1.) / uGridDim).a, uIsoValue);
+    c *= step(c, 254.);
+
+    // 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 * 255.0 * 3.0), c);
+}
\ No newline at end of file

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

+precision highp float;
+precision highp sampler2D;
+
+uniform sampler2D tActiveVoxelsPyramid;
+uniform sampler2D tActiveVoxelsBase;
+uniform sampler2D tActiveVoxelsTotal;
+uniform sampler2D tVolumeData;
+uniform sampler2D tTriIndices;
+
+uniform float uIsoValue;
+uniform float uLevels;
+uniform float uSize;
+
+uniform vec3 uGridDim;
+uniform vec3 uGridTexDim;
+uniform mat4 uGridTransform;
+
+// scale to volume data coord
+uniform vec2 uScale;
+
+varying vec2 vCoordinate;
+
+// 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.);
+// const vec3 c4  = vec3(0., 0., 1.);
+// const vec3 c5  = vec3(1., 0., 1.);
+// const vec3 c6  = vec3(1., 1., 1.);
+// const vec3 c7  = vec3(0., 1., 1.);
+
+const vec3 p0 = vec3(1., 0., 0.);
+const vec3 p1 = vec3(1., 1., 0.);
+const vec3 p2 = vec3(0., 1., 0.);
+const vec3 p3 = vec3(0., 0., 1.);
+const vec3 p4 = vec3(1., 0., 1.);
+const vec3 p5 = vec3(1., 1., 1.);
+const vec3 p6 = vec3(0., 1., 1.);
+
+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;
+}
+
+float intDiv(float a, float b) { return float(int(a) / int(b)); }
+float intMod(float a, float b) { return a - b * float(int(a) / int(b)); }
+
+vec4 texture3dFrom2dNearest(sampler2D tex, vec3 pos, vec3 gridDim, vec2 texDim) {
+    float zSlice = floor(pos.z * gridDim.z + 0.5); // round to nearest z-slice
+    float column = intMod(zSlice * gridDim.x, texDim.x) / gridDim.x;
+    float row = floor(intDiv(zSlice * gridDim.x, texDim.x));
+    vec2 coord = (vec2(column * gridDim.x, row * gridDim.y) + (pos.xy * gridDim.xy)) / texDim;
+    return texture2D(tex, coord + 0.5 / texDim);
+}
+
+vec4 voxel(vec3 pos) {
+    return texture3dFrom2dNearest(tVolumeData, pos, uGridDim, uGridTexDim.xy);
+}
+
+void main(void) {
+    // get 1D index
+    float vI = dot(floor(uSize * vCoordinate), vec2(1.0, uSize));
+
+    // ignore 1D indices outside of the grid
+    if(vI >= texture2D(tActiveVoxelsTotal, vec2(0.5)).r) discard;
+
+    float offset = uSize - 2.;
+    float k = 1. / uSize;
+
+    vec2 relativePosition = k * vec2(offset, 0.);
+    vec4 partialSums = texture2D(tActiveVoxelsPyramid, relativePosition);
+    float start = 0.;
+    vec4 starts = vec4(0.);
+    vec4 ends = vec4(0.);
+    float diff = 2.;
+    vec4 m = vec4(0.);
+    vec2 position = vec2(0.);
+    vec4 vI4 = vec4(vI);
+
+    // traverse the different levels of the pyramid
+    for(int i = 1; i < 12; i++) {
+        if(float(i) >= uLevels) break;
+
+        offset -= diff;
+        diff *= 2.;
+        relativePosition = position + k * vec2(offset, 0.);
+
+        ends = partialSums.wzyx + vec4(start);
+        starts = vec4(start, ends.xyz);
+        m = vec4(greaterThanEqual(vI4, starts)) * vec4(lessThan(vI4, ends));
+        relativePosition += m.y * vec2(k, 0.) + m.z * vec2(0., k) + m.w * vec2(k, k);
+
+        start = dot(m, starts);
+        position = 2. * (relativePosition - k * vec2(offset, 0.));
+        partialSums = texture2D(tActiveVoxelsPyramid, relativePosition);
+    }
+
+    ends = partialSums.wzyx + vec4(start);
+    starts = vec4(start, ends.xyz);
+    m = vec4(greaterThanEqual(vI4, starts)) * vec4(lessThan(vI4, ends));
+    position += m.y * vec2(k, 0.) + m.z * vec2(0., k) + m.w * vec2(k, k);
+
+    vec2 coord2d = position * uScale;
+    vec3 coord3d = floor(index3dFrom2d(coord2d) + 0.5);
+
+    float edgeIndex = texture2D(tActiveVoxelsBase, position * uScale).a;
+
+    // current vertex for the up to 15 MC cases
+    float currentVertex = vI - dot(m, starts);
+
+    // get index into triIndices table
+    float mcIndex = 16. * edgeIndex + currentVertex;
+    vec4 mcData = texture2D(tTriIndices, vec2(intMod(mcIndex, 64.), floor(mcIndex / 64.)) / 64.);
+    mcIndex = floor(mcData.a * 255.0 + 0.5);
+
+    // bit mask to avoid conditionals (see comment below) for getting MC case corner
+    vec4 m0 = vec4(mcIndex);
+
+    // get edge value masks
+    vec4 m1 = vec4(equal(m0, vec4(0., 1., 2., 3.)));
+    vec4 m2 = vec4(equal(m0, vec4(4., 5., 6., 7.)));
+    vec4 m3 = vec4(equal(m0, vec4(8., 9., 10., 11.)));
+
+    // apply bit masks
+    vec3 b0 = coord3d +
+                m1.y * p0 +
+                m1.z * p1 +
+                m1.w * p2 +
+                m2.x * p3 +
+                m2.y * p4 +
+                m2.z * p5 +
+                m2.w * p6 +
+                m3.y * p0 +
+                m3.z * p1 +
+                m3.w * p2;
+    vec3 b1 = coord3d +
+                m1.x * p0 +
+                m1.y * p1 +
+                m1.z * p2 +
+                m2.x * p4 +
+                m2.y * p5 +
+                m2.z * p6 +
+                m2.w * p3 +
+                m3.x * p3 +
+                m3.y * p4 +
+                m3.z * p5 +
+                m3.w * p6;
+
+    // vec3 b0 = coord3d;
+    // vec3 b1 = coord3d;
+    // if (mcIndex == 0.0) {
+    //     b0 += c0; b1 += c1;
+    // } else if (mcIndex == 1.0) {
+    //     b0 += c1; b1 += c2;
+    // } else if (mcIndex == 2.0) {
+    //     b0 += c2; b1 += c3;
+    // } else if (mcIndex == 3.0) {
+    //     b0 += c3; b1 += c0;
+    // } else if (mcIndex == 4.0) {
+    //     b0 += c4; b1 += c5;
+    // } else if (mcIndex == 5.0) {
+    //     b0 += c5; b1 += c6;
+    // } else if (mcIndex == 6.0) {
+    //     b0 += c6; b1 += c7;
+    // } else if (mcIndex == 7.0) {
+    //     b0 += c7; b1 += c4;
+    // } else if (mcIndex == 8.0) {
+    //     b0 += c0; b1 += c4;
+    // } else if (mcIndex == 9.0) {
+    //     b0 += c1; b1 += c5;
+    // } else if (mcIndex == 10.0) {
+    //     b0 += c2; b1 += c6;
+    // } else if (mcIndex == 11.0) {
+    //     b0 += c3; b1 += c7;
+    // }
+    // b0 = floor(b0 + 0.5);
+    // b1 = floor(b1 + 0.5);
+
+    float v0 = voxel(b0 / uGridDim).a;
+    float v1 = voxel(b1 / uGridDim).a;
+
+    float t = (uIsoValue - v0) / (v0 - v1);
+    gl_FragColor.xyz = b0 + t * (b0 - b1);
+}
\ No newline at end of file

src/mol-gl/shader/quad.vert

+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+precision mediump float;
+
+attribute vec2 aPosition;
+
+// the output UV coordinate for each fragment
+varying vec2 vCoordinate;
+
+const vec2 scale = vec2(0.5, 0.5);
+
+void main(void) {
+    vCoordinate  = aPosition * scale + scale; // scale vertex attribute to [0,1] range
+    gl_Position = vec4(aPosition, 0.0, 1.0);
+}
\ No newline at end of file