Dithering

// A picture that may only use a few colours falls apart into patches: everything that
// lies close together lands on the same colour, and soft transitions turn into hard
// steps. Dithering trades those patches for a fine grain that the eye mixes back
// together into the colour that was meant.
public static class Dithering
{
    // Floyd and Steinberg, 1976. Every pixel is set to the nearest colour that is
    // still allowed, and the difference between the wanted and the used colour is
    // handed on to the neighbours that are still to come.
    public static PngImage FloydSteinberg(PngImage image, int levels)
    {
        var channels = ToChannels(image, levels);

        for (var y = 0; y < image.Height; y++)
        {
            for (var x = 0; x < image.Width; x++)
            {
                for (var channel = 0; channel < 3; channel++)
                {
                    var index = ((y * image.Width + x) * 3) + channel;
                    var wanted = channels[index];
                    var used = Snap(wanted, levels);

                    channels[index] = used;

                    // The error is split over four neighbours. Most of it goes to the
                    // right, the rest into the line below; together the parts add up
                    // to the whole error, so nothing is lost and nothing invented.
                    var error = wanted - used;

                    Hand(channels, image, x + 1, y, channel, error * 7f / 16f);
                    Hand(channels, image, x - 1, y + 1, channel, error * 3f / 16f);
                    Hand(channels, image, x, y + 1, channel, error * 5f / 16f);
                    Hand(channels, image, x + 1, y + 1, channel, error * 1f / 16f);
                }
            }
        }

        return ToImage(channels, image);
    }

    // The other school: instead of passing the error on, every pixel is nudged up or
    // down by a fixed amount before it is snapped. The amount comes from a small
    // matrix that repeats over the picture, which is why the result looks woven.
    public static PngImage Ordered(PngImage image, int levels)
    {
        var channels = ToChannels(image, levels);
        var step = 1f / (levels - 1);

        for (var y = 0; y < image.Height; y++)
        {
            for (var x = 0; x < image.Width; x++)
            {
                // The matrix value runs from 0 to 15, the nudge from minus half a step
                // to plus half a step.
                var nudge = ((Matrix[(y % 4) * 4 + (x % 4)] + 0.5f) / 16f - 0.5f) * step;

                for (var channel = 0; channel < 3; channel++)
                {
                    var index = ((y * image.Width + x) * 3) + channel;
                    channels[index] = Snap(channels[index] + nudge, levels);
                }
            }
        }

        return ToImage(channels, image);
    }

    // The same reduction without any dithering: every pixel simply takes the nearest
    // allowed colour. This is what the two above are compared against.
    public static PngImage Quantize(PngImage image, int levels)
    {
        var channels = ToChannels(image, levels);

        for (var index = 0; index < channels.Length; index++)
            channels[index] = Snap(channels[index], levels);

        return ToImage(channels, image);
    }

    // Bayer's 4 x 4 matrix. Neighbouring numbers sit as far apart as possible, so the
    // pattern stays even and does not clump together anywhere.
    private static readonly int[] Matrix =
    [
        0, 8, 2, 10,
        12, 4, 14, 6,
        3, 11, 1, 9,
        15, 7, 13, 5,
    ];

    // The allowed values of a channel lie evenly spread between 0 and 1: with three
    // levels those are 0, 0.5 and 1. Snap picks the nearest of them.
    private static float Snap(float value, int levels)
    {
        var steps = levels - 1;

        return Math.Clamp(MathF.Round(value * steps) / steps, 0f, 1f);
    }

    // Adds a part of the error to a neighbour, as long as that neighbour exists.
    private static void Hand(float[] channels, PngImage image, int x, int y, int channel, float part)
    {
        if (x < 0 || y < 0 || x >= image.Width || y >= image.Height)
            return;

        channels[((y * image.Width + x) * 3) + channel] += part;
    }

    // Red, green and blue as values from 0 to 1. The error diffusion has to work on
    // fractions, a byte per channel would round it away immediately.
    private static float[] ToChannels(PngImage image, int levels)
    {
        ArgumentOutOfRangeException.ThrowIfLessThan(levels, 2);

        var channels = new float[image.Width * image.Height * 3];

        for (var pixel = 0; pixel < image.Width * image.Height; pixel++)
        {
            for (var channel = 0; channel < 3; channel++)
                channels[pixel * 3 + channel] = image.Pixels[pixel * 4 + channel] / 255f;
        }

        return channels;
    }

    private static PngImage ToImage(float[] channels, PngImage original)
    {
        var pixels = new byte[original.Width * original.Height * 4];

        for (var pixel = 0; pixel < original.Width * original.Height; pixel++)
        {
            for (var channel = 0; channel < 3; channel++)
                pixels[pixel * 4 + channel] = (byte)MathF.Round(channels[pixel * 3 + channel] * 255f);

            pixels[pixel * 4 + 3] = original.Pixels[pixel * 4 + 3];
        }

        return new PngImage { Width = original.Width, Height = original.Height, Pixels = pixels };
    }
}

// file points to a colour photo, 320 x 240 pixels.
var photo = PngImage.Load(file);

// Levels is how many values a channel may take. Two levels leave eight
// colours in total, three levels leave 27, four leave 64.
foreach (var levels in new[] { 2, 3, 4 })
{
    var dithered = Dithering.FloydSteinberg(photo, levels);
    var target = Path.Combine(folder, $"grand-prismatic-spring-{levels}.png");

    dithered.Save(target);

    Console.WriteLine(
        $"{levels} levels: {levels * levels * levels,3} colours allowed, " +
        $"{dithered.ColourCount(),3} used, {new FileInfo(target).Length / 1024,3} KiB");
}

var photo = PngImage.Load(file);

// Two ways of measuring the mistake: pixel by pixel, and over squares of
// 4 x 4 pixels. The second one is closer to what the eye does, because it
// averages over a small area before it judges.
Console.WriteLine("levels  colours  method      pixel   4 x 4");

foreach (var levels in new[] { 2, 3, 4 })
{
    (string Name, PngImage Result)[] ways =
    [
        ("plain", Dithering.Quantize(photo, levels)),
        ("ordered", Dithering.Ordered(photo, levels)),
        ("diffused", Dithering.FloydSteinberg(photo, levels)),
    ];

    foreach (var (name, result) in ways)
    {
        Console.WriteLine(
            $"{levels,6}  {levels * levels * levels,7}  {name,-8}  " +
            $"{Error(photo, result, 1),6:0.0}  {Error(photo, result, 4),6:0.0}");
    }
}

// The average difference to the original, measured over squares of the given
// size. Size 1 compares single pixels, size 4 compares small areas.
static double Error(PngImage original, PngImage result, int size)
{
    var sum = 0.0;
    var squares = 0;

    for (var y = 0; y + size <= original.Height; y += size)
    {
        for (var x = 0; x + size <= original.Width; x += size)
        {
            for (var channel = 0; channel < 3; channel++)
            {
                var difference = 0;

                for (var offsetY = 0; offsetY < size; offsetY++)
                {
                    for (var offsetX = 0; offsetX < size; offsetX++)
                    {
                        var pixel = (((y + offsetY) * original.Width) + x + offsetX) * 4 + channel;
                        difference += original.Pixels[pixel] - result.Pixels[pixel];
                    }
                }

                sum += Math.Abs(difference) / (double)(size * size);
                squares++;
            }
        }
    }

    return sum / squares;
}