bildbasierte algorithmen zur metall-artefakt reduktion in ... dreier.pdf · bildbasierte...
Post on 13-Aug-2019
221 Views
Preview:
TRANSCRIPT
Bildbasierte Algorithmen zur
Metall-Artefakt Reduktion in
der Computertomographie
Bachelorarbeitzur Erlangung des akademischen Grades
Bachelor of Science
Westfalische Wilhelms-Universitat Munster
Fachbereich Mathematik und Informatik
Institut fur Numerische und Angewandte Mathematik
Betreuung:
Dr. Frank Wubbeling
Prof. Dr. Martin Burger
Eingereicht von:
Nils-Arne DreierEmail: n.dreier@uni-muenster.de
Munster, September 2014
i
Inhaltsverzeichnis
1 Einleitung 1
2 Grundlagen der Computertomographie und das Problem der Metall-Arte-
fakte 3
2.1 Aufbau und Funktionsweise eines Computertomographen . . . . . . . . 3
2.2 Das Problem der Metall-Artefakte . . . . . . . . . . . . . . . . . . . . . 5
3 Mathematische Modellierung und klassische Algorithmen zur Rekonstruk-
tion 8
3.1 Radon-Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.2 Modellierung und Diskretisierung . . . . . . . . . . . . . . . . . . . . . 11
3.3 Gefilterte Ruckprojektion . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.4 Algebraic-Reconstruction-Technique . . . . . . . . . . . . . . . . . . . . 13
4 Regularisierungsmethoden 16
4.1 Tikhonov-Regularisierung . . . . . . . . . . . . . . . . . . . . . . . . . 16
4.2 Total-Variation-Regularisierung . . . . . . . . . . . . . . . . . . . . . . 21
5 Anpassung der Algorithmen auf das Problem der Metall-Artefakt Reduk-
tion 35
5.1 Lineare Interpolation und gefilterte Ruckprojektion . . . . . . . . . . . 35
5.2 Algebraic-Reconstruction-Technique . . . . . . . . . . . . . . . . . . . . 36
5.3 Tikhonov-Regularisierung . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.4 Total-Variation-Regularisierung . . . . . . . . . . . . . . . . . . . . . . 38
6 Vergleich der Methoden 40
7 Fazit und Ausblick 56
Abbildungsverzeichnis 58
Tabellenverzeichnis 59
Literaturverzeichnis 60
1
1 Einleitung
Die Computertomographie ermoglicht es, Schnittbilder eines menschlichen Korpers zu
erstellen und auf diese Weise Tumore, Knochen, Gewebe oder Blutgefaße sichtbar zu
machen sowie dessen Position und Form zu erkennen. Neben der Medizin gibt es zahl-
reiche weitere Einsatzgebiete der Computertomographie, wie etwa die Untersuchung
von Baumen oder der industriellen Materialuberprufung. Die Mathematik spielt bei
der Erstellung des Bildes eine wesentliche Rolle. Es gibt Situationen, in denen das re-
sultierende Bild Artefakte aufweist, das heißt, dass es im Bild Bereiche gibt, die gestort
sind - also nicht der Wahrheit entsprechen. Es besteht die Gefahr, dass aus dem Bild
keine, beziehungsweise eine falsche medizinische Diagnose gestellt wird.
Grunde dafur konnen Metallobjekte, wie etwa Zahnimplantate, Herzschrittmacher,
Prothesen oder Schrauben sein, die sich im Bildbereich befinden. Rontgenstrahlung
wird von Metall weitaus starker abgeschwacht als von Gewebe oder Knochen, sodass
die Intensitat der Rontgenstrahlen nach dem Durchdringen des Objekts so klein ist,
dass eine prazise Messung nicht mehr moglich ist.
Ziel dieser Arbeit ist es, verschiedene bildbasierte Methoden zur Rekonstruktion von
Daten eines Computertomographen zu erarbeiten und diese so zu verandern, dass
Metall-Artefakte reduziert werden.”Bildbasiert“ heißt dabei, dass die Reduktion der
Metall-Artefakte wahrend der Rekonstruktion auf Grundlage des Bildes geschieht. Im
Gegensatz dazu gibt es datenbasierte Algorithmen, welche die gemessenen Daten auf-
bereiten, um die Artefakte zu reduzieren.
Ein aktuelles Thema in der mathematischen Bildbearbeitung ist die Methode der Total-
Variation. Diese lasst sich dazu verwenden, Metall-Artefakte zu reduzieren, wie bei-
spielsweise von Schiffer und Bredies [SB14] aufgezeigt wird. Ein großer Teil dieser Ar-
beit wird sich damit beschaftigen, diese Methode zu untersuchen und auf das Problem
der Metall-Artefakt Reduktion anzuwenden.
Im Rahmen dieser Arbeit werden wir in Kapitel 2 die Grundlagen der Computerto-
mographie untersuchen und uns die Eigenarten der Metall-Artefakte genauer ansehen.
1 Einleitung 2
Im Anschluss daran beschaftigen wir uns mit der mathematischen Modellierung sowie
mit zwei klassischen Rekonstruktionsverfahren. Im vierten Kapitel werden wir uns wei-
tere Rekonstruktionsverfahren erarbeiten, welche das resultierende Bild glatten. Wir
erhoffen uns hiervon, die Metall-Artefakte reduzieren zu konnen. Das anschließende
Kapitel 5 untersucht, inwiefern sich die Verfahren auf das Problem der Metall-Artefakt
Reduktion anpassen lassen. Abschließend wollen wir in Kapitel 6 einen Vergleich der
Methoden aufstellen. Dabei erzeugen wir kunstlich Daten von relevanten Situationen
und vergleichen die Ergebnisse der Verfahren. Am Ende der Arbeit steht ein Fazit, wel-
ches auf die Ergebnisse des Vergleichs eingeht und die Vor- und Nachteile der Methoden
herausstellt.
3
2 Grundlagen der
Computertomographie und das
Problem der Metall-Artefakte
In diesem einfuhrenden Kapitel soll die Funktionsweise eines Computertomographen
bundig erklart und das Problem der Metall-Artefakte dargelegt werden. Außerdem wer-
den erste Modellannahmen getroffen, bevor wir im nachsten Kapitel das mathematische
Modell entwickeln.
2.1 Aufbau und Funktionsweise eines
Computertomographen
Ein Computertomograph dient der Aufnahme von Schnittbildern eines Objekts. Dazu
wird die Abschwachung von Rontgenstrahlen gemessen, die das Objekt durchleuch-
ten. Aus den gemessenen Daten kann mit einem Computer schließlich ein Schnittbild
errechnet werden. Ublicherweise werden mehrere Schnittbilder aufgenommen und zu
einem dreidimensionalen Bild zusammengesetzt. Im Gegensatz zu einer Rontgenauf-
nahme, sind die Bilder eines Computertomographen uberlagerungsfrei, das heißt, dass
sich Bereiche im Objekt auf dem Bild nicht uberlagern.
Ein Computertomograph besteht aus einer Rontgenquelle und einem Rontgendetektor,
welche sich gegenuberliegend um das zu durchleuchtende Objekt drehen. Dabei wer-
den Projektionen des Objekts erstellt, indem der Rontgendetektor die Intensitat der
Rontgenstrahlen misst, welche das Objekt durchdringen. Es gibt verschiedene Vor-
gehensweisen, die Projektionen zu erstellen – in dieser Arbeit mochten wir nur die
Parallelstrahlgeometrie (Parallel-Scanning) untersuchen. Bei anderen Geometrien sind
die Methoden ebenfalls anwendbar und wir wurden ahnliche Ergebnisse erwarten. Bei
der Parallelstrahlgeometrie verlaufen die Rontgenstrahlen, aus denen die Projektionen
2 Grundlagen der Computertomographie und das Problem der Metall-Artefakte 4
erstellt werden, parallel zueinander. Dies entspricht einer unendlich weit entfernten
Rontgenquelle, oder die Rontgenquelle und der Detektor mussten transversal in jedem
Winkelschritt verschoben werden. Dieses Vorgehen wurde bei den ersten Computerto-
mographen so praktiziert. Heutzutage wird mit mehreren Detektoren gearbeitet, welche
die Strahlung von einer Rontgenquelle messen, das heißt, dass die Strahlen facherformig
verlaufen.
Die gemessenen Daten konnen in einem sogenannten Sinogramm dargestellt werden:
Auf der X-Achse werden die Winkel und auf der Y-Achse der Radialabstand des ge-
messenen Rontgenstrahls aufgetragen (vgl. Abbildung 2.1).
(a) Shepp-Logan Phantom (b) Sinogramm des Shepp-Logan Phantoms
Abbildung 2.1: Shepp-Logan Phantom und dessen Sinogramm
Die Abschwachung von Rontgenstrahlen durch eine Substanz mit dem Abschwachungs-
koeffizienten µ wird beschrieben durch
I(s) = I(0) exp
− s∫0
µ(η) dη
. (2.1)
Dabei beschreibt I(0) die Anfangsintensitat und I(s) die am Detektor gemessene In-
tensitat. Das Ziel der Computertomographie ist es, den Abschwachungskoeffizienten
µ als Bild darzustellen. Diese Problemstellung wird auch als Inverses Problem be-
zeichnet. Es soll also von einer Wirkung (Abschwachung der Rontgenstrahlen bei dem
Durchleuchten des Objekts) auf eine Ursache (Absorptionskoeffizienten innerhalb des
Objekts) geschlossen werden. Fur die Berechnung wird hingegen das Integral uber die
2 Grundlagen der Computertomographie und das Problem der Metall-Artefakte 5
Absorptionskoeffizienten verwendet, welches durch Umstellen der Gleichung (2.1) aus
den gemessenen Daten errechnet werden kann.
s∫0
µ(η) dη = − ln(I(s)
I(0))
Dies ermoglicht uns die Modellierung mithilfe der Radon-Transformation (vgl. Kapitel
3).
2.2 Das Problem der Metall-Artefakte
Als Artefakt bezeichnet man einen Fehler im rekonstruierten Bild. Artefakte konnen zur
Folge haben, dass das Bild seine Aussagekraft verliert oder im medizinischen Bereich
falsche Diagnosen gestellt werden. Neben Metall-Artefakten gibt es viele weitere Arten
von Artefakten. Metall-Artefakte kennzeichnen sich durch helle und dunkle Streifen,
die vom Metallobjekt ausgehen, aus (vgl. Abbildung 2.2).
Metallobjekte, wie beispielsweise Zahnimplantate, Schrauben oder Herzschrittmacher,
haben einen deutlich hoheren Abschwachungskoeffizienten als ihn etwa Wasser, Kno-
chen oder Gewebe besitzen. Daher kann es ab einer bestimmten Metalldicke zu einer
sogenannten Totalabsorption kommen. Das heißt, dass die Strahlungsintensitat, die am
Detektor gemessen wird, unterhalb der Rauschschwelle liegt. Unter Betrachtung von
Gleichung (2.1) bedeutet das
I(0) exp
− s∫0
µ(η) dη
= 0
⇒s∫
0
µ(η) dη =∞.
Aus den Daten kann also nur die Information geschlossen werden, dass auf der Gera-
den ein Material mit sehr hohem Absorptionskoeffizienten liegt. Die Daten enthalten
aber keine Information uber die Absorptionskoeffizienten der Materialien, die vor so-
wie hinter dem Metallobjekt liegen. Insbesondere werden die Daten inkonsistent ge-
genuber Daten aus anderen Projektionen – dies stellt die Hauptursache fur Metall-
Artefakte dar. Andere Grunde fur die Artefakte sind Compton-Streuung, Rauschen
2 Grundlagen der Computertomographie und das Problem der Metall-Artefakte 6
oder Detektor-Unterabtastung und auch Bewegungen des Objekts wahrend der Auf-
nahme wirken sich starker aus, wenn Metall anwesend ist (vgl. [DMND+98]). Fur die
Rekonstruktion bedeutet das, dass die Methode diese Daten anders behandeln muss.
Die Standardmethode, die gefilterte Ruckprojektion, kann das nicht ohne Weiteres. Da-
her versucht man stattdessen die Daten, die durch das Metall verloren gegangen sind,
durch moglichst passende Daten zu ersetzen, um die Methode trotzdem anwenden zu
konnen. Dies geschieht etwa durch Interpolation der gestorten Daten in den einzel-
nen Projektionen oder durch Inpainting im Sinogramm. Solche Methoden nennen wir
datenbasiert. In Abschnitt 5.1 werden wir den einfachsten Fall, den der linearen In-
terpolation der Projektionen, ansehen und in Kapitel 6 mit den anderen bildbasierten
Methoden vergleichen.
2 Grundlagen der Computertomographie und das Problem der Metall-Artefakte 7
(a) Als Vergleich: Rekonstruktion vonungestorten Daten mit gefilterterRuckprojektion.
(b) Rekonstruktion von gestorten Daten mitgefilterter Ruckprojektion.
(c) Rekonstruktion mittels gefilterterRuckprojektion von gestorten Daten,die mit linearer Interpolation aufbereitetwurden.
(d) Rekonstruktion durch ein iteratives Ver-fahren, welches die gestorten Daten igno-riert.
Abbildung 2.2: Metall-Artefakte im Shepp-Logan Phantom. 180 Projektionen.
8
3 Mathematische Modellierung und
klassische Algorithmen zur
Rekonstruktion
3.1 Radon-Transformation
In diesem Abschnitt mochten wir uns kurz mit der Radon-Transformation beschaftigen.
Sie ermoglicht uns die mathematische Modellierung der Computertomographie. Wir
orientieren uns dabei an dem Buch von Natterer und Wubbeling [NW01, Kapitel 2.1]
sowie an dem Skript der Vorlesung uber Inverse Probleme von Wubbeling und Pietsch-
mann [PW12].
Einsteigen wollen wir mit der Definition eines Funktionen-Raumes, welche wir fur die
Definition der Radon-Transformation benotigen.
Definition 3.1.1 (Schwartzraum). Es sei Ω ⊂ Rd. Eine Funktion f ∈ C∞0 (Ω) heißt
Schwartzfunktion, falls gilt
supx∈Ω
∣∣∣∣xl · ∂k
∂xkf(x)
∣∣∣∣ <∞ ∀k, l ∈ N0.
Der Raum aller Schwartzfunktionen auf Ω heißt Schwartzraum und wird mit S(Ω)bezeichnet.
Definition 3.1.2 (Radon-Transformation). Sei C := S1×R der unendlich lange Ein-
heitszylinder im R3 und f ∈ S(R2) eine Schwartzfunktion auf R2. Dann ist die Radon-
Transformation Rf ∈ S(C) definiert durch
Rf(θ, s) :=∫
x·θ=s
f(x) dσ(x). (3.1)
Bemerkung 3.1.3.
• R ist linear.
3 Mathematische Modellierung und klassische Algorithmen zur Rekonstruktion 9
• Die Daten eines Computertomographen entsprechen einer Funktion g ∈ S(C).
Fur das zu errechnende Bild f ∈ S(R2) gilt also Rf = g.
• Das dσ(x) in der Definition soll andeuten, dass wir das Linienintegral meinen.
Wir werden das σ aber auch weglassen.
Als Nachstes wollen wir die L2-Adjungierte der Radon-Transformation einfuhren, wel-
che spater fur einige iterative Verfahren benotigt wird.
Lemma 3.1.4 (Die L2-Adjungierte der Radon-Transformation). Die ungefilterte Ruck-
projektion R∗ : S(C)→ S(R2)
R∗g(x) :=
∫S1
g(θ, x · θ)dθ (3.2)
ist L2-adjungiert zu der Radon-Transformation.
Beweis.
(Rf, g)L2(C) =
∫S1
∫R
Rf(θ, s)g(θ, s) ds dθ
=
∫S1
∫R
∫x·θ=s
f(x) dx g(θ, s) ds dθ
=
∫S1
∫R2
f(x) g(θ, x · θ) dx dθ
=
∫R2
f(x)
∫S1
g(θ, x · θ) dθ dx
=
∫R2
f(x)R∗g dx = (f,R∗g)L2(R2)
Der Name”Ruckprojektion“ lasst schon ahnen, was der Operator macht: Er weist jedem
Punkt die Projektionswerte zu, die durch ihn entstanden sind und integriert diese auf.
Das Adjektiv”ungefiltert“ wird verwendet, um es von der gefilterten Ruckprojektion,
welcher wir uns im Abschnitt 3.3 zuwenden wollen, abzugrenzen. Wenn klar ist, von
welcher Ruckprojektion die Rede ist, werden wir auch nur”Ruckprojektion“ sagen.
3 Mathematische Modellierung und klassische Algorithmen zur Rekonstruktion 10
An dieser Stelle wollen wir noch das Fourier-Slice-Theorem, welches eine Grundlage
fur die gefilterte Ruckprojektion bildet, anfuhren. Dafur wollen wir zunachst festlegen
was wir unter der Fourier-Transformation verstehen.
Definition 3.1.5 (Fourier-Transformation). Fur eine Funktion f ∈ S(Rn) definieren
wir die Fourier-Transformation f ∈ S(Rn) von f durch
f(ξ) := (2π)−n2
∫Rn
f(x)e−iξ·x dx.
Analog dazu definieren wir die inverse Fourier-Transformation f ∈ S(Rn) von f durch
f(x) := (2π)−n2
∫Rn
f(ξ)eiξ·x dξ.
Satz 3.1.6 (Fourier-Slice-Theorem (vgl. [NW01, Theorem 2.1])). Sei f ∈ S(R2), θ ∈S1 und σ ∈ R. Dann gilt
(Rf)(θ, σ) =√2πf(σθ).
Auf der linken Seite ist dabei a priori nicht klar, wie die Fourier-Transformation einer
Funktion auf C zu verstehen ist. Wir wollen vereinbaren, dass wir in so einem Fall die
Fourier-Transformation nur in der zweiten Variable anwenden, also Rf(θ, ·)(σ).
Beweis.
Rf(θ, σ) = 1√2π
∫R
Rf(θ, s)e−isσ ds
=1√2π
∫R
∫x·θ=s
f(x) dx e−isσ ds
=1√2π
∫R2
f(x)e−i(x·θ)σ dx
=√2πf(σθ).
Mit der Invertierbarkeit der Fourier-Transformation folgt aus diesem Satz insbesondere
die Invertierbarkeit der Radon-Transformation.
3 Mathematische Modellierung und klassische Algorithmen zur Rekonstruktion 11
3.2 Modellierung und Diskretisierung
Nachdem wir im vorherigen Abschnitt einige mathematische Grundlagen erarbeitet
haben, wollen wir in diesem Abschnitt darlegen, wie wir diese nutzen konnen. Dazu
wollen wir zunachst zwischen zwei Bildmodellen unterscheiden:
Im kontinuierlichen Bildmodell wird ein Bild durch eine Funktion f ∈ S(Ω) auf demBildgebiet Ω ⊂ R2 reprasentiert. Dabei ordnet f jedem Bildpunkt p ∈ Ω einen Grau-
wert f(p) zu. Um mit solchen Funktionen auf einem Computer arbeiten zu konnen,
muss man sie in ein diskretes Bildmodell uberfuhren. Dafur wahlen wir ein endliches,
linear unabhangiges System von Ansatzfunktionen Φ = φ1, . . . , φn ⊂ S(Ω). Der vonΦ aufgespannte Unterraum ist unser diskreter Bildraum. Nach den Erkenntnissen der
linearen Algebra lasst sich jedes Bild im diskreten Bildraum durch einen Koordinaten-
vektor zur Basis Φ darstellen. Diesen konnen wir, weil er endlich viele Eintrage hat,
auf einem Computer speichern und verarbeiten.
Fur die Wahl der Ansatzfunktionen teilen wir das Bildgebiet Ω durch ein Gitter Ω =
x1, . . . , xn × y1, . . . , ym auf und wahlen fur jeden Gitterpunkt (xi, yj) ∈ Ω eine
Ansatzfunktion φi,j ∈ S(Ω) mit
φi,j(xk, yl) =
1 , falls i = k und j = l
0 , falls i = k oder j = l,
etwa abgeschnittene Gaußglocken, Hutchenfunktionen oder Polynome. So konnen wir
eine Approximation an ein kontinuierliches Bild erstellen, indem wir sie an den Gitter-
punkten abtasten und erhalten so den Koordinatenvektor.
In unseren Implementationen der Algorithmen werden wir stets auf die in MAT-
LAB bereitgestellten Funktionen radon und iradon zuruckgreifen, um eine Radon-
Transformation oder eine ungefilterte Ruckprojektion zu berechnen. Fur eine detail-
liertere Auseinandersetzung mit kontinuierlichen und diskreten Bildern sei an dieser
Stelle auf das Buch von Bredies und Lorenz verwiesen [BL10, Kapitel 3.1].
Wir wollen davon ausgehen, dass uns bekannt ist wo sich das Metallobjekt befindet, um
unsere Algorithmen dem Problem anpassen zu konnen. Wir setzen also voraus, dass uns
eine Menge Ω0 ⊂ Ω bekannt ist, die das Metallobjekt reprasentiert. Außerdem konnen
wir damit ausrechnen, welche Daten von dem Metallobjekt gestort werden, indem wir
3 Mathematische Modellierung und klassische Algorithmen zur Rekonstruktion 12
die Menge Λ ⊂ C betrachten.
Λ := (θ, s) ∈ C∣∣RIΩ0(θ, s) = 0 mit IΩ0(x) =
1 , falls x ∈ Ω0
0 , sonst
In der Praxis kann die Region Ω0, in der sich das Metall befindet, bestimmt werden,
indem auf die Daten g ∈ S(C) die ungefilterte Ruckprojektion angewandt wird. Alle
Punkte, fur die R∗g einen Grenzwert c ∈ R uberschreiten, werden zu Ω0 gezahlt (vgl.
[PKBK09, Abschnitt 2.2]).
Ω0 =x ∈ Ω
∣∣R∗g(x) > c
3.3 Gefilterte Ruckprojektion
An dieser Stelle wollen wir uns der Standardmethode zur Rekonstruktion widmen.
Dazu fuhren wir mithilfe der Fourier-Transformation das Riesz-Potential ein.
Definition 3.3.1 (Riesz-Potential). Sei f ∈ S(Rn) eine Schwartzfunktion und α ∈ R.Dann ist das Riesz-Potential Iαf von f definiert durch
Iαf(ξ) = |ξ|−αf(ξ).
Bemerkung 3.3.2. Aus der Definition folgt, dass IβIαf = Iβ+αf und I0f = f
gilt. Außerdem definieren wir das Riesz-Potential fur Funktionen in S(C) analog zur
Fourier-Transformation, nur auf dem zweiten Argument.
Satz 3.3.3 (vgl. [NW01, Theorem 2.5]). Es sei f ∈ S(R2) und g = Rf . Dann gilt
f =1
2(2π)−1I−αR∗Iα−1g.
Beweis. Wir betrachten das Riesz-Potential von f mithilfe der inversen Fourier-Trans-
formation und wenden die Transformation auf Polarkoordinaten an.
Iαf(x) = (2π)−1
∫R2
eiξ·x|ξ|−αf(ξ) dξ
=1
2(2π)−1
∫S1
∫R
eirθ·x|r|1−αf(rθ) dr dθ
3 Mathematische Modellierung und klassische Algorithmen zur Rekonstruktion 13
Mit dem Fourier-Slice-Theorem 3.1.6 folgt nun
=1
2(2π)−1
∫S1
∫R
eirθ·x(2π)−12 |r|1−αRf(θ, r) dr dθ
=1
2(2π)−1
∫S1
∫R
eirθ·x(2π)−12 (Iα−1Rf)(θ, r) dr dθ
=1
2(2π)−1
∫S1
Iα−1Rf(θ, θ · x) dθ
=1
2(2π)−1R∗Iα−1g.
Mit der Anwendung von I−α folgt f = 12(2π)−1I−αR∗Iα−1g.
Wir erhalten hier eine ganze Familie von Rekonstruktionsformeln. Die Anwendung der
Formel im Fall α = 0 wird gefilterte Ruckprojektion genannt:
f =1
2(2π)−1R∗I−1g
3.4 Algebraic-Reconstruction-Technique
Als Algebraic-Reconstruction-Technique (ART) bezeichnet man das Losen des end-
lich-dimensionalen Gleichungssystems Rf = g anhand der Kaczmarz-Methode, wel-
che ein iteratives Losungsverfahren zur Losung linearer Gleichungssysteme darstellt.
Fur unsere Gleichung bedeutet das, dass die Projektionen sukzessive auf das Bild
zuruckprojiziert werden, sodass alle Gleichungen zu der Projektion erfullt werden. Wie
der Name schon vermuten lasst, handelt es sich dabei um eine algebraische Methode,
das heißt, dass fur die Entwicklung der Methode schon von Daten in diskreter Form
ausgegangen wird. Als Grundlage fur dieses Kapitel dient das Buch von Natterer und
Wubbeling [NW01, Kapitel 5.3.1].
Es sei H = S(R2) und Θ = θ1, . . . , θp ⊂ S1 die Menge der Projektionswinkel. Dann
betrachten wir die Radon-Transformation in eine feste Richtung θ ∈ Θ
Rθf(s) := Rf(θ, s) =∫
x·θ=s
f(x) dx.
3 Mathematische Modellierung und klassische Algorithmen zur Rekonstruktion 14
Sei eine rechte Seite g ∈ S(C) gegeben. Von unserer Losung f ∈ S(R2) erwarten wir
dann, dass fur alle θ ∈ Θ gilt
Rθf = g(θ, ·).
Fur θ ∈ Θ konnen wir die affin-linearen Unterraume Hθ := f ∈ S(R2)∣∣Rθf = g(θ, ·)
und die orthogonalen Projektionen Pθ von S(R2) auf Hθ betrachten. Fur eine Losung
f gilt dann
f ∈∩θ∈Θ
Hθ.
Wir wollen die Projektionen sukzessive anwenden. Um die Ungleichheiten bezuglich
der Reihenfolge der Anwendungen auszugleichen, fuhren wir einen Relaxionsparameter
ω ∈ (0, 2) ein und setzen
P ωθ = (1− ω)I + ωPθ.
Wir definieren dann P ω := P ωθp · · · P ω
θ1und die Kaczmarz-Methode durch
fk = P ωfk−1 mit beliebigem f 0 ∈ S(R2).
Fur die Berechnung mussen wir die P ωθ angeben. Die orthogonalen Projektionen Pθ
sind gegeben durch
Pθf = f +R∗θ(RθR∗
θ)−1(g(θ, ·)−Rθf).
Es stellt sich die Frage, wie R∗θ aussieht und ob R∗
θ berechenbar ist. Wie auch bei
der Radon-Transformation stellt sich heraus, dass die Adjungierte eine sehr einfache
Darstellung besitzt. Es gilt (R∗θφ)(x) = φ(x · θ), denn es gilt
(Rθf, φ)R =
∫R
∫x·θ=s
f(x) dxφ(s) ds
=
∫R2
f(x)φ(θ · x) dx = (f,R∗φ)R2 .
Nehmen wir an, dass Ω ein Kreis mit Radius ρ ist, so gilt
(RθR∗θ)φ(s) =
∫x·θ=s
R∗θφ(x) dx =
∫x·θ=s
φ(s) dx = 2√
ρ2 − s2φ(s).
3 Mathematische Modellierung und klassische Algorithmen zur Rekonstruktion 15
Dabei ist 2√
ρ2 − s2 die Lange der Geraden x · θ = s durch den Kreis Ω. Und somit
ergibt sich
Pθf(x) =f(x) +R∗θhf (x)
hf (s) :=1
2√
ρ2 − s2(g(θ, s)−Rθf(s)).
16
4 Regularisierungsmethoden
Regularisierungsmethoden stellen ein Werkzeug dar, um mit schlecht gestellten inversen
Problemen umgehen zu konnen.”Schlecht gestellt“ meint dabei, dass die Existenz oder
Eindeutigkeit der Losung nicht gegeben ist beziehungsweise die Losung nicht stetig von
den Anfangsdaten abhangt. Dabei lost die regularisierte Losung nicht das Problem, son-
dern ist eine Approximation der Losung, welche allerdings weniger storungsanfallig ist.
In der Bildverarbeitung bedeutet Regularisierung meist, dass das Bild geglattet wird.
Die Standardmethode zur Regularisierung eines linearen Problems ist die Tikhonov-
Regularisierung, welche wir uns im ersten Abschnitt diesen Kapitels ansehen wollen.
Im Anschluss daran untersuchen wir die Total-Variation-Regularisierung.
4.1 Tikhonov-Regularisierung
Diesem Abschnitt liegt das System Rf = g zugrunde, welches wir als lineares Glei-
chungssystem losen und mithilfe der Tikhonov-Regularisierung regularisieren wollen.
Wir wollen unsere Ergebnisse aber in allgemeiner Form fur ein System Ax = b formulie-
ren. Wir orientieren uns dabei an Muller und Siltanen [MS12] sowie an der Masterarbeit
von Koskela [Kos12].
Definition 4.1.1 (Tikhonov-Regularisierung). Sei A : Rn → Rm eine lineare Abbil-
dung, b ∈ Rm eine rechte Seite und L : Rn → Rp eine lineare Abbildung. Dann heißt
ein Minimierer x ∈ Rn von
Φα(x) = ∥Ax−m∥22 + α∥Lx∥22 (4.1)
Tikhonov-Regularisierung von Ax = b bezuglich des Regularisierungsparameters α.
Eine Tikhonov-Regularisierung lost also nicht die Gleichung Ax = b, sondern stellt
vielmehr einen Kompromiss zwischen kleinem Residuum ∥Ax − b∥ und kleiner Norm
∥Lx∥ dar. Dabei ist das L noch frei wahlbar.
4 Regularisierungsmethoden 17
Folgendes Lemma gibt uns ein notwendiges Kriterium fur eine Tikhonov-Regularisie-
rung, welches wir im Anschluss weiterentwickeln (vgl. [MS12, Kapitel 5]).
Lemma 4.1.2 (Notwendiges Kriterium fur eine Tikhonov-Regularisierung). Seien alle
Bezeichnungen wie in Def. 4.1.1. Dann gilt fur eine Tikhonov-Regularisierung x ∈ Rn
(ATA+ αLTL)x = AT b. (4.2)
Beweis. Da x ein Minimierer von Φα ist, gilt
∇Φα(x) = 0.
Also gilt fur beliebiges w ∈ Rn
0 =∂
∂t
(∥A(x+ tw)− b∥2 + α∥L(x+ tw)∥2
)∣∣∣∣t=0
.
Wir berechnen die Summanden einzeln:
∂
∂t
(∥A(x+ tw)− b∥2
)∣∣∣∣t=0
=∂
∂t⟨Ax+ Atw − b, Ax+ Atw − b⟩
∣∣∣∣t=0
=∂
∂t
(∥Ax∥2 + 2t ⟨Ax,Aw⟩+ t2∥Aw∥2 − 2t ⟨b, Aw⟩ − 2 ⟨Ax, b⟩+ ∥b∥2
)∣∣∣∣t=0
=2 ⟨Ax,Aw⟩ − 2 ⟨b, Aw⟩
und
∂
∂tα ⟨L(x+ tw), L(x+ tw)⟩
∣∣∣∣t=0
= α∂
∂t
(∥Lx∥2 + 2t ⟨Lx, Lw⟩+ t2∥Lw∥2
)∣∣∣∣t=0
=α(2 ⟨Lx, Lw⟩+ 2t∥Lw∥2
)∣∣t=0
= 2α ⟨Lx, Lw⟩ .
Zusammen gilt ⟨Ax− b, Aw⟩+ α ⟨Lx, Lw⟩ = 0, also
0 =⟨ATAx− AT b, w
⟩+ α
⟨LTLx, w
⟩=⟨(ATA+ αLTL)x− AT b, w
⟩.
Da w ∈ Rn beliebig war, muss gelten
(ATA+ αLTL)x− AT b = 0.
Lemma 4.1.3. Fur α > 0 ist R∗R+αLTL symmetrisch positiv definit, also insbeson-
dere invertierbar.
4 Regularisierungsmethoden 18
Beweis. Die Symmetrie folgt aus der Darstellung, da R∗R und αLTL symmetrisch
sind. Fur x = 0 gilt
⟨(R∗R+ αLTL)x, x
⟩= ⟨R∗Rx, x⟩+
⟨αLTLx, x
⟩= ⟨Rx,Rx⟩+ α ⟨Lx, Lx⟩ = ∥Rx∥2 + α∥Lx∥2 > 0.
Im letzten Schritt gilt =, da die Radon-Transformation nach Satz 3.1.6 invertierbar ist
und somit Rx = 0 gilt.
Lemma 4.1.4 (Hinreichendes Kriterium fur eine Tikhonov-Regularisierung). Eine
Losung von (ATA + αLTL)x = AT b ist ein Minimum von Φα aus Definition 4.1.1
und somit eine Tikhonov-Regularisierung.
Beweis. Es sei x ∈ Rn und es gelte (ATA + αLTL)x = AT b. Wir zeigen, dass fur alle
x = x ∈ Rn gilt Φα(x) ≤ Φα(x).
Φα(x)− Φα(x) =∥Ax− b∥2 + α∥Lx∥2 − ∥Ax− b∥2 − α∥Lx∥2
=⟨Ax− b, Ax− b⟩+ α⟨Lx, Lx⟩ − ⟨Ax− b, Ax− b⟩ − α⟨Lx, Lx⟩
=⟨Ax,Ax⟩ − 2⟨Ax, b⟩+ α⟨Lx,Lx⟩ − ⟨Ax,Ax⟩+ 2⟨Ax, b⟩ − α⟨Lx, Lx⟩
=2⟨Ax, b⟩ − 2⟨Ax, b⟩+ ⟨x, (ATA+ αLTL)x⟩ − ⟨x, (ATA+ αLTL)x⟩
=⟨Ax, b⟩ − ⟨Ax, b⟩+ ⟨x, (ATA+ αLTL)x− AT b⟩
− ⟨x, (ATA+ αLTL)x− AT b⟩
Einsetzen liefert:
=⟨x− x,AT b⟩+ ⟨x, (ATA+ αLTL)(x− x)⟩
=⟨x− x, (ATA+ αLTL)(x− x)⟩ = ∥A(x− x)∥2 + α∥L(x− x)∥2 ≥ 0
Bemerkung 4.1.5. Wahlen wir A = R in Lemma 4.1.4, so gilt im letzten Schritt des
Beweises sogar > 0, wie in Lemma 4.1.3. Hieraus folgt die Eindeutigkeit der Tikhonov-
Regularisierung fur unser Problem.
In den Lemmata 4.1.2 bis 4.1.4 haben wir gezeigt, dass es fur ein g ∈ S(C) genau
eine Tikhonov-Regularisierung gibt. Diese ist durch die Losung eines linearen Glei-
chungssystems gegeben, dessen Systemmatrix symmetrisch positiv definit ist. Fur ein
numerisches Losungsverfahren zeigen wir noch die Aquivalenz zu einer Least-Squares
Losung, um dann die richtige Implementation des Conjugate-Gradient-Verfahren zu
wahlen.
4 Regularisierungsmethoden 19
Korollar 4.1.6 (Aquivalenz zu einer Least-Squares Losung). Die Losung des Glei-
chungssystems (ATA + αLTL)x = AT b ist eine Least-Squares Losung des Gleichungs-
systems [A√αL
]=
[b
0
]. (4.3)
Dabei hat die 0 genau so viele Zeilen wie L.
Beweis. Wir betrachten die Normalengleichung.
[AT
√αLT
] [ A√αL
]x =
[AT
√αLT
] [ b
0
]⇔ (ATA+ αLTL)x = AT b
Wir wollen nun also eine Least-Squares Losung zu der Gleichung (4.3) finden. Wir
verwenden die CGLS-Methode aus dem Buch von Bjorck [Bjo96, Algorithmus 7.4.2]
und setzen unser Problem ein. Das L wollen wir als Differentialoperator wahlen: Sind
Dx und Dy diskrete Richtungsableitungen, so setzen wir L =
[Dx
Dy
]und es folgt
∥Lf∥22 = ∥Dxf∥22 + ∥Dyf∥22.
Dadurch wollen wir vermeiden, dass sich strahlenformige Artefakte bilden, welche einen
großen Gradienten haben.
4 Regularisierungsmethoden 20
Algorithm 1 CGLS-Methode fur Least-Squares Losungen Vgl. [Bjo96, Algorithmus7.4.2]
Initialize:
p← s←R∗(g −Rx)− αLTLx, γ ← ∥s∥22for k = 1, 2, . . . und solange γk > tol do
q1 ←Rpq2 ←
√αLp
δ ← γ/(∥q1∥22 + ∥q2∥22)x← x+ δp
s← s− δ(R∗q1 +√αLT q2)
γ ← ∥s∥22β ← γ/γ
p← s+ βp
γ ← γ
end for
(a) α = 30 (b) α = 700
Abbildung 4.1: Rekonstruktionen mit der Tikhonov-Regularisierungsmethode. Manbeachte, dass mit hohem Regularisierungparameter Kanten verlorengehen.
4 Regularisierungsmethoden 21
4.2 Total-Variation-Regularisierung
Nachdem wir uns im vorherigen Abschnitt die Tikhonov-Regularisierung angesehen
haben, werden wir uns in diesem Abschnitt mit der Total-Variation-Regularisierung
(TV-Regularisierung) beschaftigen. Wir haben festgestellt, dass bei der Tikhonov-
Regularisierung die Kanten im Bild verloren gehen. Wir erhoffen uns von der TV-
Regularisierung, dass die Kanten erhalten bleiben, wir aber trotzdem eine Glattung
erzielen. Vorgestellt wurde die TV-Regularisierung 1992 von Rudin, Osher und Fatemi
[ROF92] zur Reduzierung von Rauschen in Bildern. Dazu definieren wir die Total-
Variation eines Bildes.
Definition 4.2.1 (Total-Variation). Es sei u ∈ L1(Ω). Die Total-Variation ist definiert
durch
|u|BV (Ω) := sup
∫Ω
u divφ dx∣∣φ ∈ C1
0(Ω,Rd), ∥φ∥∞ ≤ 1
.
Bemerkung 4.2.2.
• Ist u schwach differenzierbar, so gilt nach der mehrdimensionalen partiellen In-
tegration ∫Ω
u divφ dx = −∫Ω
∇uφ dx
und damit
|u|BV (Ω) =
∫Ω
|∇u| dx. (4.4)
• Im diskreten Bildmodell sind unsere Ansatzfunktionen schwach differenzierbar,
daher konnen wir die Formulierung (4.4) verwenden.
Die regularisierte Variante v eines verrauschten Bildes v ist dann der Minimierer eines
Funktionals.
v = argminu
1
2∥u− v∥22 + λ|u|BV (Ω) (4.5)
Diese Problemstellung wird auch als ROF-Modell bezeichnet. Nun konnen wir nicht
analog zur Tikhonov-Regularisierung vorgehen, da unser Funktional nicht differenzier-
bar ist. Wir wollen deshalb das Arrow-Hurwicz-Verfahren herleiten und damit eine
4 Regularisierungsmethoden 22
TV-Regularisierung errechnen. Die Methode ist auch unter dem Namen Chambolle-
Pock bekannt (vgl. [CP11]). Wir orientieren uns an einem Paper von Sidky, Jørgensen
und Pan [SJP12]. Fur einige Grundlagen der konvexen Funktionale und der Fenchel-
Rockafellar-Dualitat ziehen wir außerdem das Buch von Bredies und Lorenz [BL10,
Kapitel 6.2] heran. Weiter sei auf das Paper von Schiffer und Bredies [SB14] verwie-
sen, in dem die TV-Regularisierung auf das Problem der Metall-Artefakte angewendet
wird.
Eine Alternative zu unserem Vorgehen ware, das Funktional durch eine differenzierbare
Approximation zu ersetzen und dieses, ahnlich wie bei der Tikhonov-Regularisierung,
durch ein Gradienten-Verfahren zu minimieren (vgl. [DV97]). Dieses Vorgehen birgt
allerdings einige Probleme, so muss man etwa zwischen einer starken Approximati-
on und einer glatten Approximation abwagen. Eine glatte Approximation fuhrt dazu,
dass Kanten verschwimmen, ahnlich wie bei der Tikhonov-Regularisierung. Eine star-
ke Approximation fuhrt zu Schwierigkeiten mit der Differenzierbarkeit, sodass es zu
Problemen mit der Konvergenz kommen kann (vgl. [BO13, Abschnitt 8.3]).
Bemerkung 4.2.3. Wir werden in diesem Abschnitt Funktionale mit Wertebereich
R∞ = R ∪ ∞ verwenden, dabei bedienen wir folgende Arithmetik:
t ≤ ∞ ∀t ∈ R∞
a+∞ =∞ ∀a ∈ R∞
a · ∞ =∞ ∀a > 0
0 · ∞ = 0.
Alle anderen Operationen gelten als nicht definiert. Ein Funktional nennen wir ei-
gentlich, falls es nicht konstant ∞ ist. Wir wollen im Folgenden annehmen, dass alle
Funktionale eigentlich sind.
Um in die Theorie der konvexen Funktionale einzusteigen, werden wir, nach Bredies
und Lorenz [BL10, Kapitel 6.2], zunachst einige Definitionen heranziehen.
Definition 4.2.4 (Konvexes Funktional). Ein Funktional F : X → R∞ auf einem
normierten Vektorraum X heißt konvex, falls fur alle u, v ∈ X und λ ∈ (0, 1) gilt
F (λu+ (1− λ)v) ≤ λF (u) + (1− λ)F (v)
und strikt konvex, falls fur u = v sogar < gilt.
4 Regularisierungsmethoden 23
Definition 4.2.5 (Unterhalbstetigkeit). Ein Funktional F : X → R∞ auf einem to-
pologischen Raum X heißt (folgen-)unterhalbstetig, falls fur jede Folge (un) und u ∈ X
mit limn→∞
un = u folgt
F (u) ≤ lim infn→∞
F (un).
Definition 4.2.6 (Duales Funktional). Es sei F : X → R∞ ein Funktional auf dem
reellen Banachraum X. Dann ist das duale Funktional F ∗ : X∗ → R∞ definiert durch
F ∗(ω) = supu∈X⟨ω, u⟩X∗×X − F (u).
Bemerkung 4.2.7.
• Das Supremum geht nur uber diejenigen Elemente u ∈ X, fur die F (u) = ∞gilt, da die Subtraktion von ∞ nicht definiert ist. Da wir aber annehmen, dass F
eigentlich ist, gibt es mindestens ein Element u ∈ X, fur welches F (u) =∞ gilt.
Somit nimmt das Supremum einen Wert in R∞ an. F ∗ ist damit wohldefiniert.
• Fur F : X → R und (u∗, ω∗) ∈ X ×X∗ mit F (u∗) =∞ gilt
F ∗(ω∗) = supu∈X⟨ω∗, u⟩X∗×X − F (u) ≥ ⟨ω∗, u∗⟩X∗×X − F (u∗) (4.6)
und damit F ∗(ω∗) + F (u∗) ≥ ⟨ω∗, u∗⟩. Diese Ungleichung nennen wir Fenchel-
Ungleichung.
Definition 4.2.8 (Subgradient, Subdifferential). Es sei X ein reeller normierter Vek-
torraum und F : X → R∞ ein konvexes Funktional. Ein ω ∈ X∗ wird ein Subgradient
von F in u ∈ X genannt, falls
F (u) + ⟨ω, v − u⟩ ≤ F (v) ∀v ∈ X (4.7)
gilt. Die Menge ∂F (u) := ω ∈ X∗|ω ist Subgradient von F in u wird Subdifferential
genannt.
Bemerkung 4.2.9.
• Geometrisch entsprechen die Subgradienten in u ∈ X den Steigungen der Gera-
den, welche durch den Punkt (u, F (u)) gehen und uberall unterhalb des Graphen
von F liegen. Das bedeutet auch, dass ∂F an differenzierbaren Stellen u maximal
einelementig ist. Dann gilt ∂F (u) = DF (u).
4 Regularisierungsmethoden 24
• Es gilt 0 ∈ ∂F (u) genau dann, wenn u ein Minimum von F ist, denn
F (u) ≤ F (v) ∀v ∈ X ⇔ F (u) + ⟨0, v − u⟩ ≤ F (v) ∀v ∈ X.
• Es gilt ∂(F + G) = ∂F + ∂G, falls F an einem Punkt u ∈ X stetig ist, fur den
G(u) = ∞ = F (u) gilt. Einen Beweis dafur finden wir in Bredies und Lorenz
[BL10, Satz 6.51].
Mit diesen Definitionen konnen wir den folgenden Satz formulieren.
Satz 4.2.10 (Fenchel-Rockafellar-Dualitat vgl. [BL10, Satz 6.68]). Es seien F : X →R∞, G : Y → R∞ zwei konvexe und unterhalbstetige Funktionale auf den reellen Ba-
nachraumen X und Y . Weiter sei K : X → Y eine lineare, stetige Abbildung und das
Minimierungsproblem
minu∈X
F (u) +G(Ku) (4.8)
besitze eine Losung. Dann gilt
maxω∈Y ∗
−F ∗(−K∗ω)−G∗(ω) = minu∈X
F (u) +G(Ku).
Einen Beweis finden wir im Buch von Bredies und Lorenz [BL10, Satz 6.68].
Das Maximierungsproblem
maxω∈Y ∗
−F ∗(−K∗ω)−G∗(ω)
nennen wir das duale Problem zu (4.8). Unter dem primal-dualen Problem, verstehen
wir die simultane Losung beider Probleme.
Lemma 4.2.11. Es seien F : X → R∞ und G : Y → R∞ konvexe unterhalbstetige
Funktionale auf reellen Banachraumen X und Y . Das Paar (u∗, ω∗) ∈ X×Y ist genau
dann Losung des primal-dualen Problems, falls (u∗, ω∗) Sattelpunkt der Funktion
L : X × Y ∗ → R L(u, ω) := ⟨ω,Ku⟩+ F (u)−G∗(ω)
ist. Das heißt, falls L(u∗, ω∗) ≥ L(u∗, ω) fur alle ω ∈ Y ∗ und L(u∗, ω∗) ≤ L(u, ω∗) fur
alle u ∈ X gilt.
4 Regularisierungsmethoden 25
Beweis. Ist u∗ ∈ X eine Losung des primalen Problems und ω∗ ∈ Y ∗ eine Losung des
dualen Problems, so gilt einerseits
L(u∗, ω∗) ≤ supω∈Y ∗
L(u∗, ω) = supω∈Y ∗⟨ω,Ku∗⟩+ F (u∗) +G∗(ω) = F (u∗) +G(Ku∗)
und andererseits
L(u∗, ω∗) ≥ infu∈X
L(u, ω∗) = − supu∈X−L(u, ω∗) = − sup
u∈X−⟨ω∗, Ku⟩ − F (u)−G∗(ω∗)
= −F ∗(−K∗ω∗)−G∗(ω∗).
Nach Satz 4.2.10 sind jedoch die rechten Seiten beider Gleichungen gleich und somit
gilt anstatt ≤ und ≥ jeweils =. Insbesondere gilt
L(u∗, ω) ≤ supω∈Y ∗
L(u∗, ω∗) = L(u∗, ω∗) = infu∈X
L(u, ω∗) ≤ L(u, ω∗) ∀u ∈ X und ω ∈ Y ∗.
Also ist (u∗, ω∗) ein Sattelpunkt von L.
Zeigen wir nun die andere Richtung: Ist (u∗, ω∗) ∈ X ×Y ein Sattelpunkt von L, dann
gilt mit der Fenchel-Ungleichung fur beliebige u ∈ X und ω ∈ Y :
F (u) +G(Ku) ≥ ⟨ω∗, Ku⟩+ F (u)−G∗(ω∗) = L(u, ω∗) ≥ L(u∗, ω∗)
= supω∈Y ∗
L(u∗, ω) = supω∈Y ∗⟨ω,Ku∗⟩+ F (u∗)−G∗(ω) = F (u∗) +G(Ku∗)
und
−F ∗(−K∗ω)−G∗(ω) ≤ ⟨ω,Ku∗⟩+ F (u∗)−G∗(ω) = L(u∗, ω) ≤ L(u∗, ω∗)
= infu∈x
L(u, ω∗) = − supu∈X−L(u, ω∗)
= − supu∈X−⟨ω∗, Ku⟩ − F (u) +G∗(ω∗) = −F ∗(−K∗ω∗)−G∗(ω∗).
Bemerkung 4.2.12. Wir benutzen in dem Beweis, dass fur das biduale Funktional
G∗∗ = G gilt. Dies folgt aus der Annahme, dass G konvex und unterhalbstetig ist (vgl.
[BL10, Korollar 6.58 und Lemma 6.63]).
Das Arrow-Hurwicz-Verfahren wird eine Folge von Punkten in X × Y ∗ berechnen, die
gegen den Sattelpunkt von L konvergieren. Dafur fuhren wir die Resolvente ein.
4 Regularisierungsmethoden 26
Definition 4.2.13 (Resolvente). Es sei F : X → R∞ ein konvexes Funktional. Die
Abbildung (id+ σ∂F )−1 : X → X, die u ∈ X den eindeutigen Minimierer
argminv∈X
∥v − u∥2X2
+ σF (v) =: (id+ σ∂F )−1(u)
zuordnet, heißt Resolvente von F zu σ > 0.
Die Wohldefiniertheit der Resolvente ist hier a priori nicht klar, da die Existenz und
Eindeutigkeit des Minimums nicht offensichtlich ist. Das zu minimierende Funktional
ist jedoch strikt konvex und koerziv, woraus die Eindeutigkeit und Existenz folgt. Fur
eine detaillierte Auseinandersetzung sei auf Bredies und Lorenz [BL10] verwiesen.
Wir wollen nun noch die Bezeichnung der Resolvente rechtfertigen.
Lemma 4.2.14. Es gilt (id+ σ∂F )−1(u) = w genau dann, wenn u ∈ w + σ∂F (w).
Beweis.
(id+ σ∂F )−1(u) = w ⇔ w = argminξ∈X
∥u− ξ∥222
+ σF (ξ)
⇔ 0 ∈ ∂
(∥u− ·∥22
2+ σF
)(w)
⇔ 0 ∈ w − u+ σ∂F (w)
⇔ u ∈ w + σ∂F (w)
Lemma 4.2.15 (Einige Rechenregeln fur die Resolvente). Es sei F1 : X → R∞ ein
konvexes unterhalbstetiges Funktional. Dann gilt
1. fur α ∈ RF2 = F1 + α⇒ (id+ σ∂F2)
−1 = (id+ σ∂F1)−1.
2. fur beliebiges ω0 ∈ X
F2(u) = F1(u) + ⟨ω0, u⟩ ⇒ (id+ σ∂F2)−1(u) = (id+ σ∂F1)
−1(u− σω0).
4 Regularisierungsmethoden 27
3. Ist F2 : Y → R∞ ebenfalls ein konvexes, unterhalbstetiges Funktional, so gilt
F3(u, v) = F1(u) + F2(v)⇒ (id+ σ∂F3)−1(u, v) =
((id+ σ∂F1)
−1(u)
(id+ σ∂F2)−1(v)
).
Beweis. 1. Folgt aus der Gleichheit argmin f(x) + α = argmin f(x).
2. Wir rechnen die Behauptung nach, dabei bedienen wir uns ebenfalls daran, Kon-
stanten zu addieren ohne das Minimum zu andern.
(id+ σ∂F2)−1(u) = argmin
v∈X
∥v − u∥222
+ σ(F1(v) + ⟨ω0, v⟩)
= argminv∈X
⟨v − u, v − u⟩+ 2⟨σω0, v⟩2
+ σF1(v)
= argminv∈X
⟨v − u, v − u⟩+ 2⟨σω0, v − u⟩+ ⟨σω0, σω0⟩2
+ σF1(v)
= argminv∈X
∥v − (u− σω0)∥222
+ σF1(v)
=(id+ σ∂F1)−1(u− σω0)
3. Bei einer Summe konnen wir die Summanden einzeln minimieren.
(id+ σ∂F3)−1(u, v) = argmin
(w,x)∈X×Y
∥u− w∥222
+ F1(w) +∥v − x∥22
2+ F2(x)
=
(argminw∈X
∥u−w∥222
+ F1(w)
argminx∈Y∥v−x∥22
2+ F2(x)
)
=
((id+ σ∂F1)
−1(u)
(id+ σ∂F2)−1(v)
)
Wir wollen nun mithilfe der Resolvente das Arrow-Hurwicz-Verfahren herleiten. Dazu
betrachten wir fur einen beliebigen Startwert (u0, ω0) ∈ X × Y ∗ die Einschrankungen
von L.
Lω0(u) = L(u, ω0) = ⟨ω0, u⟩+ F (u)−G∗(ω0)
Lu0(ω) = L(u0, ω) = ⟨ω, u0⟩+ F (u0)−G∗(ω)
4 Regularisierungsmethoden 28
Wir wollen L bezuglich der u-Variable minimieren und bezuglich der ω-Variable ma-
ximieren, daher sehen wir uns die Resolventen von Lω0 und −Lu0 an. Dabei hilft uns
Lemma 4.2.15.
(id+ σ∂Lω0)−1(u) = (id+ σ∂F )−1(u− σA∗ω0)
(id+ σ∂(−Lu0))−1(ω) = (id+ σ∂G∗)−1(ω + σAu0)
Die Aussage (u∗, ω∗) ∈ X×Y ∗ ist Sattelpunkt von L konnen wir nun durch Aquivalenz-
umformungen in eine Fixpunktgleichung umformen, um daraus das Arrow-Hurwicz-
Verfahren zu gewinnen.
(u∗, ω∗) ist Sattelpunkt von L⇔
0 ∈ ∂Lω∗(u∗)
0 ∈ ∂(−Lu∗)(ω∗)
⇔
u∗ ∈ u∗ + σ∂Lω∗(u∗)
ω∗ ∈ ω∗ + τ∂(−Lu∗)(ω∗)
⇔
u∗ = (id+ σ∂F )−1(u∗ − σA∗ω∗)
ω∗ = (id+ τ∂G∗)−1(ω∗ + τAu∗)
Aus dieser Fixpunktgleichung konnen wir nun den Algorithmus entwickeln, indem wir
zusatzlich die Variablen (u∗, ω∗) einfuhren. Außerdem fuhren wir am Ende eines Itera-
tionsschrittes noch eine Extrapolation in der primalen Variable u durch. Fur genauere
Informationen siehe Bredies und Lorenz [BL10].
Algorithm 2 Arrow-Hurwicz-Verfahren mit linearem primalen Extragradienten. (vgl.[BL10, S. 388])
u∗0 ← 0
for k = 1, 2, 3, . . . do
ω∗ ← (id+ τ∂G∗)−1(ω∗ + τKu∗)
u∗k ← (id+ σ∂F )−1(u∗ − σK∗ω∗)
u∗ ← 2u∗k − uk−1
end for
Der Algorithmus konvergiert unter der Voraussetzung τσ∥K∥2 < 1 fur k → ∞ gegen
den Sattelpunkt von L (vgl. [CP11]). Eine Frage, die sich nun noch stellt, ist, wie sich die
4 Regularisierungsmethoden 29
Resolventen berechnen lassen. Dafur schauen wir uns zunachst das ROF-Modell (4.5)
an, welches wir spater mit der ART-Rekonstruktion kombinieren werden (ART+TV),
indem wir sukzessive ART-Schritte ausfuhren und die Zwischenergebnisse mit dem
ROF-Modell regularisieren (vgl. [SP08]).
Wir wenden die Theorie nun auf das ROF-Modell an. Dafur sei X = Rm×n der diskrete
Bildraum und Ω = 1, 2, . . . ,m×1, 2, . . . , n. Weiter sei ein verrauschtes Bild u0 ∈ X
gegeben. Wir setzen dann
F : X → R∞ F (u) =1
2∥u0 − u∥22 =
∑(i,j)∈Ω
(u0i,j − ui,j)
2
und
G : X2 → R∞ G(v) = λ|v|1 = λ∑
(i,j)∈Ω
|(v1)i,j|+ |(v2)i,j|,
sowie
K : X → X2 K = ∇h =
(Dx
Dy
)wobei Dx und Dy diskrete Richtungsableitungen nach x und y sind.
Wir berechnen (id+ σ∂F )−1:
(id+ σ∂F )−1(u) = argminv∈X
∥v − u∥222
+ σ1
2∥u0 − v∥22
Der Term ist differenzierbar mit Maximum v = u+σu0
1+σund somit gilt
(id+ σ∂F )−1(u) =u+ σu0
1 + σ.
Zur Berechnung von (id+ τ∂G∗)−1 bestimmen wir zunachst G∗:
G∗(ω) = supv∈X2
⟨ω, v⟩ −G(v) = supv∈X2
⟨ω, v⟩ − λ∑
(i,j)∈Ω
|(v1)i,j|+ |(v2)i,j|
Dazu unterscheiden wir zwei Falle:
1. Ist |ωi,j|∞ ≤ λ fur alle (i, j) ∈ Ω, so gilt G∗(ω) = 0.
2. Ist |ωi,j|∞ > λ fur ein (i, j) ∈ Ω, so gilt G∗(ω) =∞.
4 Regularisierungsmethoden 30
Zusammen gilt also
G∗(ω) = I|ω|∞≤λ(ω) =
0 , falls |ω|∞ ≤ λ
∞ , sonst.
Damit konnen wir nun die Resolvente berechnen:
(id+ τ∂G∗)−1(ω) = argminξ∈X2∗
∥ξ − ω∥222
+ τG∗(ξ) = argminξ∈X2∗
∥ξ − ω∥222
+ τI∥ξ∥∞≤λ(ξ)
Der Term nimmt nur einen endlichen Wert an, falls |ξi,j|∞ ≤ λ fur alle (i, j) ∈ Ω und
ist am kleinsten, wenn ∥ξ − ω∥22 klein ist. Das heißt die Resolvente (id + τ∂G∗)−1 ist
die Projektion auf die Menge ω ∈ X2∗∣∣|ωi,j|∞ ≤ λ ∀(i, j) ∈ Ω.
(id+ τ∂G∗2)
−1(ω)i,j = P∞λ (ω)i,j =
ωi,j
max1, λ−1|ωi,j|
Der Algorithmus zur Losung des ROF-Modells sieht dann wie folgt aus:
Algorithm 3 Arrow-Hurwicz-Verfahren mit linearem primalen Extragradienten furdas ROF-Modell
Require: στ < ∥∇∥−2
for k = 1, 2, 3, . . . do
ωk+1 ← Pλ(ω + τ∇u)uk+1 ← uk+σ(−div(ω)+u0)
1+σ
u← 2uk+1 − uk
end for
4 Regularisierungsmethoden 31
(a) λ = 10−4 (b) λ = 10−1
Abbildung 4.2: ART+TV Rekonstruktionen. Man beachte, dass Kanten erhalten blei-ben, insbesondere aber kleine Details ineinander ubergehen (mittlerer,unterer Bereich).
Nachdem wir das Standardmodell betrachtet haben, werden wir noch einen zweiten
Weg wahlen und das Modell auf die Radon-Transformation anpassen. Dabei orientieren
wir uns an dem Paper von Schiffer und Bredies [SB14].
Das ROF-Modell setzt sich aus zwei Funktionalen zusammen: Auf der einen Seite
aus dem Regularisierungsterm λ|u|BV (Ω) und auf der anderen Seite aus dem Daten-
term 12∥u0 − u∥22. Nach dem Studium der Tikhonov-Regularisierung liegt es nahe, im
Datenterm den Abstand im Datenraum zu messen. Auf diese Weise wird neben der
Minimierung des Regularisierungsterms auch das Bild aus den Daten rekonstruiert.
Dafur sei Y := RM×N der diskrete Datenraum zu X. Wir betrachten also fur gegebene
Daten g ∈ Y
argminf∈X
1
2∥Rf − g∥22 + λ∥∇f∥1.
Wir legen nun wieder die Funktionale fest, um die Theorie anwenden zu konnen. Dies
wollen wir auf folgende Weise tun:
F : X → R∞ F (u) = 0
4 Regularisierungsmethoden 32
Fur η ∈ R setzen wir
G : Y ×X2 → R∞
G(p, v) =1
2∥η−1p− g∥22 + λ|v|1 =
∑(i,j)∈Ω
1
2(η−1pi,j − gi,j)
2 + λ(|(v1)i,j|+ |(v2)i,j|)
sowie
K : X → Y ×X2 K =
[ηR∇h
].
Das η ermoglicht es uns, das Verhaltnis der Normen von ηR und ∇ zu kontrollieren,
was uns Vorteile hinsichtlich der Konvergenzgeschwindigkeit bringt.
Es gilt (id+ σ∂F )−1 = id, denn
(id+ σ∂F )−1(u) = argminv∈Rm×n
∥v − u∥222
− 0 = u.
Nach Lemma 4.2.15 kann die Resolvente von G∗ komponentenweise ausgerechnet wer-
den:
(id+ τ∂G∗)−1(p, v) =
((id+ τ∂G∗
1)−1(p)
(id+ τ∂G∗2)
−1(v)
)Dabei ist G1(p) =
12∥η−1p− g∥22 und damit
G∗1(ξ) = sup
p∈Y⟨ξ, p⟩ − 1
2∥η−1p− g∥22 = (η2 − 1
2η)∥ξ∥22 + η⟨ξ, g⟩.
Die Resolvente berechnet sich nun wie folgt:
(id+ τ∂G∗1)
−1(ξ) = argminχ∈RM×N
∥χ− ξ∥222
+ τ
((η2 − 1
2η)∥χ∥22 + η⟨χ, g⟩
)und wir konnen durch Differenzieren wieder das Minimum bestimmen:
(id+ τ∂G∗1)
−1(ξ) =ξ − τηg
1 + τ(2η2 − η)
Fur G2(ω) = λ∥ω∥1 konnen wir analog zum ROF-Modell vorgehen. Die Resolvente
(id+ τ∂G∗2)
−1 ist also die Projektion auf die Menge |ω|∞ ≤ λ:
(id+ τ∂G∗2)
−1(ω)i,j = P∞λ (ω)i,j =
ωi,j
max1, λ−1|ωi,j|
4 Regularisierungsmethoden 33
Mit den Resolventen konnen wir nun das Arrow-Hurwicz-Verfahren formulieren.
Algorithm 4 Arrow-Hurwicz-Verfahren mit linearem primalen Extragradienten furdie TV-Rekonstruktion.
Require: στ < (η2∥R∥2 + ∥∇∥2)−1
for k = 1, 2, 3, . . . do
ωk+1 ← Pλ(ω + τ∇u)vk+1 ← vk + τη(Ru− g))/(1 + τ(2η2 − η))
uk+1 ← uk − σ(−div(ω) +R∗v)
u← 2uk+1 − uk
end for
(a) λ = 10−4 (b) λ = 10−1
Abbildung 4.3: TV-Regularisierungen
Der Unterschied der Bilder wird klar, wenn man die Hohenlinen eines Schnitts des
Bildes betrachtet.
4 Regularisierungsmethoden 34
0 50 100 150 200 250−0.02
−0.01
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
λ = 10−1
λ = 10−4
Abbildung 4.4: Hohenlinien von TV-Regularisierungen mit verschiedenen Regularisie-rungsparametern.
35
5 Anpassung der Algorithmen auf das
Problem der Metall-Artefakt
Reduktion
In diesem Kapitel wollen wir die von uns in Kapitel 3 und 4 untersuchten Algorithmen
auf das Problem der Metall-Artefakte anpassen. Dabei nehmen wir stets an, dass uns
die Menge Ω0 ⊂ Ω ⊂ R2, in der sich das Metall befindet, bekannt ist. Da wir Ω0
kennen, wissen wir auch, welche Daten im Sinogramm gestort sind, indem wir Λ :=
(θ, s) ∈ C∣∣RIΩ0(θ, s) = 0 betrachten. In diesem Kapitel soll Λθ = s ∈ R
∣∣(θ, s) ∈ Λdie Einschrankung von Λ auf θ × R bezeichnen. Grundidee wird es sein, im Sinne des
Themas dieser Arbeit, die gestorten Daten zu ignorieren und zu sehen, wie sich die
bildbasierten Regularisierungsmethoden auswirken. Wir wollen aber mit dem in der
Einleitung angekundigten datenbasierten Algorithmus starten.
5.1 Lineare Interpolation und gefilterte Ruckprojektion
Bei der gefilterten Ruckprojektion ersetzen wir die vom Metall gestorten Daten im
Sinogramm, indem wir sie durch umgebende Werte interpolieren. Sind also Daten g ∈S(C) gegeben, so erstellen wir g ∈ S(C), indem wir fur jedes θ ∈ S1 die Projektion
g(θ, ·) durch g(θ, ·) ersetzen. Dabei ist g(θ, ·)∣∣Λθauf jedem Intervall von Λθ linear, sodass
g(θ, ·) und g(θ, ·) auf dem Rand von Λθ ubereinstimmen. Auf g wenden wir dann die
gefilterte Ruckprojektion an.
Dieses Vorgehen ist datenbasiert, da wir die Daten andern und dann ein normales Re-
konstruktionsverfahren anwenden. Die meisten Verfahren, die in der Praxis eingesetzt
werden, sind datenbasiert.
5 Anpassung der Algorithmen auf das Problem der Metall-Artefakt Reduktion 36
(a) Gefilterte Ruckprojektion von durch lineare
Interpolation aufbereiteten Daten
200 220 240 260 2800
10
20
30
40
50
60
70
80
(b) Graph einer interpolierten Projektion
Abbildung 5.1: Lineare Interpolation und gefilterte Ruckprojektion.
5.2 Algebraic-Reconstruction-Technique
Bei der Algebraic-Reconstruction-Technique konnen wir die gestorten Daten ignorieren,
dazu streichen wir die entsprechenden Zeilen aus dem Gleichungssystem. Dies hat zur
Folge, dass die affin-linearen Unterraume Hθ, auf die wir sukzessive projizieren, großer
werden.
Hθ :=
f ∈ S(R2)
∣∣Rθf∣∣R\Λθ= g(θ, ·)∣∣R\Λθ
Auf die Projektionen Pθ wirkt sich das folgendermaßen aus:
Pθf = f +R∗θ(RθR∗
θ)−1hf
hf (s) :=
g(θ, ·)−Rθf , falls s ∈ λθ
0 , sonst
Anschaulich betrachtet ignorieren wir bei dem Zuruckprojizieren der Projektionen die
vom Metall gestorten Daten. Wir lassen also den Bereich des Bildes, auf den sich die
Daten auswirken wurden, unangetastet.
Wir erhoffen uns von der Kombination aus ART und der TV-Regularisierung eine Re-
5 Anpassung der Algorithmen auf das Problem der Metall-Artefakt Reduktion 37
(a) Ohne TV-Regularisierung (b) Mit TV-Regularisierung
Abbildung 5.2: Rekonstrunktion mit ART+TV
duktion der Artefakte. In Abbildung 5.2 sehen wir Ergebnisse der ART-Rekonstruktion
mit und ohne TV-Regularisierung.
5.3 Tikhonov-Regularisierung
Um die Tikhonov-Regularisierung anzupassen, betrachten wir das Funktional, welches
wir minimieren wollen:
Φα(x) = ∥Ax− g∥22 + α∥Lx∥22
Wir wissen nun, welche Daten wir im Funktional nicht beachten wollen, das heißt, wir
passen den Datenterm ∥Ax−g∥22 an. Schauen wir uns daher die Situation im diskreten
Fall an:
Sei Λ ⊂ 1, 2, . . . ,m × 1, 2, . . . , n. Dann ersetzen wir
∥Ax− g∥22 =∑
(i,j)∈Ω
((Ax)i,j − gi,j)2
durch
∥Ax− g∥2Ω\Λ =∑
(i,j)∈Ω\Λ
((Ax)i,j − gi,j)2.
Wir lassen die gestorten Daten in der Summierung also aus.
5 Anpassung der Algorithmen auf das Problem der Metall-Artefakt Reduktion 38
(a) α = 30 (b) α = 500
Abbildung 5.3: Tikhonov-Regularisierung mit unterschiedlichen Regularisierungspara-metern. Wir stellen fest, dass bei zu kleinem Regularisierungsparameternoch Artefakte vorhanden sind, bei zu großem das Bild unscharf wird.
Fur das Least-Squares Problem bedeutet das, dass in der Matrix A und den Daten
g entsprechende Zeilen gestrichen werden. In der Implementierung konnen wir analog
dazu in Algorithmus 1 auch q1∣∣Λ = 0 setzen. Ergebnisse der Tikhonov-Regularisierung
sind der Abbildung 5.3 zu entnehmen.
5.4 Total-Variation-Regularisierung
Bei der TV-Regularisierung konnen wir im Ansatz analog zur Tikhonov-Regularisie-
rung vorgehen. Das heißt, wir ersetzen den Datenterm ∥Rx− g∥22 durch ∥Rx− g∥2Ω\Λ.
Fur den Algorithmus mussen wir folglich die Resolvente (id+ τ∂G∗1)
−1 anpassen.
G∗1(ξ) = sup
p∈Y⟨ξ, p⟩ − ∥η−1p− g∥2Ω\Λ
= supp∈Y
∑(i,j)∈Ω
ξi,jpi,j −∑
(i,j)∈Ω\Λ
(η−1pi,j − gi,j)2
5 Anpassung der Algorithmen auf das Problem der Metall-Artefakt Reduktion 39
Falls ξi,j = 0 fur ein (i, j) ∈ Λ gilt, so giltG∗1(ξ) =∞, ansonstenG∗
1(ξ) = (η2− 12η)∥ξ∥22+
η⟨ξ, g⟩, wie in Abschnitt 4.2. Daraus konnen wir nun die Resolvente berechnen:
(id+ τ∂G∗1)
−1(ξ) = argminχ∈RM×N
∥χ− ξ∥222
+ τ
(η2 − 12η)∥χ∥22 + η⟨χ, g⟩ , falls χi,j = 0 ∀(i, j) ∈ Λ
∞ , sonst
Betrachten wir die Gleichung koeffizientenweise, so konnen wir das Minimum analytisch
berechnen.
(id+ τ∂G∗1)
−1(ξ)i,j =
(ξi,j − τηgi,j)/(1 + τ(2η2 − η)) , falls (i, j) ∈ Ω \ Λ
0 , falls (i, j) ∈ Λ
Mit dieser Resolvente wenden wir nun wieder das Arrow-Hurwicz-Verfahren an. Abbil-
dung 5.4 zeigt Ergebnisse mit und ohne Regularisierung.
(a) λ = 0 (b) λ = 10−4
Abbildung 5.4: TV-Rekonstruktion aus gestorten Daten. Das Bild kann exakt rekon-struiert werden.
40
6 Vergleich der Methoden
In diesem Kapitel wollen wir alle im Rahmen der Arbeit behandelten Methoden in
verschiedenen Situationen miteinander vergleichen. Grundlage hierfur bietet das Shepp-
Logan Phantom sowie ein Bild, welches realistischere Eigenschaften besitzt, da das
Shepp-Logan Phantom stuckweise-konstant und somit wie gemacht fur die Total-Vari-
ation ist. Wir verwenden das originale Shepp-Logan Phantom mit starken Kontrasten.
Die außere Ellipse in dem Phantom hat eine Intensitat von 1, wahrend die Intensitaten
im Rest des Bildes zwischen 0 und 0, 05 liegen. Dies wirkt sich auf die Rekonstruktionen
aus, da Bilder mit starken Kontrasten storungsanfalliger sind. Besonders bei Metall-
Artefakten spielt der Kontrast eine wesentliche Rolle, diese bilden sich starker aus, wenn
das Metall einen Kontrast abdeckt (vgl. Test 2). In unseren Darstellungen werden wir
Intensitaten, die großer als 0, 1 sind, weiß und solche, die zwischen 0 und 0, 1 liegen, in
Grauwerten darstellen.
1
0
0,03
Abbildung 6.1: Intensitaten des Shepp-Logan Phantoms
Beginnen wollen wir mit Vergleichen von verschiedenen Platzierung des Metalls inner-
halb des Bildes. Spater werden wir die Methoden auf ihr Verhalten im Falle mehrerer
Metallobjekte im Bild untersuchen und prufen wie sich Details, die nah am Metall
liegen, rekonstruieren lassen. Bei den Regularisierungsmethoden bestimmen wir den
Regularisierungsparameter moglichst so, dass sich keine Artefakte bilden, das Bild aber
6 Vergleich der Methoden 41
nicht zu stark geglattet wird.
Als Vergleichsmaß zwischen dem Originalbild u und dem Ergebnis der Methode v dient,
wie auch in der Masterarbeit von Koskela [Kos12], eine optische Begutachtung, sowie
der folgende relative Fehler:
ε2 :=∥u− v∥2Ω\Ω0
∥u∥2Ω\Ω0
Des Weiteren wollen wir den Fehler dicht am Metallobjekt beurteilen. Dafur sei eine
Umgebung ΩROI ⊂ Ω (Region of Interest) des Metallobjekts gegeben. Dann betrachten
wir
εROI2 :=
∥u− v∥2ΩROI\Ω0
∥u∥2ΩROI\Ω0
.
In den Ergebnissen wird der Bereich, in dem sich das Metall befindet, rot gekennzeich-
net. Als erstes Bild wird immer das Ausgangsbild gemeinsam mit dem Metallbereich
und ΩROI dargestellt.
6 Vergleich der Methoden 42
Test 1
Wir starten mit dem Beispiel aus Kapitel 5. Dabei liegt das Metallstuck im Inneren
des Phantoms, es stehen 180 Projektionen zur Verfugung und 2.3% der Daten werden
vom Metall verdeckt.
Die lineare Interpolation der Daten erzeugt schwache linienformige Artefakte. Bei der
Tikhonov-Regularisierung fallt auf, dass der Rand des Phantoms deutlich breiter ge-
worden ist. Als Begrundung ist hier die Intensitat des Randes zu nennen, welcher uber
dem maximal dargestellten Grauwert liegt. Die Glattung sorgt infolgedessen fur einen
breiter erscheinenden Rand. Beide TV-Methoden konnen das Ausgangsbild fast fehler-
frei wiederherstellen.
ε2 εROI2
Lineare Interpolation und
gefilterte Ruckprojektion0.0494 0.2625
Tikhonov-Regularisierung 0.1338 0.0278
ART+TV 2.0798 · 10−6 0.0166
TV-Regularisierung 8.2533 · 10−9 1.8895 · 10−4
Tabelle 6.1: Relative Fehler zu Test 1.
6 Vergleich der Methoden 43
(a) Metallobjekt und ROI
(b) Gefilterte Ruckprojektion (c) Tikhonov-Regularisierung
(d) ART+TV (e) TV-Regularisierung
Abbildung 6.2: Ergebnisse von Test 1.
6 Vergleich der Methoden 44
Test 2
Ein Fall, bei dem sich das Problem besonders auswirkt, ist wenn das Metallobjekt einen
Kontrast abdeckt beziehungsweise wenn Kontraste nah am Metallobjekt auftreten.
Auch hier fehlen 2.3% der Daten. Wir stellen fest, dass der datenbasierte Algorithmus
schlechtere Ergebnisse liefert als die bildbasierten. Bei der TV-Regularisierung wirkt
sich das Metall nur noch lokal und sehr schwach aus.
Die Ergebnisse unterscheiden sich stark von denen des ersten Tests. Besonders der
relative Fehler der gefilterten Ruckprojektion ist hier auf 45% gestiegen. Auch die TV-
Regularisierung kann das Bild nicht mehr exakt rekonstruieren.
ε2 εROI2
Lineare Interpolation und
gefilterte Ruckprojektion0.0948 0.4547
Tikhonov-Regularisierung 0.0803 0.1407
ART+TV 0.0067 0.1855
TV-Regularisierung 1.7290 · 10−3 0.1032
Tabelle 6.2: Relative Fehler zu Test 2.
6 Vergleich der Methoden 45
(a) Metallobjekt und ROI
(b) Gefilterte Ruckprojektion (c) Tikhonov-Regularisierung
(d) ART+TV (e) TV-Regularisierung
Abbildung 6.3: Ergebnisse von Test 2.
6 Vergleich der Methoden 46
Test 3
In diesem Test wollen wir uberprufen, ob sich Details, die nah am Metallobjekt liegen,
rekonstruieren lassen. In diesem Fall stehen weniger Daten uber die Bereiche nah am
Metall zur Verfugung. Daher erwarten wir, dass nicht alles rekonstruiert werden kann.
Hier liegt der Anteil an fehlenden Daten bei 4.2%.
Unsere Vermutung bestatigt sich, allerdings unterscheiden sich die Ergebnisse der Me-
thoden. Auch in diesem Test produziert die TV-Regularisierung das beste Ergebnis.
Wahrend die lineare Interpolation sowie die Tikhonov-Regularisierung Artefakte auf-
weisen, schaffen es die TV-basierten Verfahren diese vollstandig zu eliminieren. Jedoch
weisen alle Ergebnisse den Fehler auf, dass die linke Ellipse, die das Metallobjekt schnei-
det, nicht richtig rekonstruiert wird.
ε2 εROI2
Lineare Interpolation und
gefilterte Ruckprojektion0.0499 0.1211
Tikhonov-Regularisierung 0.0719 0.0282
ART+TV 8.0637 · 10−7 0.0206
TV-Regularisierung 4.8367 · 10−7 0.0100
Tabelle 6.3: Relative Fehler zu Test 3.
6 Vergleich der Methoden 47
(a) Metallobjekt und ROI
(b) Gefilterte Ruckprojektion (c) Tikhonov-Regularisierung
(d) ART+TV (e) TV-Regularisierung
Abbildung 6.4: Ergebnisse von Test 3.
6 Vergleich der Methoden 48
Test 4
Nun werden wir untersuchen wie sich die Verfahren verhalten, wenn sich mehrere
Metall-Objekte im Bild befinden. Es werden 6.7% der Daten vom Metall uberdeckt.
Besonders bei der Tikhonov-Regularisierung und der ART-Rekonstruktion mit TV-
Regularisierung bilden sich nah am Metallobjekt Artefakte in Richtung der anderen
Metallobjekte. Dies liegt an der Rekonstruktionsart: Da in diese Richtungen die Daten
fehlen, werden sie – wie etwa bei der ART-Methode – auch nicht dahin zuruckprojiziert.
Die TV-Rekonstruktionsmethode ist hier unanfallig.
ε2 εROI2
Lineare Interpolation und
gefilterte Ruckprojektion0.0958 0.4543
Tikhonov-Regularisierung 0.0817 0.1512
ART+TV 0.0101 0.2045
TV-Regularisierung 1.1193 · 10−3 0.0839
Tabelle 6.4: Relative Fehler zu Test 4.
6 Vergleich der Methoden 49
(a) Metallobjekt und ROI
(b) Gefilterte Ruckprojektion (c) Tikhonov-Regularisierung
(d) ART+TV (e) TV-Regularisierung
Abbildung 6.5: Ergebnisse von Test 4.
6 Vergleich der Methoden 50
Test 5
In Test 5 versuchen wir realistischere Bedingungen zu schaffen. Dafur verwenden wir
ein echtes Bild aus einem Computertomographen. Trotzdem sind die Ergebnisse auch
hier nur mit Einschrankungen zu bewerten. Das Bild hat nun ahnlichere Eigenschaften
wie ein echtes CT-Bild. Allerdings sind die Daten immer noch kunstlich aus dem Bild
erzeugt. Insbesondere betrachten wir das Bild auf dem gesamten Kontrastbereich, das
heißt es gibt vergleichsweise nur schwache Kontraste. 2.3% der Daten fehlen.
In der linken unteren Ecke sehen wir eine Vergroßerung von ΩROI . Wir bemerken, dass
die Details von den TV-Methoden am besten rekonstruiert werden konnen. Typischer-
weise wirkt die Tikhonov-Regularisierung unscharf.
ε2 εROI2
Lineare Interpolation und
gefilterte Ruckprojektion0.0101 0.0150
Tikhonov-Regularisierung 0.0150 0.0108
ART+TV 0.0022 0.0190
TV-Regularisierung 0.0019 0.0045
Tabelle 6.5: Relative Fehler zu Test 5.
6 Vergleich der Methoden 51
(a) Metallobjekt und ROI
(b) Gefilterte Ruckprojektion (c) Tikhonov-Regularisierung
(d) ART+TV (e) TV-Regularisierung
Abbildung 6.6: Ergebnisse von Test 5.
6 Vergleich der Methoden 52
Test 6
In den bisherigen Tests bestanden die Daten immer aus 180 Projektionen. Diese wollen
wir nun auf 360 erhohen und uns die entsprechenden Ergebnisse ansehen. Auch hier
fehlen 2.3% der Daten. Die Ergebnisse sind hier ahnlich zu denen aus Test 2, wobei
auch die Metallplatzierung die gleiche ist. Bei der gefilterten Ruckprojektion fallt auf,
dass das Rauschen verschwunden ist, daher sind die Metall-Artefakte noch deutlicher
zu erkennen.
ε2 εROI2
Lineare Interpolation und
gefilterte Ruckprojektion0.0921 0.4599
Tikhonov-Regularisierung 0.1033 0.1693
ART+TV 0.0157 0.2559
TV-Regularisierung 0.0035 0.1517
Tabelle 6.6: Relative Fehler zu Test 6.
6 Vergleich der Methoden 53
(a) Metallobjekt und ROI
(b) Gefilterte Ruckprojektion (c) Tikhonov-Regularisierung
(d) ART+TV (e) TV-Regularisierung
Abbildung 6.7: Ergebnisse von Test 6.
6 Vergleich der Methoden 54
Test 7
Der letzte Test zeigt auf, dass Metallobjekte im Bild trotz relativ guter Rekonstrukti-
onsverfahren ein Problem bleiben. Das Metallobjekt im Zentrum des Bildes wirkt sich
bis in den oberen Bereich des Bildes aus. Dies liegt vor allem darin begrundet, dass
sich die Kontraste, welche nicht rekonstruiert werden konnen, mit dem Metallobjekt
auf einer Linie befinden. Auf diese Weise ist der Teil der Projektion, in welchem der
Kontrast enthalten ist, vom Metall abgedeckt. Dadurch entsteht ein qualitativer Fehler
im Bild. Insgesamt fehlen 6% der Daten.
ε2 εROI2
Lineare Interpolation und
gefilterte Ruckprojektion0.0151 0.0529
Tikhonov-Regularisierung 0.0181 0.0520
ART+TV 6.1789 · 10−3 0.0470
TV-Regularisierung 4.5048 · 10−3 0.0383
Tabelle 6.7: Relative Fehler zu Test 7.
6 Vergleich der Methoden 55
(a) Ausgangsbild (b) Metallobjekt und ROI
(c) Gefilterte Ruckprojektion (d) Tikhonov-Regularisierung
(e) ART+TV (f) TV-Regularisierung
Abbildung 6.8: Ergebnisse von Test 7.
56
7 Fazit und Ausblick
Wir haben uns in dieser Arbeit mit der Thematik der Metall-Artefakt Reduktion aus-
einandergesetzt. Dabei haben wir zunachst die Funktionsweise eines Computertomo-
graphen studiert, um die Ursachen der Metall-Artefakte zu verstehen. Im Folgenden
haben wir uns die meist genutzten Methoden zur Rekonstruktion erarbeitet und diese
auf das Problem der Metall-Artefakte angepasst. Weiterhin wurde ein Verfahren herge-
leitet, welches die Rekonstruktion und die TV-Glattung simultan herstellt. Im Rahmen
eines Vergleichs konnten wir die Unterschiede in den Ergebnissen der Methoden her-
ausstellen.
Die Ergebnisse der TV-Regularisierung weisen die kleinsten Fehler in Bezug auf das
Ausgangsbild auf. Auch aus der optischen Begutachtung der resultierenden Bilder geht
hervor, dass die TV-basierten Methoden die besseren Ergebnisse liefern. Weiterhin hat
die TV-Regularisierung sich als robust gegenuber mehrerer Metallobjekte erwiesen.
Inwiefern die Methode allerdings Anwendung in der Praxis finden kann bleibt offen,
da die Berechnungszeiten dieser Methode mit Abstand die langsten sind. Des Weite-
ren konnten wir die Methoden nicht an echten Daten testen, was uns noch weitere
Moglichkeiten zur Beurteilung der Methoden gegeben hatte.
Eine wichtige Erkenntnis daruber hinaus ist, dass die auf die Problematik angepassten
Methoden viel bessere Ergebnisse liefern als etwa die gefilterte Ruckprojektion ohne
vorherige Aufbereitung der Daten. Das heißt, dass die Wahl der Methode zur Metall-
Artefakt Reduktion eine untergeordnete Rolle spielt – viel wichtiger ist es, dass eine
Methode angewandt wird.
Weiterhin konnten wir feststellen, dass Details, die in der Nahe des Metall-Objekts lie-
gen, mit keiner Methode wirklich sichtbar gemacht werden konnen. Dies liegt mutmaß-
lich in dem Fehlen vieler Daten uber die Details begrundet. Allerdings ist es moglich,
die strahlenformigen Artefakte, welche sich uber das gesamte Bild erstreckten, zu elimi-
nieren. Dafur eignen sich die bildbasierten Algorithmen deutlich besser als die lineare
Interpolation der gestorten Daten, wie wir es dem Test 2 entnommen haben. Der letzte
Test hat uns gezeigt, dass sich in allen Ergebnissen das Metall auch auf Regionen, die
7 Fazit und Ausblick 57
weit vom Metall entfernt liegen, auswirken und dort qualitative Unterschiede ausma-
chen konnen.
Um das Verfahren der TV-Regularisierung in der Praxis einsetzen zu konnen, mussten
weitere Tests mit echten Daten aus einem Computertomographen durchgefuhrt werden.
Insgesamt ist es uns jedoch gelungen, die Bildqualitat gegenuber der des Standardver-
fahrens deutlich zu erhohen.
58
Abbildungsverzeichnis
2.1 Shepp-Logan Phantom und dessen Sinogramm . . . . . . . . . . . . . . 4
2.2 Metall-Artefakte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
4.1 Tikhonov-Regularisierungen . . . . . . . . . . . . . . . . . . . . . . . . 20
4.2 ART-TV Rekonstruktionen . . . . . . . . . . . . . . . . . . . . . . . . 31
4.3 TV-Regularisierungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.4 Hohenlinien von TV-Regularisierungen . . . . . . . . . . . . . . . . . . 34
5.1 Lineare Interpolation und gefilterte Ruckprojektion . . . . . . . . . . . 36
5.2 Rekonstrunktion mit ART . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.3 Tikhonov-Regularisierung . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.4 TV-Rekonstruktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
6.1 Intensitaten des Shepp-Logan Phantoms . . . . . . . . . . . . . . . . . 40
6.2 Ergebnisse von Test 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
6.3 Ergebnisse von Test 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
6.4 Ergebnisse von Test 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
6.5 Ergebnisse von Test 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
6.6 Ergebnisse von Test 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
6.7 Ergebnisse von Test 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.8 Ergebnisse von Test 7. . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
59
Tabellenverzeichnis
6.1 Relative Fehler zu Test 1. . . . . . . . . . . . . . . . . . . . . . . . . . 42
6.2 Relative Fehler zu Test 2. . . . . . . . . . . . . . . . . . . . . . . . . . 44
6.3 Relative Fehler zu Test 3. . . . . . . . . . . . . . . . . . . . . . . . . . 46
6.4 Relative Fehler zu Test 4. . . . . . . . . . . . . . . . . . . . . . . . . . 48
6.5 Relative Fehler zu Test 5. . . . . . . . . . . . . . . . . . . . . . . . . . 50
6.6 Relative Fehler zu Test 6. . . . . . . . . . . . . . . . . . . . . . . . . . 52
6.7 Relative Fehler zu Test 7. . . . . . . . . . . . . . . . . . . . . . . . . . 54
60
Literaturverzeichnis
[Bjo96] Bjorck Ake: Numerical Methods for Least Squars Problems. SIAM,
1996
[BL10] Bredies, K. ; Lorenz, D.: Mathematische Bildverarbeitung: Einfuhrung
in Grundlagen und moderne Theorie. Vieweg+Teubner Verlag,
2010 http://books.google.de/books?id=0FVBkXCK3q8C. – ISBN
9783834810373
[BO13] Burger, Martin ; Osher, Standley: A Guide to the TV Zoo. In: Bur-
ger, Martin (Hrsg.) ;Osher, Standley (Hrsg.): Level Set and PDE Based
Reconstruction Methods in Imaging. Springer, 2013, S. 1–70
[CP11] Chambolle, Antonin ; Pock, Thomas: A First-Order Primal-Dual Al-
gorithm for Convex Problems with Applications to Imaging. In: Journal
of Mathematical Imaging and Vision 40 (2011), Nr. 1, 120-145. http://
dx.doi.org/10.1007/s10851-010-0251-1. – DOI 10.1007/s10851–010–
0251–1. – ISSN 0924–9907
[DMND+98] De Man, B. ; Nuyts, J. ; Dupont, P. ; Marchal, G. ; Suetens,
P.: Metal streak artifacts in X-ray computed tomography: A simulation
study. In: Nuclear Science Symposium, 1998. Conference Record. 1998
IEEE Bd. 3, 1998. – ISSN 1082–3654, S. 1860–1865 vol.3
[DV97] Dobson, D. ; Vogel, C.: Convergence of an Iterative Method for Total
Variation Denoising. In: SIAM Journal on Numerical Analysis 34 (1997),
Nr. 5, 1779-1791. http://dx.doi.org/10.1137/S003614299528701X. –
DOI 10.1137/S003614299528701X
[Kos12] Koskela, Olli: Three Ideas of Metal Artifact Reduction in CT Ima-
ging - Comparison Using Filtered Backprojection, Tikhonov regularization
Literaturverzeichnis 61
and Total Variation Regularisation Reconstruction Schemes, University of
Helsinki, Master’s Thesis, 2012
[MS12] Mueller, Jennifer L. ; Siltanen, Samuli: Linear and Nonlinear Inverse
Problems with Practical Applications. SIAM, 2012
[NW01] Natterer, Frank ;Wubbeling, Frank: Mathematical Methods in Image
Reconstruction. SIAM, 2001 (SIAM Monographs on Mathematical Mo-
deling and Computation)
[PKBK09] Prell, Daniel ; Kyriakou, Yiannis ; Beister, Marcel ; Kalender,
Willi A.: A novel forward projection-based metal artifact reduction me-
thod for flat-detector computed tomography. In: Physics in Medicine and
Biology 54 (2009), Nr. 21, 6575. http://stacks.iop.org/0031-9155/
54/i=21/a=009
[PW12] Pietschmann, Jan-F. ; Wubbeling, Frank: Inverse und schlecht
gestellte Probleme. http://wwwmath.uni-muenster.de/num/
Vorlesungen/InvProb_WS11/SkriptWS1112.pdf. Version: 2012
[ROF92] Rudin, Leonid I. ; Osher, Stanley ; Fatemi, Emad: Nonline-
ar total variation based noise removal algorithms. In: Physica D:
Nonlinear Phenomena 60 (1992), Nr. 1-4, 259 - 268. http://dx.
doi.org/http://dx.doi.org/10.1016/0167-2789(92)90242-F. – DOI
http://dx.doi.org/10.1016/0167–2789(92)90242–F. – ISSN 0167–2789
[SB14] Schiffer, Clemens ; Bredies, Kristian: Sinogram constrained TV-
minimization for metal artifact reduction in CT. 2014 (arXiv:1404.6691.
OAGM-2014-13). – Forschungsbericht. – Comments: Part of the OAGM
2014 proceedings (arXiv:1404.3538)
[SJP12] Sidky, Emil Y. ; Jørgensen, Jakob H. ; Pan, Xiaochuan: Convex
optimization problem prototyping for image reconstruction in computed
tomography with the Chambolle–Pock algorithm. In: Physics in Medicine
and Biology 57 (2012), Nr. 10, 3065. http://stacks.iop.org/0031-
9155/57/i=10/a=3065
Literaturverzeichnis 62
[SP08] Sidky, Emil Y. ; Pan, Xiaochuan: Image reconstruction in circular cone-
beam computed tomography by constrained, total-variation minimizati-
on. In: Physics in Medicine and Biology 53 (2008), Nr. 17, 4777. http://
stacks.iop.org/0031-9155/53/i=17/a=021
63
Plagiatserklarung des Studierenden
Hiermit versichere ich, dass die vorliegende Arbeit uber
Bildbasierte Algorithmen zur Metall-Artefakt Reduktion in der Computer-
tomographie
selbststandig verfasst worden ist, dass keine anderen Quellen und Hilfsmittel als die an-
gegebenen benutzt worden sind und dass die Stellen der Arbeit, die anderen Werken —
auch elektronischen Medien — dem Wortlaut oder Sinn nach entnommen wurden, auf
jeden Fall unter Angabe der Quelle als Entlehnung kenntlich gemacht worden sind.
(Datum, Unterschrift)
Ich erklare mich mit einem Abgleich der Arbeit mit anderen Texten zwecks Auffindung
von Ubereinstimmungen sowie mit einer zu diesem Zweck vorzunehmenden Speicherung
der Arbeit in eine Datenbank einverstanden.
(Datum, Unterschrift)
top related