Transformations & Homogeneous Coordinates

A linear map from ℝ³ → ℝ³ satisfies f(u + v) = f(u) + f(v) and f(αu) = αf(u). This forces f(0) = 0: every linear map fixes the origin. Rotation and uniform scaling are linear — they fix the origin. But translation moves the origin, so it is fundamentally not a linear map in the usual 3D space.

Formally, suppose we want T_t(x) = x + t for some fixed vector t ≠ 0. Then T_t(0) = t ≠ 0, violating linearity. No 3×3 matrix can represent this.

This forces us into a choice: either represent transforms as (matrix, vector) affine pairs and handle the special-case arithmetic everywhere — messy, error-prone, hard to compose — or find a new space where translation becomes linear. Homogeneous coordinates are that new space.

The practical cost of not having a uniform representation is severe: composing a rotation R followed by a translation t followed by another rotation R' would require tracking (R, t) pairs and carefully writing out (R', 0) ∘ (R, t) = (R'R, R't). As soon as you add projections this algebra explodes. The uniformity of 4×4 matrices is what makes the GPU pipeline tractable.

import numpy as np

# Attempting translation as a 3x3 matrix — impossible
R = np.eye(3)  # identity rotation
t = np.array([3.0, 4.0, 0.0])  # desired translation

# There is NO 3x3 matrix M such that M @ x == x + t for all x
# Proof: M @ [0,0,0] = [0,0,0] always, but x+t at x=0 is t != 0

The key idea: embed 3D Euclidean space into 4D projective space ℝP³. A 3D point p = (x, y, z) is represented as the 4-vector (x, y, z, 1). More generally, any 4-vector (wx, wy, wz, w) with w ≠ 0 represents the Euclidean point (x/w, y/w, z/w). The w-coordinate is the homogeneous weight.

This representation is only defined up to overall scale: (2, 4, 6, 2) and (1, 2, 3, 1) and (3, 6, 9, 3) all represent the same point (1, 2, 3). This is why it is called projective space — a point in ℝP³ is really an equivalence class (a ray through the origin of ℝ⁴).

Points at infinity: when w = 0, the expression (x/w, y/w, z/w) is undefined in Euclidean terms. We interpret (x, y, z, 0) as the point at infinity in direction (x, y, z). This is not a hack — it is geometrically meaningful: parallel lines meet at a point at infinity in projective geometry, which is exactly what perspective projection exploits.

Vectors vs points: in homogeneous coordinates, the convention is:

• Point: (x, y, z, 1)
• Direction/vector: (x, y, z, 0)

This is consistent: adding two direction vectors gives a direction; adding a point and a direction gives a point (translation); subtracting two points gives a direction. Importantly, a 4×4 matrix that translates points leaves direction vectors unchanged — which is exactly the desired behavior for surface normals and directions.

import numpy as np

def to_homogeneous_point(p):
    """3D point → homogeneous 4-vector"""
    x, y, z = p
    return np.array([x, y, z, 1.0])

def to_homogeneous_vec(v):
    """3D direction vector → homogeneous 4-vector (w=0)"""
    x, y, z = v
    return np.array([x, y, z, 0.0])

def from_homogeneous(h):
    """Homogeneous 4-vector → 3D point (perspective divide)"""
    x, y, z, w = h
    assert abs(w) > 1e-10, "Point at infinity"
    return np.array([x/w, y/w, z/w])

p = to_homogeneous_point([1, 2, 3])   # [1, 2, 3, 1]
v = to_homogeneous_vec([0, 1, 0])     # [0, 1, 0, 0]  (up direction)
print(p + v)  # [1, 3, 3, 1]  → point (1, 3, 3)  ✓

With homogeneous coordinates established, every fundamental 3D transform becomes a 4×4 matrix.

Translation by (tx, ty, tz):

T = | 1  0  0  tx |
    | 0  1  0  ty |
    | 0  0  1  tz |
    | 0  0  0   1 |

Multiplying T by (x, y, z, 1)ᵀ gives (x+tx, y+ty, z+tz, 1)ᵀ. Direction vectors (w=0) are unaffected — correct, since translating a direction is meaningless.

Uniform scale by s (or non-uniform by sx, sy, sz):

S = | sx  0   0   0 |
    |  0  sy  0   0 |
    |  0   0  sz  0 |
    |  0   0   0  1 |

Rotation about the Z-axis by angle θ:

Rz = | cosθ  -sinθ  0  0 |
     | sinθ   cosθ  0  0 |
     |    0       0  1  0 |
     |    0       0  0  1 |

Similar forms exist for Rx (rotation around X) and Ry (around Y).

Key property: the inverse of each matrix is easy:

• T(tx,ty,tz)⁻¹ = T(-tx,-ty,-tz)
• S(sx,sy,sz)⁻¹ = S(1/sx,1/sy,1/sz)
• R⁻¹ = Rᵀ (rotation matrices are orthogonal)

Normal transform: when you transform geometry by matrix M, surface normals transform by the inverse-transpose (M⁻¹)ᵀ, not by M itself. This is because normals are co-vectors (elements of the dual space) — they must stay perpendicular to the surface they describe. Non-uniform scaling will shear normals incorrectly if you apply M directly.

import numpy as np

def translation(tx, ty, tz):
    T = np.eye(4)
    T[0,3], T[1,3], T[2,3] = tx, ty, tz
    return T

def scale(sx, sy, sz):
    return np.diag([sx, sy, sz, 1.0])

def rot_z(theta):
    c, s = np.cos(theta), np.sin(theta)
    return np.array([
        [ c, -s, 0, 0],
        [ s,  c, 0, 0],
        [ 0,  0, 1, 0],
        [ 0,  0, 0, 1]
    ], dtype=float)

def normal_matrix(M):
    """Correct transform for surface normals = (M^-1)^T"""
    return np.linalg.inv(M).T

# Verify: rotation inverse = transpose
R = rot_z(np.pi / 3)
print(np.allclose(R @ R.T, np.eye(4)))  # True

The most common source of bugs in graphics is applying transforms in the wrong order. Matrix multiplication is not commutative: rotate-then-translate ≠ translate-then-rotate.

The convention: to apply transform A first, then B, the combined matrix is M = B · A (B on the left). This is because we write column vectors on the right: M·v = B·(A·v).

A vivid analogy: imagine you are in a room (local frame). If you rotate 90° then walk forward, you end up in a different place than if you walk forward then rotate 90°. The same non-commutativity applies to matrices.

The Model-World-View chain: a typical object in a scene has:

1. Model matrix M: places the object in world space (scale, rotate, translate)
2. View matrix V: transforms world space into camera space (the camera looks down -Z)
3. Projection matrix P: maps camera space to clip space

The combined MVP matrix = P · V · M transforms a model-space vertex all the way to clip space in one shot.

Reading a transform chain: The matrix M = T · Ry · S should be read right-to-left as "first scale, then rotate around Y, then translate." This corresponds to column-vector convention. Some engines (OpenGL historical convention, row-vector HLSL) use row vectors and write M = S · Ry · T instead — always check the convention of the system you are in.

Hierarchical transforms in a scene graph: a forearm is attached to an upper arm, which is attached to a torso. The forearm's world transform is M_torso · M_upper_arm · M_forearm. Changing M_torso automatically moves everything attached to it — this is why scene graphs work.

import numpy as np

def rot_y(theta):
    c, s = np.cos(theta), np.sin(theta)
    return np.array([
        [ c, 0, s, 0],
        [ 0, 1, 0, 0],
        [-s, 0, c, 0],
        [ 0, 0, 0, 1]
    ], dtype=float)

# Demonstrate non-commutativity
t = translation(1, 0, 0)
R = rot_y(np.pi / 2)

# A unit point on +Z axis
p = np.array([0, 0, 1, 1], dtype=float)

# Rotate first, then translate
result1 = (t @ R) @ p  # Rotate: (0,0,1)→(1,0,0); Translate: →(2,0,0)
print("Rotate then translate:", result1[:3])  # [2, 0, 0]

# Translate first, then rotate
result2 = (R @ t) @ p  # Translate: (0,0,1)→(1,0,1); Rotate: →(1,0,-1)
print("Translate then rotate:", result2[:3])  # [1, 0, -1]

# They differ! Order matters.

The view transform (also called the camera transform or LookAt transform) converts coordinates from world space into camera (eye) space, where:

• The camera is at the origin
• The camera looks along the −Z axis (OpenGL convention; Direct3D uses +Z into the screen)
• +Y is up, +X is to the right

We define the camera by three parameters:

• eye e: camera position in world space
• center c: the point the camera is aimed at
• up u: an approximate up-vector (typically world Y)

From these we compute an orthonormal camera frame:

forward = normalize(e - c)    # points from scene toward camera (-Z in camera space)
right   = normalize(up × forward)  # +X in camera space  
newUp   = forward × right          # +Y in camera space (orthogonalized)

The view matrix is then the inverse of the camera's world transform. Since the camera frame is orthonormal, its inverse is just its transpose (for the rotation part) composed with a negating translation:

V = R^T · T(-e)

Where R has right, newUp, forward as its rows. This matrix transforms any world-space point into camera space.

Intuition: We are not moving the camera; we are applying the inverse of the camera's transformation to everything in the world so that the camera ends up at the origin looking down −Z.

A common gotcha: if up is parallel to forward, the cross product is zero and the camera frame is degenerate. Always handle this edge case (e.g., switch to a different up vector like world X).

import numpy as np

def look_at(eye, center, up):
    """Build a view matrix (OpenGL convention: camera looks down -Z)."""
    eye, center, up = map(np.asarray, (eye, center, up))
    f = eye - center          # forward (points toward camera)
    f = f / np.linalg.norm(f)
    r = np.cross(up, f)       # right
    r = r / np.linalg.norm(r)
    u = np.cross(f, r)        # true up (orthogonalized)

    # R^T rows are the camera basis vectors
    V = np.array([
        [r[0],  r[1],  r[2], -np.dot(r, eye)],
        [u[0],  u[1],  u[2], -np.dot(u, eye)],
        [f[0],  f[1],  f[2], -np.dot(f, eye)],
        [0,     0,     0,     1             ]
    ], dtype=float)
    return V

# Camera at (0,3,5) looking at origin, Y up
V = look_at([0,3,5], [0,0,0], [0,1,0])
print(V)

The pinhole camera model is the idealized model of how a real camera (or the human eye) maps 3D space onto a 2D image. A single point of light (the pinhole / aperture) projects all rays onto an image plane behind it. By similar triangles, a 3D point at (x, y, z) in camera space (with the camera at the origin looking down -Z) projects onto the image plane at distance d as:

x' = d · x / (-z)     (−z because z is negative in front of the camera)
y' = d · y / (-z)

This is a perspective divide: distant objects appear smaller because x' and y' are divided by depth.

Now the magic of homogeneous coordinates: this nonlinear perspective divide can be encoded as a linear matrix operation if we allow the output w-coordinate to be nonzero!

The perspective projection matrix (pre-perspective-divide) is:

P = | f/a   0    0             0        |
    |  0    f    0             0        |
    |  0    0   (n+f)/(n-f)  2nf/(n-f)  |
    |  0    0   -1            0         |

where f = 1/tan(fov/2) is the focal length scaling, a = width/height is the aspect ratio, n and f are the near and far clip plane distances.

The third row maps z into a depth value in [-1,1] for NDC (this mapping is nonlinear in z — see the depth buffer discussion). The critical piece is the fourth row: it writes -z into the w-component of the output. After the matrix multiply, the vertex is at (x', y', z', -z). The perspective divide (dividing all coordinates by w = -z) then completes the projection:

NDC coordinates: (x'/-z, y'/-z, z'/-z)

This is why perspective projection is 'linear in homogeneous coordinates but nonlinear in Euclidean coordinates' — the division by depth is deferred to the perspective divide step.

import numpy as np

def perspective(fov_y, aspect, near, far):
    """
    OpenGL-style perspective matrix.
    Maps camera space into NDC clip space.
    fov_y: vertical field of view in radians
    Maps z in [-near, -far] to NDC z in [-1, 1]
    """
    f = 1.0 / np.tan(fov_y / 2)
    nf = near - far
    return np.array([
        [f/aspect, 0,              0,              0       ],
        [0,        f,              0,              0       ],
        [0,        0,  (far+near)/nf,  2*far*near/nf],
        [0,        0,             -1,              0       ]
    ], dtype=float)

# Project a point at camera-space position (1, 2, -5)
p_cam = np.array([1.0, 2.0, -5.0, 1.0])
P = perspective(np.radians(60), 16/9, 0.1, 100)
clip = P @ p_cam        # In clip space: w = -(-5) = 5
ndc  = clip / clip[3]   # Perspective divide
print("Clip:", clip)
print("NDC xy: ", ndc[:2])  # x,y in [-1,1]

Orthographic projection (also called parallel projection) maps 3D camera space to NDC by a pure scaling and translation, with no perspective divide — all parallel lines remain parallel, and size does not vary with depth.

The projection maps a view volume (an axis-aligned box defined by left, right, bottom, top, near, far) to the NDC cube [-1,1]³:

P_ortho = | 2/(r-l)     0       0     -(r+l)/(r-l) |
          |    0     2/(t-b)    0     -(t+b)/(t-b) |
          |    0        0    -2/(f-n)  -(f+n)/(f-n) |
          |    0        0       0          1        |

The w-component remains 1 after this multiply, so there is no perspective divide (dividing by 1 is a no-op). This is the essential difference from perspective: the fourth row is [0, 0, 0, 1] instead of [0, 0, -1, 0].

When to use orthographic: architectural and CAD drawings (where preserving measurements matters), isometric game views, shadow map rendering (directional lights behave like an orthographic camera), and UI/2D rendering. Wherever you need 'things look the same size regardless of distance', use orthographic.

When to use perspective: anything meant to look photorealistic, because that is how real cameras and eyes work. The perspective projection matches the physics of light converging to a point.

Note that orthographic projection makes the depth buffer (z-buffer) mapping linear in depth, which is why z-fighting is less problematic with orthographic cameras. With perspective, the depth buffer has much better precision near the near plane and poor precision far away — a fundamental issue in practice.

import numpy as np

def orthographic(left, right, bottom, top, near, far):
    """Orthographic projection matrix (OpenGL convention)."""
    rl = right - left
    tb = top - bottom
    fn = far - near
    return np.array([
        [2/rl,    0,     0,  -(right+left)/rl],
        [0,    2/tb,     0,  -(top+bottom)/tb],
        [0,       0, -2/fn,    -(far+near)/fn],
        [0,       0,     0,                1 ]
    ], dtype=float)

# A point at (1, 0, -3) with a view box [-2,2] x [-2,2] x [-1,-10]
p = np.array([1.0, 0.0, -3.0, 1.0])
P = orthographic(-2, 2, -2, 2, 1, 10)
ndc = P @ p   # w remains 1: no perspective divide needed
print("NDC:", ndc[:3])  # Should be in [-1,1]^3

Putting it all together, the vertex processing pipeline for a single vertex v in model space is:

1. Model space   : v
2. World space   : v_world = M · v           (Model matrix)
3. Camera space  : v_cam   = V · v_world     (View matrix)
4. Clip space    : v_clip  = P · v_cam       (Projection matrix)
   Combined:       v_clip  = (P·V·M) · v     (MVP matrix)
5. NDC           : v_ndc   = v_clip / v_clip.w  (Perspective divide)
6. Screen space  : v_screen = viewport_transform(v_ndc)

Clip space is where clipping (discarding geometry outside the view frustum) happens. The GPU clips against the six frustum planes, which in clip space are simply ±w planes: -w ≤ x ≤ w, -w ≤ y ≤ w, 0 ≤ z ≤ w (DirectX) or -w ≤ z ≤ w (OpenGL). Clipping in homogeneous space before the perspective divide avoids divisions by zero for geometry crossing w=0.

NDC (Normalized Device Coordinates): after the perspective divide, x and y are in [-1,1] and z is in [-1,1] (OpenGL) or [0,1] (DirectX/Vulkan/Metal). Everything outside this cube has been clipped away.

Viewport transform: maps NDC to actual pixel coordinates. If the viewport is (x0, y0) with width w and height h:

px = (ndc_x + 1) * w/2 + x0
py = (ndc_y + 1) * h/2 + y0

(Note OpenGL has Y=0 at bottom; framebuffer images often have Y=0 at top — another common source of vertically flipped images.)

The GPU's job: computes the MVP transform for every vertex in parallel (vertex shader), then interpolates vertex attributes (normals, texture coords) across triangles, applies per-fragment operations (fragment shader), and writes to the framebuffer. The entire transform chain we derived is what happens inside gl_Position = MVP * position.

import numpy as np

def viewport_transform(ndc_x, ndc_y, width, height, x0=0, y0=0):
    """Maps NDC [-1,1]^2 to pixel coordinates."""
    px = (ndc_x + 1.0) * 0.5 * width  + x0
    py = (ndc_y + 1.0) * 0.5 * height + y0
    return px, py

# Full pipeline: model -> screen
def render_vertex(v_model, M, V, P, viewport_w, viewport_h):
    MVP = P @ V @ M
    v_clip = MVP @ v_model
    # Perspective divide
    w = v_clip[3]
    if abs(w) < 1e-7:
        return None  # at infinity
    ndc = v_clip[:3] / w
    px, py = viewport_transform(ndc[0], ndc[1], viewport_w, viewport_h)
    return px, py, ndc[2]  # screen x, y, depth

# Example: place a unit cube vertex, look at it
M_model = np.eye(4)                        # identity model transform
V_view  = look_at([0,0,5],[0,0,0],[0,1,0]) # camera 5 units back
P_proj  = perspective(np.radians(60), 16/9, 0.1, 100)

vertex = np.array([1.0, 1.0, 0.0, 1.0])   # top-right-front of cube
result = render_vertex(vertex, M_model, V_view, P_proj, 1920, 1080)
print(f"Screen: ({result[0]:.1f}, {result[1]:.1f}), depth: {result[2]:.4f}")

After projection, the z-coordinate in NDC does not vary linearly with camera-space depth. Let z_cam be the camera-space depth (negative, in [-near, -far]). The NDC z after perspective projection is:

z_ndc = (f+n)/(f-n) + 2fn / ((f-n) * z_cam)

This is a rational (hyperbolic) function of z_cam, not linear. What does this mean in practice?

1. Depth precision is unevenly distributed: most of the NDC depth range is spent near the near plane. The near plane gets half the depth buffer's precision; the far half of the scene shares the other half. This is why the ratio far/near dominates depth buffer quality — setting near=0.001, far=100000 gives terrible z-fighting, while near=1, far=1000 gives acceptable precision.

2. Never linearly interpolate clip-space depth: when rasterizing a triangle, you cannot linearly interpolate z in screen space. You must interpolate 1/z (which is linear in screen space) and then recover z. This is called perspective-correct interpolation and applies to texture coordinates and all other vertex attributes too.

3. Reversed Z / floating-point depth buffers: modern engines often use reversed depth (map near → 1.0, far → 0.0) combined with a floating-point depth buffer. This dramatically improves precision for distant geometry because floating-point numbers have more precision near 0, and with reversed Z, distant geometry maps near 0. This is the "reverse Z" trick used in game engines like Frostbite and Unreal.

Orthographic depth is linear: because the orthographic projection matrix's third row maps z linearly (no perspective divide), the depth buffer values are evenly distributed over the scene depth range. This is one more reason to use orthographic for technical renderings.

import numpy as np
import matplotlib.pyplot as plt

# Visualize the nonlinear z-mapping
near, far = 1.0, 100.0
z_cam = np.linspace(-near, -far, 500)  # camera-space depth

# OpenGL perspective: maps [-near, -far] → [-1, 1]
z_ndc = (far + near) / (far - near) + (2 * far * near) / ((far - near) * z_cam)

# Notice: most precision is near z_cam = -near
print("NDC z at z_cam=-1:",   z_ndc[0])   # should be -1 (near plane)
print("NDC z at z_cam=-50:",  z_ndc[250])  # should be near +1 already!
print("NDC z at z_cam=-100:", z_ndc[-1])   # should be +1 (far plane)
# The entire far half (50 to 100) maps to only a tiny NDC range!

Three sequential Euler angle rotations (yaw, pitch, roll) are the most intuitive way to specify an orientation, but they suffer from gimbal lock: when the pitch reaches ±90°, the yaw and roll axes align, and a degree of freedom is lost. The rotation matrix formulation loses no information, but Euler angles as a parameterization become singular.

This is not just a theoretical problem — Apollo 13's guidance computer suffered gimbal lock (in its physical 3-ring gimbal mechanism), and countless 3D applications have produced animation discontinuities at ±90° pitch.

Rotation matrices represent orientations exactly (as 3×3 orthogonal matrices with determinant +1, i.e., SO(3)), but their 9-parameter redundancy makes interpolation hard.

Axis-angle representation: a rotation by angle θ around unit axis n can be encoded as just 4 numbers (nx, ny, nz, θ) or 3 numbers (θnx, θny, θnz). This is compact but interpolation is still awkward.

Quaternions: a quaternion q = (qx, qy, qz, qw) represents a rotation as:

q = (n·sin(θ/2), cos(θ/2))

Quaternions avoid gimbal lock, use only 4 floats (vs 16 for a matrix), compose efficiently via quaternion multiplication (fewer multiplications than matrix multiply), and interpolate smoothly via SLERP (Spherical Linear Interpolation). They are the standard rotation representation in game engines and animation systems.

For the pipeline, quaternions are ultimately converted to 4×4 rotation matrices for the vertex transform since the GPU works with matrices. But internally, skeletons, cameras, and physics objects store orientations as quaternions and convert only for rendering.

import numpy as np

def quat_to_mat4(q):
    """Convert unit quaternion (qx,qy,qz,qw) to 4x4 rotation matrix."""
    qx, qy, qz, qw = q / np.linalg.norm(q)
    return np.array([
        [1-2*(qy*qy+qz*qz), 2*(qx*qy-qz*qw),   2*(qx*qz+qy*qw),   0],
        [2*(qx*qy+qz*qw),   1-2*(qx*qx+qz*qz), 2*(qy*qz-qx*qw),   0],
        [2*(qx*qz-qy*qw),   2*(qy*qz+qx*qw),   1-2*(qx*qx+qy*qy), 0],
        [0,                  0,                  0,                  1]
    ], dtype=float)

def slerp(q1, q2, t):
    """Spherical linear interpolation between two quaternions."""
    q1, q2 = q1/np.linalg.norm(q1), q2/np.linalg.norm(q2)
    dot = np.clip(np.dot(q1, q2), -1, 1)
    if dot < 0:  # take the shorter arc
        q2, dot = -q2, -dot
    if dot > 0.9995:  # nearly identical: fall back to lerp
        return q1 + t*(q2-q1)
    theta0 = np.arccos(dot)
    theta  = theta0 * t
    return (np.sin(theta0-theta)*q1 + np.sin(theta)*q2) / np.sin(theta0)

# 90-degree rotation around Y axis as a quaternion
angle = np.pi / 2
q = np.array([0, np.sin(angle/2), 0, np.cos(angle/2)])
print(quat_to_mat4(q))