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.


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:


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.


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.

Thursday, December 19, 2019

2D Eddie Currents

There was a nice paper I read about exact solutions of the Navier Stokes equations in 2D for the incompressible fluid case, and for the case of zero viscosity (giving the simpler Euler equations).

The trick is that any 'stream function' Ψ(x,y) such that its mean curvature d^2Ψ(x,y)/dx^2 + d^2Ψ(x,y)/dy^2 = f(Ψ(x,y)) for some simple function f().

When we have f(Ψ) = K*Ψ then we have a nice repeating general solution:

Ψ(x,y) = sum_i k_i*sin(a_i*x + b_i*y + c_i)  such that a_i^2+b_i^2=d for all i.

giving K=-d^2.
This stream function is the sum of 2D planar waves of equal wavelength. It seems capable of approximating a wide range of functions.

The vector field is then velocities: u=dΨ/dy, v=-dΨ/dx, which is:

u = sum_i k_i*b_i*cos(a_i*x + b_i*y + c_i)
v =-sum_i k_i*a_i*cos(a_i*x + b_i*y + c_i)

The waves can then all be moved at the same rate, by replacing c_i with c_i+wt, where t is time. One could also vary the k_i parameters, but I wouldn't vary a_i or b_i over time, as it will cause very fast variation far from the origin.

This is an example vector field for:

u = 2*4*sin(3*0.1*x + 4*0.1*y + a) + 2*5*d*sin(5*0.1*y + b) + 2*3*sin(4*0.1*x + 3*0.1*y + c)
v =-2*3*sin(3*0.1*x + 4*0.1*y + a)                          - 2*4*sin(4*0.1*x + 3*0.1*y + c)

where d = 1.5.

Whenever a,b,d form a pythagorean triple, you can use combinations of these to create vector fields that tile precisely, even when the sinusoids are moving using t.

These appear to represent perfectly valid (but not necessarily stable) persistent vector fields that can change over time. Making a great tool for eddy currents around a boat, or in the atmosphere or just for wind currents.

The flexibility in k,a,b and c may allow the local neighbourhood of the currents to be somewhat interactive. For example, if you have a boat spewing out water as it drives along, maybe we can find the minimal adjustment to the set of sinusoids that matches the desired velocity at that point, or at multiple points.

If you allowed differing wavelengths then the vector field would still be divergence-free, but it would not have a valid pressure field.

Pairs of planar waves always produce rectangular 2D vortex structures:
Triplets of planar waves generally produce irregular wave structures, but you can get a nice triangular grid of vortices like so:
(Wolfram: plot sin(x) + sin(x/2 + y*sqrt(3)/2) + sin(-x/2 + y*sqrt(3)/2))

Local Eddies

Imagine you averaged all the planar cosine waves cos(ax+by) where a^2+b^2=1. The result is a radial Bessel function J0(sqrt(x^2+y^2)). Since it is a sum of sinusoids of equal wavelength, it also is a valid stream function. 

Now, if we sum together these functions (which decay with distance from centre), then our resulting velocity field decays to nothing far from the area in question. This allows local control of the vortices (and even directly specifies a main vortex), and you can add it onto a background vortex field (of planar sinusoids), or just to a constant velocity field. 

Two equal and opposite nearby Besel functions, creates a gush of water through the middle, and two vortices on either side.

There is also a dual function, which represents the linear rather than circular flow, and is also local. It is the integral of cos(r cos(angle))*cos(2*angle + k) with respect to all polar angles. Each k gives a different radial function, but they basically just blend from one Bessel function (J2) to its negative, based on the polar angle:

This gives you a local linear flow 'cross':
(Wolfram: plot cos(2*atan(x/y))*J2(sqrt(x^2+y^2))).

since eddies are made out of combinations of vortices and linear flows, it might work OK to use these two Bessel functions as the atomic generators of local eddy flows.

Tuesday, October 16, 2018

Nested Spheres

This is a new shape that I have made based on sphere inversions. It is a nested set of spheres, such that each level of spheres is a Ford-Farey sphere packing over the Riemann sphere.

The structure is a tree, it has no self contact, but is also probably the densest packing that such a tree can be.
Each level is a graphical depiction of the Eistenstein rationals. So the total shape depicts an infinite cascade of Eisenstein rationals, rather like each nested sphere depicts the rationals over the infinitesimal gap around its parent's location on the complex plane. 

As such, one might describe it as a cascade-like depiction of the complex numbers. 

Saturday, May 12, 2018

Dendrite and cross fractal families

Here are two fractal families that are probably not well known as families despite having well known members. 

The first I'll call the dendrite family because it contains the attractive Pentadendrite. The six sided version is the Hexaflake and the four sided version is the Vicsek fractal and the two sided version is just a line:

The Hausdorff dimension is log(n+1)/log(3) where n are the number of sides, for the n=even case.
The n=odd case is a little higher. The three sided fractal I haven't seen before, so I'll call it the tridendrite. here's a close up:

The fractal is built up from a three pointed star by swinging a copy of itself 180 degrees around each outer point:
For even pointed dendrites the new outer point is simply the new farthest, for odd pointed dendrites like this, there are two equidistant new outer points, we choose either the clockwise, or the anticlockwise and stick with it, repeating the process ad infinitum. 

Consequently there are two resulting fractals, one being a mirror image of the other. So while all the dendrites have n-fold rotational symmetry, the even pointed dendrites also have mirror symmetry. 

A variation is to replace the shape each iteration with n copies, rotated 180 degrees and attached by their outer points to (0,0):
and repeated, giving a denser "cross" family:

The n=3 case is sometimes called a Fudgeflake, or can be built from three terdragons. The four pointed case is just a 90 degree Koch curve or Greek cross fractal. The larger n cases include substantial self-overlap.

By the way, here is a close up of the pentacross, the density variation has some nice patterns: