Graphics Notes

Introduction

There are three main disciplines of graphics: modeling shape and appearance, animation, and rendering.

Modeling

Raster (bitmap) graphics have a fixed array of pixels with values, and a fixed resolution. Vector graphics are driven by display commands. It has a scalable resolution, which makes it ideal for typography and animation.

Rendering

The main question rendering aims to answer is: how do we make a 3D picture appear on a 2D plane?

There are two ways to think about the problem:

For every point on the paper, determine which object appears at the point
For every object in the scene, determine which points on the paper it covers

These two approaches correspond to two major classes of rendering algorithms: raytracing and rasterization.

Raytracing and rasterization

Raster Graphics

Images and Pixels

An image is a 2D grid of cells called pixels, where each pixel holds some set of values. A pixel (picture + element) is a unit with some kind of value. A pixel value can represent intensity (grayscale), color (RGB), opacity, or depth.

In the grayscale image below, each pixel requires only one value to encode its intensity. By convention, 0 is the lowest intensity (black) and 1 is the highest intensity (white).

Grayscale image

In a RGBA image, each pixel has four distinct values, which each represents red, green, blue, and alpha. Each value ranges from 0 to 255.

RGBA image

An RGB value of (255, 0, 0) gives us red, (0, 255, 0) gives us green, and (0, 0, 255) gives us blue. Surprisingly, an RGB value of (255, 255, 255) gives us pure white. This is because we are working with light from the computer screen rather than pigment from paint. This gives us additive color rather than subtractive color.

Why use RGB?

Why do we approximate the entire visible spectrum with three values?

The human eye has three types of color receptors (cones) that are each sensitive to a different range of visible wavelengths. These are sometimes called red, blue, and green cones (or short, medium, and long cones).

Notice that the medium and long cones are pretty close to each other.

Color receptors

Datatypes and Memory

Notice that intensity is represented with a float64, while color is represented with a uint8. The reason we use uint8 for RGB images is that float64 takes up a lot of storage. If we used floats for RGBA images, each pixel would require 16 bytes, and a 1000 x 1000 image would take up 16 MB of memory.

Most displays today use 24-bit true color, which is 8 bits per color channel. Since each color channel has $2^8 = 256$ distinct values, we generate $256^4 \approx 17$ million different RGB colors.

Is this enough color?

The human eye can discriminate only around 10 million colors, so this is generally sufficient. However, humans are fairly good at noticing uneven changes in color in color banding, which can be undesirable.

There are multiple ways to store images in memory:

Array of structs (AoS), where each struct represents an RGB value
Struct of arrays (SoA), where each array represents a R, G, or B channel Different memory organizations affect cache locality.

It is also common to store 2D arrays of data in 1D data structures. There are two common ways of doing this: (1) row-major ordering and (2) column-major ordering.

Row and column major ordering

References:

Vector Graphics

Pipeline

To create a vector graphic:

Define object geometry
Transform object (translate, rotate, scale)
Define "window" on the world
Render into raster graphics

In each of the steps of this pipeline, there are many different coordinate systems:

Object space is the local coordinate system in which a single object is defined
World space is the global coordinate system for the scene
Camera space is the coordinate system defined by the camera through which the scene is viewed

Why use multiple coordinate systems?

We want to keep different layers modular, so we can modify each independently. For example, we may want to scale or rotate individual objects. If all objects are defined in terms of world space, then there may not be a linear matrix transformation to do this.

Also, at higher levels of abstraction, we don't want to think about specific device coordinates to draw to. We may have to later consider the resolution of the device (normalized device coordinates) and the view plane of the OS window that the application runs in (screen space).

View Plane

A view volume contains everything the camera sees.

A conical volume approximates what the eyes can see
A frustrum is a rectangular pyramid which approximates the conical view volume
A parallel volume doesn't stimulate the eye or camera

The 3D is projected onto the view plane sitting in world space. It should be perpendicular and centered around the look vector. It should also be located somewhere between the eye and the scene.

The parameters of the view volume are:

Position
Look vector
Orientation (often represented by up vector)
Aspect ratio
View angle

A view angle¹ determines the amount of perspective distortion in the picture. It is equivalent to a photographer choosing a specific type of lens (e.g. wide-angle or telephoto lens). In a frustrum, there is both a height angle and a width angle.

Camera Space

The position, look vector, and up vector are enough to define the camera coordinate system.² The equivalent of $x$ , $y$ , and $z$ world coordinate axes in camera space are unit vectors $u$ , $v$ , and $w$ .

$w$ is in the opposite direction of the look vector
$v$ is the component of the up vector that is perpendicular to the look vector
$u$ is the vector perpendicular to both $w$ and $v$

w = \frac{-look}{\|look\|} \qquad v = \frac{up - (up \cdot w)w}{\|up - (up \cdot w) w\|} \qquad u = v \times w

Having a camera coordinate system allows us to change the orientation³ of the camera:

Roll: rotate camera around look
Yaw: rotate camera around $v$
Pitch: rotate camera around $u$

View Matrix

The canonical camera space is one where the camera sits at the origin, looks down at $-z$ , and is oriented such that $y$ is up. The view matrix transforms objects from world space into camera space.

View matrix

Suppose that the camera is in position $P$ in world space. The view matrix is made from a translation that moves the world by $-P$ , followed by a rotation that aligns $u$ , $v$ , $w$ with $x$ , $y$ , $z$ .

Note that a matrix $A$ with columns [ $u$ , $v$ , $w$ ] transforms the $x$ , $y$ , $z$ axes to $u$ , $v$ , $w$ directions. The inverse of $A$ does the opposite. As it turns out, because $A$ is a orthonormal matrix, the inverse of $A$ is just its transpose.

M_{\text{translate}} = \begin{bmatrix} 1 & 0 & 0 & -P_x \\ 0 & 1 & 0 & -P_y \\ 0 & 0 & 1 & -P_z \\ 0 & 0 & 0 & 1 \end{bmatrix}

M_{\text{rotate}} = \begin{bmatrix} u_x & u_y & u_z & 0 \\ v_x & v_y & v_z & 0 \\ w_x & w_y & w_z & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix}

By putting the rotation and translation together, we get one matrix for transforming other objects into the camera's coordinate system.⁴

M_{\text{view}} = M_{\text{rotate}}M_{\text{translate}}

Scenes

Graphics applications typically represent scenes with a DAG which stores:

Objects
Attributes (color, texture)
Transformations

To get the cumulative transformation matrix of an object, simply multiply the transforms along the path from the object node to the root. For example, to transform Object 1 to world space, first apply transformation A then C.

Note that the equivalent CTM would be $CA$ . Thus, it is fine to traverse down from the root node to the object nodes, keeping a running CTM by multiplying the current CTM with the node's CTM.

Scene graph example

To generalize a bit, each of the vertices building pointing up to the world root are model matrices which transform the object from object space to world space. The final matrix is a view matrix that transforms any object from world space to camera space. The inverses of these matrices provide the transformation of the opposite direction.

Model and view matrices

Raytracing

Raytracing takes the 3D objects we see at each point in camera space and projects them onto a view plane.

For each pixel, we need to solve two problems:

What object appears in the pixel?
What color should the pixel be?

To determine what object appears in the pixel, we trace a ray from the camera, through the middle of the pixel, and solve for the first point of intersection in the scene.

There are two types of rays:

Primary rays start from the camera and determine which parts of the scene are visible from the window.
Secondary rays starts at points in the scene, where the primary rays hit, and are traced to other parts of the scene. They are used to compute light effects.

Primary Rays

A primary ray should take the form $P + td$ , where $P$ is the origin of the ray and $d$ is the direction vector.

$P$ should be set to the camera's coordinates $(p_x, p_y, p_z)$ in world space.

We can calculate $d$ in three steps:

Given a coordinate of a pixel $(i, j)$ on the view plane, convert it into the corresponding 3D point $(u, v, k)$ on the view plane in camera space.
Transform $(u, v, k)$ into world space using the inverse view matrix, creating $(q_x, q_y, q_z)$ .
Then $d$ is simply $(p_x, p_y, p_z) - (q_x, q_y, q_z)$ .

The first step is given by

(i, j) \mapsto \left(W_v \left[\frac{i + 0.5}{W} - 0.5\right], H_v \left[\frac{j + 0.5}{H} - 0.5\right]\right)

Each coordinate is first shifted $0.5$ so that it is centered on the pixel. Then it is normalized by the dimensions of the scene plane so that it becomes a one-by-one grid. Then it is shifted by $0.5$ so that the center of the grid is at the origin $(0, 0)$ .

Finally, the $u$ and $v$ coordinates are scaled by the width or height of the view plane, $W_v$ or $H_v$ . The height and width of the view plane is given by the equations $2k\tan(\theta/2)$ , where $k$ is the distance from the camera to the view plane and $\theta$ is the height or width angle.

View plane height or width

Since the view plane sits at $z = -k$ in camera space, and $k$ is usually set to $1$ , this gives us a full 3D point $(u, v, -1)$ in camera space.

The inverse view matrix can be computed by taking the inverse of the view matrix, as described above in the View Matrix section.

Ray-Object Intersections

Ray-object intersections can easily be calculated with implicit objects.⁵ An implicit object is an object that can be definited by a function $f$ , where $f(Q) = 0$ iff $Q$ is a point on the service.

To find where a ray $P + td$ intersects an object, we just need to solve⁶ for the values of $t$ ⁷ such that $f(P + td) = 0$ .

Cubes

Consider a unit cube of side length $0$ sitting at the origin $(0, 0, 0)$ .

To determine whether a ray intersects with this cube, we should check whether it intercepts with any of its six faces:

If the ray is parallel to the face's plane, then there is no intersection
If there exists some point that intersects with one of the infinite planes, then there may be an intersection
If that point also lands within the bounds of the face, then it actually intercepts with the cube

Cube Intersection Calculations

An infinite plane with a normal vector $(A, B, C)$ that is distance $D$ away from $(0, 0, 0)$ is defined by all points $(x, y, z)$ satisfying:

Ax + By + Cz - D = 0

A unit cube can be described as six bounded planes whose normals are the unit vectors along each axis, pointing outward:

$(A, B, C) \in {(\pm1, 0, 0),\ (0, \pm1, 0),\ (0, 0, \pm1)}$

To check for ray-plane parallelism, we should check whether the normal is perpendicular to the direction vector, i.e. that $N \cdot d = 0$ .

To check for ray-plane intersection, we want to find some $t$ such that $p + dt = (x, y, z)$ satisfies the plane equation $Ax + By + Cz = 0.5$ .

By definition, this is equivalent to $N \cdot (p + dt) = 0.5$ .

The dot product is distributive, so we can derive: $t = (0.5 - N \cdot p)/(N \cdot d)$ .

Finally, we need to check that the point falls within $[-0.5, 0.5]$ on all axes.

Cylinders

Consider a cylinder of radius $0.5$ , spanning $y$ from $[-0.5, 0.5]$ . A cylinder can be conceptualized as a circle in the xz plane that is stretched along the y-axis.

To check if a ray intersects with a cylinder, we should check whether it intersects with the side of the cylinder or with either of the two end caps:

For the side, solve a quadratic for $t$ , then discard roots where $t < 0$ or where the hit's $y$ falls outside $[-0.5, 0.5]$
For each cap, intersect the ray with the plane $y = 0.5$ or $y = -0.5$ , then check that the hit point falls within the cylinder's radius

Cylinder Intersection Calculations

The side of the cylinder is the set of points satisfying:

x^2 + z^2 = 0.25, \quad y \in [-0.5, 0.5]

Substituting the ray $p + td$ for $(x, z)$ and collecting terms by power of $t$ gives a quadratic $at^2 + bt + c = 0$ , where:

$a = d_x^2 + d_z^2$
$b = 2(p_xd_x + p_zd_z)$
$c = p_x^2 + p_z^2 - 0.25$

The quadratic formula gives up to two roots for $t$ . For each root, we discard it if $t < 0$ or if the resulting $y$ falls outside $[-0.5, 0.5]$ , since that means the hit lands above or below the finite cylinder.

The surface normal at a hit point is the gradient of $x^2 + z^2 - 0.25$ , which is $(2x, 0, 2z)$ . This normalizes to the same direction as $(x, 0, z)$ , so we can use the hit point directly, with its $y$ coordinate zeroed out.

For each cap, we find where the ray intersects with the plane $y = y_\text{cap}$ , which is given by $p_y + d_yt = y_\text{cap}$ , or $t = (y_\text{cap} - p_y)/d_y$ .

We then check that the hit point at time $t$ lands within the disk $x^2 + z^2 \leq 0.25$ .

The normal is $(0, 1, 0)$ for the top cap at $y = 0.5$ , and $(0, -1, 0)$ for the bottom cap at $y = -0.5$ .

Cones

Consider a cone with an apex at $(0, 0.5, 0)$ and a base radius of $0.5$ at $y = -0.5$ .

To check if a ray intersects with a cone, we should check whether it intersects with the side of the cone or with the base disk:

For the side, solve a quadratic for $t$ , then discard roots where $t < 0$ or where the hit's $y$ falls outside $[-0.5, 0.5]$
For the base, intersect the ray with the plane $y = -0.5$ , then check that the hit point falls within the base's radius

Cone Intersection Calculations

The cone's radius shrinks linearly from $0.5$ to $0$ from the base at $y = -0.5$ to the apex at $y = 0.5$ . The radius at height $y$ is thus

r(y) = (0.5 - y)/2

The side of the cone is the set of points satisfying:

x^2 + z^2 = \left(\frac{0.5 - y}{2}\right)^2, \quad y \in [-0.5, 0.5]

Expanding the right side and moving everything to one side gives the implicit function:

f(x, y, z) = x^2 + z^2 - 0.25y^2 + 0.25y - 0.0625 = 0

Substituting the ray $p + td$ for $(x, y, z)$ and collecting terms by power of $t$ gives a quadratic $at^2 + bt + c = 0$ , where:

$a = d_x^2 + d_z^2 - 0.25d_y^2$
$b = 2p_xd_x + 2p_zd_z - 0.5p_yd_y + 0.25d_y$
$c = p_x^2 + p_z^2 - 0.25p_y^2 + 0.25p_y - 0.0625$

The surface normal at a hit point is the gradient of $f$ :

\nabla f = (2x,\ 0.5(0.5 - y),\ 2z)

which points in the same direction as the halved vector $(x,\ 0.25(0.5 - y),\ z)$ , so we can normalize that instead.

For the base, we find where the ray intersects with the plane $y = -0.5$ , which is given by $p_y + d_yt = -0.5$ , or $t = (-0.5 - p_y)/d_yt$ .

We then check that the hit point at time $t$ lands within the disk $x^2 + z^2 \leq 0.25$ .

The normal there is always $(0, -1, 0)$ , pointing straight down.

Spheres

Consider a sphere of radius $0.5$ at the origin. A sphere of radius $R$ at the origin has the equation:

x^2 + y^2 + z^2 - R^2 = 0

To check if a ray intersects with a sphere, we solve a quadratic for $t$ , then take the smallest non-negative root. Unlike the cube, cylinder, and cone, there are no bounds to check and no separate cap surfaces, since the sphere's implicit equation already describes a closed surface.

Sphere Intersection Calculations

Substituting the ray $p + td$ for $(x, y, z)$ gives:

(p + td) \cdot (p + td) - 0.25 = 0

Distributing the dot product and collecting terms by power of $t$ gives a quadratic $at^2 + bt + c = 0$ , where:

$a = d \cdot d$
$b = 2(p \cdot d)$
$c = p \cdot p - 0.25$

The quadratic formula gives up to two roots, $t_0 \leq t_1$ . We want the smallest non-negative root, since that is the closest intersection in front of the ray origin:

If $t_0 \geq 0$ , use $t_0$
Otherwise, if $t_1 \geq 0$ , use $t_1$
Otherwise, both roots are behind the ray origin, so there is no intersection

The surface normal at a hit point is the gradient of $x^2 + y^2 + z^2 - 0.25$ , which is $2(x, y, z)$ . This normalizes to the same direction as the hit point itself, so we can just normalize the hit point directly.

Lighting

Lighting is the simulation of how light interacts with objects. Because accurately simulating lighting is very hard, we usually used simplified lighting models.

Direct illumination is light striking a surface directly from a light source. Indirect illumination light bouncing off other objects toward the surface element. It is a much more difficult to simulate than direct illumination, but we can often approximate indirect illumination with direct illumination.

Light Sources

In graphics, there are three common models of light sources:

Type	Description	Parameters
Point light	A point that emits light equally in all directions	Position $p$ and light intensity $I$
Spot light	A point that emits light in a finite cone of directions	Position $p$ , light intesity $l$ , light direction $d$ , and cutoff angle $\theta$
Directional light	A light that shines in a particular direction; good for modeling the sun	Direction $d$ and light intensity $I$

Light Attenuation

Light attenuation is the idea that the amount of light illuminating an object should decrease with the square of distance. Thus, to be physically correct, when lighting a point that is distance $d$ away from a light source with intensity $I$ , we should attenuate the intensity as $I_\text{att} = \frac{I}{d^2}$ .

The problem is that this produces a very harsh lighting effect because there is no indirect illumination. Instead, we can use this hack:

I_\text{att} = f_\text{att} \cdot I = \min\left\{\frac{1}{c_1 + c_2d + c_3d^2}, 1\right\}I

Where $c_1$ , $c_2$ , and $c_3$ are experimentally defined constants.

Surface Normal Vectors

For illumination, we want to find the normal vector at the point $(x, y, z)$ . We can find this normal vector by taking the gradient at that point:

\hat{n} = \nabla(x, y, z) = \left( \frac{\partial f}{\partial x}(x, y, z), \frac{\partial f}{\partial y}(x, y, z), \frac{\partial f}{\partial z}(x, y, z) \right)

Once we have calculated the normal vector in object space, we need to transform the normal vector into world space.

Note that we cannot just apply the transformation matrix $M$ because $M\hat{n}_o \neq \hat{n}_w$ when there is scaling. Scaling distorts the normal in exactly the opposite sense of scale applied to a surface. For example, scaling $y$ by a factor of $2$ causes the normal to scale by $0.5$ .

The right equation is $\hat{n}_w = (M^t)^{-1}\hat{n}_o$ , unless the transformation includes a reflection (i.e. $\det(M) < 1$ ). The most foolproof equation is

\hat{n}_w = \text{sign}(\det(M))(M^t)^{-1}\hat{n}_o

BRDFs

Give a surface normal $\hat{n}$ , an incoming light direction $\hat{\omega}_i$ , and an outgoing light direction $\hat{\omega}_o$ , the bidirectional reflectance distribution function (BRDF) describes how much incoming light is scattered in that outgoing direction:

f_r(\hat{\omega}_i, \hat{\omega}_o) : \mathbb{R}^3 \times \mathbb{R}^3 \rightarrow [0, 1]

Some common types of BRDFs are:

Lambertian (diffuse): light scatters equally in all directions
Specular (glossy): light scatters more tightly around one direction
Mirror

The name of the Lambertian BRDF comes from Lambert's cosine law, which states that light intensity falls off with the cosine of the angle between $\hat{\omega}_i$ and $\hat{n}$ : $\cos(\theta) = \hat{n} \cdot \hat{\omega}_i$ .

Lambert's Cosine Law

Given $k_d$ as an RGB color, this gives us the Lambertian BRDF:

f_r(\hat{\omega}_i, \hat{\omega}_o) = k_d(\hat{n} \cdot \hat{\omega}_i)

Most materials are not perfectly Lambertian; most BRDFs are viewer dependent. This is where specular BRDFs. One of the most common Specular BRDFs is Phong BRDF:

f(\hat{\omega}_i, \hat{\omega}_o) = k_s(\cos \delta)^n = k_s(\hat{r} \cdot \hat{\omega}_o)^n

where:

$k_s$ is an RGB color
$\hat{\omega}_o$ is the vector to the viewer
$\hat{r} = 2(\hat{\omega}_i \cdot \hat{n})\hat{n} - \hat{\omega}_i$
$\delta$ is the angle between $\hat{\omega}_o$ and $\hat{r}$
$n$ controls how tight the specular highlight is

Phong Illumination Model

The Phong illumination model combines what we know about modeling light sources and materials to define a way to light objects.

First, it breaks lighting into three components:

Ambient stands in for the effect of indirect light; it affects all surfaces equally
Diffuse is the effect of direct light on a surface with a dull or rough finish
Specular is the effect of direct light on a shiny surface when the vector to the eye is closely aligned to the light's reflected rays

Phong Lighting Model

The Phong illumination model has three terms:

\text{Color} = \text{Ambient} + \sum_{l \in \text{lights}} [\text{Diffuse} + \text{Specular}]

Putting together the BRDFs for diffuse and specular reflection with light attenuation, the full equation, evaluated separately for each color channel $\lambda \in \{\text{red}, \text{green}, \text{blue}\}$ , is:

I_\lambda = \textcolor{#d97706}{k_a O_{a\lambda}} + \sum_{i=1}^{m} \textcolor{#e11d48}{f_\text{att}} \; \textcolor{#6b7280}{I_{\lambda,i}} \left[ \textcolor{#2563eb}{k_d O_{d\lambda} (\hat{N} \cdot \hat{L}_i)} + \textcolor{#059669}{k_s O_{s\lambda} (\hat{R}_i \cdot \hat{V})^n} \right]

where:

$\textcolor{#d97706}{k_a O_{a\lambda}}$ is the ambient term
$\textcolor{#e11d48}{f_\text{att}}$ is the attenuation factor⁸
$\textcolor{#6b7280}{I_{\lambda,i}}$ is the intensity of light
$\textcolor{#2563eb}{k_d O_{d\lambda} (\hat{N} \cdot \hat{L}_i)}$ is the diffuse term created by the Lambertian BRDF⁹
$\textcolor{#059669}{k_s O_{s\lambda} (\hat{R}_i \cdot \hat{V})^n}$ is the specular term created by the the Phong BRDF¹⁰

Here is a list of coefficients:

$k$ is a model parameter
$O$ is a material property
$\hat{L}_i$ is the unit direction from the point to the light
$\hat{R}_i = 2(\hat{N} \cdot \hat{L}_i)\hat{N} - \hat{L}_i$ is the reflection of $\hat{L}_i$ about $\hat{N}$
$\hat{V}$ is the unit direction from the point to the camera
$n$ is the shininess, a material property controlling how tight the specular highlight is

Note that the attenuation factor only multiplies the diffuse and specular terms, not the ambient term, since ambient light isn't coming from any particular light source's position. Directional lights have no position and therefore no attenuation (see Light Sources).

Note that the normal should point toward the viewer. You should check the dot product of normal with the view direction and flip the normal if it is negative.

Shading

We have to compute color for every pixel in our output image, but we don't necessarily have to evaluate the full Phong illumination equation at every equation.

Graphics systems often make the distinction between:

Lighting: what illumination model is evaluated and at what points
Shading: how every pixel is colored based on the lighting evaluations

This is mostly relevant for real-time graphics, with objects made out of polygons. With implicit objects, we usually just evaluate the exact equation at each pixel where an intersection occurs.

Common shaders include:

Flat shading: compute the normal (and lighting) once per face and apply that single color to all pixels in it
Gouraud shading: compute lighting at each vertex using the vertex normal, then linearly interpolate the resulting colors across the face
Phong shading: interpolate the normals themselves across the face, then run the lighting calculation per pixel using the interpolated normal

References:

Whitted Raytracing

There are some lighting effects that are not captured by Phong illumination, such as shadows, reflection, and refraction.¹¹ Turner Whitted introduced a recursive ray-tracing algorithm to capture these elements using shadow rays, rays of reflection, and rays of refraction.

More rays

Shadow Rays

To create shadows, we should shoot a "shadow ray" from the scene element¹² to the light source, and check if it is blocked by any other element in the scene.

If it is blocked, then do not allow the light to contribute a specular and diffuse term.¹³

Shadow rays

Reflection and Refraction

At each intersection with a surface, we perform raytracing again from that intersection point:

For a reflective surface, we trace the reflection ray about the normal vector at the surface intersection
For a refractive surface, we trace the refraction ray following Snell's law

We then add the color returned by that recursive call to the color of our intersection point, scaled by a material parameter $k_s$ or $k_r$ . We continue recursively spawning new rays until the contribution is negligible or maximum recursive depth is reached.¹⁴

Reflection Calculations

Given an incident ray direction $L$ pointing into the surface and a unit normal $N$ , decompose $L$ into a component along $N$ and a component tangent to the surface. Reflection keeps the tangential component and flips the normal component:

Reflection calculation

Refraction Calculations

Let $L$ be the incident ray direction, $N$ the surface normal oriented against $L$ , and $\eta = n_1/n_2$ the ratio of refractive indices. Then $\cos(\theta_1) = -L \cdot \hat{n}$ .

Snell's law states that

n_1 \sin(\theta_1) = n_2 \sin(\theta_2)

Dividing both sides by $n_2$ gives $\sin(\theta_2) = \eta \sin(\theta_1)$ . Squaring both sides and substituting $\sin^2(\theta_1) = 1 - \cos^2(\theta_1)$ gives:

\sin^2(\theta_2) = \eta^2(1 - \cos^2(\theta_1))

If $\sin^2(\theta_2) > 1$ , no real $\theta_2$ satisfies this equation.

Otherwise, $\cos(\theta_2) = \sqrt{1 - \sin^2(\theta_2)}$ .

To find the transmitted ray $t$ as a vector, decompose $L$ into a component along $N$ and a component tangent to the surface. Since $\cos(\theta_1) = -L \cdot N$ , the projection of $L$ onto $N$ is $(L \cdot N)N = -\cos(\theta_1)N$ , so:

L = \underbrace{L + \cos(\theta_1)N}_{\text{tangential}} - \cos(\theta_1)N

The sines in Snell's law are exactly the tangential components of these unit vectors, so the transmitted ray's tangential component is the incident ray's tangential component scaled by $\eta$ :

t_\text{tangent} = \eta(L + \cos(\theta_1)N)

The transmitted ray's normal component points further into the surface, the same way $L$ 's does, with magnitude $\cos(\theta_2)$ :

t_\text{normal} = -\cos(\theta_2)N

Adding these components together and simplifying gives the full transmitted ray:

t = t_\text{tangent} + t_\text{normal} = \eta L + (\eta \cos(\theta_1) - \cos(\theta_2))N

Extended Illumination Equation

Putting it together, this gives us the extended illumination equation:

\begin{aligned} I_\text{out} = {} & k_a + \sum_{l \in L} V(l, p)(f^{l}_\text{att} I^{l} k_d (\hat{n} \cdot \hat{\omega}^{l}_{i}) + f^{l}_\text{att} I^{l} k_s (\hat{r}^l \cdot \hat{\omega}_o)^n) \\ & + k_s I_\text{refl} + k_t I_\text{trans} \end{aligned}

where:

$V(l, p)$ is $0$ if the light ray is blocked and $1$ if the light ray is not blocked
$I_\text{refl}$ is the color returned by tracing a recursive ray reflected about the surface normal
$I_\text{trans}$ is the color returned by tracing a transmitted ray through the surface
$k_t$ is the material transmission efficiency

References: An improved illumination model for shaded display by Turner Whitted

Distributed Raytracing

Distributed raytracing¹⁵ extends Whitted raytracing by distributing ray samples across multiple domains, rather than tracing a single ray per sample point.

To achieve supersampling and antialiasing, multiple rays are shot through a single pixel into the scene.
To achieve motion blur, multiple rays are shot across time.
To achieve depth of field and defocus blur, multiple rays are shot through an area of the lens instead of just a point.
To achieve soft shadows,¹⁶ glossy reflections, and translucency, multiple rays are shot across an area.

Raytracing Pipeline

The full raytracing pipeline is as follows:

Ray tracing pipeline

Spatial Data Structures

The native ray tracing algorithm has a time complexity of O(MN), where M is the number of pixels and N is the number of objects. To reduce the time complexity, let's look at the O(N) part.

Bounding Volume Hierarchies

A bounding volume is a data structure that encloses complex objects in simpler ones. If the bounding volume isn't hit, then neither is the object inside it.

An axis aligned bounding box (AABB) is a cuboid where each edge is parallel to either the x, y, or z axis. This is convenient because we can uniquely represent the cuboid as its min and max values of x, y, and z.

A bounding volume hierarchy (BVH) arranges the bounding volumes in a hierarchical structure. There are two main approaches to building BVHs:

Find nearby objects, group them under one bounding volume. Repeat until the entire scene is under one bounding volume.¹⁷
Start with all objects in the scene, and split them into subsets using some kind of heuristic.¹⁸

BVH

Spatial Subdivisions

We can also break up space into uniform cells, and trace the ray through each cell. Compared to BVH, it takes more traversal time¹⁹ but less construction time.

We can also construct nonuniform, hierarchical cells called octrees. In an octant, each node is a AABB.

To construct an octree:

Find the bounding box of the entire scene
At each iteraction, partition the current node into 8 octants
Split the primitives at the current node into the octants. If a primitive spans both octants, put it in both.
Recurse on the octants.
Terminate on some maximum depth or maximum objects per octant

To traverse an octree:

Begin at the root node
If the node is a leaf node, test interaction with all of its primitives²⁰
Otherwise iterate through the node's children²¹ and compute intersection between ray and each child's bounding box. If it intersects, recurse on the child again.

Octrees

A kd-tree is a k-dimensional, axis-aligned binary tree. The benefit of a kd-tree over an octree is that we can be smarter about determining the size of nodes in the tree.²²

The main challenge is to decide on a split axis and split plane along the axis. Ideally, we'd like to balance the cost of entering a node (# of intersection tests) with the probability of entering the node.

This is an NP-complete problem, but a nice heuristic is the surface area heuristic:

For each split position,²³ calculate the cost cost = weight_l * cost_l + weight_r * cost_r²⁴
Use the split plane that minimizes this cost

The naive cost of constructing the kd-tree is $O(n^2 \log n)$ because for each of $O(n)$ split planes, we need to check $O(n)$ primitives to determine which side they are on. But if we sort the primitives along the split axis, we can achieve more of a $O(n \log^2 n)$ .

Summary

Structure	Build time	Traversal time	Good for dynamic scenes	Good for nonuniform density
BVH	Slow, $O(n \log n)$	$O(\log n)$	Good, it is a trie structure	Good, can adapt to clusters
Uniform grid	Fast, $O(n)$	$O(n)$ but typically faster	Great, cheap to rebuild since cells are fixed	Not good, empty or overcrowded cells waste work
Octree	Moderate, $O(n \log n)$	$O(\log n)$	Less good, objects can move between octants	Not good, can't adapt to clusters
Kd-tree	Slow, $O(n^2 \log n)$ naive, $O(n \log^2 n)$ with sorting	$O(\log n)$	Terrible, need to rebalance trees	Good, split planes can balance primitives

Texture Mapping

Other than shading, we may also want to add texture to an object. One way to do this is with texture mapping.

Pros of texture maps:

Can be stored once and reused
Can be easily compressed to reduce size
Can be rendered quickly
Fairly intuitive to work with for artists
Can be used for far-away objects (billboarding)

Cons of texture maps:

Surfaces still look smooth because geometry is unchanged

Texture mapping has two subproblems:

Surface parameterization: mapping the 3D surface to a 2D plane (UV plane)
Texturing: assigning a value to each point on the UV plane

There are two main types of texturing:

Image texturing: use $(u, v)$ to index into an image
Procedural texturing: use $(u, v)$ as input to some function which computes the texture

Surface Parameterization

The goal of surface parameterization is to find²⁵ a function $f(x, y, z) = (u, v)$ which describes the surface.

Certain 3D shapes are naturally defined by a parametric function, such as cubes, cylinders, cones, and spheres. In each case, we want a function that takes a point²⁶ $p = (p_x, p_y, p_z)$ on the surface and returns $(u, v) \in [0, 1]^2$ .

Cubes

Consider a unit cube spanning $[-0.5, 0.5]^3$ .

To find the UV coordinate of a point $p$ on a cube, given its normal $n$ :

Use the sign of each component of $n$ to determine which of the six faces $p$ lies on
Map the two coordinates of $p$ that aren't fixed by $n$ linearly into $[0, 1]$ , negating one where needed to keep the texture consistently oriented across faces

Drag the point around the shape below to see where it lands in the UV square. The checkerboard is the texture, drawn on the surface through the parameterization itself.

p(0.15, 0.10, 0.50)

u = px + 0.5 0.65

v = py + 0.5 0.60

(u, v)(0.65, 0.60)

Cube UV Mapping Calculations

Consider a point $p$ on the face whose normal points along $+z$ (i.e. $p_z = 0.5$ ). The other two coordinates, $p_x$ and $p_y$ , each range over $[-0.5, 0.5]$ and map directly and linearly into $[0, 1]$ :

u = p_x + 0.5, \qquad v = p_y + 0.5

Each coordinate is shifted by $0.5$ to map it from $[-0.5, 0.5]$ into $[0, 1]$ . Some coordinates are also negated before shifting to keep the texture's orientation consistent across faces.

Cylinders

Consider a unit cylinder of radius $0.5$ , spanning $y \in [-0.5, 0.5]$ .

To find the UV coordinate of a point $p$ on a cylinder, given its normal $n$ :

If $n_y \neq 0$ , $p$ is on a cap. Map $p$ 's $x$ and $z$ coordinates directly onto $[0, 1]$ .
Otherwise, $p$ is on the side. Find the angle $\theta$ around the circle and map it to $u$ , then normalize the height to get $v$ .

p(0.29, 0.15, 0.40)

θ = atan2(pz, px) 0.94

u = 1 - θ / 2π 0.85

v = py + 0.5 0.65

(u, v)(0.85, 0.65)

Cylinder UV Mapping Calculations

For the top cap ( $n_y > 0$ ):

(u, v) = (p_x + 0.5,\ -p_z + 0.5)

For the bottom cap ( $n_y < 0$ ):

(u, v) = (p_x + 0.5,\ p_z + 0.5)

For the side, atan produces a $\theta$ in the range $[-\pi/2, \pi/2]$ , which could map two perimeter positions to the same $\theta$ value. To fix that, use atan2, which produces a $\theta$ in the range $[-\pi, \pi]$ :

\theta = \text{atan2}(p_z, p_x)

To map $\theta$ to $u \in [0, 1]$ :

If $\theta < 0$ , $u = -\theta/(2\pi)$
If $\theta \geq 0$ , $u = 1 - \theta/(2\pi)$

These two cases agree at $\theta = \pm\pi$ , where both give $u = 0.5$ , and disagree at $\theta = 0$ , where $u$ wraps from $1$ back to $0$ . So the seam lands on the $+x$ meridian.

Since the cylinder spans $y \in [-0.5, 0.5]$ , the height maps directly into $[0, 1]$ by shifting:

v = p_y + 0.5

Cones

Consider a unit cone with apex at $(0, 0.5, 0)$ and base radius $0.5$ at $y = -0.5$ . Unlike the cylinder, the cone's side doesn't have a single normal direction, so the normal alone can't distinguish the base from the apex. Instead, we check the normal for the base, and the point's height for the apex.

To find the UV coordinate of a point $p$ on a cone, given its normal $n$ :

If $n_y = -1$ , $p$ is on the base. Map $p$ 's $x$ and $z$ coordinates into $[0, 1]$ , the same way as the cylinder's bottom cap.
Otherwise, if $p_y = 0.5$ , $p$ is at the apex. Every value of $\theta$ maps to the same point there, so $(u, v)$ defaults to a fixed value, $(0.5, 0.5)$ .
Otherwise, $p$ is on the side. Find $\theta$ and $u$ the same way as a cylinder, and let $v = p_y + 0.5$ .

p(0.10, 0.15, 0.14)

θ = atan2(pz, px) 0.94

u = 1 - θ / 2π 0.85

v = py + 0.5 0.65

(u, v)(0.85, 0.65)

Spheres

Consider a unit sphere of radius $0.5$ at the origin.

To find the UV coordinate of a point $p$ on a sphere:

Find the latitude angle $\phi$ of $p$ , then normalize it into $v$
If $v \in \{0, 1\}$ , $p$ is at a pole, where every value of $\theta$ maps to the same point, so $u$ defaults to $0.5$
Otherwise, find $\theta$ the same way as a cylinder to get $u$

p(0.25, 0.27, 0.34)

φ = asin(py / r) 0.57

v = φ / π + 0.5 0.68

θ = atan2(pz, px) 0.94

u = 1 - θ / 2π 0.85

(u, v)(0.85, 0.68)

Sphere UV Mapping Calculations

The latitude angle $\phi$ of a point at height $p_y$ on a sphere of radius $r$ is:

\phi = \sin^{-1}(p_y/r), \quad \phi \in [-\pi/2, \pi/2]

This maps to $v$ by:

v = \phi/\pi + 0.5

At $v = 0$ (south pole) and $v = 1$ (north pole), $\theta$ is undefined, so $u$ defaults to $0.5$ there instead of being computed.

Meshes

An arbitrary polygonal mesh has no closed-form parametric function, so UV coordinates must instead be assigned explicitly, as a $(u, v)$ pair stored at each vertex alongside its position and normal. There are several approaches to choosing these values:

Texture each face independently
Planar, spherical, or cylindrical mapping
Manual UV mapping

In manual UV mapping, the artist cuts seams into the object to unfold it into UV space. Seams are placed on portions of the model that are less visible or sections that could have discontinuities.

Manual UV mapping

Once each vertex of a triangle has a $(u, v)$ value, the $(u, v)$ for any point inside the triangle is found by interpolating the three vertex values, using the same barycentric-coordinate technique used to interpolate colors below. This is what lets a texture stretch smoothly across a triangle instead of being constant per-face.

Triangle UV mapping

If we can find the color on every vertex of the object, then we can find the color for any point on the surface.

First, let's consider interpolating between two values along a line. Given two colors $A$ and $B$ , we can compute any value along the line by evaluating

C(t) = tA + (1 - t)B

$t$ and $1 - t$ are the barycentric coordinates of the line segment. This technique can be generalized for triangles.

Image Texturing

In most cases, the image coordinates will be different from texture (UV) coordinates. Suppose that the image has a width $W$ and height $H$ , and the UV plane is the unit plane. Then the mapping is:

(u, v) \mapsto \left(\lfloor uW \rfloor, \lfloor (1 - v)H \rfloor\right)

For periodic textures (like brick), a better option for filling a piece of geometry is to tile the texture. Suppose that we have an image with width $W$ and height $H$ , and we want to repeat this image $m$ times horizontally and $n$ times vertically on the unit UV place. Then the mapping is:

(u, v) \mapsto \left(\lfloor umW \rfloor \bmod W,\ \lfloor (1 - v)nH \rfloor \bmod H\right)

Procedural Texturing

Procedural noise is used to create randomized data. A few examples of noise functions include white noise, value noise, perlin noise, simplex noise, and worley noise.

Types of noise

White noise:

Every pixel is given a random value independent of its neighbors
It is easy to implement and very fast
It may be "too noisy" or too discontinuous for certain use cases

Value noise:

Generate a grid of white noise values, then scale it up
It is fairly easy to implement, pretty fast, and locally continuous
The resulting noise is visibly aligned to the pixel grid

Perlin noise:

Intensity is a function of a grid of vectors instead of a grid of values
It is locally continuous and looks more organic
It is more challenging to implement and can be very inefficient if implemented poorly

An easing function has the following properties:

ease(0) is 0
ease(1) is 1
ease(x) is continuous

References: CSCI 1230 Lab 5: Terrain

Antialiasing

All images begin as a continuous signal. When a photograph is taken, this continuous signal is turned into a discrete signal in a process called sampling. Values stored in discrete image pixels are point samples of the original continuous signal.²⁷ Later on, the discrete signal can be reconstructed into a continuous signal.

When we undersample a high-frequency signal, leading to an incorrect reconstruction with a low-frequency signal, this is called aliasing. Aliasing can cause issues like jaggies²⁸ and crawlies.

We can avoid aliasing by taking some techniques from digital image processing.

Digital Image Processing

Digital image processing is the processing of images through algorithms in order to modify them or extract useful information from them. Their applications include noise removal, feature extraction, image compression, and image enhancement.

Digital image processing is really just a form of signal processing. A digital image can be viewed as a discrete 2D signal which maps a 2D coordinate to a pixel value. This gives us two ways to represent the same image:

In the spatial domain, the function $f(x, y)$ represents the value at each coordinate $(x, y)$ .
In the frequency domain, $F(u, v)$ represents the image as a sum of 2D sinusoids. Here $u$ and $v$ are spatial frequencies (cycles per unit distance) along $x$ and $y$ . Low frequencies correspond to coarse, slowly-varying features, while high frequencies correspond to fine detail and edges.

We get $F$ from $f$ with the Fourier transform.

The Nyquist Limit

The Nyquist limit says that we can avoid aliasing only if we sample at more than twice the highest frequency present in the signal. This gives us two solutions to avoid aliasing:

Raise the sample rate until it clears the signal's highest frequency. This is called supersampling.
Lower the signal's highest frequency until it fits under a sample rate we can afford. This is called prefiltering.

Note that a discontinuous signal has infinite frequency, so it is not always possible to satisfy the Nyquist limit by supersampling.

For geometry edges and shadow boundaries, we don't know where the discontinuity falls until we render, so we supersample.
For textures we have the entire signal on disk before rendering starts, so we can prefilter.

Supersampling

In supersampling, we sample at multiple places per pixel,²⁹ and then we combine these samples to get one final color for the pixel. This does not eliminate aliasing on edges, but it pushes the aliased frequencies high enough that the artifacts become small.

How should we distribute our samples for a pixel?

Regular grid sampling can still produce the stair-stepping artifacts (jaggies)
Random sampling can produce noisy artifacts
Stratified random sampling tries to combine the best of both worlds

Supersampling types

How many samples per pixel should we use?

Adaptive sampling schemes send out a small number of rays to check how much variation in color there is and check if more rays is warranted

How should we combine the samples to get the final color for the pixel?

Average the values of the samples

Supersampling a pixel is one case of a more general idea: distributed raytracing spreads samples across time, lens area, and light area as well as across the pixel.

Prefiltering

Prefiltering limits the signal, so that the frequencies we cannot represent are removed from the image. A low-pass filter keeps frequencies below some cutoff and discards everything above it.

In the frequency domain, a low-pass filter is just multiplication of the power spectrum by a unit box function.

Multiplication in the frequency domain is convolution in the spatial domain, so in practice we implement one as a convolution of the image with a blur kernel.³⁰

A box in the frequency domain is the ideal low-pass filter, but its dual in the spatial domain is $\text{sinc}$ , a continuous, oscillating function that is difficult to truncate.
The triangle filter is the linear interpolation function. The dual of the triangle filter is $\text{sinc}^2$ , since a triangle is a box convolved with itself.
The dual of the Gaussian filter is another Gaussian function. It is usually preferred over the box filter because the Gaussian decays exponentially, so it is easy to truncate.

Box and Gaussian filter

Texture Filtering

Prefiltering is especially useful for textures because we have the entire original signal.

When we shade a pixel, that pixel maps to some region of texture space. There are two cases that we could face:

Magnification: The pixel maps to less than one texel, and we are sampling the texture faster than it has detail. There is no aliasing, but there may be poor reconstruction.
Minification: The pixel maps to many texels, and we are sampling slower than the texture's detail. This causes undersampling and aliasing.

Magnification is a reconstruction problem, and minification is an aliasing problem. They need different fixes.

Bilinear Filtering

Bilinear filtering fixes magnification. The sample point almost never lands exactly on a texel center, so we find the four texels surrounding it and interpolate between them instead of snapping to the nearest one.

Let $(s, t)$ be the fractional position of the sample point within the square formed by those four texels, and let $c_{00}, c_{10}, c_{01}, c_{11}$ be their colors. Then bilinear filtering is two lerps along one axis followed by one lerp along the other:

c = \text{lerp}\big(\text{lerp}(c_{00}, c_{10}, s),\; \text{lerp}(c_{01}, c_{11}, s),\; t\big)

Up-scaling

Mipmapping

Texture pre-filtering fixes minification. Ahead of time we build a mipmap,³¹ a chain of successively halved, low-pass-filtered copies of the texture. At render time we measure how many texels the pixel covers, pick the level whose texels are about that size, and sample the texture.

We have two problems we need to solve:

How do we measure how many texels the pixel covers?
How do we determine the mip level from the texel coverage?

To measure the number of texels the pixel covers, we take a differential approach.

We want to calculate how a small shift in $x$ or $y$ on the view plane will cause a shift in the intersection point $p$ in world space. This is represented as $\frac{\partial p}{\partial x}$ and $\frac{\partial p}{\partial y}$ .

If we choose $dx$ and $dy$ to be the length and width of the pixel, then $\frac{\partial p}{\partial x} dx$ and $\frac{\partial p}{\partial y} dy$ will give us the sides of the quadrilateral covered by the pixel.

Then we want to calculate how a small shift in a point $p$ in world space will cause a shift in $u$ and $v$ on the UV map of the texture. Finally we want to calculate how a shift in $u$ and $v$ on the UV map will cause a shift in $s$ and $t$ on the actual texture map.

By chaining this together, we can calculate how a shift in $x$ or $y$ causes a change in $s$ and $t$ , giving us the amount of pixels that a pixel on the view plane covers.

Calculating region of texture map

Now, suppose we have $L$ mip levels, such that at $L = n$ , the scale factor is $1/2^n$ . If my desired scale factor is $1/s$ , then the mip level should be $L = \log_2 S$ .

Note that:

We have two scale factors $X$ and $Y$ . We can use the max of the two as the actual scale factor.
$X$ and $Y$ are not necessarily a power of $2$ , so we should also linearly interpolate between the nearest two mip levels

Because our texture lookups at each mip level are done via bilinear filtering, and then we combine them with another linear interpolation, the overall method is called trilinear filtering.

Real-Time Graphics

The real-time graphics pipeline is also called the rasterization pipeline or scanline rendering pipeline.

The main difference between raytracing and a real-time graphics pipeline is that we swap the loops.

# raytracing
for each pixel:
    for each primitive:
        trace ray 
        intersect object
        compute color

# real-time
for each primitive:
    for each pixel:
        check if primitive covers pixel
        compute color

The real-time pipeline can be parallelized better, since each object can be independently processed and each pixel can be independently shaded. However, we get less photorealistic lighting effects.

Trimeshes

Graphics platforms generally only know how to render triangles. Any shape can be decomposed into triangles; this process is called tessellation, and a shape tessellated into triangles is a trimesh.

The naive representation stores each triangle as three raw vertices, which duplicates shared vertices. The better representation is a list of vertices plus a list of triangles that index into it:

struct Triangle { int v[3]; }

struct TriMesh {
    std::vector<glm::vec3> verts;
    std::vector<Triangle> faces;
}

For each triangle:

Vertices should be specified in counter-clockwise (CCW) winding order
CCW order ensures that all triangles in a mesh have consistent normals
CCW order also enables backface culling, a rendering optimization that skips drawing triangle backs

Each vertex should store interleaved position and normal data. Normals can be:

Per-face: one normal per triangle, computed via the cross product of two edges, $\hat{n} = \frac{(v_2-v_1) \times (v_3-v_1)}{|(v_2-v_1) \times (v_3-v_1)|}$
Per-vertex: computed by interpolating the normals of the triangles touching that vertex, weighted by triangle area

Meshes are commonly stored on disk in the OBJ format. It lists vertex positions (v), vertex normals (vn), and faces that reference those by index (f 1//1 2//1 3//1).

Meshes can also be generated from:

Polygon modeling (Maya or Blender)
Curved surface modeling with splines or Bezier
Digital sculpting (ZBrush)
Procedural modeling (L-systems, shape grammars)
CAD modeling
Generative models

OpenGL

OpenGL is an industry standard API specification for real-time graphics. An implementation of OpenGL is usually implemented in C, with bindings in other languages.

There are other interfaces for programming the GPU:

Direct3D is Microsoft's graphics framework and runs only on Windows
Metal is Apple's graphics framework and optimized specifically for Apple's platforms
Vulkan is from the same group that maintains OpenGL (Khronos); it is a newer, lower-level API than OpenGL

OpenGL manages all objects and data by their ID. Every function either returns an object ID or takes in an object ID and mutates the object. Some middleware libraries (like three.js) will sit on top of the API and provide OOP support.

Furthermore, OpenGL is a state-based system. The behavior of functions depends on both the input and the state of OpenGL. OpenGL also uses the state to associate objects with one another through a process called binding. For example, to OpenGL to use a shader to render a particular mess, you need to declare the shader as the currently-bound shader to the object.

void initialize() { ... }

// Render the current scene (called every frame)
void draw() { ... }

// Free up any memory
void finalize() { ... }

// Event handlers
void onClick() { ... }

// Register all the above callbacks and start the event loop
void main() { ... }

To tell OpenGL to draw triangles, we send an array with all vertex information in a Vertex Buffer Object (VBO).

We also need to provide a Vertex Array Object (VAO) that describes the scheme of the VBO. For each attribute, the VAO takes in three parameters: size, stride, and pointer. Note that because VBO is a byte array, the parameters need to be multiplied by the size of the attribute.

OpenGL shaders are written in a DSL called GLSL. The shader programs are compiled on the CPU, and the compiled bytecode is sent to the GPU to be executed.

The vertex shader runs in parallel on every vertex contained in whatever VBO is bound when an OpenGL draw command is issued.

The inputs are the attributes contained in the VBO, as specified by the VAO, as well as the uniform variables declared and bound by the user
The outputs are the output vertex positions and and any other per-vertex outputs the user wants to compute

The fragment shader runs in parallel on every fragment produced by the rasterization stage

The inputs are the interpolated outputs from the vertex shader, as well as the uniform variables declared and bound by the user
The outputs are the depth of the fragment and any color outputs

Appendix

Vectors

A vector transforms a scalar into another scalar. Below, we transform the scalar $(x, y)$ by shifting it $a$ units right and $b$ units up.

(x, y) + \begin{bmatrix} a \\ b \end{bmatrix} = (x + a, y + b)

The dot product is a scalar representation of how similar two vectors are to each other.

\begin{bmatrix} a \\ b \end{bmatrix} \cdot \begin{bmatrix} c \\ d \end{bmatrix} = ac + bd

Geometric interpretations of the dot product:

$\begin{bmatrix} a \\ b \end{bmatrix}^2 = a^2 + b^2$ is the squared length of the vector
$\begin{bmatrix} a \\ b \end{bmatrix} \cdot \begin{bmatrix} c \\ d \end{bmatrix} = 0$ when the vectors are perpendicular
$v \cdot w = |v||w|\cos(\theta)$ , where $\theta$ is the angle between the arrows

The cross product takes in two vectors and produces a vector that is mutually perpendicular to both vectors. By convention, the direction that $a \times b$ points in is determined by the right-hand rule.

a \times b = \begin{bmatrix} a_2b_3 - a_3b_2 \\ a_3b_1 - a_1b_3 \\ a_1b_2 - a_2b_1 \end{bmatrix}

Geometric interpretations of the cross-product:

$a \times b = |a||b| \sin(\theta)$
$a \times b$ also gives the area of the parallelogram with $a$ and $b$ as edges. Half of this area gives us the triangle with $a$ and $b$ as edges.

Matrices

The identity matrix is a matrix such that $MI = IM = M$ for all matrices $M$ . The 2x2 identity matrix is:

\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}

If we have two matrices $M$ and $K$ such that $MK = I$ , then we say that $K$ is the inverse of $M$ . Not all matrices have an inverse. For a 2x2 matrix, we have:

\begin{bmatrix} a & b \\ c & d \end{bmatrix} \frac{1}{ad - bc} \begin{bmatrix} d & -c \\ -b & a \end{bmatrix} = I

Transformations

A transformation $T$ is a function that acts on space.

A linear transformation respects the property $T(a\hat{u} + b\hat{v}) = aT(\hat{u}) + bT(\hat{v})$
A matrix transformation is any transformation that can be written in terms of multiplying a matrix and a vector
Every linear transformation is a matrix transformation

Note that matrix transformations can compose. For example, applying $M$ then $K$ to vector $v$ is equivalent to $K(Mv) = (KM)v$ .

We can build matrix transformations using basis vectors. To find a matrix $M$ that sends:

\begin{equation*} e_1 = \begin{bmatrix} 1 \\ 0 \end{bmatrix} \text{to} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \quad\mathrm{and}\quad e_2 = \begin{bmatrix} 0 \\ 1 \end{bmatrix} \text{to} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} \end{equation*}

We need to find a matrix that satisfies:

\begin{equation*} \begin{bmatrix} a & b \\ c & d\end{bmatrix}\begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \quad\mathrm{and}\quad \begin{bmatrix} a & b \\ c & d\end{bmatrix}\begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} \end{equation*}

Arithmetically, this matrix is

\begin{bmatrix} x_1 & y_1 \\ x_2 & y_2\end{bmatrix}

To reverse the transformation, we can simply calculate the inverse of $M$ .

Thus, to send

\begin{equation*} \begin{bmatrix} a_1 \\ a_2 \end{bmatrix} \text{to} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \quad\mathrm{and}\quad \begin{bmatrix} b_1 \\ b_2 \end{bmatrix} \text{to} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} \end{equation*}

We can simply calculate the matrix $M$ that sends

\begin{equation*} \begin{bmatrix} 1 \\ 0 \end{bmatrix} \text{to} \begin{bmatrix} a_1 \\ a_2 \end{bmatrix} \quad\mathrm{and}\quad \begin{bmatrix} 0 \\ 1 \end{bmatrix} \text{to} \begin{bmatrix} b_1 \\ b_2 \end{bmatrix} \end{equation*}

And the matrix $K$ that sends

\begin{equation*} \begin{bmatrix} 1 \\ 0 \end{bmatrix} \text{to} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} \quad\mathrm{and}\quad \begin{bmatrix} 0 \\ 1 \end{bmatrix} \text{to} \begin{bmatrix} y_1 \\ y_2 \end{bmatrix} \end{equation*}

And then calculate $K(M^{-1})$ .

Scale, Rotate, and Translate

Both scale and rotate are linear transformations and can be represented by matrix transformations.

To scale by $s_x$ in the x-axis and scale by $s_y$ in the y-axis, apply the matrix:
$\begin{bmatrix} s_x & 0 \\ 0 & s_y \\ \end{bmatrix}$
To rotate by $\theta$ counterclockwise, apply the matrix:
$\begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \\ \end{bmatrix}$

Translation is not a linear transformation; it is an affine transformation. In order to represent it as a matrix transformation, we add another dimension to it. This is called homogeneous coordinates.

$w$ is $1$ if the value is a scalars and $0$ if the value is a vector
$w$ does not affect subtraction of two scalars (which produces a vector), addition of a vector and a scalar (which produces a scalar), or dot product (which is mostly between two vectors)
$w$ also does not affect the linearity of scaling and rotation $\begin{equation*} \begin{bmatrix} a & b \\ c & d\end{bmatrix} \quad\rightarrow\quad \begin{bmatrix} a & b & 0 \\ c & d & 0 \\ 0 & 0 & 1\end{bmatrix} \end{equation*}$

To translate $dx$ in the x-axis and $dy$ in the y-axis, apply the matrix:

\begin{bmatrix} 1 & 0 & dx \\ 0 & 1 & dy \\ 0 & 0 & 1\end{bmatrix}

Algorithmically, this works to out be:

\begin{bmatrix} 1 & 0 & dx \\ 0 & 1 & dy \\ 0 & 0 & 1\end{bmatrix} \begin{bmatrix} x \\ y \\ 1\end{bmatrix} = \begin{bmatrix} x + dx \\ y + dy \\ w\end{bmatrix}

Note that the order of the transformations matters. The most widely used ordering is TRS, which means to scale, rotate, and then translate.

References:

CSCI 1230 Lab 3: Transforms and Parsing

Non-Axis Aligned 3D Rotation

Scaling and translation in 3D is similar to that in 2D. Rotations are more complicated because we can now rotate around an arbitrary axis. There are several ways to do this: Euler angles, axis-angles, and quaternions.

Euler angles are a composition of rotations around the x, y, and z axes.

The rotation matrices for each axis, by angle $\theta$ , are:

R_x(\theta) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\theta & -\sin\theta \\ 0 & \sin\theta & \cos\theta \end{bmatrix}

R_y(\theta) = \begin{bmatrix} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \end{bmatrix}

R_z(\theta) = \begin{bmatrix} \cos\theta & -\sin\theta & 0 \\ \sin\theta & \cos\theta & 0 \\ 0 & 0 & 1 \end{bmatrix}

A full Euler angle rotation is the composition $R = R_z(\theta_z) R_y(\theta_y) R_x(\theta_x)$ , applied right to left.

There are several problems with Euler angles. First, they don't trace the shortest rotational path because they don't follow the geodesic. To fix that, we can linearly interpolate each angle independently with some correction factor. However, this does not produce an even angular velocity. Also, when two rotation axes become parallel, we lose a degree of freedom, leading to odd behavior at the point.

Axis-angle rotations involve decomposing a point $x$ into two parts relative to the axis $r$ :

$x_\parallel$ is the part parallel to $r$ ; rotation does not change it
$x_\perp$ is the part perpendicular to $r$ ; it rotates like a 2D point in a plain perpendicular to $r$

To calculate it, we need:

x' = x_\parallel + x\vdash\sin(\theta) + x_\perp\cos(\theta)

Which can be expanded to:

x' = (\hat{r}(\hat{r} \cdot x)) + (\hat{r} \times x)\sin(\theta) + (-\hat{r} \times (\hat{r} \times x))\cos(\theta)

The formula for the cross product can be represented as a matrix, so this entire equation can be rewritten in matrix form. Unfortunately, axis-angle rotations also do not sweep the shortest path between two points.

Quaternions are the only rotation method that can correctly do linear interpolation. A quaternion represents a rotation as a point on the surface of a 4D hypersphere. Because we are moving along the sphere, we get shortest path and constant angular velocity.

References: Visualizing quaternions

Convolution

Discrete 2D convolution takes in an image and kernel³², and produces an image.

(f * k)(i, j) = \sum_{m} \sum_{n} f(i - m, j - n) \, k(m, n)

To do convolution:

Flip the kernel³³³⁴
Allocate space for an output image
For each $(i, j)$ on the output image: a. Do element-wise multiplication between the rotated kernel and the region of the image centered on $(i, j)$ b. Then compute the sum of the elements of the resulting matrix c. Finally, set the value of the output image at $(i, j)$

For some intuition, consider what a few simple kernels do:

The identity kernel has a single 1 at the center and 0s elsewhere.
The right shift kernel has a 1 one position to the right of center.
The left shift kernel has a 1 one position to the left of center.

\begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix} \qquad \begin{bmatrix} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}

Near the border of the image the kernel hangs off the edge, so we have to decide what it reads there. The usual choices are to clamp to the nearest edge pixel, wrap around to the opposite side, or reflect across the edge.

Repeat convolution

Wrap convolution

Reflect convolution

Separable Kernels

Some 2D convolutions can be expressed as the outer product of two 1D kernels. A 2D convolution of an image with a separable kernel is equivalent to:

A 1D convolution of the columns of the image by the vertical 1D kernel
A 1D convolution of the rows of the image by the horizontal 1D kernel

A separable kernel can be faster than a non-separable kernel. With the 2D filter, we need to do roughly $F^2WH$ calculations. With the two 1D filters, we only need to do around $2FWH$ calculations.

Note that the view matrix does translation first and then rotation. This is because the matrix that transforms the camera to world space is composed of a rotation first and then a translation. The view matrix is the inverse of this matrix, so we must reverse the order of multiplication.
Most production raytracers model shapes using triangles, but this is fairly complicated.
Note that these objects are usually sitting in world space, but the implicit equations are defined in terms of object space. When solving for ray-object intersections, it is much easier to transform the ray into object space rather than transform the equation into world space.
There may be multiple values of $t$ that satisfies the equation, but we just care about the smallest, nonnegative value of $t$ .
See Light Attenuation
Lambertian BRDF measures the cosine similarity between the surface normal and the direction to the light. See BRDFs.
Phong BRDF meaures how closely the reflected ray $\hat{R}_i$ aligns with the unit direction $\hat{V}$ to the viewer. See BRDFs.
Light changing speeds as it enters a medium
Note that the shadow ray might intersect the starting surface at a t value close to 0. To avoid self-shadowing, we should add a small epsilon to the shadow ray's start position along the normal.
This will produce hard shadows; soft shadows are harder to compute.
Note that this only captures specular indirect light, not diffusion. A further extension of Whitted ray tracing is path tracing, which also bounces light around according to Lambert's cosine law.
"Distributed ray tracing" by Cook, Porter, and Carpenter
To implement soft shadows, we can model point lights as an area light instead. Then we can send multiple shadow rays across the area of the light, record the percentage that is blocked, and create a shadow gradient.
The bottom-up approach produces very tight bounds, but it is slow to construct because it is hard to "find nearby objects".
The top-down is easy to implement and efficient to run, but the trees it produces are generally not as good.
BVH is O(log n). Uniform cells is O(n) in the worst case, since all the objects could be lined up.
If an object is in two nodes, it may be helpful to cache the result of calculating the intersection on this object. This is known as mailboxing.
This can be optimized by picking the child that will be hit by the child first ("frontside child"). If there is an intersection in the frontside child that is earlier than the intersection to the backside child's bounding box, there is no need to traverse the backside child.
The cons of a kd-tree is that insertions and deletions are costly due to having to rebalance the tree...
We should only consider splits at the edges of primitives and edges of the bounding space.
If we assume that rays will be distributed evenly throughout the space, we can let weight_c be the ratio of the child's surface area to the parent's surface area. For simplicitly, we can also let cost_child be the number of primitives contained in the child.
Finding a good parameterization is an active area of study in graphics. Informally, a good parameterization is one that doesn't distort the texture too much when applied to the surface.
This point is usually in object space.
"A Pixel is Not A Little Square"
Here is an example of aliasing:
Most raytracing renderers will expose a "samples per pixel" (spp) parameter.
MIP stands for multum in parvo ("much in little") because it stores a lot of prefiltered texture information with little extra storage. In fact, it costs around $1/4 + 1/16 + \cdots = 1/3$ more storage.
A kernel is simply a 2D matrix of values. The type definition of image and kernel is the same.
The flipped kernel should be rotated 180 degrees, like this:

Instead of actually flipping our kernel, we can also cleverly index into the kernel "backwards" in a way that effectively flips it.
If we didn't flip the kernel, we would technically be implementing correlation.