3D math for graphics: vectors and matrices

Every renderer — a rasterizer, a path tracer, a 3D Gaussian splatting viewer — starts with the same twenty functions: dot, cross, normalize, a 4x4 matrix product, "transform this point", "transform this normal", "invert this camera". They're short, and almost every bug in them is a convention bug: a sign, an index order, a w. This lesson pins the conventions down — the ones GLM and OpenGL use, which are also what most research code (NeRF, 3DGS, Instant-NGP viewers) uses — and derives the formulas so you can check them instead of memorising them.

1. Vectors: dot, cross, length

A vec3 is three floats. Keep it a plain aggregate with constexpr operators: it's trivially copyable, passes in registers, and can be checked at compile time.

• dot: a·b = aₓbₓ + a_yb_y + a_zb_z = |a||b| cos θ. Zero means perpendicular; positive means "same side". dot(n, l) is diffuse lighting; dot(a, a) is the squared length (no square root — use it to compare distances).
• cross: a × b = (a_y b_z − a_z b_y, a_z b_x − a_x b_z, a_x b_y − a_y b_x). Perpendicular to both, with length |a||b| sin θ (the area of the parallelogram). Each component uses the other two axes in cyclic order x → y → z → x, which is how to remember it: x × y = z, y × z = x, z × x = y, and reversing the order flips the sign.
• normalize: divide by the length — unless the length is zero. 0/0 is NaN, and NaN poisons everything it touches: one degenerate triangle normal becomes a black pixel, then (through temporal accumulation) a black blob.

#include <cmath>
#include <print>

struct vec3 { float x = 0, y = 0, z = 0; };

constexpr vec3 operator-(vec3 a, vec3 b) { return {a.x - b.x, a.y - b.y, a.z - b.z}; }
constexpr vec3 operator/(vec3 a, float s) { return {a.x / s, a.y / s, a.z / s}; }
constexpr float dot(vec3 a, vec3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
constexpr vec3 cross(vec3 a, vec3 b) {
    return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
}
inline vec3 normalize(vec3 a) {
    float len = std::sqrt(dot(a, a));
    return len > 0.0f ? a / len : vec3{};   // zero stays zero: no NaN
}

static_assert(cross({1, 0, 0}, {0, 1, 0}).z == 1.0f);   // right-handed

int main() {
    // Triangle normal: counter-clockwise vertices seen from the front.
    vec3 p0{0, 0, 0}, p1{1, 0, 0}, p2{0, 1, 0};
    vec3 n = normalize(cross(p1 - p0, p2 - p0));
    std::println("n = ({}, {}, {})", n.x, n.y, n.z);
    // A degenerate (zero-area) triangle:
    vec3 d = normalize(cross(p1 - p0, p1 - p0));
    std::println("degenerate = ({}, {}, {})", d.x, d.y, d.z);
}

Handedness. The formula for the cross product is fixed; what changes between conventions is how you draw the axes. In a right-handed frame (OpenGL, Blender, PBRT, COLMAP, most of ML-land) x right, y up gives z towards you; classic Direct3D and Unity are left-handed (z away). The cross-product formula doesn't know or care: a "left-handed engine" just draws its z axis the other way, so the same numbers look mirrored. The practical consequences are the triangle winding that counts as "front" (counter-clockwise in GL by default) and the sign of z in the projection matrix.

2. Points, directions and w

A 4x4 matrix can do translation only because we pad 3D vectors with a fourth coordinate:

• a point is (x, y, z, 1) — the translation column is multiplied by 1 and added;
• a direction (a difference of two points: an edge, a velocity, a ray direction) is (x, y, z, 0) — the translation is multiplied by 0. That's right: (p + t) − (q + t) = p − q.

w is a type tag the algebra respects: point − point = direction (1 − 1 = 0), point + direction = point. After a projection matrix w becomes the depth, and dividing by it — the perspective divide — is what makes distant things small. So transform_point multiplies with w = 1 and divides by the resulting w when it isn't 1.

3. Matrices: storage order vs vector convention

Two independent choices hide behind "row-major" and "column-major", and mixing them up is the most common matrix bug there is.

1. Vector convention. Column vectors (v' = M·v, math textbooks, OpenGL, GLM, PBRT, Eigen) or row vectors (v' = v·M, classic Direct3D / DirectXMath). The matrices of the two conventions are transposes of each other: with column vectors the translation is the last column, with row vectors the last row.
2. Storage order. How 16 floats are laid out in memory: column-major (each column contiguous — OpenGL, GLM, Eigen's default) or row-major (each row contiguous — C 2D arrays, DirectXMath, NumPy's default).

Fabian Giesen's observation: column vectors stored column-major and row vectors stored row-major produce the same 16 floats — the transpose, stored transposed. Either way the translation sits at m[12], m[13], m[14]. This lesson (and the build task) uses column vectors stored column-major, like GLM: element (row r, column c) lives at m[c*4 + r].

The product is the same formula in every convention — only the indexing changes. Write one accessor and never touch m[] anywhere else:

#include <array>
#include <print>

struct mat4 {
    std::array<float, 16> m{};
    constexpr float& operator()(int r, int c) { return m[c * 4 + r]; }
    constexpr float operator()(int r, int c) const { return m[c * 4 + r]; }
    static constexpr mat4 identity() {
        mat4 I;
        for (int i = 0; i < 4; ++i) I(i, i) = 1;
        return I;
    }
};

// (A·B)(i, j) = Σₖ A(i, k)·B(k, j): row i of A against column j of B.
constexpr mat4 operator*(const mat4& a, const mat4& b) {
    mat4 r;
    for (int i = 0; i < 4; ++i)
        for (int j = 0; j < 4; ++j)
            for (int k = 0; k < 4; ++k) r(i, j) += a(i, k) * b(k, j);
    return r;
}

constexpr mat4 translate(float x, float y, float z) {
    mat4 t = mat4::identity();
    t(0, 3) = x; t(1, 3) = y; t(2, 3) = z;
    return t;
}

int main() {
    mat4 t = translate(1, 2, 3) * translate(10, 20, 30);
    std::println("m[12..14] = {} {} {}", t.m[12], t.m[13], t.m[14]);
    // transform a point (w = 1) and a direction (w = 0)
    float p[4] = {0, 0, 0, 1}, d[4] = {1, 0, 0, 0}, tp[4] = {}, td[4] = {};
    for (int i = 0; i < 4; ++i)
        for (int k = 0; k < 4; ++k) { tp[i] += t(i, k) * p[k]; td[i] += t(i, k) * d[k]; }
    std::println("point -> ({}, {}, {}), direction -> ({}, {}, {})",
                 tp[0], tp[1], tp[2], td[0], td[1], td[2]);
}

Uploading to a GPU. glUniformMatrix4fv(loc, 1, GL_FALSE, m.data()) expects column-major; HLSL's default float4x4 packing is column-major too, but HLSL code usually writes mul(v, M) (row vectors), which is the same bytes again. When a shader result is mirrored or the translation appears in w, you have one transpose too many.

4. Composing transforms: T·R·S

With column vectors (A·B)·p = A·(B·p): the matrix nearest the point acts first. A model matrix M = T·R·S therefore scales the object about its own origin, then rotates it, then moves it into place. Any other order is a different transform: scale after rotation shears, and rotation after translation swings the object around the world origin like a door.

Rotation about an axis (Rodrigues). Split v into the part along the unit axis k, which doesn't move, and the part perpendicular to it, which turns by θ in the plane spanned by v⊥ and k × v:

v' = cos θ·v + (1 − cos θ)(k·v) k + sin θ (k × v)

As a matrix: R = cos θ·I + (1 − cos θ)·k kᵀ + sin θ·[k]×, where [k]× = [[0, −k_z, k_y], [k_z, 0, −k_x], [−k_y, k_x, 0]] is the matrix of "cross with k". The sin θ part is antisymmetric — that sign pattern is the handedness. Get one sign wrong and R is no longer orthogonal; swap all of them and you have Rᵀ, the rotation by −θ.

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

struct vec3 { float x, y, z; };

// Rodrigues applied directly to a vector (k must be unit length).
vec3 rotate(vec3 v, vec3 k, float theta) {
    float c = std::cos(theta), s = std::sin(theta);
    float kv = k.x * v.x + k.y * v.y + k.z * v.z;
    vec3 kxv{k.y * v.z - k.z * v.y, k.z * v.x - k.x * v.z, k.x * v.y - k.y * v.x};
    return {c * v.x + (1 - c) * kv * k.x + s * kxv.x,
            c * v.y + (1 - c) * kv * k.y + s * kxv.y,
            c * v.z + (1 - c) * kv * k.z + s * kxv.z};
}

int main() {
    float q = std::numbers::pi_v<float> / 2;
    vec3 a = rotate({1, 0, 0}, {0, 0, 1}, q);   // x about z: should be y
    std::println("({:.3f}, {:.3f}, {:.3f})", a.x, a.y, a.z);
    float n = 1 / std::sqrt(3.0f);              // 120° about (1,1,1) cycles the axes
    vec3 b = rotate({1, 0, 0}, {n, n, n}, 2 * std::numbers::pi_v<float> / 3);
    std::println("({:.3f}, {:.3f}, {:.3f})", b.x, b.y, b.z);
}

({:.3f} hides that cos(π/2) in float is −4.4e-8, not 0 — printed unrounded you'd see -4.371139e-08. Tests must compare with a tolerance.)

5. Inverses

Rigid transforms — rotation plus translation, like a camera pose — have a free inverse. M p = R p + t, so p = Rᵀ(p' − t):

[R | t]⁻¹ = [Rᵀ | −Rᵀ t]

That's the view matrix: the inverse of the camera's pose (camera-to-world, "c2w" in NeRF datasets; COLMAP stores world-to-camera). No elimination, no division, exact up to rounding.

General 4x4 — projections, matrices with shear, anything a user hands you — needs Gauss-Jordan elimination or cofactors. Gauss-Jordan reduces the augmented matrix [A | I] to [I | A⁻¹] with row operations. The pivot for column c must be chosen by partial pivoting: the row at or below the diagonal with the largest |A(r, c)|. Without it, a perfectly invertible matrix with a zero on the diagonal (a permutation; a 90° rotation, whose cos 90° is ≈ 0) divides by zero, and small pivots blow up rounding error. If even the best pivot is (numerically) zero, the matrix squashes space onto a plane and has no inverse: return std::nullopt, don't return garbage.

The cofactor route — A⁻¹ = adj(A) / det(A) — is branch-free and is what GLM and Mesa's gluInvertMatrix do for 4x4s; both are fine at this size. Do the arithmetic in double if you can: a 4x4 inverse in float loses digits quickly when the matrix mixes large translations with small scales.

6. Normals transform by the inverse-transpose

A normal isn't a difference of two points; it's defined by being perpendicular to every tangent t of the surface: nᵀt = 0. Tangents are differences of nearby points, so they transform by M. For the new normal we need n'ᵀ(M t) = 0 whenever nᵀ t = 0; choosing n'ᵀ M = nᵀ works, i.e.

n' = (M⁻¹)ᵀ n — the inverse-transpose of the upper 3x3.

For a rotation R⁻ᵀ = R, and a uniform scale only changes the length (renormalise), which is why this bug hides until someone stretches an object.

#include <print>

int main() {
    // M = scale(2, 1, 1): M⁻ᵀ = scale(1/2, 1, 1).
    float t[3] = {1, -1, 0}, n[3] = {1, 1, 0};   // tangent and normal of a 45° face
    float t2[3] = {2 * t[0], t[1], t[2]};
    float byM[3] = {2 * n[0], n[1], n[2]};
    float byInvT[3] = {0.5f * n[0], n[1], n[2]};
    auto dot = [](const float* a, const float* b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; };
    std::println("tangent·(M n)    = {}", dot(t2, byM));
    std::println("tangent·(M⁻ᵀ n) = {}", dot(t2, byInvT));
}

The cofactor matrix is det(M)·M⁻ᵀ, so since normals are renormalised anyway, some engines use it directly: it needs no division, works for singular M, and (with the sign of det) flips normals correctly under mirroring.

7. Floats: precision and comparisons

A float has a 24-bit significand: about 7 decimal digits, and the gap between adjacent floats grows with magnitude.

• Catastrophic cancellation. Subtracting nearly equal numbers keeps only their rounding error: (1e8f + 1) − 1e8f == 0. Large world coordinates minus a large camera position is exactly this — render camera-relative (subtract the camera position in double, then convert).
• Never == on computed floats. Compare with a tolerance that scales with the magnitudes: |a − b| ≤ ε · max(1, |a|, |b|) (absolute near zero, relative elsewhere). Pick ε from the computation (1e-6 for a few float ops, 1e-4 after a chain of matrix products).
• Degenerate inputs are normal in graphics: zero-length vectors, zero-area triangles, singular matrices from a zero scale. Decide what each function does with them (return zero, nullopt, identity) and test it.

#include <algorithm>
#include <cmath>
#include <print>

bool near(float a, float b, float eps = 1e-5f) {
    return std::fabs(a - b) <= eps * std::max({1.0f, std::fabs(a), std::fabs(b)});
}

int main() {
    float a = 0.0f;
    for (int i = 0; i < 10; ++i) a += 0.1f;   // ten tenths
    std::println("{} == 1? {}   near? {}", a, a == 1.0f, near(a, 1.0f));
    float far = 16384.0f;
    std::println("16384 + 0.001 = {}", far + 0.001f);   // the millimetre is lost
}

In an interview

"Transform a normal" and "invert a camera matrix" are the two classic questions. Say why: normals keep nᵀt = 0, hence M⁻ᵀ; a rigid inverse is [Rᵀ | −Rᵀt] because Rᵀ = R⁻¹. If asked about conventions, separate "column vs row vectors" from "column- vs row-major storage" — and say GLM is column/column, DirectXMath row/row, and that both put translation in m[12..14].

In Python (NumPy / PyTorch)

NumPy arrays are row-major, but that's storage; the math is whatever you write. With an (N, 3) array of points (one per row), the column-vector transform M·p for every point is pts @ R.T + t — the transpose comes from putting points in rows. In PyTorch the same line batches over thousands of Gaussians on the GPU; 3DGS builds every Gaussian's rotation matrix this way (next skill: quaternions).