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.
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.
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:
- Nodes. Create random points in a two-dimensional space; these will be the nodes of our graph, and where the figure can branch.
- Edges. Connect the nodes to their neighbouring nodes; set the edge weight as the Euclidean distance between the nodes.
- 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.
- 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
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 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, points.shape)) # 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): # 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]) # For each neighbour, add an entry to the adjacency matrix for neighbour, distance in zip(neighbours, distances): adjacency_matrix[i, neighbour] = distance
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.
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
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...
1. The fractal dimension of a Lichtenberg figure depends on the medium it was induced in, but I've found a few sources suggesting it's in this ballpark for planar figures, and between 2 and 2.5 for figures in a solid.