Aufgabe 1 - Bildrekonstruktion

Setup:
Wir bezeichnen (als erstes, um Schreibarbeit später zu sparen) das Gebiet, auf dem das Bild der Größe \(w \times h\) lebt als \(\Omega = [0, w] \times [0, h]\).
Gegeben ist eine Bildfunktion ("\(d\) wie Daten") \[d: \Omega \to \mathbb{R}.\]
Gesucht ist eine Rekonstruktion ("\(f\) wie so-heißt-immer-unsere-Unbekannte"): \[f: \Omega \to \mathbb{R}.\]
Aufstellen eines Energiefunktionals:
Wir nehmen zur Vereinfachung an, dass unser Bild schwarz/weiß ist (Farbe hinzuzufügen ist nicht so schwer).
Das Bild wurde mit unabhängigen gleichverteilten Gauß'schem Rauschen verunstaltet (anders gesagt, die Kamera des Mobiltelefons hatte 20 Megapixel, aber es war viel zu dunkel draußen).
Um den Fehler zu bewerten formulieren wir eine least-squares Energiefunktion, die Abweichungen der Lösung bestraft \[ \begin{equation} \begin{split} E(f) & = || f - d ||^2 \\ & = \int_\Omega ( f(x) - d(x))^2 \,\, dx \\ & \to \text{min}. \\ \end{split} \end{equation} \] Diese Funktion modelliert zwar das Gauß'sche Rauschen adäquat, liefert aber an sich keine neuen Erkenntnisse.

Daher fügen wir einen Prior hinzu:
(Für Hörer der Machine-Learning-Vorlesung: statistisch formuliert nehmen wir eine Gaußverteilung auf den Gradienten als Prior Distribution an). Wir nehmen also zusätzlich an, dass starke Gradienten selten sind und ebenfalls bestraft werden müssen: \[ \begin{equation} \begin{split} E_{besser} (f) & = || f - d ||^2 + \lambda || \nabla f ||^2\\ & = \int_\Omega ( f(x) - d(x))^2 \,\, dx + \lambda \int_\Omega || \nabla f(x) ||^2 \,\, dx\\ & \to \text{min}. \\ \end{split} \end{equation} \] Der Parameter \(\lambda\) stellt dabei ein, wie stark das Bild geglättet werden soll.
Nun zur implementierbaren Lösung dieser Energie:
Wir diskretisieren diese Gleichung, indem wir \(f\) durch den Vektor aller Pixel ersetzen. Dies entspricht also einen Vektor mit \(w \times h\) Einträgen (jeder Pixel bekommt einen eigenen). Das heißt, wir nutzen die folgende Abtastung: \[ f \left( \frac{i}{w}, \frac{j}{h} \right) \approx x[i + j \times w]. \]
Der Datenterm ergibt sich dann als \[ \sum_{i=1}^{w}\sum_{j=1}^{h} \left( \, x[i + j \times w] - d[i + j \times w] \, \right) ^2. \] Den Term \(\frac{1}{wh}\) für das differentielle Flächenelement im Integral lassen wir der Einfachheit halber hier weg, und bei den Gradienten auch (in der Minimierung hebt sich das dann auf).
Weiterhin stellen wir den Gradientenvektor durch finite Differenzen dar: \[ \nabla f \left( \frac{i}{w}, \frac{j}{h} \right) = \left( \begin{array}{cc} \partial_x f( \frac{i}{w}, \frac{j}{h}) \\ \partial_y f( \frac{i}{w}, \frac{j}{h}) \\ \end{array} \right) \approx \left( \begin{array}{cc} x[(i + 1) + j \cdot w] - x[i +j \cdot w] \\ x[i + (j + 1) \cdot w] - x[i + j \cdot w] \\ \end{array} \right). \]
Und die Summe aller Quadrate dieser Gradienten approximiert das zweite Integral in \(E_{besser}\).
Den üblichen Nenner der finiten Differenzen kann man der Einfachheit halber auch weglassen; dies verändert nur den Wertebereich von \(\lambda\).
Damit lässt sich eines der üblichen multivariaten quadratischen Polynome für die Zielfunktion \(E_{besser}\) aufstellen, und dessen kritischer Punkt durch setzten der Ableitung \((= 0)\) finden.