I found a cool Voronoi variant that gives a nice fractal boundaries to each Voronoi cell, that is very efficient to compute. It seems like a good fit for coastlines, which are often similarly fractal in nature.
https://www.shadertoy.com/view/sfKSDw
Let me explain.
“Jittered” Voronoi is where you start with an infinite grid, and then for each square in the grid, you pick a random point called a site in that square. Then you make a Voronoi diagram from those sites. It’s the same thing as Worley noise, except instead of outputting the distances, the output is the choice of point you are nearest to, partitioning the plane.
It has the nice property that it’s very easy to compute on a per-point basis: for a point p you just find the nearest 25 grid cells, find each of their sites (via a hashed pseudo-random number generator), and then find which is nearest to your point. No need to actually construct the Voronoi diagram, or worry about the infinite size of the plane.
Here’s our new procedure, described mathematically:
1) Again, start with an infinite grid of squares, and pick a site for each via PRNG. These are the root, layer 0, sites.
2) Make a new grid, with squares half the size, and pick layer 1 sites for them. For each layer 1 site, find its parent site in layer 0, which is defined as the nearest site.
3) Repeat step 2 indefinitely. Each time, halve the grid size, and find parents in the layer above. The partition each site belongs to is the root site you get to when tracing the parents upwards.
As you recurse, you fill the plane with sites, which can be extended to an almost-everywhere partition of the plane.
But how do we compute this in practice? We cannot have an infinite grid, nor infinite recursion. Well the infinite grid can be dealt with like regular jittered Voronoi – for a given point p, at every level, there’s only a fixed number of cells that can actually be relevant. As for infinite recursion, the obvious solution is to just do a finite number of iterations (depending on zoom level), and make a Voronoi diagram of the final layer, which will be a good approximation of the fractal shape.
Let’s go over that in code. I code golfed the loop termination condition a bit.
from functools import cache
from math import sqrt
# Jittered Voronoi for each layer
def cell_size(layer):
return 2 ** (-layer)
@cache
def site(layer, cell):
"""
Return the site belonging to an integer grid cell (x, y).
hash2 must deterministically return two pseudorandom numbers in [0, 1),
based on (layer, cell).
"""
s = cell_size(layer)
offset = hash2(layer, cell)
return s * (cell + offset)
def nearest_site(layer, p):
"""
Find the layer-i site nearest to p.
"""
s = cell_size(layer)
centre_cell = floor(p / s)
best_cell = None
best_site = None
best_distance = infinity
# Searching the nearest 5x5 cells is sufficient:
# cells further away cannot contain the nearest site.
for dx in range(-2, 3):
for dy in range(-2, 3):
cell = centre_cell + (dx, dy)
q = site(layer, cell)
d = squared_distance(p, q)
if d < best_distance:
best_distance = d
best_cell = cell
best_site= q
return best_cell, best_site
# Parent/root
@cache
def parent(layer, cell):
"""
Return the parent of a site in layer `layer`.
"""
return nearest_site(layer - 1, site(layer, cell))[0]
@cache
def root(layer, cell):
"""
Return the layer 0 cell associated with a cell in a layer.
"""
while layer > 0:
cell = parent(layer, cell)
layer -= 1
return cell
def partition(p, depth):
"""
Finite-depth approximation to the fractal partition.
"""
cell, _ = nearest_site(depth, p)
return root(depth, cell)One variant is to have no fixed depth, but an early stopping condition. If the point p happens to be particularly near a site at layer i, we can stop, and trace parents from there. That’s sufficiently close to a site, all later layers will agree with it.
def nearest_site_safe(layer, p):
"""
Find the layer-i site nearest to p.
Also determine whether all possible layer-i
ancestors of finer sites have the same root. If so, further recursion
cannot change the answer.
"""
s = cell_size(layer)
centre_cell = floor(p / s)
best_cell = None
best_site = None
best_distance = infinity
# A site selected at an arbitrarily deeper layer can have a layer-i
# ancestor at most this far from p.
influence_radius = 2 * sqrt(2) * s
possible_roots = set()
# ±2 is enough just to find the nearest site.
#
# When checking whether it is safe to stop we use ±3, which contains
# every cell whose site could lie within influence_radius.
radius = 3 if check_safe else 2
for dx in range(-radius, radius + 1):
for dy in range(-radius, radius + 1):
cell = centre_cell + (dx, dy)
q = site(layer, cell)
d = squared_distance(p, q)
if d < best_distance:
best_distance = d
best_cell = cell
best_site = q
if check_safe and d <= influence_radius ** 2:
possible_roots.add(root(layer, cell))
safe = len(possible_roots) == 1
return best_cell, best_site, safe
def partition_adaptive(p, max_depth):
"""
Compute the partition containing p, stopping early whenever it is
guaranteed that further recursion cannot change the result.
If no such guarantee is reached, use max_depth as a finite
approximation.
"""
for layer in range(max_depth + 1):
cell, _, safe = nearest_site(layer, p, check_safe=True)
if safe:
break
return root(layer, cell)