Skip to content
Snippets Groups Projects
Commit 99f104bc authored by David Sehnal's avatar David Sehnal
Browse files

Fixed CCP4 volume parser

parent 1a041254
No related branches found
No related tags found
No related merge requests found
...@@ -109,9 +109,7 @@ async function parseInternal(file: FileHandle, ctx: RuntimeContext): Promise<Res ...@@ -109,9 +109,7 @@ async function parseInternal(file: FileHandle, ctx: RuntimeContext): Promise<Res
} }
export function parse(file: FileHandle) { export function parse(file: FileHandle) {
return Task.create<Result<Schema.Ccp4File>>('Parse CCP4', async ctx => { return Task.create<Result<Schema.Ccp4File>>('Parse CCP4', ctx => parseInternal(file, ctx));
return await parseInternal(file, ctx);
});
} }
export default parse; export default parse;
\ No newline at end of file
...@@ -212,15 +212,21 @@ export namespace Tensor { ...@@ -212,15 +212,21 @@ export namespace Tensor {
return ret; return ret;
} }
export function getCanonicalAxisIndicesFastToSlow(order: number[]) { function reorder(xs: number[], indices: number[]) {
const ret: number[] = [];
for (let i = 0; i < xs.length; i++) ret[i] = xs[indices[i]];
return ret;
}
export function convertToCanonicalAxisIndicesFastToSlow(order: number[]) {
const indices = new Int32Array(order.length) as any as number[]; const indices = new Int32Array(order.length) as any as number[];
for (let i = 0; i < order.length; i++) indices[order[i]] = i; for (let i = 0; i < order.length; i++) indices[order[i]] = i;
return indices; return (xs: number[]) => reorder(xs, indices);
} }
export function getCanonicalAxisIndicesSlowToFast(order: number[]) { export function convertToCanonicalAxisIndicesSlowToFast(order: number[]) {
const indices = new Int32Array(order.length) as any as number[]; const indices = new Int32Array(order.length) as any as number[];
for (let i = 0; i < order.length; i++) indices[order[order.length - i - 1]] = i; for (let i = 0; i < order.length; i++) indices[order[order.length - i - 1]] = i;
return indices; return (xs: number[]) => reorder(xs, indices);
} }
} }
\ No newline at end of file
...@@ -32,21 +32,42 @@ function volumeFromCcp4(source: Ccp4File): Task<VolumeData> { ...@@ -32,21 +32,42 @@ function volumeFromCcp4(source: Ccp4File): Task<VolumeData> {
const cell = SpacegroupCell.create( const cell = SpacegroupCell.create(
header.ISPG, header.ISPG,
Vec3.create(header.xLength, header.yLength, header.zLength), Vec3.create(header.xLength, header.yLength, header.zLength),
Vec3.create(degToRad(header.alpha), degToRad(header.beta), degToRad(header.gamma) Vec3.create(degToRad(header.alpha), degToRad(header.beta), degToRad(header.gamma))
)) );
const origin = Vec3.create(header.originX, header.originY, header.originZ) const axis_order_fast_to_slow = Vec3.create(header.MAPC - 1, header.MAPR - 1, header.MAPS - 1);
const dimensions = Vec3.create(header.NX, header.NY, header.NZ) const normalizeOrder = Tensor.convertToCanonicalAxisIndicesFastToSlow(axis_order_fast_to_slow);
const axisOrder = Vec3.create(header.MAPC - 1, header.MAPR - 1, header.MAPS - 1)
const space = Tensor.Space(dimensions, axisOrder, header.MODE === 0 ? Int8Array : Float32Array); const grid = [header.NX, header.NY, header.NZ];
const extent = normalizeOrder([header.NC, header.NR, header.NS]);
const gridOrigin: number[] = 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];
// }
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, Tensor.invertAxisOrder(axis_order_fast_to_slow), header.MODE === 0 ? Int8Array : Float32Array);
const data = Tensor.create(space, Tensor.Data1(values)); const data = Tensor.create(space, Tensor.Data1(values));
// TODO Calculate stats? When to trust header data? // TODO Calculate stats? When to trust header data?
// Min/max/mean are reliable (based on LiteMol/DensityServer usage)
// These, however, calculate sigma, so no data on that.
return { return {
cell, cell,
fractionalBox: Box3D.create(origin, Vec3.add(Vec3.zero(), origin, dimensions)), fractionalBox: Box3D.create(origin_frac, Vec3.add(Vec3.zero(), origin_frac, dimensions_frac)),
data, data,
dataStats: { dataStats: {
min: header.AMIN, min: header.AMIN,
......
...@@ -10,10 +10,6 @@ import { Task } from 'mol-task'; ...@@ -10,10 +10,6 @@ import { Task } from 'mol-task';
import { SpacegroupCell, Box3D } from 'mol-math/geometry'; import { SpacegroupCell, Box3D } from 'mol-math/geometry';
import { Tensor, Vec3 } from 'mol-math/linear-algebra'; import { Tensor, Vec3 } from 'mol-math/linear-algebra';
function normalizeOrder(v: number[], indices: number[]) {
return Vec3.create(v[indices[0]], v[indices[1]], v[indices[2]])
}
function parseDensityServerData(source: DensityServer_Data_Database): Task<VolumeData> { function parseDensityServerData(source: DensityServer_Data_Database): Task<VolumeData> {
return Task.create<VolumeData>('Parse Volume Data', async ctx => { return Task.create<VolumeData>('Parse Volume Data', async ctx => {
const { volume_data_3d_info: info, volume_data_3d: values } = source; const { volume_data_3d_info: info, volume_data_3d: values } = source;
...@@ -25,17 +21,17 @@ function parseDensityServerData(source: DensityServer_Data_Database): Task<Volum ...@@ -25,17 +21,17 @@ function parseDensityServerData(source: DensityServer_Data_Database): Task<Volum
const axis_order_fast_to_slow = info.axis_order.value(0); const axis_order_fast_to_slow = info.axis_order.value(0);
const indices = Tensor.getCanonicalAxisIndicesFastToSlow(axis_order_fast_to_slow); const normalizeOrder = Tensor.convertToCanonicalAxisIndicesFastToSlow(axis_order_fast_to_slow);
// sample count is in "axis order" and needs to be reordered // sample count is in "axis order" and needs to be reordered
const sample_count = normalizeOrder(info.sample_count.value(0), indices); const sample_count = normalizeOrder(info.sample_count.value(0));
const tensorSpace = Tensor.Space(sample_count, Tensor.invertAxisOrder(axis_order_fast_to_slow), Float32Array); const tensorSpace = Tensor.Space(sample_count, Tensor.invertAxisOrder(axis_order_fast_to_slow), Float32Array);
const data = Tensor.create(tensorSpace, Tensor.Data1(values.values.toArray({ array: Float32Array }))); const data = Tensor.create(tensorSpace, Tensor.Data1(values.values.toArray({ array: Float32Array })));
// origin and dimensions are in "axis order" and need to be reordered // origin and dimensions are in "axis order" and need to be reordered
const origin = normalizeOrder(info.origin.value(0), indices) const origin = Vec3.ofArray(normalizeOrder(info.origin.value(0)))
const dimensions = normalizeOrder(info.dimensions.value(0), indices); const dimensions = Vec3.ofArray(normalizeOrder(info.dimensions.value(0)));
return { return {
cell, cell,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment