add Frustum3D and Plane3D math primitives

CHANGELOG.md

@@ -18,6 +18,7 @@ Note that since we don't clearly distinguish between a public and private interf
 - Remove `JSX` reference from `loci-labels.ts`
 - Fix overpaint/transparency/substance smoothing not updated when geometry changes
 - Fix camera project/unproject when using offset viewport
+- Add `Frustum3D` and `Plane3D` math primitives
 
 ## [v3.32.0] - 2023-03-20
 

src/mol-math/geometry/_spec/frustum3d.spec.ts

@@ -0,0 +1,74 @@
+/**
+ * Copyright (c) 2023 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { Mat4, Vec3 } from '../../linear-algebra';
+import { Box3D } from '../primitives/box3d';
+import { Frustum3D } from '../primitives/frustum3d';
+import { Sphere3D } from '../primitives/sphere3d';
+
+const v3 = Vec3.create;
+const s3 = Sphere3D.create;
+
+describe('frustum3d', () => {
+    it('intersectsSphere3D', () => {
+        const f = Frustum3D();
+        const m = Mat4.perspective(Mat4(), -1, 1, 1, -1, 1, 100);
+        Frustum3D.fromProjectionMatrix(f, m);
+
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(0, 0, 0), 0))).toBe(false);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(0, 0, 0), 0.9))).toBe(false);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(0, 0, 0), 1.1))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(0, 0, -50), 0))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(0, 0, -1.001), 0))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(-1, -1, -1.001), 0))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(-1.1, -1.1, -1.001), 0))).toBe(false);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(-1.1, -1.1, -1.001), 0.5))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(1, 1, -1.001), 0))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(1.1, 1.1, -1.001), 0))).toBe(false);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(1.1, 1.1, -1.001), 0.5))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(0, 0, -99.999), 0))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(-99.999, -99.999, -99.999), 0))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(-100.1, -100.1, -100.1), 0))).toBe(false);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(-100.1, -100.1, -100.1), 0.5))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(99.999, 99.999, -99.999), 0))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(100.1, 100.1, -100.1), 0))).toBe(false);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(100.1, 100.1, -100.1), 0.2))).toBe(true);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(0, 0, -101), 0))).toBe(false);
+        expect(Frustum3D.intersectsSphere3D(f, s3(v3(0, 0, -101), 1.1))).toBe(true);
+    });
+
+    it('intersectsBox3D', () => {
+        const f = Frustum3D();
+        const m = Mat4.perspective(Mat4(), -1, 1, 1, -1, 1, 100);
+        Frustum3D.fromProjectionMatrix(f, m);
+
+        const b0 = Box3D.create(v3(0, 0, 0), v3(1, 1, 1));
+        expect(Frustum3D.intersectsBox3D(f, b0)). toBe(false);
+
+        const b1 = Box3D.create(v3(-1.1, -1.1, -1.1), v3(-0.1, -0.1, -0.1));
+        expect(Frustum3D.intersectsBox3D(f, b1)). toBe(true);
+    });
+
+    it('containsPoint', () => {
+        const f = Frustum3D();
+        const m = Mat4.perspective(Mat4(), -1, 1, 1, -1, 1, 100);
+        Frustum3D.fromProjectionMatrix(f, m);
+
+        expect(Frustum3D.containsPoint(f, v3(0, 0, 0))).toBe(false);
+        expect(Frustum3D.containsPoint(f, v3(0, 0, -50))).toBe(true);
+        expect(Frustum3D.containsPoint(f, v3(0, 0, -1.001))).toBe(true);
+        expect(Frustum3D.containsPoint(f, v3(-1, -1, -1.001))).toBe(true);
+        expect(Frustum3D.containsPoint(f, v3(-1.1, -1.1, -1.001))).toBe(false);
+        expect(Frustum3D.containsPoint(f, v3(1, 1, -1.001))).toBe(true);
+        expect(Frustum3D.containsPoint(f, v3(1.1, 1.1, -1.001))).toBe(false);
+        expect(Frustum3D.containsPoint(f, v3(0, 0, -99.999))).toBe(true);
+        expect(Frustum3D.containsPoint(f, v3(-99.999, -99.999, -99.999))).toBe(true);
+        expect(Frustum3D.containsPoint(f, v3(-100.1, -100.1, -100.1))).toBe(false);
+        expect(Frustum3D.containsPoint(f, v3(99.999, 99.999, -99.999))).toBe(true);
+        expect(Frustum3D.containsPoint(f, v3(100.1, 100.1, -100.1))).toBe(false);
+        expect(Frustum3D.containsPoint(f, v3(0, 0, -101))).toBe(false);
+    });
+});

src/mol-math/geometry/_spec/plane3d.spec.ts

@@ -0,0 +1,40 @@
+/**
+ * Copyright (c) 2023 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { Vec3 } from '../../linear-algebra';
+import { Plane3D } from '../primitives/plane3d';
+
+describe('plane3d', () => {
+    it('fromNormalAndCoplanarPoint', () => {
+        const normal = Vec3.create(1, 1, 1);
+        Vec3.normalize(normal, normal);
+        const p = Plane3D();
+        Plane3D.fromNormalAndCoplanarPoint(p, normal, Vec3.zero());
+
+        expect(p.normal).toEqual(normal);
+        expect(p.constant).toBe(-0);
+    });
+
+    it('fromCoplanarPoints', () => {
+        const a = Vec3.create(2.0, 0.5, 0.25);
+        const b = Vec3.create(2.0, -0.5, 1.25);
+        const c = Vec3.create(2.0, -3.5, 2.2);
+        const p = Plane3D();
+        Plane3D.fromCoplanarPoints(p, a, b, c);
+
+        expect(p.normal).toEqual(Vec3.create(1, 0, 0));
+        expect(p.constant).toBe(-2);
+    });
+
+    it('distanceToPoint', () => {
+        const p = Plane3D.create(Vec3.create(2, 0, 0), -2);
+        Plane3D.normalize(p, p);
+
+        expect(Plane3D.distanceToPoint(p, Vec3.create(0, 0, 0))).toBe(-1);
+        expect(Plane3D.distanceToPoint(p, Vec3.create(4, 0, 0))).toBe(3);
+        expect(Plane3D.distanceToPoint(p, Plane3D.projectPoint(Vec3(), p, Vec3.zero()))).toBe(0);
+    });
+});

src/mol-math/geometry/_spec/polygon.spec.ts

@@ -0,0 +1,21 @@
+/**
+ * Copyright (c) 2023 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ *
+ * @author Alexander Rose <alexander.rose@weirdbyte.de>
+ */
+
+import { Vec2 } from '../../linear-algebra';
+import { pointInPolygon } from '../polygon';
+
+describe('pointInPolygon', () => {
+    it('basic', () => {
+        const polygon = [
+            -1, -1,
+            1, -1,
+            1, 1,
+            -1, 1
+        ];
+        expect(pointInPolygon(Vec2.create(0, 0), polygon, 4)).toBe(true);
+        expect(pointInPolygon(Vec2.create(2, 2), polygon, 4)).toBe(false);
+    });
+});

src/mol-math/geometry/polygon.ts

@@ -0,0 +1,24 @@
+/**
+ * Copyright (c) 2023 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 { Vec2 } from '../linear-algebra';
+
+/** raycast along x-axis and apply even-odd rule */
+export function pointInPolygon(point: Vec2, polygon: NumberArray, count: number): boolean {
+    const [x, y] = point;
+    let inside = false;
+
+    for (let i = 0, j = count - 1; i < count; j = i++) {
+        const xi = polygon[i * 2], yi = polygon[i * 2 + 1];
+        const xj = polygon[j * 2], yj = polygon[j * 2 + 1];
+
+        if (((yi > y) !== (yj > y)) && (x < (xj - xi) * (y - yi) / (yj - yi) + xi)) {
+            inside = !inside;
+        }
+    }
+    return inside;
+}

src/mol-math/geometry/primitives/box3d.ts

@@ -5,10 +5,11 @@
  * @author Alexander Rose <alexander.rose@weirdbyte.de>
  */
 
-import { Vec3, Mat4 } from '../../linear-algebra';
 import { PositionData } from '../common';
 import { OrderedSet } from '../../../mol-data/int';
 import { Sphere3D } from './sphere3d';
+import { Vec3 } from '../../linear-algebra/3d/vec3';
+import { Mat4 } from '../../linear-algebra/3d/mat4';
 
 interface Box3D { min: Vec3, max: Vec3 }
 
@@ -30,26 +31,48 @@ namespace Box3D {
         return copy(zero(), a);
     }
 
+    const tmpV = Vec3();
+
     /** Get box from sphere, uses extrema if available */
     export function fromSphere3D(out: Box3D, sphere: Sphere3D): Box3D {
         if (Sphere3D.hasExtrema(sphere) && sphere.extrema.length >= 14) { // 14 extrema with coarse boundary helper
             return fromVec3Array(out, sphere.extrema);
         }
-        const r = Vec3.create(sphere.radius, sphere.radius, sphere.radius);
-        Vec3.sub(out.min, sphere.center, r);
-        Vec3.add(out.max, sphere.center, r);
+        Vec3.set(tmpV, sphere.radius, sphere.radius, sphere.radius);
+        Vec3.sub(out.min, sphere.center, tmpV);
+        Vec3.add(out.max, sphere.center, tmpV);
         return out;
     }
 
-    /** Get box from sphere, uses extrema if available */
-    export function fromVec3Array(out: Box3D, array: Vec3[]): Box3D {
-        Box3D.setEmpty(out);
+    export function addVec3Array(out: Box3D, array: Vec3[]): Box3D {
         for (let i = 0, il = array.length; i < il; i++) {
-            Box3D.add(out, array[i]);
+            add(out, array[i]);
         }
         return out;
     }
 
+    export function fromVec3Array(out: Box3D, array: Vec3[]): Box3D {
+        setEmpty(out);
+        addVec3Array(out, array);
+        return out;
+    }
+
+    export function addSphere3D(out: Box3D, sphere: Sphere3D): Box3D {
+        if (Sphere3D.hasExtrema(sphere) && sphere.extrema.length >= 14) { // 14 extrema with coarse boundary helper
+            return addVec3Array(out, sphere.extrema);
+        }
+        add(out, Vec3.subScalar(tmpV, sphere.center, sphere.radius));
+        add(out, Vec3.addScalar(tmpV, sphere.center, sphere.radius));
+        return out;
+    }
+
+    export function intersectsSphere3D(box: Box3D, sphere: Sphere3D) {
+        // Find the point on the AABB closest to the sphere center.
+        Vec3.clamp(tmpV, sphere.center, box.min, box.max);
+        // If that point is inside the sphere, the AABB and sphere intersect.
+        return Vec3.squaredDistance(tmpV, sphere.center) <= (sphere.radius * sphere.radius);
+    }
+
     export function computeBounding(data: PositionData): Box3D {
         const min = Vec3.create(Number.MAX_VALUE, Number.MAX_VALUE, Number.MAX_VALUE);
         const max = Vec3.create(-Number.MAX_VALUE, -Number.MAX_VALUE, -Number.MAX_VALUE);
@@ -139,7 +162,16 @@ namespace Box3D {
         );
     }
 
-    // const tmpTransformV = Vec3();
+    export function containsSphere3D(box: Box3D, s: Sphere3D) {
+        const c = s.center;
+        const r = s.radius;
+        return (
+            c[0] - r < box.min[0] || c[0] + r > box.max[0] ||
+            c[1] - r < box.min[1] || c[1] + r > box.max[1] ||
+            c[2] - r < box.min[2] || c[2] + r > box.max[2]
+        ) ? false : true;
+    }
+
     export function nearestIntersectionWithRay(out: Vec3, box: Box3D, origin: Vec3, dir: Vec3): Vec3 {
         const [minX, minY, minZ] = box.min;
         const [maxX, maxY, maxZ] = box.max;

src/mol-math/geometry/primitives/frustum3d.ts

@@ -0,0 +1,99 @@
+/**
+ * Copyright (c) 2022-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-2022 three.js authors. MIT License
+ */
+
+import { Mat4 } from '../../linear-algebra/3d/mat4';
+import { Vec3 } from '../../linear-algebra/3d/vec3';
+import { Box3D } from './box3d';
+import { Plane3D } from './plane3d';
+import { Sphere3D } from './sphere3d';
+
+interface Frustum3D { 0: Plane3D, 1: Plane3D, 2: Plane3D, 3: Plane3D, 4: Plane3D, 5: Plane3D; length: 6; }
+
+function Frustum3D() {
+    return Frustum3D.create(Plane3D(), Plane3D(), Plane3D(), Plane3D(), Plane3D(), Plane3D());
+}
+
+namespace Frustum3D {
+    export const enum PlaneIndex {
+        Right = 0,
+        Left = 1,
+        Bottom = 2,
+        Top = 3,
+        Far = 4,
+        Near = 5,
+    };
+
+    export function create(right: Plane3D, left: Plane3D, bottom: Plane3D, top: Plane3D, far: Plane3D, near: Plane3D): Frustum3D {
+        return [right, left, bottom, top, far, near];
+    }
+
+    export function copy(out: Frustum3D, f: Frustum3D): Frustum3D {
+        for (let i = 0 as PlaneIndex; i < 6; ++i) Plane3D.copy(out[i], f[i]);
+        return out;
+    }
+
+    export function clone(f: Frustum3D): Frustum3D {
+        return copy(Frustum3D(), f);
+    }
+
+    export function fromProjectionMatrix(out: Frustum3D, m: Mat4) {
+        const a00 = m[0], a01 = m[1], a02 = m[2], a03 = m[3];
+        const a10 = m[4], a11 = m[5], a12 = m[6], a13 = m[7];
+        const a20 = m[8], a21 = m[9], a22 = m[10], a23 = m[11];
+        const a30 = m[12], a31 = m[13], a32 = m[14], a33 = m[15];
+
+        Plane3D.setUnnormalized(out[0], a03 - a00, a13 - a10, a23 - a20, a33 - a30);
+        Plane3D.setUnnormalized(out[1], a03 + a00, a13 + a10, a23 + a20, a33 + a30);
+        Plane3D.setUnnormalized(out[2], a03 + a01, a13 + a11, a23 + a21, a33 + a31);
+        Plane3D.setUnnormalized(out[3], a03 - a01, a13 - a11, a23 - a21, a33 - a31);
+        Plane3D.setUnnormalized(out[4], a03 - a02, a13 - a12, a23 - a22, a33 - a32);
+        Plane3D.setUnnormalized(out[5], a03 + a02, a13 + a12, a23 + a22, a33 + a32);
+
+        return out;
+    }
+
+    export function intersectsSphere3D(frustum: Frustum3D, sphere: Sphere3D) {
+        const center = sphere.center;
+        const negRadius = -sphere.radius;
+
+        for (let i = 0 as PlaneIndex; i < 6; ++i) {
+            const distance = Plane3D.distanceToPoint(frustum[i], center);
+            if (distance < negRadius) return false;
+        }
+        return true;
+    }
+
+    const boxTmpV = Vec3();
+    export function intersectsBox3D(frustum: Frustum3D, box: Box3D) {
+        for (let i = 0 as PlaneIndex; i < 6; ++i) {
+            const plane = frustum[i];
+
+            // corner at max distance
+            boxTmpV[0] = plane.normal[0] > 0 ? box.max[0] : box.min[0];
+            boxTmpV[1] = plane.normal[1] > 0 ? box.max[1] : box.min[1];
+            boxTmpV[2] = plane.normal[2] > 0 ? box.max[2] : box.min[2];
+
+            if (Plane3D.distanceToPoint(plane, boxTmpV) < 0) {
+                return false;
+            }
+        }
+        return true;
+    }
+
+    export function containsPoint(frustum: Frustum3D, point: Vec3) {
+        for (let i = 0 as PlaneIndex; i < 6; ++i) {
+            if (Plane3D.distanceToPoint(frustum[i], point) < 0) {
+                return false;
+            }
+        }
+        return true;
+    }
+}
+
+export { Frustum3D };

src/mol-math/geometry/primitives/plane3d.ts

@@ -0,0 +1,93 @@
+/**
+ * Copyright (c) 2022-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-2022 three.js authors. MIT License
+ */
+
+import { NumberArray } from '../../../mol-util/type-helpers';
+import { Vec3 } from '../../linear-algebra/3d/vec3';
+import { Sphere3D } from './sphere3d';
+
+interface Plane3D { normal: Vec3, constant: number }
+
+function Plane3D() {
+    return Plane3D.create(Vec3.create(1, 0, 0), 0);
+}
+
+namespace Plane3D {
+    export function create(normal: Vec3, constant: number): Plane3D { return { normal, constant }; }
+
+    export function copy(out: Plane3D, p: Plane3D): Plane3D {
+        Vec3.copy(out.normal, p.normal);
+        out.constant = p.constant;
+        return out;
+    }
+
+    export function clone(p: Plane3D): Plane3D {
+        return copy(Plane3D(), p);
+    }
+
+    export function normalize(out: Plane3D, p: Plane3D): Plane3D {
+        // Note: will lead to a divide by zero if the plane is invalid.
+        const inverseNormalLength = 1.0 / Vec3.magnitude(p.normal);
+        Vec3.scale(out.normal, p.normal, inverseNormalLength);
+        out.constant = p.constant * inverseNormalLength;
+        return out;
+    }
+
+    export function negate(out: Plane3D, p: Plane3D): Plane3D {
+        Vec3.negate(out.normal, p.normal);
+        out.constant = -p.constant;
+        return out;
+    }
+
+    export function toArray<T extends NumberArray>(p: Plane3D, out: T, offset: number) {
+        Vec3.toArray(p.normal, out, offset);
+        out[offset + 3] = p.constant;
+        return out;
+    }
+
+    export function fromArray(out: Plane3D, array: NumberArray, offset: number) {
+        Vec3.fromArray(out.normal, array, offset);
+        out.constant = array[offset + 3];
+        return out;
+    }
+
+    export function fromNormalAndCoplanarPoint(out: Plane3D, normal: Vec3, point: Vec3) {
+        Vec3.copy(out.normal, normal);
+        out.constant = -Vec3.dot(out.normal, point);
+        return out;
+    }
+
+    export function fromCoplanarPoints(out: Plane3D, a: Vec3, b: Vec3, c: Vec3) {
+        const normal = Vec3.triangleNormal(Vec3(), a, b, c);
+        fromNormalAndCoplanarPoint(out, normal, a);
+        return out;
+    }
+
+    const unnormTmpV = Vec3();
+    export function setUnnormalized(out: Plane3D, nx: number, ny: number, nz: number, constant: number) {
+        Vec3.set(unnormTmpV, nx, ny, nz);
+        const inverseNormalLength = 1.0 / Vec3.magnitude(unnormTmpV);
+        Vec3.scale(out.normal, unnormTmpV, inverseNormalLength);
+        out.constant = constant * inverseNormalLength;
+        return out;
+    }
+
+    export function distanceToPoint(plane: Plane3D, point: Vec3) {
+        return Vec3.dot(plane.normal, point) + plane.constant;
+    }
+
+    export function distanceToSpher3D(plane: Plane3D, sphere: Sphere3D) {
+        return distanceToPoint(plane, sphere.center) - sphere.radius;
+    }
+
+    export function projectPoint(out: Vec3, plane: Plane3D, point: Vec3) {
+        return Vec3.scaleAndAdd(out, out, plane.normal, -distanceToPoint(plane, point));
+    }
+}
+
+export { Plane3D };

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

@@ -1,5 +1,5 @@
 /**
- * Copyright (c) 2017-2021 mol* contributors, licensed under MIT, See LICENSE file for more info.
+ * Copyright (c) 2017-2022 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>
@@ -18,7 +18,7 @@
  */
 
 import { Mat4 } from './mat4';
-import { spline as _spline, quadraticBezier as _quadraticBezier, clamp } from '../../interpolate';
+import { spline as _spline, quadraticBezier as _quadraticBezier, clamp as _clamp } from '../../interpolate';
 import { NumberArray } from '../../../mol-util/type-helpers';
 import { Mat3 } from './mat3';
 import { Quat } from './quat';
@@ -246,6 +246,16 @@ namespace Vec3 {
         return out;
     }
 
+    /**
+     * Assumes min < max, componentwise
+     */
+    export function clamp(out: Vec3, a: Vec3, min: Vec3, max: Vec3) {
+        out[0] = Math.max(min[0], Math.min(max[0], a[0]));
+        out[1] = Math.max(min[1], Math.min(max[1], a[1]));
+        out[2] = Math.max(min[2], Math.min(max[2], a[2]));
+        return out;
+    }
+
     export function distance(a: Vec3, b: Vec3) {
         const x = b[0] - a[0],
             y = b[1] - a[1],
@@ -341,7 +351,7 @@ namespace Vec3 {
 
     const slerpRelVec = zero();
     export function slerp(out: Vec3, a: Vec3, b: Vec3, t: number) {
-        const d = clamp(dot(a, b), -1, 1);
+        const d = _clamp(dot(a, b), -1, 1);
         const theta = Math.acos(d) * t;
         scaleAndAdd(slerpRelVec, b, a, -d);
         normalize(slerpRelVec, slerpRelVec);
@@ -429,6 +439,14 @@ namespace Vec3 {
         return out;
     }
 
+    export function transformDirection(out: Vec3, a: Vec3, m: Mat4) {
+        const x = a[0], y = a[1], z = a[2];
+        out[0] = m[0] * x + m[4] * y + m[8] * z;
+        out[1] = m[1] * x + m[5] * y + m[9] * z;
+        out[2] = m[2] * x + m[6] * y + m[10] * z;
+        return normalize(out, out);
+    }
+
     /**
      * Like `transformMat4` but with offsets into arrays
      */
@@ -477,7 +495,7 @@ namespace Vec3 {
         const denominator = Math.sqrt(squaredMagnitude(a) * squaredMagnitude(b));
         if (denominator === 0) return Math.PI / 2;
         const theta = dot(a, b) / denominator;
-        return Math.acos(clamp(theta, -1, 1)); // clamp to avoid numerical problems
+        return Math.acos(_clamp(theta, -1, 1)); // clamp to avoid numerical problems
     }
 
     const tmp_dh_ab = zero();