Weighted Voronoi Stippling

using MagicBook.Algorithms.Randomness;

// Ein Punkt des fertigen Bildes. Der Radius zählt in Pixeln der Vorlage.
public readonly record struct Stipple(float X, float Y, float Radius);

public static class WeightedVoronoiStippling
{
    // Macht aus einem Bild Punkte: dunkle Bereiche bekommen viele, helle wenige.
    // brightness enthält width * height Werte von 0 (schwarz) bis 1 (weiß).
    public static Stipple[] Distribute(float[] brightness, int width, int height, int dotCount, int rounds, ulong seed)
    {
        // Tinte ist der Stoff, aus dem die Punkte sind, also genau andersherum als Helligkeit.
        var ink = new float[brightness.Length];
        for (var i = 0; i < ink.Length; i++)
            ink[i] = 1f - brightness[i];

        if (ink.Sum() <= 0)
            throw new ArgumentException("The picture is completely white, there is nothing to draw.", nameof(brightness));

        var points = Scatter(ink, width, height, dotCount, seed);

        // Lloyds Algorithmus: Jede Runde schiebt jeden Punkt in den Schwerpunkt der
        // Fläche, die ihm näher liegt als jedem anderen Punkt. Weil der Schwerpunkt mit
        // Tinte gewichtet wird, wandern die Punkte ins Dunkle und verteilen sich zugleich
        // gleichmäßig. Eine Handvoll Runden genügt.
        for (var round = 0; round < rounds; round++)
            Relax(points, ink, width, height);

        return Measure(points, ink, width, height);
    }

    // Wirft Punkte ins Bild und behält die, die auf genug Tinte landen.
    // Dieser Start stimmt schon ungefähr und spart dadurch später ein paar Runden.
    private static (float X, float Y)[] Scatter(float[] ink, int width, int height, int count, ulong seed)
    {
        var random = new Pcg(seed);
        var points = new (float X, float Y)[count];

        for (var i = 0; i < count; i++)
        {
            while (true)
            {
                var x = random.Next(width);
                var y = random.Next(height);

                // Je heller die Stelle, desto wahrscheinlicher wird der Wurf verworfen.
                if (random.NextDouble() > ink[y * width + x])
                    continue;

                points[i] = (x, y);
                break;
            }
        }

        return points;
    }

    // Eine Runde von Lloyds Algorithmus.
    private static void Relax((float X, float Y)[] points, float[] ink, int width, int height)
    {
        var sumX = new float[points.Length];
        var sumY = new float[points.Length];
        var sumInk = new float[points.Length];

        for (var y = 0; y < height; y++)
        {
            for (var x = 0; x < width; x++)
            {
                var weight = ink[y * width + x];

                // Weiße Pixel tragen keine Tinte und ziehen deshalb auch nichts.
                if (weight <= 0)
                    continue;

                var nearest = Nearest(points, x, y);

                sumX[nearest] += x * weight;
                sumY[nearest] += y * weight;
                sumInk[nearest] += weight;
            }
        }

        for (var i = 0; i < points.Length; i++)
        {
            // Ein Punkt, der gar keine Tinte erwischt hat, bleibt einfach liegen.
            if (sumInk[i] > 0)
                points[i] = (sumX[i] / sumInk[i], sumY[i] / sumInk[i]);
        }
    }

    // Stumpfe Suche über alle Punkte. Quadrierte Abstände genügen hier, weil die
    // Wurzel nichts daran ändert, welcher davon der kleinste ist.
    private static int Nearest((float X, float Y)[] points, int x, int y)
    {
        var best = 0;
        var bestDistance = float.MaxValue;

        for (var i = 0; i < points.Length; i++)
        {
            var dx = points[i].X - x;
            var dy = points[i].Y - y;
            var distance = dx * dx + dy * dy;

            if (distance >= bestDistance)
                continue;

            bestDistance = distance;
            best = i;
        }

        return best;
    }

    // Zum Schluss bekommt jeder Punkt eine Größe. Die Regel ist einfach und sorgt
    // dafür, dass das Ergebnis stimmt: Ein Punkt bedeckt genau so viel Fläche, wie in
    // seiner Zelle Tinte steckt. Alle Punkte zusammen tragen also die Tinte des Bildes.
    private static Stipple[] Measure((float X, float Y)[] points, float[] ink, int width, int height)
    {
        var sumInk = new float[points.Length];

        for (var y = 0; y < height; y++)
        {
            for (var x = 0; x < width; x++)
                sumInk[Nearest(points, x, y)] += ink[y * width + x];
        }

        var stipples = new Stipple[points.Length];

        for (var i = 0; i < points.Length; i++)
        {
            // Die Kreisfläche ist Pi mal Radius im Quadrat, das hier ist der Weg zurück.
            var radius = MathF.Sqrt(sumInk[i] / MathF.PI);

            stipples[i] = new Stipple(points[i].X, points[i].Y, radius);
        }

        return stipples;
    }
}

// file zeigt auf ein Graustufen-Portrait, 260 x 300 Pixel.
var image = PngImage.Load(file);

var dots = WeightedVoronoiStippling.Distribute(
    image.ToBrightness(), image.Width, image.Height, dotCount: 6000, rounds: 12, seed: 4711);

// Die Punkte doppelt so groß zeichnen wie die Vorlage, damit sie rund bleiben.
StippleRenderer.Render(dots, image.Width, image.Height, scale: 2).Save(result);

Console.WriteLine($"{dots.Length} dots, average radius {dots.Average(dot => dot.Radius):0.00} pixels");
Console.WriteLine($"written to {Path.GetFileName(result)}");