Principled Math for Graphics
In graphics, if things are upside down, inside out, backwards, and reversed, it is often
tempting to start flipping random signs until all the errors cancel out.
This is an experience I want to remove from my life, so with this page, I am hoping to
lay out all the math I've used for the various kinds of renderers I've devised over the years.
By building up every equation from first principles, any mistakes in an implementation can also be
derived using those same principles, and overall code comprehensibility improves.
Table of contents
- 1. Basic ray-plane raycasting
- 1.1 Start of intersection
- 1.2 Actual intersection
- 1.3 Further optimization
- 1.4 Closing remarks, implementation details
- 2. Basic ray-sphere raycasting
1. Basic ray-plane raycasting
There are multiple ways to put textures on the screen, and many of them take an approach of plane-to-screen projection. Instead of that, this method renders from the screen outwards.
The ray \(\vec{r}\) in world space of every pixel can be described as a single-variable parametric equation: \begin{align} \vec{r}(\lambda) = \vec{p} + \lambda\vec{l} \end{align} Where \(\vec{p}\) is the position of the camera, and \(\vec{l}\) is the direction of that pixel's ray. \(\lambda\) expresses the amount of \(\vec{l}\) traveled by the ray; it is proportional to distance from the camera.
A plane can be similarly described with a double-variable parametric equation: \begin{align} \vec{t}(\alpha, \beta) = \vec{q} + \alpha\vec{v} + \beta\vec{w} \end{align} \(\vec{v}\) and \(\vec{w}\) are the plane's basis vectors, and \(\vec{q}\) is the position of the plane's origin. \(\vec{t}\) is the position in world space of a given normalized texture coordinate \((\alpha, \beta)\).
1.1 Start of intersection
With these basic definitions out of the way, we can already define what a ray-plane intersection is:
\begin{align}
\vec{r}(\lambda) &= \vec{t}(\alpha, \beta)\\
\vec{p} + \lambda\vec{l} &= \vec{q} + \alpha\vec{v} + \beta\vec{w}
\end{align}
This equation can be rearranged such that all unknowns (\(\alpha, \beta, \lambda\)) are on one side.
\begin{align}
\vec{p}-\vec{q} = \alpha\vec{v} + \beta\vec{w} - \lambda\vec{l}
\end{align}
So far these definitions have been dimensionally agnostic, but from here on, we assume that world space is 3-dimensional.
The acquired equation can be shown to be a system with 3 unknowns and 3 equations, meaning that it may be fruitful to try and solve it.
\begin{align}
p_x - q_x &= v_x \alpha + w_x \beta - l_x \lambda\\
p_y - q_y &= v_y \alpha + w_y \beta - l_y \lambda\\
p_z - q_z &= v_z \alpha + w_z \beta - l_z \lambda\\
\end{align}
This, in turn, can be translated to a matrix-based equation.
\begin{align}
\begin{bmatrix}p_x - q_x\\ p_y - q_y\\ p_z - q_z\end{bmatrix} &=
\begin{bmatrix}v_x& w_x& -l_x\\v_y& w_y& -l_y\\ v_z& w_z& -l_z\end{bmatrix} \cdot
\begin{bmatrix}\alpha\\ \beta\\ \lambda\end{bmatrix}\\\\
\vec{n} &= M\cdot\vec{g}
\end{align}
Now we can say: "Oh, the solution is trivial! It's simply:
\begin{align}
\vec{g} = M^{-1}\cdot\vec{n}
\end{align}
Where \(\vec{g}\) holds the unknowns!", and call it a day.
However, this method is likely to perform many redundant calculations, and is
very dependent on the implementation of the 3x3 matrix inverse.
1.2 Actual intersection
As basis for solving this system we'll use Cramer's Rule,
where the determinants of various carefully constructed matrices are combined together to calculate the unknowns.
For our case, that looks like this (put on appropriate welding goggles for the amount of variables from here on):
\begin{align}
d &= \begin{vmatrix}v_x& w_x& -l_x\\v_y& w_y& -l_y\\ v_z& w_z& -l_z\end{vmatrix}\\\\
\alpha = \frac{\begin{vmatrix}n_x& w_x& -l_x\\n_y& w_y& -l_y\\ n_z& w_z& -l_z\end{vmatrix}}{d} &&
\beta = \frac{\begin{vmatrix}v_x& n_x& -l_x\\v_y& n_y& -l_y\\ v_z& n_z& -l_z\end{vmatrix}}{d}\\
\lambda &= \frac{\begin{vmatrix}v_x& w_x& n_x\\v_y& w_y& n_y\\ v_z& w_z& n_z\end{vmatrix}}{d}
\end{align}
A reader familiar with determinants may already be able to tell that there is a lot of oppertunity
for factoring out repeated calculations here, so let's explore further what an implementation of this would look like.
For ease of verifying that the determinant is correctly implemented, we'll define a new vector \(\vec{l}'\): \begin{align} \vec{l}' = -\vec{l} \end{align} With this, we can start cooking the spaghetti: \begin{align} d &= v_x(w_yl'_z - w_zl'_y) - w_x(v_yl'_z - v_zl'_y) + l'_x(v_yw_z - v_zw_y)\\ \alpha &= \frac{n_z(w_yl'_z - w_zl'_y) - w_x(n_yl'_z - n_zl'_y) + l'_x(n_yw_z - n_zw_y)}{d}\\ \beta &= \frac{v_x(n_yl'_z - n_zl'_y) - n_x(v_yl'_z - v_zl'_y) + l'_x(v_yn_z - v_zn_y)}{d}\\ \lambda &= \frac{v_x(w_yn_z - w_zn_y) - w_x(v_yn_z - v_zn_y) + n_x(v_yw_z - v_zw_y)}{d} \end{align} We can now also start pulling out common factors. \begin{align} C_{wl'} &= w_yl'_z - w_zl'_y\\ C_{vl'} &= v_yl'_z - v_zl'_y\\ C_{nl'} &= n_yl'_z - n_zl'_y\\ C_{vw} &= v_yw_z - v_zw_y\\ C_{vn} &= v_yn_z - v_zn_y\\ C_{wn} &= w_yn_z - w_zn_y\\ \end{align} \begin{align} d &= v_xC_{wl'} - w_xC_{vl'} + l'_xC_{vw}\\ \alpha &= \frac{n_zC_{wl'} - w_xC_{nl'} - l'_xC_{wn}}{d}\\ \beta &= \frac{v_xC_{nl'} - n_xC_{vl'} + l'_xC_{vn}}{d}\\ \lambda &= \frac{v_xC_{wn} - w_xC_{vn} + n_xC_{vw}}{d} \end{align} Note the flipped sign on the last term of the calculation for \(\alpha\). Finally, we'll undo \(\vec{l}'\), back to \(\vec{l}\). \begin{align} C_{wl} &= w_zl_y - w_yl_z\\ C_{vl} &= v_zl_y - v_yl_z\\ C_{nl} &= n_zl_y - n_yl_z\\ C_{vw} &= v_yw_z - v_zw_y\\ C_{vn} &= v_yn_z - v_zn_y\\ C_{wn} &= w_yn_z - w_zn_y\\ \end{align} \begin{align} d &= v_xC_{wl} - w_xC_{vl} - l_xC_{vw}\\ \alpha &= \frac{n_zC_{wl} - w_xC_{nl} + l_xC_{wn}}{d}\\ \beta &= \frac{v_xC_{nl} - n_xC_{vl} - l_xC_{vn}}{d}\\ \lambda &= \frac{v_xC_{wn} - w_xC_{vn} + n_xC_{vw}}{d} \end{align}
1.3 Further optimization
From here on, we might even consider what values change more often than others. We know, for example,
that \(\vec{l}\) changes for every pixel, so we can segregate calculations based on this.
Cold path:
\begin{align}
C_{vw} &= v_yw_z - v_zw_y\\
C_{vn} &= v_yn_z - v_zn_y\\
C_{wn} &= w_yn_z - w_zn_y\\
\end{align}
\begin{align}
\lambda' &= v_xC_{wn} - w_xC_{vn} + n_xC_{vw}
\end{align}
Hot path:
\begin{align}
C_{wl} &= w_zl_y - w_yl_z\\
C_{vl} &= v_zl_y - v_yl_z\\
C_{nl} &= n_zl_y - n_yl_z\\
\end{align}
\begin{align}
d &= v_xC_{wl} - w_xC_{vl} - l_xC_{vw}\\
\alpha &= \frac{n_zC_{wl} - w_xC_{nl} + l_xC_{wn}}{d}\\
\beta &= \frac{v_xC_{nl} - n_xC_{vl} - l_xC_{vn}}{d}\\
\lambda &= \frac{\lambda'}{d}
\end{align}
Of course, we will want to render more than one plane in a scene, so the cold path will need to memoize
the pre-calculated values per plane. Lots of clever caching, SIMD, and threading can be thrown at this approach.
1.4 Closing remarks, implementation details
Many types of camera topology can be simulated with this approach. The standard way (that is elaborated upon most here) is rays eminating from a point, but one could have rays eminating from a plane, in parallel, or perhaps rays pointed inwards, eminating from a half-sphere, simulating the behaviour of a lens.
Similarly, many different types of texture sampling can be calculated using \(\alpha\) and \(\beta\). They can be clipped between 0 and 1, or infinitely repeat with a modulo, or even make a fading circle, mixing with the colors from other plane intersections.