"""Pure NumPy matrix utilities — replaces GLM dependency.
All matrices use NumPy's native row-major layout. When sending to GLSL:
- Matrices must be transposed (GLSL expects column-major)
- This is handled at GPU boundary points only
This eliminates the entire class of row-major vs column-major bugs.
"""
from __future__ import annotations
import logging
import numpy as np
log = logging.getLogger(__name__)
__all__ = [
"identity",
"perspective",
"look_at",
"translate",
"rotate",
"scale",
"orthographic",
"quat_to_mat4",
]
[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. glm.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 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