8. Orthonormal Bases

Ray Tracing: The Rest of Your Life (v3.2.3): 8 Orthonormal Bases / 3.8 Orthonormal Bases

The previous chapter generated random directions around the $z$ axis. Actual rendering requires scattering directions to be generated around an arbitrary normal vector $\mathbf{n}$. The tool for doing so is an orthonormal basis (ONB).

This chapter covers three topics:

• Defining and constructing an ONB: deriving mutually orthogonal unit vectors $\mathbf{u},\mathbf{v},\mathbf{w}$ from a normal $\mathbf{n}$
• Transforming local coordinates into world coordinates: $x\mathbf{u}+y\mathbf{v}+z\mathbf{w}$
• Verifying cosine-weighted sampling through an ONB transform with ${\int }_{\text{hemi}}{\cos }^3\theta d\omega \approx \pi /2$

Relative Coordinate Systems

In a coordinate system with origin ${\mathbf{O}}^{\prime }$ and orthonormal basis vectors $\mathbf{u},\mathbf{v},\mathbf{w}$, a position is expressed as

$${\mathbf{O}}^{\prime }+u\mathbf{u}+v\mathbf{v}+w\mathbf{w}$$

The three ONB vectors are mutually orthogonal unit vectors:

$$\mathbf{u}\cdot \mathbf{v}=\mathbf{u}\cdot \mathbf{w}=\mathbf{v}\cdot \mathbf{w}=0,|\mathbf{u}|=|\mathbf{v}|=|\mathbf{w}|=1$$

Constructing an ONB from a Normal

After fixing $\mathbf{w}=\text{normalize}(\mathbf{n})$, the tangent-plane vectors $\mathbf{v}$ and $\mathbf{u}$ follow from these cross products:

$$\mathbf{v}=\text{normalize}(\mathbf{w}\times \mathbf{a}),\mathbf{u}=\mathbf{w}\times \mathbf{v}$$

Any vector not parallel to $\mathbf{w}$ can serve as $\mathbf{a}$. For numerical stability, choose

$$\mathbf{a}=\{\begin{array}{ll}(0,1,0)& |w_x|>0.9\\ (1,0,0)& \text{otherwise}\end{array}$$

When $|w_x|>0.9$, the vector $(1,0,0)$ is nearly parallel to $\mathbf{w}$, making their cross product close to zero, so the construction switches to the $y$ axis.

Transforming Local Coordinates into World Coordinates

The random_cosine_direction() function from the previous chapter returns local coordinates $(x,y,z)$ around the $z$ axis. To transform them into world coordinates, use

$$\text{direction}=x\mathbf{u}+y\mathbf{v}+z\mathbf{w}$$

local(a, b, c) evaluates this linear combination, while local_vec(d) is a thin wrapper that accepts an array $\mathbf{d}=(d_0,d_1,d_2)$.

Rust Implementation

The r308-onb crate implements the Onb structure and random_cosine_direction function.

r308-onb/src/lib.rs

/// Orthonormal basis built from a single surface normal.
///
/// `w` is aligned to the normal; `u` and `v` span the tangent plane.
struct Onb {
    u: [f64; 3],
    v: [f64; 3],
    w: [f64; 3],
}

impl Onb {
    /// Build an ONB whose `w` axis is aligned to `n`.
    fn build_from_w(n: [f64; 3]) -> Self {
        let w = normalize(n);
        // Pick an arbitrary vector not parallel to w.
        let a = if w[0].abs() > 0.9 {
            [0.0, 1.0, 0.0]
        } else {
            [1.0, 0.0, 0.0]
        };
        let v = normalize(cross(w, a));
        let u = cross(w, v);
        Self { u, v, w }
    }

    /// Transform local coordinates `(a, b, c)` to world space.
    ///
    /// Returns `a·u + b·v + c·w`.
    fn local(&self, a: f64, b: f64, c: f64) -> [f64; 3] {
        [
            a * self.u[0] + b * self.v[0] + c * self.w[0],
            a * self.u[1] + b * self.v[1] + c * self.w[1],
            a * self.u[2] + b * self.v[2] + c * self.w[2],
        ]
    }

    /// Transform a local direction vector `d` to world space.
    fn local_vec(&self, d: [f64; 3]) -> [f64; 3] {
        self.local(d[0], d[1], d[2])
    }
}

Lambertian scattering uses the ONB as follows. Because cos_theta comes from the dot product with $\mathbf{w}$, no explicit angle calculation is needed.

r308-onb/src/lib.rs

let local_dir = random_cosine_direction(&mut rng);
let world_dir = onb.local_vec(local_dir);
// cos(θ) relative to normal w = dot(world_dir, w)
let cos_theta = dot(world_dir, onb.w);
let pdf = cos_theta / PI;

generate_report produces two sections. It first constructs an ONB around the normal normalize([1, 2, 3]), then displays $\mathbf{u},\mathbf{v},\mathbf{w}$ and checks their orthogonality and unit lengths. It next estimates ${\int }_{\text{hemi}}{\cos }^3\theta d\omega \approx \pi /2$ with $N={10}^6$ samples. Unlike Section 3.7, it transforms the result of random_cosine_direction into world coordinates with onb.local_vec before computing $\cos \theta ={\mathbf{d}}_{\text{world}}\cdot \mathbf{w}$. An ONB transform preserves dot products, so $\cos \theta$ remains correct after the transformation.

Differences Between C++ and Rust

The C++ onb class stores an axis[3] array of vec3 values and provides indexed access through operator[]. The Rust version uses a structure with explicit u, v, and w fields.

In C++, build_from_w is a void method that mutates axis. In Rust, build_from_w follows the constructor pattern and returns Self, since returning a value is more idiomatic than mutating through a reference.

r308-onb/src/lib.rs

use std::f64::consts::PI;

struct XorShift64 {
    state: u64,
}

impl XorShift64 {
    fn new(seed: u64) -> Self {
        let state = if seed == 0 { 0x9e3779b97f4a7c15 } else { seed };
        Self { state }
    }

    fn next_u64(&mut self) -> u64 {
        let mut x = self.state;
        x ^= x << 13;
        x ^= x >> 7;
        x ^= x << 17;
        self.state = x;
        x
    }

    fn next_f64(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 * (1.0 / (1u64 << 53) as f64)
    }
}

fn dot(a: [f64; 3], b: [f64; 3]) -> f64 {
    a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
}

fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
    [
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0],
    ]
}

fn normalize(v: [f64; 3]) -> [f64; 3] {
    let len = dot(v, v).sqrt();
    [v[0] / len, v[1] / len, v[2] / len]
}

/// Orthonormal basis built from a single surface normal.
///
/// `w` is aligned to the normal; `u` and `v` span the tangent plane.
struct Onb {
    u: [f64; 3],
    v: [f64; 3],
    w: [f64; 3],
}

impl Onb {
    /// Build an ONB whose `w` axis is aligned to `n`.
    fn build_from_w(n: [f64; 3]) -> Self {
        let w = normalize(n);
        // Pick an arbitrary vector not parallel to w.
        let a = if w[0].abs() > 0.9 {
            [0.0, 1.0, 0.0]
        } else {
            [1.0, 0.0, 0.0]
        };
        let v = normalize(cross(w, a));
        let u = cross(w, v);
        Self { u, v, w }
    }

    /// Transform local coordinates `(a, b, c)` to world space.
    ///
    /// Returns `a·u + b·v + c·w`.
    fn local(&self, a: f64, b: f64, c: f64) -> [f64; 3] {
        [
            a * self.u[0] + b * self.v[0] + c * self.w[0],
            a * self.u[1] + b * self.v[1] + c * self.w[1],
            a * self.u[2] + b * self.v[2] + c * self.w[2],
        ]
    }

    /// Transform a local direction vector `d` to world space.
    fn local_vec(&self, d: [f64; 3]) -> [f64; 3] {
        self.local(d[0], d[1], d[2])
    }
}

fn random_cosine_direction(rng: &mut XorShift64) -> [f64; 3] {
    let r1 = rng.next_f64();
    let r2 = rng.next_f64();
    let z = (1.0 - r2).sqrt();
    let phi = 2.0 * PI * r1;
    let r = r2.sqrt();
    [phi.cos() * r, phi.sin() * r, z]
}

pub fn report_onb() -> String {
    generate_report()
}

pub fn generate_report() -> String {
    let mut out = String::new();

    // --- Section 1: ONB construction and orthogonality check ----------------
    let n = normalize([1.0, 2.0, 3.0]);
    let onb = Onb::build_from_w(n);

    out.push_str("=== ONB Construction: n = normalize(1, 2, 3) ===\n");
    out.push_str(&format!(
        "  w (normal):    [{:.6}, {:.6}, {:.6}]\n",
        onb.w[0], onb.w[1], onb.w[2]
    ));
    out.push_str(&format!(
        "  v (tangent):   [{:.6}, {:.6}, {:.6}]\n",
        onb.v[0], onb.v[1], onb.v[2]
    ));
    out.push_str(&format!(
        "  u (bitangent): [{:.6}, {:.6}, {:.6}]\n",
        onb.u[0], onb.u[1], onb.u[2]
    ));
    out.push_str("\n  Orthogonality check (should be ≈ 0):\n");
    out.push_str(&format!("  |u·v| = {:.2e}\n", dot(onb.u, onb.v).abs()));
    out.push_str(&format!("  |u·w| = {:.2e}\n", dot(onb.u, onb.w).abs()));
    out.push_str(&format!("  |v·w| = {:.2e}\n", dot(onb.v, onb.w).abs()));
    out.push_str("\n  Unit-length check (should be ≈ 1):\n");
    out.push_str(&format!("  |u| = {:.12}\n", dot(onb.u, onb.u).sqrt()));
    out.push_str(&format!("  |v| = {:.12}\n", dot(onb.v, onb.v).sqrt()));
    out.push_str(&format!("  |w| = {:.12}\n\n", dot(onb.w, onb.w).sqrt()));

    let n2 = 1_000_000usize;

    // --- Section 2: MC estimate via ONB-transformed cosine direction ---------
    out.push_str("=== Estimate ∫_hemi cos³θ dA with ONB-Transformed Cosine Directions ===\n");
    out.push_str("    (θ is the angle relative to normal w)\n");
    out.push_str(&format!("  π/2 (analytic solution) = {:.12}\n", PI / 2.0));

    let onb2 = Onb::build_from_w(n);
    let mut rng2 = XorShift64::new(0x1a2b_3c4d_0001);
    let sum2: f64 = (0..n2)
        .map(|_| {
            let local_dir = random_cosine_direction(&mut rng2);
            let world_dir = onb2.local_vec(local_dir);
            // cos(θ) relative to normal w = dot(world_dir, w)
            let cos_theta = dot(world_dir, onb2.w);
            // p = cos(θ)/π (cosine PDF), f = cos³θ → f/p = π·cos²θ
            cos_theta * cos_theta * cos_theta / (cos_theta / PI)
        })
        .sum();
    out.push_str(&format!(
        "  estimate (N={}) = {:.12}\n",
        n2,
        sum2 / n2 as f64
    ));

    out
}

#[cfg(test)]
mod tests {
    use super::*;

    /// The three ONB vectors are mutually orthogonal (dot products ≈ 0).
    #[test]
    fn onb_vectors_are_orthogonal() {
        let normals: &[[f64; 3]] = &[
            normalize([1.0, 0.0, 0.0]),
            normalize([0.0, 1.0, 0.0]),
            normalize([0.0, 0.0, 1.0]),
            normalize([1.0, 2.0, 3.0]),
            normalize([-1.0, 0.5, 0.2]),
        ];
        for &n in normals {
            let onb = Onb::build_from_w(n);
            assert!(dot(onb.u, onb.v).abs() < 1e-12);
            assert!(dot(onb.u, onb.w).abs() < 1e-12);
            assert!(dot(onb.v, onb.w).abs() < 1e-12);
        }
    }

    /// The three ONB vectors are unit length.
    #[test]
    fn onb_vectors_are_unit_length() {
        let onb = Onb::build_from_w(normalize([1.0, 2.0, 3.0]));
        assert!((dot(onb.u, onb.u).sqrt() - 1.0).abs() < 1e-12);
        assert!((dot(onb.v, onb.v).sqrt() - 1.0).abs() < 1e-12);
        assert!((dot(onb.w, onb.w).sqrt() - 1.0).abs() < 1e-12);
    }

    /// ONB-transformed cosine sampling estimates ∫_hemi cos³θ dA ≈ π/2.
    #[test]
    fn onb_cosine_sampling_estimates_pi_over_2() {
        let onb = Onb::build_from_w(normalize([1.0, 2.0, 3.0]));
        let n = 500_000;
        let mut rng = XorShift64::new(0xdead_beef_0001);
        let sum: f64 = (0..n)
            .map(|_| {
                let local_dir = random_cosine_direction(&mut rng);
                let world_dir = onb.local_vec(local_dir);
                let cos_theta = dot(world_dir, onb.w);
                cos_theta * cos_theta * cos_theta / (cos_theta / PI)
            })
            .sum();
        let est = sum / n as f64;
        let exact = PI / 2.0;
        assert!(
            (est - exact).abs() < 0.005,
            "ONB cosine estimate {est:.6} too far from π/2={exact:.6}"
        );
    }
}

Results in the Browser

The orthogonality check confirms that the absolute dot products are at the level of floating-point error, approximately ${10}^{-16}$. The estimate is close to the analytic solution $\pi /2\approx 1.5708$.

Summary

• I established a procedure for constructing an ONB from a normal $\mathbf{n}$: $\mathbf{w}=\text{normalize}(\mathbf{n})$, $\mathbf{v}=\text{normalize}(\mathbf{w}\times \mathbf{a})$, and $\mathbf{u}=\mathbf{w}\times \mathbf{v}$.
• local(a, b, c) transforms a local direction around the $z$ axis into world coordinates.
• I numerically confirmed that ONB-transformed cosine-direction sampling still produces ${\int }_{\text{hemi}}{\cos }^3\theta d\omega =\pi /2$.

The next chapter introduces PDF classes and a mechanism for sampling from mixtures of multiple PDFs.