mat3 symmetricEigenvalues and eigenvectors

src/mol-math/linear-algebra/3d/mat3.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>
@@ -17,10 +17,11 @@
  * furnished to do so, subject to the following conditions:
  */
 
-import { Mat4 } from '../3d'
+import { Mat4, Vec3 } from '../3d'
 import { NumberArray } from '../../../mol-util/type-helpers';
 
 interface Mat3 extends Array<number> { [d: number]: number, '@type': 'mat3', length: 9 }
+interface ReadonlyMat3 extends Array<number> { readonly [d: number]: number, '@type': 'mat3', length: 9 }
 
 function Mat3() {
     return Mat3.zero();
@@ -102,6 +103,20 @@ namespace Mat3 {
         return out;
     }
 
+    export function create(a00: number, a01: number, a02: number, a10: number, a11: number, a12: number, a20: number, a21: number, a22: number): Mat3 {
+        const out = zero();
+        out[0] = a00
+        out[1] = a01
+        out[2] = a02
+        out[3] = a10
+        out[4] = a11
+        out[5] = a12
+        out[6] = a20
+        out[7] = a21
+        out[8] = a22
+        return out;
+    }
+
     export function hasNaN(m: Mat3) {
         for (let i = 0; i < 9; i++) if (isNaN(m[i])) return true
         return false
@@ -114,6 +129,13 @@ namespace Mat3 {
         return Mat3.copy(Mat3.zero(), a);
     }
 
+    export function areEqual(a: Mat3, b: Mat3, eps: number) {
+        for (let i = 0; i < 9; i++) {
+            if (Math.abs(a[i] - b[i]) > eps) return false;
+        }
+        return true;
+    }
+
     export function setValue(a: Mat3, i: number, j: number, value: number) {
         a[3 * j + i] = value;
     }
@@ -210,6 +232,172 @@ namespace Mat3 {
         // Calculate the determinant
         return a00 * b01 + a01 * b11 + a02 * b21;
     }
+
+    export function trace(a: Mat3) {
+        return a[0] + a[4] + a[8]
+    }
+
+    export function sub(out: Mat3, a: Mat3, b: Mat3) {
+        out[0] = a[0] - b[0]
+        out[1] = a[1] - b[1]
+        out[2] = a[2] - b[2]
+        out[3] = a[3] - b[3]
+        out[4] = a[4] - b[4]
+        out[5] = a[5] - b[5]
+        out[6] = a[6] - b[6]
+        out[7] = a[7] - b[7]
+        out[8] = a[8] - b[8]
+        return out
+    }
+
+    export function add(out: Mat3, a: Mat3, b: Mat3) {
+        out[0] = a[0] + b[0]
+        out[1] = a[1] + b[1]
+        out[2] = a[2] + b[2]
+        out[3] = a[3] + b[3]
+        out[4] = a[4] + b[4]
+        out[5] = a[5] + b[5]
+        out[6] = a[6] + b[6]
+        out[7] = a[7] + b[7]
+        out[8] = a[8] + b[8]
+        return out
+    }
+
+    export function mul(out: Mat3, a: Mat3, b: Mat3) {
+        const a00 = a[0], a01 = a[1], a02 = a[2],
+            a10 = a[3], a11 = a[4], a12 = a[5],
+            a20 = a[6], a21 = a[7], a22 = a[8];
+
+        const b00 = b[0], b01 = b[1], b02 = b[2],
+            b10 = b[3], b11 = b[4], b12 = b[5],
+            b20 = b[6], b21 = b[7], b22 = b[8];
+
+        out[0] = b00 * a00 + b01 * a10 + b02 * a20;
+        out[1] = b00 * a01 + b01 * a11 + b02 * a21;
+        out[2] = b00 * a02 + b01 * a12 + b02 * a22;
+
+        out[3] = b10 * a00 + b11 * a10 + b12 * a20;
+        out[4] = b10 * a01 + b11 * a11 + b12 * a21;
+        out[5] = b10 * a02 + b11 * a12 + b12 * a22;
+
+        out[6] = b20 * a00 + b21 * a10 + b22 * a20;
+        out[7] = b20 * a01 + b21 * a11 + b22 * a21;
+        out[8] = b20 * a02 + b21 * a12 + b22 * a22;
+        return out;
+    }
+
+    export function subScalar(out: Mat3, a: Mat3, s: number) {
+        out[0] = a[0] - s
+        out[1] = a[1] - s
+        out[2] = a[2] - s
+        out[3] = a[3] - s
+        out[4] = a[4] - s
+        out[5] = a[5] - s
+        out[6] = a[6] - s
+        out[7] = a[7] - s
+        out[8] = a[8] - s
+        return out
+    }
+
+    export function addScalar(out: Mat3, a: Mat3, s: number) {
+        out[0] = a[0] + s
+        out[1] = a[1] + s
+        out[2] = a[2] + s
+        out[3] = a[3] + s
+        out[4] = a[4] + s
+        out[5] = a[5] + s
+        out[6] = a[6] + s
+        out[7] = a[7] + s
+        out[8] = a[8] + s
+        return out
+    }
+
+    export function mulScalar(out: Mat3, a: Mat3, s: number) {
+        out[0] = a[0] * s
+        out[1] = a[1] * s
+        out[2] = a[2] * s
+        out[3] = a[3] * s
+        out[4] = a[4] * s
+        out[5] = a[5] * s
+        out[6] = a[6] * s
+        out[7] = a[7] * s
+        out[8] = a[8] * s
+        return out
+    }
+
+    const piThird = Math.PI / 3
+    const tmpB = Mat3()
+    /**
+     * Given a real symmetric 3x3 matrix A, compute the eigenvalues
+     *
+     * From https://en.wikipedia.org/wiki/Eigenvalue_algorithm#3.C3.973_matrices
+     */
+    export function symmetricEigenvalues(out: Vec3, a: Mat3) {
+        const p1 = a[1] * a[1] + a[2] * a[2] + a[5] * a[5]
+        if (p1 === 0) {
+            out[0] = a[0]
+            out[1] = a[4]
+            out[2] = a[8]
+        } else {
+            const q = trace(a) / 3
+            const a1 = a[0] - q
+            const a2 = a[4] - q
+            const a3 = a[8] - q
+            const p2 = a1 * a1 + a2 * a2 + a3 * a3 + 2 * p1
+            const p = Math.sqrt(p2 / 6)
+            mulScalar(tmpB, Identity, q)
+            sub(tmpB, a, tmpB)
+            mulScalar(tmpB, tmpB, (1 / p))
+            const r = determinant(tmpB) / 2
+            // In exact arithmetic for a symmetric matrix  -1 <= r <= 1
+            // but computation error can leave it slightly outside this range.
+            const phi = r <= -1 ? piThird : r >= 1 ?
+                0 : Math.acos(r) / 3
+            // the eigenvalues satisfy eig3 <= eig2 <= eig1
+            out[0] = q + 2 * p * Math.cos(phi)
+            out[2] = q + 2 * p * Math.cos(phi + (2 * piThird))
+            out[1] = 3 * q - out[0] - out[2] // since trace(A) = eig1 + eig2 + eig3
+        }
+        return out
+    }
+
+    const tmpR0 = [0.1, 0.0, 0.0] as Vec3
+    const tmpR1 = [0.1, 0.0, 0.0] as Vec3
+    const tmpR2 = [0.1, 0.0, 0.0] as Vec3
+    const tmpR0xR1 = [0.1, 0.0, 0.0] as Vec3
+    const tmpR0xR2 = [0.1, 0.0, 0.0] as Vec3
+    const tmpR1xR2 = [0.1, 0.0, 0.0] as Vec3
+    /**
+     * Calculates the eigenvector for the given eigenvalue `e` of matrix `a`
+     */
+    export function eigenvector(out: Vec3, a: Mat3, e: number) {
+        Vec3.set(tmpR0, a[0] - e, a[1], a[2])
+        Vec3.set(tmpR1, a[1], a[4] - e, a[5])
+        Vec3.set(tmpR2, a[2], a[5], a[8] - e)
+        Vec3.cross(tmpR0xR1, tmpR0, tmpR1)
+        Vec3.cross(tmpR0xR2, tmpR0, tmpR2)
+        Vec3.cross(tmpR1xR2, tmpR1, tmpR2)
+        const d0 = Vec3.dot(tmpR0xR1, tmpR0xR1)
+        const d1 = Vec3.dot(tmpR0xR2, tmpR0xR2)
+        const d2 = Vec3.dot(tmpR1xR2, tmpR1xR2)
+        let dmax = d0
+        let imax = 0
+        if (d1 > dmax) {
+            dmax = d1
+            imax = 1
+        }
+        if (d2 > dmax) imax = 2
+        if (imax === 0) {
+            Vec3.scale(out, tmpR0xR1, 1 / Math.sqrt(d0))
+        } else if (imax === 1) {
+            Vec3.scale(out, tmpR0xR2, 1 / Math.sqrt(d1))
+        } else {
+            Vec3.scale(out, tmpR1xR2, 1 / Math.sqrt(d2))
+        }
+        return out
+    }
+
+    export const Identity: ReadonlyMat3 = identity()
 }
 
 export default Mat3

src/mol-math/linear-algebra/_spec/mat3.spec.ts

@@ -0,0 +1,35 @@
+/**
+ * Copyright (c) 2019 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { Mat3, Vec3 } from '../3d'
+
+describe('Mat3', () => {
+    it('symmetricEigenvalues', () => {
+        const m = Mat3.create(
+            0.1945, -0.0219, -0.0416,
+            -0.0219, 0.1995, -0.0119,
+            -0.0416, -0.0119, 0.3673
+        )
+        const e = Vec3.create(0.377052701425898, 0.21713981522725134, 0.1671074833468507)
+        expect(Vec3.equals(e, Mat3.symmetricEigenvalues(Vec3(), m))).toBe(true);
+    });
+
+    it('eigenvectors', () => {
+        const m = Mat3.create(
+            0.1945, -0.0219, -0.0416,
+            -0.0219, 0.1995, -0.0119,
+            -0.0416, -0.0119, 0.3673
+        )
+        const e = Vec3.create(0.377052701425898, 0.21713981522725134, 0.1671074833468507)
+        const v0 = Vec3.create(-0.2176231019882068, -0.038522620041966125, 0.9752723687391808)
+        const v1 = Vec3.create(-0.5905636938047126, 0.8007524989198634, -0.10014968314142503)
+        const v2 = Vec3.create(0.7770937582036648, 0.5977553372576602, 0.19701230352667118)
+
+        expect(Vec3.equals(v0, Mat3.eigenvector(Vec3(), m, e[0]))).toBe(true);
+        expect(Vec3.equals(v1, Mat3.eigenvector(Vec3(), m, e[1]))).toBe(true);
+        expect(Vec3.equals(v2, Mat3.eigenvector(Vec3(), m, e[2]))).toBe(true);
+    });
+});