Monday, August 3, 2020

Fourier-Mellin Fractals

The Fourier-Mellin Transform (FMT) seems to be fundamental in encoding 2D fractals.

Here are example kernels of which the FMT composes all shapes:

Notice that these are standard Fourier transformation planar waves, but in log-polar space. Also notice that these are basic components of shapes with single scale symmetry (an image that is invariant to a single similarity transformation), you can build any sort of scale-symmetric spiral from summing these.

If we take a 2D image I:
and get its Fourier transformation F(I) then all shapes in any location will produce the same abs(F(I)):   (zoom in on centre to see):
Now, if we take the log of this log(abs(F(I))) (as though it were on a complex plane), i.e. convert to log polar coordinates: (horiz = angle, vertical = log radius from bottom to top)
The result is vertically periodic only if it is a broad class of scale-symmetric shape. (Above we see the lines getting softer towards the bottom, that is just because there aren't enough pixels close to the abs(F(I)) centre to represent the low spatial frequencies of the red squares shape. The precise result would be periodic).

The height offset of the vertical translation symmetry represents the scale involved in the fractal. In this case it is scale=2.

Next, we take the Fourier transformation of this, and get its absolute, to give the Fourier Mellin Transformation: FMT = abs(F(log(abs(F(I))))):
This gives an encoding of the shape that is now invariant to rotations and scalings.

Let's do it again, this time with a Menger sponge...
                      I                                      abs(F(I))                       log(abs(F(I)))      
                                               FMT = abs(F(log(abs(F(I)))))
and with the 2D Cantor set:
                                   I                                                 FMT(I)
As you can see, the Fourier-Mellin transformation is different to the squares examples. The intermediate images are hard to see, so you'll need to zoom in or brighten them to see better.

Now, this FMT is invariant to translation, rotation and scale. To illustrate, here is the FMT calculated for the Viscek fractal at two different scales and rotations:

Not only does the FMT represent the shape independent of its translation, rotation and scale, but for a large class of fractals, they can represent their infinitely detailed structure.

This is the class of fractal sets S of the form S=X u Ui(TiS)  where X is a set, U is union and Ti are a set of transformations whose scale are mutually rational, and whose rotations are rational multiples of 2PI. The class includes most simple fractals, including the ones above.

For such fractals, the log(abs(F(I))) image is vertically periodic (and horizontally periodic as it represents angle), this means its abs(FFT()) (the FMT of I) fully captures the entire shape.

However, we still need to quantise this FMT to some finite resolution, but the important thing here is that a finite resolution of the FMT does not inverse transform into a finite resolution image, instead the inverse FMT generates an unbounded fractal, with infinite detail as you zoom in to the centre. Further from the centre gets more blurry (faster if the FMT is lower res), but you can move the fractal so any different parts can go under the detailed centre. So there is access to the whole fractal in high detail, you just have a foviated view of it.

Thursday, February 20, 2020

Laplacian Trees

As a follow on to the previous post about Laplacian growth. The conclusion was really that the majority of techniques tend to the same structure at large scales, namely a void-tree.

There is an interesting extension to bidirectional Laplacian growth, which should generate 'wide' tree structures (a tree-tree to use my terminology).

The idea is to have two sources of chemical (A and B) which both diffuse with equal diffusion rate from their sources (in our example, a line at the top for A and a line at the bottom for B). There is then a single scalar field C that is both grown and shrunk by interacting with these A and B chemicals (in our example the boundary starts midway between the two lines). A and B are both set to 0 at the boundary of C, so the growth is based on the gradient at the boundary dA and dB, using the formula:

growth rate of C = g((dA^a)(dB^b)-s)(dA - dB).

Note that the growth rate of C is signed, so the boundary can grow or reduce.

The g value just dictates the rate of growth. For good results it should be slower than the diffusion of A and B, if it isn't then you can get dynamic behaviour with pulsating growths.

The s value represents smoothness, when it is smaller then thinner trees are accepted.

The a and b values represent the relative important of the two chemicals in building the tree. When a is small the branches filled by chemical A are thinner, and vice versa.

Implementation

As seen in the previous post, there are multiple ways to implement this idea:

1. Using a shader, run the diffusion, and allow a slight gradient to C, with some careful interaction to achieve the formula above. Conceptually simple, but slow (particularly in 3D).

2. Use the fuse-breakdown paper as a model to grow it stochastically on a grid (but without the diffusion component or the random walk).

3. Model the boundary as a polyline with one 1/x drop off per vertex, using the fuse-breakdown model, but not stochastic, growth the whole line.

4. Alternates: we might do a combination of 2,3, or a hybrid of 1,3 like the coral paper, basically a shader for the diffusion and a curve for the growth.

I have used method 1, and preliminary results are as follows:
Chemical A is diffused from the bottom as 5 red sources (to give a variation). Chemical B is the cyan bar on top. a=b=0.5, s=0.0025


If we remove the middle red source then the bounary drops down in the middle, and gets smoother. That is a consequence of the method, the boundary is smoother where the chemical concentrations are less:


If we modify it to a=0.375,b=0.675 then the red-filled areas are slightly thinner, looking more like upwards trees and less like downwards trees:


Here is a more extreme assymetry: a=0.3,b=0.7:

Update:

Using proper gradients (not just neighbour pixel values)
Here are progressively smoother versions from s=0.05 up to s=0.8 (a,b=0.5):

The power of a,b seems to change it from super-fractal (rough at small scales) to sub-fractal (smooth at small scales), from 0.125 to 2:

We can also define a power on the difference, i.e. g((dA^a)(dB^b)-s)(dA^n - dB^n). We might expect larger n to be more directed as it is with DBM, here it is from n=0.5 to n=2:


As you can see, there is not a lot of difference in the results here. This could be that the actual gradients are too approximate to yield a strong difference. But we can see something... the boundary starts as a slightly sine wave, and the n=0.5 has held and exaggerated that broad wave more then n=2. This seems more the opposite of what I'd expect. 

The current algorithm (s=0.2,a,b=0.5,n=1)  does grow towards targets... but seemingly only slightly:


A consequence of the model is that the growth is faster and more erratic closer to the source, and really starts to move like something on fire, as the growth balances with the diffusion of the chemicals (like air and the combustible compound).

Here's another asymmetrical one (a=0.45,b=0.55), notice the branching is more upwards:


It currently suffers from the upwards branches coming and going in repeated manner, and not as much height variation as for equal a,b for some reason.

Update

Here I have made one source the centre and another source the surrounds:
This is a little different, I have made the initial sources self-perpetuate by scaling up the chemicals everywhere by a small amount, to counter the boundary which is a sink for both chemicals. The consequence is that the blob doesn't have to remain around a spherical source. But if the scaling is too great, then large branches can literally bud off into new blobs.

Instead of grown the shape pixel by pixel stochastically, we can grow it using the value of the boundary channel to give sub-pixel accuracy. The result is quite similar:
One difference is that the stochastic one doesn't symmetry-break as easily, so it takes longer for smooth surfaces to become rough. On the other hand, in cases where the result should be smooth, it models the smoothness better. Using stochastic or not is just a boolean flag. Since this version isn't stochastic, it is symmetric if the original circle boundary (and chemical distribution) is exactly symmetric:
The pattern here is just a function of the imperfect disk rendered to pixels. 

Wednesday, January 29, 2020

The multiple guises of Laplacian growth

There is an interesting 3-parameter family of aggregation structures, which seems to encompass a lot of different disparate ideas, so I'll lay out what I've learnt here, to try and get it clear.

The idea of aggregation structures is that a seed shape then grows by accumulating a dissolved resource from the medium that it is in. An example is a sponge or other coral, growing from the accumulation of calcium in the sea water.

More geometrically, we can think of an initial set in 2D or 3D space, which is a sink, and a surrounding distant source continues to replenish the space within it, modelling the effect of a large sea. Point sources or other shaped sources are possible, each producing the resource at various rates.

The resource is usually modelled as diffusing in the medium, the simplest diffusion is Laplacian, where the curvature of resource concentration is everywhere zero. This form of aggregation is therefore generally termed Laplacian growth (or Laplacian instability).

Now we have three parameters describing the process. This is best seen by thinking of the resource as millions of individual grains. Each grain is emitted from the source, and follows a random walk until it touches the aggregation structure, where it can stick and form part of the new structure. I'll look at each parameter in turn:

1. Determinism n

This describes how stochastic or deterministic the system is, with two extremes:
  1. Stochastic: The particles can be released one at a time, and so there is always a larger structure for each subsequent particle. This is the Diffusion Limited Aggregation method:   
    A similar method is called the Dialectric Breakdown Model, and treat the problem as an equivalent electrical one. It constrains the charge on the aggregation structure to be 0 (and positive on the source structures), then solves the Laplacian field for the (conductive) space in between. The probability of the particle attaching to any point on the structure is proportional to the gradient in the Laplacian at that point.
    The results of DBM look very similar to DLA and so they should, as far as I can see this is the same algorithm. The Laplacian field is exactly the probability density of particles in the medium due to random walks. 
  2. Deterministic: The particles are all laid down at once, layer by layer. This means that in any one layer, the aggregation for each particle is independent of the other particles. We can implement it like the DBM: instead of using the gradient of the Laplacian on the border as a probability distribution, you use it as the layer's thickness. This second extreme gives smooth results. The mathematics apparently gets unsolvable in finite time, but this can be fixed by adding a tiny amount of surface tension (next).
Figure: Laplacian growth: source is the two red dots, the start seed shape is the green half-plane.
These extremes also represent the difference between sequential and concurrent aggregation. We can find any process in between these two extreme options by choosing the number of points n released for each layer. With 1 release being option 1 and infinitely many being option 2.

Here is an image of option 2 iterated to high precision, with a disk as the seed and large ring as the source. You might expect the deterministic approach to add a uniform skin to the seed each iteration, and this is true, but any small inconsistencies escalate in what is called Laplacian instability, leading to a pattern like below:

 
Figure: showing the 'Laplacian instability' inherent in Laplacian growth.

Now this is starting to look a bit more like the stochastic option, and indeed I suspect that they tend to the same thing. As each is computed at higher and higher precision (and with smaller and smaller particles in the stochastic case), it seems they should converge to the same structure. Indeed this paper suggests that they are, here is a higher resolution comparison:
Notice that the high resolution DLA (50 million particles) starts to mimic the radially outward branching of the Laplacian growth. 

Given that these two (probably) tend to the same structure, we could replace the parameter n with a particle size parameter. When it is 0, the result is equal to the deterministic case (with no surface tension).

2. Surface tension p   (scale-dependent)

For the deterministic case we can specify a surface tension by setting the charge equal to p times the curvature of the shape boundary. Larger values of p show smoother, fatter fingers in what is now called viscous fingering.
Figure: Viscous fingering by squeezing a viscous fluid between two plates.

We can mimic the same surface tension by taking the charge some small skin width away from the surface, as the probability density in the stochastic case and as the growth speed in the deterministic case.
Figure: increasing surface tension

On a stochastic (grid-based) aggregation, this is mimicked by making the probability of particle attachment proportional to p^(3-N) where N is the number of neighbour cells that are part of the fixed structure. (The 3 changes depending on your neighbourhood size).

Figure: 3D DLA, left: low surface tension (sticks when hits one), right: higher surface tension (images from syntopia's blog).

The surface tension parameter p has units of distance, and fixes a feature scale on the structure, giving non-scale symmetry. However, for larger and larger structures it (should) tend towards the zero surface tension appearance. 

3. Directedness d

Directedness appears in the DBM model by raising the distribution of the Laplacian gradient along the border to some power d. For n=1, the effect of increasing d above 1 is that the dendritic aggregation structure has longer lines between branches, and looks more like lightning, and ultimately cracks and then just a straight line. It takes the fractal dimension down from about 1.72 (for normal DBM/DLA) down to 1 as d approaches infinity. On the other extreme, when d is lower than one, the DBM gets denser, ultimately leading to an Eden cluster when d=0.






Figure: left d=0, middle d=1, right d=2

For the n=infinity (deterministic) case, large values of d tend to pointier and less branching fingers in the structure. As p tends upwards, the result seems to tend towards a greedy shortest path from the seed to the nearest point in the source, then onwards to the next nearest point. Giving a sort of spanning tree.

Figure: left d=1, right d=2 gives thinner branches

Figure: d=2, when the right source point is slightly higher, it grows to the nearest point first.

This d parameter has some sort of meaning under the dialectric breakdown model, but it also has meaning in terms of resource particles. It means that aggregation only occurs when p particles hit the surface point at the same time. This is reminiscent of the Reaction Diffusion Systems (particularly the Gray-Scott model), where powers of the density at each pixel correspond to number of simultaneous particles in the reaction. Here it means that you need d particles and one border point in the same pixel to create a new piece of the structure.