VolumeData, mesh, tweaks

src/apps/render-test/utils/mcubes.ts

@@ -5,7 +5,7 @@
  */
 
 import { Run } from 'mol-task'
-import { compute } from 'mol-geo/util/marching-cubes/algorithm'
+import { computeMarchingCubes } from 'mol-geo/util/marching-cubes/algorithm'
 import { Mesh } from 'mol-geo/shape/mesh'
 import { Tensor, Mat4, Vec3 } from 'mol-math/linear-algebra'
 
@@ -38,7 +38,7 @@ export default async function computeSurface(f: (x: number, y: number, z: number
     const min = Vec3.create(-1.1, -1.1, -1.1), max = Vec3.create(1.1, 1.1, 1.1);
 
     fillField(field, f, min, max);
-    const surface = await Run(compute({
+    const surface = await Run(computeMarchingCubes({
         scalarField: field,
         isoLevel: 0,
         oldSurface: data ? data.surface : void 0

src/apps/structure-info/volume.ts

@@ -1,35 +1,77 @@
-
 /**
  * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
  */
 
+import * as fs from 'fs'
 import * as argparse from 'argparse'
-import { VolumeData, parseDensityServerData } from 'mol-model/volume'
+import * as util from 'util'
+
+import { VolumeData, parseDensityServerData, VolumeIsoValue } from 'mol-model/volume'
 import { downloadCif } from './helpers'
 import CIF from 'mol-io/reader/cif'
 import { DensityServer_Data_Database } from 'mol-io/reader/cif/schema/density-server';
 import { Run } from 'mol-task';
 import { Table } from 'mol-data/db';
+import { computeVolumeMesh } from 'mol-geo/representation/volume/Mesh';
+import { StringBuilder } from 'mol-util';
+
+require('util.promisify').shim();
+const writeFileAsync = util.promisify(fs.writeFile);
 
 type Volume = { source: DensityServer_Data_Database, volume: VolumeData }
 
 async function getVolume(url: string): Promise<Volume> {
-    const cif = await downloadCif(url, false);
+    const cif = await downloadCif(url, true);
     const data = CIF.schema.densityServer(cif.blocks[1]);
     return { source: data, volume: await Run(parseDensityServerData(data)) };
 }
 
-function print(volume: Volume) {
-    const { volume_data_3d_info } = volume.source;
+function print(data: Volume) {
+    const { volume_data_3d_info } = data.source;
     const row = Table.getRow(volume_data_3d_info, 0);
     console.log(row);
+    console.log(data.volume.cell);
+    console.log(data.volume.dataStats);
+    console.log(data.volume.fractionalBox);
+}
+
+async function doMesh(data: Volume, filename: string) {
+    const mesh = await Run(computeVolumeMesh(data.volume, VolumeIsoValue.relative(data.volume.dataStats, 1.5)));
+    console.log({ vc: mesh.vertexCount, tc: mesh.triangleCount });
+
+    // Export the mesh in OBJ format.
+    const { vertexCount, triangleCount } = mesh;
+
+    const vs = mesh.vertexBuffer.ref.value;
+    const ts = mesh.indexBuffer.ref.value;
+
+    const obj = StringBuilder.create();
+    for (let i = 0; i < vertexCount; i++) {
+        StringBuilder.write(obj, 'v ');
+        StringBuilder.writeFloat(obj, vs[3 * i + 0], 100);
+        StringBuilder.whitespace1(obj);
+        StringBuilder.writeFloat(obj, vs[3 * i + 1], 100);
+        StringBuilder.whitespace1(obj);
+        StringBuilder.writeFloat(obj, vs[3 * i + 2], 100);
+        StringBuilder.newline(obj);
+    }
+    for (let i = 0; i < triangleCount; i++) {
+        StringBuilder.write(obj, 'f ');
+        StringBuilder.writeIntegerAndSpace(obj, ts[3 * i + 0] + 1);
+        StringBuilder.writeIntegerAndSpace(obj, ts[3 * i + 1] + 1);
+        StringBuilder.writeInteger(obj, ts[3 * i + 2] + 1);
+        StringBuilder.newline(obj);
+    }
+
+    await writeFileAsync(filename, StringBuilder.getString(obj));
 }
 
-async function run(url: string) {
+async function run(url: string, meshFilename: string) {
     const volume = await getVolume(url);
     print(volume);
+    await doMesh(volume, meshFilename);
 }
 
 const parser = new argparse.ArgumentParser({
@@ -39,9 +81,14 @@ description: 'Info about VolumeData from mol-model module'
 parser.addArgument([ '--emdb', '-e' ], {
     help: 'EMDB id, for example 8116',
 });
+parser.addArgument([ '--mesh' ], {
+    help: 'Mesh filename',
+    required: true
+});
 interface Args {
-    emdb?: string
+    emdb?: string,
+    mesh: string
 }
 const args: Args = parser.parseArgs();
 
-run(`https://webchem.ncbr.muni.cz/DensityServer/em/emd-${args.emdb}/cell?detail=1&encoding=cif`);
+run(`https://webchem.ncbr.muni.cz/DensityServer/em/emd-${args.emdb}/cell?detail=4`, args.mesh);

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

@@ -0,0 +1,29 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import { VolumeData, VolumeIsoValue } from 'mol-model/volume'
+import { Task } from 'mol-task'
+import { Mesh } from '../../shape/mesh';
+import { computeMarchingCubes } from '../../util/marching-cubes/algorithm';
+
+function computeVolumeMesh(volume: VolumeData, isoValue: VolumeIsoValue) {
+    return Task.create<Mesh>('Volume Surface', async ctx => {
+        ctx.update({ message: 'Marching cubes...' });
+
+        const mesh = await ctx.runChild(computeMarchingCubes({
+            isoLevel: VolumeIsoValue.toAbsolute(isoValue).absoluteValue,
+            scalarField: volume.data
+        }));
+
+        const transform = VolumeData.getGridToCartesianTransform(volume);
+        ctx.update({ message: 'Transforming mesh...' });
+        Mesh.transformImmediate(mesh, transform);
+
+        return mesh;
+    });
+}
+
+export { computeVolumeMesh }

src/mol-geo/util/marching-cubes/algorithm.ts

@@ -26,7 +26,7 @@ export interface MarchingCubesParameters {
     oldSurface?: Mesh
 }
 
-export function compute(parameters: MarchingCubesParameters) {
+export function computeMarchingCubes(parameters: MarchingCubesParameters) {
     return Task.create('Marching Cubes', async ctx => {
         let comp = new MarchingCubesComputation(parameters, ctx);
         return await comp.run();

src/mol-math/geometry/primitives/box3d.ts

@@ -12,7 +12,7 @@ interface Box3D { min: Vec3, max: Vec3 }
 
 namespace Box3D {
     export function create(min: Vec3, max: Vec3): Box3D { return { min, max }; }
-    
+
     export function computeBounding(data: PositionData): Box3D {
         const min = [Number.MAX_VALUE, Number.MAX_VALUE, Number.MAX_VALUE];
         const max = [-Number.MAX_VALUE, -Number.MAX_VALUE, -Number.MAX_VALUE];

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

@@ -15,8 +15,7 @@ interface SpacegroupCell {
     /** Transfrom cartesian -> fractional coordinates within the cell */
     readonly toFractional: Mat4,
     /** Transfrom fractional coordinates within the cell -> cartesian */
-    readonly fromFractional: Mat4,
-    readonly fractionalBasis: Vec3[]
+    readonly fromFractional: Mat4
 }
 
 interface Spacegroup {
@@ -59,14 +58,8 @@ namespace SpacegroupCell {
         ]);
         const toFractional = Mat4.invert(Mat4.zero(), fromFractional)!;
 
-        const fractionalBasis = [
-            Vec3.transformMat4(Vec3.zero(), Vec3.create(1, 0, 0), toFractional),
-            Vec3.transformMat4(Vec3.zero(), Vec3.create(0, 1, 0), toFractional),
-            Vec3.transformMat4(Vec3.zero(), Vec3.create(0, 0, 1), toFractional)
-        ];
-
         const num = typeof nameOrNumber === 'number' ? nameOrNumber : SpacegroupNumbers[nameOrNumber];
-        return { number: num, size, anglesInRadians, toFractional, fromFractional, fractionalBasis };
+        return { number: num, size, anglesInRadians, toFractional, fromFractional };
     }
 }
 

src/mol-model/volume/_spec/data.spec.ts

@@ -1,6 +0,0 @@
-/**
- * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
- *
- * @author David Sehnal <david.sehnal@gmail.com>
- */
-

src/mol-model/volume/data.ts

@@ -30,4 +30,36 @@ namespace VolumeData {
     }
 }
 
-export { VolumeData }
+type VolumeIsoValue = VolumeIsoValue.Absolute | VolumeIsoValue.Relative
+
+namespace VolumeIsoValue {
+    export type Relative = Readonly<{ kind: 'relative', stats: VolumeData['dataStats'], relativeValue: number }>
+    export type Absolute = Readonly<{ kind: 'absolute', stats: VolumeData['dataStats'], absoluteValue: number }>
+
+    export function absolute(stats: VolumeData['dataStats'], value: number): Absolute { return { kind: 'absolute', stats, absoluteValue: value }; }
+    export function relative(stats: VolumeData['dataStats'], value: number): Relative { return { kind: 'relative', stats, relativeValue: value }; }
+
+    export function toAbsolute(value: VolumeIsoValue): Absolute {
+        if (value.kind === 'absolute') return value;
+
+        const { mean, sigma } = value.stats
+        return {
+            kind: 'absolute',
+            stats: value.stats,
+            absoluteValue: value.relativeValue * sigma + mean
+        }
+    }
+
+    export function toRelative(value: VolumeIsoValue): Relative {
+        if (value.kind === 'relative') return value;
+
+        const { mean, sigma } = value.stats
+        return {
+            kind: 'relative',
+            stats: value.stats,
+            relativeValue: (mean - value.absoluteValue) / sigma
+        }
+    }
+}
+
+export { VolumeData, VolumeIsoValue }

src/mol-model/volume/formats/density-server.ts

@@ -13,7 +13,9 @@ import { Tensor, Vec3 } from 'mol-math/linear-algebra';
 function parseDensityServerData(source: DensityServer_Data_Database): Task<VolumeData> {
     return Task.create<VolumeData>('Parse Volume Data', async ctx => {
         const { volume_data_3d_info: info, volume_data_3d: values } = source;
-        const cell = SpacegroupCell.create(info.spacegroup_number.value(0), Vec3.ofArray(info.spacegroup_cell_size.value(0)), Vec3.ofArray(info.spacegroup_cell_angles.value(0)));
+        const cell = SpacegroupCell.create(info.spacegroup_number.value(0),
+            Vec3.ofArray(info.spacegroup_cell_size.value(0)),
+            Vec3.scale(Vec3.zero(), Vec3.ofArray(info.spacegroup_cell_angles.value(0)), Math.PI / 180));
 
         const tensorSpace = Tensor.Space(info.sample_count.value(0), info.axis_order.value(0), Float32Array);
         const data = Tensor.create(tensorSpace, Tensor.Data1(values.values.toArray()));