Cameras and projection

Every renderer — rasterizer, ray tracer, NeRF, 3D Gaussian splatting — starts with the same question: which pixel does this 3D point land on, and which 3D ray does this pixel see? The answer is the pinhole camera: a rigid transform into the camera's frame, a divide by depth, and a scale-and-shift into pixels. The maths fits on an index card. The bugs come from conventions: which way the axes point, which direction a matrix maps, where a pixel's centre is. In radiance-field work you load poses from COLMAP, Blender and phone captures in the same week, and a single flipped axis turns a scene inside out.

1. The pinhole model

Put the centre of projection at the origin, looking down +z. A point at (x, y, z) and its image on a plane at distance f form similar triangles: the image point is (f·x/z, f·y/z). That divide by z is perspective: twice as far away, half as big.

Three coordinate systems are involved, and every camera bug is a mix-up between two of them:

2. Extrinsics: world-to-camera vs camera-to-world

A camera's pose is a rotation and a translation. It can be written in two directions, and both are in common use:

• world-to-camera (w2c, the extrinsics / view matrix): x_c = R·x_w + t. COLMAP's images.txt, OpenCV's solvePnP and the OpenGL view matrix are w2c.
• camera-to-world (c2w, the pose): x_w = R_c2w·x_c + C. Its columns are the camera's x, y, z axes expressed in world space, and C is the camera centre. NeRF's transforms.json, Blender and most SLAM trajectories are c2w.

They are inverses. For a rotation the inverse is the transpose, so R_c2w = Rᵀ and the centre is the point with camera coordinates zero: R·C + t = 0 ⇒ C = −Rᵀt. Note what t is not: it is not the camera position (it's where the world origin lands in camera coordinates). Treating COLMAP's tvec as the camera position is a rite of passage.

#include <array>
#include <print>

struct V3 { float x, y, z; };
using M3 = std::array<std::array<float, 3>, 3>;   // m[row][col]

V3 mul(const M3& m, V3 v) {
    return {m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z,
            m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z,
            m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z};
}
M3 transpose(const M3& m) {
    M3 t{};
    for (int r = 0; r < 3; ++r)
        for (int c = 0; c < 3; ++c) t[r][c] = m[c][r];
    return t;
}

int main() {
    // A camera at C = (1, 2, 3), turned 90° about y: its forward axis (+z,
    // third column of R_c2w) points along world +x.
    M3 R_c2w{{{0, 0, 1}, {0, 1, 0}, {-1, 0, 0}}};
    V3 C{1, 2, 3};

    M3 R = transpose(R_c2w);                        // world-to-camera rotation
    V3 RC = mul(R, C);
    V3 t{-RC.x, -RC.y, -RC.z};                      // t = -R C

    V3 p{4, 2, 3};                                  // 3 units in front of the camera
    V3 c = mul(R, p);
    c = {c.x + t.x, c.y + t.y, c.z + t.z};
    std::println("t = ({}, {}, {})   p in camera = ({}, {}, {})", t.x, t.y, t.z, c.x, c.y, c.z);

    V3 back = mul(transpose(R), t);                 // centre = -R^T t
    std::println("centre = ({}, {}, {})", -back.x, -back.y, -back.z);
}

t = (3, −2, −1) looks nothing like the position (1, 2, 3) — yet −Rᵀt recovers it. Name your variables R_w2c / R_c2w (or world_from_cam, as Eigen-style code does) and the type system of your head will catch most of these bugs.

3. Intrinsics: focal length and principal point

After the divide, (x/z, y/z) is in units of "tan of the angle off-axis". The intrinsic matrix turns that into pixels:

    ⎡ fx  0  cx ⎤        u = fx · x/z + cx
K = ⎢  0 fy  cy ⎥        v = fy · y/z + cy
    ⎣  0  0   1 ⎦

fx, fy are the focal length measured in pixels (equal for square pixels), (cx, cy) is the principal point where the optical axis hits the image — about the image centre, but calibrated cameras put it a few pixels off. Real lenses add a skew term (≈ 0) and radial/tangential distortion; COLMAP's OPENCV camera model has four distortion parameters, and datasets are usually undistorted first so the pure pinhole model applies.

Graphics APIs talk about a field of view instead. Half the image width subtends half the angle: tan(fovx/2) = (W/2) / fx, so

fx = W / (2 · tan(fovx / 2))        fovx = 2 · atan(W / (2 fx))

#include <cmath>
#include <numbers>
#include <print>

float focal_from_fov(int width, float fov_degrees) {
    float fov = fov_degrees * std::numbers::pi_v<float> / 180.0f;   // tan wants radians
    return width / (2.0f * std::tan(fov / 2.0f));                     // and the HALF angle
}

int main() {
    for (float deg : {30.0f, 60.0f, 90.0f, 120.0f})
        std::println("fov {:>5}°  ->  fx = {:7.1f} px (W = 800)", deg, focal_from_fov(800, deg));
    std::println("wrong (degrees to tan): {:.1f}", 800 / (2.0f * std::tan(90.0f / 2)));
    std::println("wrong (full angle)    : {:.1f}", 800 / (2.0f * std::tan(std::numbers::pi_v<float> / 2)));
}

Both wrong versions produce a number, and the image even looks like a scene — just at the wrong zoom. Unit-test fx(W, 90°) == W/2. Also watch whether a file gives fovx or fovy (OpenGL's gluPerspective takes the vertical one; NeRF's Blender transforms.json stores camera_angle_x).

4. Axis conventions: OpenCV vs OpenGL

Everyone agrees x points right. They disagree about y and z:

Both camera frames sit at the same place looking at the same thing; only the labels of two axes differ. So converting a camera-to-world matrix flips its y and z columns:

c2w_opencv = c2w_opengl · diag(1, −1, −1, 1)       (and w2c: flip rows)

This is exactly what the 3DGS and nerfstudio loaders do when reading Blender scenes (c2w[:3, 1:3] *= -1). Flip only one axis and the rotation becomes a reflection (determinant −1): the scene renders mirrored, and quaternion conversions return garbage. Flip rows instead of columns and you move the world instead of relabelling the camera. Check det(R) = +1 when loading poses; it's one line and catches half these bugs.

look_at in each convention, from fwd = normalize(target − eye):

…and t = −R·eye in both. If up is parallel to fwd the cross product is zero and the basis is NaN — production code switches to another up vector.

5. Projecting a point

#include <optional>
#include <print>

struct V3 { float x, y, z; };
struct Pixel { float u, v, depth; };

// Camera-space point (OpenCV axes) to pixel. Check BEFORE dividing.
std::optional<Pixel> project(V3 c, float fx, float fy, float cx, float cy) {
    if (!(c.z > 0)) return std::nullopt;         // behind, on the plane, or NaN
    return Pixel{fx * c.x / c.z + cx, fy * c.y / c.z + cy, c.z};
}

int main() {
    for (V3 c : {V3{1, -0.5f, 10}, V3{1, -0.5f, -10}, V3{0, 0, 0}}) {
        if (auto p = project(c, 100, 100, 64, 48))
            std::println("({}, {}, {}) -> u = {}, v = {}, depth {}", c.x, c.y, c.z, p->u, p->v, p->depth);
        else
            std::println("({}, {}, {}) -> not visible", c.x, c.y, c.z);
    }
    // What a divide-first version would compute for the point behind the camera:
    std::println("mirrored: u = {}, v = {}", 100 * 1 / -10.0f + 64, 100 * -0.5f / -10.0f + 48);
}

The point behind the camera, divided anyway, lands at (54, 53) — a perfectly plausible pixel, mirrored through the centre. After the divide nothing tells you it was behind: the sign of z is the only witness, so test it first. Write the test as !(z > 0), which also rejects z = 0 (±inf) and NaN; a point in front but outside the image still projects, and the caller decides (splatting code keeps a guard band so Gaussians whose centre is just off-screen still contribute).

In a full pipeline you don't stop at z > 0: points closer than the near plane are clipped, and primitives that straddle the camera plane must be clipped in 3D before the divide (see the rasterization lesson).

6. Pixel centres and primary rays

Pixel (i, j) is a little square [i, i+1) × [j, j+1); the sample that represents it is its centre, (i + 0.5, j + 0.5). Ray tracers and NeRF shoot a ray through that point; rasterizers test coverage at it.

A primary ray runs the projection backwards: undo the intrinsics at the pixel centre to get a camera-space direction with z = 1, rotate it into the world with the camera-to-world rotation Rᵀ, normalize, and start it at the centre −Rᵀt:

Ray pixel_ray(int i, int j) const {
    float u = i + 0.5f, v = j + 0.5f;
    Vec3 d_cam{(u - K.cx) / K.fx, (v - K.cy) / K.fy, 1.0f};
    return {centre(), normalize(transpose(R) * d_cam)};
}

Unprojection is the same with a known depth: the camera-space point is d_cam · depth (so its z is depth), then camera_to_world. That's how a depth map becomes a point cloud, and how 3DGS initialises Gaussians from COLMAP points' reprojections. Mind which "depth" a dataset means: camera-space z (what project returns) or distance along the ray (|x_c|). They differ by the factor |d_cam|, which is 1 only at the principal point.

In Python (NumPy / PyTorch). NeRF's get_rays builds all rays at once: i, j = meshgrid(arange(W), arange(H), indexing="xy"), dirs = stack([(i + 0.5 - cx) / fx, -(j + 0.5 - cy) / fy, -ones_like(i)], -1) (OpenGL axes: y and z negated), then rays_d = dirs @ c2w[:3, :3].T, rays_o = c2w[:3, 3]. The original code omitted the + 0.5; later implementations (Mip-NeRF, nerfstudio) put it back.

7. The OpenGL projection matrix, NDC and depth

A rasterizer wants more than (u, v): it needs a depth per pixel for the z-buffer and a box to clip against. The projection matrix maps the view frustum (OpenGL view space: looking down −z) into clip space; the GPU then divides by w to get normalized device coordinates (NDC), the cube [−1, 1]³; finally the viewport transform maps NDC x, y to pixels (flipping y, because NDC y is up and pixel rows go down).

          ⎡ f/aspect  0      0               0          ⎤      f = 1 / tan(fovy / 2)
P_gl =    ⎢ 0         f      0               0          ⎥
          ⎢ 0         0  (far+near)/(near−far)  2·far·near/(near−far) ⎥
          ⎣ 0         0     −1               0          ⎦

The last row makes w_clip = −z_view — the depth — so the divide by w is the pinhole divide. With aspect = W/H, fy = H/(2 tan(fovy/2)) and square pixels, P_gl · view followed by the viewport lands on exactly the pixel that K gives; the build task checks this.

Clipping happens in clip space, before the divide: a vertex is inside when −w ≤ x, y, z ≤ w. Dividing first would lose the sign of w, just as in §5.

The third row is the interesting one. z can't be kept linear — after the divide it must be of the form A + B/depth — so solving for −near → −1, −far → +1 gives a hyperbola:

With a fixed number of depth values, surfaces far away get only a few, and nearby-but-distinct surfaces z-fight. Two remedies:

• Push the near plane out. Precision depends on far/near and is dominated by near; moving near from 0.01 to 0.1 buys more than moving far from 1000 to 100.
• Reverse-Z. Map near → 1 and far → 0 (often with far at infinity: z = n/d) into a 32-bit float depth buffer and test with GREATER. Floats are dense near zero, exactly where the hyperbola is flat, and the two cancel.

#include <bit>
#include <cstdint>
#include <cstdlib>
#include <print>

// Depth in [0, 1] (Direct3D / Vulkan style) for a surface at distance d.
float standard_z(float d, float n, float f) { return f / (f - n) * (1.0f - n / d); }
float reverse_z(float d, float n) { return n / d; }   // reverse-Z, infinite far plane

int ulps_apart(float a, float b) {   // how many floats lie between a and b
    return std::abs(static_cast<int>(std::bit_cast<std::uint32_t>(a)) -
                    static_cast<int>(std::bit_cast<std::uint32_t>(b)));
}

int main() {
    const float n = 0.1f, f = 1000.0f;
    for (float d : {1.0f, 100.0f, 500.0f}) {
        float d2 = d * 1.0001f;   // a surface 0.01 % further away
        std::println("d = {:>5}: standard {:>4} ulps apart, reverse-Z {:>4} ulps apart", d,
                     ulps_apart(standard_z(d, n, f), standard_z(d2, n, f)),
                     ulps_apart(reverse_z(d, n), reverse_z(d2, n)));
    }
}

At 500 units the standard mapping can't tell the two surfaces apart at all (0 ulps); reverse-Z still has over a thousand float values between them. This is the default in modern engines (Unreal, most Vulkan/D3D12 renderers).

Note the split: NDC depth is affine in screen space and is interpolated linearly across a triangle, while colours and UVs need perspective-correct interpolation — the next lesson.

8. Where this shows up in neural rendering

• 3D Gaussian Splatting (Kerbl et al. 2023) projects each Gaussian's mean with exactly project (view matrix, then the projection matrix) and its covariance with the Jacobian of the perspective map, J = [[fx/z, 0, −fx·x/z²], [0, fy/z, −fy·y/z²]] (Zwicker et al. 2001, EWA splatting). Gaussians with z below a small near threshold are culled before that divide.
• NeRF (Mildenhall et al. 2020) and Instant-NGP (Müller et al. 2022) generate one ray per pixel centre and march it; ray directions are kept in world space and the sampled depths along them are what the network sees.

In an interview. "Given a 4×4 view matrix and a projection matrix, where is the camera, and what's the vertical FOV?" The camera centre is −Rᵀt from the view matrix's top-left 3×3 and last column; fovy = 2·atan(1/P[1][1]). Then expect "why is depth precision non-linear and what do you do about it?" (§7) and "a point behind the camera shows up on screen — why?" (§5).

Summary

• x_c = R·x_w + t, then u = fx·x/z + cx, v = fy·y/z + cy; check z > 0 before dividing.
• c2w and w2c are inverses: R_c2w = Rᵀ, centre = −Rᵀt. t is not the camera position.
• fx = W / (2 tan(fovx/2)), radians, half angle.
• Pixel (i, j) is sampled at (i + 0.5, j + 0.5).
• OpenCV: y down, +z forward. OpenGL: y up, −z forward. Convert c2w by flipping the y and z columns; check det(R) = +1.
• NDC depth is hyperbolic in distance; reverse-Z with float depth fixes the precision.