diff --git a/src/mol-gl/compute/histogram-pyramid/reduction.ts b/src/mol-gl/compute/histogram-pyramid/reduction.ts
new file mode 100644
index 0000000000000000000000000000000000000000..c8cb905906c1fe089254abb5305a46e72a330cd5
--- /dev/null
+++ b/src/mol-gl/compute/histogram-pyramid/reduction.ts
@@ -0,0 +1,165 @@
+/**
+ * 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
diff --git a/src/mol-gl/compute/histogram-pyramid/sum.ts b/src/mol-gl/compute/histogram-pyramid/sum.ts
new file mode 100644
index 0000000000000000000000000000000000000000..64d377238054b5a0f49c10584c25408b74e0057c
--- /dev/null
+++ b/src/mol-gl/compute/histogram-pyramid/sum.ts
@@ -0,0 +1,66 @@
+/**
+ * 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
diff --git a/src/mol-gl/compute/marching-cubes/active-voxels.ts b/src/mol-gl/compute/marching-cubes/active-voxels.ts
new file mode 100644
index 0000000000000000000000000000000000000000..3c84c4d7080976e2e24495b23aa7cdfd5468b51b
--- /dev/null
+++ b/src/mol-gl/compute/marching-cubes/active-voxels.ts
@@ -0,0 +1,90 @@
+/**
+ * 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
diff --git a/src/mol-gl/compute/marching-cubes/isosurface.ts b/src/mol-gl/compute/marching-cubes/isosurface.ts
new file mode 100644
index 0000000000000000000000000000000000000000..03f6ec0e7c8cf864afcbe6dcf536afa65f17df00
--- /dev/null
+++ b/src/mol-gl/compute/marching-cubes/isosurface.ts
@@ -0,0 +1,183 @@
+/**
+ * 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
diff --git a/src/mol-gl/compute/marching-cubes/tables.ts b/src/mol-gl/compute/marching-cubes/tables.ts
new file mode 100644
index 0000000000000000000000000000000000000000..3fafe4de31e5b1d5c7d6120b468b19ba483210e6
--- /dev/null
+++ b/src/mol-gl/compute/marching-cubes/tables.ts
@@ -0,0 +1,36 @@
+/**
+ * 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
diff --git a/src/mol-gl/compute/util.ts b/src/mol-gl/compute/util.ts
new file mode 100644
index 0000000000000000000000000000000000000000..2c3da2710b00a08295b0bad717d3aeef3d444d7e
--- /dev/null
+++ b/src/mol-gl/compute/util.ts
@@ -0,0 +1,32 @@
+/**
+ * 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
diff --git a/src/mol-gl/shader/histogram-pyramid/reduction.frag b/src/mol-gl/shader/histogram-pyramid/reduction.frag
new file mode 100644
index 0000000000000000000000000000000000000000..1307e20e3a6d189af5f3968024af96acd78a5565
--- /dev/null
+++ b/src/mol-gl/shader/histogram-pyramid/reduction.frag
@@ -0,0 +1,23 @@
+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
diff --git a/src/mol-gl/shader/histogram-pyramid/sum.frag b/src/mol-gl/shader/histogram-pyramid/sum.frag
new file mode 100644
index 0000000000000000000000000000000000000000..9e53bc2b5c4994c251bcc1adf5baa17b6a65c5d8
--- /dev/null
+++ b/src/mol-gl/shader/histogram-pyramid/sum.frag
@@ -0,0 +1,19 @@
+/**
+ * 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
diff --git a/src/mol-gl/shader/marching-cubes/active-voxels.frag b/src/mol-gl/shader/marching-cubes/active-voxels.frag
new file mode 100644
index 0000000000000000000000000000000000000000..5c5cffdc8296bea39398f56af318b9f3af0dab4f
--- /dev/null
+++ b/src/mol-gl/shader/marching-cubes/active-voxels.frag
@@ -0,0 +1,54 @@
+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
diff --git a/src/mol-gl/shader/marching-cubes/isosurface.frag b/src/mol-gl/shader/marching-cubes/isosurface.frag
new file mode 100644
index 0000000000000000000000000000000000000000..65a167edf4513369dec51c3d561b5c646d6aeb02
--- /dev/null
+++ b/src/mol-gl/shader/marching-cubes/isosurface.frag
@@ -0,0 +1,188 @@
+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
diff --git a/src/mol-gl/shader/quad.vert b/src/mol-gl/shader/quad.vert
new file mode 100644
index 0000000000000000000000000000000000000000..d8392e124de9b066efd63b30be5896818a050ccc
--- /dev/null
+++ b/src/mol-gl/shader/quad.vert
@@ -0,0 +1,19 @@
+/**
+ * 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