Fix surface normals

src/apps/render-test/mcubes.ts

@@ -7,14 +7,14 @@
 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'
+import { Tensor, Mat4, Vec3 } 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;
+
+    const dx = (max[0] - min[0]) / (ii - 1);
+    const dy = (max[1] - min[1]) / (jj - 1);
+    const dz = (max[2] - min[2]) / (kk - 1);
 
     for (let i = 0, x = min[0]; i < ii; i++, x += dx) {
         for (let j = 0, y = min[1]; j < jj; j++, y += dy) {
@@ -35,15 +35,23 @@ export default async function computeSurface(f: (x: number, y: number, z: number
         field = Tensor.create(space, space.create(Float32Array));
     }
 
-    fillField(field, f, [-1.1, -1.1, -1.1], [1.1, 1.1, 1.1]);
+    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({
         scalarField: field,
         isoLevel: 0,
-        buffers: data ? {
-            vertex: data.surface.vertexBuffer,
-            index: data.surface.indexBuffer,
-            normal: data.surface.normalBuffer
-        } : void 0
+        oldSurface: data ? data.surface : void 0
     }));
+
+    const translation = Mat4.fromTranslation(Mat4.zero(), min);
+    const grid = Vec3.zero();
+    Vec3.fromArray(grid, field.space.dimensions as any, 0);
+    const size = Vec3.sub(Vec3.zero(), max, min);
+    const scale = Mat4.fromScaling(Mat4.zero(), Vec3.create(size[0] / (grid[0] - 1), size[1] / (grid[1] - 1), size[2] / (grid[2] - 1)));
+
+    const transform = Mat4.mul(Mat4.zero(), translation, scale);
+    Surface.transformImmediate(surface, transform);
+    Surface.computeNormalsImmediate(surface);
     return { surface, field };
 }

src/apps/render-test/state.ts

@@ -110,7 +110,7 @@ export default class State {
             ...createTransformAttributes(regl, transformArray1)
         })
 
-        let rr = 1.0;
+        let rr = 1;
         function cubesF(x: number, y: number, z: number) {
             return x * x + y * y + z * z - rr * rr;
             // const a = ca;
@@ -122,8 +122,8 @@ export default class State {
 
         const makeCubesMesh = () => MeshRenderable.create(regl,
             {
-                position: Attribute.create(regl, cubes.surface.vertexBuffer, { size: 3 }),
-                normal: Attribute.create(regl, cubes.surface.vertexBuffer, { size: 3 }),
+                position: Attribute.create(regl, cubes.surface.vertexBuffer.value, { size: 3 }),
+                normal: Attribute.create(regl, cubes.surface.normalBuffer.value!, { size: 3 }),
                 ...createTransformAttributes(regl, transformArray1)
             },
             {
@@ -135,18 +135,18 @@ export default class State {
                 'light.falloff': 0,
                 'light.radius': 500
             },
-            cubes.surface.indexBuffer
+            cubes.surface.indexBuffer.value
         );
 
         let mesh2 = makeCubesMesh();
 
-        const makeCubes = async () => {
-            rr = Math.random();
-            cubes = await mcubes(cubesF, cubes);
-            mesh2 = makeCubesMesh();
-            setTimeout(makeCubes, 1000 / 15);
-        };
-        makeCubes();
+        // const makeCubes = async () => {
+        //     rr = Math.random();
+        //     cubes = await mcubes(cubesF, cubes);
+        //     mesh2 = makeCubesMesh();
+        //     setTimeout(makeCubes, 1000 / 15);
+        // };
+        // makeCubes();
 
         // const mesh2 = MeshRenderable.create(regl,
         //     {

src/mol-geo/shape/surface.ts

@@ -5,16 +5,18 @@
  */
 
 import { Task } from 'mol-task'
+import { ValueBox } from 'mol-util'
+import { Vec3, Mat4 } from 'mol-math/linear-algebra'
 
 export interface Surface {
     vertexCount: number,
     triangleCount: number,
-    vertexBuffer: Float32Array,
-    indexBuffer: Uint32Array,
-    normalBuffer?: Float32Array,
+    vertexBuffer: ValueBox<Float32Array>,
+    indexBuffer: ValueBox<Uint32Array>,
+    normalBuffer: ValueBox<Float32Array | undefined>,
     normalsComputed: boolean,
 
-    vertexAnnotation?: ArrayLike<number>[]
+    vertexAnnotation?: ValueBox<ArrayLike<number>>
     //boundingSphere?: { center: Geometry.LinearAlgebra.Vector3, radius: number };
 }
 
@@ -22,22 +24,25 @@ 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 normals = surface.normalBuffer.value && surface.normalBuffer.value!.length >= surface.vertexCount * 3
+            ? surface.normalBuffer.value : new Float32Array(surface.vertexBuffer.value!.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]);
+        const v = surface.vertexBuffer.value, triangles = surface.indexBuffer.value;
 
-            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;
+        const x = Vec3.zero(), y = Vec3.zero(), z = Vec3.zero(), d1 = Vec3.zero(), d2 = Vec3.zero(), n = Vec3.zero();
+        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];
+
+            Vec3.fromArray(x, v, a);
+            Vec3.fromArray(y, v, b);
+            Vec3.fromArray(z, v, c);
+            Vec3.sub(d1, z, y);
+            Vec3.sub(d2, y, x);
+            Vec3.cross(n, d1, d2);
+
+            normals[a] += n[0]; normals[a + 1] += n[1]; normals[a + 2] += n[2];
+            normals[b] += n[0]; normals[b + 1] += n[1]; normals[b + 2] += n[2];
+            normals[c] += n[0]; normals[c + 1] += n[1]; normals[c + 2] += n[2];
         }
 
         for (let i = 0, ii = 3 * surface.vertexCount; i < ii; i += 3) {
@@ -46,8 +51,10 @@ export namespace Surface {
             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;
+
+           // console.log([normals[i], normals[i + 1], normals[i + 2]], [v[i], v[i + 1], v[i + 2]])
         }
-        surface.normalBuffer = normals;
+        surface.normalBuffer = ValueBox(surface.normalBuffer, normals);
         surface.normalsComputed = true;
     }
 
@@ -60,6 +67,22 @@ export namespace Surface {
             return surface;
         });
     }
+
+    export function transformImmediate(surface: Surface, t: Mat4) {
+        const p = Vec3.zero();
+        const vertices = surface.vertexBuffer.value;
+        for (let i = 0, _c = surface.vertexCount * 3; i < _c; i += 3) {
+            p[0] = vertices[i];
+            p[1] = vertices[i + 1];
+            p[2] = vertices[i + 2];
+            Vec3.transformMat4(p, p, t);
+            vertices[i] = p[0];
+            vertices[i + 1] = p[1];
+            vertices[i + 2] = p[2];
+        }
+        surface.normalsComputed = false;
+        //surface.boundingSphere = void 0;
+    }
 }
 
 //     function addVertex(src: Float32Array, i: number, dst: Float32Array, j: number) {
@@ -177,22 +200,6 @@ export namespace Surface {
 //         });
 //     }
 
-//     export function transformImmediate(surface: Surface, t: number[]) {
-//         const p = LinearAlgebra.Vector3.zero();
-//         const m = LinearAlgebra.Vector3.transformMat4;
-//         const vertices = surface.vertices;
-//         for (let i = 0, _c = surface.vertices.length; i < _c; i += 3) {
-//             p[0] = vertices[i];
-//             p[1] = vertices[i + 1];
-//             p[2] = vertices[i + 2];
-//             m(p, p, t);
-//             vertices[i] = p[0];
-//             vertices[i + 1] = p[1];
-//             vertices[i + 2] = p[2];
-//         }
-//         surface.normals = void 0;
-//         surface.boundingSphere = void 0;
-//     }
 
 //     export function transform(surface: Surface, t: number[]): Computation<Surface> {
 //         return computation<Surface>(async ctx => {

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

@@ -9,6 +9,7 @@ import { ChunkedArray } from 'mol-data/util'
 import { Tensor } from 'mol-math/linear-algebra'
 import { Surface } from '../../shape/surface'
 import { Index, EdgeIdInfo, CubeEdges, EdgeTable, TriTable } from './tables'
+import { ValueBox } from 'mol-util'
 
 /**
  * The parameters required by the algorithm.
@@ -22,12 +23,7 @@ export interface MarchingCubesParameters {
 
     annotationField?: Tensor,
 
-    buffers?: {
-        vertex?: Float32Array,
-        index?: Uint32Array,
-        normal?: Float32Array,
-        annotation?: ArrayLike<number>
-    }
+    oldSurface?: Surface
 }
 
 export function compute(parameters: MarchingCubesParameters) {
@@ -40,6 +36,7 @@ export function compute(parameters: MarchingCubesParameters) {
 class MarchingCubesComputation {
     private size: number;
     private sliceSize: number;
+    private parameters: MarchingCubesParameters;
 
     private minX = 0; private minY = 0; private minZ = 0;
     private maxX = 0; private maxY = 0; private maxZ = 0;
@@ -68,18 +65,23 @@ class MarchingCubesComputation {
     }
 
     private finish() {
-        const vertexBuffer = ChunkedArray.compact(this.state.vertexBuffer, true) as Float32Array;
-        const indexBuffer = ChunkedArray.compact(this.state.triangleBuffer, true) as Uint32Array;
+        const vb = ChunkedArray.compact(this.state.vertexBuffer, true) as Float32Array;
+        const ib = ChunkedArray.compact(this.state.triangleBuffer, true) as Uint32Array;
 
         this.state.vertexBuffer = <any>void 0;
         this.state.verticesOnEdges = <any>void 0;
 
         let ret: Surface = {
-            vertexCount: this.state.vertexCount,
+            vertexCount:  this.state.vertexCount,
             triangleCount: this.state.triangleCount,
-            vertexBuffer,
-            indexBuffer,
-            // vertexAnnotation: this.state.annotate ? ChunkedArray.compact(this.state.annotationBuffer) : void 0,
+            vertexBuffer: this.parameters.oldSurface ? ValueBox(this.parameters.oldSurface.vertexBuffer, vb) : ValueBox(vb),
+            indexBuffer: this.parameters.oldSurface ? ValueBox(this.parameters.oldSurface.indexBuffer, ib) : ValueBox(ib),
+            normalBuffer: this.parameters.oldSurface ? this.parameters.oldSurface.normalBuffer : ValueBox(void 0),
+            vertexAnnotation: this.state.annotate
+                ? this.parameters.oldSurface && this.parameters.oldSurface.vertexAnnotation
+                    ? ValueBox(this.parameters.oldSurface.vertexAnnotation, ChunkedArray.compact(this.state.annotationBuffer))
+                    : ValueBox(ChunkedArray.compact(this.state.annotationBuffer))
+                : void 0,
             normalsComputed: false
         }
 
@@ -98,6 +100,7 @@ class MarchingCubesComputation {
         private ctx: RuntimeContext) {
 
         let params = { ...parameters };
+        this.parameters = params;
 
         if (!params.bottomLeft) params.bottomLeft = [0, 0, 0];
         if (!params.topRight) params.topRight = params.scalarField.space.dimensions;
@@ -195,8 +198,10 @@ 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, params.buffers && params.buffers.vertex);
-        this.triangleBuffer = ChunkedArray.create<number>(s => new Uint32Array(s), 3, triangleBufferSize, params.buffers && params.buffers.index);
+        this.vertexBuffer = ChunkedArray.create<number>(s => new Float32Array(s), 3, vertexBufferSize,
+            params.oldSurface && params.oldSurface.vertexBuffer.value);
+        this.triangleBuffer = ChunkedArray.create<number>(s => new Uint32Array(s), 3, triangleBufferSize,
+            params.oldSurface && params.oldSurface.indexBuffer.value);
 
         this.annotate = !!params.annotationField;
         if (this.annotate) this.annotationBuffer = ChunkedArray.create(s => new Int32Array(s), 1, vertexBufferSize);