ccp4/mrc handling improvements

- read mode 6 (uint16)
- correctly normalize mrc origin offset
- use heuristic to determine endianess if not given
- handle erroneous spacegroup
- calculate data stats if not available
- add offset param

src/mol-io/common/typed-array.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2018-2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * Taken/adapted from DensityServer (https://github.com/dsehnal/DensityServer)
  *
@@ -10,15 +10,16 @@
 import { FileHandle } from '../../mol-io/common/file-handle';
 import { SimpleBuffer } from '../../mol-io/common/simple-buffer';
 
-export type TypedArrayValueType = 'float32' | 'int8' | 'int16'
+export type TypedArrayValueType = 'float32' | 'int8' | 'int16' | 'uint16'
 
 export namespace TypedArrayValueType {
     export const Float32: TypedArrayValueType = 'float32';
     export const Int8: TypedArrayValueType = 'int8';
     export const Int16: TypedArrayValueType = 'int16';
+    export const Uint16: TypedArrayValueType = 'uint16';
 }
 
-export type TypedArrayValueArray = Float32Array | Int8Array | Int16Array
+export type TypedArrayValueArray = Float32Array | Int8Array | Int16Array | Uint16Array
 
 export interface TypedArrayBufferContext {
     type: TypedArrayValueType,
@@ -31,12 +32,14 @@ export interface TypedArrayBufferContext {
 export function getElementByteSize(type: TypedArrayValueType) {
     if (type === TypedArrayValueType.Float32) return 4;
     if (type === TypedArrayValueType.Int16) return 2;
+    if (type === TypedArrayValueType.Uint16) return 2;
     return 1;
 }
 
 export function makeTypedArray(type: TypedArrayValueType, buffer: ArrayBuffer, byteOffset = 0, length?: number): TypedArrayValueArray {
     if (type === TypedArrayValueType.Float32) return new Float32Array(buffer, byteOffset, length);
     if (type === TypedArrayValueType.Int16) return new Int16Array(buffer, byteOffset, length);
+    if (type === TypedArrayValueType.Uint16) return new Uint16Array(buffer, byteOffset, length);
     return new Int8Array(buffer, byteOffset, length);
 }
 
@@ -45,6 +48,7 @@ export function createTypedArray(type: TypedArrayValueType, size: number) {
         case TypedArrayValueType.Float32: return new Float32Array(new ArrayBuffer(4 * size));
         case TypedArrayValueType.Int8: return new Int8Array(new ArrayBuffer(1 * size));
         case TypedArrayValueType.Int16: return new Int16Array(new ArrayBuffer(2 * size));
+        case TypedArrayValueType.Uint16: return new Uint16Array(new ArrayBuffer(2 * size));
     }
     throw Error(`${type} is not a supported value format.`);
 }

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

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2018-2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
@@ -27,10 +27,14 @@ export async function readCcp4Header(file: FileHandle): Promise<{ header: Ccp4He
     // 54  MACHST      Machine stamp indicating machine type which wrote file
     //                 17 and 17 for big-endian or 68 and 65 for little-endian
     const MACHST = [ buffer.readUInt8(53 * 4), buffer.readUInt8(53 * 4 + 1) ]
-    let littleEndian = true
-    // found MRC files that don't have the MACHST stamp set and are big-endian
-    if (MACHST[0] !== 68 && MACHST[1] !== 65) {
+    let littleEndian = false
+    if (MACHST[0] === 68 && MACHST[1] === 65) {
+        littleEndian = true;
+    } else if (MACHST[0] === 17 && MACHST[1] === 17) {
         littleEndian = false;
+    } else {
+        const modeLE = buffer.readInt32LE(3 * 4)
+        if (modeLE <= 16) littleEndian = true;
     }
 
     const readInt = littleEndian ? (o: number) => buffer.readInt32LE(o * 4) : (o: number) => buffer.readInt32BE(o * 4)
@@ -117,11 +121,15 @@ export async function readCcp4Slices(header: Ccp4Header, buffer: TypedArrayBuffe
 
 function getCcp4DataType(mode: number) {
     switch (mode) {
-        case 2: return TypedArrayValueType.Float32
-        case 1: return TypedArrayValueType.Int16
         case 0: return TypedArrayValueType.Int8
+        case 1: return TypedArrayValueType.Int16
+        case 2: return TypedArrayValueType.Float32
+        case 3: throw new Error('mode 3 unsupported, complex 16-bit integers')
+        case 4: throw new Error('mode 4 unsupported, complex 32-bit reals')
+        case 6: TypedArrayValueType.Uint16
+        case 16: throw new Error('mode 16 unsupported, unsigned char * 3 (for rgb data, non-standard)')
     }
-    throw new Error(`ccp4 mode '${mode}' unsupported`);
+    throw new Error(`unknown mode '${mode}'`);
 }
 
 /** check if the file was converted by mapmode2to0, see https://github.com/uglymol/uglymol */

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

@@ -115,5 +115,5 @@ export interface Ccp4Header {
  */
 export interface Ccp4File {
     header: Ccp4Header
-    values: Float32Array | Int16Array | Int8Array
+    values: Float32Array | Int16Array | Int8Array | Uint16Array
 }

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

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
@@ -12,16 +12,19 @@ import { Ccp4File, Ccp4Header } from '../../mol-io/reader/ccp4/schema';
 import { degToRad } from '../../mol-math/misc';
 import { getCcp4ValueType } from '../../mol-io/reader/ccp4/parser';
 import { TypedArrayValueType } from '../../mol-io/common/typed-array';
+import { arrayMin, arrayRms, arrayMean, arrayMax } from '../../mol-util/array';
 
 /** When available (e.g. in MRC files) use ORIGIN records instead of N[CRS]START */
-export function getCcp4Origin(header: Ccp4Header) {
-    let gridOrigin: number[]
+export function getCcp4Origin(header: Ccp4Header): Vec3 {
     if (header.originX === 0.0 && header.originY === 0.0 && header.originZ === 0.0) {
-        gridOrigin = [header.NCSTART, header.NRSTART, header.NSSTART];
+        return Vec3.create(header.NCSTART, header.NRSTART, header.NSSTART)
     } else {
-        gridOrigin = [header.originX, header.originY, header.originZ];
+        return Vec3.create(
+            header.originX / (header.xLength / header.NX),
+            header.originY / (header.yLength / header.NY),
+            header.originZ / (header.zLength / header.NZ)
+        )
     }
-    return gridOrigin
 }
 
 function getTypedArrayCtor(header: Ccp4Header) {
@@ -30,25 +33,28 @@ function getTypedArrayCtor(header: Ccp4Header) {
         case TypedArrayValueType.Float32: return Float32Array;
         case TypedArrayValueType.Int8: return Int8Array;
         case TypedArrayValueType.Int16: return Int16Array;
+        case TypedArrayValueType.Uint16: return Uint16Array;
     }
     throw Error(`${valueType} is not a supported value format.`);
 }
 
-export function volumeFromCcp4(source: Ccp4File, params?: { voxelSize?: Vec3 }): Task<VolumeData> {
+export function volumeFromCcp4(source: Ccp4File, params?: { voxelSize?: Vec3, offset?: Vec3 }): Task<VolumeData> {
     return Task.create<VolumeData>('Create Volume Data', async ctx => {
         const { header, values } = source;
-        console.log({ header, values })
         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 spacegroup = header.ISPG > 65536 ? 0 : header.ISPG
+        const cell = SpacegroupCell.create(spacegroup || '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);
 
         const grid = [header.NX, header.NY, header.NZ];
         const extent = normalizeOrder([header.NC, header.NR, header.NS]);
-        const gridOrigin = normalizeOrder(getCcp4Origin(header));
+        const origin = getCcp4Origin(header)
+        if (params?.offset) Vec3.add(origin, origin, params.offset)
+        const gridOrigin = normalizeOrder(origin);
 
         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]);
@@ -65,10 +71,10 @@ export function volumeFromCcp4(source: Ccp4File, params?: { voxelSize?: Vec3 }):
             fractionalBox: Box3D.create(origin_frac, Vec3.add(Vec3.zero(), origin_frac, dimensions_frac)),
             data,
             dataStats: {
-                min: header.AMIN,
-                max: header.AMAX,
-                mean: header.AMEAN,
-                sigma: header.ARMS
+                min: isNaN(header.AMIN) ? arrayMin(values) : header.AMIN,
+                max: isNaN(header.AMAX) ? arrayMax(values) : header.AMAX,
+                mean: isNaN(header.AMEAN) ? arrayMean(values) : header.AMEAN,
+                sigma: (isNaN(header.ARMS) || header.ARMS === 0) ? arrayRms(values) : header.ARMS
             }
         };
     });

src/mol-plugin/state/actions/volume.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2018-2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
@@ -28,7 +28,7 @@ export const Ccp4Provider: DataFormatProvider<any> = {
     },
     getDefaultBuilder: (ctx: PluginContext, data: StateBuilder.To<PluginStateObject.Data.Binary>, options: DataFormatBuilderOptions, state: State) => {
         return Task.create('CCP4/MRC/BRIX default builder', async taskCtx => {
-            let tree: StateBuilder.To<any> = data.apply(StateTransforms.Data.ParseCcp4)
+            let tree: StateBuilder.To<any> = data.apply(StateTransforms.Data.ParseCcp4, {}, { state: { isGhost: true } })
                 .apply(StateTransforms.Volume.VolumeFromCcp4)
             if (options.visuals) {
                 tree = tree.apply(StateTransforms.Representation.VolumeRepresentation3D)
@@ -48,7 +48,7 @@ export const Dsn6Provider: DataFormatProvider<any> = {
     },
     getDefaultBuilder: (ctx: PluginContext, data: StateBuilder.To<PluginStateObject.Data.Binary>, options: DataFormatBuilderOptions, state: State) => {
         return Task.create('DSN6/BRIX default builder', async taskCtx => {
-            let tree: StateBuilder.To<any> = data.apply(StateTransforms.Data.ParseDsn6)
+            let tree: StateBuilder.To<any> = data.apply(StateTransforms.Data.ParseDsn6, {}, { state: { isGhost: true } })
                 .apply(StateTransforms.Volume.VolumeFromDsn6)
             if (options.visuals) {
                 tree = tree.apply(StateTransforms.Representation.VolumeRepresentation3D)

src/mol-plugin/state/transforms/volume.ts

@@ -25,7 +25,8 @@ const VolumeFromCcp4 = PluginStateTransform.BuiltIn({
     to: SO.Volume.Data,
     params(a) {
         return {
-            voxelSize: PD.Vec3(Vec3.create(1, 1, 1))
+            voxelSize: PD.Vec3(Vec3.create(1, 1, 1)),
+            offset: PD.Vec3(Vec3.create(0, 0, 0))
         };
     }
 })({