Theoretische und methodische Grundlagen des Visual Computing – Interpolation unregelmäßig verteilter Daten

Voronoi-Diagramm und Delaunay-Triangulierung

Voronoi-Diagramm: Gegeben sei eine Menge \(\S := \{\vecs {x}{1}, \dotsc , \vecs {x}{N}\} \subset \real ^n\) von \(N\) Punkten und eine Metrik \(d\) auf \(\real ^n\). Die Voronoi-Zelle von \(\vecs {x}{i}\) ist definiert als
\(V(\vecs {x}{i}) := \{\vec {x} \in \real ^n \;|\; \forall _{j\not =i}\; d(\vec {x}, \vecs {x}{i}) < d(\vec {x}, \vecs {x}{j})\}\).
Das Voronoi-Diagramm von \(\S \) ist \(\Vor (\S ) := \{V(\vecs {x}{i})\}_{i=1}^N\).

Delaunay-Triangulierung: Die Delaunay-Triangulierung \(\Del (\S )\) ist dual zu \(\Vor (\S )\).
Zwei Punkte \(\vecs {x}{i}\), \(\vecs {x}{j}\) werden verbunden, falls \(V(\vecs {x}{i})\) und \(V(\vecs {x}{j})\) eine gemeinsame Kante besitzen. \(\vecs {x}{i}, \vecs {x}{j}, \vecs {x}{k}\) bilden ein Dreieck genau dann, wenn kein anderer Punkt aus \(\S \) im Umkreis des Dreiecks \(\Delta \vecs {x}{i}\vecs {x}{j}\vecs {x}{k}\) liegt. \(\Del (\S )\) maximiert lexikografisch die minimalen Innenwinkel der Triang. und außerdem das Verhältnis \(a := \frac {r_\text {Inkreis}}{r_\text {Umkreis}}\).

Bowyer-Watson-Algorithmus: Der Bowyer-Watson-Algorithmus konstruiert \(\Del (\S )\) inkrementell. Um einen Punkt \(\vecs {x}{i}\) in die aktuelle Delaunay-Triang. der Punkte \(\vecs {x}{1}, \dotsc , \vecs {x}{i-1}\) einzufügen, wird zunächst betrachtet, in welchen Dreiecksumkreisen \(\vecs {x}{i}\) liegt. Diese Dreiecke werden dann aus der aktuellen Delaunay-Triang. gelöscht, was ein sternförmiges Gebiet erzeugt. Verbinde nun den Punkt \(\vecs {x}{i}\) mit allen Ecken, um die neue Delaunay-Triang. zu erhalten.

Seien nun Datenwerte \(f_1, \dotsc , f_N\) an den \(\vecs {x}{1}, \dotsc , \vecs {x}{N}\) gegeben.

Interpolation mit Voronoi-Zerlegung: Mit der Voronoi-Zerlegung kann man interpolieren, indem man jeder Zelle \(V(\vecs {x}{i})\) den Wert \(f_i\) zuweist (nächster Nachbar).

Interpolation mit Delaunay-Triangulierung: Mit der Delaunay-Triang. kann man interpolieren, indem man lineare Interpolation in jedem Dreieck durchführt.

Shepard-Interpolation

Shepard-Interpolation: Die Shepard-Interpolation (oder inverse Distanzwichtung) interpoliert \(f_i \in \real \), \(\vecs {x}{i} := (x_i, y_i) \in \real ^2\) (\(i = 1, \dotsc , N\)) mit dem Ansatz \(f(x, y) := \sum _{j=1}^N f_j w_j(x, y)\). Damit erhält man \(\forall _{i=1,\dotsc ,N}\; f_i = f(x_i, y_i) \iff \forall _{i,j=1,\dotsc ,N}\; w_j(x_i, y_i) = \delta _{j,i}\).
Dies motiviert den Ansatz \(w_j(x, y) := \frac {\varphi _j(x, y)^\mu }{\sum _{k=1}^N \varphi _k(x, y)^\mu }\) mit \(\varphi _j(x, y) := \frac {1}{|\vec {x} - \vecs {x}{j}|}\) und \(\mu \in (0, \infty )\).
Es gilt \(w_j(x, y) \in [0, 1]\), \(\sum _{j=1}^N w_j(x, y) = 1\) und \(w_j(x_i, y_i) = \delta _{j,i}\)
(wenn man \(\varphi _j(x_j, y_j) := \infty \) und \(w_j(x_j, y_j) := 1\) setzt).

lokale Shepard-Interpolation: Bei der lokalen Shepard-Interpolation verwendet man
Gewichtsfunktionen \(w_j\), die nur von einer lokalen Teilmenge von Samplepunkten abhängen. Definiert man \(r_j = r_j(x, y) := |\vec {x} - \vecs {x}{j}|\), dann kann man z. B.

\(\varphi _j(x, y) := \left (\frac {R - r_j}{R \cdot r_j}\right )^2\) für \(r_j \in [0, R]\) (modifizierte Shepard-Gewichte),
\(\varphi _j(x, y) := \frac {1}{r_j}\) für \(r_j \in [0, \frac {R}{3}]\) und \(\varphi _j(x, y) := \frac {27}{4R} \left (\frac {r_j}{R} - 1\right )^2\) für \(r_j \in [\frac {R}{3}, R]\) oder alternativ
\(\varphi _j(x, y) := 1 - \frac {r_j}{R}\) für \(r_j \in [0, R]\) (Franke-Little-Gewichte) verwenden.

Methode der radialen Basisfunktionen

Methode der radialen Basisfunktionen: Seien \(N\) paarw. disj. Datenpunkte \(\vecs {x}{i} \in \real ^n\) und Werte \(f_i \in \real \) für \(i = 1, \dotsc , N\) gegeben. Die Methode der radialen Basisfunktionen arbeitet ähnlich wie die Shepard-Intp., mit dem Unterschied, dass hier \(f(\vec {x}) := \sum _{i=1}^N \lambda _i \phi (|\vec {x} - \vecs {x}{i}|)\) mit einer radialen Basisfunktion (RBF) \(\phi (r)\colon [0, \infty ) \to \real \) angesetzt wird (d. h. die Koeffizienten sind jetzt allgemein und die Funktionen sind zwingend radialsymmetrisch).

gebräuchliche Basisfunktionen: \(\phi (r) := e^{-(\varepsilon r)^2}\) (Gauß (GA)),
\(\phi (r) := \frac {1}{1 + (\varepsilon r)^2}\) (invers quadratisch (IQ)), \(\phi (r) := \frac {1}{\sqrt {1 + (\varepsilon r)^2}}\) (invers multiquadratisch (IMQ)),
\(\phi (r) := \sqrt {1 + (\varepsilon r)^2}\) (Hardy-multiquadratisch (MQ)), \(\phi (r) := r\) (linear), \(\phi (r) := r^3\) (kubisch),
\(\phi (r) := r^2 \ln r\) (Thin Plate Spline (TPS))

Ermitteln der Koeffizienten: Die Koeffizienten bekommt man aus dem LGS
\(\forall _{i=1,\dotsc ,N}\; f(\vecs {x}{i}) = f_i \iff \Phi \vec {\lambda } = \vec {f}\) mit \(\Phi := (\phi _{i,j})_{i,j=1}^N\), \(\phi _{i,j} := \phi (|\vecs {x}{i} - \vecs {x}{j}|)\), \(\vec {\lambda } := (\lambda _j)_{j=1}^N\) und \(\vec {f} := (f_i)_{i=1}^N\). Ist \(\Phi \) invertierbar, so gilt \(\vec {\lambda } = \Phi ^{-1} \vec {f}\).

Im Folgenden werden hinreichende Bedingungen für die Invertierbarkeit von \(\Phi \) angegeben.

vollständig monoton: Eine Funktion \(\psi \colon [0, \infty ) \to \real \) heißt vollständig monoton, falls

\(\psi \in \C ^0([0, \infty ))\),
\(\psi \in \C ^\infty ((0, \infty ))\) und
\(\forall _{\ell \in \natural _0} \forall _{r > 0}\; (-1)^\ell \psi ^{(\ell )}(r) \ge 0\).

Satz von Schoenberg: Sei \(\psi (r) := \phi (\sqrt {r})\). Ist \(\psi \) vollständig monoton, aber nicht konstant auf \([0, \infty )\), dann ist \(\Phi := (\phi _{i,j})_{i,j=1}^N\) für paarw. disjunkte Punkte \(\vecs {x}{i}\) p. d. und insb. invertierbar.

Beispiel: Für die Gauß-RBF ist \(\psi (r) = e^{-\varepsilon ^2 r}\) und Schoenberg ist anwendbar (\((-1)^\ell \psi ^{(\ell )}(r) > 0\)).
Für die MQ-RBF ist \(\psi (r) = \sqrt {1 + \varepsilon ^2 r}\) und Schoenberg ist nicht anwendbar (\(-\psi ’(r) < 0\)).

Satz 1 von Micchelli: Sei \(\psi (r) := \phi (\sqrt {r})\). Ist \(\psi \in \C ^0([0, \infty ))\), \(\forall _{r>0}\; \psi (r) > 0\) und \(\psi ’\) vollständig monoton, aber nicht konstant auf \([0, \infty )\), dann ist \(\Phi := (\phi _{i,j})_{i,j=1}^N\) für paarw. disjunkte Punkte \(\vecs {x}{i}\) invertierbar.

Beispiel: Für die MQ-RBF ist der Satz anwendbar. Allerdings kann man den Satz für TPS nicht anwenden, weil für \(\psi (r) = r \ln \sqrt {r}\) gilt, dass \(\psi (r) < 0\) für \(r > 0\) klein (Schoenberg geht auch nicht, weil \(-\psi ’(r) = -\frac {1}{2} - \ln \sqrt {r} < 0\) für \(r\) groß).

augmentierte RBF-Interpolation: Sei \(\{p_k\}_{k=1}^M\) eine Basis von \(\PP _m(\real ^d)\) (\(d\)-variate Polynome vom Grad \(\le m\), \(M := \binom {m+d}{m}\)) und \(m \in \natural _0\). Dann ist der Ansatz für die augmentierte RBF-Interpolation \(f(\vec {x}) = \sum _{i=1}^N \lambda _j \phi (|\vec {x} - \vecs {x}{j}|) + \sum _{k=1}^M \gamma _k p_k(\vec {x})\). Damit die Interpolation nicht unterbestimmt ist, legt man die zusätzlichen Bedingungen \(\forall _{k=1,\dotsc ,M}\; \sum _{j=1}^N \lambda _j p_k(\vecs {x}{j}) = 0\) fest.
Man erhält das Interpolations-LGS \(A\smallpmatrix {\vec {\lambda }\\\vec {\gamma }} = \smallpmatrix {\vec {f}\\\vec {0}}\) mit \(A := \smallpmatrix {\Phi &P\\P^\tp &0}\), \(\Phi := (\phi (|\vecs {x}{j} - \vecs {x}{k}|))_{j,k=1}^N\) und \(P := (p_k(\vecs {x}{j}))_{j,k=1}^{N,M}\).

Satz 2 von Micchelli: Seien \(\psi (r) := \phi (\sqrt {r})\) und \(m \in \natural _0\). Ist \(\psi \in \C ^0([0, \infty ))\), \(\psi ^{(m+1)}\) vollständig monoton, aber nicht konstant auf \([0, \infty )\), und \(\Rang (P) = M\) (d. h. \(P\) hat vollen Spaltenrang), dann ist \(A\) für paarw. disjunkte Punkte \(\vecs {x}{i}\) invertierbar.