An introduction to Voronoi

The one dimensional case

Voronoy tesselation is named after the russian mathmatician Georgy Voronoy and is also often written as Voronoi. As a tesselation, it is a method of dividing a given space into geometric shapes. In this case convex polygons. Often this is done in two dimensions, though it generalizes to as many dimensions as you want. It should be noted that the Voronoi tesselation is dual to the delaunay triangulation which divides the space into trianlges, hence triangulation. A dual graph is what you get if you convert each polygon into a vertex (i.e. dot) and then connect the dots of adjacent polygons. Especially triangulations are interesting for computer graphics since they allow arbitrarily complex shapes while only relying on a common geometric primitive: the triangle. Which intern makes it not only much easier to implement various algorithms (such as collision detection) but also allows for much greater optimization.

Even though in the technical sense Voronoi only refers to the Voronoi tesselation it is also often used to refer to other algorithms using the same underlying principle. And the core of Voronoi is rather simple. Let's play around with the simplest case: one dimension and only two points. Move your mouse around and see how the graph changes.


Falloff Steepness
Line Density

Note that the green line indicates the y-axis and red the x axis. There are no conrete numbers on the axes but for demonstration puposes this makes no difference. You probably also tried various settings on the right but not all of them are as easy to understand. So let me begin by explaining what you are seeing in the first place.

There are two vertical lines. One in the middle and one beneath your cursor. For every point along the x-axis we check the distance to the closest line and plot that either as a graph or as a grayscale value. The distance metric here determines how we measure that distance.

Euclidean is probably the one you are already familiar with. Just subtract your point from the other point. So 5 would be 2 units apart from 3 because 5-3=2. Since we don't care about in which direction the distance goes we simply ignore the sign. Which would be noted as abs(5-3). Or more formally: sqrt((5-3)²) as the square root is defined to always be positive. You might already suspect how this works in two dimensions. If we have two points a and b, then their distance is sqrt((a.x - b.x)²+(a.y - b.y)²). Sidenote: a.x notates the x-coordinate of the point a etc. See how this is exacly the same just with two dimensions? And since squaring always gives a positive result we also always get a positive distance.
Now Manhattan seems just the same as Euclidean. And that is because for one dimension, it is. But as we soon will see this changes if we go into two dimensions.
Hyperbolic is far to complicated to explain here but just play around and you might see what kind of effects you get.
At last: Spherical. This is just if you take two points on the globe and draw the shortest line between them. And if you go in a straight line on the globe you end up where you started. So it should be no suprise that if your mouse leaves the right edge it comes right back in at the left side.

Now that we got a distance we need a way to represent that graphically. Here in the one dimensional case that was easy by just printing a graph. However as we get to two dimensions we would need to show a three dimensional graph which is not as intuitive to understand. So I made various ways to display that distance.
The default one (absolute) is just that. If the graph reaches 1 (the top) we are directly at the line/point. If we go further away the graph will go down as much as we are away from the line.
Absolue inverted is just 1 - absolute. So if the graph is at y=0 the distance is also 0. If the graph is at 1 we are one unit apart. You get the idea.
Inverse will take 1 / distance. When you change the falloff steepness you change how quickly the graph goes towards 0. When we get to two dimensions you will notice that this looks almost how light intensity works.
And finally Equidistant Lines. This looks really weird in one dimension but I think you will understand immediatly when we get to the two dimensional case.

One last technical note here: I only draw the line for each pixel in the x-axis. This means if the graph becomes steep the points will seem detached and not as a line anymore. To circumvent that I would need to do multisampling. Basically check how the adjacent pixels look and factor that in. Which then smoothes out the curve. This is also why multisampling makes games beautiful but require much more computaional power as for every additional sample you take you have to do squared many computations. Meaning if I only check two neighbours I need to do 4 times the amount of calculation. If you want to play around more and even zoom in, give yourself the challenge and try to create the graphs at desmos, an online graph calculator, which does the multi-sampling for you.

Now I've been talking a lot about the case in two dimensions, so without further ado, here it is.

Two dimensions

As always move your mouse around and see what happens.

Falloff Steepness
Line Density

I hope you have really played around because there is already some crazy stuff going on. Lets begin with the distance metrics again.

When using Euclidean distance this time you will notice it forms circles. That is because all the points on a circle are the same distance to the center. In fact this is the very definition of a circle: The set of points that have the exact distance of the radius to the center.
Manhattan looks very different now. That is because it doesn't add the distances under the root. So it is caluclated as abs(a.x- b.x) + abs(a.y - b.y). Now pure fomulas aren't intuitive to most so let's take a real-life example. Imagine you're in a big american city build on a grid. If you ask someone how far something is they'll probably give you a distance in blocks. Which is not only much easier to keep track of but also makes sense, since you can't go through buildings. And this is precicely where the name comes from: the borough of Manhattan. So instead of taking the shortest direct path you take the difference in each axis and then add them up.
There isn't much to say about hyperbolic other than it looks really weird. That is also the main reason I added it.
When you select spherical you might already have seen this kind of pattern somewhere before. If you ever saw a map of the day/night border on the globe you will now see where that comes from. This doesn't have anything to do with spheres per se but rather how we present spheres on maps. In fact there is no way to draw a map of a sphere without stretching or displacement. You can see this when you enable axes and then take a look at the y-axis. As we get closer to the poles the axis spreads out and eventually the whole top is just a single point on the globe. If you are interested in further reading: I used a simple equirectangular projection.

This is also a fantastic time to revisit the equidistant lines. These are exactly the same as the isobars you see on weather maps or the height lines on mountain maps. Each ring represents the same distance to the closest point.
If you play around like I do, you might have pulled the line density all the way down and noticed these crazy patterns come in. There patterns are called Moiré patterns. Though they have nothing to to with Voronoi, they do look cool so I decided to keep them in.
Equidistant lines also allow us to inspect the internal structure of each of the two cells. Cells? Remember: all we do is seperate the plane into two parts. That's what tesselation is. Now you might say that the two cells don't exacly look like polygons right now. And in a sense that is true but see what happens when we add more points.

Many points at once

As we have many points at once, click and drag each point around.

Falloff Steepness
Line Density

Mix UV

Now that are polygons.
Also make sure to check out the "use colour" option. As you can see we can also identify in which cell we are.

At this point you should have gotten the gist of what is going on here. Though if you don't quite understand how hyperbolic or spherical distance work that is fine. The main point is: Voronoi doesn't care about which distance metric you use as long as it is differentiable. I.e. when you travel a small amount, the distance should also only differ by a small amount.

Though there isn't much to explain there still is a few things to mention:

  1. When you select Euclidean distance with "absolute inverted" and no color you get a common variation of Worley-Noise.
  2. The slider Mix UV lets you see the coordinates used. This is done with a common technique in shaders where you color-code information you want to see. In this case the more green a pixel is the more it is away from the y-axes. Same with red but using the x-axis.
  3. Since for every additional point you have to check the distance to each point this scales very poorly. For example if I double the amount of points I have to do four times as many checks. In computer science you'd say it has a complexity of O(n²) (more on what that means later).
  4. Since we divide the space into cells we always get hard borders. Though sometimes it is desired to have fuzzy borders for artistic purposes. You can achieve that many different ways and through different means. A fantastic (and shader oriented) explaination by Inigo Quilez can be found here.
  5. Voronoi usually only uses the closest point. But you could also take the second closest or even the furthest point with interesting effects.
  6. With the same train of thought you can also make certain points "stronger" in that they have a bigger influence. Those are then called Weighted Voronoi Diagram.
  7. You can also calculate the exact border and have certain effects applied to it too. Though the computations are more complex. A description by Inigo Quilez can be found here

Optimization

When you use many points at once, Vornoi tesselation can become quite expensive to calculate. However there is an easy trick to solve that with an added benefit. Beware that if you increase the scale the FPS (Frames per second) can drop quite fast, so especially on weaker devices don't go to the biggest scale immediately.


Falloff Steepness
Line Density

Mix UV

Scale
Animation Speed

The basic idea is to devide the plane into a grid. You can see the grid if you enable axes. Each grid box then only has one point. This has two advantages:
Firstly, all the dots are roughly equally spaced, so they don't bunch up together as much and you always get a nice distribution of points.
Secondly (and this is where the optimization comes in), we only need to check the eight surrounding boxes of a point to check for the closest neighbour. As all other dots must be further away the closest dot can only be in the surrounding grid cells.
So instead of checking every dot, we only need to perform eight checks. In theorecial computer science you would say it changes the complexity from O(n²) to O(n). This is called Big O notation and simply puts it mathematically what we already discovered. If we double the amount of points without the trick, then we need to do four (2²=4) times as many calculations, whereas with the trick we only need to do twice the amount. When we need a few points this may be acceptable but what if we need 30 times as many. Then suddenly we need 30²=900 times the calculations. So without the optimization Voronoi becomes really expensive really fast.

This may bring some advantages but as you maybe also noticed, comes with a few drawbacks. For exampe a sphere has only limited space so we can not arbitraliy add more cells as needed. As a result the border become jagged as each cell has it's own sphere it lives on. The same goes for hyperbolic space. Even though there are ways around this issue they become even more complex and aren't as trivial to implement.

As always feel free to play around with the various paramters as there are many cool effects you can achieve. One reccomendation from myself is euclidean distance with equidistant lines of density ~0.8. No colour, no axes and no UV mix and an animation speed of around 0.61. Best viewed at a big scale (of course with optimization).

But how do I do that in a GLSL shader?

Depending on your desired effect or purpose (You can also use Voronoi in a Vertex Shader!) your shader will be implemented differently. But I am going over the muliple points demo and how to implement the optimization. Though I try to keep it as simple as possible I will not explain how to write GLSL in general so you should be familiar with the basics of GLSL but I think you don't need to know every little detail to understand what the code does.

This is an implementation of the multi shader in a slightly modified manner. I will be including code snippets but if you want to see the whole code or just play around with it you can see it at the Shadertoy website. Don't worry about modifying the code at the Shadertoy Website as your canges will only be local to your machine.

As you might have noticed, you can only move one point. This is due to shadertoy only hosting pure shaders. Shaders themselves can't remember anything so they need a companion programm to tell them the info they need. On here I do that with Javascript and P5 (Processing) but on shadertoy I can not do that. However Shadertoy provides the current mouse position when the left mouse button is pressed so I use that instead. If you implement Voronoi yourself, you can pass Uniforms to your shader instead. This also applies to the constants.

Speaking of which, here they are. These are used to configure the internal workings and you are already familiar with them through the examples above.


			const float uvMix = 0.;
			const int distanceMetric = EUCLIDEAN;
			const int distanceVisualization = ABSOLUTE;
			const float steepness = .2;
			const float lineDensity = .5;
			const bool showAxes = true;
			const bool showColor = true;
		

Distance plays a crucial part in Voronoi so I defined myself a function, that calculates the distance beween two points p1 and p2 in the given distance metric. Most implementations will only use Euclidean distance and therefore just do length(p1 - p2). Since the distance will need to be calculated many times, there is usually some low-level optimization done which I skipped for readability purposes. However you will notice that generally it is easier to calculate the distance in Euclidean space rather than spherical or hyperbolic.


			float distanceFunc(vec2 p1, vec2 p2) {
				switch (distanceMetric) {
				case EUCLIDEAN:
					return length(p1 - p2);
				case MANHATTAN:
					return abs(p1.x - p2.x) + abs(p1.y - p2.y);
				case HYPERBOLIC:
					//Distance in the Poincaré disk model https://en.wikipedia.org/wiki/Poincar%C3%A9_disk_model#Distance
					return acosh(1.
						+ 2. * pow(length(p1 - p2), 2.)
						/ (
							(1.-pow(length(p1), 2.))
							*(1.-pow(length(p2), 2.))
						)
					);
				case SPHERICAL:
					//Shift and scale coordinates to align sphere with uniform coordinate space
					p1 = (p1 + 1.) * PI;
					p2 = (p2 + 1.) * PI;
					p1.y = (p1.y + PI) / 2.;
					p2.y = (p2.y + PI) / 2.;
					//Distance between two points on a sphere https://en.wikipedia.org/wiki/Great-circle_distance
					return return asin(sqrt(haversine(abs(p1.y - p2.y))+cos(p1.y)*cos(p2.y)*haversine(abs(p1.x - p2.x)))) / PI;;
				}
			}
		
This is a whole lot of code for such a simple task. But remember that this implements various distance metrics and for a single effect you probably will only need one. Also as we will soon see, we will only call the function at one place so in most cases you wouldn't even make a seperate function.

This is already all we need. Just a distance function. So let's go ahead and write the tesselation.


			//Cell id of the voronoi cell (or point if you want)
			float cellID;

			//Set the minimum distance to something impossibly large
			float minDist = 1000.0;

			//Find point with the smallest distance
			for (int i = 0; i < POINTS; i++) {
				//Distance from pixel to point
				float distanceToPoint = distanceFunc(uv, point(i));
				if (distanceToPoint < minDist) {
					minDist = distanceToPoint;
					//Normalize point/cell id to range from 0-1 
					cellID = float(i) / float(POINTS);
				} 
			}
		
In case we want to colour each Voronoi cell afterwards we will need to differentiate them. So we define an ID for the closest point and of course the distance to that point. Since we compare each point to the current smallest distance we need to initialize the minimal distance to something impossibly large.
Now we go through every point. I used a function called point to get the current point but you can use an array or any other way to get each point. Then we calculate the distance from the pixel (in OpenGL called fragment) to the point. However we don't know where pixel is in relation to the points. To find that out you usually transform each pixel into a world-space coordinate (or each world-coordinate into a pixel position). There are many ways to do this especially in 3D but for us a simple shift and scale is sufficient. This also takes care of different resolutions and therefore normalizes the pixel coordinate to a given space. This is then also called the Uniform Vector or UV in short. I left that part out as the actual calculation is not that important; just be aware that our world coordinates are called UV.
Then as a last step we check if the current point is closer than the current closest point and if so, set the new clostest distance and update the cell ID. I divide the cell ID by the total amount of points so that the cell ID is in the range of 0 to 1, which makes calculations later on easier.

That at least is the naive implementation that scales very badly. So how do you go about the optimization. Well quite easy, take a look:


			//Divide world into 1x1 grid cells
			vec2 gridCell = floor(uv);
			vec2 gridUV = fract(uv);

			float cellID;
			float minDist = 1000.0;
			//Find point with the smallest distance of the
			// surrounding 3x3 grid
			for (int xOffset = -1; xOffset <= 1; xOffset++) {
				for (int yOffset = -1; yOffset <= 1; yOffset++) {
					vec2 offset = vec2(xOffset, yOffset);
					vec2 neighbourGridCell = gridCell + offset;
					float distanceToPoint = distanceFunc(gridUV, point(neighbourGridCell) + offset);

					if (distanceToPoint < minDist) {
						minDist = distanceToPoint;
						cellID = float(random(neighbourGridCell));
					} 
				}
			}
		
This is what it looks like in a shader. I've also taken the liberty of adding a small animation in the point function:
Firstly we split the world into a grid. We do that by taking the fractional part (i.e. the part after the decimal point) as our new UV and then the floor (i.e. the value before the decimal point) as a vector to get to that cell. The concept to create the grid may not be intuitive but to explain how and why it works would make this tutorial unneccesarily long.
The actual calculation is rather simple again. We take a three by three field of cells and check for each of the neighbours (and of course the current cell) which has the closest point. As we need (theoretically) infinite points I randomly generate the coordinates for each point in the point function. Only one thing left to explain: The points aren't a simple list of points anymore but instead infinite in both axes, so we can no longer get a unique identifier for each Voronoi cell. Though there is a way to cheat. We generate a random value based on which grid cell we are in and use that instead. Since each grid cells has a unique position and ideally the random function never spits out the same value for two different inputs we can assign it as a pseudo-identifier. That is, it is very unlikely for a different cell to have the same identifier and even if so, as long as it is not closeby no one would notice.


And that's it already. Well not quite. We only get a distance and a cell ID at the end. What do we do now? That is a good question as your creativity is basically unlimited. I used the two values to create the interactive sections above but you can use them for any purpose. CGMatter used it to create procedural crystals. Flopine made lava inside a volcano at the Revision Shader Showdown 2020 Final. You can see the basis at 6:42 and the full effect at 29:21. And there is so much more you can do. Just take a look at these.

Now I just want to leave you with a really bare-bones implementation that you can use wherever you want. This uses the grid approach with the optimization and generates the points on-the-fly randomly using the random function. Do whatever you want with it but if you make something awsome toot me at pixdigit@layer8.space


			vec2 random(vec2 seed) {
				float x = fract(sin(dot(seed.yx ,vec2(12.9898,78.233))) * 43758.5453);
				float y = fract(sin(dot(seed.xy ,vec2(19.1235,57.764))) * 82653.5789);
				return vec2(x,y);
			}

			vec2 voronoi(vec2 uv) {
				//Divide world into 1x1 grid cells
				vec2 gridCell = floor(uv);
				vec2 gridUV = fract(uv);

				float cellID;
				float minDist = 1000.0;
				//Find point with the smallest distance of the
				// surrounding 3x3 grid
				for (int xOffset = -1; xOffset <= 1; xOffset++) {
					for (int yOffset = -1; yOffset <= 1; yOffset++) {
						vec2 offset = vec2(xOffset, yOffset);
						vec2 neighbourGridCell = gridCell + offset;
						float distanceToPoint = length(gridUV - (random(neighbourGridCell) + offset));

						if (distanceToPoint < minDist) {
							minDist = distanceToPoint;
							cellID = float(random(neighbourGridCell));
						}
					}
				}
				return vec2(minDist, cellID);
			}