volume related fixes

src/apps/canvas/index.ts

@@ -26,8 +26,11 @@ if (pdbId) app.loadPdbIdOrMmcifUrl(pdbId, { assemblyId })
 // app.loadPdbIdOrMmcifUrl('3pqr')
 // app.loadCcp4Url('http://localhost:8091/ngl/data/3pqr-mode0.ccp4')
 
-app.loadPdbIdOrMmcifUrl('6DRV')
-app.loadCcp4Url('http://localhost:8091/ngl/data/betaGal.mrc')
+// app.loadPdbIdOrMmcifUrl('1lee')
+// app.loadCcp4Url('http://localhost:8091/ngl/data/1lee.ccp4')
+
+// app.loadPdbIdOrMmcifUrl('6DRV')
+// app.loadCcp4Url('http://localhost:8091/ngl/data/betaGal.mrc')
 
 // app.loadPdbIdOrMmcifUrl('3pqr')
 // app.loadVolCifUrl('https://webchem.ncbr.muni.cz/DensityServer/x-ray/3pqr/cell?space=fractional', true)
@@ -41,52 +44,4 @@ app.loadCcp4Url('http://localhost:8091/ngl/data/betaGal.mrc')
 // app.loadPdbIdOrMmcifUrl('http://localhost:8091/test/pdb-dev/carb/1B5F-carb.cif')
 // app.loadPdbIdOrMmcifUrl('http://localhost:8091/test/pdb-dev/carb/2HYV-carb.cif')
 // app.loadPdbIdOrMmcifUrl('http://localhost:8091/test/pdb-dev/carb/2WMG-carb.cif')
-// app.loadPdbIdOrMmcifUrl('http://localhost:8091/test/pdb-dev/carb/5KDS-carb.cif')
-
-const basisX = [
-    h.xlen,
-    0,
-    0
-  ]
-
-  const basisY = [
-    h.ylen * Math.cos(Math.PI / 180.0 * h.gamma),
-    h.ylen * Math.sin(Math.PI / 180.0 * h.gamma),
-    0
-  ]
-
-  const basisZ = [
-    h.zlen * Math.cos(Math.PI / 180.0 * h.beta),
-    h.zlen * (
-      Math.cos(Math.PI / 180.0 * h.alpha) -
-      Math.cos(Math.PI / 180.0 * h.gamma) *
-      Math.cos(Math.PI / 180.0 * h.beta)
-    ) / Math.sin(Math.PI / 180.0 * h.gamma),
-    0
-  ]
-  basisZ[ 2 ] = Math.sqrt(
-    h.zlen * h.zlen * Math.sin(Math.PI / 180.0 * h.beta) *
-    Math.sin(Math.PI / 180.0 * h.beta) - basisZ[ 1 ] * basisZ[ 1 ]
-  )
-
-  const basis = [ 0, basisX, basisY, basisZ ]
-  const nxyz = [ 0, h.MX, h.MY, h.MZ ]
-  const mapcrs = [ 0, h.MAPC, h.MAPR, h.MAPS ]
-
-  const matrix = new Matrix4()
-
-  matrix.set(
-    basis[ mapcrs[1] ][0] / nxyz[ mapcrs[1] ],
-    basis[ mapcrs[2] ][0] / nxyz[ mapcrs[2] ],
-    basis[ mapcrs[3] ][0] / nxyz[ mapcrs[3] ],
-    0,
-    basis[ mapcrs[1] ][1] / nxyz[ mapcrs[1] ],
-    basis[ mapcrs[2] ][1] / nxyz[ mapcrs[2] ],
-    basis[ mapcrs[3] ][1] / nxyz[ mapcrs[3] ],
-    0,
-    basis[ mapcrs[1] ][2] / nxyz[ mapcrs[1] ],
-    basis[ mapcrs[2] ][2] / nxyz[ mapcrs[2] ],
-    basis[ mapcrs[3] ][2] / nxyz[ mapcrs[3] ],
-    0,
-    0, 0, 0, 1
-  )
+// app.loadPdbIdOrMmcifUrl('http://localhost:8091/test/pdb-dev/carb/5KDS-carb.cif')

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

@@ -114,41 +114,44 @@ function createVolumeTexture3d(volume: VolumeData) {
 
     console.log('stats', stats)
 
+    let i = 0
+    for (let z = 0; z < depth; ++z) {
+        for (let y = 0; y < height; ++y) {
+            for (let x = 0; x < width; ++x) {
+                if (i < 100) {
+                    console.log(get(data, x, y, z), ((get(data, x, y, z) - stats.min) / (stats.max - stats.min)) * 255)
+                }
+                array[i + 3] = ((get(data, x, y, z) - stats.min) / (stats.max - stats.min)) * 255
+                // array[i + 3] = ((get(data, x, z, y) - stats.min) / (stats.max - stats.min)) * 255
+                // array[i + 3] = ((get(data, y, x, z) - stats.min) / (stats.max - stats.min)) * 255
+                // array[i + 3] = ((get(data, y, z, x) - stats.min) / (stats.max - stats.min)) * 255
+                // array[i + 3] = ((get(data, z, y, x) - stats.min) / (stats.max - stats.min)) * 255
+                // array[i + 3] = ((get(data, z, x, y) - stats.min) / (stats.max - stats.min)) * 255
+                i += 4
+            }
+        }
+    }
+
     // let i = 0
     // for (let z = 0; z < depth; ++z) {
-    //     for (let y = 0; y < height; ++y) {
-    //         for (let x = 0; x < width; ++x) {
+    //     for (let x = 0; x < width; ++x) {
+    //         for (let y = 0; y < height; ++y) {
     //             array[i + 3] = ((get(data, x, y, z) - stats.min) / (stats.max - stats.min)) * 255
-    //             // array[i + 3] = ((get(data, x, z, y) - stats.min) / (stats.max - stats.min)) * 255
-    //             // array[i + 3] = ((get(data, y, x, z) - stats.min) / (stats.max - stats.min)) * 255
-    //             // array[i + 3] = ((get(data, y, z, x) - stats.min) / (stats.max - stats.min)) * 255
-    //             // array[i + 3] = ((get(data, z, y, x) - stats.min) / (stats.max - stats.min)) * 255
-    //             // array[i + 3] = ((get(data, z, x, y) - stats.min) / (stats.max - stats.min)) * 255
     //             i += 4
     //         }
     //     }
     // }
 
     // let i = 0
-    // for (let z = 0; z < depth; ++z) {
-    //     for (let x = 0; x < width; ++x) {
-    //         for (let y = 0; y < height; ++y) {
+    // for (let x = 0; x < width; ++x) {
+    //     for (let y = 0; y < height; ++y) {
+    //         for (let z = 0; z < depth; ++z) {
     //             array[i + 3] = ((get(data, x, y, z) - stats.min) / (stats.max - stats.min)) * 255
     //             i += 4
     //         }
     //     }
     // }
 
-    let i = 0
-    for (let x = 0; x < width; ++x) {
-        for (let y = 0; y < height; ++y) {
-            for (let z = 0; z < depth; ++z) {
-                array[i + 3] = ((get(data, x, y, z) - stats.min) / (stats.max - stats.min)) * 255
-                i += 4
-            }
-        }
-    }
-
     return textureVolume
 }
 

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

@@ -164,9 +164,9 @@ vec4 raymarch(vec3 startLoc, vec3 step, vec3 viewDir) {
                     src.rgb = finalColor;
                     src.a = uAlpha;
 
-                    float marker = readFromTexture(tMarker, instance * float(uGroupCount) + group, uMarkerTexDim).a * 256.0;
+                    float marker = readFromTexture(tMarker, instance * float(uGroupCount) + group, uMarkerTexDim).a * 255.0;
                     if (marker > 0.1) {
-                        if (mod(marker, 2.0) < 0.1) {
+                        if (mod(marker, 2.0) > 0.1) {
                             src.rgb = mix(uHighlightColor, src.rgb, 0.3);
                         } else {
                             src.rgb = mix(uSelectColor, src.rgb, 0.3);

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

@@ -26,8 +26,7 @@ uniform mat4 uProjection;
 
 void main() {
     unitCoord = aPosition + vec3(0.5);
-    vec4 mvPosition = uModelView * uTransform * vec4(unitCoord * uGridDim + uBboxMin, 1.0);
-    // vec4 mvPosition = uModelView * uTransform * vec4(unitCoord, 1.0);
+    vec4 mvPosition = uModelView * uTransform * vec4(unitCoord * uGridDim, 1.0);
     origPos = unitCoord * uBboxSize + uBboxMin;
     instance = aInstance;
     gl_Position = uProjection * mvPosition;

src/mol-io/reader/ccp4/parser.ts

@@ -104,6 +104,18 @@ async function parseInternal(file: FileHandle, ctx: RuntimeContext): Promise<Res
         return Result.error(`ccp4 mode '${header.MODE}' unsupported`);
     }
 
+    // if the file was converted by mapmode2to0 - scale the data
+    // based on uglymol (https://github.com/uglymol/uglymol) by Marcin Wojdyr (wojdyr)
+    if (intView[39] === -128 && intView[40] === 127) {
+        values = new Float32Array(values)
+        // scaling f(x)=b1*x+b0 such that f(-128)=min and f(127)=max
+        const b1 = (header.AMAX - header.AMIN) / 255.0
+        const b0 = 0.5 * (header.AMIN + header.AMAX + b1)
+        for (let j = 0, jl = values.length; j < jl; ++j) {
+            values[j] = b1 * values[j] + b0
+        }
+    }
+
     const result: Schema.Ccp4File = { header, values };
     return Result.success(result);
 }

src/mol-math/geometry/spacegroup/construction.ts

@@ -26,10 +26,10 @@ interface Spacegroup {
 }
 
 namespace SpacegroupCell {
-    // Create a 'P 1' with cellsize [1, 1, 1]
+    /** Create a 'P 1' with cellsize [1, 1, 1] */
     export const Zero: SpacegroupCell = create('P 1', Vec3.create(1, 1, 1), Vec3.create(Math.PI / 2, Math.PI / 2, Math.PI / 2));
 
-    // True if 'P 1' with cellsize [1, 1, 1]
+    /** True if 'P 1' with cellsize [1, 1, 1] */
     export function isZero(cell: SpacegroupCell) {
         return cell.index === 0 && cell.size[0] === 1 && cell.size[1] === 1 && cell.size[1] === 1;
     }
@@ -37,7 +37,10 @@ namespace SpacegroupCell {
     // returns Zero cell if the spacegroup does not exist
     export function create(nameOrNumber: number | string | SpacegroupName, size: Vec3, anglesInRadians: Vec3): SpacegroupCell {
         const index = getSpacegroupIndex(nameOrNumber);
-        if (index < 0) return Zero;
+        if (index < 0) {
+            console.warn(`Unknown spacegroup '${nameOrNumber}', returning a 'P 1' with cellsize [1, 1, 1]`)
+            return Zero;
+        }
 
         const alpha = anglesInRadians[0];
         const beta = anglesInRadians[1];

src/mol-model/structure/model/formats/mmcif.ts

@@ -27,7 +27,7 @@ import { ChemicalComponent, ChemicalComponentMap } from '../properties/chemical-
 import { ComponentType, getMoleculeType } from '../types';
 
 import mmCIF_Format = Format.mmCIF
-import { SaccharideComponentMap, SaccharideComponent, SaccharidesSnfgMap, UnknownSaccharideComponent, SaccharideCompIdMap } from 'mol-model/structure/structure/carbohydrates/constants';
+import { SaccharideComponentMap, SaccharideComponent, SaccharidesSnfgMap, SaccharideCompIdMap } from 'mol-model/structure/structure/carbohydrates/constants';
 
 type AtomSite = mmCIF_Database['atom_site']
 

src/mol-model/structure/structure/carbohydrates/compute.ts

@@ -17,7 +17,7 @@ import { getAtomicMoleculeType, getPositionMatrix } from '../../util';
 import StructureElement from '../element';
 import Structure from '../structure';
 import Unit from '../unit';
-import { SaccharideCompIdMap, UnknownSaccharideComponent, SaccharideComponent, SaccharidesSnfgMap } from './constants';
+import { UnknownSaccharideComponent, SaccharideComponent } from './constants';
 import { CarbohydrateElement, CarbohydrateLink, Carbohydrates, CarbohydrateTerminalLink, PartialCarbohydrateElement } from './data';
 import { UnitRings, UnitRing } from '../unit/rings';
 import { ElementIndex } from '../../model/indexing';

src/mol-model/volume/data.ts

@@ -25,10 +25,7 @@ namespace VolumeData {
     export function getGridToCartesianTransform(volume: VolumeData) {
         const { data: { space } } = volume;
         const scale = Mat4.fromScaling(_scale, Vec3.div(Vec3.zero(), Box3D.size(Vec3.zero(), volume.fractionalBox), Vec3.ofArray(space.dimensions)));
-        const scale1 = Mat4.getScaling(Vec3.zero(), scale)
-        console.log('scale1', scale1)
         const translate = Mat4.fromTranslation(_translate, volume.fractionalBox.min);
-        console.log('translate1', translate)
         return Mat4.mul3(Mat4.zero(), volume.cell.fromFractional, translate, scale);
     }
 }

src/mol-model/volume/formats/ccp4.ts

@@ -7,33 +7,17 @@
 import { VolumeData } from '../data'
 import { Task } from 'mol-task';
 import { SpacegroupCell, Box3D } from 'mol-math/geometry';
-import { Tensor, Vec3, Mat4 } from 'mol-math/linear-algebra';
+import { Tensor, Vec3 } from 'mol-math/linear-algebra';
 import { Ccp4File } from 'mol-io/reader/ccp4/schema';
 import { degToRad } from 'mol-math/misc';
 
-// TODO implement voxelSize parameter
-
-// TODO support for rescaling of values
-// based on uglymol (https://github.com/uglymol/uglymol) by Marcin Wojdyr (wojdyr)
-// if the file was converted by mapmode2to0 - scale the data
-// if (intView[ 39 ] === -128 && intView[ 40 ] === 127) {
-//     // scaling f(x)=b1*x+b0 such that f(-128)=min and f(127)=max
-//     const b1 = (header.DMAX - header.DMIN) / 255.0
-//     const b0 = 0.5 * (header.DMIN + header.DMAX + b1)
-//     for (let j = 0, jl = data.length; j < jl; ++j) {
-//         data[ j ] = b1 * data[ j ] + b0
-//     }
-// }
-
-function volumeFromCcp4(source: Ccp4File): Task<VolumeData> {
+function volumeFromCcp4(source: Ccp4File, params?: { voxelSize?: Vec3 }): Task<VolumeData> {
     return Task.create<VolumeData>('Parse Volume Data', async ctx => {
         const { header, values } = source;
-        console.log('header', header)
-        const cell = SpacegroupCell.create(
-            header.ISPG,
-            Vec3.create(header.xLength, header.yLength, header.zLength),
-            Vec3.create(degToRad(header.alpha), degToRad(header.beta), degToRad(header.gamma))
-        );
+        const size = Vec3.create(header.xLength, header.yLength, header.zLength)
+        if (params && params.voxelSize) Vec3.mul(size, size, params.voxelSize)
+        const angles = Vec3.create(degToRad(header.alpha), degToRad(header.beta), degToRad(header.gamma))
+        const cell = SpacegroupCell.create(header.ISPG || 'P 1', size, angles);
 
         const axis_order_fast_to_slow = Vec3.create(header.MAPC - 1, header.MAPR - 1, header.MAPS - 1);
         const normalizeOrder = Tensor.convertToCanonicalAxisIndicesFastToSlow(axis_order_fast_to_slow);
@@ -41,30 +25,17 @@ function volumeFromCcp4(source: Ccp4File): Task<VolumeData> {
         const grid = [header.NX, header.NY, header.NZ];
         const extent = normalizeOrder([header.NC, header.NR, header.NS]);
 
-        const gridOrigin = normalizeOrder([header.NCSTART, header.NRSTART, header.NSSTART]);
-
-        // TODO: this is code from LiteMol, but not sure about the part based on originXYZ. The
-        // if (header.originX === 0.0 && header.originY === 0.0 && header.originZ === 0.0) {
-        //     gridOrigin = normalizeOrder([header.NCSTART, header.NRSTART, header.NSSTART]);
-        // } else {
-        //     // Use ORIGIN records rather than old n[xyz]start records
-        //     //   http://www2.mrc-lmb.cam.ac.uk/image2000.html
-        //     // XXX the ORIGIN field is only used by the EM community, and
-        //     //     has undefined meaning for non-orthogonal maps and/or
-        //     //     non-cubic voxels, etc.
-        //     gridOrigin = [header.originX, header.originY, header.originZ];
-        // }
+        let gridOrigin: number[]
+        if (header.originX === 0.0 && header.originY === 0.0 && header.originZ === 0.0) {
+            gridOrigin = normalizeOrder([header.NCSTART, header.NRSTART, header.NSSTART]);
+        } else {
+            // When available (e.g. in MRC files) use ORIGIN records instead of N[CRS]START
+            gridOrigin = [header.originX, header.originY, header.originZ];
+        }
 
         const origin_frac = Vec3.create(gridOrigin[0] / grid[0], gridOrigin[1] / grid[1], gridOrigin[2] / grid[2]);
         const dimensions_frac = Vec3.create(extent[0] / grid[0], extent[1] / grid[1], extent[2] / grid[2]);
 
-        console.log('length', header.xLength, header.yLength, header.zLength)
-        console.log('grid', grid)
-        console.log('extent', extent)
-        console.log('gridOrigin', gridOrigin)
-        console.log('origin_frac', origin_frac)
-        console.log('dimensions_frac', dimensions_frac)
-
         const space = Tensor.Space(extent, Tensor.invertAxisOrder(axis_order_fast_to_slow), header.MODE === 0 ? Int8Array : Float32Array);
         const data = Tensor.create(space, Tensor.Data1(values));
 
@@ -72,7 +43,7 @@ function volumeFromCcp4(source: Ccp4File): Task<VolumeData> {
         // Min/max/mean are reliable (based on LiteMol/DensityServer usage)
         // These, however, calculate sigma, so no data on that.
 
-        const vd = {
+        return {
             cell,
             fractionalBox: Box3D.create(origin_frac, Vec3.add(Vec3.zero(), origin_frac, dimensions_frac)),
             data,
@@ -83,12 +54,6 @@ function volumeFromCcp4(source: Ccp4File): Task<VolumeData> {
                 sigma: header.ARMS
             }
         };
-
-        const transform = VolumeData.getGridToCartesianTransform(vd);
-        const scale = Mat4.getScaling(Vec3.zero(), transform)
-        console.log('transform', transform)
-        console.log('scale', scale)
-        return vd
     });
 }