"""Pure NumPy matrix utilities.
All matrices use NumPy's native row-major layout. When sending to GLSL,
matrices must be transposed (GLSL expects column-major); this transpose
happens only at GPU boundary points, keeping the rest of the engine free
of row-major vs column-major ambiguity.
"""
import logging
import numpy as np
log = logging.getLogger(__name__)
__all__ = [
"identity",
"perspective",
"oblique_near_plane",
"look_at",
"translate",
"rotate",
"scale",
"orthographic",
"quat_to_mat4",
"mat4_from_trs",
"batch_mat4_from_trs",
"mat4_to_bytes",
"halton_jitter",
"apply_jitter",
]
[docs]
def quat_to_mat4(q, dtype: np.dtype = np.float32) -> np.ndarray:
"""Convert a quaternion to a 4x4 rotation matrix.
Args:
q: Quaternion: any object with .w, .x, .y, .z attributes
(e.g. Quat) or a 4-element array [w, x, y, z].
dtype: NumPy data type (default float32)
Returns:
4x4 rotation matrix as numpy array
"""
if hasattr(q, "w"):
w, x, y, z = float(q.w), float(q.x), float(q.y), float(q.z)
else:
w, x, y, z = float(q[0]), float(q[1]), float(q[2]), float(q[3])
xx, yy, zz = x * x, y * y, z * z
xy, xz, yz = x * y, x * z, y * z
wx, wy, wz = w * x, w * y, w * z
result = np.eye(4, dtype=dtype)
result[0, 0] = 1.0 - 2.0 * (yy + zz)
result[0, 1] = 2.0 * (xy - wz)
result[0, 2] = 2.0 * (xz + wy)
result[1, 0] = 2.0 * (xy + wz)
result[1, 1] = 1.0 - 2.0 * (xx + zz)
result[1, 2] = 2.0 * (yz - wx)
result[2, 0] = 2.0 * (xz - wy)
result[2, 1] = 2.0 * (yz + wx)
result[2, 2] = 1.0 - 2.0 * (xx + yy)
return result
[docs]
def identity(dtype: np.dtype = np.float32) -> np.ndarray:
"""Create 4x4 identity matrix.
Args:
dtype: NumPy data type (default float32)
Returns:
4x4 identity matrix as numpy array
"""
return np.eye(4, dtype=dtype)
[docs]
def perspective(
fov: float,
aspect: float,
near: float,
far: float,
dtype: np.dtype = np.float32,
) -> np.ndarray:
"""Create perspective projection matrix.
Args:
fov: Field of view in radians
aspect: Aspect ratio (width / height)
near: Near clipping plane
far: Far clipping plane
dtype: NumPy data type (default float32)
Returns:
4x4 perspective projection matrix
Note:
Does NOT include Y-flip for Vulkan. Caller must do:
proj[1, 1] *= -1 # Flip Y-axis for Vulkan
"""
# Guard against degenerate inputs
if aspect < 1e-6:
aspect = 1.0
if abs(far - near) < 1e-10:
far = near + 1.0
fov = max(np.radians(1.0), min(fov, np.radians(179.0)))
f = 1.0 / np.tan(fov / 2.0)
result = np.zeros((4, 4), dtype=dtype)
result[0, 0] = f / aspect
result[1, 1] = f
result[2, 2] = (far + near) / (near - far)
result[2, 3] = (2.0 * far * near) / (near - far)
result[3, 2] = -1.0
return result
[docs]
def oblique_near_plane(
proj: np.ndarray,
view: np.ndarray,
plane: tuple[float, float, float, float] | np.ndarray,
dtype: np.dtype = np.float32,
) -> np.ndarray:
"""Replace *proj*'s near clip plane with an arbitrary world-space plane.
The oblique-frustum technique (Lengyel, "Oblique View Frustum Depth
Projection and Clipping"): the projection's z row is rewritten so the GPU's
fixed ``0 <= z_clip <= w_clip`` near clip coincides with *plane*, scaled so
the new far plane still touches the corner of the original frustum that
lies deepest on the kept side. Used by planar reflections: the mirrored
camera must not render geometry behind the reflection plane.
Args:
proj: 4x4 projection matrix (row-major, any convention whose hardware
clip volume is ``0 <= z_clip <= w_clip``; both backends qualify).
view: 4x4 view matrix of the camera the projection belongs to.
plane: World-space plane ``(nx, ny, nz, d)`` with the equation
``n . x + d = 0``; geometry on the ``n . x + d > 0`` side is kept.
Returns:
A new 4x4 projection matrix (input is not modified). Depth values are
redistributed (inherent to the technique); depth testing within the
pass stays consistent. Returns *proj* unchanged when the plane does not
face any part of the view frustum (degenerate input).
"""
p = np.array(proj, dtype=np.float64)
c_cam = np.linalg.inv(np.asarray(view, dtype=np.float64)).T @ np.asarray(plane, dtype=np.float64)
# Far-plane corner of the original frustum deepest on the plane's kept
# side, found in camera space by inverse-projecting all four clip corners
# (robust to Y-flips and axis sign conventions baked into ``proj``).
inv_p = np.linalg.inv(p)
corners = (
inv_p
[docs]
@ np.array(
[[-1.0, -1.0, 1.0, 1.0], [1.0, -1.0, 1.0, 1.0], [-1.0, 1.0, 1.0, 1.0], [1.0, 1.0, 1.0, 1.0]],
dtype=np.float64,
).T
) # (4, 4): one camera-space corner per column
dots = c_cam @ corners
best = int(np.argmax(dots))
if dots[best] <= 1e-12:
log.warning("oblique_near_plane: plane faces away from the whole frustum; projection unchanged")
return np.array(proj, dtype=dtype)
q = corners[:, best]
p[2] = c_cam * (float(p[3] @ q) / float(dots[best]))
return p.astype(dtype)
def look_at(
eye: np.ndarray | tuple[float, float, float],
center: np.ndarray | tuple[float, float, float],
up: np.ndarray | tuple[float, float, float],
dtype: np.dtype = np.float32,
) -> np.ndarray:
"""Create view matrix using look-at vectors.
Args:
eye: Camera position
center: Point to look at
up: Up vector (should be normalized)
dtype: NumPy data type (default float32)
Returns:
4x4 view matrix
"""
eye = np.asarray(eye, dtype=dtype)
center = np.asarray(center, dtype=dtype)
up = np.asarray(up, dtype=dtype)
# Compute orthonormal basis
forward = center - eye
fwd_len = np.linalg.norm(forward)
if fwd_len < 1e-10:
return np.eye(4, dtype=dtype)
forward = forward / fwd_len
right = np.cross(forward, up)
right_len = np.linalg.norm(right)
if right_len < 1e-6:
# forward is nearly parallel to up: pick a fallback up vector
# Use the axis most perpendicular to forward for a stable fallback
abs_fwd = np.abs(forward)
if abs_fwd[0] <= abs_fwd[1] and abs_fwd[0] <= abs_fwd[2]:
fallback = np.array([1, 0, 0], dtype=dtype)
elif abs_fwd[1] <= abs_fwd[2]:
fallback = np.array([0, 1, 0], dtype=dtype)
else:
fallback = np.array([0, 0, 1], dtype=dtype)
right = np.cross(forward, fallback)
right_len = np.linalg.norm(right)
if right_len < 1e-10:
return np.eye(4, dtype=dtype)
right = right / right_len
up_new = np.cross(right, forward)
# Verify right-handed basis: det(right, up_new, -forward) must be positive.
# Equivalent to checking that the triple product right . (up_new x -forward) > 0.
det = np.dot(right, np.cross(up_new, -forward))
if det < 0:
right = -right
up_new = -up_new
# Build view matrix
result = np.eye(4, dtype=dtype)
result[0, :3] = right
result[1, :3] = up_new
result[2, :3] = -forward
result[0, 3] = -np.dot(right, eye)
result[1, 3] = -np.dot(up_new, eye)
result[2, 3] = np.dot(forward, eye)
return result
[docs]
def translate(
t: np.ndarray | tuple[float, float, float],
dtype: np.dtype = np.float32,
) -> np.ndarray:
"""Create translation matrix.
Args:
t: Translation vector (x, y, z)
dtype: NumPy data type (default float32)
Returns:
4x4 translation matrix
"""
t = np.asarray(t, dtype=dtype)
result = np.eye(4, dtype=dtype)
result[0, 3] = t[0]
result[1, 3] = t[1]
result[2, 3] = t[2]
return result
[docs]
def rotate(
axis: np.ndarray | tuple[float, float, float],
angle: float,
dtype: np.dtype = np.float32,
) -> np.ndarray:
"""Create rotation matrix using axis-angle representation.
Args:
axis: Rotation axis (should be normalized)
angle: Rotation angle in radians
dtype: NumPy data type (default float32)
Returns:
4x4 rotation matrix
Uses Rodrigues' rotation formula.
"""
axis = np.asarray(axis, dtype=dtype)
axis = axis / np.linalg.norm(axis) # Normalize
cos_a = np.cos(angle)
sin_a = np.sin(angle)
# Rodrigues' formula: R = I + sin(θ)K + (1-cos(θ))K²
# where K is the skew-symmetric cross-product matrix of axis
# Skew-symmetric matrix
K = np.array(
[
[0, -axis[2], axis[1]],
[axis[2], 0, -axis[0]],
[-axis[1], axis[0], 0],
],
dtype=dtype,
)
# Compute rotation matrix
R = np.eye(3, dtype=dtype) + sin_a * K + (1 - cos_a) * (K @ K)
# Embed in 4x4 matrix
result = np.eye(4, dtype=dtype)
result[:3, :3] = R
return result
[docs]
def scale(
s: np.ndarray | tuple[float, float, float] | float,
dtype: np.dtype = np.float32,
) -> np.ndarray:
"""Create scale matrix.
Args:
s: Scale factors (x, y, z) or uniform scale factor
dtype: NumPy data type (default float32)
Returns:
4x4 scale matrix
"""
if isinstance(s, int | float):
# Uniform scale
scale_vec = np.array([s, s, s], dtype=dtype)
else:
scale_vec = np.asarray(s, dtype=dtype)
result = np.eye(4, dtype=dtype)
result[0, 0] = scale_vec[0]
result[1, 1] = scale_vec[1]
result[2, 2] = scale_vec[2]
return result
[docs]
def orthographic(
left: float,
right: float,
bottom: float,
top: float,
near: float,
far: float,
dtype: np.dtype = np.float32,
) -> np.ndarray:
"""Create orthographic projection matrix.
Args:
left: Left plane
right: Right plane
bottom: Bottom plane
top: Top plane
near: Near plane
far: Far plane
dtype: NumPy data type (default float32)
Returns:
4x4 orthographic projection matrix
Note:
For Vulkan, you may need to flip Y: proj[1, 1] *= -1
"""
# Guard against degenerate inputs
if abs(right - left) < 1e-10:
right = left + 1.0
if abs(top - bottom) < 1e-10:
top = bottom + 1.0
if abs(far - near) < 1e-10:
far = near + 1.0
result = np.zeros((4, 4), dtype=dtype)
result[0, 0] = 2.0 / (right - left)
result[1, 1] = 2.0 / (top - bottom)
result[2, 2] = -2.0 / (far - near)
result[0, 3] = -(right + left) / (right - left)
result[1, 3] = -(top + bottom) / (top - bottom)
result[2, 3] = -(far + near) / (far - near)
result[3, 3] = 1.0
return result
[docs]
def mat4_from_trs(
pos: tuple[float, float, float] | np.ndarray,
rot,
scl: tuple[float, float, float] | np.ndarray,
) -> np.ndarray:
"""Build model matrix from position, rotation quaternion, and scale.
Args:
pos: Position (x, y, z)
rot: Rotation quaternion: Quat or any object with .w/.x/.y/.z
scl: Scale (x, y, z)
Returns:
4x4 model matrix as numpy array (Translate * Rotate * Scale)
"""
return translate(pos) @ quat_to_mat4(rot) @ scale(scl)
[docs]
def batch_mat4_from_trs(
positions: np.ndarray,
rotations: np.ndarray,
scales: np.ndarray,
) -> np.ndarray:
"""Build N model matrices from arrays of positions, quaternions, and scales.
Args:
positions: (N, 3) float32 positions
rotations: (N, 4) float32 quaternions [w, x, y, z]
scales: (N, 3) float32 scale factors
Returns:
(N, 4, 4) float32 model matrices (Translate * Rotate * Scale)
"""
n = positions.shape[0]
w, x, y, z = rotations[:, 0], rotations[:, 1], rotations[:, 2], rotations[:, 3]
xx, yy, zz = x * x, y * y, z * z
xy, xz, yz = x * y, x * z, y * z
wx, wy, wz = w * x, w * y, w * z
sx, sy, sz = scales[:, 0], scales[:, 1], scales[:, 2]
out = np.zeros((n, 4, 4), dtype=np.float32)
out[:, 0, 0] = (1.0 - 2.0 * (yy + zz)) * sx
out[:, 0, 1] = (2.0 * (xy - wz)) * sy
out[:, 0, 2] = (2.0 * (xz + wy)) * sz
out[:, 1, 0] = (2.0 * (xy + wz)) * sx
out[:, 1, 1] = (1.0 - 2.0 * (xx + zz)) * sy
out[:, 1, 2] = (2.0 * (yz - wx)) * sz
out[:, 2, 0] = (2.0 * (xz - wy)) * sx
out[:, 2, 1] = (2.0 * (yz + wx)) * sy
out[:, 2, 2] = (1.0 - 2.0 * (xx + yy)) * sz
out[:, 0, 3] = positions[:, 0]
out[:, 1, 3] = positions[:, 1]
out[:, 2, 3] = positions[:, 2]
out[:, 3, 3] = 1.0
return out
[docs]
def mat4_to_bytes(m: np.ndarray) -> bytes:
"""Convert mat4 to bytes for GPU upload (64 bytes, row-major float32)."""
return np.ascontiguousarray(m, dtype=np.float32).tobytes()
def _halton(index: int, base: int) -> float:
"""Radical-inverse of ``index`` in ``base`` (Halton low-discrepancy sequence)."""
f, r, i = 1.0, 0.0, index
while i > 0:
f /= base
r += f * (i % base)
i //= base
return r
[docs]
def halton_jitter(index: int, base_x: int = 2, base_y: int = 3) -> tuple[float, float]:
"""Sub-pixel jitter offset in [-0.5, 0.5] from the Halton(base_x, base_y) sequence.
Used for TAA: each frame samples a different sub-pixel location so the temporal
accumulation reconstructs detail below one pixel. ``index`` is the frame counter
(callers typically cycle it modulo the history length, e.g. 8). Mirrors the web
runtime's ``TAAPass.haltonJitter`` exactly so desktop and web jitter identically.
Returns:
``(jx, jy)`` offsets in pixels, each in ``[-0.5, 0.5]``.
"""
return _halton(index + 1, base_x) - 0.5, _halton(index + 1, base_y) - 0.5
[docs]
def apply_jitter(proj: np.ndarray, jx: float, jy: float, width: float, height: float) -> np.ndarray:
"""Return a copy of ``proj`` translated by a sub-pixel jitter in clip space.
``proj`` is a row-major Vulkan projection matrix (post Y-flip). The pixel-space
jitter ``(jx, jy)`` is converted to NDC (``2 / dim`` per pixel) and added to the
clip-x/clip-y rows of the third column (``proj[0, 2]`` / ``proj[1, 2]``), the
same slots the renderer transposes to the GPU's column-2 rows 0/1. The original
matrix is left untouched so callers keep an unjittered copy for culling and
motion-vector reprojection.
"""
out = proj.astype(np.float32, copy=True)
if width > 0:
out[0, 2] += jx * (2.0 / width)
if height > 0:
out[1, 2] += jy * (2.0 / height)
return out