Reaction-Diffusion

using MagicBook.Algorithms.Randomness;

// Alan Turing asked in 1952 how a uniform blob of cells can end up as a leopard, and
// answered it with two substances that spread and react with each other. The model
// used here is Gray-Scott: A is fed in from outside, B eats A and turns into more B,
// and B is drained off. Both spread, but A twice as fast as B.
//
// Everything depends on two numbers. feed is how much A is fed in, kill how much B is
// drained. A hundredth in either direction turns stripes into dots or lets the whole
// pattern die.
public static class ReactionDiffusion
{
    private const float SpreadA = 1.0f;
    private const float SpreadB = 0.5f;

    public static float[] Simulate(int width, int height, int steps, float feed, float kill, ulong seed)
    {
        var a = new float[width * height];
        var b = new float[width * height];

        Array.Fill(a, 1f);
        Splash(b, width, height, seed);

        var nextA = new float[a.Length];
        var nextB = new float[b.Length];

        for (var step = 0; step < steps; step++)
        {
            for (var cell = 0; cell < a.Length; cell++)
            {
                // How much A and B meet here. Two B are needed to turn one A into
                // another B, which is why the square is in there.
                var reaction = a[cell] * b[cell] * b[cell];

                nextA[cell] = a[cell] + SpreadA * Spread(a, width, height, cell) - reaction + feed * (1 - a[cell]);
                nextB[cell] = b[cell] + SpreadB * Spread(b, width, height, cell) + reaction - (kill + feed) * b[cell];

                nextA[cell] = Math.Clamp(nextA[cell], 0f, 1f);
                nextB[cell] = Math.Clamp(nextB[cell], 0f, 1f);
            }

            (a, nextA) = (nextA, a);
            (b, nextB) = (nextB, b);
        }

        return b;
    }

    // How much flows in from the neighbours: the average around the cell minus the
    // cell itself. Diagonal neighbours count less, because they are further away.
    // This is the Laplacian, the same operator the edge filters use.
    private static float Spread(float[] values, int width, int height, int cell)
    {
        var x = cell % width;
        var y = cell / width;
        var sum = -values[cell];

        for (var dy = -1; dy <= 1; dy++)
        {
            for (var dx = -1; dx <= 1; dx++)
            {
                if (dx == 0 && dy == 0)
                    continue;

                // The edges are glued together, so the pattern has no border.
                var neighbour = ((y + dy + height) % height) * width + (x + dx + width) % width;

                sum += values[neighbour] * (dx == 0 || dy == 0 ? 0.2f : 0.05f);
            }
        }

        return sum;
    }

    // Without a disturbance nothing would ever happen: A is everywhere, B nowhere,
    // and that state stays as it is. So a few drops of B are put in.
    private static void Splash(float[] b, int width, int height, ulong seed)
    {
        var random = new Pcg(seed);

        for (var drop = 0; drop < 14; drop++)
        {
            var centreX = random.Next(width);
            var centreY = random.Next(height);

            for (var y = -4; y <= 4; y++)
            {
                for (var x = -4; x <= 4; x++)
                    b[((centreY + y + height) % height) * width + (centreX + x + width) % width] = 1f;
            }
        }
    }
}

const int width = 200;
const int height = 150;
const int steps = 6000;

// The same two rules, the same starting drops, only feed and kill differ -
// and out comes something completely different every time.
(string Name, float Feed, float Kill)[] recipes =
[
    ("coral", 0.0545f, 0.0620f),
    ("stripes", 0.0580f, 0.0650f),
    ("spots", 0.0300f, 0.0620f),
];

Console.WriteLine("pattern   feed    kill    covered by B");

foreach (var (name, feed, kill) in recipes)
{
    var b = ReactionDiffusion.Simulate(width, height, steps, feed, kill, seed: 20240921);

    // B never gets anywhere near 1, so for the picture the values are pulled
    // apart until the strongest spot is white.
    var strongest = b.Max();
    var shown = new float[b.Length];

    for (var i = 0; i < b.Length; i++)
        shown[i] = b[i] / strongest;

    Channels.ToImage(shown, width, height).Save(Path.Combine(folder, $"reaction-{name}.png"));

    Console.WriteLine($"{name,-8}  {feed:0.0000}  {kill:0.0000}  " +
        $"{b.Count(value => value > 0.25f) * 100.0 / b.Length,10:0.0} %");
}