Open any radiance-field codebase — bmild/nerf, yenchenlin/nerf-pytorch,
graphdeco-inria/gaussian-splatting, nerfstudio-project/gsplat — and the
first hundred lines are the same geometry you wrote in C++ for the gfx-*
skills: homogeneous coordinates, 4×4 camera transforms, a pinhole projection,
nearest neighbours. The difference is that there is no for loop over points.
A million Gaussians, a batch of cameras and every pixel of an image are
handled by one array expression, and the thing you reason about is not the
loop body but the shape of every array. In research Python, shapes are the
type system: get them right and the code is short, fast and differentiable;
get one axis wrong and it usually still runs — and computes something else.
This lesson is the NumPy half of that toolkit. py-torch-autograd carries
the same ideas to PyTorch (the broadcasting rules and einsum are identical;
torch.Tensor adds a device and a gradient).
1. Shapes are the contract
The conventions every codebase above shares:
- Points are rows.
Npoints in 3D are an(N, 3)array, not(3, N). A batch ofCcameras' world-to-camera matrices is(C, 4, 4), an image of rays is(H, W, 3), Gaussians' quaternions are(N, 4). - Batch dimensions lead, the "vector" dimension is last. A function that
works on one point set
(N, 3)should work unchanged on(B, N, 3)— write shapes as(..., N, 3)and index from the right (p[..., 2],axis=-1). - No per-element Python loops. A Python loop costs ~50–100 ns per iteration before doing any work; NumPy's inner loops are C. The gap is routinely 100× — and in PyTorch, a loop also builds a huge autograd graph.
// The C++ version (gfx-vectors-matrices): one point at a time.
for (std::size_t i = 0; i < pts.size(); ++i) out[i] = T * vec4{pts[i], 1.0f};# The NumPy version: the loop is inside the library.
out = to_homogeneous(pts) @ T.T # (N, 4) @ (4, 4) -> (N, 4)2. Broadcasting: align from the right, stretch the ones
Every elementwise operation (+ - * / **, comparisons, np.where, ufuncs)
broadcasts its operands to a common shape:
- Line the shapes up from the right; pad the shorter one with 1s on the left.
- In each position the sizes must be equal, or one of them must be 1.
- A size-1 axis is stretched (without copying — stride 0) to match.
The rule is mechanical, which is exactly why it bites: a shape that happens
to be compatible broadcasts silently. The classic case is a reduction that
drops its axis. np.linalg.norm(v, axis=1) of an (N, 3) array is (N,);
dividing (N, 3) by (N,) aligns that with the last axis — the
coordinates — and fails for N ≠ 3, or for N = 3 divides column j by the
norm of row j. keepdims=True keeps a size-1 axis where the reduction was,
so the result lines up with the rows:
import numpy as np
a = np.zeros((5, 1, 3))
b = np.zeros((4, 3))
print("(5,1,3) + (4,3) ->", (a + b).shape)
# The classic silent bug: normalising rows without keepdims.
v = np.array([[3.0, 4.0, 0.0], [0.0, 0.0, 2.0], [1.0, 0.0, 0.0]])
n = np.linalg.norm(v, axis=1) # shape (3,): one norm per row
wrong = v / n # (3,3) / (3,) divides COLUMN j by norm j
right = v / np.linalg.norm(v, axis=1, keepdims=True) # (3,3) / (3,1)
print("norms ", n)
print("wrong row 0 ", wrong[0], "length", np.linalg.norm(wrong[0]).round(3))
print("right row 0 ", right[0], "length", np.linalg.norm(right[0]).round(3))
try:
np.zeros((4, 3)) / np.zeros(4)
except ValueError as e:
print("(4,3) / (4,) ->", e)The tools for steering broadcasting are keepdims=True, x[..., None] (or
np.newaxis: insert a size-1 axis), slicing with a range (q[..., 3:] keeps
a length-1 axis where q[..., 3] drops it) and np.broadcast_to /
np.broadcast_shapes to check what you'll get. The 3DGS code is full of
r / r.norm(dim=1, keepdim=True) (PyTorch spells it keepdim) for exactly
this reason.
3. Homogeneous coordinates and batched transforms
A rigid transform (rotation R, translation t) and a perspective
projection both become a single 4×4 matrix product once a point gets a
fourth coordinate w = 1: [R t; 0 1] · [p; 1] = [R p + t; 1]. A direction
gets w = 0, so translation doesn't touch it. After a projective matrix,
divide by the resulting w.
Here is the trap: the maths writes column vectors (p' = T p), but the
code stores points as rows. For an (N, 4) stack the column-vector
product is ph @ T.T — transposing the matrix, not the points. Forget the
transpose and you apply Tᵀ: the rotation runs backwards and the translation
lands in w. np.einsum sidesteps the question by naming the axes:
"ij,nj->ni" reads "for each point n, contract T's column index j with the
point's coordinate j".
import numpy as np
rng = np.random.default_rng(0)
c, s = np.cos(np.pi / 2), np.sin(np.pi / 2)
T = np.array([[c, -s, 0, 1],
[s, c, 0, 2],
[0, 0, 1, 3],
[0, 0, 0, 1]]) # rotate 90° about z, then translate (1, 2, 3)
p = rng.normal(size=(1000, 3))
ph = np.concatenate([p, np.ones((len(p), 1))], axis=-1) # (N, 4)
loop = np.array([T @ x for x in ph]) # one 4x4 @ 4 per point
colvec = (T @ ph.T).T # (4,4) @ (4,N), transposed back
rowvec = ph @ T.T # rows: x' = x Tᵀ
einsum = np.einsum("ij,nj->ni", T, ph) # say exactly which axes contract
print("all agree:", np.allclose(loop, colvec) and np.allclose(loop, rowvec) and np.allclose(loop, einsum))
print("+x ->", (np.array([1.0, 0, 0, 1]) @ T.T)[:3].round(6))
print("ph @ T (forgot the transpose), +x ->", (np.array([1.0, 0, 0, 1]) @ T)[:3].round(6),
" w =", (np.array([1.0, 0, 0, 1]) @ T)[3])
# Batches: C camera transforms against the same N points, with '...'
Ts = np.stack([T, np.eye(4), T @ T]) # (3, 4, 4)
out = np.einsum("...ij,...nj->...ni", Ts, ph) # (3, N, 4): '...' broadcasts (3,) with ()
print("batched:", out.shape, "| matmul agrees:", np.allclose(out, ph @ Ts.swapaxes(-1, -2)))@ (np.matmul) treats the last two axes as matrices and broadcasts
everything before them, so (C, 4, 4) against (N, 4, 4) would fail but
(C, 1, 4, 4) against (N, 4, 4) gives (C, N, 4, 4). einsum's ...
broadcasts the same way. Use .swapaxes(-1, -2) (PyTorch: .transpose(-1, -2)
or .mT), never .T, on batched matrices: .T reverses all axes, so a
(C, 4, 4) stack becomes (4, 4, C).
Where codebases differ. The 3DGS reference stores its matrices
transposed: getWorld2View2 builds a column-vector world-to-camera matrix
and Camera keeps world_view_transform = ….transpose(0, 1), so its CUDA and
Python code multiply row vectors, p_h @ M. gsplat takes viewmats of
shape (C, 4, 4) as ordinary column-vector world-to-camera matrices. Before
reusing a matrix from another codebase, push one known point through it.
4. Projecting N points through K[R | t]
The pinhole camera in OpenCV's (and COLMAP's, and gsplat's) convention:
R, tare world-to-camera:x_cam = R X + t. The camera centre in world coordinates is−Rᵀ t.- Camera axes: x right, y down, z forward (into the scene). A point is in
front of the camera iff
z_cam > 0. K = [[fx, 0, cx], [0, fy, cy], [0, 0, 1]]maps camera coordinates to pixels:(u, v) = (fx·x/z + cx, fy·y/z + cy), i.e.(K x_cam)[:2] / z.
NeRF's Blender data uses OpenGL's axes instead (y up, z backwards, camera
looks down −z) and stores camera-to-world matrices — py-camera-rays is about
converting. Here we stay in OpenCV.
The divide by z is where points behind the camera go wrong. A point at
z = −1 divides to a perfectly finite pixel — the mirror image through the
camera centre, often right in the middle of the frame. Rasterisers clip such
points (3DGS's in_frustum culls p_view.z <= 0.2; gsplat has near_plane);
array code must mask them explicitly, and must not "fix" the sign with
abs(z):
import numpy as np
K = np.array([[500.0, 0, 320],
[ 0, 500, 240],
[ 0, 0, 1]]) # fx, fy, cx, cy for a 640x480 image
R = np.eye(3) # world-to-camera rotation
t = np.array([0.0, 0, 2]) # the world origin is 2 m in front
X = np.array([[0.0, 0, 0], # straight ahead
[0.5, -0.25, 0], # right and up (OpenCV y points down)
[0.0, 0, -3], # 1 m BEHIND the camera
[-0.5, 0.25, -4]]) # behind, mirrored
Xc = X @ R.T + t # (N, 3) camera coordinates
z = Xc[:, 2]
naive = (Xc @ K.T)[:, :2] / z[:, None] # what you get with no mask
valid = z > 1e-6
uv = np.where(valid[:, None], naive, np.nan)
for i in range(len(X)):
print(f"z = {z[i]:5.1f} naive = {naive[i].round(1)} masked = {uv[i].round(1)}")The last point is behind the camera, yet lands on exactly the same pixel as
the second. To avoid divide-by-zero warnings for points on the plane z = 0,
divide by a "safe" np.where(valid, z, 1) and mask afterwards. Batching over
cameras is the same code with ...: R (C, 3, 3), t (C, 3) and one cloud
X (N, 3) give x_cam (C, N, 3) with
np.einsum("...ij,...nj->...ni", R, X) + t[..., None, :] — note the
[..., None, :], without which t (C, 3) would try to broadcast against the
point axis.
5. Pairwise distances without loops
Nearest neighbours appear everywhere: 3DGS initialises each Gaussian's scale
from the mean distance to its three nearest neighbours (distCUDA2 from the
simple-knn extension), point-cloud losses (Chamfer) need all-pairs
distances, densification needs neighbours. With broadcasting, all-pairs is one
line: insert axes so A becomes (N, 1, D) and B becomes (1, M, D);
their difference is (N, M, D), and summing squares over the last axis gives
the (N, M) matrix.
The other classic formulation expands the square,
‖a − b‖² = ‖a‖² + ‖b‖² − 2 a·b, so the (N, M) matrix comes from one matrix
product A @ B.T — O(N·M) memory instead of O(N·M·D), and BLAS-fast. It is
what torch.cdist does by default for large inputs
(compute_mode="use_mm_for_euclid_dist_if_necessary"). But it subtracts two
large numbers to get a small one, and in float32 that is catastrophic
cancellation: far from the origin ‖a‖² is huge, its rounding error is
bigger than the distance you want, and results go negative.
import numpy as np
rng = np.random.default_rng(0)
for offset in (0.0, 1e2, 1e3, 1e4):
A = (offset + rng.normal(scale=0.1, size=(200, 3))).astype(np.float32)
B = (offset + rng.normal(scale=0.1, size=(300, 3))).astype(np.float32)
exact = ((A.astype(np.float64)[:, None] - B.astype(np.float64)[None]) ** 2).sum(-1)
diff = A[:, None, :] - B[None, :, :] # (N, M, D): subtract first
d_diff = np.einsum("nmd,nmd->nm", diff, diff)
d_gram = (A * A).sum(1)[:, None] + (B * B).sum(1)[None, :] - 2 * A @ B.T # (N, M) via matmul
print(f"offset {offset:7.0f}: max error diff {np.abs(d_diff - exact).max():.1e}"
f" expansion {np.abs(d_gram - exact).max():.1e} negatives {int((d_gram < 0).sum())}")
print("typical squared distance:", round(float(np.median(exact)), 3))Real scenes are exactly this case: a COLMAP reconstruction in metres with the
origin wherever the first camera was, points thousands of units out. The
fixes, in order of preference: subtract first (exact, more memory — chunk over
N if needed); or centre a copy of the points (A - c, B - c with a
shared c) and use the expansion in float64; and always np.maximum(d2, 0).
6. Views and copies
A NumPy array is a pointer, a dtype, a shape and strides (bytes to step per axis). Many operations only make a new header over the same memory — a view — and writing through a view writes the original:
| view (shares memory) | copy |
|---|---|
basic slicing P[:, :2], P[::2], P[..., None] |
fancy indexing P[:, [0, 1]], P[idx] |
.T, .swapaxes, np.broadcast_to (stride 0, read-only) |
boolean masks P[valid] |
.reshape / .view when the strides allow it |
.reshape when they don't, .astype, .copy() |
np.asarray(x) of an array (no copy at all) |
np.array(x) (copies by default) |
import numpy as np
P = np.arange(12.0).reshape(4, 3) # 4 points, xyz rows
xyz_view = P[:, :2] # basic slice: a view
fancy = P[:, [0, 1]] # fancy index: a copy
mask = P[P[:, 2] > 4] # boolean mask: a copy
print("slice shares memory:", np.shares_memory(P, xyz_view),
"| fancy:", np.shares_memory(P, fancy), "| mask:", np.shares_memory(P, mask))
xyz_view *= 10 # writes through into P
fancy[:] = -1 # writes into the copy only
print(P[0], "<- P[0] after scaling the view")
PT = P.T # transpose: a view with swapped strides
print("P strides", P.strides, "| P.T strides", PT.strides, "| P.T contiguous:", PT.flags.c_contiguous)
flat = PT.reshape(-1) # can't be a view of a non-contiguous array: copies
print("reshape of P.T shares memory:", np.shares_memory(PT, flat))
def center_in_place(pts): # looks harmless...
pts -= pts.mean(axis=0)
return pts
Q = np.array([[1.0, 1, 1], [3, 3, 3]])
center_in_place(Q)
print("caller's Q after center_in_place:", Q.tolist())Two rules follow. Never mutate an argument you didn't create: x -= c
and x[...] = … modify the caller's array (and np.asarray won't protect
you, it returns the same object); write x = x - c. And when you want an
in-place update for speed, make sure the target isn't a broadcast or a
shared view. PyTorch has the same split (view fails on non-contiguous
tensors where reshape copies; expand is stride-0 like broadcast_to)
plus a third hazard, in-place ops on tensors autograd saved — that's
py-torch-autograd.
7. float32 vs float64
NumPy defaults to float64; PyTorch, GPUs and every rendering codebase to
float32 (half the memory and bandwidth, and GPU float64 throughput is a
small fraction of float32). float32 has a 24-bit significand: about 7
significant digits, a relative step of 1.19e-7, and an absolute step of
~0.001 at 10⁴. Two consequences: keep coordinates centred and scaled (NeRF
pipelines normalise scenes into a unit box for exactly this reason), and
don't let float64 sneak in by accident — it doubles memory and, when the
result goes back into PyTorch, raises a dtype mismatch.
NumPy 2 promotes by NEP 50 rules: Python scalars (2.0, 1, np.nan)
are "weak" and take the array's dtype; NumPy scalars and arrays are not.
import numpy as np
x = np.ones(3, dtype=np.float32)
print("f32 * 2.0 ->", (x * 2.0).dtype) # Python float: weak, stays float32
print("f32 * np.float64(2) ->", (x * np.float64(2)).dtype) # NumPy float64 scalar: promotes
print("f32 + np.ones(3) ->", (x + np.ones(3)).dtype) # np.ones defaults to float64
print("f32 + 1 (int) ->", (x + 1).dtype)
print("int arange / 2 ->", (np.arange(4) / 2).dtype, " // ->", (np.arange(4) // 2).dtype)
# float32 has a 24-bit significand: ~7 decimal digits.
print("eps32 =", np.finfo(np.float32).eps, " eps64 =", np.finfo(np.float64).eps)
far = np.float32(16_777_216.0) # 2**24
print("2**24 + 1 in float32:", far + np.float32(1), "(the +1 is lost)")
print("spacing of float32 at 1e4:", np.spacing(np.float32(1e4)))The usual culprits are np.ones(shape) / np.zeros(shape) / np.eye(n)
without dtype=, np.linspace, constants wrapped in np.float64(...), and loading
JSON (Python floats → float64 arrays). Write helpers
that create arrays with dtype=p.dtype and test them with float32 inputs.
Where you need float64 is in checks: torch.autograd.gradcheck and
reference comparisons run in float64 so that rounding doesn't drown the error
you're looking for.
In an interview
"Vectorise this loop" is the standard Python question for a graphics role:
expect to project N points through K[R|t] with a mask, compute all-pairs
distances, or apply a batch of transforms, and to say the output shape of
each line out loud. Good things to volunteer: column vectors in the maths are
rows in the array, so it's X @ R.T + t; keepdims and [..., None] to
steer broadcasting; masking z ≤ 0 before (not after) trusting a pixel; the
memory/precision trade-off between (N, M, D) differences and the
‖a‖² + ‖b‖² − 2a·b expansion; and float32 by default with float64 only for
verification.