Growing a lightning tree

There's something deeply beautiful about the fractal, branching processes in Lichtenberg figures. They're reminiscent of lightning, fern fronds, and river networks viewed from above, all things that show some degree of fractal geometry.

These tree-like patterns are caused by subjecting an insulating medium to high voltages to the point where dielectric breakdown occurs. Lichtenberg figures have fractal dimensions varying between 1.5 and 2 for those formed on a two-dimensional surface.^[1]

Short of getting on my hands on a particle accelerator, I wanted to try producing something like these figures: something that exhibited elaborate branching and filled a two-dimensional space to some extent, despite being a one-dimensional object.

A quick-and-dirty physics approach

We can simulate Lichtenberg figures, or something that looks like them, in a variety of ways. A pretty standard approach is that of simulating a diffusion-limited aggregation process, where we simulate particles randomly walking around until they meet a "seed" or aggregate, at which point they become part of the static aggregate. A very minimal implementation may look something like this.

# Set the size our figure will be
fig_size = 301

# Create an empty array where our figure will form
lichtenberg = np.zeros((fig_size, fig_size), dtype=bool)

# Set a "seed" near the bottom middle of the figure
lichtenberg[-10:, int(np.ceil(fig_size / 2))] = 1


def particle_is_adjacent_to_figure(lichtenberg, particle):
    
    # Determine the particle's neighbourhood
    min_x = (particle[0] - 1) % fig_size
    min_y = (particle[1] - 1) % fig_size
    max_x = (particle[0] + 2) % fig_size
    max_y = (particle[1] + 2) % fig_size
    
    # If there are any non-zeros here, we're adjacent
    return lichtenberg[min_x:max_x, min_y:max_y].sum() > 0


# Let's simulate 25,000 random particles, one at a time
for _ in trange(25000):

    # Place the particle randomly on the figure
    particle = np.random.randint(0, fig_size, size=2)
    
    # Have it perform a random walk until it comes into contact with the figure
    while not particle_is_adjacent_to_figure(lichtenberg, particle):
        particle = (particle + np.random.randint(-1, 2, size=2)) % fig_size
        
    # When the particle comes adjacent to the figure, it becomes part of it
    lichtenberg[particle[0], particle[1]] = 1

This gives us a result that looks something vaguely like a Lichtenberg figure.

However, this approach is slow – very slow. This code is not optimised in any way, to keep it as clear and minimal as possible, but it still takes several minutes to generate a fairly small figure – this one took an hour, at 301 by 301 pixels. A key limitation is that we can't simulate very many particles undergoing Brownian motion at once, because this would change the dynamics of the diffusion, which is central to the formation of these patterns.

Walking along branches

Instead, we can turn to graph theory to build something similar (although far less representative of the underlying physical process). Here's the general procedure we'll use:

1. Nodes. Create random points in a two-dimensional space; these will be the nodes of our graph, and where the figure can branch.
2. Edges. Connect the nodes to their neighbouring nodes; set the edge weight as the Euclidean distance between the nodes.
3. Paths. Fix one node as the root of the tree; then, select a random set of other nodes, and find the shortest paths between each of these nodes and the root.
4. Collect. Aggregate all of these shortest paths together, and visualise.

Let's go over each piece in turn.

from itertools import chain

import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
from scipy.sparse import lil_matrix
from scipy.spatial.distance import cdist
from scipy.sparse.csgraph import shortest_path
import networkx as nx

Nodes

For the nodes, we'll use uniformly-distributed points in 2D.

# Set the number of points in the space
n_points = 30000

# Generate this many randomly-distributed points
# between 0 and 1 in two dimensions
points = np.random.uniform(0, 1, size=(n_points, 2))

While we're here, we'll find the node we can call the origin, or root of the tree.

# Find the point nearest the point (0.5, 0)
origin = np.argmin(np.sum((points - [0.5, 0]) ** 2, axis=1))

Edges

Edges go between neighbouring nodes, but how do you define that ? I opted to use a triangulation to determine which nodes should have edges between them. I then iterated over all points and, for each neighbour, added an entry into a sparse adjacency matrix with the distance between those points.

# First, instantiate a sparse matrix
adjacency_matrix = lil_matrix((points.shape[0], points.shape[0]))

# Triangulate between our points
tri = Delaunay(points)

# Which nodes are connected to which ?
indptr, indices = tri.vertex_neighbor_vertices

# Iterate over each node
for i in range(points.shape[0]):

    # Grab the neighbours of this node
    neighbours = indices[indptr[i] : indptr[i+1]]

    # Calculate the distance between this node and its neighbours
    distances = cdist(points[i].reshape(1, -1), points[neighbours])[0]

    # For each neighbour, add an entry to the adjacency matrix
    for neighbour, distance in zip(neighbours, distances):
        adjacency_matrix[i, neighbour] = distance

Paths

The shortest path between any two nodes is a very well-studied problem in graph theory, and scipy offers a number of implementations of algorithms designed to find it. However, the perhaps not-optimally-named shortest_path function returns the length of the shortest path, not the path itself, and so we have to build it ourselves using the predecessors of each node along the path (which scipy does kindly provide).

# Find the predecessors of every node, along
# every shortest path from the origin
_, predecessor = shortest_path(
    adjacency_matrix, 
    return_predecessors=True,
    indices=origin
)

# Then, we can build a path in this manner
def build_path(predecessor, node):
    
    # Start at the node; move towards the root
    path = [node]

    # scipy denotes "no predecessors" as -9999
    # so while there are predecessors to the node
    # we're at, add them to the path list
    while predecessor[path[-1]] != -9999:
        path.append(predecessor[path[-1]])

    return path

At this point, we can find the shortest path between any node and our root node. However, which nodes are we interested in ? On one hand, we note that many paths in the Lichtenberg figures make it to the edges of the medium, or close. On the other hand, there are many paths that follow major channels and bifurcate off of them, never reaching the edges.

Collect

Our choice of nodes dictates the general morphology of the end result. We can select only nodes on the edges of the space, or we can select nodes randomly. We can even do both and combine the two, for some interesting results, visualised in networkx.

I think the combination is best, but it's not quite what we're after. It's perhaps a little too homogeneous, and it doesn't really give us a few paths that dominate, from which many smaller branches bifurcate. In an attempt at replicating some of this behaviour, I tried finding a small number of paths all the way to the edge of the space, and then finding the nodes that neighboured those on those paths; then I added their neighbours, and so on.

Despite not having the major channels we were hoping for, this seems a lot more natural! It's worth remembering that those major channels are electrically-conductive trajectories in an otherwise insulating medium (due to dielectric breakdown), and that branching is a dynamic process; by simply creating points randomly and connecting them to their neighbours, we're going to lose some of that behaviour.

Lichtenberg figures are fascinating, and there's a lot of opportunity to tweak this algorithm, from the distribution of points (Poisson point process?) to how edges are defined (with probability decreasing with distance?), to building up the branching structures. This could be a rabbit hole...