Generating Periodic Voronoi Tilings

I recently wanted to simulate fluid flowing through a periodic 2D packed-bed, and decided to procedurally generate the packing using Voronoi diagrams.

I thought that Voronoi diagrams would be a good starting point for several reasons;

• All edges are connected at both ends, meaning the fluid can flow through the domain.
• The positions of the cells can be easily set.
• The widths of the resulting flow-channels can be well controlled.

Since this was to be a lattice-Boltzmann simulation I only needed a binary image to serve as my boundary conditions, not some description of the edges & their connectivity or other mesh based data structure. This means I could use image processing techniques to achieve the desired effect. I came up with what I think is a relatively neat way of generating images that contain flow-channels of the specified width.

Before going over the method, below is a snapshot of the output from the resulting fluid simulation. I’m pretty happy with how it turned out.

Okay, let’s dive into the method.

Step 0: Set-Up the Domain

The very first step is to set up the domain and create some arrays we will use later.

import numpy as np

# Set-up the domain.
n = 1000  # The grid cell count along each axis.
x = np.linspace(0, 1, n)
y = np.linspace(0, 1, n)
xx, yy = np.meshgrid(x, y)
ii, jj = np.meshgrid(np.arange(n), np.arange(n), indexing="ij")

Step 1: Choosing Seed Points & Tiling

The next step is to choose some seed points, around which to grow our Voronoi cells. The method will work for any seed points, but we would like them to be both somewhat randomized, and approximately evenly spaced over the domain. To achieve this I used Poisson disk sampling. Conveniently SciPy has an implementation we can use without having to write our own: scipy.stats.qmc.PoissonDisk.

from scipy.stats.qmc import PoissonDisk

# Sample the seed points.
rng = np.random.default_rng(42)
poisson_disk = PoissonDisk(d=2, radius=1 / 6, rng=rng)
points = poisson_disk.fill_space()

The variable points contains an array, whose last axis represents the x & y coordinates for the points. However, we want an image to work with. What we do is create an image which is the same size as the domain, and is initially zero everywhere. Then, we set only the pixels containing our seed points to be one. We do this by converting the locations (points) into indices (idx) and updating the array using these:

# Place points into a 2d image array.
idx = (n * points).astype(np.int64)
pp = np.zeros_like(xx)
pp[idx[:, 0], idx[:, 1]] = 1

We want the final result to be periodic. Points near, say, the left-hand edge should see virtual “images” of the points that are near the right-hand edge. If we do nothing there will be edge effects because these virtual points are not present. We have two options; make sure all operators we apply are intrinsically aware of the periodicity, or we can create tiled copies of the domain. Creating a 3x3 tiled copy of the domain means that any edge effects are guaranteed not to affect the middle copy. At the end we can then chop out just the middle copy to get our final result.

# Tile so the middle chunk looks periodic.
pp = np.tile(pp, (3, 3))

Once done, we have an image which looks like the following:

This plot, and subsequent plots, show the middle tile with a small part of the extended 3x3 domain, in order to highlight the periodicity.

Step 2: Cell Labelling & Image Partitioning

The next step is to divide the domain into separate Voronoi cells for each point. This can be achieved in two (sub-)steps.

Firstly, use scipy.ndimage.label to give each seed pixel a unique ’label’ (i.e. a unique integer value).

from scipy.ndimage import label

# Give each seed pixel a unique non-zero value.
ll, *_ = label(pp, structure=np.ones((3, 3)))

Next, we use scipy.interpolate.griddata to expand these unique labels out from the seed points to those pixels which have a value of zero.

from scipy.interpolate import griddata

# For zero pixels, get the value of the nearest seed pixel.
mask = (ll != 0)
f_in = ll[mask].flatten()
x_in = np.stack([ii[mask], jj[mask]], axis=-1)
x_out = (ii, jj)
ll = griddata(x_in, f_in, x_out, method="nearest")

The idea behind this is to say that we have samples of some function’s (completely arbitrary) value at the seed points only, and we want to interpolate them to the rest of the domain. If we interpolate them using a neartest neighbour search then we get a piecewise-constant function over our domain. Moreover the constant values correspond exactly to the Voronoi cells. This is by definition; a Voronoi cell is all the points which are closer to that seed point than to any other, which is what our nearest neighbour interpolation yields.

Once we have run the interpolation we get back the following image:

Great, we can see the Voronoi cells here. However, there’s no way to get the final result we want directly from this image.

Step 3: Cell Edge Detection

To make progress we want to detect the edges of the cells. Because the pixels comprising each cell have one value, we can check if any neighbouring pixels have a different value. If so, we know we must be at a cell edge.

A quick way to check this is to calculate the minimum and maximum in a rolling 3x3 window over our image; if they differ we know we are at a cell edge.

from scipy.ndimage import maximum_filter, minimum_filter

# Identify cell edges.
ee = maximum_filter(ll, size=3) - minimum_filter(ll, size=3)
ee = (ee != 0)

This works great:

Step 4: Edge Thickening

We know that any pixels in the flow-channel must be within half of the desired width of the cell edge, which lies at the middle of the channel. This implies that for each pixel, we want to know the shortest distance to the nearest edge.

There exists an image transform to calculate just this: the Euclidean distance transform. Again SciPy has our backs, and provides a function to calculate it; scipy.ndimage.distance_transform_edt.

from scipy.ndimage import distance_transform_edt

# Calculate distance to the nearest cell edge.
dd = distance_transform_edt(1 - ee)

This yields the following plot:

Now we can threshold this distance field at the desired half-width to get the final boundary-condition image.

# Clip distance field to get channels of the desired width.
width = 40  # Units of pixels.
tt = (dd <= width / 2)

Giving

Now all that remains is to select out the central tile from our 3x3 tesselation:

tt = tt[n:2*n, n:2*n]

Extension - Noisy Boundaries

The method above works great and is completely generic over the seed points & channel widths. However, it definitely looks very artificial given every channel has exactly the same width! We can easily remedy this. Note that we could have used a spatially varying half-width, which will give channels of differing widths throughout the domain.

You can use any sort of field you like for the half-width, e.g. Perlin noise, but there are two things to note:

• High frequency components in the width-field can lead to disconnected wall boundaries.
• The field should be periodic with the same size as the domain, else we will see discontinuities at the edges.

Here I have used a simple linear combination of some plane-waves, which allows me to both avoid high-frequency components and ensure periodicity (by using integer multiples of the domain size for the wave-numbers’ components.)

# Randomize the width field.
ww = np.ones_like(xx)
ks_ = [[3, 5], [4, -1], [7, 15], [3, 1]]  # Wave-number
as_ = [0.3, 0.2, 0.1, 0.1]                # Amplitude
for k, a in zip(ks_, as_):
    ww += a * np.sin(2 * np.pi * (k[0] * xx / size_x + k[1] * yy / size_y))

ww = np.tile(ww, (3, 3)) * (width / 2)

This gives the following periodic “width-field”:

which we can then use to threshold the edge-distance field we calculated above:

tt = (dd <= ww)

This yields the final boundary condition image:

It reminds me of a Giraffe!

Extension - Spatially Extended Seeds

I mentioned above that, for the point-like seeds, we needn’t have done the double-step of placing the points into the image as pixels of one colour, only to then label and extract them to feed to the nearest-neighbour interpolation. So why do it?

It means the method will work for any binary image you choose to feed it, not just one containing single isolated seed pixels. Even concave shapes are fine.

For example, here are the steps of the method applied to some random squiggles I drew with my mouse. All of the labelling, edge detection and thickening work just fine. Nice!

Code

The code that produced the plots for this article is available here, and the code for the fluid simulation is available on GitHub.