Kantenerkennung

// Eine Kante ist eine Stelle, an der sich die Helligkeit schnell ändert. Kanten zu
// suchen heißt also, die größte Steigung zu suchen - die Ableitung des Bildes. Weil
// ein Bild keine glatte Funktion ist, sondern ein Raster aus Zahlen, wird sie mit
// einer kleinen Matrix geschätzt, und die Wahl dieser Matrix ist alles, was die
// klassischen Filter unterscheidet.
public static class EdgeDetection
{
    // Sobel gewichtet die mittlere Zeile dreifach, was das Ergebnis quer zur Kante
    // glättet und den Filter deutlich unempfindlicher gegen Rauschen macht.
    public static readonly float[] SobelX = [-1, 0, 1, -2, 0, 2, -1, 0, 1];
    public static readonly float[] SobelY = [-1, -2, -1, 0, 0, 0, 1, 2, 1];

    // Prewitt ist dieselbe Idee ohne Gewichtung: drei schlichte Differenzen
    // nebeneinander. Etwas schärfer, etwas unruhiger.
    public static readonly float[] PrewittX = [-1, 0, 1, -1, 0, 1, -1, 0, 1];
    public static readonly float[] PrewittY = [-1, -1, -1, 0, 0, 0, 1, 1, 1];

    // Der Laplace-Operator fragt etwas anderes: nicht wie steil es ist, sondern wo die
    // Steigung kippt. Er braucht keine Richtung, eine Matrix genügt - und weil er die
    // zweite Ableitung ist, reagiert er doppelt so stark auf Rauschen.
    public static readonly float[] Laplace = [0, 1, 0, 1, -4, 1, 0, 1, 0];

    // Wie stark sich die Helligkeit an jedem Pixel ändert, und in welche Richtung.
    // Die beiden Matrizen messen die Steigung quer und längs; zusammen bilden sie
    // einen kleinen Pfeil je Pixel, und dessen Länge ist die Stärke.
    public static (float[] Strength, float[] Direction) Gradient(
        float[] gray, int width, int height, float[] horizontal, float[] vertical)
    {
        var alongX = Convolution.Apply(gray, width, height, horizontal);
        var alongY = Convolution.Apply(gray, width, height, vertical);

        var strength = new float[gray.Length];
        var direction = new float[gray.Length];

        for (var i = 0; i < gray.Length; i++)
        {
            strength[i] = MathF.Sqrt(alongX[i] * alongX[i] + alongY[i] * alongY[i]);
            direction[i] = MathF.Atan2(alongY[i], alongX[i]);
        }

        return (strength, direction);
    }

    public static float[] Sobel(float[] gray, int width, int height) =>
        Gradient(gray, width, height, SobelX, SobelY).Strength;

    public static float[] Prewitt(float[] gray, int width, int height) =>
        Gradient(gray, width, height, PrewittX, PrewittY).Strength;

    // Der Laplace-Operator antwortet mit einem Vorzeichen, es zählt also, wie weit die
    // Antwort von null entfernt ist.
    public static float[] Laplacian(float[] gray, int width, int height)
    {
        var answer = Convolution.Apply(gray, width, height, Laplace);
        var strength = new float[answer.Length];

        for (var i = 0; i < answer.Length; i++)
            strength[i] = MathF.Abs(answer[i]);

        return strength;
    }
}

// Fast alles, was lokal mit einem Bild geschieht, funktioniert gleich: eine kleine
// Matrix über jeden Pixel legen, mit dem multiplizieren, was darunter liegt, und
// aufsummieren. Weichzeichnen, Schärfen, Prägen und jeder Kantenfilter sind
// dieselbe Rechnung mit einer anderen Matrix.
public static class Convolution
{
    // Wendet eine quadratische Matrix auf ein Graubild an. Am Rand wird der äußerste
    // Pixel wiederholt, damit es keinen Sonderfall gibt.
    public static float[] Apply(float[] values, int width, int height, float[] kernel)
    {
        var size = (int)Math.Sqrt(kernel.Length);

        if (size * size != kernel.Length || size % 2 == 0)
            throw new ArgumentException("A kernel has to be square and of an odd size.", nameof(kernel));

        var reach = size / 2;
        var result = new float[values.Length];

        for (var y = 0; y < height; y++)
        {
            for (var x = 0; x < width; x++)
            {
                var sum = 0f;

                for (var ky = 0; ky < size; ky++)
                {
                    for (var kx = 0; kx < size; kx++)
                    {
                        var sampleX = Math.Clamp(x + kx - reach, 0, width - 1);
                        var sampleY = Math.Clamp(y + ky - reach, 0, height - 1);

                        sum += values[sampleY * width + sampleX] * kernel[ky * size + kx];
                    }
                }

                result[y * width + x] = sum;
            }
        }

        return result;
    }

    // Ein Gauß-Weichzeichner, und ein Trick dabei: Statt eine n-mal-n-Matrix über das
    // Bild zu legen, wird dieselbe Glockenkurve einmal waagerecht und einmal senkrecht
    // angewandt. Das Ergebnis ist dasselbe, der Aufwand fällt von n² auf 2n.
    public static float[] Blur(float[] values, int width, int height, float radius)
    {
        if (radius <= 0)
            return (float[])values.Clone();

        var bell = Bell(radius);

        return Sweep(Sweep(values, width, height, bell, horizontal: true), width, height, bell, horizontal: false);
    }

    // Die Glockenkurve, an ganzen Pixeln abgetastet und so skaliert, dass ihre Werte
    // sich zu eins addieren. Sonst würde das Bild heller oder dunkler.
    public static float[] Bell(float radius)
    {
        var reach = Math.Max(1, (int)MathF.Ceiling(radius * 3));
        var weights = new float[reach * 2 + 1];
        var sum = 0f;

        for (var i = -reach; i <= reach; i++)
        {
            var weight = MathF.Exp(-(i * i) / (2 * radius * radius));

            weights[i + reach] = weight;
            sum += weight;
        }

        for (var i = 0; i < weights.Length; i++)
            weights[i] /= sum;

        return weights;
    }

    private static float[] Sweep(float[] values, int width, int height, float[] weights, bool horizontal)
    {
        var reach = weights.Length / 2;
        var result = new float[values.Length];

        for (var y = 0; y < height; y++)
        {
            for (var x = 0; x < width; x++)
            {
                var sum = 0f;

                for (var i = 0; i < weights.Length; i++)
                {
                    var sampleX = horizontal ? Math.Clamp(x + i - reach, 0, width - 1) : x;
                    var sampleY = horizontal ? y : Math.Clamp(y + i - reach, 0, height - 1);

                    sum += values[sampleY * width + sampleX] * weights[i];
                }

                result[y * width + x] = sum;
            }
        }

        return result;
    }
}

var photo = PngImage.Load(file);

// Die Filter schauen nur auf die Helligkeit. Farbe sagt nichts über Kanten;
// zwei völlig verschiedene Farben können exakt gleich hell sein.
var gray = photo.ToBrightness();

(string Name, float[] Strength)[] filters =
[
    ("sobel", EdgeDetection.Sobel(gray, photo.Width, photo.Height)),
    ("prewitt", EdgeDetection.Prewitt(gray, photo.Width, photo.Height)),
    ("laplace", EdgeDetection.Laplacian(gray, photo.Width, photo.Height)),
];

Console.WriteLine("filter    strongest  clearly an edge");

foreach (var (name, strength) in filters)
{
    // Die Antworten werden umgedreht: Eine starke Kante soll dunkel auf Weiß sein,
    // so wie eine Bleistiftzeichnung aussieht.
    var drawn = new float[strength.Length];

    for (var i = 0; i < strength.Length; i++)
        drawn[i] = 1 - Math.Min(strength[i], 1f);

    Channels.ToImage(drawn, photo.Width, photo.Height)
        .Save(Path.Combine(folder, $"spring-{name}.png"));

    var strongest = strength.Max();
    var clear = strength.Count(value => value > 0.25f) * 100.0 / strength.Length;

    Console.WriteLine($"{name,-9} {strongest,9:0.00}  {clear,14:0.0} %");
}

// John Canny hat 1986 gefragt, was ein guter Kantenfinder überhaupt sein soll, drei
// Forderungen aufgeschrieben - jede Kante finden, sie richtig verorten, nur einmal -
// und daraus ein Verfahren abgeleitet. Es ist keine Matrix, sondern eine Kette aus
// vier Schritten, und es ist bis heute der Standard.
public static class Canny
{
    public static float[] Detect(float[] gray, int width, int height, float blur, float low, float high)
    {
        // 1. Weichzeichnen. Der Gradient reagiert auf jedes Korn, also wird das Bild
        // zuerst beruhigt. Wie stark, ist der eine echte Regler des Verfahrens.
        var calm = Convolution.Blur(gray, width, height, blur);

        // 2. Gradient, genau wie bei Sobel: wie steil, und in welche Richtung.
        var (strength, direction) = EdgeDetection.Gradient(
            calm, width, height, EdgeDetection.SobelX, EdgeDetection.SobelY);

        // 3. Grate ausdünnen. Ein Anstieg ist mehrere Pixel breit, eine Kante aber
        // eine Linie. Nur der Pixel oben auf dem Grat überlebt.
        var thin = Thin(strength, direction, width, height);

        // 4. Zwei Schwellen. Was über der oberen liegt, ist sicher eine Kante.
        // Was dazwischen liegt, zählt nur, wenn es an so einer sicheren Kante
        // hängt - das hält eine schwächer werdende Linie ganz, statt sie zu punkten.
        return Follow(thin, width, height, low, high);
    }

    // Non-Maximum Suppression: Ein Pixel bleibt nur, wenn er mindestens so stark ist wie
    // seine beiden Nachbarn quer zur Kante. Die Richtung wird auf eine von vieren
    // gerundet, denn andere Nachbarn gibt es auf einem Raster nicht.
    private static float[] Thin(float[] strength, float[] direction, int width, int height)
    {
        var thin = new float[strength.Length];

        for (var y = 1; y < height - 1; y++)
        {
            for (var x = 1; x < width - 1; x++)
            {
                var here = y * width + x;
                var angle = direction[here] * 180 / MathF.PI;

                if (angle < 0)
                    angle += 180;

                var (dx, dy) = angle switch
                {
                    < 22.5f or >= 157.5f => (1, 0),
                    < 67.5f => (1, 1),
                    < 112.5f => (0, 1),
                    _ => (-1, 1),
                };

                var before = strength[(y - dy) * width + x - dx];
                var after = strength[(y + dy) * width + x + dx];

                if (strength[here] >= before && strength[here] >= after)
                    thin[here] = strength[here];
            }
        }

        return thin;
    }

    // Hysterese: bei den sicheren Kanten anfangen und allem folgen, was wenigstens
    // halbwegs stark ist und sie berührt.
    private static float[] Follow(float[] thin, int width, int height, float low, float high)
    {
        var edges = new float[thin.Length];
        var pending = new Stack<int>();

        for (var i = 0; i < thin.Length; i++)
        {
            if (thin[i] >= high)
            {
                edges[i] = 1;
                pending.Push(i);
            }
        }

        while (pending.Count > 0)
        {
            var pixel = pending.Pop();
            var x = pixel % width;
            var y = pixel / width;

            for (var dy = -1; dy <= 1; dy++)
            {
                for (var dx = -1; dx <= 1; dx++)
                {
                    var nx = x + dx;
                    var ny = y + dy;

                    if (nx < 0 || ny < 0 || nx >= width || ny >= height)
                        continue;

                    var neighbour = ny * width + nx;

                    if (edges[neighbour] > 0 || thin[neighbour] < low)
                        continue;

                    edges[neighbour] = 1;
                    pending.Push(neighbour);
                }
            }
        }

        return edges;
    }
}

var photo = PngImage.Load(file);
var gray = photo.ToBrightness();

// Drei Einstellungen desselben Verfahrens. Das Weichzeichnen entscheidet, wie viel
// Detail überlebt, die beiden Schwellen, was noch als Kante zählt.
(string Name, float Blur, float Low, float High)[] settings =
[
    ("fine", 1.0f, 0.08f, 0.20f),
    ("calm", 2.0f, 0.08f, 0.20f),
    ("strict", 2.0f, 0.15f, 0.35f),
];

Console.WriteLine("setting  blur   low  high  edge pixels");

foreach (var (name, blur, low, high) in settings)
{
    var edges = Canny.Detect(gray, photo.Width, photo.Height, blur, low, high);
    var drawn = new float[edges.Length];

    for (var i = 0; i < edges.Length; i++)
        drawn[i] = 1 - edges[i];

    Channels.ToImage(drawn, photo.Width, photo.Height)
        .Save(Path.Combine(folder, $"spring-canny-{name}.png"));

    Console.WriteLine(
        $"{name,-7}  {blur,4:0.0}  {low,4:0.00}  {high,4:0.00}  " +
        $"{edges.Sum() * 100 / edges.Length,10:0.0} %");
}