Drei Punkte, eine falsche Antwort: Shewchuks robuster Orientierungstest in Grasshopper C#
Warum 'links der Linie?' in Floating Point falsch sein kann und wie man Shewchuks robustes Orientierungstest in Grasshopper C# baut.
Drei Punkte auf dem Tisch. Drehen sie nach links, nach rechts oder liegen sie auf einer Linie? Es ist die kleinste Frage in der Geometrie, und deine Software stellt sie häufiger als jede andere. Jede boolesche Operation, konvexe Hülle, Mesh-Schnitt und Delaunay-Triangulation wird aus Millionen dieser Ja/Nein-Antworten gebaut. Computergeomet nennen sie Prädikate. Der Orientierungstest ist der einfachste: das Vorzeichen einer 2×2-Determinante, (qx−px)·(ry−py) − (qy−py)·(rx−px). Positiv heisst links, negativ heisst rechts, null heisst kollinear.
Und hier die Überraschung. Auf einem Computer kommt die Antwort falsch zurück. Nicht knapp daneben. Falsch. Koordinaten sind als IEEE-754-Doubles gespeichert, eine 53-Bit-Mantisse mal eine Zweierpotenz, und jede Subtraktion und Multiplikation rundet. Wenn die Punkte fast kollinear sind, sind die zwei Produkte fast gleich, und der Rundungsfehler kann grösser sein als der echte Unterschied. Das Vorzeichen ist die ganze Antwort, und es kommt als Null zurück oder kippt das Vorzeichen.
Die Treppe, wo eine Linie sein sollte
Lutz Kettner, Kurt Mehlhorn, Sylvain Pion, Stefan Schirra und Chee Yap machten das Problem in einem Lehrbeitrag sichtbar: Classroom Examples of Robustness Problems in Geometric Computations (Computational Geometry: Theory and Applications 40(1), 2008), vom Max-Planck-Institut für Informatik, INRIA Sophia Antipolis, Otto-von-Guericke-Universität Magdeburg und New York University. Sie nahmen die Diagonale durch (12, 12) und (24, 24), testeten Abfragepunkte nahe (0.5, 0.5) in Schritten von 2−53 über ein 256×256-Gitter und färbten jeden Punkt nach der Floating-Point-Antwort.
Kettner et al. schreiben, dass sie “ein gelbes Band um die Diagonale mit fast geraden Grenzen” erwarteten. Sie fanden eine Treppe aus Blöcken, mit “Punkten, die die Linienseite wechseln”. Der Grund: 12 hat vier Binärstellen vor dem Komma, also wirft die Subtraktion die letzten vier Bits des Schritts weg, und das Ergebnis bleibt für 24 aufeinanderfolgende Schritte gleich (25 für 24). Niemand hat Treppen bestellt.
PAZ wiederholte das Setup am 2026-09-28, einmal in Python mit exakten Brüchen und einmal mit dem C# unten; beide gaben identische Zählungen. Der naive Double-Test unterscheidet sich vom exakten Ergebnis bei 11.972 von 65.536 Punkten, etwa 18%: 11.300 falsch als kollinear bezeichnet, und 672 auf die falsche Seite gestellt. Nur 256 Gitterpunkte liegen wirklich auf der Linie. Das sind PAZs Zahlen, nicht Zahlen aus dem Papier.
Ein falsches Vorzeichen ist eine falsche Antwort. Das Problem ist, dass Algorithmen viele Antworten kombinieren und erwarten, dass sie übereinstimmen. Wenn sie nicht übereinstimmen, ist das Ergebnis topologisch unmöglich. Das Lehrbeitrag zeigt Double-Genauigkeits-Konvexhüllen, die einen extremen Punkt auslassen oder sich selbst schneiden, und Implementierungen, die “möglicherweise sogar für immer laufen”. Nach dem CGAL-Kernel-Handbuch “kapseln Prädikate Entscheidungen”, und ungenaue Arithmetik “kann zu inkonsistenten Entscheidungen führen, was zu unerwarteten Fehlschlägen bei einigen korrekten Eingabedaten führt”.
Die Lösung: Billig zuerst, exakt nur wenn es zählt
Jonathan Richard Shewchuks Antwort, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates (Discrete & Computational Geometry 18(3):305–363, 1997, geschrieben an der Carnegie Mellon; er ist jetzt Professor an der UC Berkeley), passt in einen Satz. Führe die billige Berechnung durch, arbeite aus, wie falsch sie möglicherweise sein könnte, und mache mehr Arbeit nur, wenn dieser Fehler das Vorzeichen kippen könnte.
Zwei Ideen tragen sie. Eine exakte Zahl kann als Expansion gespeichert werden, eine Summe nicht überlappender Doubles: der Rundungsfehler von a + b ist selbst genau ein Double, das man in wenigen Operationen zurückgewinnen kann (die TWO-SUM-Abstammungslinie von Dekker, Knuth und Priest). Und die Berechnung läuft in vier Stufen. Stufe A ist die einfache Determinante plus eine Fehlergrenze-Prüfung; wenn das Ergebnis die Grenze schlägt, ist das Vorzeichen garantiert und die Arbeit stoppt. Stufen B und C fügen Korrektionen hinzu, Stufe D berechnet den exakten Wert, jede nutzt die letzte.
Shewchuk mass die Auszahlung auf einer 2D-Delaunay-Triangulation eines geneigten Millionen-Punkte-Gitters, einer seiner schlimmsten Eingaben: etwa 9,4 Millionen Orientierungsaufrufe, 9.318.610 in Stufe A settled, 121.081 in B, 118 in C, und nur 3 mit voller exakter Arithmetik. Nach Shewchuk beträgt die Kosten etwa 8% bei Zufallspunkten und bis zu 30% bei schwierigen, “genau für die Punktsätze, die am wahrscheinlichsten Schwierigkeiten verursachen”. Der Kompromiss ist real: In 3D kosten robuste Prädikate etwa 35% bei Zufallseingaben und einen Faktor von elf bei Punkten nahe einer Kugel. Aber sein 3D-Mesher Pyramid, auf dem geneigten Gitter ohne robuste Arithmetik ausgeführt, “konnte nicht beendet werden”, gefangen in einer Endlosschleife durch eine Rundungsinkonsistenz. Sein Urteil: “Robuste Arithmetik ist schliesslich nicht immer langsamer.”
Schau, wer den Kredit erhält. Shewchuks Danksagungen nennen Steven Fortune, Douglas Priest und Christopher Van Wyk als die Grundlagen, und sagen, Fortune “hat diese Forschung Mitte 1994 unwissentlich mit ein paar kurzen E-Mail-Antworten ausgelöst.” Ein paar E-Mails wurden zu Code, der jetzt, grösstenteils unsichtbar, unter einer Menge Geometrie-Software sitzt.
←HEUTE: Shewchuks öffentlich zugängliche predicates.c von 1997 ist Gemeingut und wird immer noch ausgeliefert, direkt oder als Ports, einschliesslich Vladimir Agafonkins JavaScript robust-predicates. →3012: Zürichs gedrehte Gitter und ausgerichtete Fassaden hängen immer noch von der gleichen vierstufigen Vorzeichenentscheidung ab. Drehpunkt: Code, den jeder noch kompilieren kann, überlebt jeden Anbieter.
Warum ein Gebäude die schlimmstmögliche Eingabe ist
Gebäude sind voller degenerierter Fälle mit Absicht: ausgerichtete Wände, koplanare Platten, wiederholte Module, Teile, die genau aufeinandertreffen. Ein strukturelles Gitter, das gedreht ist, um einer Grundstücksgrenze zu folgen, ist das geneigte Gitter aus beiden Papieren. Rhino und Grasshopper handhaben das mit der absoluten Toleranz des Dokuments: Vereinbare im Voraus, welche Distanz als “derselbe Punkt” zählt. Toleranz verwaltet Konstruktionen, wo ein Punkt landet. Robuste Prädikate halten Entscheidungen konsistent, welche Seite er landet. Zwei Schichten eines Problems, und deine eigenen Skripte brauchen beide.
Das Werkzeug: Shewchuks adaptive Prädikate (orient2d, orient3d, incircle, insphere), öffentlich zugängliches C in predicates.c, portiert zu JavaScript von Vladimir Agafonkin als robust-predicates (Unlicense), mit der gleichen Filter-First-Idee in CGALs gefilterten Kerneln. Heute bauen wir eine kleine C#-Lehrbeitrag-Version: ein naiver Fast-Test und ein Brute-Force-Exact-Test. Exact ist nur Shewchuks Stufe D; sein echter Beitrag ist, dass du sie fast nie laufen musst.
Setup:
dotnet new console -n OrientLab
cd OrientLab
# paste the Orient class and the loop into Program.cs
dotnet run
# expected: 11972 wrong, 672 on the wrong sideusing System;
using System.Numerics;
static class Orient
{
public static int Fast(double px, double py, double qx, double qy, double rx, double ry)
=> Math.Sign((qx - px) * (ry - py) - (qy - py) * (rx - px));
// Every double is m * 2^e: scale all six to integers, use BigInteger.
public static int Exact(double px, double py, double qx, double qy, double rx, double ry)
{
double[] v = { px, py, qx, qy, rx, ry };
var m = new BigInteger[6]; var e = new int[6]; int minE = int.MaxValue;
for (int i = 0; i < 6; i++)
{
long bits = BitConverter.DoubleToInt64Bits(v[i]);
int exp = (int)((bits >> 52) & 0x7FF);
long man = bits & 0xFFFFFFFFFFFFFL;
if (exp == 0) exp = 1; else man |= 1L << 52;
m[i] = bits < 0 ? -man : man; e[i] = exp - 1075;
if (man != 0) minE = Math.Min(minE, e[i]);
}
if (minE == int.MaxValue) return 0;
for (int i = 0; i < 6; i++) m[i] <<= e[i] - minE;
BigInteger d = (m[2] - m[0]) * (m[5] - m[1]) - (m[3] - m[1]) * (m[4] - m[0]);
return d.Sign;
}
}Erste Schritte:
- In Rhino 8 lege eine C#-Skript-Komponente auf die Grasshopper-Leinwand und füge die
Orient-Klasse unter der Skript-Klasse nebenRunScriptein. - Füge die Schleife in
RunScriptein:double u = Math.Pow(2, -53); int wrong = 0, flipped = 0; for (int X = 0; X < 256; X++) for (int Y = 0; Y < 256; Y++) { double px = 0.5 + X * u, py = 0.5 + Y * u; int f = Orient.Fast(px, py, 12, 12, 24, 24); int x = Orient.Exact(px, py, 12, 12, 24, 24); if (f != x) wrong++; if (f != 0 && f == -x) flipped++; } Print($"{wrong} wrong, {flipped} on the wrong side"); - Lies das Ergebnis: 11972 falsch, 672 auf der falschen Seite.
- Dehnziel: gib einen gefärbten
Point3d(X, Y, 0)pro Zelle nachFast-Antwort aus und backe. Gitter-Indizes plotten, nicht echte Koordinaten, die 2−53 auseinander liegen, weit unter jeder Modellierungstoleranz. Die Treppe erscheint.
Atelier: Ateliers, die ihre eigenen Grasshopper-C#- oder Python-Komponenten schreiben, tragen normalerweise ein paar Zeilen wie if (cross == 0) herum, die entscheiden, ob eine Spalte auf einer Gitterlinie sitzt oder eine Raumkante schliesst. Sie übergeben das Test-Modell und brechen auf dem gedrehten Flügel. Der Montagzug: Deine Skript-Bibliothek nach Null- oder Vorzeichen-Tests auf rohen Kreuzprodukten durchsuchen und jeden durch ein gefiltertes Prädikat wie den Hack unten leiten, bevor das nächste Tender-Stage-Modell läuft.
Hack: Gatter den Brute-Force-Test hinter Shewchuks Stufe A, so dass du für Genauigkeit nur zahlst, wenn die Grenze es verlangt. Die Konstante kommt aus seiner öffentlich zugänglichen predicates.c, mit Epsilon 2−53. Wenn die billige Determinante den schlimmsten Rundungsfehler schlägt, ist sein Vorzeichen sicher.
double l = (qx - px) * (ry - py), r2 = (qy - py) * (rx - px), det = l - r2;
double eps = Math.Pow(2, -53), bound = (3.0 + 16.0 * eps) * eps * (Math.Abs(l) + Math.Abs(r2));
if (Math.Abs(det) > bound) return Math.Sign(det); // stage A: sign is certain
return Orient.Exact(px, py, qx, qy, rx, ry); // rare: go exactVon wo aus ich in den späten 2070ern sitze, war die Geometrie, die arbeitete, die Geometrie, deren Quelle jeder noch kompilieren konnte. Schreibe Atelier-Skripte auf die gleiche Weise: lesbar, offen, reparierbar von einem 25-Jährigen. Zuversicht entscheidet nicht das Vorzeichen.
Lerne es:
- Shewchuk, Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates, DCG 18(3), 1997 (CMU-CS-96-140R). Abschnitt 4 behandelt die Stufen.
- Kettner, Mehlhorn, Pion, Schirra, Yap, Classroom Examples of Robustness Problems in Geometric Computations, CGTA 40(1):61–78, 2008. Die Treppe ist in Abschnitt 3.
- Vladimir Agafonkins robust-predicates: ein lesbarer JavaScript-Port mit schnellen Varianten.
- CGAL-Kernel-Handbuch,
Exact_predicates_inexact_constructions_kernel. - PAZ-Notiz: PAZ Academy unterrichtet seit der Alpha Grasshopper 2 und C#-Skripting. Bring die Geometrieskripte deines Ateliers mit und wir führen sie zusammen gegen ein geneigtes Gitter aus.
Füge die Klasse ein, führe die 256×256-Schleife aus, und beobachte 672 Punkte auf der falschen Seite landen, bevor du dem nächsten == 0 vertraust.
PAZ Kaffi · interdisziplinäre Redaktionsarbeit, geleitet von der PAZ Academy