3D Vision and Neural Rendering

This guide covers the mathematical and practical foundations unique to 3D vision — concepts not found in 2D computer vision or general deep learning. Every technique includes real computations, memory requirements, and actual performance metrics from papers.

Video Version: https://www.youtube.com/watch?v=CoF_ca5NvD4

Part I: 3D Data Representations

Chapter 1: Voxels — 3D Pixels

Real-World Analogy: Like Minecraft blocks — space divided into a 3D grid where each cube (voxel) can be empty or filled, possibly with properties like color or material type.

Basic Voxel Grid

Mathematical Representation:

Voxel grid V ∈ {0,1}^(D×H×W) for binary occupancy

or V ∈ ℝ^(D×H×W×C) for C-channel features

Example: 32×32×32 binary grid = 32,768 voxels
Memory: 32^3 × 1 bit = 4KB (binary)
        32^3 × 32 bits (float) = 128KB

Indexing a Voxel:

# World coordinates to voxel index
def world_to_voxel(point, voxel_size=0.1, grid_origin=[0,0,0]):
    # point = [x, y, z] in meters
    # voxel_size = 0.1 means each voxel is 10cm³
    voxel_idx = np.floor((point - grid_origin) / voxel_size)
    return voxel_idx.astype(int)
  
# Example: point at [1.25, 2.37, 0.84]m
# With 10cm voxels: idx = [12, 23, 8]

Memory Scaling Problem

Resolution vs Memory:

Resolution    Voxels        Memory (float32)
32³          32,768         128 KB
64³          262,144        1 MB
128³         2,097,152      8 MB
256³         16,777,216     64 MB
512³         134,217,728    512 MB
1024³        1,073,741,824  4 GB

Why It Matters: Doubling resolution = 8× memory. This cubic scaling makes high-resolution voxels prohibitive.

Sparse Voxel Representations

Real-World Analogy: Instead of storing every block in a Minecraft world, only store the non-empty ones with their positions.

Hash Table Approach:

class SparseVoxelGrid:
    def __init__(self):
        self.voxels = {}  # (i,j,k) -> value
    def set_voxel(self, i, j, k, value):
        if value != 0:  # Only store non-empty
            self.voxels[(i, j, k)] = value
        elif (i, j, k) in self.voxels:
            del self.voxels[(i, j, k)]
    def get_voxel(self, i, j, k):
        return self.voxels.get((i, j, k), 0)
  
# Memory: Only O(occupied voxels), not O(grid size)
# Typical scene: 5% occupancy → 20× memory savings

Chapter 2: Octrees — Hierarchical Voxels

Real-World Analogy: Like organizing a library — instead of checking every shelf, you first go to the right floor, then section, then shelf. Octrees divide 3D space hierarchically, allowing quick spatial queries.

Octree Structure

Recursive Subdivision:

Level 0: Entire space (1 node)
Level 1: 8 octants (2×2×2)
Level 2: 64 possible nodes (4×4×4)
Level n: 8^n possible nodes
  
Each node either:
- Is a leaf (contains data)
- Has 8 children (subdivided)
- Is empty (pruned)

Memory Efficient Construction:

class OctreeNode:
    def __init__(self, center, size):
        self.center = center  # [x, y, z]
        self.size = size      # Half-width
        self.children = None  # None or list of 8
        self.value = None     # Leaf node data
    def subdivide(self):
        if self.children is not None:
            return
        self.children = []
        s = self.size / 2  # Child size
        # 8 octants: ±x, ±y, ±z combinations
        for dx in [-1, 1]:
            for dy in [-1, 1]:
                for dz in [-1, 1]:
                    child_center = self.center + s * np.array([dx, dy, dz])
                    self.children.append(OctreeNode(child_center, s))

Octree Traversal for Ray Marching

Finding voxels along a ray:

def traverse_octree(node, ray_origin, ray_dir, t_min, t_max):
    # Check if ray intersects node's bounding box
    t_enter, t_exit = ray_box_intersection(
        ray_origin, ray_dir, 
        node.center - node.size,
        node.center + node.size
    )
    if t_exit < t_min or t_enter > t_max:
        return []  # No intersection
    if node.children is None:
        # Leaf node - return if occupied
        return [(node, t_enter)] if node.value else []
    # Recursively check children
    intersections = []
    for child in node.children:
        intersections.extend(
            traverse_octree(child, ray_origin, ray_dir, 
                          max(t_enter, t_min), min(t_exit, t_max))
        )
    return sorted(intersections, key=lambda x: x[1])

Performance:

• Brute force voxel check: O(n) for n voxels along ray
• Octree traversal: O(log n) average case

Chapter 3: Point Clouds — Unstructured 3D

Real-World Analogy: Like a swarm of fireflies in 3D space — each point has a position and possibly color/normal, but no fixed grid structure.

Point Cloud Representation

Data Structure:

Points P ∈ ℝ^(N×3) for positions
Features F ∈ ℝ^(N×C) for colors/normals/descriptors
  
Example: 100K points with RGB
P: 100,000 × 3 × 4 bytes = 1.2 MB (positions)
F: 100,000 × 3 × 1 byte = 300 KB (colors)
Total: ~1.5 MB

K-Nearest Neighbors (KNN) in 3D

Finding nearby points efficiently:

class KDTree:
    def build(self, points, depth=0):
        if len(points) == 0:
            return None
        # Choose axis based on depth
        axis = depth % 3
        # Sort and split at median
        points = sorted(points, key=lambda p: p[axis])
        median = len(points) // 2
        node = {
            'point': points[median],
            'left': self.build(points[:median], depth + 1),
            'right': self.build(points[median + 1:], depth + 1)
        }
        return node
    def nearest(self, root, target, k=5):
        # Returns k nearest neighbors
        # O(log n) average, O(n) worst case
        pass

Point Cloud Convolution

PointNet++ Approach:

def pointnet_conv(points, features, radius=0.1, k=32):
    """
    points: [N, 3] positions
    features: [N, C] input features
    radius: neighborhood radius
    k: max neighbors
    """
    # 1. Sample centroids (furthest point sampling)
    centroids = furthest_point_sample(points, n_samples=N//2)
    # 2. Group neighbors within radius
    groups = []
    for centroid in centroids:
        dists = np.linalg.norm(points - centroid, axis=1)
        neighbors = np.where(dists < radius)[0][:k]
        # 3. Extract local features
        local_points = points[neighbors] - centroid  # Relative positions
        local_features = features[neighbors]
        # 4. Apply PointNet to local region
        aggregated = max_pool(mlp(local_points, local_features))
        groups.append(aggregated)
    return np.array(groups)

Furthest Point Sampling (FPS)

Real-World Analogy: Placing Wi-Fi routers in a large building — install the first router anywhere, then always place the next one at the spot with the worst signal, farthest from every router installed so far. A handful of routers ends up covering the whole building evenly.

Definition: FPS picks a subset of S points from a cloud of N by repeatedly selecting the point with the maximum distance to all points already selected. It is the furthest_point_sample call in the PointNet++ code above—the standard way to choose centroids that cover the shape uniformly.

Get Daniel Stallworth’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

Step-by-Step Code:

def furthest_point_sample(points, n_samples):
    """points: [N, 3] -> indices of n_samples points, O(N * n_samples)."""
    N = len(points)
    if n_samples > N:
        raise ValueError(f"Requested {n_samples} samples from {N} points")
    selected = [0]  # Start from any point (index 0 by convention)
    # Track each point's distance to its nearest selected point
    min_dists = np.linalg.norm(points - points[0], axis=1)
    for _ in range(n_samples - 1):
        next_idx = int(np.argmax(min_dists))  # Worst-covered point
        selected.append(next_idx)
        new_dists = np.linalg.norm(points - points[next_idx], axis=1)
        min_dists = np.minimum(min_dists, new_dists)
    return np.array(selected)

How It Is Used in ML: Every PointNet++-style hierarchy uses FPS to pick the centroids of local neighborhoods, because random sampling oversamples dense regions and can miss thin structures entirely (a chair’s legs vanish, its seat gets 90% of the samples). Edge cases: FPS always selects outliers early — they are by definition far from everything — so denoise first; the result depends on the seed point (fix index 0 or the centroid for reproducibility); and the O(N·S) loop is the sampling bottleneck on large clouds, so implementations run it on GPU.

Chapter 4: Mesh Representations

Real-World Analogy: Like a wireframe model made of triangles — vertices connected by edges form faces that define surfaces.

Mesh Data Structure

Components:

Vertices V ∈ ℝ^(N×3): 3D positions
Faces F ∈ ℤ^(M×3): Triangle indices
Edges E ∈ ℤ^(K×2): Vertex pairs (optional)
Normals N ∈ ℝ^(N×3): Per-vertex normals
UVs ∈ ℝ^(N×2): Texture coordinates
  
Example mesh (cube):
V = [[0,0,0], [1,0,0], [1,1,0], [0,1,0],  # Bottom
     [0,0,1], [1,0,1], [1,1,1], [0,1,1]]   # Top
F = [[0,1,2], [0,2,3],  # Bottom face
     [4,6,5], [4,7,6],   # Top face
     ...]  # 12 triangles total

Mesh Convolution (Graph Neural Networks)

Operating on irregular mesh structure:

def mesh_convolution(vertices, faces, features):
    # Build adjacency from faces
    adjacency = build_adjacency(faces)
    # Graph convolution
    new_features = []
    for v_idx, v_feat in enumerate(features):
        # Aggregate from neighbors
        neighbor_idxs = adjacency[v_idx]
        neighbor_feats = features[neighbor_idxs]
        # Weighted by edge lengths
        edge_vectors = vertices[neighbor_idxs] - vertices[v_idx]
        edge_weights = 1.0 / (np.linalg.norm(edge_vectors, axis=1) + 1e-8)
        # Normalized aggregation
        aggregated = np.average(neighbor_feats, weights=edge_weights, axis=0)
        new_feat = mlp(np.concatenate([v_feat, aggregated]))
        new_features.append(new_feat)
    return np.array(new_features)

Part II: 3D Reconstruction Fundamentals

Chapter 5: Homogeneous Coordinates and Rigid-Body Transformations

Real-World Analogy: Imagine moving a picture frame around a room. The frame can slide across the wall (translation) and tilt (rotation), but the picture inside never stretches or shears — the corners always stay the same distance apart. Homogeneous coordinates are like putting a “1” sticker on every corner of the frame so a single rubber stamp operation can slide and tilt it at once, instead of needing one tool for sliding and a separate tool for turning.

Definition

A rigid-body transformation moves points without changing distances or angles: only rotation and translation, no scaling or shearing. The set of all such transforms in 3D is called SE(3). The problem is that rotation is a matrix multiplication while translation is an addition:

p' = R p + t    (two different operations)

Homogeneous coordinates fix this by appending a 1 to every point, so both operations become one 4×4 matrix multiplication:

3D point p = [x, y, z]  →  homogeneous p̃ = [x, y, z, 1]

Rigid transform as a single 4×4 matrix:
T = [ R  t ]    R ∈ ℝ³ˣ³ (rotation), t ∈ ℝ³ (translation)
    [ 0  1 ]
  
p̃' = T p̃   performs   p' = R p + t   in one multiplication
Chaining: T_total = T₂ @ T₁   (apply T₁ first, then T₂)

Step-by-Step Calculation

Rotate 90° about the z-axis, then translate by [2, 0, 1], applied to the point [1, 0, 0]:

R (90° about z) = [ 0 -1  0 ]      t = [ 2 ]
                  [ 1  0  0 ]          [ 0 ]
                  [ 0  0  1 ]          [ 1 ]

T = [ 0 -1  0  2 ]        p̃ = [ 1 ]
    [ 1  0  0  0 ]            [ 0 ]
    [ 0  0  1  1 ]            [ 0 ]
    [ 0  0  0  1 ]            [ 1 ]
  
T p̃ = [ 0·1 + (-1)·0 + 0·0 + 2·1 ]   [ 2 ]
      [ 1·1 +   0·0  + 0·0 + 0·1 ] = [ 1 ]
      [ 0·1 +   0·0  + 1·0 + 1·1 ]   [ 1 ]
      [ 0·1 +   0·0  + 0·0 + 1·1 ]   [ 1 ]
  
→ p' = [2, 1, 1]: the point rotated onto the y-axis, then shifted by t

Coordinate Frames

Every 3D system juggles several frames: the world frame (fixed reference), body frame (attached to a moving object like a drone), and camera frame (attached to a camera, z pointing forward). A transform T_world_body converts body-frame points to world-frame points; its inverse goes the other way.

import numpy as np
def make_transform(R, t):
    """Build a 4x4 rigid transform from rotation R and translation t."""
    R, t = np.asarray(R, dtype=float), np.asarray(t, dtype=float)
    if R.shape != (3, 3):
        raise ValueError(f"R must be 3x3, got {R.shape}")
    # Guard against non-rigid matrices (scale/shear corrupt everything downstream)
    if not np.allclose(R @ R.T, np.eye(3), atol=1e-5):
        raise ValueError("R is not orthonormal - not a pure rotation")
    T = np.eye(4)
    T[:3, :3] = R
    T[:3, 3] = t
    return T
  
def invert_rigid(T):
    """Closed-form inverse: [R, t]^-1 = [R^T, -R^T t].
    Exact and ~3x faster than np.linalg.inv, with no numerical fill-in."""
    R, t = T[:3, :3], T[:3, 3]
    T_inv = np.eye(4)
    T_inv[:3, :3] = R.T
    T_inv[:3, 3] = -R.T @ t
    return T_inv
  
def transform_points(T, points):
    """points: [N, 3] -> [N, 3]. Appends the homogeneous 1 internally."""
    points = np.atleast_2d(points)
    homogeneous = np.hstack([points, np.ones((len(points), 1))])
    return (T @ homogeneous.T).T[:, :3]
  
# Drone example: body frame -> world frame
T_world_body = make_transform(R_attitude, drone_position)
lidar_hits_world = transform_points(T_world_body, lidar_hits_body)

How It Is Used in ML

• Camera extrinsics: the view_matrix in the Gaussian splatting projection code (Chapter 12) is exactly a world-to-camera rigid transform, and NeRF datasets ship one camera-to-world 4×4 matrix per training image.
• Robot and drone world models: a quadrotor’s pose at any instant is a time-varying T ∈ SE(3). World models like SkyJEPA decode their latent state into exactly this interpretable pose (position + attitude) via a physics prober, and model-based controllers score candidate trajectories as sequences of these transforms.
• Point cloud registration: ICP and learned registration networks output a single rigid T aligning two scans.
• Edge cases: (1) Repeatedly composing transforms accumulates floating-point drift, so R slowly stops being orthonormal — re-orthonormalize periodically with SVD (U, _, Vt = svd(R); R = U @ Vt). (2) Directions (normals, ray directions) use w = 0, not w = 1: they must rotate but never translate. Transforming a normal with w = 1 is a classic silent bug. (3) Use the closed-form inverse above, not a general matrix inverse, to preserve rigidity.

Chapter 6: Rotation Representations — Quaternions

Real-World Analogy: Think of a potato on a skewer. Any orientation of the potato can be reached by choosing the direction the skewer points (an axis) and twisting around it (an angle). A quaternion stores exactly that skewer-plus-twist recipe in four numbers. Euler angles are instead like a camera tripod with three stacked dials (pan, tilt, roll) — and when the camera points straight up, two dials suddenly do the same thing and you lose a direction of control. That jam is gimbal lock, and it is why 3D vision prefers quaternions.

Definition

A unit quaternion encodes a rotation of angle θ around unit axis a:

q = [w, x, y, z] = [cos(θ/2), sin(θ/2)·aₓ, sin(θ/2)·a_y, sin(θ/2)·a_z]
Constraint: w² + x² + y² + z² = 1   (unit norm)

Four numbers with one constraint match the three degrees of freedom of 3D rotation, with no jamming configurations and a smooth optimization landscape.

Step-by-Step Calculation

Encode the 90° z-axis rotation from Chapter 5 as a quaternion, then convert it back to a matrix:

Axis a = [0, 0, 1], angle θ = 90°

q = [cos(45°), 0, 0, sin(45°)] = [0.7071, 0, 0, 0.7071]
  
Conversion formula:
R = [ 1-2(y²+z²)   2(xy-wz)     2(xz+wy)  ]
    [ 2(xy+wz)     1-2(x²+z²)   2(yz-wx)  ]
    [ 2(xz-wy)     2(yz+wx)     1-2(x²+y²)]
  
Plugging in w = z = 0.7071, x = y = 0:
R = [ 1-2(0.5)   -2(0.5)   0 ]   [ 0 -1  0 ]
    [ 2(0.5)     1-2(0.5)  0 ] = [ 1  0  0 ]
    [ 0          0         1 ]   [ 0  0  1 ]
  
→ Exactly the 90° z-rotation matrix from Chapter 5

def quaternion_to_matrix(q, eps=1e-8):
    """q = [w, x, y, z] -> 3x3 rotation matrix."""
    q = np.asarray(q, dtype=float)
    norm = np.linalg.norm(q)
    if norm < eps:
        raise ValueError("Zero-norm quaternion has no defined rotation")
    w, x, y, z = q / norm  # Always renormalize: gradient steps break unit norm
    return np.array([
        [1 - 2*(y*y + z*z),     2*(x*y - w*z),     2*(x*z + w*y)],
        [    2*(x*y + w*z), 1 - 2*(x*x + z*z),     2*(y*z - w*x)],
        [    2*(x*z - w*y),     2*(y*z + w*x), 1 - 2*(x*x + y*y)]
    ])
  
def quaternion_multiply(q1, q2):
    """Compose rotations: apply q2 first, then q1."""
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    return np.array([
        w1*w2 - x1*x2 - y1*y2 - z1*z2,
        w1*x2 + x1*w2 + y1*z2 - z1*y2,
        w1*y2 - x1*z2 + y1*w2 + z1*x2,
        w1*z2 + x1*y2 - y1*x2 + z1*w2
    ])

How It Is Used in ML

• 3D Gaussian splatting (Chapter 12) stores every Gaussian’s orientation as a quaternion precisely because it is the friendliest rotation parameterization for gradient descent: 4 unconstrained numbers that are renormalized after each optimizer step, versus 9 matrix entries that would need 6 orthogonality constraints.
• Drone and robot state: attitude in flight controllers and world models (SkyJEPA’s interpretable state includes the quadrotor’s orientation) is almost universally a quaternion, avoiding gimbal lock during aggressive maneuvers where Euler angles fail.
• Pose estimation networks regress quaternions for object and camera orientation.
• Edge cases: (1) Double cover — q and -q encode the same rotation, so a naive L2 loss between quaternions can report a huge error for identical orientations; use min(||q - q'||, ||q + q'||) or the distance 1 - |q·q'|. (2) Always renormalize after optimization updates, or the "rotation" silently scales geometry. (3) Interpolate with slerp (spherical interpolation), not linear blending—linear interpolation cuts through the sphere and speeds up mid-way. (4) A near-zero quaternion (e.g., from an untrained network head) must be caught before normalization divides by ~0.

Chapter 7: The Pinhole Camera Model and Ray Generation

Real-World Analogy: A shoebox with a pinhole pricked in one end is a working camera: light from the scene passes through the hole and paints an upside-down image on the far wall. The depth of the box is the focal length — a deeper box spreads the same rays over a bigger picture, zooming in. Every camera equation in 3D vision is just bookkeeping for this shoebox: where the hole is, how deep the box is, and where the box stands in the room.

Definition

The pinhole model splits a camera into intrinsics K (properties of the box: focal length and where the optical axis hits the sensor) and extrinsics [R|t] (where the box is in the world — a rigid transform from Chapter 5):

Projection of world point P to pixel (u, v):
p̃ ∝ K [R | t] P̃
  
K = [ fx  0  cx ]    fx, fy: focal length in pixels
    [ 0  fy  cy ]    cx, cy: principal point (image center)
    [ 0   0   1 ]
  
In camera coordinates (X, Y, Z):
u = fx · X/Z + cx
v = fy · Y/Z + cy       (the divide by Z is the perspective effect)

Step-by-Step Calculation

Camera: fx = fy = 500 pixels, cx = cy = 320 (a 640×640 sensor)
Point in camera frame: [0.2, -0.1, 2.0] meters
  
u = 500 × 0.2 / 2.0 + 320 = 50 + 320 = 370
v = 500 × (-0.1) / 2.0 + 320 = -25 + 320 = 295
  
→ Pixel (370, 295). Move the point twice as far (Z = 4) and it lands
  at (345, 307.5)-closer to center, exactly how distant objects shrink.

Projection and Ray Generation Code

def project_point(P_world, K, T_world_to_cam):
    """World point -> pixel coordinates."""
    P_cam = transform_points(T_world_to_cam, P_world)[0]
    if P_cam[2] <= 1e-6:
        return None  # Behind the camera: projecting would flip/explode
    u = K[0, 0] * P_cam[0] / P_cam[2] + K[0, 2]
    v = K[1, 1] * P_cam[1] / P_cam[2] + K[1, 2]
    return np.array([u, v])
  
def generate_rays(K, T_cam_to_world, height, width):
    """One ray per pixel: the exact inverse of projection.
    This is the generate_rays used by every NeRF training loop."""
    u, v = np.meshgrid(np.arange(width), np.arange(height))
    # Unproject pixels to directions in the camera frame (Z = 1 plane)
    dirs_cam = np.stack([
        (u - K[0, 2]) / K[0, 0],
        (v - K[1, 2]) / K[1, 1],
        np.ones_like(u)
    ], axis=-1)
    # Rotate directions to world frame (w = 0: rotate, never translate)
    R = T_cam_to_world[:3, :3]
    ray_dirs = dirs_cam @ R.T
    ray_dirs /= np.linalg.norm(ray_dirs, axis=-1, keepdims=True)
    # All rays start at the camera center
    ray_origins = np.broadcast_to(T_cam_to_world[:3, 3], ray_dirs.shape)
    return ray_origins, ray_dirs

How It Is Used in ML

• NeRF training (Chapter 10): every batch starts by calling exactly this generate_rays—pixels become rays, rays get sampled, samples get rendered, and the loss compares rendered color to the pixel color.
• Multi-view stereo (Chapter 9): the K₁[R₁|t₁] and K₂[R₂|t₂] in the epipolar constraint are these same intrinsics and extrinsics; calibration (finding K) is the prerequisite for any metric reconstruction.
• Gaussian splatting (Chapter 12): the perspective divide is nonlinear, so covariances cannot be projected exactly — the projection Jacobian J in the splatting code is the local linearization of this camera model.
• Edge cases: (1) Points at or behind Z = 0 must be culled before the divide — they otherwise project to infinity or flip across the image. (2) Real lenses distort: apply undistortion before using pinhole math, or reconstructions warp. (3) Convention mismatches are the top source of silent bugs: OpenCV puts +y down and +z forward, OpenGL +y up and -z forward; mixing them flips images and breaks ray directions. (4) Resizing an image scales fx, fy, cx, cy proportionally — forgetting this breaks depth scale.

Chapter 8: Signed Distance Functions (SDFs)

Real-World Analogy: Imagine every point in space knows its distance to the nearest surface, with negative values inside objects and positive outside — like a 3D weather map showing distance to the storm front.

SDF Definition

Mathematical Formulation:

SDF(x) = {
    -d(x, ∂Ω)  if x inside object Ω
     0         if x on surface ∂Ω
    +d(x, ∂Ω)  if x outside object
}
  
where d(x, ∂Ω) = min ||x - y|| for y ∈ ∂Ω

Computing SDF for Simple Shapes:

def sphere_sdf(point, center, radius):
    # Distance to sphere surface
    return np.linalg.norm(point - center) - radius
    # < 0: inside, = 0: on surface, > 0: outside
  
def box_sdf(point, corner_min, corner_max):
    center = (corner_min + corner_max) / 2
    half_size = (corner_max - corner_min) / 2
    # Distance to box surface
    q = np.abs(point - center) - half_size
    return np.linalg.norm(np.maximum(q, 0)) + min(np.max(q), 0)

Neural SDF (DeepSDF)

Learning implicit surfaces:

class NeuralSDF(nn.Module):
    def __init__(self, latent_dim=256):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(3 + latent_dim, 512),
            nn.ReLU(),
            nn.Linear(512, 512) + Skip(),  # Skip connection
            nn.ReLU(),
            nn.Linear(512, 512),
            nn.ReLU(),
            nn.Linear(512, 512) + Skip(),
            nn.ReLU(),
            nn.Linear(512, 1)  # Output SDF value
        )
    def forward(self, points, latent_code):
        # points: [B, N, 3]
        # latent_code: [B, 256] shape encoding
        N = points.shape[1]
        latent_expanded = latent_code.unsqueeze(1).expand(-1, N, -1)
        features = torch.cat([points, latent_expanded], dim=-1)
        sdf_values = self.net(features)
        return sdf_values

Marching Cubes — SDF to Mesh

Extracting the zero level set:

def marching_cubes_simple(sdf_grid, threshold=0):
    """
    Convert SDF grid to mesh
    sdf_grid: [D, H, W] signed distances
    """
    vertices = []
    faces = []
    # Check each voxel
    for i in range(D-1):
        for j in range(H-1):
            for k in range(W-1):
                # Get 8 corner values
                corners = [
                    sdf_grid[i,j,k],     sdf_grid[i+1,j,k],
                    sdf_grid[i+1,j+1,k], sdf_grid[i,j+1,k],
                    sdf_grid[i,j,k+1],   sdf_grid[i+1,j,k+1],
                    sdf_grid[i+1,j+1,k+1], sdf_grid[i,j+1,k+1]
                ]
                # Determine cube configuration (256 cases)
                cube_index = sum(2**i if c > threshold else 0 
                               for i, c in enumerate(corners))
                # Look up triangulation from table
                triangulation = MARCHING_CUBES_TABLE[cube_index]
                # Add vertices and faces...
                # (Complex implementation details omitted)
    return vertices, faces

Chapter 9: Multi-View Stereo (MVS)

Real-World Analogy: Like how your two eyes create depth perception — multiple camera views of the same scene allow 3D reconstruction by finding corresponding points.

Epipolar Geometry

Fundamental Constraint:

For point P seen in cameras C1 and C2:
p1 = K1 × [R1|t1] × P  (projection in camera 1)
p2 = K2 × [R2|t2] × P  (projection in camera 2)
  
Epipolar constraint: p2^T × F × p1 = 0
where F = fundamental matrix

Depth from Stereo Matching

Computing depth from disparity:

def stereo_depth(left_img, right_img, baseline=0.1, focal=500):
    """
    baseline: distance between cameras (meters)
    focal: focal length (pixels)
    """
    # 1. Rectify images (align epipolar lines)
    left_rect, right_rect = rectify_stereo_pair(left_img, right_img)
    # 2. Compute disparity (pixel offset between matches)
    disparity_map = np.zeros_like(left_rect[:,:,0])
    for y in range(height):
        for x in range(width):
            # Search along epipolar line (same y)
            patch_left = left_rect[y-w:y+w+1, x-w:x+w+1]
            min_ssd = float('inf')
            best_disp = 0
            # Search in right image
            for x_right in range(max(0, x-max_disp), min(width, x+1)):
                patch_right = right_rect[y-w:y+w+1, x_right-w:x_right+w+1]
                ssd = np.sum((patch_left - patch_right)**2)
                if ssd < min_ssd:
                    min_ssd = ssd
                    best_disp = x - x_right
            disparity_map[y, x] = best_disp
    # 3. Convert disparity to depth
    # depth = baseline × focal / disparity
    depth_map = np.where(disparity_map > 0, 
                         baseline * focal / disparity_map,
                         0)
    return depth_map

Practical numbers:

Baseline = 10cm, Focal = 500 pixels
Disparity = 50 pixels
→ Depth = 0.1 × 500 / 50 = 1 meter
  
Disparity = 5 pixels  
→ Depth = 0.1 × 500 / 5 = 10 meters

Part III: Neural Rendering

Chapter 10: Neural Radiance Fields (NeRF)

Real-World Analogy: Instead of storing a 3D model, train a neural network that acts like a “3D photograph” — ask it what color/density exists at any 3D point from any viewing angle.

Volume Rendering Equation

Mathematical Foundation:

Color at pixel = ∫ T(t) × σ(r(t)) × c(r(t), d) dt

where:
- r(t) = o + t×d  (ray from origin o in direction d)
- T(t) = exp(-∫σ(s)ds)  (transmittance)
- σ = density
- c = color
  
Discretized:
C = Σᵢ Tᵢ × αᵢ × cᵢ
where:
- αᵢ = 1 - exp(-σᵢ × δᵢ)  (opacity)
- Tᵢ = Πⱼ<ᵢ (1 - αⱼ)  (accumulated transmittance)

Positional Encoding

Why it’s crucial:

def positional_encoding(x, L=10):
    """
    Map x to higher dimensional space
    γ(x) = [sin(2^0×π×x), cos(2^0×π×x), ..., sin(2^L×π×x), cos(2^L×π×x)]
    """
    freqs = 2**torch.arange(L) * np.pi
    # x: [..., 3]
    # freqs: [L]
    x_freq = x[..., None, :] * freqs[:, None]  # [..., L, 3]
    # Compute sin and cos
    encoding = torch.cat([
        x,  # Original input
        torch.sin(x_freq).flatten(-2),  # Sin components
        torch.cos(x_freq).flatten(-2)   # Cos components
    ], dim=-1)
    # Output dimension: 3 + 2×L×3 = 3 + 60 = 63 for L=10
    return encoding

Why it works: Neural networks struggle with high-frequency details. Positional encoding maps low-dimensional inputs to high-dimensional space where linear layers can learn high frequencies.

Hierarchical Sampling

Coarse-to-Fine Strategy:

def hierarchical_sample(ray_o, ray_d, near=2., far=6., n_coarse=64, n_fine=128):
    # 1. Coarse sampling - uniform
    t_coarse = torch.linspace(near, far, n_coarse)
    points_coarse = ray_o + ray_d * t_coarse[:, None]
    # 2. Query network for density
    density_coarse = nerf_coarse(points_coarse, ray_d)
    # 3. Compute weights (importance)
    weights = compute_volume_rendering_weights(density_coarse, t_coarse)
    # 4. Sample more points where density is high
    # Inverse transform sampling from weight distribution
    t_fine = sample_pdf(t_coarse, weights, n_fine)
    # 5. Combine and sort all samples
    t_combined = torch.sort(torch.cat([t_coarse, t_fine]))[0]
    points_combined = ray_o + ray_d * t_combined[:, None]
    return points_combined, t_combined

Performance Impact:

• Uniform sampling only: 256 samples needed
• Hierarchical: 64 + 128 = 192 samples
• Quality improvement: 2–3 dB PSNR

Chapter 11: Spherical Harmonics for View-Dependent Color

Real-World Analogy: Like a graphic equalizer, but for directions instead of sound frequencies. The first slider sets the overall volume (the object’s base color, the same from every viewpoint). The next sliders add gentle variation — brighter from the left, darker from below. Higher sliders add sharper and sharper directional detail, like a glossy highlight that only appears when you look from one narrow angle. A few sliders reproduce most real materials; only a mirror needs infinitely many.

Definition

Spherical harmonics (SH) are an orthogonal set of basis functions defined on the sphere — the directional analogue of Fourier series. Any function of viewing direction, such as the color an object shows from each angle, is approximated as a weighted sum:

color(d) = Σ_ℓ Σ_m  c_ℓm · Y_ℓm(d)      d = unit view direction

Degree ℓ band has 2ℓ+1 functions:
ℓ = 0: 1 function  (constant - base color)
ℓ = 1: 3 functions (linear variation across the sphere)
ℓ = 2: 5 functions (quadratic lobes)
ℓ = 3: 7 functions (sharper lobes)
Total through degree 3: 16 functions per color channel

The first bands, for d = (x, y, z):

Y₀₀   = 0.2821
Y₁,₋₁ = 0.4886·y     Y₁,₀ = 0.4886·z     Y₁,₁ = 0.4886·x
Y₂ band: 5 quadratic terms like 1.0925·xy, 0.3154·(3z²-1), ...

Step-by-Step Calculation

A degree-1 red channel with coefficients c₀₀ = 1.0, c₁,₋₁ = 0, c₁,₀ = 0.4, c₁,₁ = 0, viewed from straight above (d = [0, 0, 1]):

red(d) = 1.0 × 0.2821 + 0.4 × 0.4886 × z
       = 0.2821 + 0.4 × 0.4886 × 1.0 = 0.4775

Viewed from below (d = [0, 0, -1]):
red(d) = 0.2821 - 0.1954 = 0.0867
  
→ One extra coefficient made the object bright from above
  and dark from below, at the cost of a single multiply-add.

def evaluate_sh_degree1(sh_coeffs, view_dir, eps=1e-8):
    """sh_coeffs: [4, 3] (4 basis functions x RGB). view_dir: [3]."""
    norm = np.linalg.norm(view_dir)
    if norm < eps:
        raise ValueError("View direction must be non-zero")
    x, y, z = view_dir / norm  # SH assumes unit directions
    basis = np.array([0.2821, 0.4886 * y, 0.4886 * z, 0.4886 * x])
    color = basis @ sh_coeffs
    # Coefficients are unconstrained during optimization: clamp output
    return np.clip(color, 0.0, 1.0)

How It Is Used in ML

• 3D Gaussian splatting (Chapter 12): every Gaussian’s sh_coefficients are exactly this—48 numbers (16 basis functions × RGB at degree 3) that make color view-dependent. Evaluating a 16-term polynomial per Gaussian is thousands of times cheaper than querying an MLP per sample, which is a large part of why splatting renders at 130+ FPS while NeRF takes seconds. Those 48 of ~59 floats per Gaussian also dominate the 500 MB memory figure in the performance table.
• Plenoxels and baked NeRF variants store SH coefficients on voxel grids to remove the neural network entirely.
• Graphics pipelines use SH for environment lighting and ambient occlusion — precomputing how a scene’s illumination varies by direction.
• Edge cases: (1) Low degrees cannot represent mirror-sharp reflections — highlights smear or “ring” (oscillating halos), the same Gibbs ringing as truncated Fourier series. (2) Coefficients are unconstrained, so evaluated colors can go negative or above 1 at extreme angles — clamp them. (3) Directions must be normalized before evaluation. (4) Degree 3 is the practical sweet spot: degree 4+ quadruples memory for marginal quality on typical scenes.

Chapter 12: 3D Gaussian Splatting

Real-World Analogy: Instead of pixels or voxels, represent the scene as millions of semi-transparent 3D Gaussian “blobs” that can be quickly projected and blended to create images from any viewpoint.

3D Gaussian Representation

Each Gaussian has:

class Gaussian3D:
    def __init__(self):
        # Position in 3D
        self.position = np.array([x, y, z])  # μ ∈ ℝ³
        # Covariance (defines ellipsoid shape)
        # Factorized as RS S^T R^T for optimization
        self.rotation = np.array([qw, qx, qy, qz])  # Quaternion
        self.scale = np.array([sx, sy, sz])  # 3D scale
        # Appearance  
        self.opacity = α  # ∈ [0, 1]
        self.sh_coefficients = np.array([...])  # Spherical harmonics for view-dependent color
        # Computed covariance
        R = quaternion_to_matrix(self.rotation)
        S = np.diag(self.scale)
        self.covariance = R @ S @ S.T @ R.T  # Σ ∈ ℝ³ˣ³

Gaussian Projection to 2D

World to screen transformation:

def project_gaussian(gaussian_3d, view_matrix, projection_matrix):
    """
    Project 3D Gaussian to 2D screen space
    """
    # 1. Transform center to camera space
    p_cam = view_matrix @ [*gaussian_3d.position, 1]
    # 2. Compute Jacobian of projection
    J = compute_projection_jacobian(p_cam, projection_matrix)
    # 3. Transform covariance to 2D
    # Σ' = J W Σ W^T J^T
    W = view_matrix[:3, :3]  # Rotation part
    Σ_2d = J @ W @ gaussian_3d.covariance @ W.T @ J.T
    # 4. Project center to screen
    p_screen = projection_matrix @ p_cam
    p_screen = p_screen[:2] / p_screen[2]  # Perspective divide
    return {
        'center_2d': p_screen,
        'covariance_2d': Σ_2d,
        'opacity': gaussian_3d.opacity,
        'color': evaluate_sh(gaussian_3d.sh_coefficients, view_dir),
        'depth': p_cam[2]
    }

Efficient Rasterization

Tile-based rendering:

def rasterize_gaussians(gaussians_2d, image_width, image_height, tile_size=16):
    """
    Fast GPU-friendly rasterization
    """
    # 1. Divide screen into tiles
    n_tiles_x = (image_width + tile_size - 1) // tile_size
    n_tiles_y = (image_height + tile_size - 1) // tile_size
    # 2. Assign Gaussians to tiles
    tiles = [[[] for _ in range(n_tiles_x)] for _ in range(n_tiles_y)]
    for g in gaussians_2d:
        # Compute 99% confidence bounds
        radius = 3.0 * np.sqrt(np.max(np.linalg.eigvals(g['covariance_2d'])))
        # Find affected tiles
        min_tile_x = int((g['center_2d'][0] - radius) // tile_size)
        max_tile_x = int((g['center_2d'][0] + radius) // tile_size)
        min_tile_y = int((g['center_2d'][1] - radius) // tile_size)
        max_tile_y = int((g['center_2d'][1] + radius) // tile_size)
        # Add to tiles
        for ty in range(max(0, min_tile_y), min(n_tiles_y, max_tile_y + 1)):
            for tx in range(max(0, min_tile_x), min(n_tiles_x, max_tile_x + 1)):
                tiles[ty][tx].append(g)
    # 3. Render each pixel
    image = np.zeros((image_height, image_width, 3))
    for y in range(image_height):
        for x in range(image_width):
            tile_y, tile_x = y // tile_size, x // tile_size
            # Sort Gaussians in tile by depth
            gaussians = sorted(tiles[tile_y][tile_x], key=lambda g: g['depth'])
            # Alpha blending
            color = np.zeros(3)
            accumulated_alpha = 0
            for g in gaussians:
                if accumulated_alpha > 0.99:
                    break
                # Evaluate Gaussian at pixel
                diff = np.array([x, y]) - g['center_2d']
                Σ_inv = np.linalg.inv(g['covariance_2d'])
                power = -0.5 * diff @ Σ_inv @ diff
                if power > -5:  # Threshold for numerical stability
                    alpha = g['opacity'] * np.exp(power)
                    weight = alpha * (1 - accumulated_alpha)
                    color += weight * g['color']
                    accumulated_alpha += weight
            image[y, x] = color
    return image

Adaptive Density Control

Split and prune Gaussians:

def adaptive_control(gaussians, gradients, τ_pos=0.002, τ_scale=0.05):
    """
    Split large Gaussians, remove small ones
    """
    new_gaussians = []
    for i, g in enumerate(gaussians):
        grad_pos = gradients['position'][i]
        # Split if gradient too large (under-reconstruction)
        if np.linalg.norm(grad_pos) > τ_pos:
            # Create two smaller Gaussians
            for offset in [-1, 1]:
                g_new = copy(g)
                # Offset along largest scale axis
                axis = np.argmax(g.scale)
                offset_vec = np.zeros(3)
                offset_vec[axis] = offset * g.scale[axis] * 0.5
                g_new.position = g.position + offset_vec
                g_new.scale = g.scale * 0.8  # Reduce scale
                new_gaussians.append(g_new)
        # Clone if scale too large
        elif np.max(g.scale) > τ_scale:
            g_new = copy(g)
            g_new.scale = g.scale * 0.7
            new_gaussians.append(g)
            new_gaussians.append(g_new)
        # Prune if opacity too low
        elif g.opacity < 0.01:
            continue  # Remove
        else:
            new_gaussians.append(g)
    return new_gaussians

Performance Comparison

Rendering Speed (1920×1080):

NeRF:              ~30 seconds per frame
Instant-NGP:       ~50 ms per frame
3D Gaussian Splatting: ~5 ms per frame (>130 FPS)
  
Training time (300K iterations):
NeRF:              ~24 hours
Instant-NGP:       ~10 minutes  
3D Gaussian Splatting: ~30 minutes
  
Memory:
NeRF:              ~5 MB (network weights)
3D Gaussian Splatting: ~500 MB (2M Gaussians × 250 bytes each)

Chapter 13: Instant Neural Graphics (Instant-NGP)

Real-World Analogy: Like a filing system that uses hash tables instead of regular folders — randomly distributing information but with instant lookup, allowing massive resolution without massive memory.

Multi-Resolution Hash Encoding

Core Innovation:

class HashGrid:
    def __init__(self, resolution, n_features=2, table_size=2**19):
        self.resolution = resolution
        self.table_size = table_size
        # Hash table of learnable features
        self.table = nn.Parameter(torch.randn(table_size, n_features))
    def hash_function(self, coords):
        """
        Spatial hash: map integer coords to table index
        Using prime numbers for better distribution
        """
        primes = [1, 2654435761, 805459861]
        result = 0
        for i, c in enumerate(coords):
            result ^= c * primes[i % 3]
        return result % self.table_size
    def forward(self, positions):
        # positions: [..., 3] in [0, 1]³
        # Scale to grid
        scaled = positions * self.resolution
        # Get integer coordinates of voxel corners
        corners_min = torch.floor(scaled).long()
        corners_max = corners_min + 1
        # Interpolation weights
        weights = scaled - corners_min
        # Hash all 8 corners
        features = 0
        for dx in [0, 1]:
            for dy in [0, 1]:
                for dz in [0, 1]:
                    corner = corners_min + torch.tensor([dx, dy, dz])
                    hash_idx = self.hash_function(corner)
                    # Trilinear weight
                    w = ((1-weights[0]) if dx==0 else weights[0]) * \
                        ((1-weights[1]) if dy==0 else weights[1]) * \
                        ((1-weights[2]) if dz==0 else weights[2])
                    features += w * self.table[hash_idx]
        return features

Multi-Level Feature Grid

class MultiResHashEncoding:
    def __init__(self, n_levels=16, base_resolution=16, max_resolution=2048):
        self.n_levels = n_levels
        # Geometric progression of resolutions
        b = np.exp((np.log(max_resolution) - np.log(base_resolution)) / (n_levels-1))
        self.grids = []
        for level in range(n_levels):
            resolution = int(base_resolution * b**level)
            self.grids.append(HashGrid(resolution))
    def forward(self, positions):
        # Concatenate features from all levels
        features = []
        for grid in self.grids:
            features.append(grid(positions))
        return torch.cat(features, dim=-1)
        # Output: n_levels × n_features per level = 16 × 2 = 32 dims

Memory vs Resolution:

Traditional voxel grid at 2048³:
  Memory = 2048³ × 4 bytes = 32 GB

Hash encoding with 2^19 entries:
  Memory = 2^19 × 8 bytes = 4 MB
  Effective resolution: Up to 2048³ with collisions

Chapter 14: Differentiable Rendering

Real-World Analogy: Like having a camera that not only takes photos but can tell you exactly how to move objects to make the photo look more like a target image — enabling optimization through the rendering process.

Soft Rasterization

Making rendering differentiable:

def soft_rasterize(vertices, faces, image_size=256, sigma=1e-4):
    """
    Differentiable rasterization using soft assignments
    """
    # For each pixel, compute distance to all faces
    pixels = []
    for y in range(image_size):
        for x in range(image_size):
            pixel_pos = np.array([x, y]) / image_size
            min_dist = float('inf')
            closest_face = None
            for face_idx, face in enumerate(faces):
                # Get triangle vertices
                v0, v1, v2 = vertices[face]
                # Compute barycentric coordinates
                bary = compute_barycentric(pixel_pos, v0[:2], v1[:2], v2[:2])
                # Check if inside triangle
                if all(b >= 0 for b in bary):
                    # Distance to triangle plane
                    z = bary[0]*v0[2] + bary[1]*v1[2] + bary[2]*v2[2]
                    if z < min_dist:
                        min_dist = z
                        closest_face = face_idx
                else:
                    # Distance to triangle edges
                    edge_dist = distance_to_triangle_edge(pixel_pos, v0[:2], v1[:2], v2[:2])
                    # Soft assignment based on distance
                    prob = torch.exp(-edge_dist**2 / (2*sigma**2))
                    # Aggregate probabilities for antialiasing
            pixels.append(rendered_value)
    return torch.stack(pixels).reshape(image_size, image_size, -1)

Gradient Flow Through Rendering

def rendering_gradient(loss, vertices):
    """
    How vertex positions affect rendered image
    """
    # Automatic differentiation gives us:
    # ∂L/∂v = ∂L/∂I × ∂I/∂p × ∂p/∂v
    # where:
    #   L = loss
    #   I = rendered image
    #   p = projected 2D positions
    #   v = 3D vertices
    grad_vertices = torch.autograd.grad(loss, vertices)[0]
    # Gradient indicates how to move vertices to minimize loss
    return grad_vertices

Part IV: Advanced Topics

Chapter 15: Neural Surface Reconstruction

Real-World Analogy: Like a detective reconstructing a crime scene from photos — using multiple clues (views, lighting, shadows) to figure out the exact 3D shape.

NeuS — Neural Implicit Surfaces

Volume rendering with SDF:

def neus_rendering(sdf_net, ray_o, ray_d, n_samples=64):
    """
    Convert SDF to density for volume rendering
    """
    # Sample points along ray
    t = torch.linspace(0, 1, n_samples)
    points = ray_o + ray_d * t[:, None]
    # Query SDF
    sdf_values = sdf_net(points)
    # Convert SDF to density using logistic distribution
    # φ_s(x) = s × Φ_s(x) × (1 - Φ_s(x))
    # where Φ_s(x) = sigmoid(s × x)
    s = 100.0  # Sharpness parameter
    density = s * torch.sigmoid(s * sdf_values) * (1 - torch.sigmoid(s * sdf_values))
    # Volume rendering
    delta = t[1] - t[0]
    alpha = 1 - torch.exp(-density * delta)
    transmittance = torch.cumprod(1 - alpha[:-1], dim=0)
    transmittance = torch.cat([torch.ones(1), transmittance])
    weights = transmittance * alpha
    # Render color
    colors = sdf_net.get_color(points, ray_d)
    rendered_color = torch.sum(weights[:, None] * colors, dim=0)
    return rendered_color, weights

DMTet — Deep Marching Tetrahedra

Hybrid explicit-implicit representation:

class DMTet:
    def __init__(self, grid_res=128):
        # Tetrahedral grid
        self.vertices = create_tet_grid(grid_res)  # [N_v, 3]
        self.tets = create_tet_topology(grid_res)  # [N_t, 4]
        # Learnable SDF values at vertices
        self.sdf_values = nn.Parameter(torch.zeros(len(self.vertices)))
        # Learnable vertex offsets
        self.vertex_offsets = nn.Parameter(torch.zeros_like(self.vertices))
    def extract_mesh(self):
        # Deform vertices
        deformed_vertices = self.vertices + self.vertex_offsets
        # Extract surface where SDF = 0
        triangles = []
        for tet in self.tets:
            v_sdf = self.sdf_values[tet]  # 4 SDF values
            # Check sign changes (surface crosses tet)
            pos_mask = v_sdf > 0
            if pos_mask.sum() in [1, 2, 3]:  # Surface crosses
                # Marching tetrahedra table lookup
                config = pos_mask.long()
                tri_indices = MARCHING_TET_TABLE[config]
                # Interpolate vertices on edges
                for tri in tri_indices:
                    triangle_vertices = []
                    for edge in tri:
                        v1, v2 = tet[edge[0]], tet[edge[1]]
                        s1, s2 = v_sdf[edge[0]], v_sdf[edge[1]]
                        # Linear interpolation
                        t = s1 / (s1 - s2)
                        v_interp = (1-t) * deformed_vertices[v1] + t * deformed_vertices[v2]
                        triangle_vertices.append(v_interp)
                    triangles.append(triangle_vertices)
        return triangles

Chapter 16: 3D Loss Functions — Chamfer Distance and Eikonal Regularization

Real-World Analogy: Two rival scout troops are scattered across a park and you want to measure how similarly they are spread out. Ask every scout in the red troop to walk to the nearest blue scout and write down the distance, then ask every blue scout to walk to the nearest red scout. The average of all those walks is the Chamfer distance — small only when the two troops occupy the park the same way. The Eikonal constraint is a different kind of check: on a proper trail-distance map, walking one meter straight away from the trail must raise the number on the map by exactly one meter. Any map where the number climbs faster or slower than you walk is lying about distance.

Chamfer Distance

Definition: The standard loss for comparing two point clouds P and Q with no point-to-point correspondence:

CD(P, Q) = (1/|P|) Σ_{p∈P} min_{q∈Q} ||p - q||²
         + (1/|Q|) Σ_{q∈Q} min_{p∈P} ||q - p||²

Both directions are required: one-sided distance is blind to
missing geometry on the other side.

Step-by-Step Calculation:

P = {(0,0,0), (1,0,0)}          Q = {(0,0,0), (2,0,0)}

P → Q:  (0,0,0) → nearest is (0,0,0), dist² = 0
        (1,0,0) → nearest is (0,0,0) or (2,0,0), dist² = 1
        average = (0 + 1) / 2 = 0.5
  
Q → P:  (0,0,0) → nearest is (0,0,0), dist² = 0
        (2,0,0) → nearest is (1,0,0), dist² = 1
        average = (0 + 1) / 2 = 0.5
  
CD(P, Q) = 0.5 + 0.5 = 1.0

from scipy.spatial import cKDTree
def chamfer_distance(P, Q):
    """P: [N, 3], Q: [M, 3]. KD-tree makes this O(N log M) instead of O(N*M)."""
    P, Q = np.atleast_2d(P), np.atleast_2d(Q)
    if len(P) == 0 or len(Q) == 0:
        raise ValueError("Chamfer distance is undefined for empty point sets")
    d_pq, _ = cKDTree(Q).query(P)  # Nearest Q-neighbor for each p
    d_qp, _ = cKDTree(P).query(Q)  # Nearest P-neighbor for each q
    return np.mean(d_pq**2) + np.mean(d_qp**2)

Eikonal Regularization

Definition: A true signed distance function must satisfy the eikonal equation — its gradient has unit length everywhere, because moving one meter perpendicular to the surface changes the distance by one meter:

||∇SDF(x)|| = 1  for all x

Eikonal loss:  L_eik = E_x [ (||∇f(x)|| - 1)² ]

Without this term, a network trained only to be zero on the surface can invent any function that happens to cross zero there — valid-looking surface, garbage distances everywhere else.

def eikonal_loss(sdf_net, points):
    """points: [N, 3] sampled in space (uniform + near-surface)."""
    points = points.clone().requires_grad_(True)
    sdf = sdf_net(points)
    grad = torch.autograd.grad(
        outputs=sdf, inputs=points,
        grad_outputs=torch.ones_like(sdf),
        create_graph=True  # Needed: this loss itself gets backpropagated
    )[0]
    return ((grad.norm(dim=-1) - 1.0) ** 2).mean()

How It Is Used in ML

• Chamfer distance is the default training loss and benchmark metric for point cloud autoencoders, single-image shape prediction, and mesh reconstruction (meshes are compared by sampling points on both surfaces). One-sided variants handle partial scans: a LiDAR sweep should match part of the model, not all of it.
• Eikonal loss is the standard regularizer whenever a network represents geometry as an SDF — DeepSDF-style models, NeuS (Chapter 15), and IGR all add it so the learned field behaves like a real distance function, which sphere-tracing renderers and marching cubes both rely on.
• Edge cases: (1) The squared form of Chamfer lets a single outlier dominate the loss — robust variants use unsquared distance or trim the largest matches. (2) Chamfer is fooled by density mismatch: a cloud with all points bunched on one wall can score deceptively well against a uniformly sampled room; Earth Mover’s Distance fixes this with a one-to-one matching but costs O(N³) versus Chamfer’s O(N log N). (3) Chamfer is not a true metric (no triangle inequality), so don’t use it to build metric data structures. (4) Eikonal weight is a balance: too high oversmooths and erases thin details, too low leaves invalid distances that break sphere tracing with overshoots.

Chapter 17: Dynamic 3D Scenes

Real-World Analogy: Like stop-motion animation in 3D — representing how objects move and deform over time, not just static shapes.

4D Neural Fields

class DynamicNeRF(nn.Module):
    def __init__(self):
        super().__init__()
        # Canonical space network
        self.canonical_net = NeRF()
        # Deformation network
        self.deformation_net = nn.Sequential(
            nn.Linear(3 + 1, 256),  # xyz + time
            nn.ReLU(),
            nn.Linear(256, 256),
            nn.ReLU(),
            nn.Linear(256, 3)  # Displacement
        )
    def forward(self, x, t):
        # x: [N, 3] positions
        # t: scalar time
        # Compute deformation
        t_embed = torch.ones(len(x), 1) * t
        deform_input = torch.cat([x, t_embed], dim=-1)
        displacement = self.deformation_net(deform_input)
        # Warp to canonical space
        x_canonical = x + displacement
        # Query canonical network
        density, color = self.canonical_net(x_canonical)
        return density, color

Gaussian Splatting for Dynamic Scenes

class Dynamic4DGaussians:
    def __init__(self, n_frames=100):
        # Gaussians with time-varying properties
        self.base_positions = []  # Canonical positions
        self.base_rotations = []
        self.base_scales = []
        # Learned deformations per frame
        self.position_deltas = nn.Parameter(torch.zeros(n_frames, n_gaussians, 3))
        self.rotation_deltas = nn.Parameter(torch.zeros(n_frames, n_gaussians, 4))
        self.scale_deltas = nn.Parameter(torch.zeros(n_frames, n_gaussians, 3))
    def get_gaussians_at_time(self, frame_idx):
        gaussians = []
        for i in range(len(self.base_positions)):
            g = Gaussian3D()
            g.position = self.base_positions[i] + self.position_deltas[frame_idx, i]
            g.rotation = quaternion_multiply(
                self.base_rotations[i], 
                self.rotation_deltas[frame_idx, i]
            )
            g.scale = self.base_scales[i] * torch.exp(self.scale_deltas[frame_idx, i])
            gaussians.append(g)
        return gaussians

Chapter 18: Large-Scale Scene Representation

Real-World Analogy: Like Google Earth — you can’t load the entire planet at full resolution, so you use levels of detail and only load what’s visible.

Octree-based NeRF

class OctreeNeRF:
    def __init__(self, max_depth=10):
        self.root = OctreeNode([0, 0, 0], size=1.0)
        self.max_depth = max_depth
        # Feature grid per octree node
        self.features = {}
    def adaptive_sampling(self, ray_o, ray_d):
        """
        Sample more in detailed regions
        """
        samples = []
        # Traverse octree along ray
        intersected_nodes = self.root.ray_intersection(ray_o, ray_d)
        for node, t_enter, t_exit in intersected_nodes:
            # Sample density based on node level
            level = node.get_depth()
            n_samples = 2 ** (self.max_depth - level)
            t_samples = torch.linspace(t_enter, t_exit, n_samples)
            samples.extend(t_samples)
        return torch.sort(torch.tensor(samples))[0]

Block-NeRF for City-Scale

class BlockNeRF:
    def __init__(self, block_size=100):  # 100 meter blocks
        self.blocks = {}
        self.block_size = block_size
    def get_block_id(self, position):
        return tuple(np.floor(position / self.block_size).astype(int))
    def render_ray(self, ray_o, ray_d):
        # Find blocks intersected by ray
        blocks_on_ray = self.get_intersected_blocks(ray_o, ray_d)
        colors = []
        transmittances = []
        for block_id, t_enter, t_exit in blocks_on_ray:
            if block_id not in self.blocks:
                # Skip empty blocks
                continue
            block_nerf = self.blocks[block_id]
            # Transform ray to block coordinates
            block_origin = np.array(block_id) * self.block_size
            ray_o_local = ray_o - block_origin
            # Render within block
            color, trans = block_nerf.render(ray_o_local, ray_d, t_enter, t_exit)
            colors.append(color)
            transmittances.append(trans)
        # Composite across blocks
        return composite_colors(colors, transmittances)

Chapter 19: 3D Generation and Synthesis

Real-World Analogy: Like a 3D printer that works from imagination — generating new 3D objects from text descriptions or single images.

3D Diffusion Models

class DiffusionSDF:
    """
    Generate 3D shapes by denoising SDF grids
    """
    def __init__(self, grid_size=64):
        self.grid_size = grid_size
        self.denoiser = UNet3D(in_channels=1, out_channels=1)
    def forward_diffusion(self, sdf_grid, t):
        """Add noise to SDF"""
        noise = torch.randn_like(sdf_grid)
        alpha = self.noise_schedule(t)
        noisy_sdf = np.sqrt(alpha) * sdf_grid + np.sqrt(1 - alpha) * noise
        return noisy_sdf, noise
    def reverse_diffusion(self, noisy_sdf, t):
        """Predict and remove noise"""
        predicted_noise = self.denoiser(noisy_sdf, t)
        alpha = self.noise_schedule(t)
        denoised = (noisy_sdf - np.sqrt(1 - alpha) * predicted_noise) / np.sqrt(alpha)
        return denoised
    def generate(self, n_steps=1000):
        """Generate SDF from noise"""
        # Start with random noise
        sdf = torch.randn(1, 1, self.grid_size, self.grid_size, self.grid_size)
        # Gradually denoise
        for t in reversed(range(n_steps)):
            sdf = self.reverse_diffusion(sdf, t)
        return sdf

Text-to-3D with DreamFusion

class DreamFusion:
    def __init__(self):
        self.nerf = NeRF()
        self.text_encoder = CLIP()  # Pretrained
        self.diffusion_model = StableDiffusion()  # Pretrained
    def score_distillation_loss(self, rendered_image, text_prompt, t):
        """
        Use 2D diffusion model to guide 3D generation
        """
        # Encode text
        text_features = self.text_encoder(text_prompt)
        # Add noise to rendered image
        noise = torch.randn_like(rendered_image)
        noisy_image = self.add_noise(rendered_image, noise, t)
        # Predict noise with diffusion model
        predicted_noise = self.diffusion_model(noisy_image, t, text_features)
        # Score Distillation Sampling loss
        # Gradient of log-probability
        loss = (predicted_noise - noise).detach() * rendered_image
        return loss.mean()
    def optimize(self, text_prompt, n_iters=10000):
        optimizer = torch.optim.Adam(self.nerf.parameters(), lr=1e-3)
        for i in range(n_iters):
            # Random camera
            camera_pos = sample_camera_on_sphere()
            rays = generate_rays(camera_pos)
            # Render view
            rendered = self.nerf.render(rays)
            # Random timestep
            t = torch.randint(0, 1000, (1,))
            # Compute SDS loss
            loss = self.score_distillation_loss(rendered, text_prompt, t)
            optimizer.zero_grad()
            loss.backward()
            optimizer.step()

Part V: Performance and Optimization

Chapter 20: Acceleration Structures

Real-World Analogy: Like a courier company sorting packages — parcels are grouped by region, then city, then street, so no driver ever scans every address in the country. A BVH wraps triangles in nested boxes so a ray only ever tests the boxes it actually enters.

Bounding Volume Hierarchies (BVH)

class BVHNode:
    def __init__(self, primitives):
        # Compute bounding box
        self.bbox_min = np.min([p.bbox_min for p in primitives], axis=0)
        self.bbox_max = np.max([p.bbox_max for p in primitives], axis=0)
        if len(primitives) <= 4:
            # Leaf node
            self.primitives = primitives
            self.left = None
            self.right = None
        else:
            # Split along longest axis
            axis = np.argmax(self.bbox_max - self.bbox_min)
            primitives.sort(key=lambda p: p.center[axis])
            mid = len(primitives) // 2
            self.left = BVHNode(primitives[:mid])
            self.right = BVHNode(primitives[mid:])
    def intersect_ray(self, ray_o, ray_d):
        # Check bbox intersection
        t_min, t_max = ray_bbox_intersect(ray_o, ray_d, self.bbox_min, self.bbox_max)
        if t_max < 0 or t_min > t_max:
            return None
        if self.is_leaf():
            # Test primitives
            best_hit = None
            for prim in self.primitives:
                hit = prim.intersect(ray_o, ray_d)
                if hit and (not best_hit or hit.t < best_hit.t):
                    best_hit = hit
            return best_hit
        # Recurse
        hit_left = self.left.intersect_ray(ray_o, ray_d) if self.left else None
        hit_right = self.right.intersect_ray(ray_o, ray_d) if self.right else None
        if hit_left and hit_right:
            return hit_left if hit_left.t < hit_right.t else hit_right
        return hit_left or hit_right

GPU Memory Management

def optimize_gaussian_memory():
    """
    Pack Gaussian data for GPU efficiency
    """
    # Structure of Arrays (SoA) instead of Array of Structures (AoS)
    # Bad: Each Gaussian is an object
    gaussians = [Gaussian() for _ in range(1000000)]
    # Random memory access, cache misses
    # Good: Separate arrays for each property
    positions = torch.zeros(1000000, 3)      # Coalesced access
    rotations = torch.zeros(1000000, 4)      # Better cache usage
    scales = torch.zeros(1000000, 3)         # SIMD friendly
    opacities = torch.zeros(1000000, 1)
    sh_coeffs = torch.zeros(1000000, 48)     # (3 color × 16 coefficients)
    # Quantization for memory reduction
    positions_quantized = (positions * 1000).short()  # mm precision
    # 6 bytes instead of 12 per position

Conclusion: Choosing the Right Representation

Decision Tree

When to use what:

Need real-time rendering?
├─ Yes → 3D Gaussian Splatting
│        (130+ FPS, good quality)
└─ No → Need perfect quality?
        ├─ Yes → NeRF variants
        │        (Slow but photorealistic)
        └─ No → Need to edit/manipulate?
                ├─ Yes → Mesh or SDF
                │        (Explicit geometry)
                └─ No → Point clouds
                        (Simple, efficient)
  
Memory constrained?
├─ Yes → Instant-NGP (hash encoding)
└─ No → Storage constrained?
        ├─ Yes → Neural fields (compact)
        └─ No → Gaussian Splatting (fast)
  
Dynamic scenes?
├─ Yes → 4D representations
│        (Deformation fields, dynamic Gaussians)
└─ No → Static methods
  
Large scale?
├─ Yes → Hierarchical (Octree, Block-NeRF)
└─ No → Single model

Performance Summary

| Method                | Training | Rendering    | Memory       | Quality   |
| - - - -               | - - - - -| - - - - - -  | - - - - -    | - - - - - |
| Voxels                | Instant  | Real-time    | High (O(n³)) | Low       |
| Point Clouds          | Instant  | Fast         | Moderate     | Moderate  |
| Meshes                | Varies   | Real-time    | Low          | Good      |
| NeRF                  | 24 hours | 30 sec/frame | Low (5MB)    | Excellent |
| Instant-NGP           | 5 min    | 50 ms/frame  | Moderate     | Very Good |
| 3D Gaussian Splatting | 30 min   | 5 ms/frame   | High (500MB) | Very Good |

The future of 3D vision combines these approaches: Gaussian Splatting for rendering, neural fields for compression, and hierarchical structures for scale — each solving different pieces of the 3D puzzle.

Highly Cited Paper Check: Attention Is All You Need

“Attention Is All You Need” is one of the most-cited papers of the twenty-first century, and although it is a language paper, its machinery shows up all over modern 3D vision — most directly through positional encoding, which NeRF lifted almost verbatim to make MLPs see high-frequency geometry. Here is where each 3D-relevant concept from the paper lives in this post:

Transformer paper concept                          Covered where in this post
----------------------------------------------------------------------------------
Sinusoidal positional encoding sin/cos(2^k · x)    Chapter 10: NeRF - Positional Encoding
Frequency encodings unlock high-frequency detail   Chapter 10: NeRF - Positional Encoding ("Why it works")
Learned embeddings for discrete entities           Chapter 8: Neural SDF - DeepSDF latent shape codes
Position-wise feed-forward (MLP) layers            Chapters 8 & 10: DeepSDF and NeRF MLP backbones
Residual / skip connections                        Chapter 8: Neural SDF skip connections
Normalized weights as a distribution over inputs   Chapter 10: Hierarchical sampling - volume rendering
                                                   weights used as an importance distribution
Encoder-decoder split (representation → output)    Chapter 19: DreamFusion - frozen CLIP text encoder
                                                   guiding a NeRF "decoder"

The remaining transformer concepts — scaled dot-product attention, multi-head attention, tokenization, softmax, cross-entropy, label smoothing, BLEU, Adam, learning-rate warmup, and dropout — are not 3D vision or neural rendering concepts, so they are intentionally not covered in this post.

Recent Paper Check: SkyJEPA (arXiv:2606.23444)

SkyJEPA is a recent JEPA-style latent world model for real-time quadrotor control: temporal convolutional encoders compress sensor history into a latent state, a recurrent predictor rolls that latent forward, a physics-inspired prober decodes frozen latents into interpretable physical state (position, velocity, attitude), and a sampling-based controller plans over the model’s predictions. Its 3D content is exactly the rigid-body geometry this post now covers:

SkyJEPA concept                                    Covered where in this post
----------------------------------------------------------------------------------
Quadrotor pose: 3D position + attitude             Chapter 5: Homogeneous Coordinates and
                                                   Rigid-Body Transformations (SE(3))
World frame vs body frame conventions              Chapter 5: Coordinate Frames
Attitude as a rotation (quaternion / SO(3))        Chapter 6: Rotation Representations - Quaternions
Prober output: interpretable pose state            Chapters 5-6: the pose and rotation
(position, velocity, orientation)                  representations the prober decodes into
Composing pose changes along a trajectory          Chapter 5: chaining rigid transforms
Predicting how 3D state evolves over time          Chapter 17: Dynamic 3D Scenes (deformation
                                                   fields and time-varying Gaussians)

Everything else in SkyJEPA — TCN encoders, GRU predictors, the JEPA objective, the SIGReg collapse regularizer, MPPI sampling-based control, reward design, and sim-to-real transfer — is not a 3D vision or neural rendering concept, so it is intentionally not covered in this post.