Marching Squares

// Ein Stück einer Höhenlinie, in den Koordinaten des Feldes, in dem es gefunden wurde.
public readonly record struct Segment(float X1, float Y1, float X2, float Y2);

// Ein Feld aus Zahlen - Höhen, Temperaturen, Helligkeiten - und die Frage, wo genau
// ein bestimmter Wert hindurchläuft. Die Antwort ist eine Linie, und der Trick, sie
// zu finden, ist, immer nur vier benachbarte Werte anzuschauen. Jeder davon liegt
// entweder über dem Wert oder darunter, das ergibt sechzehn Eckmuster,
// und für jedes davon steht fest, wie die Linie durch dieses Quadrat läuft.
public static class MarchingSquares
{
    public static Segment[] Contour(float[] values, int width, int height, float level)
    {
        var segments = new List<Segment>();

        for (var y = 0; y < height - 1; y++)
        {
            for (var x = 0; x < width - 1; x++)
            {
                var topLeft = values[y * width + x];
                var topRight = values[y * width + x + 1];
                var bottomRight = values[(y + 1) * width + x + 1];
                var bottomLeft = values[(y + 1) * width + x];

                // Ein Bit je Ecke: Liegt sie über dem Wert oder nicht?
                var corners = (topLeft > level ? 8 : 0)
                    | (topRight > level ? 4 : 0)
                    | (bottomRight > level ? 2 : 0)
                    | (bottomLeft > level ? 1 : 0);

                if (corners is 0 or 15)
                    continue;

                // Wo die Linie eine Kante schneidet, wird nicht geraten, sondern interpoliert:
                // Stehen an den Ecken 0,2 und 0,6 und gesucht ist 0,3, dann liegt der Schnitt
                // bei einem Viertel. Genau das macht die Höhenlinie glatt
                // statt treppig.
                var top = new Segment(x + Between(topLeft, topRight, level), y, 0, 0);
                var bottom = new Segment(x + Between(bottomLeft, bottomRight, level), y + 1, 0, 0);
                var left = new Segment(x, y + Between(topLeft, bottomLeft, level), 0, 0);
                var right = new Segment(x + 1, y + Between(topRight, bottomRight, level), 0, 0);

                switch (corners)
                {
                    case 1 or 14: Add(segments, left, bottom); break;
                    case 2 or 13: Add(segments, bottom, right); break;
                    case 3 or 12: Add(segments, left, right); break;
                    case 4 or 11: Add(segments, top, right); break;
                    case 6 or 9: Add(segments, top, bottom); break;
                    case 7 or 8: Add(segments, left, top); break;

                    // Zwei diagonal gegenüberliegende Ecken liegen oben, die anderen beiden
                    // unten. Beide Arten, sie zu verbinden, sind zulässig, und der
                    // Mittelwert der vier entscheidet, welche genommen wird - so können
                    // sich benachbarte Zellen nicht widersprechen.
                    case 5 or 10:
                        var middle = (topLeft + topRight + bottomRight + bottomLeft) / 4;

                        // Im Fall 5 liegen die oberen Ecken rechts oben und links unten,
                        // im Fall 10 andersherum, deshalb fällt dieselbe Entscheidung
                        // spiegelverkehrt aus.
                        var joinLeftToTop = middle > level == (corners == 5);

                        if (joinLeftToTop)
                        {
                            Add(segments, left, top);
                            Add(segments, bottom, right);
                        }
                        else
                        {
                            Add(segments, left, bottom);
                            Add(segments, top, right);
                        }

                        break;
                }
            }
        }

        return [.. segments];
    }

    // Wie weit zwischen zwei Ecken der Wert liegt, von 0 bis 1.
    private static float Between(float from, float to, float level) =>
        Math.Abs(to - from) < 1e-6f ? 0.5f : Math.Clamp((level - from) / (to - from), 0f, 1f);

    private static void Add(List<Segment> segments, Segment from, Segment to)
    {
        // Läuft der Wert genau durch eine Ecke, können beide Schnitte auf derselben
        // Ecke landen. So ein Stück hat keine Länge und bleibt weg.
        if (from.X1 == to.X1 && from.Y1 == to.Y1)
            return;

        segments.Add(new Segment(from.X1, from.Y1, to.X1, to.Y1));
    }
}

const int width = 120;
const int height = 90;

// Eine Landschaft aus Perlin Noise und die Frage, wo bestimmte Höhen
// hindurchlaufen. Genau das sind die Höhenlinien einer Karte.
var random = new Pcg(20240921);
var seeds = new float[width * height];

for (var i = 0; i < seeds.Length; i++)
    seeds[i] = (float)random.NextDouble();

var land = PerlinNoise.Generate(seeds, width, height, octaves: 5, persistence: 0.5f);

// Perlin Noise nutzt nicht den ganzen Bereich von 0 bis 1, deshalb werden die
// Höhen zuerst gestreckt. Sonst fänden die unteren Höhenlinien nichts.
var lowest = land.Min();
var span = land.Max() - lowest;

for (var i = 0; i < land.Length; i++)
    land[i] = (land[i] - lowest) / span;

var lines = new List<Segment>();

Console.WriteLine("height  pieces of line");

foreach (var level in new[] { 0.2f, 0.35f, 0.5f, 0.65f, 0.8f })
{
    var contour = MarchingSquares.Contour(land, width, height, level);

    lines.AddRange(contour);

    Console.WriteLine($"{level,6:0.00}  {contour.Length,14}");
}

// Jedes Stück wird für sich gezeichnet. Zu geschlossenen Ringen sortiert werden
// sie hier nicht - für ein Bild braucht es das nicht, für einen Plotter schon.
ContourRenderer.Render([.. lines], width, height, scale: 4, background: land)
    .Save(Path.Combine(folder, "contour-lines.png"));

Channels.ToImage(land, width, height).Save(Path.Combine(folder, "contour-field.png"));

Console.WriteLine($"{lines.Count} pieces of line altogether");