diff --git a/src/mol-math/linear-algebra/3d/mat3.ts b/src/mol-math/linear-algebra/3d/mat3.ts
index 56459c9b01435f4ea17251f58d24452cae9eef2a..d37d4c18cd6debddd01c43f74dbe5e54b2b362f3 100644
--- a/src/mol-math/linear-algebra/3d/mat3.ts
+++ b/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
\ No newline at end of file
diff --git a/src/mol-math/linear-algebra/_spec/mat3.spec.ts b/src/mol-math/linear-algebra/_spec/mat3.spec.ts
new file mode 100644
index 0000000000000000000000000000000000000000..35901555c356fc166634682b23c77ff0012fb3f6
--- /dev/null
+++ b/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);
+ });
+});
\ No newline at end of file