rndr

What it is

rndr is a small 3D engine in Rust. It draws triangles on the CPU and has simple physics: ray casting, collisions and rigid bodies.

I used no math library and did almost no research. I wrote the vectors, the 3×3 matrices and the matrix inverse myself. SDL2 only shows the final pixels in a window. Meshes load from STL files.

Coordinates and math

The world has $x$ forward, $y$ to the left and $z$ up. Angles are in degrees.

V3 is a 3D vector with the dot product, the cross product, the length and the element-wise (Hadamard) product. M3x3 stores its three columns, so a matrix times a vector is a sum of the columns:

$$M\vec{v} = \vec{c}_0\,v_x + \vec{c}_1\,v_y + \vec{c}_2\,v_z$$

The inverse comes from the cofactors and the determinant. If $|\det M|$ is smaller than the float epsilon, the matrix has no inverse and the caller skips the case:

$$M^{-1} = \frac{1}{\det M}\,\operatorname{adj}(M)$$

Transforms

An object has a position $\vec{t}$ and a rotation $(\gamma, \beta, \alpha)$ around $x$, $y$ and $z$. The rotation matrix applies roll, then pitch, then yaw:

$$R = R_z(\alpha)\,R_y(\beta)\,R_x(\gamma) = \begin{pmatrix} c_\alpha c_\beta & c_\alpha s_\beta s_\gamma - s_\alpha c_\gamma & c_\alpha s_\beta c_\gamma + s_\alpha s_\gamma \\ s_\alpha c_\beta & s_\alpha s_\beta s_\gamma + c_\alpha c_\gamma & s_\alpha s_\beta c_\gamma - c_\alpha s_\gamma \\ -s_\beta & c_\beta s_\gamma & c_\beta c_\gamma \end{pmatrix}$$

where $s$ and $c$ are the sine and the cosine. A vertex of a mesh moves into the world with $\vec{p}_w = R\,\vec{p} + \vec{t}$.

Camera

The camera has the same transform. Its rotation turns the base vectors into its forward $\hat{f}$, right $\hat{r}$ and up $\hat{u}$ vectors. These three vectors are the rows of the camera matrix. Thus one matrix product gives the coordinates of a point relative to the camera:

$$\vec{p}_c = \begin{pmatrix} \hat{r}^{\mathsf T} \\ \hat{u}^{\mathsf T} \\ \hat{f}^{\mathsf T} \end{pmatrix} (\vec{p}_w - \vec{c}) = \begin{pmatrix} \hat{r}\cdot\vec{d} \\ \hat{u}\cdot\vec{d} \\ \hat{f}\cdot\vec{d} \end{pmatrix}, \qquad \vec{d} = \vec{p}_w - \vec{c}$$

The perspective projection is a pinhole. The screen is at an offset $\vec{s} = (0, 0, -100)$ from the pinhole, and similar triangles give the screen position:

$$x_s = \frac{s_z}{z_c}\,x_c + s_x, \qquad y_s = \frac{s_z}{z_c}\,y_c + s_y$$

$z_c$ is kept as the depth. If $z_c$ is smaller than 0.01, the division is skipped, because a very small $z_c$ makes lines millions of pixels long. If one vertex of a triangle is nearer than the near plane ($z_c \le 0.1$), the whole triangle is rejected. There is no clipping.

Rasterization

The triangles of a mesh are projected and rasterized in parallel with Rayon. Each triangle returns a list of pixel writes, and the main thread applies them.

For a triangle $abc$ on the screen, the rasterizer visits each pixel $p$ in the bounding box of the triangle, clamped to the screen. The area of a triangle comes from its corners:

$$A(p_1, p_2, p_3) = \tfrac{1}{2}\,\big|\,x_1(y_2 - y_3) + x_2(y_3 - y_1) + x_3(y_1 - y_2)\,\big|$$

$p$ splits $abc$ into three smaller triangles. If $p$ is inside, their areas add up to the full area:

$$\big|\,A(p,b,c) + A(p,a,c) + A(p,a,b) - A(a,b,c)\,\big| \le 0.01$$

The same areas give the barycentric weights of $p$:

$$\lambda_a = \frac{A(p,b,c)}{A(a,b,c)}, \quad \lambda_b = \frac{A(p,a,c)}{A(a,b,c)}, \quad \lambda_c = \frac{A(p,a,b)}{A(a,b,c)}$$

The weights interpolate the position, the normal, the color and the depth of the three vertices, for example $z = \lambda_a z_a + \lambda_b z_b + \lambda_c z_c$. The weights are calculated on the screen, so the interpolation is not perspective-correct. Large triangles near the camera show it.

The interpolated values go to the fragment shader. A shader is a trait with one method. It gets the position of the pixel in the world and in the mesh, and it changes the output depth and color. A new shader does not change the rasterizer.

The pixel buffer keeps a depth for each pixel. A write happens only if its depth is not larger than the stored depth.

Ray casting

A point on the triangle $abc$ is $a + v(b - a) + w(c - a)$. A point on the ray is $\vec{o} + t\,\vec{d}$. When they are equal, the result is a 3×3 linear system:

$$\begin{pmatrix} b - a & c - a & -\vec{d} \end{pmatrix} \begin{pmatrix} v \\ w \\ t \end{pmatrix} = \vec{o} - a$$

I invert the matrix and get $v$, $w$ and $t$ in one product. The ray hits the triangle if $t \ge 0$, $v \ge 0$, $w \ge 0$ and $v + w \le 1$. With $u = 1 - v - w$, the same interpolation as in the rasterizer gives the position and the normal of the hit.

For a sphere with center $\vec{m}$ and radius $r$, the condition $|\vec{o} + t\,\vec{d} - \vec{m}|^2 = r^2$ is a quadratic in $t$:

$$(\vec{d}\cdot\vec{d})\,t^2 + 2(\vec{n}\cdot\vec{d})\,t + (\vec{n}\cdot\vec{n} - r^2) = 0, \qquad \vec{n} = \vec{o} - \vec{m}$$

The normal at the hit is the direction from the center to the hit. A ray against the scene returns the hit with the smallest $t$.

Collision detection

Mesh against mesh uses the ray cast. The center $\vec{m}$ of a mesh is the average of its vertices. For each vertex $\vec{v}$, in parallel, a ray starts at $\vec{m}$ in the direction of $\vec{v}$, with a maximum distance of $|\vec{v} - \vec{m}|$. If the ray hits the other mesh before it gets to the vertex, part of the other mesh is inside this mesh. The hit is the contact point, and the inverse of the hit normal is the contact normal.

Two spheres touch if the distance between the centers is at most $r_1 + r_2$. The contact point is halfway between the centers.

Each pair of objects is checked once per tick, so the check is $O(n^2)$ in the number of objects.

Rigid bodies

A rigid body has a mass $m$, a diagonal inertia $\vec{I}$, a linear velocity $\vec{v}$ and an angular velocity $\vec{\omega}$. Each tick is one explicit Euler step:

$$\vec{x} \leftarrow \vec{x} + \vec{v}\,\Delta t, \qquad \vec{\theta} \leftarrow \vec{\theta} + \vec{\omega}\,\Delta t, \qquad v_z \leftarrow v_z + g\,\Delta t, \quad g = -5$$

The angular velocity is added directly to the rotation angles.

For a collision with normal $\hat{n}$, each body has an offset $\vec{r}_i$ from its center to the contact point. Its momentum at the contact point is

$$\vec{p}_i = m_i\,(\vec{v}_i + \vec{\omega}_i \times \vec{r}_i)$$

This splits into a normal part and a tangent part:

$$\vec{n}_i = (\vec{p}_i \cdot \hat{n})\,\hat{n}, \qquad \vec{t}_i = \vec{p}_i - \vec{n}_i, \qquad N = |\vec{n}_1| + |\vec{n}_2|$$

Friction uses the Coulomb model. If the tangent momentum is small, static friction stops it. Otherwise dynamic friction slows it:

$$\vec{f}_i = \begin{cases} -\vec{t}_i & |\vec{t}_i| \le \mu_s N \\ -\mu_d N\,\hat{t}_i & \text{otherwise} \end{cases}$$

The impulses exchange the normal momentum of the two bodies and apply the friction:

$$\vec{J}_1 = -\vec{n}_1 + \vec{n}_2 + \vec{f}_1 - \vec{f}_2, \qquad \vec{J}_2 = -\vec{n}_2 + \vec{n}_1 + \vec{f}_2 - \vec{f}_1$$

Each impulse changes the linear velocity, and its torque changes the angular velocity. The division by $\vec{I}_i$ is element-wise:

$$\Delta\vec{v}_i = \frac{\vec{J}_i}{m_i}, \qquad \Delta\vec{\omega}_i = \frac{\vec{r}_i \times \vec{J}_i}{\vec{I}_i}$$

A body can lock its movement or its rotation. A locked body ignores the matching change.

What was hard

The camera matrix. I tried to derive it without the internet and used brute force for a long time. The answer was simple: the rows of the matrix are the forward, right and up vectors of the camera.

That was a good lesson for my ego. I wanted to derive everything myself. Most knowledge is a compression that other people already did, and it is rational to use it.