marching cubes

src/apps/render-test/mcubes.ts

@@ -0,0 +1,49 @@
+/**
+ * Copyright (c) 2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author David Sehnal <david.sehnal@gmail.com>
+ */
+
+import { Run } from 'mol-task'
+import { compute } from 'mol-geo/util/marching-cubes/algorithm'
+import { Surface } from 'mol-geo/shape/surface'
+import { Tensor } from 'mol-math/linear-algebra'
+
+function fillField(tensor: Tensor, f: (x: number, y: number, z: number) => number, min: number[], max: number[]): Tensor {
+    const { space: { set, dimensions: [ii, jj, kk] }, data } = tensor;
+ 
+    const dx = (max[0] - min[0]) / ii;
+    const dy = (max[1] - min[1]) / jj;
+    const dz = (max[2] - min[2]) / kk;
+
+    for (let i = 0, x = min[0]; i < ii; i++, x += dx) {
+        for (let j = 0, y = min[1]; j < jj; j++, y += dy) {
+            for (let k = 0, z = min[2]; k < kk; k++, z += dz) {
+                set(data, i, j, k, f(x, y, z));
+            }
+        }
+    }
+
+    return tensor
+}
+
+export default async function computeSurface(f: (x: number, y: number, z: number) => number, data?: { field: Tensor, surface: Surface }) {
+    let field: Tensor;
+    if (data) field = data.field;
+    else {
+        const space = Tensor.Space([30, 30, 30], [0, 1, 2]);
+        field = Tensor.create(space, space.create(Float32Array));
+    }
+
+    fillField(field, f, [-1.1, -1.1, -1.1], [1.1, 1.1, 1.1]);
+    const surface = await Run(compute({
+        scalarField: field,
+        isoLevel: 0,
+        buffers: data ? {
+            vertex: data.surface.vertexBuffer,
+            index: data.surface.indexBuffer,
+            normal: data.surface.normalBuffer
+        } : void 0
+    }));
+    return { surface, field };
+}

src/apps/render-test/state.ts

@@ -16,10 +16,12 @@ import { calculateTextureInfo } from 'mol-gl/util';
 import Icosahedron from 'mol-geo/primitive/icosahedron'
 import Box from 'mol-geo/primitive/box'
 
+import mcubes from './mcubes'
+
 export default class State {
     regl: REGL.Regl
 
-    initRegl (container: HTMLDivElement) {
+    async initRegl (container: HTMLDivElement) {
         const regl = glContext.create({
             container,
             extensions: [
@@ -108,11 +110,21 @@ export default class State {
             ...createTransformAttributes(regl, transformArray1)
         })
 
-        const mesh2 = MeshRenderable.create(regl,
+        let rr = 1.0;
+        function cubesF(x: number, y: number, z: number) {
+            return x * x + y * y + z * z - rr * rr;
+            // const a = ca;
+            // const t = (x + y + z + a);
+            // return x * x * x + y * y * y + z * z * z + a * a * a - t * t * t;
+        }
+
+        let cubes = await mcubes(cubesF);
+
+        const makeCubesMesh = () => MeshRenderable.create(regl,
             {
-                position: Attribute.create(regl, new Float32Array(box.vertices), { size: 3 }),
-                normal: Attribute.create(regl, new Float32Array(box.normals), { size: 3 }),
-                ...createTransformAttributes(regl, transformArray2)
+                position: Attribute.create(regl, cubes.surface.vertexBuffer, { size: 3 }),
+                normal: Attribute.create(regl, cubes.surface.vertexBuffer, { size: 3 }),
+                ...createTransformAttributes(regl, transformArray1)
             },
             {
                 colorTex,
@@ -123,8 +135,36 @@ export default class State {
                 'light.falloff': 0,
                 'light.radius': 500
             },
-            box.indices
-        )
+            cubes.surface.indexBuffer
+        );
+
+        let mesh2 = makeCubesMesh();
+
+        const makeCubes = async () => {
+            rr = Math.random();
+            cubes = await mcubes(cubesF, cubes);
+            mesh2 = makeCubesMesh();
+            setTimeout(makeCubes, 1000 / 15);
+        };
+        makeCubes();
+
+        // const mesh2 = MeshRenderable.create(regl,
+        //     {
+        //         position: Attribute.create(regl, new Float32Array(box.vertices), { size: 3 }),
+        //         normal: Attribute.create(regl, new Float32Array(box.normals), { size: 3 }),
+        //         ...createTransformAttributes(regl, transformArray1)
+        //     },
+        //     {
+        //         colorTex,
+        //         colorTexSize: [ colorTexInfo.width, colorTexInfo.height ],
+        //         'light.position': Vec3.create(0, 0, -20),
+        //         'light.color': Vec3.create(1.0, 1.0, 1.0),
+        //         'light.ambient': Vec3.create(0.5, 0.5, 0.5),
+        //         'light.falloff': 0,
+        //         'light.radius': 500
+        //     },
+        //     box.indices
+        // )
 
         const baseContext = regl({
             context: {
@@ -143,12 +183,12 @@ export default class State {
                 baseContext(() => {
                     // console.log(ctx)
                     regl.clear({color: [0, 0, 0, 1]})
-                    position.update(array => { array[0] = Math.random() })
+                    //position.update(array => { array[0] = Math.random() })
                     // points.update(a => { a.position[0] = Math.random() })
                     // mesh.draw()
                     // points.draw()
                     mesh2.draw()
-                    points2.draw()
+                    // points2.draw()
                     // model1({}, ({ transform }) => {
                     //     points.draw()
                     // })

src/mol-data/util/chunked-array.ts

@@ -75,6 +75,10 @@ namespace ChunkedArray {
     }
 
     export function compact<T>(array: ChunkedArray<T>, doNotResizeSingleton = false): ArrayLike<T> {
+        return _compact(array, doNotResizeSingleton);
+    }
+
+    export function _compact<T>(array: ChunkedArray<T>, doNotResizeSingleton: boolean): ArrayLike<T> {
         const { ctor, chunks, currentIndex } = array;
 
         if (!chunks.length) return ctor(0);

src/mol-geo/shape/surface.ts

@@ -4,6 +4,8 @@
  * @author David Sehnal <david.sehnal@gmail.com>
  */
 
+import { Task } from 'mol-task'
+
 export interface Surface {
     vertexCount: number,
     triangleCount: number,
@@ -16,47 +18,49 @@ export interface Surface {
     //boundingSphere?: { center: Geometry.LinearAlgebra.Vector3, radius: number };
 }
 
-// export namespace Surface {
-//     export function computeNormalsImmediate(surface: Surface) {
-//         if (surface.normals) return;
-
-//         const normals = new Float32Array(surface.vertices.length),
-//             v = surface.vertices, triangles = surface.triangleIndices;
-//         for (let i = 0; i < triangles.length; i += 3) {
-//             const a = 3 * triangles[i],
-//                 b = 3 * triangles[i + 1],
-//                 c = 3 * triangles[i + 2];
-
-//             const nx = v[a + 2] * (v[b + 1] - v[c + 1]) + v[b + 2] * v[c + 1] - v[b + 1] * v[c + 2] + v[a + 1] * (-v[b + 2] + v[c + 2]),
-//                 ny = -(v[b + 2] * v[c]) + v[a + 2] * (-v[b] + v[c]) + v[a] * (v[b + 2] - v[c + 2]) + v[b] * v[c + 2],
-//                 nz = v[a + 1] * (v[b] - v[c]) + v[b + 1] * v[c] - v[b] * v[c + 1] + v[a] * (-v[b + 1] + v[b + 1]);
-
-//             normals[a] += nx; normals[a + 1] += ny; normals[a + 2] += nz;
-//             normals[b] += nx; normals[b + 1] += ny; normals[b + 2] += nz;
-//             normals[c] += nx; normals[c + 1] += ny; normals[c + 2] += nz;
-//         }
-
-//         for (let i = 0; i < normals.length; i += 3) {
-//             const nx = normals[i];
-//             const ny = normals[i + 1];
-//             const nz = normals[i + 2];
-//             const f = 1.0 / Math.sqrt(nx * nx + ny * ny + nz * nz);
-//             normals[i] *= f; normals[i + 1] *= f; normals[i + 2] *= f;
-//         }
-//         surface.normals = normals;
-//     }
-
-//     export function computeNormals(surface: Surface): Computation<Surface> {
-//         return computation<Surface>(async ctx => {
-//             if (surface.normals) {
-//                 return surface;
-//             };
-
-//             await ctx.updateProgress('Computing normals...');
-//             computeNormalsImmediate(surface);
-//             return surface;
-//         });
-//     }
+export namespace Surface {
+    export function computeNormalsImmediate(surface: Surface) {
+        if (surface.normalsComputed) return;
+
+        const normals = surface.normalBuffer && surface.normalBuffer.length >= surface.vertexCount * 3
+            ? surface.normalBuffer : new Float32Array(surface.vertexBuffer.length);
+
+        const v = surface.vertexBuffer, triangles = surface.indexBuffer;
+        for (let i = 0, ii = 3 * surface.triangleCount; i < ii; i += 3) {
+            const a = 3 * triangles[i],
+                b = 3 * triangles[i + 1],
+                c = 3 * triangles[i + 2];
+
+            const nx = v[a + 2] * (v[b + 1] - v[c + 1]) + v[b + 2] * v[c + 1] - v[b + 1] * v[c + 2] + v[a + 1] * (-v[b + 2] + v[c + 2]),
+                ny = -(v[b + 2] * v[c]) + v[a + 2] * (-v[b] + v[c]) + v[a] * (v[b + 2] - v[c + 2]) + v[b] * v[c + 2],
+                nz = v[a + 1] * (v[b] - v[c]) + v[b + 1] * v[c] - v[b] * v[c + 1] + v[a] * (-v[b + 1] + v[b + 1]);
+
+            normals[a] += nx; normals[a + 1] += ny; normals[a + 2] += nz;
+            normals[b] += nx; normals[b + 1] += ny; normals[b + 2] += nz;
+            normals[c] += nx; normals[c + 1] += ny; normals[c + 2] += nz;
+        }
+
+        for (let i = 0, ii = 3 * surface.vertexCount; i < ii; i += 3) {
+            const nx = normals[i];
+            const ny = normals[i + 1];
+            const nz = normals[i + 2];
+            const f = 1.0 / Math.sqrt(nx * nx + ny * ny + nz * nz);
+            normals[i] *= f; normals[i + 1] *= f; normals[i + 2] *= f;
+        }
+        surface.normalBuffer = normals;
+        surface.normalsComputed = true;
+    }
+
+    export function computeNormals(surface: Surface): Task<Surface> {
+        return Task.create<Surface>('Surface (Compute Normals)', async ctx => {
+            if (surface.normalsComputed) return surface;
+
+            await ctx.update('Computing normals...');
+            computeNormalsImmediate(surface);
+            return surface;
+        });
+    }
+}
 
 //     function addVertex(src: Float32Array, i: number, dst: Float32Array, j: number) {
 //         dst[3 * j] += src[3 * i];

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

@@ -68,15 +68,15 @@ class MarchingCubesComputation {
     }
 
     private finish() {
-        const vertexBuffer = ChunkedArray.compact(this.state.vertexBuffer) as Float32Array;
-        const indexBuffer = ChunkedArray.compact(this.state.triangleBuffer) as Uint32Array;
+        const vertexBuffer = ChunkedArray.compact(this.state.vertexBuffer, true) as Float32Array;
+        const indexBuffer = ChunkedArray.compact(this.state.triangleBuffer, true) as Uint32Array;
 
         this.state.vertexBuffer = <any>void 0;
         this.state.verticesOnEdges = <any>void 0;
 
         let ret: Surface = {
-            triangleCount: 0,
-            vertexCount: 0,
+            vertexCount: this.state.vertexCount,
+            triangleCount: this.state.triangleCount,
             vertexBuffer,
             indexBuffer,
             // vertexAnnotation: this.state.annotate ? ChunkedArray.compact(this.state.annotationBuffer) : void 0,
@@ -128,6 +128,8 @@ class MarchingCubesState {
     vertexBuffer: ChunkedArray<number>;
     annotationBuffer: ChunkedArray<number>;
     triangleBuffer: ChunkedArray<number>;
+    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);
@@ -173,6 +175,8 @@ class MarchingCubesState {
             ChunkedArray.add(this.annotationBuffer, a);
         }
 
+        this.vertexCount++;
+
         return id;
     }
 
@@ -191,8 +195,8 @@ class MarchingCubesState {
             vertexBufferSize = Math.min(262144, Math.max(dX * dY * dZ / 16, 1024) | 0),
             triangleBufferSize = Math.min(1 << 16, vertexBufferSize * 4);
 
-        this.vertexBuffer = ChunkedArray.create<number>(s => new Float32Array(s), 3, vertexBufferSize);
-        this.triangleBuffer = ChunkedArray.create<number>(s => new Uint32Array(s), 3, triangleBufferSize);
+        this.vertexBuffer = ChunkedArray.create<number>(s => new Float32Array(s), 3, vertexBufferSize, params.buffers && params.buffers.vertex);
+        this.triangleBuffer = ChunkedArray.create<number>(s => new Uint32Array(s), 3, triangleBufferSize, params.buffers && params.buffers.index);
 
         this.annotate = !!params.annotationField;
         if (this.annotate) this.annotationBuffer = ChunkedArray.create(s => new Int32Array(s), 1, vertexBufferSize);
@@ -203,7 +207,7 @@ class MarchingCubesState {
 
     get(i: number, j: number, k: number) {
         return this.scalarFieldGet(this.scalarField, i, j, k);
-    } 
+    }
 
     processCell(i: number, j: number, k: number) {
         let tableIndex = 0;
@@ -236,6 +240,7 @@ class MarchingCubesState {
 
         let 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]]);
         }
     }

src/mol-math/linear-algebra/tensor.ts

@@ -41,6 +41,8 @@ export namespace Tensor {
         return { dimensions, axisOrderFastToSlow, axisOrderSlowToFast, accessDimensions, defaultCtor: ctor || Float64Array }
     }
 
+    export function create(space: Space, data: Data): Tensor { return { space, data }; }
+
     export function Space(dimensions: number[], axisOrderSlowToFast: number[], ctor?: ArrayCtor): Space {
         const layout = Layout(dimensions, axisOrderSlowToFast, ctor);
         const { get, set } = accessors(layout);