The matrix rain effect is a really common one among creative coders/shader enthusiasts, and for very good reason: it’s a cool-looking effect and it’s not incredibly complicated to achieve.

Because I’m a fan of retro-futuristic visuals and because it felt like a rite of passage as a shader programmer to create this effect, I decided to code my own version with a unique twist to it – instead of being the traditional matrix rain where the letters move downwards on the screen in 2D, mine would appear to scroll in 3D.

Here’s a gif of the result (click on the gif to see the code):

The characters themselves were rendered by sampling a glyph-containing texture provided by shadertoy.

To achieve the 3D scrolling effect, two primary techniques were used:

1. Faking a pair of parallel 2D planes.
2. Using a modified version of parallax scrolling to simulate depth.

We’ll go over each of these techniques in turn. Like my other shader tutorials, this one assumes that you have some prior experience writing shaders.

Faking Parallel 2D Planes

The effect can, in a nutshell, be thought of as two parallel textured planes, each split up into columns of text that are scrolling at different speeds and which have different opacity values. The first step is to take the above texture and get it to render on these planes:

This is actually a pseudo-3D effect – there’s no raymarching, matrix multiplication, or any traditional 3D math being used. Instead, we just modify the x and y values of the sampled coordinate based on the vertical distance of the pixel from the origin. Here’s the responsible code snippet:

vec2 uv = fragCoord.xy / iResolution.xy - vec2(.5);

// in the above example, the vec2 "lookup" defined below is used to
// sample the glyph texture directly 
vec2 lookup = vec2(uv.x / uv.y, .5 / uv.y);

Dividing uv.x by uv.y causes the x value to to increase more slowly as uv.y, which is the vertical offset from the origin, increases. This is why the grid squares at the top and bottom edges of the screen are wider than those closer to the center.

The calculation of lookup.y uses the function .5 / uv.y, which looks like this when plotted on a graph:

The x-axis is uv.y. The y-axis is lookup.y.

As you can see, the values of lookup.y change very quickly near the origin and slow down as uv.y increases. Observe the absolute value of the derivative of the above function to get another look into how the uv-value changes with respect to vertical distance from center:

The value of lookup.y changes more slowly as vertical UV distance increases, thus causing the squares closer to the top and bottom edges to be more vertically stretched out (in addition to being wider, as we illustrated previously). These distortions of lookup.x and lookup.y account for the planar appearance of the rendered result.

Parallax Effect

If you’ve dabbled in gamedev or web design, you’re likely familiar with the concept of parallax scrolling. Although in practice the technique is almost exclusively applied to scrolling backgrounds, the concept can actually be used for other interesting pseudo-3D effects. The multi-tiered appearance of each of the flat parallel planes in my matrix rain effect is created by two simple tricks:

1. Each column within each plane scrolls at a different pseudo-random speed.
2. The opacity of each column depends on its scroll speed.

Trick #1 uses a very similar principle to standard parallax scrolling: things that scroll more slowly appear to be further away. It turns out that this sort of applies even if the things that are scrolling are tightly packed columns. However, it’s really trick #2 that sells the illusion. Things that are fainter in color appear further away (interestingly, this is a technique that is commonly used in traditional landscape painting – desaturated colors are frequently used to convey distance from the viewer).

By making the slower-scrolling columns also less opaque, these columns get pushed even further into the background.

Donezo

That’s basically all there is to it. In my explanation of the effect, I purposely refrained from explaining certain aspects of the shader (coloring the text green, generating a random speed per column, etc.) in order to avoid over-explaining things that aren’t particularly complex. Here’s a link to the code again though, in case you’d like to take a closer look at the implementation details. It’s a short read – only 23 lines.

I’ve noticed recently that a lot of the really interesting 2D effects on shadertoy.com involve repeating patterns and grids of some kind. Standard Cartesian grids, triangle grids, and hexagon grids seem to be the most common. Hexagons in particular have a very cool and futuristic look to them, resulting in hexagonal patterns appearing in many sci-fi wallpapers and designs.

Here are a few effects that I created recently using this technique (click the gifs to see the code):

We’ll be exploring two particularly elegant ways to render hexagons. The first technique is better suited for rendering a single hexagon. The second technique generates a tiled hexagonal grid. Both of these techniques are useful if you want to achieve a wide range of hex effects and both of them are used in the effects that I showed above. They’re far from being the only hex rendering techniques out there, but they’re clean, fast, and the first ones that I came across. Before diving in, I need to give credit where it’s due: I learned these mathematical techniques by studying this shader by the shadertoy user Shane (I encourage you to check out the rest of his work by the way – the guy is a shader wizard).

For both algorithms below, assume that our goal is to render a hexagon (or a grid of hexagons) where each hexagon is exactly 1 unit wide. We will be dealing with pointy-topped hexagons rather than flat-topped hexagons, and we will be dealing with hexagonal grids which are staggered by row rather than by column.

Before proceeding, it will also be useful to know the height of the hexagon, which we can easily calculate based on its width:

Hex height is 2 / sqrt(3) because it is 2 * the length of the hypotenuse calculated above.

Method 1: Hexagonal Distance Function

To render a single hexagon, we’ll be using a signed distance function (SDF for short). An SDF is a way to describe a shape with a mathematical function. Typically, an SDF will return the distance of a given point from the boundary of the shape. Positive values indicate that the point that is outside of the shape, negative values indicate that the point is inside of the shape, and a value of zero indicates that the point is on the boundary of the shape.

If we wanted, for instance, to draw a red hexagon against a black background, we would implement a distance function for the hexagon and then assign a red color to pixels that are inside of the shape and a black color to pixels that are on the outside of the shape. Here’s the distance function for a hexagon.

float calcHexDistance(vec2 p)
{
     const float hexHalfWidth= .5;
     // 1.7320508 is sqrt(3)
     const float s = vec2(1, 1.7320508);
     p = abs(p);
     return max(dot(p, s * .5), p.x) - hexHalfWidth;
}

Here’s a shader which uses this formula to render a single hexagon. It’s just calculating the distance using the above function and then coloring the pixel based on that distance.

So how does all of this stuff work?

First, let’s just consider the top-right quadrant of the hexagon. Consider that for a hexagon of width 1, any point more than .5 units away from the center along the x-axis will lie outside of the hexagon. That’s why the p.x appears inside of the max function in the return statement. A point also is outside of the hexagon if it is more than .5 units along the diagonal line depicted below (don’t worry about the points in the other quadrants for now – that’s what the abs function call is for):

To see how far a point is along this diagonal line, you can simply take the dot product of the point and a unit vector pointing in the same direction as the diagonal. You can find this diagonal by taking advantage of the fact that hexagons are composed of equilateral triangles joined at the center. We’ll do this with a generic hexagon whose aforementioned diagonal is of length 1, since we need it to be a unit vector for the dot product to give us the displacement:

The leftmost angle in the red circle is 60 degrees because the red vectors radiating from the center both bisect 60 degree angles. 30 + 30 = 60.

As pictured above, to find the relevant diagonal we form a 30-60-90 triangle from the perpendicular bisector of the diagonal edge and then calculate the length of the sides to create a vec2 (by using half the width of the hexagon, .5, and taking advantage of the relationship between the lengths of the sides of 30-60-90 triangles which always have the same ratios with respect to each other).

The constant “s”, multiplied by .5, represents this diagonal (1.7320508 is the square root of 3). The max function call thus gives us the maximum of the displacements of the current point along the two axes that we care about (the horizontal axis and the diagonal).

You’ll notice that we’re operating on the absolute value of p rather than p. This is because, without the “p = abs(p)” line, we would get the following shape:

We don’t want this.

To get it to be an actual hexagon, we need to apply the logic to all four quadrants. Taking the absolute value of p gets all points into quadrant 1 so that we can use the same calculation. That way we don’t have to worry about checking the horizontal displacement of the point in the negative direction, nor do we have to find the displacement along the other three diagonals.

Finally, we subtract .5 from the resulting value to complete the distance function. With this subtraction, any fragment that either extends further than .5 in any horizontal direction or extends further than .5 in the diagonal direction (in any of the four quadrants) will be considered to be outside of the ‘gon. Those that don’t extend .5 or more in those directions will be considered to be inside of the ‘gon.

The value returned from this function gives you the distance of the fragment from the hexagon boundary. You can use this to shade the hexagon (by coloring any pixels with a non-positive distance value) a certain color, as is done in the example. With a bit of tweaking,  you can also use it to draw concentric hexagons/isolines, as Shane does in his shader. I’ll leave the implementation of that as an exercise to the reader.

Method 2: Hexagonal Grid

Interestingly, this next method of rendering a grid that I’m about to discuss doesn’t rely at all on the distance field formula that I described above. This is because we can take advantage of an interesting property of hexagons: if you have a staggered grid of points, as shown below, then the set of all pixels for whom a given point is the closest point forms a hexagon!

If you simply select a color for each pixel based on which of the staggered points is the closest, you automatically get a colored hex grid for free:

Here’s an example shader which does exactly that. Let’s talk about how it works.

At a high level, this is what we’re doing:

1. First, we divide the space into two different Cartesian grids – one whose cell centers represent the unstaggered points described above and another whose cell centers represent the staggered points (recall that these points represent the hexagon centers). The result of this is two staggered Cartesian grids who, if you collectively look at each cell center (the center of each “grid square” in the Cartesian grid), gives us the dotted pattern shown above in the image with the red and black dots.
2. Next, we figure out the closest hexagon center to the fragment by comparing the distance of the nearest cell center in each of the two grids (the staggered grid and the unstaggered grid).
3. FInally, we return the unique ID (in this case, we’re just using the position) of the closest hexagon center in the .zw components of the returned vec4 and returning the distance from the closest hexagon center in the .xy components of the returned vec4. We can later use this unique ID to uniquely color each hexagon and we can use the distance to draw isolines and render smaller hexagons within each cell.

The magic happens in the calcHexInfo function, displayed below.

vec4 calcHexInfo(vec2 uv)
{
     // remember, s is vec2(1, sqrt(3))
     vec4 hexCenter = round(vec4(uv, uv - vec2(.5, 1.)) / s.xyxy);
     vec4 offset = vec4(uv - hexCenter.xy * s, uv - (hexCenter.zw + .5) * s);
     return dot(offset.xy, offset.xy) < dot(offset.zw, offset.zw) ? 
          vec4(offset.xy, hexCenter.xy) : vec4(offset.zw, hexCenter.zw);
}

Let’s go over it line-by-line.

The first line essentially splits the space into a grid of squares of size (1, sqrt(3)), by dividing the UV coordinate by s.xy, and using round to figure out which grid square is the closest. The resulting value, stored in hexCenter, is the index of the nearest grid square. But why does the grid consist of squares of size (1, sqrt(3))?

Think of the grid squares as containing the non-staggered rows of the hexagonal grid (in our hexagonal grid, every other row will be staggered). We want these squares to be sized such that the hexagons are completely packed together horizontally, but leave enough space between the hexagons to fit in the staggered row. The cell width of 1 is obvious – our hexagons (as we decided early on in the article) are of width 1.  We use a height of sqrt(3) because sqrt(3) is exactly 3/2 of the height of each hexagon that we’re rendering (recall that the height of each hexagon is 2 / sqrt(3)). This allows us to stagger the rows perfectly, since the height 3 / 2 * hexHeight gives us exactly enough space between each vertically aligned hexagon to squeeze in a staggered hexagon. See the diagram below to see why this is the case:

.1/2 * h needs to equal n in order to perfectly fit the staggered hex row. We know that n is 1 / sqrt(3) – see the first diagram in this article to see why (simply double the .5 / sqrt(3) side). We also know that h is 2 / sqrt(3). So 1/2 * h = n.

Now let’s talk about the code. Notice that we pass a vec4 into the round function. Think of this vec4 as two vec2’s concatenated together – we’re sticking them in a single vec4 in a single line just for convenience. In the xy components, we simply store the UV value of the point. The zw components store a different set of coordinates which are used for the staggered rows. You’ll see how exactly these are used in the next step, but for now it’ll suffice to realize that we’re subtracting (.5, 1) in order to make sure that the rounded result of the line contains, in the zw component, the 1 by sqrt(3) cell immediately to the lower-left of any point that lies within any of the staggered hexagons. The number (.5, 1) is somewhat arbitrary – it could’ve been any vec2 that would be guaranteed to bring all points in the staggered hexagon (labeled “B”) shown below into the region indicated by the green rectangle after subtracting it from the fragment’s position:

For any point in a non-staggered hexagon, you’ll see that the zw value is irrelevant and will later be discarded.

The division by s.xyxy and the rounding are basically just taking each point (the original UV and the UV value offset to the lower left) and giving you the closest cell center.

In the next line, we calculate the offset of the fragment from the nearest cell center and from the staggered cell coordinate immediately to the top-right of the closest cell center to “fragment minus (.5, 1)”. This latter coordinate will always give you the correct cell center for fragments in staggered hex cells – for those that aren’t in staggered hex cells, this coordinate will be meaningless and will end up being ignored because we compare the offsets and discard the further one. We get these offsets by simply subtracting the positions of the respective hex centers from the uv coordinates. The cell center stored in hexCenter.xy is easy to calculate – simply multiply it by the cell size s (the round function gave us the index of the hex cell, not the actual position). The staggered cell center is calculated by offsetting the cell by half of the cell size s:

We then check which cell is closer to the point (the staggered cell or the unstaggered cell) by comparing the squared distances. Finally, we return the offset and unique ID (basically just the rounded index) of the closer hexagon.

Using these Functions

Armed with these functions, you now have everything you need to create some cool and unique hex-based effects. For instance, the sci-fi hexagon effect that I showed at the beginning of this article used the hex grid function to randomly color the hexagons a different shade of blue. The rotating lines inside of each hexagon are just hexagonal outlines, rendered using the distance function algorithm, that are clipped by rotating masks whose rotations are offset by a random number and by the current time.

Happy shader-writing! The only limit is your imagination. And your coding ability.

One of the earliest features that I implemented for Gridpulse Legions was fog of war. My intention was for the feature to work similarly to how it works in the original Starcraft; that is to say, it should satisfy the following criteria:

• Areas in the line of sight of the player’s units are not covered by fog.
• Areas that are not in the line of sight of the player’s units, but which have been explored previously, are covered by semi-transparent fog (you can see the terrain and any static structures/elements that were in that are last time you were there, but you cannot see any new non-allied units that are there).
• Fog covering unexplored areas is completely opaque – you cannot see anything there except for a thick black shroud.

At a high level, the solution I came up with was to render the visibility information into a RenderTexture, which is then drawn onto the terrain using a projector through a shader that takes care of rendering and blending the shroud.

Here’s a video of the feature in action:

Let’s go over each of the components in turn.

Recording Visible Areas

The first piece of information that we’ll need is the visible areas – that is, the areas that are currently within the line of sight of the player’s units. My solution involved two steps:

1. Generate a mesh representing the visible range of each unit.
2. Render these visibility meshes to a RenderTexture.

The visibility mesh generation can be as simple as a static mesh with its parent set to the unit whose visibility it is representing (so that it automatically moves with that unit). If your line of sight algorithm is more complex and takes elevation and obstacles into account, you may needed to procedurally update your mesh on the fly using a 2D visibility algorithm. The specifics of these algorithms are beyond the scope of this article and will differ per game, but for simplicity’s sake, you can just assume that this mesh is simply a static flat mesh floating above the unit’s head. It makes no difference to the rest of the implementation how the mesh is generated.

Place these meshes on their own layer. We’ll call this layer FOV and we’ll use it in the next step. Make sure that the main camera isn’t rendering this layer. Here’s what it’ll look like in the scene:

This mesh represents the visible range for a single unit.

Camera Setup

Once you have your visibility meshes, you’ll need to set up two cameras to render these meshes, each to their own render texture. Why two cameras? Because we’ll need a separate texture to display unexplored areas (completely opaque black) from the one used to display discovered areas that are not currently within the line of sight (semi-transparent black). Ultimately, we want to render textures that look like this

, which we will then project onto the terrain using a Projector component. The pink areas in the image above correspond to the visibility meshes.

We’ll render these textures by creating dedicated cameras for fog of war and configuring them to only render the visibility mesh layer. Our shader will look at this image and will project a black shroud onto the areas where the pink is not present and will ignore the areas where the pink is present.

First, create two new RenderTextures. You can actually get away with using fairly low-res RenderTextures if you turn anti-aliasing on. In Gridpulse Legions, the RenderTexture for the dark shroud (the shroud for unexplored areas) is 64×64. The RenderTexture for the normal shroud is 128×128.

Set up the camera size/position so that it covers the entire map and configure it to only render the FOV layer. Also set the Target Texture to the corresponding RenderTexture.

Both cameras will be identical in setup except for two things:

1. The camera for the dark shroud will have its clear mode set to Don’t Clear. As you’ll see a little bit later on, the openings in the fog will correspond to the locations of the visibility meshes within the camera’s viewport. As the visibility meshes move around, the RenderTextures will change to reflect the new locations of the meshes. But for the dark shroud, we actually want areas that have been explored to remain uncovered. See this image for clarification – notice that the dark shroud shows a superset of the visible areas that the normal shroud shows
   • 
2. The target RenderTextures will be different. The camera for the dark shroud will render to the dark shroud RenderTexture and vice versa for the other camera.

Rendering the Fog

With the above setup, you should now have two RenderTextures accessible during runtime whose pixels directly represent the areas of your map that need to be covered in fog. The next step is to actually draw the fog. For this, we’ll need a few different things:

1. A shader for rendering the fog.
2. Two materials, one for each fog type. Each material will have a different opacity value that will be used when rendering the fog (100% for the dark shroud and some lower value for the normal shroud).
3. A script, which will control the blending of the fog.
4. Two projectors, one for each fog type. Each projector will project its fog onto the map through the material described in step #2.

Let’s take a look at that shader. I’ll show just the fragment shader below, since the rest of it is just boilerplate shader code for projectors (see the full code snippet here).

fixed4 frag (v2f i) : SV_Target
{
    float aPrev = tex2Dproj(_PrevTexture, i.uv).a;
    float aCurr = tex2Dproj(_CurrTexture, i.uv).a;

    _Color.a = lerp(aPrev, aCurr, _Blend);
    return _Color;
}

As you can see, it’s using the alpha values of the shroud textures to figure out where to render the fog. When the pixel corresponds to an area of the shroud texture where the visibility mesh was rendered, the alpha will be 1 and the _Color variable will become transparent, resulting in the fog not being visible at that pixel.

You’re probably wondering now why we have _PrevTexture and a _CurrTexture and are lerping between the two instead of just using a single shroud texture for the lookup. This is a special technique for making the actual revealing/hiding of the fog look very smooth. This, plus the anti-aliasing, lets you achieve the illusion of smooth fog despite your fog texture not having nearly enough pixels to represent every visible spot on the map. As you’ll see in the next code snippet, each projector actually creates two textures internally – one for the previous visibility state and one for the current state. The script will lerp between the two and, upon completing the lerp, will blit the texture for the current state into the texture for the previous state and will blit the texture from the corresponding fog camera’s target texture into the texture for the current state. After doing this swap, it will repeat the process.

using System.Collections;
using UnityEngine;

public class FogProjector : MonoBehaviour
{
    public Material projectorMaterial;
    public float blendSpeed;
    public int textureScale;

    public RenderTexture fogTexture;

    private RenderTexture prevTexture;
    private RenderTexture currTexture;
    private Projector projector;

    private float blendAmount;

    private void Awake ()
    {
        projector = GetComponent();
        projector.enabled = true;

        prevTexture = GenerateTexture();
        currTexture = GenerateTexture();

        // Projector materials aren't instanced, resulting in the material asset getting changed.
        // Instance it here to prevent us from having to check in or discard these changes manually.
        projector.material = new Material(projectorMaterial);

        projector.material.SetTexture("_PrevTexture", prevTexture);
        projector.material.SetTexture("_CurrTexture", currTexture);

        StartNewBlend();
    }

    RenderTexture GenerateTexture()
    {
        RenderTexture rt = new RenderTexture(
            fogTexture.width * textureScale,
            fogTexture.height * textureScale,
            0,
            fogTexture.format) { filterMode = FilterMode.Bilinear };
        rt.antiAliasing = fogTexture.antiAliasing;

        return rt;
    }

    public void StartNewBlend()
    {
        StopCoroutine(BlendFog());
        blendAmount = 0;
        // Swap the textures
        Graphics.Blit(currTexture, prevTexture);
        Graphics.Blit(fogTexture, currTexture);

        StartCoroutine(BlendFog());
    }

    IEnumerator BlendFog()
    {
        while (blendAmount < 1)
        {
            // increase the interpolation amount
            blendAmount += Time.deltaTime * blendSpeed;
            // Set the blend property so the shader knows how much to lerp
            // by when checking the alpha value
            projector.material.SetFloat("_Blend", blendAmount);
            yield return null;
        }
        // once finished blending, swap the textures and start a new blend
        StartNewBlend();
    }
}

Both projectors will have this component. The script just creates two inner textures and does a continuous blend between the two textures by increasing the interpolation amount in a coroutine and feeding this value to the shader, which uses this float to lerp between the alphas from the previous and current textures.

Here’s what the actual projector component will look like in the inspector (make sure you configure it to cover your map and to project onto the correct layers based on your game’s specific needs):

Here’s what the two shroud materials look like. Notice that they differ only in the color (notice the white bars under the color, which represent the alpha, are different lengths):

The texture fields are empty because those are created in script.

With these components in your project/scene, you now have everything you need to get the fog working your game. The fog cameras will render the visibility meshes to the shroud textures. The projectors will then use these textures to render the shroud onto the terrain, blending the fog for smooth animation.

And there you have it. Attractive fog of war in Unity.

If you’re a game developer then you’ve probably, at some point, had to make an object spin. In your quest to rotate said object, you’ve undoubtedly come across the following matrix:

In this article, I’ll be diving into how exactly this matrix functions. I’ll be discussing the 2D rotation matrix but, as you’ll see, the 3D rotation matrix works in the exact same way and doesn’t require any new mathematical knowledge to be able to create or understand it.

To understand what the matrix is doing, we first need to understand what an actual point in 2D space represents. Consider the point (5, 3). As we know, this point represents a displacement of 5 units in the x-direction and 3 units in the y-direction.

But what exactly are the x and y directions?

Whenever you specify a point in space, it is always in the context of a coordinate system whose axes are represented by basis vectors. A 2D coordinate system is represented by two basis vectors (one for x and one for y), a 3D coordinate system is represented by three basis vectors (one for x, one for y, and one for z), etc. Typically when dealing with coordinates, we don’t think about the basis at all. This is because we usually default to using the standard basis – the basis where the x-axis is represented by the vector (1, 0) and the y-axis is represented by (0, 1).

To get the actual location of a point in a given coordinate system, we multiply each component of the point by the basis vector representing the corresponding axis. So for the point (5, 3) in the standard basis, we’d multiply 5 by (1, 0) and 3 by (0, 1). This gives us (5, 3) – the exact same point we had before.

But what if we were using a different basis? What if, say, the x and y axes were rotated by 30 degrees counter-clockwise?

In this case, the x and y components of (5, 3) represent a displacement of 5 units along the rotated x axis and 3 units along the rotated y-axis. To calculate the resulting point, we multiply the components of (5, 3) by the new basis vectors. This would give us the value of (5, 3) rotated 30 degrees around the origin!

Why? Because if you rotate a coordinate system, everything inside of it gets rotated as well:

Multiplying the components of (5, 3) by the rotated basis vectors tells you what (5, 3) in the new, rotated coordinate system resolves to in the original coordinate system with the standard basis.

This is essentially what the rotation matrix is doing – it’s rotating the x and y axes that you are using to interpret the point. When you perform the multiplication, you’re using the rotated basis vectors to calculate what the same coordinates in the rotated coordinate system represent in the original coordinate system (via displacement of x and y units along the respective basis vectors).

Let’s take a look at the rotation matrix again:

The columns of the matrix represent the x and y axes, respectively, of the rotated space. When you multiply the matrix by a point, the result is actually equivalent to multiplying each component by the corresponding basis vector (in the formula below is the rotated basis vector for the x-axis and is the rotated basis vector for the y-axis):

This is equivalent to moving 5 units in the new x-direction and 3 units in the new y-direction.

Anyway, all of this begs the question: why does the rotation matrix represent the axes of the rotated space?

This question can be answered by some simple trigonometry. Recall these trigonometric identities for right triangles:

We can apply these identities to figure out the new vectors for the x and y axes. As you can see below, the rotated vector for the x-axis is at and the rotated vector for the y-axis is at .

o is opposite, a is adjacent. The denominator is 1 because the basis vectors are unit vectors (they have a magnitude of 1).

We now know the basis vectors for the rotated coordinate system. As we saw before, these vectors become the columns of the rotation matrix. And when you multiply the rotation matrix by a vector representing a point that you want to rotate, it is equivalent to multiplying each component of that point by the corresponding basis vector (due to the properties of matrix multiplication), giving you the rotated point.

So your final matrix for rotating by 30 degrees counter-clockwise around the origin is:

Multiplying by this matrix is equivalent to interpreting the multiplied point in a rotated coordinate system. You’re basically answering the question, “what would the point (5, 3) in this rotated coordinate system represent in the standard basis?”

3D rotation works the exact same way, except you also factor in the z-axis when creating your matrix and figuring out your basis vectors. And other transformations, such as scaling and translation, can be thought of in a similar way – you’re essentially modifying the coordinate system that you’re using to interpret your positions. So scaling by 2 is the same as doubling the basis vectors, resulting in x and y vectors of (2, 0) and (0, 2),  and then interpreting all points as displacements of (x, y) units along these new axes.

Page curl is a fairly common effect in e-reader applications. It looks like this:

The effect simulates turning a page of a paperback book, usually in an interactive way where the deformation of the page is affected by the location of your finger on the screen and the direction of your drag. I’m a big fan of skeumorphism in UI’s (contrary to modern design trends), which makes this particular effect very delightful to me.

Too skeumorphic? Nahhh

From the moment I first saw this effect, I’ve been curious about the math behind it. I finally decided recently to dive in and attempt to implement it as a learning exercise. As per my usual strategy for learning a new topic for which there is likely an abundance of pre-existing literature, I turned to Google for my initial research. I came across a few writeups on the technique (including this research paper) but found them to be lacking. Most of them were extremely vague, sometimes with code that was obfuscated to an obnoxious degree. My initial reaction was that the math was simply beyond my current level of understanding. I couldn’t have been more mistaken.

I quickly grew frustrated with trying to reverse-engineer write-only code and decided to try to figure the whole thing out from scratch. It turned out to be far less complicated than I expected and didn’t involve any math beyond geometry and trig. Here’s the end result as a GLSL shader:

Here’s my implementation on shadertoy. Click and drag on the image to see the effect. I’m not going to go too deep into the code itself in the following explanation since it handles a lot of shader-related minutiae (aspect ratio, clamping, etc) and I’d rather focus on the math, but feel free to leave a comment or reach out to me if you’d like me to elaborate further on the GLSL code.

How the Effect Works

High Level Summary

The basic idea behind the effect is that the curl of the page is represented as a cylinder whose left and right edges are parallel to what I’ll call the “curl axis”. The curl axis is the crease that you would get if you folded the page along the base of the curl:

Each page (the current page and the next page) is represented as a different texture passed into the shader. In the fragment shader, the fragment’s position is analyzed to determine where it lies with respect to the cylinder and the curl axis. This information is used to determine both which texture to use and which UV coordinate on that texture to display at the fragment’s position.

This gives us two sub-problems to solve:

1. Finding the location and orientation of the cylinder.
2. Determining the texture and UV, based on #1 and the current fragment position.

Sub-Problem #1: Finding the Curl Axis and the Cylinder

To deform the quad around a cylinder, we need to know the actual location and orientation of the cylinder. This information can be represented in many forms, so it is helpful to consider how we will use this information in sub-problem #2 to determine which form we need it in. The key piece of information we will need is how far the current fragment is from the curl axis. This distance alone gives us most of the information we need to choose the correct texture and map the UV in sub-problem #2, so we don’t need to worry about finding the angle of rotation of the cylinder or anything like that. Here’s an illustration of the piece of information that we are trying to find in this step (we’ll call it d from this point onward):

d is the distance we want to find

Finding d involves a few different steps. Let’s take it one step at a time.

First, we’ll want to represent the orientation of the curl axis in a way that makes solving for d easy. To do this, we can simply look at the direction in which the user is dragging the screen. If we find the vector between the point where the user placed his finger on the screen (we’ll call this clickPos) and the point to which he has currently dragged his finger (we’ll call this dragPos), we can treat the curl axis as a line that is perpendicular to this vector and passes through the dragPos. See these images to gain an intuition for why this makes a good curl axis:

A is clickPos, B is dragPos

We’ll calculate the direction vector (dir) as follows:

vec2 dir = normalize(clickPos - dragPos);

Now we need to figure out how to represent the fragment (f) as a vector so that the dot product actually gives us the length of the projection along the vector to the curl axis. In other words, we need to find the origin point depicted below:

f is an arbitrary fragment

To find that point, we can calculate the intersection of the lefthand side of the page with the direction vector from dragPos:

vec2 origin = clamp(dragPos - dir * dragPos.x / dir.x, 0., 1.);

Basically, we’re finding out how many instances of dir fit between dragPos and the lefthand side of the screen. We then subtract by that many instances of dir, bringing us right up to that edge. The point that we arrive at is the origin point that we’re looking for. We can then find a vector representing f by simply subtracting the origin point from the fragment position.

Now that we have the fragment and the direction vector, we can take the dot product to get the distance along the direction vector.

float distOfFragmentAlongDirectionVector = dot(fragVec, dir);

To get the critical piece of information that we need for part #2 (the distance of fragment from the curl axis), we can do the following calculation:

float distOfFragmentFromCurlAxis = distOfFragmentAlongDirectionVector - distOfCurlAxisAlongDirectionVector;

Since we’re using dragPos as the curl axis position, the distance between dragPos and the origin calculated previously gives us distOfCurlAxisAlongDirectionVector. distOfFragmentAlongDirectionVector is the projected distance that we just calculated by taking the dot product. We now have both terms that we need to execute the subtraction and find the distance of the fragment from the curl axis, which is what we need to proceed to sub-problem #2.

Sub-Problem #2: Mapping the Point to the Cylinder

Having obtained the distance between the current fragment and the curl axis, we can now properly deform the image along the cylinder. There are three scenarios now that will determine how this deformation calculation is performed. These scenarios depend on the distance calculated in step 1, as well as a pre-configured radius for the page curl cylinder. I’ll quickly summarize these three scenarios below and then go into more detail:

1. The fragment is ahead of the curl axis and not within the curl radius
2. The fragment is ahead of the curl axis but within the curl radius
3. The fragment is behind the curl axis

Scenario 1: Fragment is ahead of curl axis and not within radius

Viewed from the side, this scenario looks like this:

In this case, it’s clear that the point doesn’t lie on the curl at all and is completely on the second page. So we can simply sample the texture of the second page without deforming the UV coordinates at all.

Scenario 2: Fragment is ahead of curl axis and within curl radius

If we’re ahead of the curl axis but within the curl radius, then that means we’re on the curl itself. In this case, we’re definitely on the first page but we could either be on the front side or back side of the page. Here’s what this scenario looks like:

p1 and p2 represent the possible positions on the page for the fragment. p1 is on the front of the page and p2 is on the back of the page, since it’s curled back around and since the viewer is looking down on the page.

To find the actual UV coordinates that p1 and p2 represent, we can essentially “unroll” the curl and see where p1 and p2 would be if the page was laid flat. We can unroll it by calculating the distance along the circumference from the curl axis to the point. Here’s where the geometry and trig come in.

We know that the circumference of a circle is and we know that there are radians in a circle. If we can calculate in the diagram above, then we can calculate the distance to p1 along the circumference by taking the proportion of to and multiplying this by the circumference. This works because rotating of the way around a circle is equivalent to traveling of the distance around its circumference.

We can calculate fairly easily by using some high school trigonometry. We know the H and the O in SOH CAH TOA – H is the radius, which is simply a pre-determined value. O is the distance that we found in sub-problem #1.

We can thus find the distance to p1 as follows:

The distance to p2 is calculated in a similar way, only we replace theta by :

Now that we’ve solved for d1 and d2, we can find the values of the unrolled p1 and p2 by multiplying d1 and d2 by dir and adding each product to the point where the direction vector to the fragment would’ve intersected with the curl axis (we’ll call this point linePoint – you can find it by moving the fragment position towards the curlAxis by dir * distOfFragmentFromCurlAxis, where the latter variable is the critical piece of info we found in sub-problem #1).

vec2 linePoint = fragmentPos - distOfFragmentFromCurlAxis * dir;
vec2 p1 = linePoint + dir * d1;
vec2 p2 = linePoint + dir * d2;

To determine whether to use p1 (the front side) or p2 (the back side), we simply check to see whether or not the UV coordinate at p2 is within the UV bounds of [0, 1]. If it is outside of these boundaries, then p2 doesn’t actually exist at all on the first page and we therefore use p1 as the UV coordinate. See the following diagram to gain intuition for why this makes sense:

In the above image, we look at the unrolled p2 values for several fragments. Because the page is a rectangle, some of the unrolled points will lie beyond the page boundaries. For these fragments, we’ll use p1 to get the UV coordinates rather than p2, since p2 is no longer on the page.

Scenario 3: The fragment is behind the curl axis

Here’s a visualization of scenario 3, assuming that the curled page is straight once it is behind the curl axis:

We can find the UV here using a similar technique to scenario #2 – unroll the page and see where the point would lay on the unrolled page. If the unrolled point is within the UV bounds [0, 1], then we use that as the UV. Otherwise we just use the original UV of the fragment. Finding the unrolled UV is a little bit easier in this case since we don’t have to do any trig – we can get the distance to p by just adding half of the circumference to the distance of the fragment from the curl axis (that wonderful number we calculated in the very first sub-problem).

Like the previous scenario, we ignore the UV coord for the backside of the page if it is out of bounds and instead just use the original UV.

Miscellaneous Tweaks

And there you have it. To see the nitty gritty implementation details, feel free to check out my shadertoy implementation. As far as I know, this is about as simple as this effect can get (< 50 lines of code for the whole thing). No crazy mesh deformations or insanely complex math – just a single quad with a bit of trig. There are a few things going on in there that I didn’t go over in this post, since they’re not fundamental to the primary deformation effect:

• Adding a pseudo-shadow cast by the curled page.
• Doing some clamping to make sure the page always stays attached to the “spine” of the book on the lefthand side.

I’ve pretty much satisfied my curiosity with this effect, but there are a few small additions that I may add to the shader in the future:

• Anti-aliasing of page edges.
• Turning to the previous page (rather than only the next page).

Splines are incredibly useful in game development. Anytime you want to generate a smooth curve that goes through (or “interpolates”) any arbitrary number of points, you can easily do so by using splines. Maybe you want to procedurally generate a wacky Sonic the Hedgehog-style level, complete with crazy loops and turns. Or maybe you want an enemy AI in your game to smoothly move through some waypoints. Or maybe you just want to understand how that weird pen tool in Photoshop works. The secret to all of these things is – you guessed it – splines.

Quick warning: this article assumes some degree of mathematical fluency and it’s not at all necessary to understand all of the math behind splines to be able to make use of them. If you’re on this page because you’re just looking for a simple way to generate a curve that interpolates a list of points, then scroll down to the bottom of this article for a code sample that you can plug right into your project (well, as long as it’s a C#/Unity project).

Before we dive into the math, we’ll go over some practical examples of what you can do with splines. Both of these examples are from games I worked on in the past.

Procedural Mesh Generation

This is a video from a game I worked on at a jam last year.

The tail motion is implemented via a chain of sprites connected with springs. One way to make the tail smooth would be to increase the number of joints, but this would get expensive fairly quickly. Instead, I created a curve that smoothly interpolates all of the joints and procedurally generated a mesh that conforms to that curve. Here’s the source code for the spline generation and for the tail mesh generation.

Movement Along a Curve

My friend Ed Lu also utilized splines in a game that we worked on together for Global Game Jam 2018 called Dog Walk. In Dog Walk, the player travels between neurons along axons of varying degrees of curvature. Using splines, Ed was able to ensure that the character moved smoothly along the path.

The player gliding along a curved axon in Dog Walk.

Spline Internals

So… what exactly is a spline, anyway? According to Wikipedia, it’s a “numeric function that is piecewise-defined by polynomial functions, and which possesses a high degree of smoothness at the places where the polynomial pieces connect.” I’ll elaborate since it might not make sense without some context. For our purposes, you can think of a spline as a bunch of curves that are stitched together at the endpoints in a way that allows the transitions between the individual curves to be smooth. The end result is that all of these stitched-together curves actually look like one big continuous curve. That’s where the “piecewise” part of the wikipedia definition comes in – a spline is a bunch of “pieces”, each of which is an individual curve. The “functions” that are mentioned in the definition are just the curves, each of which is defined by a function.

Hermite Curves

Before diving into how splines themselves work, it’s important to understand the individual curves that make up the splines. These curves are “third-degree polynomials specified in Hermite form”, which is kind of a mouthful, so I’m going to call them “Hermite curves”. Anyway, here’s what your typical third degree (aka cubic) polynomial looks like:

We use p to emphasize that we are calculating position – that’s what we care about when drawing our curve, since we’re plotting points with x, y, and possibly z coordinates. You’ll also notice that this formula is a function of t. Hermite curves (and other related curves, such as Bézier curves) are defined by parametric equations – both x and y are a function of some independent variable t, rather than depending on each other. To demonstrate this concept using the cubic polynomial equation from above, we define x in terms of t:

y(t) will look identical but most likely with different values for the coefficients. t is basically a value that represents how far along the curve you are, normalized between 0 and 1. t is 0 at the beginning of the curve and is 1 at the end of the curve. If you’ve ever used the lerp function of Unity/GLSL/etc, then you’ve used  t.

So what does it mean to express a cubic curve in Hermite form?  In a nutshell, it just means to provide the four Hermite values, which I will explain in more depth shortly. These four values tell you everything you need to know to generate a curve. They also give you enough information to get the coefficients of t in the standard cubic polynomial equation, and with these coefficients we have the entire formula for a polynomial curve. The reason we prefer Hermite form is because these values are more intuitive to reason about than the coefficients. You can kind of predict how a cubic curve will look if you know the Hermite values. The Hermite values are:

– the position of the starting point of the curve.
– the position of the end point of the curve
– the rate of change at p0
– the rate of change at p1

To get a bit of intuition around how this works, think of it this way – p0 and p1 give you the endpoints of the curve. v0 and v1, which represent the rate of change around p0 and p1, tell you how the curve changes as it goes away from p0 and p1. See the following diagram to get a sense for how the  Hermite values affect the final resulting curve:

A key thing to understand is that a Hermite curve is a cubic polynomial curve. Any Hermite curve can be described by the formula . It’s just a different way of expressing the same thing – you’re supplying Hermite values instead of coefficients of the powers of t. Bézier curves are also just another form of cubic curves, except instead of using Hermite values they use “control points”. All three of these forms are capable of generating the exact same set of curves. To make this more clear – cubic curves can only generate a certain type of curve. For example it’s impossible to trace the outline of Abraham Lincoln’s face with a cubic curve, since a cubic curve has a degree of three and can therefore only have two “wiggles” (just like how a parabola, which is of degree 2, can only have one “wiggle” or “turn”). Any curve that can be described by a third degree polynomial can also be described by the Hermite or Bézier form.

Moving on. Now that we know what the Hermite values are, we can actually solve for them by plugging them into the standard equation for a cubic curve. p(0) and p(1) are pretty straightforward:

To find v0 and v1, recall that the velocity is just the derivative of the position. So the velocity of a cubic curve at any value of t is:

To solve for v(0) and v(1), we just plug in 0 and 1. Easy peasy:

Now that we have the Hermite values, let’s try to actually get a formula for the curve that we can plug numbers into and get position values out of. The first thing we’ll do is solve for the coefficients of t. We already have enough information to do this. We have all four Hermite values and there are exactly four unknowns. We can find these coefficients by solving a system of equations, which I won’t go over here since the details are pretty mundane and straightforward. Anyway, you end up with these coefficients for your cubic polynomial:

Since these are the coefficients for the various powers of t, we can actually rewrite this as a matrix product (multiply it out and you’ll see that we end up with a polynomial where each of the above coefficients is multiplied by the corresponding t):

A quick side note on this – you don’t actually have to use matrix form. I’m using matrices because that’s what most of the literature on this topic uses and one of my goals is for you to be able to easily and intuitively digest that sort of material after reading this article. But if you don’t like dealing with matrices, you can just keep the coefficients in their present form and just replace p0, p1, v0, and v1 with the new values that we come up with in the next section.

Anyway, this formula is already starting to become pretty useful. It’s a formula for a Hermite curve – we can supply any two endpoints (p0 and p1) and the tangents v0 and v1 and this formula will give us points along that curve as we plug in values of t. But this isn’t what we set out to do. We don’t want to supply tangents. That’s a pain in the ass. We just want to supply a bunch of points and magically generate a curve that goes through them. So how do we do that?

Creating our Spline

Let’s start with a small example. Let’s say we have four points and we want to generate a spline that goes through all of these points:

We’re going to do this by making a bunch of Hermite curves and sticking them together. Yup, that’s right. A Catmull-Rom spline is nothing more than a bunch of cubic curves joined together at their endpoints. The first curve will go from p0 to p1. The second one will go from p1 to p2. The final one will go from p2 to p3. So you’ll apply that formula we derived above for p(t), but in three different variants, each variant with a different pair of endpoints. The formula for the first curve will have endpoint p0 and p1, the second one will have p1 and p2, and the third will have p2 and p3:

But wait, don’t we also need the tangents v0 and v1? How do we decide what those should be?

Generally, two curves joined together at some point will be smooth if the tangent at the end of the first curve (which is analogous to v1 in our formula) is equal to the tangent at the beginning of the second curve (which is analogous to v0 in our formula). So all we need to do is make sure those are equal and we’re golden. But what should those tangents actually be? A long time ago, two computer scientists named Edwin Catmull and Raphael Rom discovered that a good way to pick tangents is to just subtract the previous point from the next point. So if you’re calculating the tangent for point p2, you’d take the result of (p3 – p1) since p3 is the next point and p1 is the previous point.  It’s also common to divide these tangents by 2, not for any mathematically rigorous reason but just because doing so tends to give curves that are aesthetically pleasing. The key insight to gain from this is that the only information that you need in order to calculate tangents is the adjacent points. So you can actually rewrite the above equation in terms of four points p0 through p3, instead of in terms of p0 p1 v0 and v1. We’ll do this shortly.

Another tricky thing that comes up is that if you need to use the previous and next point to calculate the tangent, it raises the question of what the previous point of the first point will be and what the next point of the last point will be. One solution to this is to add a “ghost” start point and another ghost end point, which aren’t actually part of the curve but are used to calculate the tangent:

The ghost points are used to calculate the tangents at p0 and p3. The yellow lines are there to indicate that the tangents at the endpoints are calculated by finding the vector difference of the adjacent points (one of which is a ghost vector in both cases) and then dividing by 2.

Another common technique is to just assume some pre-determined default value for the tangent at the start and endpoints instead of using ghost endpoints. This is less flexible but slightly simpler. For the rest of this tutorial, we’ll be using the ghost endpoint technique.

Anyway, let’s modify our formula from above to reflect this previously mentioned method of calculating the tangents. Keep in mind that p0 and p1 below are not the same p0 and p1 from the previous matrix. We are simply rewriting the above matrix as a function of four arbitrary control points:

In this form, p1 is analogous to p0 in the previous matrix (it’s the first interpolated point). p2 is the same as p1 in the previous matrix (it’s the last interpolated point). The next two elements of the column matrix are the tangents, which are calculated using the adjacent points as we just discussed. We introduced two other points p0 and p3 to use for tangent calculation. The symbol is tension, which you can tweak to change how sharply your spline changes direction at the knots between curves. This is commonly set to 1. You can play around with tension values to change the appearance of your curve. Just to give you an idea of how this looks, here’s a picture of a spline with the tension bumped up from 1 to 10:

Anyway, the equation we derived above technically gives us everything we need in order to draw a curve from four points, but the equation is a bit ugly. You can actually do a bit of matrix manipulation (not shown, since it’s not really relevant) to get it into the form below, which is the form you’d probably actually want to use:

At this point, you’re basically done. Believe it or not, the formula above gives you everything you need to create a spline that goes through any arbitrarily-sized list of points (minus the first and last points in the list, which will be the ghost points). p0 and p3 are used for tangent calculation, and p1 and p2 are actually passed through by the curve.

From here, you can create a list of points that you want your spline to go through. For each piece of the spline, you can get the formula for that piece by plugging in values to the formula and getting a Hermite curve equation. For the first curve you’d plug in the first, second, third, and fourth points. This curve will go through the second and third points, using the first and fourth points just for tangent calculation. For the second curve you’d plug in the second, third, fourth, and fifth points. This curve will go through the third and fourth points, using the second and fifth points just for tangent calculation. And so on and so forth. Note that all of the points used for tangent calculation will eventually become an actual interpolated endpoint of another curve, except for the very first and last points in the list (which is why I call them “ghost” endpoints).

You can evaluate these functions at various values of t, spaced out however you like (depending on how many points you want to evaluate on that curve). As I mentioned previously, the spline won’t go through the first and last points of your list, since they’re being used just for calculating tangents, but it will go through every single other point.

Also don’t forget that since these curves are parametric, you’re actually solving for x(t) and y(t) separately. This doesn’t really complicate things though – after solving individually for x(t) and y(t), you have the coordinates you need to plot an (x, y) point in 2D space. Similarly, if you’re doing a 3D game, you could just do it again for z and get a 3D point. No biggie.

So there you go. Create a big list of as many control points as you want. Use the above template to get the actual formulas for each of the curves that make up your spline. Plug in a bunch of values of t spaced out however you like. And you’re done.

Example Code

In practice, you generally won’t want to actually put the matrices in your code. If you multiply out the matrices that we derived above, you get a very simple equation that you can then plug into a function:

public static List GenerateSpline(List points, int stepsPerCurve = 3, float tension = 1)
    {
        List result = new List();

        for (int i = 0; i < points.Count - 1; i++)
        {
            Vector3 prev = i == 0 ? points[i] : points[i - 1];
            Vector3 currStart = points[i];
            Vector3 currEnd = points[i + 1];
            Vector3 next = i == points.Count - 2 ? points[i + 1] : points[i + 2];

            for (int step = 0; step <= stepsPerCurve; step++)
            {
                float t = (float)step / stepsPerCurve;
                float tSquared = t * t;
                float tCubed = tSquared * t;

                Vector3 interpolatedPoint =
                    (-.5f * tension * tCubed + tension * tSquared - .5f * tension * t) * prev +
                    (1 + .5f * tSquared * (tension - 6) + .5f * tCubed * (4 - tension)) * currStart +
                    (.5f * tCubed * (tension - 4) + .5f * tension * t - (tension - 3) * tSquared) * currEnd +
                    (-.5f * tension * tSquared + .5f * tension * tCubed) * next;

                result.Add(interpolatedPoint);
            }
        }

        return result;
    }

Basically all this is doing is taking a list of points and returning a list of points along a spline that smoothly interpolates all of the provided points. The stepsPerCurve parameter specifies the number of points to generate between each of the original points. The greater the value, the smoother the spline. The tension parameter should remain at the default value in most cases.

References

I came across quite a few online resources while I was learning this material but these two were by far the best:

• Bézier Curves from the Ground Up – Jamie Wong – A very clear and intuitive explanation of Bézier curves with great animations. Doesn’t go into too much detail on the math, but it’ll give you a very solid intuition of how Bézier curves work and how to construct them.
• A Primer on Bézier Curves – This is an amazing resource and will give you a really solid mathematical foundation for understanding all sorts of splines and curves. I used this resource heavily during my own learning. I haven’t finished the whole book yet and I’m pretty eager to see what other cool knowledge awaits me in there.