wip, dsn6/brix parsing

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

@@ -121,10 +121,9 @@ async function parseInternal(file: FileHandle, ctx: RuntimeContext): Promise<Res
 }
 
 export function parseFile(file: FileHandle) {
-    return Task.create<Result<Ccp4File>>('Parse CCP4', ctx => parseInternal(file, ctx));
+    return Task.create<Result<Ccp4File>>('Parse CCP4/MRC', ctx => parseInternal(file, ctx));
 }
 
 export function parse(buffer: Uint8Array) {
-    const file = FileHandle.fromBuffer(buffer)
-    return Task.create<Result<Ccp4File>>('Parse CCP4', ctx => parseInternal(file, ctx));
+    return parseFile(FileHandle.fromBuffer(buffer))
 }

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

@@ -0,0 +1,126 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { Task, RuntimeContext } from 'mol-task';
+import { Dsn6File, Dsn6Header } from './schema'
+import Result from '../result'
+import { FileHandle } from '../../common/file-handle';
+
+function parseBrixHeader(str: string): Dsn6Header {
+    return {
+        xStart: parseInt(str.substr(10, 5)),
+        yStart: parseInt(str.substr(15, 5)),
+        zStart: parseInt(str.substr(20, 5)),
+        xExtent: parseInt(str.substr(32, 5)),
+        yExtent: parseInt(str.substr(38, 5)),
+        zExtent: parseInt(str.substr(42, 5)),
+        xRate: parseInt(str.substr(52, 5)),
+        yRate: parseInt(str.substr(58, 5)),
+        zRate: parseInt(str.substr(62, 5)),
+        xlen: parseFloat(str.substr(73, 10)),
+        ylen: parseFloat(str.substr(83, 10)),
+        zlen: parseFloat(str.substr(93, 10)),
+        alpha: parseFloat(str.substr(103, 10)),
+        beta: parseFloat(str.substr(113, 10)),
+        gamma: parseFloat(str.substr(123, 10)),
+        divisor: parseFloat(str.substr(138, 12)) / 100,
+        summand: parseInt(str.substr(155, 8)),
+        sigma: parseFloat(str.substr(170, 12)) * 100
+    }
+}
+
+function parseDsn6Header(int: Int16Array): Dsn6Header {
+    const factor = 1 / int[ 17 ]
+    return {
+        xStart: int[ 0 ],
+        yStart: int[ 1 ],
+        zStart: int[ 2 ],
+        xExtent: int[ 3 ],
+        yExtent: int[ 4 ],
+        zExtent: int[ 5 ],
+        xRate: int[ 6 ],
+        yRate: int[ 7 ],
+        zRate: int[ 8 ],
+        xlen: int[ 9 ] * factor,
+        ylen: int[ 10 ] * factor,
+        zlen: int[ 11 ] * factor,
+        alpha: int[ 12 ] * factor,
+        beta: int[ 13 ] * factor,
+        gamma: int[ 14 ] * factor,
+        divisor: int[ 15 ] / 100,
+        summand: int[ 16 ],
+        sigma: undefined
+    }
+}
+
+async function parseInternal(file: FileHandle, ctx: RuntimeContext): Promise<Result<Dsn6File>> {
+    await ctx.update({ message: 'Parsing CCP4 file...' });
+
+    const { buffer } = await file.readBuffer(0, file.length)
+    const bin = buffer.buffer
+
+    const intView = new Int16Array(bin)
+    const byteView = new Uint8Array(bin)
+    const brixStr = String.fromCharCode.apply(null, byteView.subarray(0, 512))
+    const isBrix = brixStr.startsWith(':-)')
+
+    const header = isBrix ? parseBrixHeader(brixStr) : parseDsn6Header(intView)
+    const { divisor, summand } = header
+
+    if (!isBrix) {
+        // for DSN6, swap byte order when big endian
+        if (intView[ 18 ] !== 100) {
+            for (let i = 0, n = intView.length; i < n; ++i) {
+              const val = intView[ i ]
+                intView[ i ] = ((val & 0xff) << 8) | ((val >> 8) & 0xff)
+            }
+        }
+    }
+    const values = new Float32Array(header.xExtent * header.yExtent * header.zExtent)
+
+    let offset = 512
+    const xBlocks = Math.ceil(header.xExtent / 8)
+    const yBlocks = Math.ceil(header.yExtent / 8)
+    const zBlocks = Math.ceil(header.zExtent / 8)
+
+    // loop over blocks
+    for (let zz = 0; zz < zBlocks; ++zz) {
+        for (let yy = 0; yy < yBlocks; ++yy) {
+            for (let xx = 0; xx < xBlocks; ++xx) {
+                // loop inside block
+                for (let k = 0; k < 8; ++k) {
+                    const z = 8 * zz + k
+                    for (let j = 0; j < 8; ++j) {
+                        const y = 8 * yy + j
+                        for (let i = 0; i < 8; ++i) {
+                            const x = 8 * xx + i
+                            // check if remaining slice-part contains values
+                            if (x < header.xExtent && y < header.yExtent && z < header.zExtent) {
+                                const idx = ((((x * header.yExtent) + y) * header.zExtent) + z)
+                                values[ idx ] = (byteView[ offset ] - summand) / divisor
+                                ++offset
+                            } else {
+                                offset += 8 - i
+                                break
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    const result: Dsn6File = { header, values };
+    return Result.success(result);
+}
+
+export function parseFile(file: FileHandle) {
+    return Task.create<Result<Dsn6File>>('Parse DSN6/BRIX', ctx => parseInternal(file, ctx));
+}
+
+export function parse(buffer: Uint8Array) {
+    return parseFile(FileHandle.fromBuffer(buffer))
+}

src/mol-io/reader/dsn6/schema.ts

@@ -0,0 +1,44 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+export interface Dsn6Header {
+    /** the origin of the map in grid units, along X, Y, Z */
+    xStart: number
+    yStart: number
+    zStart: number
+    /** the extent (size) of the map, along X, Y, Z, in grid units */
+    xExtent: number
+    yExtent: number
+    zExtent: number
+    /** number of grid points along the whole unit cell, along X, Y, Z */
+    xRate: number
+    yRate: number
+    zRate: number
+    /** Unit cell parameters */
+    xlen: number
+    ylen: number
+    zlen: number
+    alpha: number
+    beta: number
+    gamma: number
+    /**
+     * Constants that bring the electron density from byte to normal scale.
+     * They are calculated like this: prod = 255.0/(rhomax-rhomin), plus = -rhomin*prod.
+     */
+    divisor: number
+    summand: number
+    /** Rms deviation of electron density map (only given in BRIX but not in DSN6) */
+    sigma: number | undefined
+}
+
+/**
+ * DSN6 http://www.uoxray.uoregon.edu/tnt/manual/node104.html
+ * BRIX http://svn.cgl.ucsf.edu/svn/chimera/trunk/libs/VolumeData/dsn6/brix-1.html
+ */
+export interface Dsn6File {
+    header: Dsn6Header
+    values: Float32Array | Int8Array
+}

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

@@ -12,7 +12,7 @@ import { Ccp4File } from 'mol-io/reader/ccp4/schema';
 import { degToRad } from 'mol-math/misc';
 
 function volumeFromCcp4(source: Ccp4File, params?: { voxelSize?: Vec3 }): Task<VolumeData> {
-    return Task.create<VolumeData>('Parse Volume Data', async ctx => {
+    return Task.create<VolumeData>('Create Volume Data', async ctx => {
         const { header, values } = source;
         const size = Vec3.create(header.xLength, header.yLength, header.zLength)
         if (params && params.voxelSize) Vec3.mul(size, size, params.voxelSize)

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

@@ -0,0 +1,48 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { VolumeData } from '../data'
+import { Task } from 'mol-task';
+import { SpacegroupCell, Box3D } from 'mol-math/geometry';
+import { Tensor, Vec3 } from 'mol-math/linear-algebra';
+import { degToRad } from 'mol-math/misc';
+import { Dsn6File } from 'mol-io/reader/dsn6/schema';
+import { arrayMin, arrayMax, arrayMean, arrayRms } from 'mol-util/array';
+
+function volumeFromDsn6(source: Dsn6File, params?: { voxelSize?: Vec3 }): Task<VolumeData> {
+    return Task.create<VolumeData>('Create Volume Data', async ctx => {
+        const { header, values } = source;
+        const size = Vec3.create(header.xlen, header.ylen, header.zlen)
+        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('P 1', size, angles)
+
+        const grid = [header.xRate, header.yRate, header.zRate];
+        const extent = [header.xExtent, header.yExtent, header.zExtent];
+
+        const gridOrigin = [header.xStart, header.yStart, header.zStart];
+
+        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]);
+
+        const space = Tensor.Space(extent, [2, 1, 0], Float32Array);
+        const data = Tensor.create(space, Tensor.Data1(values));
+
+        return {
+            cell,
+            fractionalBox: Box3D.create(origin_frac, Vec3.add(Vec3.zero(), origin_frac, dimensions_frac)),
+            data,
+            dataStats: {
+                min: arrayMin(values),
+                max: arrayMax(values),
+                mean: arrayMean(values),
+                sigma: header.sigma !== undefined ? header.sigma : arrayRms(values)
+            }
+        };
+    });
+}
+
+export { volumeFromDsn6 }