"""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__)
_warned_zero_axis = False
def _warn_zero_axis() -> None:
"""Report a zero-length rotation axis once, then stay quiet.
The cause is almost always a cross product of two directions that turned out
parallel, which happens for a frame at a time in ordinary play, so this must
not fill the log.
"""
global _warned_zero_axis
if _warned_zero_axis:
return
_warned_zero_axis = True
log.warning(
"rotate() was given a zero-length axis and returned identity. This usually means "
"an axis came from a cross product of two parallel directions. "
"Further occurrences are not logged."
)
__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 a reverse-Z 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:
Depth maps to ``z_ndc`` in **[0, 1]** with the **near plane at 1 and the
far plane at 0**, matching :func:`orthographic`. Two things follow, and
the engine's pipelines depend on both: the depth attachment clears to
**0.0**, and every camera-space depth test compares ``GREATER`` or
``GREATER_OR_EQUAL``.
The reversal is what makes a float depth buffer usable over a wide
range. A float's exponent crowds its precision near zero, and a
perspective divide crowds depth near the near plane; forward-mapping
puts both concentrations in the same place and wastes them, while
reversing lines the float's fine steps up with the far field where the
projection's own steps are coarsest.
The older OpenGL ``[-1, 1]`` form is not merely a different spelling
here. The hardware clip volume is ``0 <= z_clip <= w_clip`` on both
backends, so that form put the effective near plane at
``2 f n / (f + n)``: a camera asking for ``near = 0.1`` really clipped
at ``0.1998``, and nothing said so.
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] = near / (far - near)
result[2, 3] = (near * far) / (far - near)
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 one of
the GPU's two fixed depth clip planes coincides with *plane*, scaled so the
opposite 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.
The engine is reverse-Z, so the near clip is ``z_clip <= w_clip`` rather
than ``z_clip >= 0`` and the substituted row is ``w_row - alpha * plane``
rather than Lengyel's ``alpha * plane``. Both forms clip identically; only
the reversed one leaves depth running near-to-far the way the pipelines'
``GREATER`` comparison expects, so the textbook form would silently invert
depth testing inside every planar-reflection pass.
Args:
proj: 4x4 reverse-Z projection matrix (row-major): hardware clip volume
``0 <= z_clip <= w_clip`` with the near plane at ``z_clip = w_clip``.
Both backends and both of this module's projections 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``). Under
# reverse-Z the far plane is ``z_ndc = 0``.
inv_p = np.linalg.inv(p)
corners = (
inv_p
[docs]
@ np.array(
[[-1.0, -1.0, 0.0, 1.0], [1.0, -1.0, 0.0, 1.0], [-1.0, 1.0, 0.0, 1.0], [1.0, 1.0, 0.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] = p[3] - 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)
# 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)
norm = float(np.linalg.norm(axis))
# A zero-length axis names no rotation. Dividing by its length would put NaN
# in every element, and NaN spreads: it would reach a transform, a bounding
# box and a GPU buffer before anything looked wrong. Identity is what an
# axis of no length means, and it is what the quaternion path returns.
if norm < 1e-10:
_warn_zero_axis()
return np.eye(4, dtype=dtype)
axis = axis / norm
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:
Depth maps to ``z_ndc`` in **[0, 1]**, the clip volume both backends
(Vulkan, WebGPU) actually enforce, with the **near plane at 1 and the
far plane at 0** so an orthographic camera shares the reverse-Z
convention :func:`perspective` uses and the same clear value and depth
comparison serve both.
An OpenGL-style [-1, 1] ortho would put half the depth range at
negative ``z_clip`` and the hardware would clip away half the scene:
unlike the perspective case there is no non-linearity to hide it.
Reversal buys an orthographic camera no precision, since its depth is
linear in distance. It is here for uniformity: one clear value and one
compare op across both projections is what keeps the two from needing
separate pipelines.
Does NOT include the Y-flip for Vulkan. Caller must do:
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] = 1.0 / (far - near)
result[0, 3] = -(right + left) / (right - left)
result[1, 3] = -(top + bottom) / (top - bottom)
result[2, 3] = far / (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.
Under an orthographic projection ``w`` is a constant 1 rather than ``-z_eye``,
so a third-column offset would shear with depth instead of translating. The
offset moves to the translation column, negated to land on the same NDC shift
the perspective path produces (``+p02`` divided by ``w = -z_eye``); the TAA
resolve inverts the same signed jitter for both.
"""
out = proj.astype(np.float32, copy=True)
orthographic_w = abs(float(proj[3, 2])) < 1e-9
col = 3 if orthographic_w else 2
sign = -1.0 if orthographic_w else 1.0
if width > 0:
out[0, col] += sign * jx * (2.0 / width)
if height > 0:
out[1, col] += sign * jy * (2.0 / height)
return out