Weighted Voronoi Stippling

using MagicBook.Algorithms.Randomness;

// One dot of the finished picture. The radius is measured in pixels of the source.
public readonly record struct Stipple(float X, float Y, float Radius);

public static class WeightedVoronoiStippling
{
    // Turns a picture into dots: dark areas get many, bright areas few.
    // brightness holds width * height values from 0 (black) to 1 (white).
    public static Stipple[] Distribute(float[] brightness, int width, int height, int dotCount, int rounds, ulong seed)
    {
        // Ink is what the dots are made of, so it is the other way round than brightness.
        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);

        // Lloyd's algorithm: every round moves each dot into the centre of gravity of
        // the area that is closer to it than to any other dot. Because the centre is
        // weighted by ink, the dots drift towards the dark parts and spread out evenly
        // at the same time. A handful of rounds is enough.
        for (var round = 0; round < rounds; round++)
            Relax(points, ink, width, height);

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

    // Throws dots at the picture and keeps the ones that land on enough ink.
    // This start is already roughly right, which saves a couple of rounds later.
    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);

                // The brighter the spot, the more likely the throw is rejected.
                if (random.NextDouble() > ink[y * width + x])
                    continue;

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

        return points;
    }

    // One round of Lloyd's algorithm.
    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];

                // White pixels carry no ink and therefore pull nothing.
                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++)
        {
            // A dot that caught no ink at all simply stays where it is.
            if (sumInk[i] > 0)
                points[i] = (sumX[i] / sumInk[i], sumY[i] / sumInk[i]);
        }
    }

    // Plain search over all dots. Squared distances are enough here, because the
    // square root would not change which one comes out smallest.
    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;
    }

    // Finally each dot gets a size. The rule is simple and it is what makes the
    // result look right: a dot covers exactly as much area as there is ink in its
    // cell. All dots together therefore carry exactly the ink of the picture.
    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++)
        {
            // Area of a circle is pi times radius squared, so this is the way back.
            var radius = MathF.Sqrt(sumInk[i] / MathF.PI);

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

        return stipples;
    }
}

// file points to a grayscale portrait, 260 x 300 pixels.
var image = PngImage.Load(file);

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

// Draw the dots twice as large as the original, so they stay round.
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)}");