add euler math primitive

CHANGELOG.md

@@ -8,6 +8,7 @@ Note that since we don't clearly distinguish between a public and private interf
 
 - Fix display issue with SIFTS mapping
 - Properly switch-off fog
+- Add `Euler` math primitive
 
 ## [v3.37.1] - 2023-06-20
 

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

@@ -0,0 +1,167 @@
+/**
+ * Copyright (c) 2023 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ *
+ * This code has been modified from https://github.com/mrdoob/three.js/,
+ * copyright (c) 2010-2023 three.js authors. MIT License
+ */
+
+import { Mat4 } from './mat4';
+import { assertUnreachable, NumberArray } from '../../../mol-util/type-helpers';
+import { Quat } from './quat';
+import { Vec3 } from './vec3';
+import { clamp } from '../../interpolate';
+
+interface Euler extends Array<number> { [d: number]: number, '@type': 'euler', length: 3 }
+
+function Euler() {
+    return Euler.zero();
+}
+
+namespace Euler {
+    export type Order = 'XYZ' | 'YXZ' | 'ZXY' | 'ZYX' | 'YZX' | 'XZY'
+
+    export function zero(): Euler {
+        // force double backing array by 0.1.
+        const ret = [0.1, 0, 0];
+        ret[0] = 0.0;
+        return ret as any;
+    }
+
+    export function create(x: number, y: number, z: number): Euler {
+        const out = zero();
+        out[0] = x;
+        out[1] = y;
+        out[2] = z;
+        return out;
+    }
+
+    export function set(out: Euler, x: number, y: number, z: number) {
+        out[0] = x;
+        out[0] = y;
+        out[0] = z;
+        return out;
+    }
+
+    export function clone(a: Euler): Euler {
+        const out = zero();
+        out[0] = a[0];
+        out[1] = a[1];
+        out[2] = a[2];
+        return out;
+    }
+
+    export function copy(out: Euler, a: Euler) {
+        out[0] = a[0];
+        out[1] = a[1];
+        out[2] = a[2];
+        return out;
+    }
+
+    /**
+     * Assumes the upper 3x3 of m is a pure rotation matrix (i.e, unscaled)
+     */
+    export function fromMat4(out: Euler, m: Mat4, order: Order): Euler {
+        const m11 = m[0], m12 = m[4], m13 = m[8];
+        const m21 = m[1], m22 = m[5], m23 = m[9];
+        const m31 = m[2], m32 = m[6], m33 = m[10];
+
+        switch (order) {
+            case 'XYZ':
+                out[1] = Math.asin(clamp(m13, -1, 1));
+                if (Math.abs(m13) < 0.9999999) {
+                    out[0] = Math.atan2(-m23, m33);
+                    out[2] = Math.atan2(-m12, m11);
+                } else {
+                    out[0] = Math.atan2(m32, m22);
+                    out[2] = 0;
+                }
+                break;
+            case 'YXZ':
+                out[0] = Math.asin(-clamp(m23, -1, 1));
+                if (Math.abs(m23) < 0.9999999) {
+                    out[1] = Math.atan2(m13, m33);
+                    out[2] = Math.atan2(m21, m22);
+                } else {
+                    out[1] = Math.atan2(-m31, m11);
+                    out[2] = 0;
+                }
+                break;
+            case 'ZXY':
+                out[0] = Math.asin(clamp(m32, -1, 1));
+                if (Math.abs(m32) < 0.9999999) {
+                    out[1] = Math.atan2(-m31, m33);
+                    out[2] = Math.atan2(-m12, m22);
+                } else {
+                    out[1] = 0;
+                    out[2] = Math.atan2(m21, m11);
+                }
+                break;
+            case 'ZYX':
+                out[1] = Math.asin(-clamp(m31, -1, 1));
+                if (Math.abs(m31) < 0.9999999) {
+                    out[0] = Math.atan2(m32, m33);
+                    out[2] = Math.atan2(m21, m11);
+                } else {
+                    out[0] = 0;
+                    out[2] = Math.atan2(-m12, m22);
+                }
+                break;
+            case 'YZX':
+                out[2] = Math.asin(clamp(m21, -1, 1));
+                if (Math.abs(m21) < 0.9999999) {
+                    out[0] = Math.atan2(-m23, m22);
+                    out[1] = Math.atan2(-m31, m11);
+                } else {
+                    out[0] = 0;
+                    out[1] = Math.atan2(m13, m33);
+                }
+                break;
+            case 'XZY':
+                out[2] = Math.asin(-clamp(m12, -1, 1));
+                if (Math.abs(m12) < 0.9999999) {
+                    out[0] = Math.atan2(m32, m22);
+                    out[1] = Math.atan2(m13, m11);
+                } else {
+                    out[0] = Math.atan2(-m23, m33);
+                    out[1] = 0;
+                }
+                break;
+            default:
+                assertUnreachable(order);
+        }
+
+        return out;
+    }
+
+    const _mat4 = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] as Mat4;
+    export function fromQuat(out: Euler, q: Quat, order: Order) {
+        Mat4.fromQuat(_mat4, q);
+        return fromMat4(out, _mat4, order);
+    }
+
+    export function fromVec3(out: Euler, v: Vec3) {
+        return set(out, v[0], v[1], v[2]);
+    }
+
+    export function exactEquals(a: Euler, b: Euler) {
+        return a[0] === b[0] && a[1] === b[1] && a[2] === b[2];
+    }
+
+    export function fromArray(e: Euler, array: ArrayLike<number>, offset: number) {
+        e[0] = array[offset + 0];
+        e[1] = array[offset + 1];
+        e[2] = array[offset + 2];
+        return e;
+    }
+
+    export function toArray<T extends NumberArray>(e: Euler, out: T, offset: number) {
+        out[offset + 0] = e[0];
+        out[offset + 1] = e[1];
+        out[offset + 2] = e[2];
+        return out;
+    }
+}
+
+export { Euler };

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

@@ -23,6 +23,7 @@ import { Quat } from './quat';
 import { degToRad } from '../../misc';
 import { NumberArray } from '../../../mol-util/type-helpers';
 import { Mat3 } from './mat3';
+import { Euler } from './euler';
 
 interface Mat4 extends Array<number> { [d: number]: number, '@type': 'mat4', length: 16 }
 interface ReadonlyMat4 extends Array<number> { readonly [d: number]: number, '@type': 'mat4', length: 16 }
@@ -717,6 +718,82 @@ namespace Mat4 {
         return out;
     }
 
+    export function compose(out: Mat4, position: Vec3, quaternion: Quat, scale: Vec3) {
+        const [x, y, z, w] = quaternion;
+        const x2 = x + x,	y2 = y + y, z2 = z + z;
+        const xx = x * x2, xy = x * y2, xz = x * z2;
+        const yy = y * y2, yz = y * z2, zz = z * z2;
+        const wx = w * x2, wy = w * y2, wz = w * z2;
+
+        const [sx, sy, sz] = scale;
+
+        out[0] = (1 - (yy + zz)) * sx;
+        out[1] = (xy + wz) * sx;
+        out[2] = (xz - wy) * sx;
+        out[3] = 0;
+
+        out[4] = (xy - wz) * sy;
+        out[5] = (1 - (xx + zz)) * sy;
+        out[6] = (yz + wx) * sy;
+        out[7] = 0;
+
+        out[8] = (xz + wy) * sz;
+        out[9] = (yz - wx) * sz;
+        out[10] = (1 - (xx + yy)) * sz;
+        out[11] = 0;
+
+        out[12] = position[0];
+        out[13] = position[1];
+        out[14] = position[2];
+        out[15] = 1;
+
+        return out;
+    }
+
+    const _v3 = [0, 0, 0] as Vec3;
+    const _m4 = zero();
+    export function decompose(m: Mat4, position: Vec3, quaternion: Quat, scale: Vec3) {
+
+        let sx = Vec3.magnitude(Vec3.set(_v3, m[0], m[1], m[2]));
+        const sy = Vec3.magnitude(Vec3.set(_v3, m[4], m[5], m[6]));
+        const sz = Vec3.magnitude(Vec3.set(_v3, m[8], m[9], m[10]));
+
+        // if determine is negative, we need to invert one scale
+        const det = determinant(m);
+        if (det < 0) sx = -sx;
+
+        position[0] = m[12];
+        position[1] = m[13];
+        position[2] = m[14];
+
+        // scale the rotation part
+        copy(_m4, m);
+
+        const invSX = 1 / sx;
+        const invSY = 1 / sy;
+        const invSZ = 1 / sz;
+
+        _m4[0] *= invSX;
+        _m4[1] *= invSX;
+        _m4[2] *= invSX;
+
+        _m4[4] *= invSY;
+        _m4[5] *= invSY;
+        _m4[6] *= invSY;
+
+        _m4[8] *= invSZ;
+        _m4[9] *= invSZ;
+        _m4[10] *= invSZ;
+
+        getRotation(quaternion, _m4);
+
+        scale[0] = sx;
+        scale[1] = sy;
+        scale[2] = sz;
+
+        return m;
+    }
+
     export function makeTable(m: Mat4) {
         let ret = '';
         for (let i = 0; i < 4; i++) {
@@ -851,6 +928,94 @@ namespace Mat4 {
         return out;
     }
 
+    export function fromEuler(out: Mat4, euler: Euler, order: Euler.Order) {
+        const x = euler[0], y = euler[1], z = euler[2];
+        const a = Math.cos(x), b = Math.sin(x);
+        const c = Math.cos(y), d = Math.sin(y);
+        const e = Math.cos(z), f = Math.sin(z);
+
+        if (order === 'XYZ') {
+            const ae = a * e, af = a * f, be = b * e, bf = b * f;
+            out[0] = c * e;
+            out[4] = - c * f;
+            out[8] = d;
+            out[1] = af + be * d;
+            out[5] = ae - bf * d;
+            out[9] = - b * c;
+            out[2] = bf - ae * d;
+            out[6] = be + af * d;
+            out[10] = a * c;
+        } else if (order === 'YXZ') {
+            const ce = c * e, cf = c * f, de = d * e, df = d * f;
+            out[0] = ce + df * b;
+            out[4] = de * b - cf;
+            out[8] = a * d;
+            out[1] = a * f;
+            out[5] = a * e;
+            out[9] = - b;
+            out[2] = cf * b - de;
+            out[6] = df + ce * b;
+            out[10] = a * c;
+        } else if (order === 'ZXY') {
+            const ce = c * e, cf = c * f, de = d * e, df = d * f;
+            out[0] = ce - df * b;
+            out[4] = - a * f;
+            out[8] = de + cf * b;
+            out[1] = cf + de * b;
+            out[5] = a * e;
+            out[9] = df - ce * b;
+            out[2] = - a * d;
+            out[6] = b;
+            out[10] = a * c;
+        } else if (order === 'ZYX') {
+            const ae = a * e, af = a * f, be = b * e, bf = b * f;
+            out[0] = c * e;
+            out[4] = be * d - af;
+            out[8] = ae * d + bf;
+            out[1] = c * f;
+            out[5] = bf * d + ae;
+            out[9] = af * d - be;
+            out[2] = - d;
+            out[6] = b * c;
+            out[10] = a * c;
+        } else if (order === 'YZX') {
+            const ac = a * c, ad = a * d, bc = b * c, bd = b * d;
+            out[0] = c * e;
+            out[4] = bd - ac * f;
+            out[8] = bc * f + ad;
+            out[1] = f;
+            out[5] = a * e;
+            out[9] = - b * e;
+            out[2] = - d * e;
+            out[6] = ad * f + bc;
+            out[10] = ac - bd * f;
+        } else if (order === 'XZY') {
+            const ac = a * c, ad = a * d, bc = b * c, bd = b * d;
+            out[0] = c * e;
+            out[4] = - f;
+            out[8] = d * e;
+            out[1] = ac * f + bd;
+            out[5] = a * e;
+            out[9] = ad * f - bc;
+            out[2] = bc * f - ad;
+            out[6] = b * e;
+            out[10] = bd * f + ac;
+        }
+
+        // bottom row
+        out[3] = 0;
+        out[7] = 0;
+        out[11] = 0;
+
+        // last column
+        out[12] = 0;
+        out[13] = 0;
+        out[14] = 0;
+        out[15] = 1;
+
+        return out;
+    }
+
     /**
      * Generates a perspective projection (frustum) matrix with the given bounds
      */

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

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2017-2021 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2017-2023 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,7 +25,8 @@
 import { Mat3 } from './mat3';
 import { Vec3 } from './vec3';
 import { EPSILON } from './common';
-import { NumberArray } from '../../../mol-util/type-helpers';
+import { assertUnreachable, NumberArray } from '../../../mol-util/type-helpers';
+import { Euler } from './euler';
 
 interface Quat extends Array<number> { [d: number]: number, '@type': 'quat', length: 4 }
 interface ReadonlyQuat extends Array<number> { readonly [d: number]: number, '@type': 'quat', length: 4 }
@@ -238,6 +239,10 @@ namespace Quat {
         return out;
     }
 
+    export function dot(a: Quat, b: Quat) {
+        return a[0] * b[0] + a[1] * b[1] + a[2] * b[2] + a[3] * b[3];
+    }
+
     /**
      * Creates a quaternion from the given 3x3 rotation matrix.
      *
@@ -277,6 +282,63 @@ namespace Quat {
         return out;
     }
 
+    export function fromEuler(out: Quat, euler: Euler, order: Euler.Order) {
+        const [x, y, z] = euler;
+
+        // http://www.mathworks.com/matlabcentral/fileexchange/20696-function-to-convert-between-dcm-euler-angles-quaternions-and-euler-vectors/content/SpinCalc.m
+
+        const c1 = Math.cos(x / 2);
+        const c2 = Math.cos(y / 2);
+        const c3 = Math.cos(z / 2);
+
+        const s1 = Math.sin(x / 2);
+        const s2 = Math.sin(y / 2);
+        const s3 = Math.sin(z / 2);
+
+        switch (order) {
+            case 'XYZ':
+                out[0] = s1 * c2 * c3 + c1 * s2 * s3;
+                out[1] = c1 * s2 * c3 - s1 * c2 * s3;
+                out[2] = c1 * c2 * s3 + s1 * s2 * c3;
+                out[3] = c1 * c2 * c3 - s1 * s2 * s3;
+                break;
+            case 'YXZ':
+                out[0] = s1 * c2 * c3 + c1 * s2 * s3;
+                out[1] = c1 * s2 * c3 - s1 * c2 * s3;
+                out[2] = c1 * c2 * s3 - s1 * s2 * c3;
+                out[3] = c1 * c2 * c3 + s1 * s2 * s3;
+                break;
+            case 'ZXY':
+                out[0] = s1 * c2 * c3 - c1 * s2 * s3;
+                out[1] = c1 * s2 * c3 + s1 * c2 * s3;
+                out[2] = c1 * c2 * s3 + s1 * s2 * c3;
+                out[3] = c1 * c2 * c3 - s1 * s2 * s3;
+                break;
+            case 'ZYX':
+                out[0] = s1 * c2 * c3 - c1 * s2 * s3;
+                out[1] = c1 * s2 * c3 + s1 * c2 * s3;
+                out[2] = c1 * c2 * s3 - s1 * s2 * c3;
+                out[3] = c1 * c2 * c3 + s1 * s2 * s3;
+                break;
+            case 'YZX':
+                out[0] = s1 * c2 * c3 + c1 * s2 * s3;
+                out[1] = c1 * s2 * c3 + s1 * c2 * s3;
+                out[2] = c1 * c2 * s3 - s1 * s2 * c3;
+                out[3] = c1 * c2 * c3 - s1 * s2 * s3;
+                break;
+            case 'XZY':
+                out[0] = s1 * c2 * c3 - c1 * s2 * s3;
+                out[1] = c1 * s2 * c3 - s1 * c2 * s3;
+                out[2] = c1 * c2 * s3 + s1 * s2 * c3;
+                out[3] = c1 * c2 * c3 + s1 * s2 * s3;
+                break;
+            default:
+                assertUnreachable(order);
+        }
+
+        return out;
+    }
+
     const fromUnitVec3Temp = [0, 0, 0] as Vec3;
     /** Quaternion from two normalized unit vectors. */
     export function fromUnitVec3(out: Quat, a: Vec3, b: Vec3) {

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

@@ -0,0 +1,35 @@
+/**
+ * Copyright (c) 2023 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { Mat4 } from '../3d/mat4';
+import { Euler } from '../3d/euler';
+import { Quat } from '../3d/quat';
+
+const t = [
+    [Euler.create(0, 0, 0), 'XYZ'],
+    [Euler.create(1, 0, 0), 'XYZ'],
+    [Euler.create(0, 1, 0), 'ZYX'],
+] as const;
+
+describe('Euler', () => {
+    it('fromMat4', () => {
+        for (const [e, o] of t) {
+            const m = Mat4.fromEuler(Mat4(), e, o);
+            const e2 = Euler.fromMat4(Euler(), m, o);
+            const m2 = Mat4.fromEuler(Mat4(), e2, o);
+            expect(Mat4.areEqual(m, m2, 0.0001)).toBe(true);
+        }
+    });
+
+    it('fromQuat', () => {
+        for (const [e, o] of t) {
+            const q = Quat.fromEuler(Quat(), e, o);
+            const e2 = Euler.fromQuat(Euler(), q, o);
+            const q2 = Quat.fromEuler(Quat(), e2, o);
+            expect(Quat.equals(q, q2)).toBe(true);
+        }
+    });
+});