Reaction-Diffusion

using MagicBook.Algorithms.Randomness;

// Alan Turing hat 1952 gefragt, wie aus einem gleichförmigen Klumpen Zellen ein
// Leopard wird, und mit zwei Stoffen geantwortet, die sich ausbreiten und miteinander
// reagieren. Hier steht Gray-Scott: A kommt von außen nach, B frisst A und wird dabei
// zu mehr B, und B wird abgeführt. Beide streuen, A aber doppelt so schnell wie B.
//
// Alles hängt an zwei Zahlen. feed ist, wie viel A nachkommt, kill, wie viel B
// abgeführt wird. Ein Hundertstel in eine der beiden Richtungen macht aus Streifen
// Punkte oder lässt das ganze Muster absterben.
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++)
            {
                // Wie viel A und B hier zusammentreffen. Es braucht zwei B, um ein A in ein
                // weiteres B zu verwandeln, daher das Quadrat.
                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;
    }

    // Wie viel von den Nachbarn hereinfließt: der Durchschnitt ringsum minus dem
    // eigenen Wert. Diagonale Nachbarn zählen weniger, weil sie weiter weg sind.
    // Das ist der Laplace-Operator, derselbe, den die Kantenfilter verwenden.
    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;

                // Die Ränder sind zusammengeklebt, deshalb hat das Muster keinen Rand.
                var neighbour = ((y + dy + height) % height) * width + (x + dx + width) % width;

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

        return sum;
    }

    // Ohne Störung passiert nie etwas: A überall, B nirgends, und dieser Zustand
    // bleibt, wie er ist. Also werden ein paar Tropfen B hineingesetzt.
    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;

// Dieselben zwei Regeln, dieselben Starttropfen, nur feed und kill sind anders -
// und heraus kommt jedes Mal etwas völlig anderes.
(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 kommt nie in die Nähe von 1, deshalb werden die Werte für das Bild
    // auseinandergezogen, bis die stärkste Stelle weiß ist.
    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} %");
}