gaussian density gpu calculation fixes

src/mol-geo/representation/structure/visual/gaussian-density-volume.ts

@@ -46,7 +46,10 @@ export function GaussianDensityVolumeVisual(): UnitsVisual<GaussianDensityVolume
         setUpdateState: (state: VisualUpdateState, newProps: GaussianDensityVolumeProps, currentProps: GaussianDensityVolumeProps) => {
             if (newProps.resolution !== currentProps.resolution) state.createGeometry = true
             if (newProps.radiusOffset !== currentProps.radiusOffset) state.createGeometry = true
-            if (newProps.smoothness !== currentProps.smoothness) state.createGeometry = true
+            if (newProps.smoothness !== currentProps.smoothness) {
+                state.createGeometry = true
+                newProps.isoValueAbsolute = Math.exp(-newProps.smoothness)
+            }
             if (newProps.useGpu !== currentProps.useGpu) state.createGeometry = true
             if (newProps.ignoreCache !== currentProps.ignoreCache) state.createGeometry = true
         }

src/mol-gl/renderable/direct-volume.ts

@@ -34,15 +34,15 @@ export const DirectVolumeBaseSchema = {
 export type DirectVolumeBaseSchema = typeof DirectVolumeBaseSchema
 export type DirectVolumeBaseValues = Values<DirectVolumeBaseSchema>
 
-function getInternalValues(ctx: Context, id: number, version: '100es' | '300es'): InternalValues {
+function getInternalValues(ctx: Context, id: number): InternalValues {
     return {
         uObjectId: ValueCell.create(id)
     }
 }
 
-function DirectVolumeRenderable<T extends DirectVolumeBaseValues, S extends DirectVolumeBaseSchema>(ctx: Context, id: number, values: T, state: RenderableState, schema: S, version: '100es' | '300es'): Renderable<T> {
+function DirectVolumeRenderable<T extends DirectVolumeBaseValues, S extends DirectVolumeBaseSchema>(ctx: Context, id: number, values: T, state: RenderableState, schema: S): Renderable<T> {
     const fullSchema = Object.assign({}, GlobalUniformSchema, InternalSchema, schema)
-    const internalValues = getInternalValues(ctx, id, version)
+    const internalValues = getInternalValues(ctx, id)
     const fullValues = Object.assign({}, values, internalValues)
     const shaderCode = DirectVolumeShaderCode
     const renderItem = createRenderItem(ctx, 'triangles', shaderCode, fullSchema, fullValues)
@@ -65,7 +65,7 @@ export type DirectVolume2dSchema = typeof DirectVolume2dSchema
 export type DirectVolume2dValues = Values<DirectVolume2dSchema>
 
 export function DirectVolume2dRenderable(ctx: Context, id: number, values: DirectVolume2dValues, state: RenderableState): Renderable<DirectVolume2dValues> {
-    return DirectVolumeRenderable(ctx, id, values, state, DirectVolume2dSchema, '100es')
+    return DirectVolumeRenderable(ctx, id, values, state, DirectVolume2dSchema)
 }
 
 // via 3d texture
@@ -79,5 +79,5 @@ export type DirectVolume3dSchema = typeof DirectVolume3dSchema
 export type DirectVolume3dValues = Values<DirectVolume3dSchema>
 
 export function DirectVolume3dRenderable(ctx: Context, id: number, values: DirectVolume3dValues, state: RenderableState): Renderable<DirectVolume3dValues> {
-    return DirectVolumeRenderable(ctx, id, values, state, DirectVolume3dSchema, '300es')
+    return DirectVolumeRenderable(ctx, id, values, state, DirectVolume3dSchema)
 }
\ No newline at end of file

src/mol-gl/shader/direct-volume.frag

@@ -69,7 +69,7 @@ const vec3 color = vec3(0.45, 0.55, 0.8);
 vec4 raymarch(vec3 startLoc, vec3 step, vec3 viewDir) {
     vec3 scaleVol = vec3(1.0) / uGridDim;
     vec3 pos = startLoc + scaleVol * 0.5;
-    float prevValue = -127.0;
+    float prevValue = -1.0;
     float value = 0.0;
     vec4 src = vec4(0.0);
     vec4 dst = vec4(0.0);

src/mol-gl/shader/gaussian-density.frag

@@ -5,7 +5,7 @@
  * @author Michael Krone <michael.krone@uni-tuebingen.de>
  */
 
-precision mediump float;
+precision highp float;
 
 varying vec3 position;
 varying float radius;
@@ -31,27 +31,28 @@ uniform float uAlpha;
     layout(location = 7) out vec4 out7;
 #endif
 
-float calcDensity(float x, float y, float z) {
+float calcDensity(float x, float y, float z, float radiusSq) {
     vec3 fragPos = vec3(x, y, z) / uGridDim;
-    float dist = length(fragPos * uBboxSize - position * uBboxSize);
-    float density = exp(-uAlpha * ((dist * dist) / (radius * radius)));
+    float dist = distance(fragPos * uBboxSize, position * uBboxSize);
+    float density = exp(-uAlpha * ((dist * dist) / radiusSq));
     return density;
 }
 
 const vec3 color = vec3(1.0, 1.0, 1.0);
 
 void main() {
+    float radiusSq = radius * radius;
     vec2 v = gl_FragCoord.xy - vec2(uCurrentX, uCurrentY) - 0.5;
-    gl_FragColor = vec4(color, calcDensity(v.x, v.y, uCurrentSlice));
+    out_FragColor = vec4(color, calcDensity(v.x, v.y, uCurrentSlice, radiusSq));
     #if dDrawBuffers >= 4
-        out1 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 1.0));
-        out2 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 2.0));
-        out3 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 3.0));
+        out1 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 1.0, radiusSq));
+        out2 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 2.0, radiusSq));
+        out3 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 3.0, radiusSq));
     #endif
     #if dDrawBuffers >= 8
-        out4 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 4.0));
-        out5 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 5.0));
-        out6 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 6.0));
-        out7 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 7.0));
+        out4 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 4.0, radiusSq));
+        out5 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 5.0, radiusSq));
+        out6 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 6.0, radiusSq));
+        out7 = vec4(color, calcDensity(v.x, v.y, uCurrentSlice + 7.0, radiusSq));
     #endif
 }
\ No newline at end of file

src/mol-gl/shader/gaussian-density.vert

@@ -5,7 +5,7 @@
  * @author Michael Krone <michael.krone@uni-tuebingen.de>
  */
 
-precision mediump float;
+precision highp float;
 
 attribute vec3 aPosition;
 attribute float aRadius;

src/mol-gl/webgl/texture.ts

@@ -198,6 +198,7 @@ export function createTexture(ctx: Context, kind: TextureKind, _format: TextureF
             // unpack alignment of 1 since we use textures only for data
             gl.pixelStorei(gl.UNPACK_ALIGNMENT, 1);
             gl.pixelStorei(gl.UNPACK_COLORSPACE_CONVERSION_WEBGL, gl.NONE);
+            gl.pixelStorei(gl.UNPACK_PREMULTIPLY_ALPHA_WEBGL, 0);
             if (target === gl.TEXTURE_2D) {
                 const { array, width: _width, height: _height } = data as TextureImage<any>
                 width = _width, height = _height;

src/mol-math/geometry/gaussian-density/gpu.ts

@@ -25,11 +25,9 @@ export async function GaussianDensityGPU(ctx: RuntimeContext, position: Position
 
     const { transform, texture, gridDimension } = await GaussianDensityTexture(ctx, webgl, position, box, radius, props)
 
-    console.time('gpu gaussian density read')
     const field = webgl.maxDrawBuffers > 0 ?
         fieldFromTexture3d(webgl, texture, gridDimension) :
         fieldFromTexture2d(webgl, texture, gridDimension)
-    console.timeEnd('gpu gaussian density read')
 
     const idData = field.space.create()
     const idField = Tensor.create(field.space, idData)
@@ -38,9 +36,11 @@ export async function GaussianDensityGPU(ctx: RuntimeContext, position: Position
 }
 
 export async function GaussianDensityTexture(ctx: RuntimeContext, webgl: Context, position: PositionData, box: Box3D, radius: (index: number) => number, props: GaussianDensityProps, oldTexture?: Texture): Promise<DensityTextureData> {
+    console.time(`GaussianDensityTexture, ${webgl.maxDrawBuffers > 0 ? 'multi' : 'single'}`)
     const { texture, scale, bbox, dim } = webgl.maxDrawBuffers > 0 ?
         await GaussianDensityMultiDrawBuffer(ctx, webgl, position, box, radius, props, oldTexture) :
         await GaussianDensitySingleDrawBuffer(ctx, webgl, position, box, radius, props, oldTexture)
+    console.timeEnd(`GaussianDensityTexture, ${webgl.maxDrawBuffers > 0 ? 'multi' : 'single'}`)
 
     const transform = Mat4.identity()
     Mat4.fromScaling(transform, scale)
@@ -283,12 +283,13 @@ function setRenderingDefaults(gl: GLRenderingContext) {
     gl.frontFace(gl.CCW)
     gl.cullFace(gl.BACK)
 
-    gl.blendFunc(gl.SRC_ALPHA, gl.ONE_MINUS_SRC_ALPHA)
+    gl.blendFunc(gl.ONE, gl.ONE)
     gl.blendEquation(gl.FUNC_ADD)
     gl.enable(gl.BLEND)
 }
 
 function fieldFromTexture2d(ctx: Context, texture: Texture, dim: Vec3) {
+    console.time('fieldFromTexture2d')
     const { gl } = ctx
     const [ dx, dy, dz ] = dim
     const { width, height } = texture
@@ -316,7 +317,7 @@ function fieldFromTexture2d(ctx: Context, texture: Texture, dim: Vec3) {
         }
         for (let iy = 0; iy < dy; ++iy) {
             for (let ix = 0; ix < dx; ++ix) {
-                data[idx] = image[4 * (tmpCol * dx + (iy + tmpRow) * width + ix)] / 255
+                data[idx] = image[4 * (tmpCol * dx + (iy + tmpRow) * width + ix) + 3] / 255
                 idx++
             }
         }
@@ -324,11 +325,13 @@ function fieldFromTexture2d(ctx: Context, texture: Texture, dim: Vec3) {
     }
 
     framebuffer.destroy()
+    console.timeEnd('fieldFromTexture2d')
 
     return field
 }
 
 function fieldFromTexture3d(ctx: Context, texture: Texture, dim: Vec3) {
+    console.time('fieldFromTexture3d')
     const { gl } = ctx
     const { width, height, depth } = texture
 
@@ -347,13 +350,14 @@ function fieldFromTexture3d(ctx: Context, texture: Texture, dim: Vec3) {
         gl.readPixels(0, 0, width, height, gl.RGBA, gl.UNSIGNED_BYTE, slice)
         for (let iy = 0; iy < height; ++iy) {
             for (let ix = 0; ix < width; ++ix) {
-                data[j] = slice[4 * (iy * width + ix)] / 255
+                data[j] = slice[4 * (iy * width + ix) + 3] / 255
                 ++j
             }
         }
     }
 
     framebuffer.destroy()
+    console.timeEnd('fieldFromTexture3d')
 
     return field
 }
\ No newline at end of file