5.4 Delaunay-Algorithmus

Das ist der Inhalt von Praktikumsaufgabe 4 (randomisierte Delaunay-Triangulierung).

Idee (Folie 40)

Ein einfacher Algorithmus in drei Schritten:

1. Die konvexe Hülle $\mathcal{CH}(P)$ berechnen (Kapitel 2).
2. Einen zufälligen Punkt $p\in P\setminus\mathcal{CH}(P)$ ins Innere einfügen: ihn mit allen Ecken der Hülle verbinden (Fächer) und die Delaunay-Bedingung wiederherstellen.
3. Die übrigen Punkte in zufälliger Reihenfolge einfügen und jeweils die Delaunay-Bedingung wiederherstellen.

Liegen mehrere Punkte auf einem Kreis, kann die entsprechende Fläche beliebig trianguliert werden. Dann gibt es mehr als eine Delaunay-Triangulierung.

Algorithmus 4 (Folie 41)

Algorithmus 4: Delaunay-Triangulierung
Eingabe: Menge P von n verschiedenen Punkten in der Ebene
Ausgabe: Delaunay-Triangulierung D
1: Berechne die konvexe Hülle CH(P);
2: Berechne eine zufällige Permutation p_1, …, p_m
   der übrigen Punkte P \ CH(P);
3: Initiale Triangulierung D: p_1 mit allen Ecken von CH(P) verbinden
   und die Delaunay-Bedingung wiederherstellen;
4: for (r = 2, …, m) {
5:   Finde das Dreieck p_i p_j p_k, das p_r enthält;
6:   Finde alle Dreiecke, deren Delaunay-Bedingung verletzt ist
     (über die Nachbarschaft der Dreiecke, ausgehend vom Fund aus 5);
7:   Entferne diese Dreiecke aus D → polygonales Loch Δ;
8:   Verbinde alle Ecken von Δ mit p_r (Loch neu triangulieren);
9: }

Das ist die Bowyer–Watson-Variante: Statt einzeln Kanten zu flippen, werden alle „schlechten“ Dreiecke auf einmal entfernt, und das Loch wird sternförmig neu gefüllt.

• Schlechtes Dreieck = eines, dessen Umkreis $p_r$ enthält.
• Die schlechten Dreiecke bilden immer eine zusammenhängende Menge um das Dreieck, das $p_r$ enthält. Deshalb reicht in Zeile 6 eine Breitensuche über die Nachbarn: Man geht nur zu Nachbarn weiter, die selbst schlecht sind.
• Das Loch $\Delta$ ist sternförmig bezüglich $p_r$: Von $p_r$ aus sieht man jede Ecke von $\Delta$. Deshalb sind die Fächerkanten in Zeile 8 alle gültig.

Schritt 1/16: Zeile 1: konvexe Hülle CH(P) berechnen (6 Ecken). Die übrigen 4 Punkte werden zufällig permutiert (Zeile 2).

langsamnormalschnellsehr schnell

neuer Punkt / enthaltendes Dreieckverletzt Delaunayneue Dreiecke

Satz 11: Korrektheit (Folien 42–43)

Satz 11

Algorithmus 4 berechnet eine korrekte Delaunay-Triangulierung.

Beweis. Zu zeigen ist, dass beim Einfügen von $p_r$ (Zeilen 7 und 8)

1. alle neuen Kanten im Inneren des Lochs $\Delta$ liegen und
2. alle neuen Dreiecke die Delaunay-Bedingung erfüllen.

Zu (1) (Folien: „Tafel“, hier ausgeführt): Für jede Randkante $\overline{p_ap_b}$ von $\Delta$ gehörte das angrenzende entfernte Dreieck $t$ zu den schlechten, $p_r$ liegt also im Umkreis von $t$. Das Dreieck auf der anderen Seite der Kante ist gut, $p_r$ liegt also nicht in dessen Umkreis. Daraus folgt, dass $p_r$ auf derselben Seite von $\overline{p_ap_b}$ liegt wie $t$. Gilt das für alle Randkanten, sieht $p_r$ jede Kante von innen: $\Delta$ ist sternförmig bzgl. $p_r$, und die Strecken $\overline{p_rp_a}$ liegen in $\Delta$.

Zu (2) (Widerspruch): Angenommen, ein neues Dreieck $t$ mit Kante $\overline{p_rp_j}$ (aus Zeile 8) verletzt die Delaunay-Bedingung. Dann müsste $\overline{p_rp_j}$ geflippt werden. Dabei entsteht ein Dreieck $\Gamma = p_ip_jp_k$ aus drei aufeinanderfolgenden Ecken von $\Delta$, dessen Umkreis $p_r$ nicht enthält.

• a) War $\Gamma\in D$ vor dem Einfügen, wurde es entfernt, also enthielt sein Umkreis $p_r$. Widerspruch.
• b) War $\Gamma\notin D$ vor dem Einfügen, hat es die Delaunay-Bedingung schon vorher nicht erfüllt (sonst wäre es ja ein Delaunay-Dreieck gewesen). Ein Punkt außer $p_r$ liegt in seinem Umkreis. Widerspruch dazu, dass der Flip es legal machen soll.

Beide Fälle sind unmöglich, also ist jedes neue Dreieck Delaunay. Per Induktion über $r$ ist $D$ nach jedem Schritt eine Delaunay-Triangulierung. □

Laufzeit (Folie 45)

	Aufwand
naiv (Dreieck in Zeile 5 durch lineare Suche)	$O(r)$ pro Schritt ⇒ $O(n^2)$
mit Point-Location-Struktur (vgl. Kap. 6, Trapezzerlegung bzw. History-DAG)	erwartet $O(\log n)$ pro Suche
Änderungen pro Schritt (entfernte/neue Dreiecke)	erwartet konstant ($\le 6$)
gesamt	erwartet $O(n\log n)$
Speicher	$O(n)$ in jedem Fall

Die Kantenflips lassen sich in der Point-Location-Struktur lokal in erwartet $O(\log n)$ nachführen.

Warum erwartet $\le 6$? (Rückwärtsanalyse, Ergänzung)

Nach dem Einfügen von $p_r$ hat $p_r$ in der Triangulierung von $r$ Punkten einen gewissen Grad $\deg(p_r)$. Die Zahl neuer Dreiecke ist genau $\deg(p_r)$. Stellt man sich vor, der zuletzt eingefügte Punkt ist ein zufälliger der $r$ Punkte (wegen der zufälligen Permutation), ist sein erwarteter Grad der Durchschnittsgrad:

$\mathbb E[\deg(p_r)] = \frac{2\cdot\#\text{Kanten}}{r} \le \frac{2(3r-6)}{r} < 6.$

Deshalb braucht man die zufällige Reihenfolge: Für eine böse Reihenfolge (z.B. Punkte auf einer Parabel von links nach rechts) kann jeder Schritt $\Theta(r)$ Dreiecke ändern, insgesamt $\Theta(n^2)$.

Variante mit Super-Dreieck und Flips (Ergänzung, so in der Demo)

Die Folien starten mit der konvexen Hülle. Häufiger in Implementierungen:

• Super-Dreieck: Ein riesiges Dreieck um alle Punkte. Dann liegt jeder Punkt in einem Dreieck, ohne Hülle vorab. Am Ende alle Dreiecke mit einer Super-Ecke löschen. Vorsicht: Ist das Super-Dreieck zu klein, fehlen am Rand Kanten der Hülle.
• Lawson-Flips statt Loch: $p_r$ ins Dreieck einfügen (in 3 Dreiecke teilen) und dann die drei gegenüberliegenden Kanten prüfen. Ist eine illegal, flippen und die zwei neuen gegenüberliegenden Kanten auf einen Stack legen. Ergibt dieselbe Triangulierung wie Bowyer–Watson.

C++-Skizze (Bowyer–Watson, naiv O(n²))

struct P { double x, y; };
struct Tri { int a, b, c; };   // Indizes in pts, gegen den Uhrzeigersinn

// > 0  ⇔  d liegt im Umkreis des ccw-Dreiecks (a, b, c)
double inCircle(P a, P b, P c, P d) {
    double ax = a.x - d.x, ay = a.y - d.y;
    double bx = b.x - d.x, by = b.y - d.y;
    double cx = c.x - d.x, cy = c.y - d.y;
    return (ax*ax + ay*ay) * (bx*cy - cx*by)
         - (bx*bx + by*by) * (ax*cy - cx*ay)
         + (cx*cx + cy*cy) * (ax*by - bx*ay);
}

std::vector<Tri> delaunay(std::vector<P> pts) {
    std::shuffle(pts.begin(), pts.end(), std::mt19937{42});
    int n = pts.size();
    const double R = 1e6;                       // Super-Dreieck
    pts.push_back({-R, -R}); pts.push_back({R, -R}); pts.push_back({0, R});
    std::vector<Tri> T{{n, n + 1, n + 2}};

    for (int r = 0; r < n; ++r) {
        std::vector<Tri> keep;
        std::map<std::pair<int,int>, int> edgeCount;  // Randkanten des Lochs
        for (auto t : T) {
            if (inCircle(pts[t.a], pts[t.b], pts[t.c], pts[r]) > 0) {
                for (auto [u, v] : {std::pair{t.a, t.b}, {t.b, t.c}, {t.c, t.a}})
                    edgeCount[{std::min(u, v), std::max(u, v)}]++;
            } else keep.push_back(t);
        }
        // Kanten, die nur zu einem schlechten Dreieck gehören, bilden den Rand
        for (auto t : T)
            if (inCircle(pts[t.a], pts[t.b], pts[t.c], pts[r]) > 0)
                for (auto [u, v] : {std::pair{t.a, t.b}, {t.b, t.c}, {t.c, t.a}})
                    if (edgeCount[{std::min(u, v), std::max(u, v)}] == 1)
                        keep.push_back({u, v, r});   // ccw bleibt erhalten
        T = std::move(keep);
    }
    std::erase_if(T, [n](Tri t) { return t.a >= n || t.b >= n || t.c >= n; });
    return T;
}

Für $O(n\log n)$ braucht man statt der Schleife über alle Dreiecke eine Point-Location (History-DAG) und die Breitensuche über Dreiecks-Nachbarn.