matrix and axes improvements

src/mol-math/linear-algebra/3d.ts

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2017-2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2017-2019 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>
@@ -25,4 +25,10 @@ import Vec4 from './3d/vec4'
 import Quat from './3d/quat'
 import { EPSILON } from './3d/common'
 
-export { Mat4, Mat3, Vec2, Vec3, Vec4, Quat, EPSILON }
+export { Mat4, Mat3, Vec2, Vec3, Vec4, Quat, EPSILON }
+
+export type Vec<T> =
+    T extends 4 ? Vec4 :
+    T extends 3 ? Vec3 :
+    T extends 2 ? Vec2 :
+    number[]

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

@@ -1,10 +1,11 @@
 /**
- * Copyright (c) 2018-2018 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2018-2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
  *
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
 
 import { NumberArray } from '../../../mol-util/type-helpers';
+import { Vec } from '../3d';
 
 interface Matrix<N extends number = number, M extends number = number> {
     data: NumberArray,
@@ -20,14 +21,30 @@ namespace Matrix {
     }
 
     /** Get element assuming data are stored in column-major order */
-    export function get(m: Matrix, i: number, j: number) { return m.data[m.rows * j + i]; }
+    export function get(m: Matrix, i: number, j: number) {
+        return m.data[m.rows * j + i];
+    }
+
     /** Set element assuming data are stored in column-major order */
-    export function set(m: Matrix, i: number, j: number, value: number) { m.data[m.rows * j + i] = value; }
+    export function set<N extends number, M extends number>(m: Matrix<N, M>, i: number, j: number, value: number) {
+        m.data[m.rows * j + i] = value;
+        return m
+    }
+
     /** Add to element assuming data are stored in column-major order */
-    export function add(m: Matrix, i: number, j: number, value: number) { m.data[m.rows * j + i] += value; }
+    export function add<N extends number, M extends number>(m: Matrix<N, M>, i: number, j: number, value: number) {
+        m.data[m.rows * j + i] += value;
+        return m
+    }
+
     /** Zero out the matrix */
-    export function makeZero(m: Matrix) {
-        for (let i = 0, _l = m.data.length; i < _l; i++) m.data[i] = 0.0;
+    export function makeZero<N extends number, M extends number>(m: Matrix<N, M>) {
+        m.data.fill(0.0)
+        return m
+    }
+
+    export function clone<N extends number, M extends number>(m: Matrix<N, M>): Matrix<N, M> {
+        return { data: m.data.slice(), size: m.size, cols: m.cols, rows: m.rows }
     }
 
     export function fromArray<N extends number, M extends number>(data: NumberArray, cols: N, rows: M): Matrix<N, M> {
@@ -36,8 +53,9 @@ namespace Matrix {
 
     export function transpose<N extends number, M extends number>(out: Matrix<M, N>, mat: Matrix<N, M>): Matrix<M, N> {
         const nrows = mat.rows, ncols = mat.cols
+        // TODO add in-place transpose
+        if (out as any === mat) mat = clone(mat)
         const md = mat.data, mtd = out.data
-
         for (let i = 0, mi = 0, mti = 0; i < nrows; mti += 1, mi += ncols, ++i) {
             let ri = mti
             for (let j = 0; j < ncols; ri += nrows, j++) mtd[ri] = md[mi + j]
@@ -46,7 +64,7 @@ namespace Matrix {
     }
 
     /** out = matA * matB' */
-    export function multiplyABt (out: Matrix, matA: Matrix, matB: Matrix) {
+    export function multiplyABt<NA extends number, NB extends number, M extends number>(out: Matrix<M, M>, matA: Matrix<NA, M>, matB: Matrix<NB, M>): Matrix<M, M> {
         const ncols = matA.cols, nrows = matA.rows, mrows = matB.rows
         const ad = matA.data, bd = matB.data, cd = out.data
 
@@ -64,10 +82,10 @@ namespace Matrix {
     }
 
     /** Get the mean of rows in `mat` */
-    export function meanRows (mat: Matrix) {
+    export function meanRows<N extends number, M extends number, V extends Vec<N>>(mat: Matrix<N, M>): V {
         const nrows = mat.rows, ncols = mat.cols
         const md = mat.data
-        const mean = new Array(ncols)
+        const mean = new Array(ncols) as V
 
         for (let j = 0; j < ncols; ++j) mean[ j ] = 0.0
         for (let i = 0, p = 0; i < nrows; ++i) {
@@ -79,7 +97,7 @@ namespace Matrix {
     }
 
     /** Subtract `row` from all rows in `mat` */
-    export function subRows (mat: Matrix, row: NumberArray) {
+    export function subRows<N extends number, M extends number>(mat: Matrix<N, M>, row: NumberArray) {
         const nrows = mat.rows, ncols = mat.cols
         const md = mat.data
 

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

@@ -37,7 +37,6 @@ namespace PrincipalAxes {
     export function ofPoints(points: Matrix<3, number>): PrincipalAxes {
         const n = points.rows
         const n3 = n / 3
-        const pointsT = Matrix.create(n, 3)
         const A = Matrix.create(3, 3)
         const W = Matrix.create(1, 3)
         const U = Matrix.create(3, 3)
@@ -45,8 +44,8 @@ namespace PrincipalAxes {
 
         // calculate
         const mean = Matrix.meanRows(points)
-        Matrix.subRows(points, mean)
-        Matrix.transpose(pointsT, points)
+        const pointsM = Matrix.subRows(Matrix.clone(points), mean)
+        const pointsT = Matrix.transpose(pointsM as Matrix<number, 3>, pointsM)
         Matrix.multiplyABt(A, pointsT, pointsT)
         svd(A, W, U, V)
 
@@ -59,9 +58,9 @@ namespace PrincipalAxes {
         const normVecC = Vec3.create(U.data[2], U.data[5], U.data[8])
 
         // scaled
-        const vecA = Vec3.scale(Vec3.zero(), normVecA, Math.sqrt(W.data[0] / n3))
-        const vecB = Vec3.scale(Vec3.zero(), normVecB, Math.sqrt(W.data[1] / n3))
-        const vecC = Vec3.scale(Vec3.zero(), normVecC, Math.sqrt(W.data[2] / n3))
+        const vecA = Vec3.scale(Vec3(), normVecA, Math.sqrt(W.data[0] / n3))
+        const vecB = Vec3.scale(Vec3(), normVecB, Math.sqrt(W.data[1] / n3))
+        const vecC = Vec3.scale(Vec3(), normVecC, Math.sqrt(W.data[2] / n3))
 
         // points
         const begA = Vec3.sub(Vec3.clone(center), center, vecA)
@@ -92,7 +91,7 @@ namespace PrincipalAxes {
      * 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) {
+    export function getProjectedScale(positions: NumberArray, principalAxes: PrincipalAxes) {
         let d1a = -Infinity
         let d1b = -Infinity
         let d2a = -Infinity

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

@@ -138,8 +138,8 @@ export async function getStructureOrientationRepresentation(ctx: RuntimeContext,
     const repr = prev || ShapeRepresentation(getOrientationShape, Mesh.Utils);
     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 principalAxes = PrincipalAxes.ofPoints(Matrix.fromArray(positions, 3, structure.elementCount))
+    const projectedScale = PrincipalAxes.getProjectedScale(positions, principalAxes)
     const data: OrientationData = { label, principalAxes, projectedScale }
     await repr.createOrUpdate(params, data).runInContext(ctx);
     return repr;