add oriented-box visual to structure-orientation

src/mol-math/linear-algebra/matrix/principal-axes.ts

@@ -5,8 +5,9 @@
  */
 
 import Matrix from './matrix';
-import { Vec3 } from '../3d';
+import { Vec3, Mat4 } from '../3d';
 import { svd } from './svd';
+import { NumberArray } from '../../../mol-util/type-helpers';
 
 export { PrincipalAxes }
 
@@ -78,95 +79,68 @@ namespace PrincipalAxes {
         }
     }
 
-    // TODO
-    // /**
-    //  * Get the basis matrix descriping the axes
-    //  * @param  {Matrix4} [optionalTarget] - target object
-    //  * @return {Matrix4} the basis
-    //  */
-    // getBasisMatrix(optionalTarget = new Matrix4()) {
-    //     const basis = optionalTarget
-
-    //     basis.makeBasis(this.normVecB, this.normVecA, this.normVecC)
-    //     if (basis.determinant() < 0) {
-    //         basis.scale(negateVector)
-    //     }
-
-    //     return basis
-    // }
-
-    // TODO
-    // /**
-    //  * Get a quaternion descriping the axes rotation
-    //  * @param  {Quaternion} [optionalTarget] - target object
-    //  * @return {Quaternion} the rotation
-    //  */
-    // getRotationQuaternion(optionalTarget = new Quaternion()) {
-    //     const q = optionalTarget
-    //     q.setFromRotationMatrix(this.getBasisMatrix(tmpMatrix))
-
-    //     return q.inverse()
-    // }
-
-    // TODO
-    // /**
-    //  * Get the scale/length for each dimension for a box around the axes
-    //  * to enclose the atoms of a structure
-    //  * @param  {Structure|StructureView} structure - the structure
-    //  * @return {{d1a: Number, d2a: Number, d3a: Number, d1b: Number, d2b: Number, d3b: Number}} scale
-    //  */
-    // getProjectedScaleForAtoms(structure: Structure) {
-    //     let d1a = -Infinity
-    //     let d1b = -Infinity
-    //     let d2a = -Infinity
-    //     let d2b = -Infinity
-    //     let d3a = -Infinity
-    //     let d3b = -Infinity
-
-    //     const p = new Vector3()
-    //     const t = new Vector3()
-
-    //     const center = this.center
-    //     const ax1 = this.normVecA
-    //     const ax2 = this.normVecB
-    //     const ax3 = this.normVecC
-
-    //     structure.eachAtom(function (ap: AtomProxy) {
-    //         projectPointOnVector(p.copy(ap as any), ax1, center)  // TODO
-    //         const dp1 = t.subVectors(p, center).normalize().dot(ax1)
-    //         const dt1 = p.distanceTo(center)
-    //         if (dp1 > 0) {
-    //             if (dt1 > d1a) d1a = dt1
-    //         } else {
-    //             if (dt1 > d1b) d1b = dt1
-    //         }
-
-    //         projectPointOnVector(p.copy(ap as any), ax2, center)
-    //         const dp2 = t.subVectors(p, center).normalize().dot(ax2)
-    //         const dt2 = p.distanceTo(center)
-    //         if (dp2 > 0) {
-    //             if (dt2 > d2a) d2a = dt2
-    //         } else {
-    //             if (dt2 > d2b) d2b = dt2
-    //         }
-
-    //         projectPointOnVector(p.copy(ap as any), ax3, center)
-    //         const dp3 = t.subVectors(p, center).normalize().dot(ax3)
-    //         const dt3 = p.distanceTo(center)
-    //         if (dp3 > 0) {
-    //             if (dt3 > d3a) d3a = dt3
-    //         } else {
-    //             if (dt3 > d3b) d3b = dt3
-    //         }
-    //     })
-
-    //     return {
-    //         d1a: d1a,
-    //         d2a: d2a,
-    //         d3a: d3a,
-    //         d1b: -d1b,
-    //         d2b: -d2b,
-    //         d3b: -d3b
-    //     }
-    // }
+    /**
+     * Set basis matrix for given axes
+     */
+    export function setBasisMatrix(out: Mat4, principalAxes: PrincipalAxes) {
+        Mat4.setAxes(out, principalAxes.normVecB, principalAxes.normVecA, principalAxes.normVecC)
+        if (Mat4.determinant(out) < 0) Mat4.scaleUniformly(out, out, -1)
+        return out
+    }
+
+    /**
+     * Get the scale/length for each dimension for a box around the axes
+     * to enclose the given positions
+     */
+    export function getProjectedScaleForPositions(positions: NumberArray, principalAxes: PrincipalAxes) {
+        let d1a = -Infinity
+        let d1b = -Infinity
+        let d2a = -Infinity
+        let d2b = -Infinity
+        let d3a = -Infinity
+        let d3b = -Infinity
+
+        const p = Vec3()
+        const t = Vec3()
+
+        const { center, normVecA, normVecB, normVecC } = principalAxes
+
+        for (let i = 0, il = positions.length; i < il; i += 3) {
+            Vec3.projectPointOnVector(p, Vec3.fromArray(p, positions, i), normVecA, center)
+            const dp1 = Vec3.dot(normVecA, Vec3.normalize(t, Vec3.sub(t, p, center)))
+            const dt1 = Vec3.distance(p, center)
+            if (dp1 > 0) {
+                if (dt1 > d1a) d1a = dt1
+            } else {
+                if (dt1 > d1b) d1b = dt1
+            }
+
+            Vec3.projectPointOnVector(p, Vec3.fromArray(p, positions, i), normVecB, center)
+            const dp2 = Vec3.dot(normVecB, Vec3.normalize(t, Vec3.sub(t, p, center)))
+            const dt2 = Vec3.distance(p, center)
+            if (dp2 > 0) {
+                if (dt2 > d2a) d2a = dt2
+            } else {
+                if (dt2 > d2b) d2b = dt2
+            }
+
+            Vec3.projectPointOnVector(p, Vec3.fromArray(p, positions, i), normVecC, center)
+            const dp3 = Vec3.dot(normVecC, Vec3.normalize(t, Vec3.sub(t, p, center)))
+            const dt3 = Vec3.distance(p, center)
+            if (dp3 > 0) {
+                if (dt3 > d3a) d3a = dt3
+            } else {
+                if (dt3 > d3b) d3b = dt3
+            }
+        }
+
+        return {
+            d1a: d1a,
+            d2a: d2a,
+            d3a: d3a,
+            d1b: -d1b,
+            d2b: -d2b,
+            d3b: -d3b
+        }
+    }
 }

src/mol-plugin/util/structure-orientation.ts

@@ -4,7 +4,7 @@
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
 
-import { Structure } from '../../mol-model/structure';
+import { Structure, StructureElement } from '../../mol-model/structure';
 import { ShapeRepresentation } from '../../mol-repr/shape/representation';
 import { Shape } from '../../mol-model/shape';
 import { ColorNames } from '../../mol-util/color/names';
@@ -17,25 +17,28 @@ import Matrix from '../../mol-math/linear-algebra/matrix/matrix';
 import { PrincipalAxes } from '../../mol-math/linear-algebra/matrix/principal-axes';
 import { createCage } from '../../mol-geo/primitive/cage';
 import { stringToWords } from '../../mol-util/string';
+import { structureElementStatsLabel } from '../../mol-theme/label';
 
 const tmpMatrixPos = Vec3.zero()
-export function getPositionMatrix(structure: Structure) {
-    const mat = Matrix.create(3, structure.elementCount)
+function getPositions(structure: Structure) {
+    const positions = new Float32Array(structure.elementCount * 3)
     for (let i = 0, m = 0, il = structure.units.length; i < il; ++i) {
         const unit = structure.units[i]
         const { elements } = unit
         const pos = unit.conformation.position
         for (let j = 0, jl = elements.length; j < jl; ++j) {
             pos(elements[j], tmpMatrixPos)
-            Vec3.toArray(tmpMatrixPos, mat.data, m + j * 3)
+            Vec3.toArray(tmpMatrixPos, positions, m + j * 3)
         }
         m += elements.length * 3
     }
-    return mat
+    return positions
 }
 
 interface OrientationData {
+    label: string
     principalAxes: PrincipalAxes
+    projectedScale: { d1a: number, d2a: number, d3a: number, d1b: number, d2b: number, d3b: number }
 }
 
 const OrientationVisuals = { 'principal-axes': '', 'oriented-box': '' }
@@ -44,7 +47,7 @@ const OrientationVisualOptions = Object.keys(OrientationVisuals).map(name => [na
 
 export const OrientationParams = {
     ...Mesh.Params,
-    visuals: PD.MultiSelect<OrientationVisualName>(['principal-axes'], OrientationVisualOptions),
+    visuals: PD.MultiSelect<OrientationVisualName>(['oriented-box'], OrientationVisualOptions),
     orientationColor: PD.Color(ColorNames.orange),
     orientationScale: PD.Numeric(2, { min: 0.1, max: 5, step: 0.1 })
 }
@@ -56,6 +59,11 @@ enum VisualGroup {
     OrientedBox = 2
 }
 
+function getVolume(data: OrientationData) {
+    const { d1a, d2a, d3a, d1b, d2b, d3b } = data.projectedScale
+    return (d1a - d1b) * (d2a - d2b) * (d3a - d3b)
+}
+
 function buildPrincipalAxes(state: MeshBuilder.State, data: OrientationData, props: OrientationProps) {
     const vertices = new Float32Array(6 * 3)
     const edges = new Uint8Array([0, 1, 2, 3, 4, 5])
@@ -69,24 +77,70 @@ function buildPrincipalAxes(state: MeshBuilder.State, data: OrientationData, pro
     const matrix = Mat4.fromTranslation(Mat4(), Vec3.inverse(Vec3(), data.principalAxes.center))
 
     const cage = createCage(vertices, edges)
-    const radius = 0.3 * props.orientationScale
+    const radius = (Math.cbrt(getVolume(data)) / 300) * props.orientationScale
     state.currentGroup = VisualGroup.PrincipalAxes
     MeshBuilder.addCage(state, matrix, cage, radius, 2, 20)
 }
 
+const tmpBoxVec = Vec3()
+function buildOrientedBox(state: MeshBuilder.State, data: OrientationData, props: OrientationProps) {
+    const { center, normVecA, normVecB, normVecC } = data.principalAxes
+    const { d1a, d2a, d3a, d1b, d2b, d3b } = data.projectedScale
+
+    const vertices = new Float32Array(8 * 3)
+    const edges = new Uint8Array([
+        0, 1, 0, 3, 0, 6, 1, 2, 1, 7, 2, 3,
+        2, 4, 3, 5, 4, 5, 4, 7, 5, 6, 6, 7
+    ])
+
+    let offset = 0
+    const addCornerHelper = function (d1: number, d2: number, d3: number) {
+        Vec3.copy(tmpBoxVec, center)
+        Vec3.scaleAndAdd(tmpBoxVec, tmpBoxVec, normVecA, d1)
+        Vec3.scaleAndAdd(tmpBoxVec, tmpBoxVec, normVecB, d2)
+        Vec3.scaleAndAdd(tmpBoxVec, tmpBoxVec, normVecC, d3)
+        Vec3.toArray(tmpBoxVec, vertices, offset)
+        offset += 3
+    }
+    addCornerHelper(d1a, d2a, d3a)
+    addCornerHelper(d1a, d2a, d3b)
+    addCornerHelper(d1a, d2b, d3b)
+    addCornerHelper(d1a, d2b, d3a)
+    addCornerHelper(d1b, d2b, d3b)
+    addCornerHelper(d1b, d2b, d3a)
+    addCornerHelper(d1b, d2a, d3a)
+    addCornerHelper(d1b, d2a, d3b)
+
+    const matrix = Mat4.fromTranslation(Mat4(), Vec3.inverse(Vec3(), data.principalAxes.center))
+
+    const cage = createCage(vertices, edges)
+    const radius = (Math.cbrt(getVolume(data)) / 300) * props.orientationScale
+    state.currentGroup = VisualGroup.OrientedBox
+    MeshBuilder.addCage(state, matrix, cage, radius, 2, 20)
+}
+
 function getOrientationMesh(data: OrientationData, props: OrientationProps, mesh?: Mesh) {
     const state = MeshBuilder.createState(256, 128, mesh)
 
     if (props.visuals.includes('principal-axes')) buildPrincipalAxes(state, data, props)
+    if (props.visuals.includes('oriented-box')) buildOrientedBox(state, data, props)
 
     return MeshBuilder.getMesh(state)
 }
 
+function getLabel(structure: Structure) {
+    const loci = Structure.toStructureElementLoci(structure)
+    const remappedLoci = StructureElement.Loci.remap(loci, structure.root)
+    return structureElementStatsLabel(StructureElement.Stats.ofLoci(remappedLoci), true)
+}
+
 export async function getStructureOrientationRepresentation(ctx: RuntimeContext, structure: Structure, params: OrientationProps, prev?: ShapeRepresentation<OrientationData, Mesh, Mesh.Params>) {
     const repr = prev || ShapeRepresentation(getOrientationShape, Mesh.Utils);
-    const data = {
-        principalAxes: PrincipalAxes.ofPoints(getPositionMatrix(structure))
-    }
+    const label = getLabel(structure)
+    const positions = getPositions(structure)
+    const principalAxes = PrincipalAxes.ofPoints(Matrix.fromArray(positions.slice(), 3, structure.elementCount))
+    const projectedScale = PrincipalAxes.getProjectedScaleForPositions(positions, principalAxes)
+    const data: OrientationData = { label, principalAxes, projectedScale }
     await repr.createOrUpdate(params, data).runInContext(ctx);
     return repr;
 }
@@ -102,7 +156,7 @@ function getOrientationLabel(data: OrientationData, groupId: number): string {
 function getOrientationShape(ctx: RuntimeContext, data: OrientationData, props: OrientationProps, shape?: Shape<Mesh>) {
     const geo = getOrientationMesh(data, props, shape && shape.geometry);
     const getLabel = function (groupId: number ) {
-        return getOrientationLabel(data, groupId)
+        return `${getOrientationLabel(data, groupId)} of ${data.label}`
     }
     return Shape.create('Principal Axes', data, geo, () => props.orientationColor, () => 1, getLabel)
 }