Diffusion Limited Aggregation is a growing object consisting of particles that diffuse and aggregate to the object, much like how flocks of snow are created. Particles are introduced one after the other and made to undergo a random walk. Eventually, a newly introduced particle reaches an existing particle and remains there. This clustering of particles leads to a formation of aggregates of particles. This theory, proposed by T.A. Witten Jr. and L.M. Sander in 1981 is applicable to aggregation in any system where diffusion is the primary means of transport in the system. (André Offringa's report).
Another interesting description (and my favourite) involves a city square surrounded by taverns. Drunks leave the taverns and stagger randomly around the square until they finally trip over one their insensate companions at which time, lulled by the sounds of peaceful snoring, they lie down and fall asleep. The tendril-like structure is an aerial view of the sleeping crowd in the morning (Paul Bourke's blog).
P.S. Code for the simulations is available in this repository.
This article is the first of two parts on DLA. It's mostly about my implementation of a 2D DLA simulation with a point attractor, and an optimization I use to make these simulations faster. I plan to follow this up with another one on the stickiness factor estimation (more on stickiness below).
'Brute force' DLA simulations
Generating these images is pretty straightforward.
- We begin by taking a blank square grid (all white pixels) with one black pixel in the centre
- Then we introduce a new particle (a black pixel) randomly at any point on the boundary.
- This new particle undergoes Brownian motion till it comes in contact with another black pixel on the canvas.
- As soon as the new pixel comes in contact with another black pixel, it stops moving and stays there forever.
The pixel undergoing Brownian motion can move to any of its eight neighbours with equal probability. This version of the DLA which we discussed has a point attractor (i.e. new particles get attracted towards a central point). This attractor itself can have a different locus (it could be a line or a circle), but we just focus on point attractors for the moment.
Another attribute which I incorporate in this model is called 'stickiness'. As soon as the new particle undergoing Brownian motion comes in contact with an existing particle at rest, instead of stopping, it stops with some pre-decided probability (stickiness factor) k.
A python notebook of this simulation can be found here. It has a DLA class which can be initialised with the size of the square grid and the stickiness factor k. The addPoint method can be used to add any number of particles to the canvas, the printState method can be used to print the grid out (uses matplotlib.pyplot.imshow).
Making the simulations run faster
As you would have guessed, adding a new particle is a very time-consuming process, especially as the number of particles increases. To accommodate a large number of particles, we need a larger grid. As a result, the area which a new particle explores before coming to rest increases. Of course, an easy way to quicken the whole process is to not use Python. Let's discuss a more mathematical approach.
At each iteration of DLA, we introduce a particle \(C_1\) at a randomly selected position \((x, y)\) on the border of our m x m canvas. This particle follows Brownian motion in the canvas until it meets the termination criterion. As m increases, the area in which \(C_1\) can move also increases, making simulations slow. In order to quicken the process, we would like to decrease the search space. Intuitively, we would like to spawn as close as possible to the current group of particles in the canvas.
We can do this by finding a group of points \(L_2\) such that the probability that a point that spawned on the square border reaches any of the points in \(L_2\) is the same. Formally, this would be a locus such that any two points \(l_{21}\) and \(l_{22}\) in it satisfy the condition
\[W_{reach}(a, b) = \frac{n(a, b)}{\sum_{c \in L_1}n(c, b)}\]
\(L_1\) is the set of points which make up the square boundary. \(W_{reach}(a, b)\) is the ratio of the number of ways of reaching \(b\) from \(a\) to the total number of ways of reaching \(b\) from all points in \(L_1\). We take a ratio since there are infinitely many ways of reaching \(b\) from any point in \(L_1\).
If the area covered by locus \(L_2\) is less than the area of the square boundary, then we can hope to see a faster DLA. Intuitively it seems \(L_2\) can be a circle. We shall be proving it below.
Minimum bounding circle
\(W_{total}\) of a point \(C_2\) on the circle is calculated by summing \(W_{reach}(c, C_2)\) for all points \(c\) on the square. \[ W_{total}(C_2) = \sum_{c \in side1}W_{reach}(c, C_2) \quad + \sum_{c \in side2}W_{reach}(c, C_2) \] \[ \quad \quad \quad \quad \; \; + \sum_{c \in side3}W_{reach}(c, C_2) \quad + \sum_{c \in side4}W_{reach}(c, C_2) \]
Now we shall show that \(W_{total}\) is equal for all points on the circle. There are two assumptions we make to simplify our analysis
- Instead of considering discrete coordinates, we shall consider all points on the 2D plane.
- We define \(W_{reach}(a, b) = k(|a_x - b_x| + |a_y - b_y|)\), where \(k\) is a constant.
Using our two assumptions,
Similarly,
Therefore,
\[ W_{total}(C_2) = \sum_{c \in side1}W_{reach}(c, C_2) \quad + \sum_{c \in side2}W_{reach}(c, C_2) \] \[ \quad \quad \quad \quad \; \; + \sum_{c \in side3}W_{reach}(c, C_2) \quad + \sum_{c \in side4}W_{reach}(c, C_2) \] \[ = 4\frac{S^2}{2} - kSR\cos\theta + kSR\sin\theta + kSR\cos\theta - kSR\sin\theta \] \[ \quad \quad \quad \quad + 2R^2\cos^2\theta + 2R^2\sin^2\theta + 2R^2\cos^2\theta + 2R^2\sin^2\theta \] \[ = 2S^2 + 4R^2 \]
Clearly \(W_{total}(C_2)\) is a constant, for every \(C_2\) on the circle with radius \(R\). Hence the probability of reaching any point on a circle with radius \(R\) is the same. As a result, we can choose to spawn at any point on the circle with a radius \(R\) instead of the square boundary.
The only constraint is that all the points outside the circle should be free. And obviously, we would like \(R\) to be as small as possible. Hence while spawning a new point, we chose the smallest circle that encloses the current image (minimum bounding circle).
A python notebook of DLA simulations using this optimization can be found here. Unfortunately I haven't been able to do any kind of benchmarking comparing the running times, however, I was able to go up to 10,000 particles (compared to 500 in the previous version) in a reasonable amount of time.