Browse Source

added mc that outputs contour lines

Alexander Rose 6 years ago
parent
commit
ed31e06c06

+ 23 - 0
src/mol-geo/geometry/lines/lines.ts

@@ -16,6 +16,8 @@ import { TransformData } from '../transform-data';
 import { LocationIterator } from '../../util/location-iterator';
 import { SizeThemeProps } from 'mol-view/theme/size';
 import { LinesValues } from 'mol-gl/renderable/lines';
+import { Mesh } from '../mesh/mesh';
+import { LinesBuilder } from './lines-builder';
 
 /** Wide line */
 export interface Lines {
@@ -52,6 +54,27 @@ export namespace Lines {
         }
     }
 
+    export function fromMesh(mesh: Mesh, lines?: Lines) {
+        const vb = mesh.vertexBuffer.ref.value
+        const ib = mesh.indexBuffer.ref.value
+        const gb = mesh.groupBuffer.ref.value
+
+        const builder = LinesBuilder.create(mesh.triangleCount * 3, mesh.triangleCount / 10, lines)
+
+        // TODO avoid duplicate lines
+        for (let i = 0, il = mesh.triangleCount * 3; i < il; i += 3) {
+            const i0 = ib[i], i1 = ib[i + 1], i2 = ib[i + 2];
+            const x0 = vb[i0 * 3], y0 = vb[i0 * 3 + 1], z0 = vb[i0 * 3 + 2];
+            const x1 = vb[i1 * 3], y1 = vb[i1 * 3 + 1], z1 = vb[i1 * 3 + 2];
+            const x2 = vb[i2 * 3], y2 = vb[i2 * 3 + 1], z2 = vb[i2 * 3 + 2];
+            builder.add(x0, y0, z0, x1, y1, z1, gb[i0])
+            builder.add(x0, y0, z0, x2, y2, z2, gb[i0])
+            builder.add(x1, y1, z1, x2, y2, z2, gb[i1])
+        }
+
+        return builder.getLines();
+    }
+
     export function transformImmediate(line: Lines, t: Mat4) {
         transformRangeImmediate(line, t, 0, line.lineCount)
     }

+ 3 - 1
src/mol-geo/representation/structure/representation/surface.ts

@@ -13,12 +13,14 @@ import { Loci } from 'mol-model/loci';
 import { PickingId } from '../../../geometry/picking';
 import { Task } from 'mol-task';
 import { DefaultGaussianWireframeProps, GaussianWireframeVisual } from '../visual/gaussian-surface-wireframe';
+import { ColorThemeProps } from 'mol-view/theme/color';
 
 export const DefaultSurfaceProps = {
     ...DefaultGaussianSurfaceProps,
     ...DefaultGaussianWireframeProps,
 
-    visuals: { surface: true, wireframe: false }
+    visuals: { surface: true, wireframe: false },
+    colorTheme: { name: 'element-symbol' } as ColorThemeProps,
 }
 export type SurfaceProps = typeof DefaultSurfaceProps
 

+ 5 - 5
src/mol-geo/representation/structure/visual/gaussian-surface-mesh.ts

@@ -10,19 +10,19 @@ import { RuntimeContext } from 'mol-task'
 import { Mesh } from '../../../geometry/mesh/mesh';
 import { UnitsMeshVisual, DefaultUnitsMeshProps } from '../units-visual';
 import { StructureElementIterator, getElementLoci, markElement } from './util/element';
-import { computeMarchingCubes } from '../../../util/marching-cubes/algorithm';
+import { computeMarchingCubesMesh } from '../../../util/marching-cubes/algorithm';
 import { DefaultGaussianDensityProps, GaussianDensityProps } from 'mol-model/structure/structure/unit/gaussian-density';
 
 async function createGaussianSurfaceMesh(ctx: RuntimeContext, unit: Unit, structure: Structure, props: GaussianDensityProps, mesh?: Mesh): Promise<Mesh> {
     const { smoothness } = props
     const { transform, field, idField } = await unit.computeGaussianDensity(props, ctx)
 
-    const surface = await computeMarchingCubes({
+    const params = {
         isoLevel: Math.exp(-smoothness),
         scalarField: field,
-        idField,
-        oldSurface: mesh
-    }).runAsChild(ctx)
+        idField
+    }
+    const surface = await computeMarchingCubesMesh(params, mesh).runAsChild(ctx)
 
     Mesh.transformImmediate(surface, transform)
     Mesh.computeNormalsImmediate(surface)

+ 6 - 26
src/mol-geo/representation/structure/visual/gaussian-surface-wireframe.ts

@@ -7,46 +7,26 @@
 import { Unit, Structure } from 'mol-model/structure';
 import { UnitsVisual, VisualUpdateState } from '..';
 import { RuntimeContext } from 'mol-task'
-import { Mesh } from '../../../geometry/mesh/mesh';
 import { UnitsLinesVisual, DefaultUnitsLinesProps } from '../units-visual';
 import { StructureElementIterator, getElementLoci, markElement } from './util/element';
-import { computeMarchingCubes } from '../../../util/marching-cubes/algorithm';
+import { computeMarchingCubesLines } from '../../../util/marching-cubes/algorithm';
 import { Lines } from '../../../geometry/lines/lines';
-import { LinesBuilder } from '../../../geometry/lines/lines-builder';
 import { GaussianDensityProps, DefaultGaussianDensityProps } from 'mol-model/structure/structure/unit/gaussian-density';
 
 async function createGaussianWireframe(ctx: RuntimeContext, unit: Unit, structure: Structure, props: GaussianDensityProps, lines?: Lines): Promise<Lines> {
     const { smoothness } = props
     const { transform, field, idField } = await unit.computeGaussianDensity(props, ctx)
 
-    const surface = await computeMarchingCubes({
+    const params = {
         isoLevel: Math.exp(-smoothness),
         scalarField: field,
         idField
-    }).runAsChild(ctx)
-
-    Mesh.transformImmediate(surface, transform)
-
-    const vb = surface.vertexBuffer.ref.value
-    const ib = surface.indexBuffer.ref.value
-    const gb = surface.groupBuffer.ref.value
-
-    const builder = LinesBuilder.create(surface.triangleCount * 3, surface.triangleCount / 10, lines)
-
-    // TODO avoid duplicate lines and move to '../../../geometry/lines/lines' as Lines.fromMesh()
-    for (let i = 0, il = surface.triangleCount * 3; i < il; i += 3) {
-        const i0 = ib[i], i1 = ib[i + 1], i2 = ib[i + 2];
-        const x0 = vb[i0 * 3], y0 = vb[i0 * 3 + 1], z0 = vb[i0 * 3 + 2];
-        const x1 = vb[i1 * 3], y1 = vb[i1 * 3 + 1], z1 = vb[i1 * 3 + 2];
-        const x2 = vb[i2 * 3], y2 = vb[i2 * 3 + 1], z2 = vb[i2 * 3 + 2];
-        builder.add(x0, y0, z0, x1, y1, z1, gb[i0])
-        builder.add(x0, y0, z0, x2, y2, z2, gb[i0])
-        builder.add(x1, y1, z1, x2, y2, z2, gb[i1])
     }
+    const wireframe = await computeMarchingCubesLines(params, lines).runAsChild(ctx)
+
+    Lines.transformImmediate(wireframe, transform)
 
-    const l = builder.getLines();
-    console.log(l)
-    return l
+    return wireframe
 }
 
 export const DefaultGaussianWireframeProps = {

+ 2 - 2
src/mol-geo/representation/volume/surface.ts

@@ -7,7 +7,7 @@
 
 import { VolumeData, VolumeIsoValue } from 'mol-model/volume'
 import { Task, RuntimeContext } from 'mol-task'
-import { computeMarchingCubes } from '../../util/marching-cubes/algorithm';
+import { computeMarchingCubesMesh } from '../../util/marching-cubes/algorithm';
 import { Mesh } from '../../geometry/mesh/mesh';
 import { VolumeVisual } from '.';
 import { createMeshRenderObject, MeshRenderObject } from 'mol-gl/render-object';
@@ -23,7 +23,7 @@ export function computeVolumeSurface(volume: VolumeData, isoValue: VolumeIsoValu
     return Task.create<Mesh>('Volume Surface', async ctx => {
         ctx.update({ message: 'Marching cubes...' });
 
-        const mesh = await computeMarchingCubes({
+        const mesh = await computeMarchingCubesMesh({
             isoLevel: VolumeIsoValue.toAbsolute(isoValue).absoluteValue,
             scalarField: volume.data
         }).runAsChild(ctx);

+ 66 - 86
src/mol-geo/util/marching-cubes/algorithm.ts

@@ -2,41 +2,78 @@
  * Copyright (c) 2018 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>
  */
 
 import { Task, RuntimeContext } from 'mol-task'
-import { ChunkedArray } from 'mol-data/util'
 import { Tensor } from 'mol-math/linear-algebra'
 import { Mesh } from '../../geometry/mesh/mesh'
 import { Index, EdgeIdInfo, CubeEdges, EdgeTable, TriTable } from './tables'
-import { ValueCell } from 'mol-util'
+import { defaults } from 'mol-util'
+import { MarchinCubesBuilder, MarchinCubesMeshBuilder, MarchinCubesLinesBuilder } from './builder';
+import { Lines } from '../../geometry/lines/lines';
+// import { Lines } from '../../geometry/lines/lines';
 
 /**
  * The parameters required by the algorithm.
  */
-export interface MarchingCubesParameters {
+export interface MarchingCubesParams {
     isoLevel: number,
     scalarField: Tensor,
+    bottomLeft?: ReadonlyArray<number>,
+    topRight?: ReadonlyArray<number>,
+    idField?: Tensor,
+}
 
-    bottomLeft?: ArrayLike<number>,
-    topRight?: ArrayLike<number>,
+interface MarchingCubesInputParams extends MarchingCubesParams {
+    bottomLeft: ReadonlyArray<number>,
+    topRight: ReadonlyArray<number>,
+}
 
-    idField?: Tensor,
+function getInputParams(params: MarchingCubesParams): MarchingCubesInputParams {
+    return {
+        ...params,
+        bottomLeft: defaults(params.bottomLeft, [0, 0, 0] as ReadonlyArray<number>),
+        topRight: defaults(params.topRight, params.scalarField.space.dimensions)
+    }
+}
 
-    oldSurface?: Mesh
+function getExtent(inputParams: MarchingCubesInputParams) {
+    return {
+        dX: inputParams.topRight[0] - inputParams.bottomLeft[0],
+        dY: inputParams.topRight[1] - inputParams.bottomLeft[1],
+        dZ: inputParams.topRight[2] - inputParams.bottomLeft[2]
+    }
+}
+
+export function computeMarchingCubesMesh(params: MarchingCubesParams, mesh?: Mesh) {
+    return Task.create('Marching Cubes Mesh', async ctx => {
+        const inputParams = getInputParams(params)
+        const { dX, dY, dZ } = getExtent(inputParams)
+        // TODO should it be configurable? Scalar fields can produce meshes with vastly different densities.
+        const vertexChunkSize = Math.min(262144, Math.max(dX * dY * dZ / 32, 1024))
+        const builder = MarchinCubesMeshBuilder(vertexChunkSize, mesh)
+        await (new MarchingCubesComputation(ctx, builder, inputParams)).run()
+        return builder.get()
+    });
 }
 
-export function computeMarchingCubes(parameters: MarchingCubesParameters) {
-    return Task.create('Marching Cubes', async ctx => {
-        const comp = new MarchingCubesComputation(parameters, ctx);
-        return await comp.run();
+export function computeMarchingCubesLines(params: MarchingCubesParams, lines?: Lines) {
+    return Task.create('Marching Cubes Lines', async ctx => {
+        const inputParams = getInputParams(params)
+        const { dX, dY, dZ } = getExtent(inputParams)
+        // TODO should it be configurable? Scalar fields can produce meshes with vastly different densities.
+        const vertexChunkSize = Math.min(262144, Math.max(dX * dY * dZ / 32, 1024))
+        const builder = MarchinCubesLinesBuilder(vertexChunkSize, lines)
+        await (new MarchingCubesComputation(ctx, builder, inputParams)).run()
+        return builder.get()
     });
 }
 
 class MarchingCubesComputation {
     private size: number;
     private sliceSize: number;
-    private parameters: MarchingCubesParameters;
+    private edgeFilter: number
 
     private minX = 0; private minY = 0; private minZ = 0;
     private maxX = 0; private maxY = 0; private maxZ = 0;
@@ -45,7 +82,8 @@ class MarchingCubesComputation {
     private async doSlices() {
         let done = 0;
 
-        for (let k = this.minZ; k < this.maxZ; k++) {
+        this.edgeFilter = 15
+        for (let k = this.minZ; k < this.maxZ; k++, this.edgeFilter &= ~4) {
             this.slice(k);
 
             done += this.sliceSize;
@@ -56,52 +94,24 @@ class MarchingCubesComputation {
     }
 
     private slice(k: number) {
-        for (let j = this.minY; j < this.maxY; j++) {
-            for (let i = this.minX; i < this.maxX; i++) {
-                this.state.processCell(i, j, k);
+        this.edgeFilter |= 2
+        for (let j = this.minY; j < this.maxY; j++, this.edgeFilter &= ~2) {
+            this.edgeFilter |= 1
+            for (let i = this.minX; i < this.maxX; i++, this.edgeFilter &= ~1) {
+                this.state.processCell(i, j, k, this.edgeFilter);
             }
         }
         this.state.clearEdgeVertexIndexSlice(k);
     }
 
-    private finish(): Mesh {
-        const vb = ChunkedArray.compact(this.state.vertexBuffer, true) as Float32Array;
-        const ib = ChunkedArray.compact(this.state.triangleBuffer, true) as Uint32Array;
-        const gb = ChunkedArray.compact(this.state.groupBuffer, true) as Float32Array;
-
-        this.state.vertexBuffer = <any>void 0;
-        this.state.verticesOnEdges = <any>void 0;
-        this.state.groupBuffer = <any>void 0;
-
-        const os = this.parameters.oldSurface
-
-        return {
-            kind: 'mesh',
-            vertexCount:  this.state.vertexCount,
-            triangleCount: this.state.triangleCount,
-            vertexBuffer: os ? ValueCell.update(os.vertexBuffer, vb) : ValueCell.create(vb),
-            groupBuffer: os ? ValueCell.update(os.groupBuffer, gb) : ValueCell.create(gb),
-            indexBuffer: os ? ValueCell.update(os.indexBuffer, ib) : ValueCell.create(ib),
-            normalBuffer: os ? os.normalBuffer : ValueCell.create(new Float32Array(0)),
-            normalsComputed: false
-        }
-    }
-
     async run() {
         await this.ctx.update({ message: 'Computing surface...', current: 0, max: this.size });
         await this.doSlices();
         await this.ctx.update('Finalizing...');
-        return this.finish();
     }
 
-    constructor(parameters: MarchingCubesParameters, private ctx: RuntimeContext) {
-        const params = { ...parameters };
-        this.parameters = params;
-
-        if (!params.bottomLeft) params.bottomLeft = [0, 0, 0];
-        if (!params.topRight) params.topRight = params.scalarField.space.dimensions;
-
-        this.state = new MarchingCubesState(params);
+    constructor(private ctx: RuntimeContext, builder: MarchinCubesBuilder<any>, params: MarchingCubesInputParams) {
+        this.state = new MarchingCubesState(builder, params);
         this.minX = params.bottomLeft[0];
         this.minY = params.bottomLeft[1];
         this.minZ = params.bottomLeft[2];
@@ -121,19 +131,12 @@ class MarchingCubesState {
     scalarField: Tensor.Data;
     idFieldGet?: Tensor.Space['get'];
     idField?: Tensor.Data;
-    assignIds: boolean;
 
     // two layers of vertex indices. Each vertex has 3 edges associated.
     verticesOnEdges: Int32Array;
     vertList: number[] = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
     i: number = 0; j: number = 0; k: number = 0;
 
-    vertexBuffer: ChunkedArray<number, 3>;
-    groupBuffer: ChunkedArray<number, 1>;
-    triangleBuffer: ChunkedArray<number, 3>;
-    vertexCount = 0;
-    triangleCount = 0;
-
     private get3dOffsetFromEdgeInfo(index: Index) {
         return (this.nX * (((this.k + index.k) % 2) * this.nY + this.j + index.j) + this.i + index.i);
     }
@@ -153,7 +156,7 @@ class MarchingCubesState {
         const edgeId = 3 * this.get3dOffsetFromEdgeInfo(info) + info.e;
 
         const ret = this.verticesOnEdges[edgeId];
-        if (ret > 0) return (ret - 1) | 0;
+        if (ret > 0) return ret - 1;
 
         const edge = CubeEdges[edgeNum];
         const a = edge.a, b = edge.b;
@@ -163,13 +166,7 @@ class MarchingCubesState {
         const v1 = this.scalarFieldGet(this.scalarField, hi, hj, hk);
         const t = (this.isoLevel - v0) / (v0 - v1);
 
-        const id = ChunkedArray.add3(
-            this.vertexBuffer,
-            li + t * (li - hi),
-            lj + t * (lj - hj),
-            lk + t * (lk - hk)
-        ) | 0;
-
+        const id = this.builder.addVertex(li + t * (li - hi), lj + t * (lj - hj), lk + t * (lk - hk));
         this.verticesOnEdges[edgeId] = id + 1;
 
         if (this.idField) {
@@ -177,17 +174,15 @@ class MarchingCubesState {
             const v = this.idFieldGet!(this.idField, hi, hj, hk)
             let a = t < 0.5 ? u : v;
             if (a < 0) a = t < 0.5 ? v : u;
-            ChunkedArray.add(this.groupBuffer, a);
+            this.builder.addGroup(a);
         } else {
-            ChunkedArray.add(this.groupBuffer, 0);
+            this.builder.addGroup(0);
         }
 
-        this.vertexCount++;
-
         return id;
     }
 
-    constructor(params: MarchingCubesParameters) {
+    constructor(private builder: MarchinCubesBuilder<any>, params: MarchingCubesInputParams) {
         const dims = params.scalarField.space.dimensions;
         this.nX = dims[0]; this.nY = dims[1]; this.nZ = dims[2];
         this.isoLevel = params.isoLevel;
@@ -198,20 +193,6 @@ class MarchingCubesState {
             this.idFieldGet = params.idField.space.get;
         }
 
-        const dX = params.topRight![0] - params.bottomLeft![0]
-        const dY = params.topRight![1] - params.bottomLeft![1]
-        const dZ = params.topRight![2] - params.bottomLeft![2]
-        // TODO should it be configurable? Scalar fields can produce meshes with vastly different densities.
-        const vertexBufferSize = Math.min(262144, Math.max(dX * dY * dZ / 32, 1024) | 0)
-        const triangleBufferSize = Math.min(1 << 16, vertexBufferSize * 4)
-
-        this.vertexBuffer = ChunkedArray.create(Float32Array, 3, vertexBufferSize,
-            params.oldSurface && params.oldSurface.vertexBuffer.ref.value);
-        this.groupBuffer = ChunkedArray.create(Float32Array, 1, vertexBufferSize,
-            params.oldSurface && params.oldSurface.groupBuffer.ref.value);
-        this.triangleBuffer = ChunkedArray.create(Uint32Array, 3, triangleBufferSize,
-            params.oldSurface && params.oldSurface.indexBuffer.ref.value);
-
         // two layers of vertex indices. Each vertex has 3 edges associated.
         this.verticesOnEdges = new Int32Array(3 * this.nX * this.nY * 2);
     }
@@ -220,7 +201,7 @@ class MarchingCubesState {
         return this.scalarFieldGet(this.scalarField, i, j, k);
     }
 
-    processCell(i: number, j: number, k: number) {
+    processCell(i: number, j: number, k: number, edgeFilter: number) {
         let tableIndex = 0;
 
         if (this.get(i, j, k) < this.isoLevel) tableIndex |= 1;
@@ -235,7 +216,7 @@ class MarchingCubesState {
         if (tableIndex === 0 || tableIndex === 255) return;
 
         this.i = i; this.j = j; this.k = k;
-        let edgeInfo = EdgeTable[tableIndex];
+        const edgeInfo = EdgeTable[tableIndex];
         if ((edgeInfo & 1) > 0) this.vertList[0] = this.interpolate(0); // 0 1
         if ((edgeInfo & 2) > 0) this.vertList[1] = this.interpolate(1); // 1 2
         if ((edgeInfo & 4) > 0) this.vertList[2] = this.interpolate(2); // 2 3
@@ -251,8 +232,7 @@ class MarchingCubesState {
 
         const triInfo = TriTable[tableIndex];
         for (let t = 0; t < triInfo.length; t += 3) {
-            this.triangleCount++;
-            ChunkedArray.add3(this.triangleBuffer, this.vertList[triInfo[t]], this.vertList[triInfo[t + 1]], this.vertList[triInfo[t + 2]]);
+            this.builder.addTriangle(this.vertList, triInfo[t], triInfo[t + 1], triInfo[t + 2], edgeFilter)
         }
     }
 }

+ 111 - 0
src/mol-geo/util/marching-cubes/builder.ts

@@ -0,0 +1,111 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ * @author David Sehnal <david.sehnal@gmail.com>
+ * @author Fred Ludlow <fred.ludlow@gmail.com>
+ */
+
+import { ChunkedArray } from '../../../mol-data/util';
+import { ValueCell } from 'mol-util';
+import { Mesh } from '../../geometry/mesh/mesh';
+import { AllowedContours } from './tables';
+import { LinesBuilder } from '../../geometry/lines/lines-builder';
+import { Lines } from '../../geometry/lines/lines';
+
+ export interface MarchinCubesBuilder<T> {
+    addVertex(x: number, y: number, z: number): number
+    addGroup(group: number): void
+    addTriangle(vertList: number[], a: number, b: number, c: number, edgeFilter: number): void
+    get(): T
+}
+
+export function MarchinCubesMeshBuilder(vertexChunkSize: number, mesh?: Mesh): MarchinCubesBuilder<Mesh> {
+    const triangleChunkSize = Math.min(1 << 16, vertexChunkSize * 4)
+
+    const vertices = ChunkedArray.create(Float32Array, 3, vertexChunkSize, mesh && mesh.vertexBuffer.ref.value);
+    const groups = ChunkedArray.create(Float32Array, 1, vertexChunkSize, mesh && mesh.groupBuffer.ref.value);
+    const indices = ChunkedArray.create(Uint32Array, 3, triangleChunkSize, mesh && mesh.indexBuffer.ref.value);
+
+    let vertexCount = 0
+    let triangleCount = 0
+
+    return {
+        addVertex: (x: number, y: number, z: number) => {
+            ++vertexCount
+            return ChunkedArray.add3(vertices, x, y, z );
+        },
+        addGroup: (group: number) => {
+            ChunkedArray.add(groups, group);
+        },
+        addTriangle: (vertList: number[], a: number, b: number, c: number) => {
+            ++triangleCount
+            ChunkedArray.add3(indices, vertList[a], vertList[b], vertList[c]);
+        },
+        get: () => {
+            const vb = ChunkedArray.compact(vertices, true) as Float32Array;
+            const ib = ChunkedArray.compact(indices, true) as Uint32Array;
+            const gb = ChunkedArray.compact(groups, true) as Float32Array;
+
+            return {
+                kind: 'mesh',
+                vertexCount,
+                triangleCount,
+                vertexBuffer: mesh ? ValueCell.update(mesh.vertexBuffer, vb) : ValueCell.create(vb),
+                groupBuffer: mesh ? ValueCell.update(mesh.groupBuffer, gb) : ValueCell.create(gb),
+                indexBuffer: mesh ? ValueCell.update(mesh.indexBuffer, ib) : ValueCell.create(ib),
+                normalBuffer: mesh ? mesh.normalBuffer : ValueCell.create(new Float32Array(0)),
+                normalsComputed: false
+            }
+        }
+    }
+}
+
+export function MarchinCubesLinesBuilder(vertexChunkSize: number, lines?: Lines): MarchinCubesBuilder<Lines> {
+    const vertices = ChunkedArray.create(Float32Array, 3, vertexChunkSize);
+    const groups = ChunkedArray.create(Float32Array, 1, vertexChunkSize);
+    const indices = ChunkedArray.create(Float32Array, 2, vertexChunkSize);
+
+    let linesCount = 0
+
+    return {
+        addVertex: (x: number, y: number, z: number) => {
+            return ChunkedArray.add3(vertices, x, y, z);
+        },
+        addGroup: (group: number) => {
+            ChunkedArray.add(groups, group);
+        },
+        addTriangle: (vertList: number[], a: number, b: number, c: number, edgeFilter: number) => {
+            if (AllowedContours[a][b] & edgeFilter) {
+                ++linesCount
+                ChunkedArray.add2(indices, vertList[a], vertList[b])
+            }
+            if (AllowedContours[b][c] & edgeFilter) {
+                ++linesCount
+                ChunkedArray.add2(indices, vertList[b], vertList[c])
+            }
+            if (AllowedContours[a][c] & edgeFilter) {
+                ++linesCount
+                ChunkedArray.add2(indices, vertList[a], vertList[c])
+            }
+        },
+        get: () => {
+            const vb = ChunkedArray.compact(vertices, true) as Float32Array;
+            const ib = ChunkedArray.compact(indices, true) as Uint32Array;
+            const gb = ChunkedArray.compact(groups, true) as Float32Array;
+
+            const builder = LinesBuilder.create(linesCount, linesCount / 10, lines)
+
+            for (let i = 0; i < linesCount; ++i) {
+                const la = ib[i * 2], lb = ib[i * 2 + 1]
+                builder.add(
+                    vb[la * 3], vb[la * 3 + 1], vb[la * 3 + 2],
+                    vb[lb * 3], vb[lb * 3 + 1], vb[lb * 3 + 2],
+                    gb[la]
+                )
+            }
+
+            return builder.getLines()
+        }
+    }
+}

+ 35 - 0
src/mol-geo/util/marching-cubes/tables.ts

@@ -2,6 +2,7 @@
  * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author David Sehnal <david.sehnal@gmail.com>
+ * @author Fred Ludlow <fred.ludlow@gmail.com>
  */
 
 export interface Index { i: number, j: number, k: number }
@@ -417,4 +418,38 @@ export const TriTable = [
     [0, 9, 1],
     [0, 3, 8],
     []
+];
+
+/**
+ * Triangles are constructed between points on cube edges.
+ * AllowedContours[edge1][edge1] indicates which lines from a given
+ * triangle should be shown in line mode.
+ *
+ * Values are bitmasks:
+ * In loop over cubes we keep another bitmask indicating whether our current
+ * cell is the first x-value (1),
+ * first y-value (2) or first z-value (4) of the current loop.
+ * We draw all lines on leading faces but only draw trailing face lines the first
+ * time through the loop
+ * A value of 8 below means the edge is always drawn (leading face)
+ *
+ * E.g. the first row, lines between edge0 and other edges in the bottom
+ * x-y plane are only drawn for the first value of z, edges in the
+ * x-z plane are only drawn for the first value of y. No other lines
+ * are drawn as they're redundant
+ * The line between edge 1 and 5 is always drawn as it's on the leading edge
+ */
+export const AllowedContours = [
+    [ 0, 4, 4, 4, 2, 0, 0, 0, 2, 2, 0, 0 ], // 1 2 3 4 8 9
+    [ 4, 0, 4, 4, 0, 8, 0, 0, 0, 8, 8, 0 ], // 0 2 3 5 9 10
+    [ 4, 4, 0, 4, 0, 0, 8, 0, 0, 0, 8, 8 ], // 0 1 3 6 10 11
+    [ 4, 4, 4, 0, 0, 0, 0, 1, 1, 0, 0, 1 ], // 0 1 2 7 8 11
+    [ 2, 0, 0, 0, 0, 8, 8, 8, 2, 2, 0, 0 ], // 0 5 6 7 8 9
+    [ 0, 8, 0, 0, 8, 0, 8, 8, 0, 8, 8, 0 ], // And rotate it
+    [ 0, 0, 8, 0, 8, 8, 0, 8, 0, 0, 8, 8 ],
+    [ 0, 0, 0, 1, 8, 8, 8, 0, 1, 0, 0, 1 ],
+    [ 2, 0, 0, 1, 2, 0, 0, 1, 0, 2, 0, 1 ], // 0 3 4 7 9 11
+    [ 2, 8, 0, 0, 2, 8, 0, 0, 2, 0, 8, 0 ], // And rotate some more
+    [ 0, 8, 8, 0, 0, 8, 8, 0, 0, 8, 0, 8 ],
+    [ 0, 0, 8, 1, 0, 0, 8, 1, 1, 0, 8, 0 ]
 ];