Merge pull request #244 from molstar/meshproc

Mesh processing: border smoothing

CHANGELOG.md

@@ -6,6 +6,10 @@ Note that since we don't clearly distinguish between a public and private interf
 
 ## [Unreleased]
 
+- Add ``invertCantorPairing`` helper function
+- Add ``Mesh`` processing helper ``.smoothEdges``
+- Smooth border of molecular-surface with ``includeParent`` enabled
+- Hide ``includeParent`` option from gaussian-surface visuals (not particularly useful)
 
 ## [v2.2.2] - 2021-08-11
 

src/mol-data/util/hash-functions.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2017-2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2017-2021 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>
@@ -70,6 +70,15 @@ export function sortedCantorPairing(a: number, b: number) {
     return a < b ? cantorPairing(a, b) : cantorPairing(b, a);
 }
 
+export function invertCantorPairing(out: [number, number], z: number) {
+    const w = Math.floor((Math.sqrt(8 * z + 1) - 1) / 2);
+    const t = (w * w + w) / 2;
+    const y = z - t;
+    out[0] = w - y;
+    out[1] = y;
+    return out;
+}
+
 /**
  * 32 bit FNV-1a hash, see http://isthe.com/chongo/tech/comp/fnv/
  */

src/mol-geo/geometry/mesh/mesh.ts

@@ -8,13 +8,13 @@
 import { ValueCell } from '../../../mol-util';
 import { Vec3, Mat4, Mat3, Vec4 } from '../../../mol-math/linear-algebra';
 import { Sphere3D } from '../../../mol-math/geometry';
-import { transformPositionArray, transformDirectionArray, computeIndexedVertexNormals, GroupMapping, createGroupMapping} from '../../util';
+import { transformPositionArray, transformDirectionArray, computeIndexedVertexNormals, GroupMapping, createGroupMapping } from '../../util';
 import { GeometryUtils } from '../geometry';
 import { createMarkers } from '../marker-data';
 import { TransformData } from '../transform-data';
 import { LocationIterator, PositionLocation } from '../../util/location-iterator';
 import { createColors } from '../color-data';
-import { ChunkedArray, hashFnv32a } from '../../../mol-data/util';
+import { ChunkedArray, hashFnv32a, invertCantorPairing, sortedCantorPairing } from '../../../mol-data/util';
 import { ParamDefinition as PD } from '../../../mol-util/param-definition';
 import { calculateInvariantBoundingSphere, calculateTransformBoundingSphere } from '../../../mol-gl/renderable/util';
 import { Theme } from '../../../mol-theme/theme';
@@ -25,6 +25,8 @@ import { createEmptyOverpaint } from '../overpaint-data';
 import { createEmptyTransparency } from '../transparency-data';
 import { createEmptyClipping } from '../clipping-data';
 import { RenderableState } from '../../../mol-gl/renderable';
+import { arraySetAdd } from '../../../mol-util/array';
+import { degToRad } from '../../../mol-math/misc';
 
 export interface Mesh {
     readonly kind: 'mesh',
@@ -332,10 +334,10 @@ export namespace Mesh {
         mesh.vertexCount = newVertexCount;
         mesh.triangleCount = newTriangleCount;
 
-        ValueCell.update(vertexBuffer, newVb) as ValueCell<Float32Array>;
-        ValueCell.update(groupBuffer, newGb) as ValueCell<Float32Array>;
-        ValueCell.update(indexBuffer, newIb) as ValueCell<Uint32Array>;
-        ValueCell.update(normalBuffer, newNb) as ValueCell<Float32Array>;
+        ValueCell.update(vertexBuffer, newVb);
+        ValueCell.update(groupBuffer, newGb);
+        ValueCell.update(indexBuffer, newIb);
+        ValueCell.update(normalBuffer, newNb);
 
         // keep some original data, e.g., for geometry export
         (mesh.meta.originalData as OriginalData) = { indexBuffer: ib, vertexCount, triangleCount };
@@ -345,6 +347,276 @@ export namespace Mesh {
 
     //
 
+    function getNeighboursMap(mesh: Mesh) {
+        const { vertexCount, triangleCount } = mesh;
+        const elements = mesh.indexBuffer.ref.value;
+
+        const neighboursMap: number[][] = [];
+        for (let i = 0; i < vertexCount; ++i) {
+            neighboursMap[i] = [];
+        }
+
+        for (let i = 0; i < triangleCount; ++i) {
+            const v1 = elements[i * 3];
+            const v2 = elements[i * 3 + 1];
+            const v3 = elements[i * 3 + 2];
+            arraySetAdd(neighboursMap[v1], v2);
+            arraySetAdd(neighboursMap[v1], v3);
+            arraySetAdd(neighboursMap[v2], v1);
+            arraySetAdd(neighboursMap[v2], v3);
+            arraySetAdd(neighboursMap[v3], v1);
+            arraySetAdd(neighboursMap[v3], v2);
+        }
+        return neighboursMap;
+    }
+
+    function getEdgeCounts(mesh: Mesh) {
+        const { triangleCount } = mesh;
+        const elements = mesh.indexBuffer.ref.value;
+
+        const edgeCounts = new Map<number, number>();
+        const add = (a: number, b: number) => {
+            const z = sortedCantorPairing(a, b);
+            const c = edgeCounts.get(z) || 0;
+            edgeCounts.set(z, c + 1);
+        };
+
+        for (let i = 0; i < triangleCount; ++i) {
+            const a = elements[i * 3];
+            const b = elements[i * 3 + 1];
+            const c = elements[i * 3 + 2];
+            add(a, b); add(a, c); add(b, c);
+        }
+        return edgeCounts;
+    }
+
+    function getBorderVertices(edgeCounts: Map<number, number>) {
+        const borderVertices = new Set<number>();
+        const pair: [number, number] = [0, 0];
+        edgeCounts.forEach((c, z) => {
+            if (c === 1) {
+                invertCantorPairing(pair, z);
+                borderVertices.add(pair[0]);
+                borderVertices.add(pair[1]);
+            }
+        });
+
+        return borderVertices;
+    }
+
+    function getBorderNeighboursMap(neighboursMap: number[][], borderVertices: Set<number>, edgeCounts: Map<number, number>) {
+        const borderNeighboursMap = new Map<number, number[]>();
+        const add = (v: number, nb: number) => {
+            if (borderNeighboursMap.has(v)) arraySetAdd(borderNeighboursMap.get(v)!, nb);
+            else borderNeighboursMap.set(v, [nb]);
+        };
+
+        borderVertices.forEach(v => {
+            const neighbours = neighboursMap[v];
+            for (const nb of neighbours) {
+                if (borderVertices.has(nb) && edgeCounts.get(sortedCantorPairing(v, nb)) === 1) {
+                    add(v, nb);
+                }
+            }
+        });
+
+        return borderNeighboursMap;
+    }
+
+    function trimEdges(mesh: Mesh, neighboursMap: number[][]) {
+        const { indexBuffer, triangleCount } = mesh;
+        const ib = indexBuffer.ref.value;
+
+        // new
+        const index = ChunkedArray.create(Uint32Array, 3, 1024, triangleCount);
+
+        let newTriangleCount = 0;
+        for (let i = 0; i < triangleCount; ++i) {
+            const a = ib[i * 3];
+            const b = ib[i * 3 + 1];
+            const c = ib[i * 3 + 2];
+            if (neighboursMap[a].length === 2 ||
+                neighboursMap[b].length === 2 ||
+                neighboursMap[c].length === 2) continue;
+
+            ChunkedArray.add3(index, a, b, c);
+            newTriangleCount += 1;
+        }
+
+        const newIb = ChunkedArray.compact(index);
+        mesh.triangleCount = newTriangleCount;
+        ValueCell.update(indexBuffer, newIb);
+
+        return mesh;
+    }
+
+    function fillEdges(mesh: Mesh, neighboursMap: number[][], borderNeighboursMap: Map<number, number[]>, maxLengthSquared: number) {
+        const { vertexBuffer, indexBuffer, normalBuffer, triangleCount } = mesh;
+        const vb = vertexBuffer.ref.value;
+        const ib = indexBuffer.ref.value;
+        const nb = normalBuffer.ref.value;
+
+        // new
+        const index = ChunkedArray.create(Uint32Array, 3, 1024, triangleCount);
+
+        let newTriangleCount = 0;
+        for (let i = 0; i < triangleCount; ++i) {
+            ChunkedArray.add3(index, ib[i * 3], ib[i * 3 + 1], ib[i * 3 + 2]);
+            newTriangleCount += 1;
+        }
+
+        const vA = Vec3();
+        const vB = Vec3();
+        const vC = Vec3();
+        const vD = Vec3();
+        const vAB = Vec3();
+        const vAC = Vec3();
+        const vAD = Vec3();
+        const vABC = Vec3();
+
+        const vAN = Vec3();
+        const vN = Vec3();
+
+        const AngleThreshold = degToRad(120);
+        const added = new Set<number>();
+
+        const indices = Array.from(borderNeighboursMap.keys())
+            .filter(v => borderNeighboursMap.get(v)!.length < 2)
+            .map(v => {
+                const bnd = borderNeighboursMap.get(v)!;
+
+                Vec3.fromArray(vA, vb, v * 3);
+                Vec3.fromArray(vB, vb, bnd[0] * 3);
+                Vec3.fromArray(vC, vb, bnd[1] * 3);
+                Vec3.sub(vAB, vB, vA);
+                Vec3.sub(vAC, vC, vA);
+
+                return [v, Vec3.angle(vAB, vAC)];
+            });
+
+        // start with the smallest angle
+        indices.sort(([, a], [, b]) => a - b);
+
+        for (const [v, angle] of indices) {
+            if (added.has(v) || angle > AngleThreshold) continue;
+
+            const nbs = borderNeighboursMap.get(v)!;
+            if (neighboursMap[nbs[0]].includes(nbs[1]) &&
+                !borderNeighboursMap.get(nbs[0])?.includes(nbs[1])
+            ) continue;
+
+            Vec3.fromArray(vA, vb, v * 3);
+            Vec3.fromArray(vB, vb, nbs[0] * 3);
+            Vec3.fromArray(vC, vb, nbs[1] * 3);
+            Vec3.sub(vAB, vB, vA);
+            Vec3.sub(vAC, vC, vA);
+            Vec3.add(vABC, vAB, vAC);
+
+            if (Vec3.squaredDistance(vA, vB) >= maxLengthSquared) continue;
+
+            let add = false;
+            for (const nb of neighboursMap[v]) {
+                if (nbs.includes(nb)) continue;
+
+                Vec3.fromArray(vD, vb, nb * 3);
+                Vec3.sub(vAD, vD, vA);
+                if (Vec3.dot(vABC, vAD) < 0) {
+                    add = true;
+                    break;
+                }
+            }
+            if (!add) continue;
+
+            Vec3.fromArray(vAN, nb, v * 3);
+            Vec3.triangleNormal(vN, vA, vB, vC);
+            if (Vec3.dot(vN, vAN) > 0) {
+                ChunkedArray.add3(index, v, nbs[0], nbs[1]);
+            } else {
+                ChunkedArray.add3(index, nbs[1], nbs[0], v);
+            }
+            added.add(v); added.add(nbs[0]); added.add(nbs[1]);
+            newTriangleCount += 1;
+        }
+
+        const newIb = ChunkedArray.compact(index);
+        mesh.triangleCount = newTriangleCount;
+        ValueCell.update(indexBuffer, newIb);
+
+        return mesh;
+    }
+
+    function laplacianEdgeSmoothing(mesh: Mesh, borderNeighboursMap: Map<number, number[]>, options: { iterations: number, lambda: number }) {
+        const { iterations, lambda } = options;
+
+        const a = Vec3();
+        const b = Vec3();
+        const c = Vec3();
+        const t = Vec3();
+
+        const mu = -lambda;
+
+        let dst = new Float32Array(mesh.vertexBuffer.ref.value.length);
+
+        const step = (f: number) => {
+            const pos = mesh.vertexBuffer.ref.value;
+            dst.set(pos);
+
+            borderNeighboursMap.forEach((nbs, v) => {
+                if (nbs.length !== 2) return;
+
+                Vec3.fromArray(a, pos, v * 3);
+                Vec3.fromArray(b, pos, nbs[0] * 3);
+                Vec3.fromArray(c, pos, nbs[1] * 3);
+
+                const wab = 1 / Vec3.distance(a, b);
+                const wac = 1 / Vec3.distance(a, c);
+                Vec3.scale(b, b, wab);
+                Vec3.scale(c, c, wac);
+
+                Vec3.add(t, b, c);
+                Vec3.scale(t, t, 1 / (wab + wac));
+                Vec3.sub(t, t, a);
+
+                Vec3.scale(t, t, f);
+                Vec3.add(t, a, t);
+
+                Vec3.toArray(t, dst, v * 3);
+            });
+
+            const tmp = mesh.vertexBuffer.ref.value;
+            ValueCell.update(mesh.vertexBuffer, dst);
+            dst = tmp;
+        };
+
+        for (let k = 0; k < iterations; ++k) {
+            step(lambda);
+            step(mu);
+        }
+    }
+
+    export function smoothEdges(mesh: Mesh, options: { iterations: number, maxNewEdgeLength: number }) {
+        trimEdges(mesh, getNeighboursMap(mesh));
+
+        for (let k = 0; k < 10; ++k) {
+            const oldTriangleCount = mesh.triangleCount;
+            const edgeCounts = getEdgeCounts(mesh);
+            const neighboursMap = getNeighboursMap(mesh);
+            const borderVertices = getBorderVertices(edgeCounts);
+            const borderNeighboursMap = getBorderNeighboursMap(neighboursMap, borderVertices, edgeCounts);
+            fillEdges(mesh, neighboursMap, borderNeighboursMap, options.maxNewEdgeLength * options.maxNewEdgeLength);
+            if (mesh.triangleCount === oldTriangleCount) break;
+        }
+
+        const edgeCounts = getEdgeCounts(mesh);
+        const neighboursMap = getNeighboursMap(mesh);
+        const borderVertices = getBorderVertices(edgeCounts);
+        const borderNeighboursMap = getBorderNeighboursMap(neighboursMap, borderVertices, edgeCounts);
+        laplacianEdgeSmoothing(mesh, borderNeighboursMap, { iterations: options.iterations, lambda: 0.5 });
+        return mesh;
+    }
+
+    //
+
     export const Params = {
         ...BaseGeometry.Params,
         doubleSided: PD.Boolean(false, BaseGeometry.CustomQualityParamInfo),

src/mol-repr/structure/visual/gaussian-density-volume.ts

@@ -44,6 +44,7 @@ export const GaussianDensityVolumeParams = {
     ...ComplexDirectVolumeParams,
     ...GaussianDensityParams,
     ignoreHydrogens: PD.Boolean(false),
+    includeParent: PD.Boolean(false, { isHidden: true }),
 };
 export type GaussianDensityVolumeParams = typeof GaussianDensityVolumeParams
 
@@ -99,6 +100,7 @@ export const UnitsGaussianDensityVolumeParams = {
     ...UnitsDirectVolumeParams,
     ...GaussianDensityParams,
     ignoreHydrogens: PD.Boolean(false),
+    includeParent: PD.Boolean(false, { isHidden: true }),
 };
 export type UnitsGaussianDensityVolumeParams = typeof UnitsGaussianDensityVolumeParams
 

src/mol-repr/structure/visual/gaussian-surface-mesh.ts

@@ -30,6 +30,7 @@ const SharedParams = {
     ...ColorSmoothingParams,
     ignoreHydrogens: PD.Boolean(false),
     tryUseGpu: PD.Boolean(true),
+    includeParent: PD.Boolean(false, { isHidden: true }),
 };
 type SharedParams = typeof SharedParams
 

src/mol-repr/structure/visual/gaussian-surface-wireframe.ts

@@ -42,6 +42,7 @@ export const GaussianWireframeParams = {
     sizeFactor: PD.Numeric(3, { min: 0, max: 10, step: 0.1 }),
     lineSizeAttenuation: PD.Boolean(false),
     ignoreHydrogens: PD.Boolean(false),
+    includeParent: PD.Boolean(false, { isHidden: true }),
 };
 export type GaussianWireframeParams = typeof GaussianWireframeParams
 

src/mol-repr/structure/visual/molecular-surface-mesh.ts

@@ -11,7 +11,7 @@ import { VisualContext } from '../../visual';
 import { Unit, Structure } from '../../../mol-model/structure';
 import { Theme } from '../../../mol-theme/theme';
 import { Mesh } from '../../../mol-geo/geometry/mesh/mesh';
-import { computeUnitMolecularSurface, MolecularSurfaceProps } from './util/molecular-surface';
+import { computeUnitMolecularSurface } from './util/molecular-surface';
 import { computeMarchingCubesMesh } from '../../../mol-geo/util/marching-cubes/algorithm';
 import { ElementIterator, getElementLoci, eachElement } from './util/element';
 import { VisualUpdateState } from '../../util';
@@ -29,6 +29,7 @@ export const MolecularSurfaceMeshParams = {
     ...ColorSmoothingParams,
 };
 export type MolecularSurfaceMeshParams = typeof MolecularSurfaceMeshParams
+export type MolecularSurfaceMeshProps = PD.Values<MolecularSurfaceMeshParams>
 
 type MolecularSurfaceMeta = {
     resolution?: number
@@ -37,7 +38,7 @@ type MolecularSurfaceMeta = {
 
 //
 
-async function createMolecularSurfaceMesh(ctx: VisualContext, unit: Unit, structure: Structure, theme: Theme, props: MolecularSurfaceProps, mesh?: Mesh): Promise<Mesh> {
+async function createMolecularSurfaceMesh(ctx: VisualContext, unit: Unit, structure: Structure, theme: Theme, props: MolecularSurfaceMeshProps, mesh?: Mesh): Promise<Mesh> {
     const { transform, field, idField, resolution } = await computeUnitMolecularSurface(structure, unit, props).runInContext(ctx.runtime);
 
     const params = {
@@ -47,6 +48,11 @@ async function createMolecularSurfaceMesh(ctx: VisualContext, unit: Unit, struct
     };
     const surface = await computeMarchingCubesMesh(params, mesh).runAsChild(ctx.runtime);
 
+    if (props.includeParent) {
+        const iterations = Math.ceil(2 / props.resolution);
+        Mesh.smoothEdges(surface, { iterations, maxNewEdgeLength: Math.sqrt(2) });
+    }
+
     Mesh.transform(surface, transform);
     if (ctx.webgl && !ctx.webgl.isWebGL2) Mesh.uniformTriangleGroup(surface);
 
@@ -71,6 +77,7 @@ export function MolecularSurfaceMeshVisual(materialId: number): UnitsVisual<Mole
             if (newProps.ignoreHydrogens !== currentProps.ignoreHydrogens) state.createGeometry = true;
             if (newProps.traceOnly !== currentProps.traceOnly) state.createGeometry = true;
             if (newProps.includeParent !== currentProps.includeParent) state.createGeometry = true;
+
             if (newProps.smoothColors.name !== currentProps.smoothColors.name) {
                 state.updateColor = true;
             } else if (newProps.smoothColors.name === 'on' && currentProps.smoothColors.name === 'on') {