Theoretische und methodische Grundlagen des Visual Computing – Interpolation auf Gittern

1D-Polynom-Interpolation

Satz (Polynom-Interpolation): Seien $(x_i,y_i)\in {\mathbb{R}}^2$ für $i=0,\dots ,n$ gegeben, wobei $x_i\ne x_j$ für $i\ne j$. Dann gibt es genau ein interpolierendes Polynom $P_n(x)=\sum _{j=0}^n{a}_j{x}^j$ vom Grad $\le n$.

Vandermonde-Matrix: Die Koeffizienten $a_0,\dots ,a_n$ lassen sich als Lösung des LGS
$\left(\begin{array}{cccc}1& x_0& \dots & x_0^n\\ ⋮& ⋮& & ⋮\\ 1& x_n& \dots & x_n^n\end{array}\right)\left(\begin{array}{c}{a}_0\\ ⋮\\ a_n\end{array}\right)=\left(\begin{array}{c}{y}_0\\ ⋮\\ y_n\end{array}\right)⟺V\overrightarrow{a}=\overrightarrow{y}$ bestimmen. $V$ heißt Vandermonde-Matrix. Allerdings ist $V$ schwierig und teuer zu berechnen und weitere Interpolationspkt.e sind aufwendig.

Lagrange-Interpolation: Die Lagrange-Interpolation bestimmt $P_n$ durch den Ansatz $P_n(x)=\sum _{i=0}^n{y}_i{L}_i(x)$ mit den Lagrange-Polynomen $L_i(x):=\prod _{j\ne i}\frac{x-x_j}{x_i-x_j}$. $L_i$ ist ein Polynom vom Grad $n$ mit $L_i(x_k)={\delta }_{i,k}$.

Newton-Interpolation: Die Newton-Interpolation bestimmt $P_n$ durch den Ansatz
$P_n(x)=\sum _{i=0}^n{a}_i{N}_i(x)$ mit den Newton-Polynomen $N_0(x):=1$ und $N_i(x):=\prod _{j=0}^{i-1}(x-x_j)$ für $i=1,\dots ,n$. Wegen der rekursiven Definition ($N_i(x)=(x-x_{i-1})N_{i-1}(x)$) gilt $N_i(x_k)=0$ für $i>k$. Damit ist das resultierende LGS für die $a_i$ in unterem Dreiecksformat und ein zusätzlicher Samplepunkt $x_{n+1}$ verändert $a_0,\dots ,a_n$ nicht. Außerdem ist dieser Ansatz numerisch stabiler.

Kubische 1D-Interpolation

kubische 1D-Interpolation: Seien $y_0,y_0^{\prime },y_1,y_1^{\prime }\in \mathbb{R}$ gegeben. Gesucht ist ein höchstens kubisches Polynom $f(x)=a{x}^3+b{x}^2+cx+d$ mit $f(i)=y_i$ und $f^{\prime }(i)=y_i^{\prime }$ für $i=0,1$. Man erhält das reguläre LGS $A(a,b,c,d{)}^{\mathrm{T}}=\overrightarrow{v}$ mit $A:=\left(\begin{array}{cccc}0& 0& 0& 1\\ 0& 0& 1& 0\\ 3& 2& 1& 0\\ 1& 1& 1& 1\end{array}\right)$ und $\overrightarrow{v}:=(y_0,y_0^{\prime },y_1^{\prime },y_1{)}^{\mathrm{T}}$. Damit bekommt man $f(x)=(a,b,c,d)(x^3,x^2,x,1{)}^{\mathrm{T}}={\overrightarrow{v}}^{\mathrm{T}}{A}^{-\mathrm{T}}(x^3,x^2,x,1{)}^{\mathrm{T}}$ mit den
kubischen Hermite-Polynomen $\left(\begin{array}{c}{H}_0^3(x)\\ H_1^3(x)\\ H_2^3(x)\\ H_3^3(x)\end{array}\right)=A^{-\mathrm{T}}(x^3,x^2,x,1{)}^{\mathrm{T}}=\left(\begin{array}{c}(2x+1)(1-x{)}^2\\ x(1-x{)}^2\\ -x^2(1-x)\\ (3-2x)x^2\end{array}\right)$

Bikubische Interpolation

bikubische Interpolation: Gegeben seien ${\overrightarrow{f}}_{i,j}=(x_i,y_j{)}^{\mathrm{T}}$ für $x_i,y_j\in \{0,1\}$ und Werte
$z_{i,j},{\partial }_x{z}_{i,j},{\partial }_y{z}_{i,j},{\partial }_x{\partial }_y{z}_{i,j}\in \mathbb{R}$. Gesucht ist bei der bikubischen Interpolation das bikubische Polynom $f(x,y)=\sum _{n=0}^3\sum _{m=0}^3{a}_{n,m}{x}^n{y}^m$ mit $f({\overrightarrow{p}}_{i,j})=z_{i,j}$, ${\partial }_{x}f({\overrightarrow{p}}_{i,j})={\partial }_x{z}_{i,j}$, ${\partial }_{y}f({\overrightarrow{p}}_{i,j})={\partial }_y{z}_{i,j}$ und ${\partial }_x{\partial }_{y}f({\overrightarrow{p}}_{i,j})={\partial }_x{\partial }_y{z}_{i,j}$. Schreibt man $\overrightarrow{\alpha }:=(a_{0,0},a_{1,0},a_{2,0},a_{3,0},a_{0,1},\dots ,a_{3,3}{)}^{\mathrm{T}}$ und
$\overrightarrow{z}=(z_{0,0},z_{1,0},z_{0,1},z_{1,1},{\partial }_x{z}_{0,0},\dots ,{\partial }_x{\partial }_y{z}_{1,1}{)}^{\mathrm{T}}$, dann erhält man ein LGS $B\overrightarrow{\alpha }=\overrightarrow{z}$ mit einer $(16\times 16)$-Matrix $B$. Es gilt dann $\overrightarrow{\alpha }=B^{-1}\overrightarrow{z}$.

bikubische Interpolation mit Hermite-Polynomen: Man kann den Ansatz auch schreiben als $f(x,y)={\overrightarrow{y}}^{\mathrm{T}}A\overrightarrow{x}$ mit $\overrightarrow{x}:=(1,x,x^2,x^3{)}^{\mathrm{T}}=C^{\mathrm{T}}{\overrightarrow{h}}_x$, ${\overrightarrow{h}}_x:=(H_0^3(x),\dots ,H_3^3(x){)}^{\mathrm{T}}$ und einer bestimmten Basiswechsel-Matrix $C=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 1& 2& 3\\ 1& 1& 1& 1\end{array}\right)$, die man aus $C^{-1}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ -3& -2& -1& 3\\ 2& 1& 1& -2\end{array}\right)$ erhält (analog $\overrightarrow{y}=C^{\mathrm{T}}{\overrightarrow{h}}_y$). Damit bekommt man $f(x,y)={\overrightarrow{h}}_y^{\mathrm{T}}CA{C}^{\mathrm{T}}{\overrightarrow{h}}_x$. Man kann zeigen, dass
$CA{C}^{\mathrm{T}}=F:=\left(\begin{array}{cccc}f(0,0)& f_x(0,0)& f_x(1,0)& f(1,0)\\ f_y(0,0)& f_{xy}(0,0)& f_{xy}(1,0)& f_y(1,0)\\ f_y(0,1)& f_{xy}(0,1)& f_{xy}(1,1)& f_y(1,1)\\ f(0,1)& f_x(0,1)& f_x(1,1)& f(1,1)\end{array}\right)$ bzw. $A=C^{-1}F{C}^{-\mathrm{T}}$.

Interpolation auf Dreiecken

baryzentrische Koordinaten: Seien $\overrightarrow{a},\overrightarrow{b},\overrightarrow{c}\in {\mathbb{R}}^2$ nicht kollinear. Dann gibt es für $\overrightarrow{p}\in {\mathbb{R}}^2$ baryzentrische Koordinaten ${\alpha }_1,{\alpha }_2,{\alpha }_3\in \mathbb{R}$ mit $\overrightarrow{p}={\alpha }_1\overrightarrow{a}+{\alpha }_2\overrightarrow{b}+{\alpha }_3\overrightarrow{c}$ und ${\alpha }_1+{\alpha }_2+{\alpha }_3=1$.
Durch Umschreiben $\overrightarrow{p}-\overrightarrow{c}={\alpha }_1(\overrightarrow{a}-\overrightarrow{c})+{\alpha }_2(\overrightarrow{b}-\overrightarrow{c})$ erhält man $({\alpha }_1,{\alpha }_2{)}^{\mathrm{T}}=T^{-1}(\overrightarrow{p}-\overrightarrow{c})$ mit $T:=\left(\begin{array}{cc}\overrightarrow{a}-\overrightarrow{c}& \overrightarrow{b}-\overrightarrow{c}\end{array}\right)$. Durch Ausschreiben der Inverse und Multiplikation bekommt man
${\alpha }_1=\frac{1}{D}det(\left(\begin{array}{cc}\overrightarrow{p}-\overrightarrow{c}& \overrightarrow{b}-\overrightarrow{c}\end{array}\right))$, ${\alpha }_2=\frac{1}{D}det(\left(\begin{array}{cc}\overrightarrow{a}-\overrightarrow{c}& \overrightarrow{p}-\overrightarrow{c}\end{array}\right))$ und ${\alpha }_3=\frac{1}{D}det(\left(\begin{array}{cc}\overrightarrow{b}-\overrightarrow{a}& \overrightarrow{p}-\overrightarrow{a}\end{array}\right))$ mit der Determinate $D:=det(T)$ (doppelter Flächeninhalt des Dreiecks).

geometrische Interpretation: ${\alpha }_1$ ist der Anteil der Fläche des Dreiecks mit der zu $\overrightarrow{a}$ gegenüber liegender Kante $\overrightarrow{b}-\overrightarrow{c}$ als Kante und Ecke $\overrightarrow{p}$ am gesamten Dreieck (falls $\overrightarrow{p}$ im Dreieck liegt). Ob $\overrightarrow{p}$ im Dreieck, auf einer Ecke/Kante oder außerhalb liegt, kann man leicht an den baryzentrischen Koordinaten ablesen.

Transformation zu Einheitsdreieck: Mit der Transformation $x=a_x+(b_x-a_x)\xi +(c_x-a_x)\eta$ und $y=a_y+(b_y-a_y)\xi +(c_y-a_y)\eta$ transformiert man ein Dreieck mit den Ecken $\overrightarrow{a},\overrightarrow{b},\overrightarrow{c}$ auf das Einheitsdreieck. Im Einheitsdreieck gilt ${\alpha }_1=1-\xi -\eta$, ${\alpha }_2=\xi$ und ${\alpha }_3=\eta$.

baryzentrische Interpolation: Seien Werte $f_1,f_2,f_3\in \mathbb{R}$ an den Ecken des Dreiecks gegeben. Dann erhält man baryzentrische Interpolation durch $f(\overrightarrow{p})=\sum _{i=1}^3{\alpha }_i(\overrightarrow{p})f_i$.

quadratische Interpolation: Für quadr. Interpolation im Standarddreieck sei der Ansatz $\overrightarrow{f}(\xi ,\eta )=\sum _{i=1}^6{f}_i{N}_i(\xi ,\eta )$ gegeben mit $N_i({\xi }_j,{\eta }_j)={\delta }_{i,j}$ für $({\xi }_1,{\eta }_1):=(0,0)$, $({\xi }_2,{\eta }_2):=(1,0)$, $({\xi }_3,{\eta }_3):=(0,1)$, $({\xi }_4,{\eta }_4):=(1/2,0)$, $({\xi }_5,{\eta }_5):=(1/2,1/2)$ und $({\xi }_6,{\eta }_6):=(0,1/2)$.
Man erhält $N_1(\xi ,\eta )=(1-\xi -\eta )(1-2\xi -2\eta )$, $N_2(\xi ,\eta )=\xi (2\xi -1)$, $N_3(\xi ,\eta )=\eta (2\eta -1)$, $N_4(\xi ,\eta )=4\xi (1-\xi -\eta )$, $N_5(\xi ,\eta )=4\xi \eta$ und $N_6(\xi ,\eta )=4\eta (1-\xi -\eta )$.

kubische Interpolation: Man kann den Ansatz auch für höhere Ordnungen verallgemeinern. Zum Beispiel wäre der kubische Ansatz
$u(\xi ,\eta )=c_1+c_2\xi +c_3\eta +c_4{\xi }^2+c_5\xi \eta +c_6{\eta }^2+c_7{\xi }^3+c_8{\xi }^2\eta +c_9\xi {\eta }^2+c_{10}{\eta }^3$.

Bikubische Interpolation auf krummlinigen Gittern

krummlinige Gitter: Zur Interpolation auf krummlinigen Gittern muss man die Ableitungen, die im physischen Raum (P-Raum) mit Koord.en $(x,y)$ gegeben sind, umrechnen auf den Berechnungsraum (C-Raum) mit Koord.en $(\xi ,\eta )$. Im C-Raum sollte das Gitter ein uniformes 2D-Rechtecksgitter mit Gitterabständen $\mathrm{\Delta }\xi =1=\mathrm{\Delta }\eta$ sein.

Umrechnung der Ableitungen: Für ${\overrightarrow{\partial }}_{xy}:=({\partial }_x,{\partial }_y,{\partial }_x^2,{\partial }_x{\partial }_y,{\partial }_y^2{)}^{\mathrm{T}}$ und
${\overrightarrow{\partial }}_{\xi \eta }:=({\partial }_{\xi },{\partial }_{\eta },{\partial }_{\xi }^2,{\partial }_{\xi }{\partial }_{\eta },{\partial }_{\eta }^2{)}^{\mathrm{T}}$ erhält man ${\overrightarrow{\partial }}_{xy}=B{\overrightarrow{\partial }}_{\xi \eta }$ mit $B:=\left(\begin{array}{cc}{B}_1& 0\\ B_2& B_3\end{array}\right)$ und
$B_1:=\left(\begin{array}{cc}{\partial }_x\xi & {\partial }_x\eta \\ {\partial }_y\xi & {\partial }_y\eta \end{array}\right)$, $B_2:=\left(\begin{array}{cc}{\partial }_x^2\xi & {\partial }_x^2\eta \\ {\partial }_x{\partial }_y\xi & {\partial }_x{\partial }_y\eta \\ {\partial }_y^2\xi & {\partial }_y^2\eta \end{array}\right)$ sowie $B_3:=\left(\begin{array}{ccc}({\partial }_x\xi {)}^2& 2{\partial }_x\xi \cdot {\partial }_x\eta & ({\partial }_x\eta {)}^2\\ {\partial }_x\xi \cdot {\partial }_y\xi & {\partial }_x\xi \cdot {\partial }_y\eta +{\partial }_x\eta \cdot {\partial }_y\xi & {\partial }_x\eta \cdot {\partial }_y\eta \\ ({\partial }_y\xi {)}^2& 2{\partial }_y\xi \cdot {\partial }_y\eta & ({\partial }_y\eta {)}^2\end{array}\right)$.
Umgekehrt gilt ${\overrightarrow{\partial }}_{\xi \eta }=C{\overrightarrow{\partial }}_{xy}$ mit $C:=B^{-1}=\left(\begin{array}{cc}{C}_1& 0\\ C_2& C_3\end{array}\right)$ und
$C_1:=B_1^{-1}=\left(\begin{array}{cc}{\partial }_{\xi }x& {\partial }_{\xi }y\\ {\partial }_{\eta }x& {\partial }_{\eta }y\end{array}\right)$, $C_2:=-B_3^{-1}{B}_2{B}_1^{-1}=\left(\begin{array}{cc}{\partial }_x^{2}x& {\partial }_{\xi }^{2}y\\ {\partial }_{\xi }{\partial }_{\eta }x& {\partial }_{\xi }{\partial }_{\eta }y\\ {\partial }_{\eta }^{2}x& {\partial }_{\eta }^{2}y\end{array}\right)$ sowie
$C_3:=B_3^{-1}=\left(\begin{array}{ccc}({\partial }_{\xi }x{)}^2& 2{\partial }_{\xi }x\cdot {\partial }_{\xi }y& ({\partial }_{\xi }y{)}^2\\ {\partial }_{\xi }x\cdot {\partial }_{\eta }x& {\partial }_{\xi }x\cdot {\partial }_{\eta }y+{\partial }_{\xi }y\cdot {\partial }_{\eta }x& {\partial }_{\xi }y\cdot {\partial }_{\eta }y\\ ({\partial }_{\eta }x{)}^2& 2{\partial }_{\eta }x\cdot {\partial }_{\eta }y& ({\partial }_{\eta }y{)}^2\end{array}\right)$.