SPECT-Dampfungskorrektur
ohne zusatzliche Messungen
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:
Carolin Roßmanith
Munster, November 2011
i
Zusammenfassung
Die Bildgebung mittels SPECT-Messungen spielt in der modernen Medizin eine im-
mer großer werdende Rolle. Dabei besteht ein Bedarf an Algorithmen, die moglichst
detailgenaue Ergebnisse liefern. Zu diesem Zweck ist bei der SPECT-Rekonstruktion
nicht nur die Verteilungsfunktion von Interesse, sondern auch die Abschwachung der
Gamma-Strahlung auf dem Weg durch das betrachtete Gewebe. Um diese mit moglichst
wenig zusatzlichem Aufwand miteinzubeziehen, gibt es viele verschiedene Ansatze.
Diese Arbeit beschaftigt sich mit den Grundlagen der SPECT-Rekonstruktion mit be-
sonderem Blick auf die verschiedenen Methoden zur Bestimmung der Abschwachung.
Zunachst wird ein Uberblick uber die Problemstellung gegeben und ein mathematisches
Modell fur SPECT vorgestellt. Die Grundzuge der gebrauchlichsten Losungsalgorithmen,
die gefilterte Ruckprojektion und der EM-Algorithmus, werden ebenfalls erlautert.
Anschließend folgen zwei unterschiedliche Modelle zur Ermittlung der auftretenden
Dampfung, die keine zusatzlichen Messungen wahrend der eigentlichen SPECT-Messung
benotigen: Das erste Modell baut auf einer affin-linearen Transformation eines Proto-
typs fur die Abschwachung auf, wahrend beim zweiten eine nichtlineare Transformation
verwendet wird. Zum Abschluss der Arbeit wird eine kurze Auswertung und ein Ver-
gleich der Methoden vorgenommen.
ii
Eidesstattliche Erklarung
Hiermit versichere ich, Carolin Roßmanith, dass ich die vorliegende Arbeit selbststandig
verfasst und keine anderen als die angegebenen Quellen und Hilfsmittel verwendet
habe. Gedanklich, inhaltlich oder wortlich ubernommenes habe ich durch Angabe von
Herkunft und Text oder Anmerkung belegt bzw. kenntlich gemacht. Dies gilt in gleicher
Weise fur Bilder, Tabellen, Zeichnungen und Skizzen, die nicht von mir selbst erstellt
wurden.
Munster, 25. Juni 2010
Carolin Roßmanith
iii
Inhaltsverzeichnis
1 Was ist Emissionstomographie? 1
1.1 Einleitung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Mathematisches Modell fur SPECT . . . . . . . . . . . . . . . . . . . . 2
1.3 Bildbeeinflussende Faktoren . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 8
2.1 Einleitung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Die gefilterte Ruckprojektion . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Der EM-Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3.1 Diskretisierung des Modells . . . . . . . . . . . . . . . . . . . . 13
2.3.2 Statistische Grundlagen . . . . . . . . . . . . . . . . . . . . . . 16
2.3.3 Der ML-EM-Algorithmus . . . . . . . . . . . . . . . . . . . . . . 18
2.3.4 Maximierung der Log-Likelihood-Funktion . . . . . . . . . . . . 19
3 Dampfungskorrektur 22
3.1 Problemstellungen und Losungsansatze . . . . . . . . . . . . . . . . . . 22
3.2 Bestimmung mittels affin-linearer Transformation eines vorhandenen Pro-
totyps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.1 Grundidee und Rekonstruktion . . . . . . . . . . . . . . . . . . 24
3.2.2 Beweis der Konsistenzbedingung . . . . . . . . . . . . . . . . . . 27
3.2.3 Numerische Ergebnisse . . . . . . . . . . . . . . . . . . . . . . . 35
3.3 Bestimmung mittels nichtlinearer Transformation . . . . . . . . . . . . 36
3.3.1 Grundidee und Rekonstruktion . . . . . . . . . . . . . . . . . . 36
3.3.2 Der Algorithmus . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.3.3 Numerische Ergebnisse . . . . . . . . . . . . . . . . . . . . . . . 42
3.4 Vergleich und Auswertung . . . . . . . . . . . . . . . . . . . . . . . . . 42
Abbildungsverzeichnis 45
Literaturverzeichnis 46
1
1 Was ist Emissionstomographie?
1.1 Einleitung
Als Emissions-Tomographie (ET) bezeichnet man eine Form der medizinischen Bild-
gebung, die sich die Eigenschaften radioaktiver Substanzen zunutze macht. Die zwei
Haupttechniken, Positron-Emission-Tomography (PET) und Single-Photon-Emission-
Computed-Tomography (SPECT), finden in der Praxis weitlaufige Anwendung. Sie
dienen dazu, die Verteilung eines radioaktiven Stoffes im Korper des Patienten zu vi-
sualisieren. Dies ermoglicht beispielsweise die Erkennung von Tumorzellen oder der
Wirkung von Medikamenten auf bestimmte Hirnregionen.
Die Grundidee der Emissionstomographie lasst sich in wenigen Worten beschreiben:
Zunachst wird dem Organismus eine radioaktiv angereicherte Substanz, ein sogenann-
tes Radiopharmazeutikum, injiziert, das sich dann an bestimmten Stellen im Korper
ablagert. Von dieser Position aus werden Photonen emittiert, die von außerhalb detek-
tiert werden konnen. Um beispielsweise Tumorzellen zu erkennen, wird dem Korper eine
radioaktiv angereicherte Zuckerart ubergeben. Da Tumorzellen viel Energie (in Form
von Zucker) benotigen, wird dieser dorthin transportiert und lagert sich an den Zellen
ab. Von dort aus werden dann Photonen emittiert. Um die Zahl der emittierten Pho-
tonen zu messen, befinden sich außerhalb des Patienten eine oder mehrere sogenannte
Gamma-Kameras, die um den Korper herum rotieren und Messungen aus unterschied-
lichen Raumrichtungen vornehmen. Aus diesen Messungen wird schließlich mit Hilfe
mathematischer Methoden ein zwei- oder dreidimensionales Bild der entsprechenden
Korperregion rekonstruiert, das medizinisch ausgewertet werden kann.
Der Unterschied zwischen PET und SPECT besteht im Wesentlichen in der Art der ver-
wendeten radioaktiven Substanz. PET verwendet sogenannte Positronen-Emitter: Bei
jedem radioaktiven Zerfall eines Teilchens wird ein Positron aus dem Kern emittiert.
Beim Auftreffen auf ein Elektron tritt Annihilation ein: Beide Teilchen verwandeln sich
1 Was ist Emissionstomographie? 2
je in ein Photon, diese werden dann in entgegengesetzte Richtung emittiert. Bei SPECT
werden dagegen Radiopharmazeutika verwendet, bei denen bei jedem radioaktiven Zer-
fall jeweils nur ein Photon emittiert wird. Diese Photonen, die von der Objektregion
ausgesandt werden, treffen auf Gamma-Kameras außerhalb des Korpers. In diesen be-
finden sich sogenannte Kollimatoren. Sie dienen dazu, einen Großteil der emittierten
Strahlung abzublocken und nur Photonen”durchzulassen“, die in eine bestimmte Rich-
tung abgestrahlt werden. In der Praxis werden großtenteils Parallel-Hole-Kollimatoren
verwendet. Sie bestehen aus parallel angeordneten Kanalen, die nur Photonen mit
gleicher Strahlrichtung wie die Ausrichtung der Kanale berucksichtigen. Eine graphi-
sche Darstellung der Parallel-Hole-Kollimatoren ist in Bild 1.1 zu sehen. Mit diesen
entstehen Messungen, die als Parallelprojektionen entlang einer bestimmten Achse be-
schrieben werden konnen. Um ein vollstandiges Bild des Objektes zu erhalten, behalten
die Kollimatoren nicht eine Position bei, sondern bewegen sich um das Objekt herum
und messen Photonen aus verschiedenen Winkeln. Hierbei gilt bis zu einer bestimm-
ten Grenze: Je mehr Positionen die Kollimatoren einnehmen, desto genauer wird das
entstandene Bild, jedoch ist dies auch mit mehr Aufwand und langerer Dauer des Vor-
gangs verbunden. Ab einer Grenze verbessert sich die Genauigkeit des Bildes jedoch
nicht mehr [18].
1.2 Mathematisches Modell fur SPECT
Das mathematische Modell fur die SPECT-Rekonstruktion richtet sich weitestgehend
nach [18], [12] und [15]. Um die mathematische Funktionsweise bei der ET deutlich zu
machen, beschranken wir uns in dieser Arbeit ausschließlich auf den zweidimensionalen
Fall (Details zum dreidimensionalen Modell sind in [18] nachzulesen). Die Punkte im
zweidimensionalen betrachteten Bereich identifizieren wir mit Koordinaten (x, y). Die
Verteilung der radioaktiven Substanz in jedem Punkt ist dann gegeben durch f (x, y).
Der zweidimensionale Bereich wird aus verschiedenen Blickwinkeln betrachtet, indem
die Koordinatenachsen um einen Winkel ϕ gedreht werden. Die neuen, mit Winkel
ϕ entstandenen Koordinaten bezeichnen wir mit (xr, yr). Der betrachtete Querschnitt
wird in Abbildung 1.2 dargestellt.
Die Messungen der Emissionen aus einem bestimmten Winkel ϕ sind dann, im mathe-
1 Was ist Emissionstomographie? 3
Abbildung 1.1: Photonenemission aus einem dreidimensionalen Objekt mit umgeben-den Kollimatoren, aus [18].
Nur einige passieren die Kollimatorkanale, wahrend andere sie verfehlen oderblockiert werden.
Abbildung 1.2: Zweidimensionales Objekt im Koordinatensystem, aus [18]
1 Was ist Emissionstomographie? 4
matischen Sinne, anschaulich in etwa die Linienintegrale uber die Verteilung f entlang
einer zur yr-Achse parallelen Linie L, also die Parallelprojektionen auf die xr-Achse
[18]:
P (xr, ϕ) =
∫L
f (x) dyr
Wir wollen dies noch klarer formulieren. Dazu benotigen wir die folgenden Definitionen
[15].
Definition und Satz: Schwartzraum
Der Schwartzsche Raum oder Raum der schnell fallenden Funktionen S (R2) auf R2 ist
definiert durch
S(R2)
:= f ∈ C∞(R2)| supx∈R2
∣∣xkDlf (x)∣∣ <∞ ∀ k, l ∈ Z2
+
Der Schwartzraum ist mit der Verkettung von Abbildungen ein Vektorraum. Fur 1 ≤p <∞ liegt S (Rn) dicht in Lp (Rn) bezuglich der Lp-Norm.
Definition: Radontransformation
Sei θ ∈ S1, der Einheitssphare im R2, s ∈ R und f ∈ S (R2). Dann ist die Radontrans-
formation Rf auf S × R definiert durch
Rf (θ, s) :=
∫x·θ=s
f (x) dx
Man integriert hierbei uber die Linie, die senkrecht auf θ steht und Abstand s zum Ko-
ordinatenursprung hat. Die oben genannten Parallelprojektionen lassen sich auf diese
Weise durch die Radontransformation ausdrucken. Gegeben die Radontransformation
fur verschiedene Werte θ und s soll nun die Verteilungsfunktion f aus diesen Daten
bestimmt werden. Wir stehen hier vor einem schlecht gestellten inversen Problem. Ein
inverses Problem heißt gut gestellt, wenn folgende drei Bedingungen erfullt sind:
1. Existenz: Das Problem besitzt eine Losung.
2. Eindeutigkeit: Die Losung ist eindeutig.
3. Stabilitat: Die Losung hangt stetig von den Eingabedaten ab.
1 Was ist Emissionstomographie? 5
Das Problem heißt schlecht gestellt, wenn mindestens eine der Bedingungen verletzt
ist. In der Theorie besitzt das Problem der Inversion der Radontransformation nicht
immer eine Losung, daher sprechen wir von einem schlecht gestellten Problem.
Zu beachten ist hierbei, dass sich die Messungen in der Praxis nicht exakt als die
Parallelprojektionen ergeben. Die Starke der emittierten Gamma-Strahlen wird auf
dem Weg bis zu den Kollimatoren abgeschwacht, die Intensitat der Abschwachung va-
riiert, je nachdem welche Art von Gewebe die Strahlung passiert. Aus diesem Grund
muss in der Berechnung der Verteilung f noch eine zusatzliche Dampfungsfunktion µ
miteinbezogen werden. Wir betrachten dann nicht mehr die gewohnliche, sondern die
gedampfte Radontransformation [15]:
Definition: Gedampfte Radontransformation
Seien f, µ ∈ S (R2). Dann ist die gedampfte Radontransformation Rµf auf S × Rdefiniert durch
Rµf (θ, s) =
∫xθ=s
f (x) e−∫∞0 µ(x+tθ⊥)dtdx
mit θ⊥ =( −sinϕcosϕ
). µ ist die sogenannte Dampfungsfunktion.
Die Abschwachung der Gammastrahlung wird durch den Exponentialterm, der die
Dampfungsfunktion µ enthalt, miteinbezogen. Wir konnen nun annehmen, dass die
Messungen des SPECT-Gerates von einer Funktion g reprasentiert werden. Es gilt
dann also:
Rµf (θ, s) = g (θ, s)
Naturlich sind diese Daten in der Praxis diskret, um das Problem jedoch mit Hilfe
analytischer Methoden zu losen, geht man zunachst von einem kontinuierlichen Modell
mit einer Funktion g aus. Man trifft hierbei die (modelltheoretisch sinnvolle) Annahme,
dass die Exponentialfunktion den zu berucksichtigenden Dampfungsterm reprasentiert.
Das zu losende Problem bei der SPECT-Rekonstruktion ergibt sich also als die Inver-
sion der gedampften Radontransformation.
Bei der PET-Rekonstruktion steht man vor einem etwas anderen Problem: Da bei je-
1 Was ist Emissionstomographie? 6
dem radioaktiven Zerfall zwei Photonen in genau entgegengesetzte Richtung emittiert
werden, wird die Dampfungsfunktion µ nicht nur uber einen Teil der Messlinie inte-
griert, sondern uber die gesamte Linie. Die Messdaten werden dann von einer Funktion
der Form
Iµ (θ, s) =
∫x·θ=s
f (x) e−∫y·θ=s µ(y)dydx
= e−∫y·θ=s µ(y)dy
∫x·θ=s
f (x) dx
reprasentiert [12].
Grundsatzlich ist man an der Bestimmung der Dampfungsfunktion µ nicht interes-
siert. Ist diese allerdings nicht vernachlassigbar klein, so ist es notwendig, auch sie
mit Hilfe eines Algorithmus zu ermitteln. Dazu gibt es in der Literatur eine Vielzahl
moglicher Ansatze, auf die in Kapitel 3 genauer eingegangen wird.
1.3 Bildbeeinflussende Faktoren
Bei SPECT gibt es neben der bereits erwahnten Dampfung noch einige weitere Fakto-
ren, welche das rekonstruierte Bild und dessen Qualitat beeinflussen und bei einer rea-
listischen Modellierung miteinbezogen werden sollten. Ein wahrend des Messvorgangs
auftretendes Problem ist Streuung, eine Storung der Messdaten durch eine Interaktion
zwischen den Gammastrahlen und dem betrachteten Objekt. Eine Streuung entsteht
beispielsweise durch eine Bewegung des Patienten wahrend der Messung und verursacht
eine Abweichung der angezeigten Strahlung. Dies wird im rekonstruierten Bild in Form
von Verwischung oder Unscharfe sichtbar. Desweiteren entsteht eine tiefenabhangige
Verwischung dadurch, dass Kollimatoren in der Praxis nicht nur Messungen uber eine
einzige Linie vornehmen konnen. Waren die Kanale der Kollimatoren unendlich dunn,
so wurde tatsachlich immer nur eine Linie, entlang der Photonen laufen, berucksichtigt.
So allerdings muss man in Kauf nehmen, dass eine Ungenauigkeit durch Messungen
uber mehrere Linien entsteht [18].
Ein weiterer Faktor von Bedeutung, der in [18] naher erlautert wird, ist die Storung
der Bildqualitat durch Bildrauschen. Ist I das Originalbild und X eine Zufallsvariable,
1 Was ist Emissionstomographie? 7
so ergibt sich das rekonstruierte Bild als [9]
IX = I +X
also als eine leicht abgewandelte Version des Originals. Diese Art von Bildrauschen wird
als additives Rauschen bezeichnet. Der Grund fur diese Abweichung liegt in der Detek-
tion der Gamma-Strahlen. Gamma-Strahlen sind Photonen und unterliegen den Geset-
zen der Quantenphysik und der zufalligen Zerfalle. Die Zahl der gemessenen Zerfalle in
einer bestimmten Region (z.B. in einem Pixel) in einem festen Zeitintervall ist zufallig:
Wurde man dieselbe Region immer wieder im gleichen Zeitintervall untersuchen, so
erhielte man jedes Mal ein etwas anderes Ergebnis, unter Umstanden also auch ein
von den anderen Messungen stark abweichendes. Diese Variation verursacht eine”ge-
sprenkelte“ Erscheinung im rekonstruierten Bild und wird als Bildrauschen bezeichnet.
Fur die Wahrscheinlichkeit dieses Zufallsvorgangs gibt es unterschiedliche Modelle. Die
gangigste Methode ist die Annahme der Poisson-Verteilung fur die Zufallsvariable X,
da physikalische Zerfallsprozesse im Allgemeinen dieser Verteilung unterliegen. Diese
Idee kommt in der Rekonstruktion mittels des EM-Algorithmus zur Anwendung.
8
2 Standard-Algorithmen bei der
SPECT-Rekonstruktion
2.1 Einleitung
Bei den Standard-Algorithmen zur Bestimmung der Verteilung f unterscheidet man
zwischen zwei Typen, den analytischen und den iterativen Methoden. Vertreter beider
Arten sind mittlerweile in der Praxis haufig zu finden.
Ziel der analytischen Methoden ist es, eine explizite Inversionsformel fur die Integral-
transformation anzugeben und eine diskrete Version dieser zu implementieren. Ty-
pischerweise werden hierbei beeinflussende Faktoren wie Rauschen missachtet, um die
mathematische Grundstruktur zu erhalten, welche die Angabe einer solchen Inversions-
formel erlaubt. Analytische Methoden ergeben demnach Losungen, die relativ einfach
und mit wenig Aufwand zu berechnen sind. Die meistvertretene Methode ist die gefil-
terte Ruckprojektion, auf die wir in Kapitel 2.2 genauer eingehen werden.
Aufgrund ihres hoheren numerischen Aufwands wurden iterative Methoden in der Ver-
gangenheit seltener praktisch angewandt. Da sich die Rechenleistung moderner Com-
puter in den letzten Jahren jedoch drastisch verbessert hat, erfreuen sie sich wachsender
Beliebtheit. Iterative Methoden verwenden eine diskrete Version der Integraltransfor-
mation, um mit Hilfe statistischer Schatzverfahren eine Losung zu konstruieren. Dies
erlaubt auch die Miteinbeziehung realistischer Faktoren wie Rauschen oder Streuung.
Der Hauptkonflikt zwischen analytischen und iterativen Verfahren besteht demnach in
der Entscheidung zwischen hoherer Genauigkeit und hoherer Effizienz. Eine bekannte
iterative Methode ist der EM-Algorithmus, der in Kapitel 2.3 vorgestellt werden soll
[18].
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 9
2.2 Die gefilterte Ruckprojektion
Eine analytische Methode mit haufiger Anwendung ist die gefilterte Ruckprojektion
(FBP). Sie basiert auf dem idealisierten Modell von SPECT, was bedeutet, dass meh-
rere realistische Faktoren außer Acht gelassen werden. Ein Vorteil der Methode ergibt
sich dabei aus der Einfachheit der Berechnung. Der folgende Abschnitt orientiert sich
großtenteils an [15], [12] und [14].
Die Funktionsweise des Algorithmus lasst sich im dreidimensionalen Fall anhand ei-
nes einfachen Beispiels veranschaulichen (dies ist auch auf den zweidimensionalen Fall
ubertragbar). Man stelle sich ein Aquarium vor, in dem sich an festen Positionen einige
Fische befinden. Nun positioniert man vor das Aquarium eine Lichtquelle und dahinter
eine Leinwand. Auf dieser sind dann die Schatten der Fische zu sehen, eine zweidimen-
sionale Abbildung der dreidimensionalen Objekte. Dreht man jetzt Lichtquelle und
Leinwand parallel zueinander um das Aquarium herum, erhalt man die zweidimensio-
nalen Bilder der Fische aus verschiedenen Blickwinkeln. Die Vorgehensweise ist dann,
die einzelnen Projektionen auf das dreidimensionale Aquarium zuruckzuprojizieren.
Da man jedoch keine Informationen daruber besitzt, in welcher Tiefe sich die Fische
befinden, werden die Projektionen orthogonal zur jeweiligen Position der Leinwand
”verschmiert“. Dadurch entsteht ein Fehler, der wiederum bei der Ruckprojektion mit
einem geeigneten Bildfilter unterdruckt werden muss.
Im Zweidimensionalen wird dieser Vorgang in Abbildung 2.1 deutlich.
Abbildung 2.1: Parallel- und Ruckprojektion, aus [18](a) Parallelprojektion entlang der zur yr-Achse parallelen Linien, (b)Ruckverschiebung einer einzelnen Projektion in den Objektbereich
Die Kollimatoren bewegen sich um das Objekt und nehmen Messungen aus verschie-
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 10
Abbildung 2.2: (a) Ursprungliches zu rekonstruierendes Bild, (b) Aus Projektionen ent-standenes Sinogramm, aus [18]
denen Blickwinkeln ϕ vor. Diese lassen sich als Projektionen auf die jeweilige Koordi-
natenachse xr veranschaulichen. Diese Projektionen erzeugen wiederum ein neues Bild
mit Koordinaten (xr, ϕ), ein sogenanntes Sinogramm (s. Abb. 2.2).
Durch”Zuruckschieben“ der Projektionen in die Objektregion entsteht dann eine An-
naherung des ursprunglichen Bildes, wobei auch hier zu beachten ist, dass beeinflus-
sende Faktoren die Bildqualitat verschlechtern. Aus diesem Grund ist die Einbringung
eines geeigneten Filters notwendig. Ungefilterte Ruckprojektion fuhrt zu keinem gu-
ten Ergebnis. Sie fuhrt beim rekonstruierten Bild zu einer starken”Verwischung“ des
Originalbildes. Um dies zu verdeutlichen, betrachten wir die Abbildung 2.2(a) als zu
rekonstruierendes Bild. Wir wollen nun einen Punkt rekonstruieren, der außerhalb der
beiden weißen Kreise liegt, also in dem die tatsachliche Verteilung des Radiotracers
gleich Null ist. Der rekonstruierte Wert ergibt sich als die Summe aller gemessenen
Verteilungen uber alle Linien, die durch diesen Punkt laufen. Je naher der Punkt sich
an einem der weißen Kreise befindet, desto hoher ist die Wahrscheinlichkeit, dass eine
der Messlinien durch diesen Kreis lauft (s. Abb. 2.3). Der uber diese Linie gemessene
Wert ist demnach nicht Null, wird aber dennoch zu unserem außerhalb liegenden Punkt
hinzugezahlt. Dadurch ergibt sich dort letztendlich ein von Null abweichender Wert,
obwohl die Verteilung im Originalbild exakt Null war. Im rekonstruierten Bild werden
also ohne eine Filterung keine scharfen Kanten erzielt, es entsteht eine”verschmierte“
Version des ursprunglichen Bildes.
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 11
Abbildung 2.3: Rekonstruktion eines Punktes mit tatsachlicher Verteilung 0, aus [18]
Die mathematische Idee der gefilterten Ruckprojektion ist es, eine explizite Inversions-
formel der gedampften Radontransformation zu finden und eine diskrete Version dieser
zu implementieren. Das mathematische Problem besteht also zunachst in der Inver-
sion der gedampften Radontransformation. Dieses Problem ist in der Theorie bereits
vollstandig gelost. Um den Vorgang der gefilterten Ruckprojektion zu verstehen, zeigen
wir diesen nur fur die ungedampfte Form
(Rf) (θ, s) =
∫xθ=s
f (x) dx
und losen das Problem im zweidimensionalen Raum. Wir benotigen zunachst die zu R
adjungierte Radontransformation.
Satz:
Sei S1 die Einheitssphare im R2, g : S1 × R→ R und R∗g auf R2 gegeben durch
(R∗g) (x) :=
∫S1
g (θ, x · θ) dθ
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 12
Abbildung 2.4: Ergebnis der gefilterten Ruckprojektion mit verschiedenen Anzahlenvon Blickwinkeln, aus [18]
R∗ ist die zu R adjungierte Radontransformation bezuglich des L2-Skalarproduktes.
Die adjungierte Radontransformation ist der sogenannte Ruckprojektionsoperator. An-
gewandt auf die Funktion g entspricht er der ungefilterten Ruckprojektion. Es handelt
sich also nicht um eine inverse Radontransformation, da wir bereits gesehen hatten,
dass ungefilterte Ruckprojektion nicht die tatsachliche Funktion f ergibt, sondern eine
gestorte Version von f. Das bedeutet [12]:
Satz:
Fur f ∈ S (R2) gilt:
R∗ (Rf) (x) = 2f (x) ∗ 1
|x|
Beweis: [12]
Allgemeiner besteht zwischen der Radontransformation und ihrer Adjungierten folgen-
der Zusammenhang:
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 13
Satz:
Seien f ∈ S (R2), υ ∈ S (S1 × R). Dann gilt:
(R∗υ) ∗ f (x) = R∗ (υ ∗Rf) (x)
Beweis: [15]
In der Theorie gibt es eine explizite Inversionsformel fur die Radontransformation, f
kann also explizit angegeben werden. Fur den Ubergang vom kontinuierlichen zum dis-
kreten Problem macht sich die gefilterte Ruckprojektion die Eigenschaften der Faltung
zunutze: Man wahlt υ so, dass R∗υ eine Approximation der Dirac-Funktion δ darstellt.
(R∗υ) ∗ f ist dann aufgrund der Faltungseigenschaft der Dirac-Funktion eine Approxi-
mation von f. Damit lost man die rechte Seite der Gleichung: Man”filtert“ zunachst
die Datenfunktion g mit υ und wendet anschließend die Ruckprojektion R∗ darauf an.
Daraus leitet sich der Name des Algorithmus ab, die gefilterte Ruckprojektion. Fur die
Approximation von δ gibt es verschiedene Moglichkeiten, demnach existieren verschie-
dene geeignete Filter, die unterschiedliche Ergebnisse liefern.
Kurz zusammengefasst enthalt der Algorithmus der gefilterten Ruckprojektion also
folgende Schritte:
1. Wahle einen geeigneten Filter υ mit R∗υ ≈ δ
2. Berechne υ ∗ g3. Berechne R∗ (υ ∗ g)
2.3 Der EM-Algorithmus
2.3.1 Diskretisierung des Modells
Wir wollen im Folgenden die Funktionsweise des EM-Algorithmus, einer iterativen Me-
thode, erlautern. Naheres findet man in [15], [18] und [17]. Der Expectation-Maximization-
Algorithmus (EM) ist ein iteratives Verfahren zur Berechnung der Verteilung f, das
sich statistische Methoden zunutze macht. Zur Erinnerung: Wir wollen die gedampfte
Radontransformation (auf diskrete Weise) invertieren. Um die Vorgehensweise zu ver-
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 14
stehen, betrachten wir zunachst die ungedampfte Radontransformation
(Rf) (θ, s) =
∫x·θ=s
f (x) dx
Wir nehmen an, dass der Trager von f im Einheitskreis des R2 liegt, das heißt wir
betrachten als zweidimensionalen Querschnitt des Objektes den Einheitskreis. Diesen
uberdecken wir nun mit endlich vielen Quadraten (Pixeln) P1, . . . , Pn und setzen vor-
aus, dass f in jedem dieser Pixel einen konstanten Wert annimmt. Gesucht ist also der
Vektor F = (f1, . . . , fn)t, wobei fj den Wert von f in Pixel Pj darstellt. Desweiteren
betrachten wir Geraden L1, . . . , Lm, entlang denen von den Kollimatoren Strahlungs-
intensitaten g1, . . . , gm gemessen werden (s. Abb. 2.5). Es ergibt sich also fur jedes
i = 1, . . . ,m:
gi =
∫Li
f (x) dx
Wir wollen nun eine diskrete Version dieses Integrals bestimmen. Zu diesem Zweck defi-
nieren wir aij als die Lange des Schnitts der Linie Li mit Pixel Pj. Eine Approximation
des Integrals ergibt sich damit als Summe der Intensitaten fj, jeweils multipliziert mit
der Lange des Schnitts der Gerade mit dem entsprechenden Pixel. Das bedeutet:
(ai1, . . . , ain)F = gi ⇔ ai1f1 + . . .+ ainfn = gi
Daraus ergibt sich direkt, durch Untereinanderschreiben der einzelnen Zeilen, das li-
neare Gleichungssystem
a11 . . . a1n...
...
am1 . . . amn
f1...
fn
=
g1...
gm
F ist also die gesuchte Losung unseres Problems.
Wir wollen nun zur gedampften Radontransformation ubergehen, die die von einem
SPECT-Gerat gemessenen Daten beschreibt. Dazu mussen wir nur eine Abanderung
der Matrixelemente der zuvor gegebenen Matrix vornehmen. Wir gehen hier davon aus,
dass die Dampfungsfunktion µ bereits bekannt ist. Die gemessenen Intensitaten gi sind
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 15
Abbildung 2.5: Diskretisiertes Modell der Objektregion aus [18]Bildbereich bestehend aus n Pixeln und Kollimatoren mit m Kanalen
gegeben durch:
gi =
∫Li
f (x) e−
∫Liµ(y)dy
dx
Li steht fur den Abschnitt von Linie Li, der zwischen Pixel Pj und Detektor i liegt,
also genau das Wegstuck, in dem Dampfung auftritt. Analoges Vorgehen zum un-
gedampften Fall fuhrt auf ein ahnliches Gleichungssystem: Die fj werden nun nicht
nur mit den Langen der jeweiligen Schnitte multipliziert, sondern zusatzlich noch mit
dem Damfungsterm. Es entsteht also eine Anderung in den Matrixwerten. Fur die
neuen Werte definieren wir:
aij := aije−
∫Liµ(y)dy
und losen das Gleichungssystem AF = g. Im Folgenden schreiben wir der Einfachheit
halber weiter aij statt aij.
Bei der Rekonstruktion von PET-Bildern mit Hilfe des EM-Algorithmus werden die-
selben Methoden angewandt. Der Unterschied besteht lediglich in der Art, in der
Dampfung der Gammastrahlung auftritt und demnach in den Werten der konstruierten
Matrix. Wir konnen fur SPECT-Rekonstruktionen also die Ideen von PET weiterver-
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 16
wenden.
Die Losbarkeit des aus obiger Konstruktion resultierenden Gleichungssystems ist abhan-
gig von der Anzahl der Gleichungen und Variablen. Die Zahl der Gleichungen ist durch
die Anzahl der Messwerte fest vorgeschrieben, bei der Anzahl der Pixel und somit der
Variablen hat man dagegen freie Wahl. Die einfachste Idee ist daher, die Pixelzahl
n kleiner oder gleich der Anzahl der Messwerte zu wahlen. Eine eindeutige Losung
kann dann mit Hilfe der Kleinste-Quadrate-Methode konstruiert werden. Der EM-
Algorithmus geht von einem anderen Ansatz aus: Hier wird die Pixelanzahl deutlich
großer gewahlt als die Zahl der Gleichungen, das Gleichungssystem ist daher unterbe-
stimmt und es existieren in der Regel unendlich viele Losungen. Um eine Losung zu
finden, ist also ein Verfahren notwendig, welches aus allen moglichen Losungen dieje-
nige auswahlt, die am besten zu den gegebenen Daten passt. Dies geschieht mit Hilfe
eines statistischen Schatzverfahrens, des EM-Algorithmus.
2.3.2 Statistische Grundlagen
Wir definieren eine Familie von Zufallsvariablen Xij, die die Anzahl der Photonen re-
prasentiert, die aus Pixel Pj stammen und auf Linie Li gemessen werden. Man kann
annehmen, dass diese Zufallsvariablen unabhangig und Poissonverteilt sind. Als Erin-
nerung:
Definition und Satz: Poissonverteilte Zufallsvariable
Eine Zufallsvariable Z heißt Poissonverteilt zum Parameter λ, wenn gilt:
P (Z = k) = e−λλk
k!
Fur Poissonverteilte Zufallsvariablen gilt:
E (Z) = V ar (Z) = λ
Die Voraussetzung der Poissonverteilung ist eine sinnvolle und in der Praxis haufig
angewandte Vorgehensweise und grundet sich darin, dass die Anzahl der emittierten
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 17
Photonen bei einem Zerfallsprozess dieser Verteilung unterliegt. Es sind allerdings auch
andere Modelle fur die Verteilung denkbar, wie beispielsweise ein Modell mit Normal-
verteilung [18].
Weiter definieren wir Zufallsvariablen γi, ϕj durch:
γi :=n∑j=1
Xij und ϕj :=m∑i=1
Xij
γi steht fur die Anzahl aller auf Linie Li gemessenen Photonen, ϕj misst die Gesamtzahl
der aus Pixel Pj stammenden Photonen. Mit den Xij sind demnach auch die γi und
ϕj unabhangige Poissonverteilte Zufallsvariablen. Wir werden nun die gi nicht mehr
als Integrale uber die Verteilungsfunktion f mit Dampfungsterm betrachten, sondern
als Realisierungen einer Poissonverteilten Zufallsvariable, deren Mittelwert das Integral
darstellt.
Sei aij nun die bedingte Wahrscheinlichkeit, dass ein Photon auf Linie Li gemessen
wird unter der Voraussetzung, dass es aus Pj stammt. Es muss an dieser Stelle gelten:
m∑i=1
aij ≤ 1
Man kann zur Vereinfachung annehmen, dass ein Photon sicher auf irgendeiner Linie
gemessen wird, also die Spaltensummen der Matrix normiert sind:
m∑i=1
aij = 1
Mit diesen Voraussetzungen ergibt sich, in der Sprache der Wahrscheinlichkeitstheorie,
dass die gi Realisierungen der Zufallsvariablen γi darstellen. Wir haben also eine Stich-
probe g1, . . . , gm der Zufallsvariablen γ1, . . . , γm vorliegen.
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 18
2.3.3 Der ML-EM-Algorithmus
Um nun den Wert der Verteilungsfunktion f in jedem Pixel naherungsweise zu bestim-
men, der am besten zu der Stichprobe g1, . . . , gm passt, verwenden wir die Maximum-
Likelihood-Methode, die auch zu dem Begriff Maximum-Likelihood-Expectation-Maxi-
mization (ML-EM) fuhrt. Diese sucht sich unter allen moglichen Verteilungen fur f
diejenige heraus, fur die die gezogene Stichprobe g unter allen moglichen Stichproben
die großte Wahrscheinlichkeit hat, das heißt genau diese Stichprobe g zu ziehen ist
wahrscheinlicher als jede andere. In Zeichen:
P (γ1 = g1, . . . , γm = gm) = maxh
P (γ1 = h1, . . . , γm = hm)
h durchlauft alle moglichen Daten, die das Gerat messen konnte. Wir berechnen zunachst
die Wahrscheinlichkeit der einzelnen γi fur i = 1, . . . ,m. Da diese aus den Summen uber
die Zufallsvariablen Xij bestehen und somit ebenso der Poissonverteilung unterliegen,
ergibt sich
P (γi = gi) = P
(n∑j=1
Xij = gi
)= e−λ
λgi
gi!
mit einem unbekannten Parameter λ. Fur Poissonverteilte Zufallsvariablen gilt, dass
ihr Erwartungswert genau ihrem Parameter entspricht. Wir benotigen also den Er-
wartungswert der Variablen. γi misst die Anzahl aller Photonen, die auf einer Linie
Li registriert werden. Demnach erhalt man als erwartete Anzahl die Summe uber alle
Werte fj mal die Wahrscheinlichkeit, dass ein Photon aus eben diesem Pixel Pj stammt
und auf Li gemessen wird, also:
λ = E (γi) =n∑j=1
aijfj = (Af)i
Damit sind wir in der Lage, die Wahrscheinlichkeit unserer gezogenen Stichprobe zu
bestimmen:
P (γ1 = g1, . . . , γm = gm) = P (γ1 = g1) . . . P (γm = gm) =m∏i=1
P (γi = gi)
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 19
=m∏i=1
(Af)giigi!
e−(Af)i := L (f)
wobei wir uns die Unabhangigkeit der γi zunutze machen. Dies ist also die Wahrschein-
lichkeit, die wir maximieren wollen. Wir definieren diese Funktion als die Maximum-
Likelihood-Funktion L (f) und bilden, zur Vereinfachung des Maximierungsvorgangs,
die Log-Likelihood-Funktion durch Logarithmieren:
` (f) := logL (f) =m∑i=1
(gilog (Af)i − (Af)i)
Den konstanten Term −log (gi!) konnen wir vernachlassigen, er spielt fur die Extrem-
stellen der Funktion keine Rolle. Die Anwendung des Logarithmus wird durch seine
Monotonie ermoglicht: Die Extremstellen werden auch dadurch nicht beeinflusst.
2.3.4 Maximierung der Log-Likelihood-Funktion
Um ein Maximum der Log-Likelihood-Funktion zu ermitteln, benotigen wir ein geeig-
netes Verfahren. Zunachst stellen wir fest, dass ` (f) konkav ist, jedes lokale Maximum
ist demnach auch ein globales Maximum. Betrachten wir die Funktion auf einem un-
beschrankten Gebiet, so erhalten wir als notwendige und hinreichende Bedingung fur
die Existenz eines Extremums, dass die partiellen Ableitungen nach allen Kompo-
nenten von f verschwinden. Da wir in der Praxis jedoch nur Werte von f innerhalb
eines beschrankten Gebietes beobachten, benotigen wir die sogenannte Kuhn-Tucker-
Bedingung [15]:
Satz: Kuhn-Tucker-Bedingung
Sei f ∈ Rn+, ` : Rn
+ → R eine konkave Funktion. An der Stelle f befindet sich ein
globales Maximum von ` genau dann, wenn gilt:
∂`
∂fk(f) = 0 fur fk > 0
∂`
∂fk(f) ≥ 0 fur fk = 0
Anschaulich bedeutet das: Wenn das lokale Maximum der Funktion innerhalb des be-
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 20
Abbildung 2.6: Maxima konkaver Funktionen(a) Maximum bei Betrachtung eines unbeschrankten Gebietes, (b) Maximum bei
Betrachtung nichtnegativer Werte
trachteten Bereichs, also der positiven f, liegt, so mussen dort alle Ableitungen gleich
Null sein. Ist jedoch das lokale Maximum außerhalb dieser Region, so ist die Funktion
im relevanten Bereich monoton fallend. Das Maximum befindet sich in diesem Fall ge-
nau auf dem Rand (s. Abb. 2.6).
Mit Hilfe der partiellen Ableitungen konnen wir nun den Gradienten von ` berechnen.
Nach kurzer Rechnung erhalten wir:
∇` (f) = AT(g
Af− 1
)Die Rechenoperationen geschehen hier komponentenweise, die 1 steht anstelle des m-
dimensionalen Einsvektors. Jedes globale Maximum erfullt nach der Kuhn-Tucker-
Bedingung
fAT(g
Af− 1
)= 0
Dies konnen wir, unter Verwendung der Normierungsbedingung AT1 = 1, umschreiben:
fAT(g
Af− 1
)= 0 ⇔ fAT
g
Af= fAT1 ⇔ fAT
g
Af= f
Wir haben somit das Problem in eine Fixpunktaufgabe transformiert. Um sie zu losen,
verwendet der EM-Algorithmus die Standard-Fixpunktiteration
2 Standard-Algorithmen bei der SPECT-Rekonstruktion 21
fk+1 = fkATg
Afk
Man kann zeigen, dass die Iteration fur einen beliebigen Startwert f 0 > 0 tatsachlich
gegen einen Fixpunkt und somit gegen das Maximum der Log-Likelihood-Funktion
` (f) konvergiert [15].
22
3 Dampfungskorrektur
3.1 Problemstellungen und Losungsansatze
Unser oberstes Ziel bei der Verarbeitung von SPECT-Daten besteht in einer moglichst
originalgetreuen Ermittlung der Verteilung der radioaktiven Substanz im Organismus.
An der Abschwachung des betrachteten Gewebes besteht kein direktes Interesse. Den-
noch ist ihre (zumindest naherungsweise) Bestimmung fur die Erzielung eines”schonen“
Ergebnisses unumganglich. In fruheren praktischen Anwendungen wurde die benotigte
Information uber die Dampfung haufig durch einen zusatzlichen Messvorgang wahrend
der eigentlichen SPECT-Messung gegeben. Dieses Vorgehen hat einige Nachteile: Es
entstehen sowohl hohere Kosten als auch eine langere Dauer der Behandlung, was wie-
derum mit mehr Stress sowie hoherer Strahlenbelastung fur den Patienten verbunden
ist. Aus diesen Grunden wird nach Moglichkeiten gesucht, die Notwendigkeit einer sol-
chen Extramessung zu umgehen [3].
In der Vergangenheit gab es verschiedene Ansatze zur Losung des Problems. Eine auf
Kosten der Genauigkeit relativ einfache Moglichkeit ist die Idee, die Abschwachung als
eine stuckweise konstante Funktion auf einer bekannten konvexen Menge zu betrachten
[7]. Bronnikov verwendete diskrete Konsistenzbedingungen, um die Dampfungsfunktion
ohne zusatzliche Messungen zu ermitteln [6]. Zudem existieren verschiedene Wege einer
simultanen Bestimmung von Verteilung und Abschwachung allein aus SPECT-Daten
[5], [8], [10].
Eine relativ neue Idee, die in der Praxis haufige Anwendung findet, besteht darin, eine
Naherung der Abschwachungsfunktion mit Hilfe einer Computertomographie (CT) vor
der eigentlichen SPECT-Messung zu bestimmen. Dadurch kann diese bei der Rekon-
struktion von f als bekannt angenommen werden. Die zusatzliche Messung wahrend der
SPECT-Aufnahme entfallt hierbei - wir sprechen daher von Dampfungskorrektur ohne
zusatzliche Messungen. Problematisch ist hierbei, dass die durch eine CT-Aufnahme
3 Dampfungskorrektur 23
gegebene Dampfung eventuell keine brauchbare Approximation des tatsachlichen µ
darstellt. Sie ist eine Momentaufnahme und berucksichtigt nicht die unvermeidbare
Bewegung des Patienten zwischen CT-Aufnahme und SPECT-Messung und wahrend
dieser, wie beispielsweise Atmung oder Herzschlag [3]. Eine Losung, wie der Herzschlag
mit berucksichtigt werden kann, liefert das sogenannte”Gated SPECT“. Die eigentli-
che SPECT-Messung wird hier von einem EKG (Elektrokardiogramm) begleitet, um
den Herzschlag des Patienten zu messen. Dabei werden die verschiedenen Phasen des
Herzschlags grob unterteilt und die Messdaten jeweils einer Phase zugeordnet. Dabei
misst das Gerat, in welcher Phase sich das Herz jeweils beim Auftreffen eines Photons
befindet und ordnet es der entsprechenden Phase zu. Jede dieser Phasen setzt eine
andere Abschwachungsfunktion voraus, die Rekonstruktion wird dann fur jede dieser
Herzphasen einzeln durchgefuhrt. Dieses Verfahren hat jedoch einen entscheidenden
Nachteil: Die ohnehin bereits geringe gegebene Datenmenge, die das SPECT-Gerat lie-
fert, wird noch weiter unterteilt, sodass die Rekonstruktion zwangslaufig noch um ein
Vielfaches ungenauer ausfallen muss. Die Aufteilung der Daten sollte also bei genaue-
ren Rekonstruktionsverfahren umgangen werden [3].
Um die Probleme entstehend durch Bewegung des Patienten dennoch miteinzubezie-
hen, prasentierte Natterer einen einfachen Ansatz [13]: Anstatt die zuvor ermittelte
Abschwachung beizubehalten, verwendet man eine affin-lineare Transformation dieser
und ermittelt die benotigten Zusatzinformationen mit Hilfe sogenannter Konsistenzbe-
dingungen fur die Radontransformation. Die Grundzuge dieser Idee werden in Kapitel
3.2 dargelegt.
Aufbauend auf der Idee Natterers entstand der folgende Ansatz [3], [4], [2]: Anstelle der
affin-linearen wird nun eine nichtlineare Transformation durchgefuhrt und die Konsis-
tenzbedingungen werden durch Regularisierung ersetzt. Verteilung und Abschwachung
werden dann in einem Algorithmus simultan ermittelt. In Kapitel 3.3 werden wir diese
Methode weiter veranschaulichen.
3 Dampfungskorrektur 24
3.2 Bestimmung mittels affin-linearer Transformation
eines vorhandenen Prototyps
3.2.1 Grundidee und Rekonstruktion
In diesem Kapitel wollen wir die Grundidee sowie den Algorithmus zur Bestimmung
der Abschwachungsfunktion µ darstellen, wie er in [13] zu finden ist. Grundlage ist die
Annahme, dass µ durch eine affin-lineare Transformation eines gegebenen Prototyps
µ0 approximiert werden kann. Das bedeutet, wir gehen von einer Dampfungsfunktion
der folgenden Gestalt aus:
µ (x) = µ0 (Ax+ a)
Hierbei ist A eine 2× 2-Matrix und a ein Vektor der Lange 2. Die Idee dabei ist, dass
mit einer solchen Transformation sowohl Verschiebungen als auch Drehungen darstell-
bar sind. A und a sind zunachst unbekannt und sollen nun mit Hilfe eines Algorithmus
ermittelt werden. Dabei ist im Allgemeinen nicht davon auszugehen, dass das Problem
der Bestimmung von f und µ aus den Messdaten g eindeutig losbar ist. Allerdings
konnen wir unter Verwendung zusatzlicher sogenannter”A-priori-Informationen“ eine
numerische Losung entwickeln. Bei diesen Extra-Informationen handelt es sich um so-
genannte Konsistenzbedingungen, die fur die gedampfte Radontransformation gelten
[12]. Basis fur den Algorithmus zur Bestimmung der Abschwachung ist die folgende
Bedingung:
Satz: Konsistenzbedingung fur die gedampfte Radontransformation
Seien f, µ ∈ S (R2). Dann gilt fur ganze Zahlen k,m mit 0 ≤ m < k:
∫ 2π
0
∫ ∞−∞
sme12(I+iH)Rµ(θ,s)+ikϕRµf (θ, s) dsdϕ = 0
Beweis: [12], siehe Kapitel 3.2.2.
Hierbei steht I fur die Identitat, ϕ ist der Grad der Drehung im Vektor θ =( cosϕsinϕ
)und
R im Exponentialterm steht fur die gewohnliche, ungedampfte Radontransformation
der Abschwachung µ. H bezeichnet die sogenannte Hilberttransformation:
3 Dampfungskorrektur 25
Definition: Hilberttransformation
Sei g : R→ R eine Funktion. Die Hilberttransformation von g ist definiert durch
(Hg) (s) :=1
π
∫ ∞−∞
g (t)
s− tdt
Die Konsistenzbedingung der gedampften Radontransformation ist unter der Annahme,
dass die Messdaten durch Rµf = g reprasentiert werden, nicht mehr von f abhangig.
Die Idee des Algorithmus ist es nun, die Gleichung unter der Voraussetzung µ (x) =
µ0 (Ax+ a) nach µ aufzulosen und somit A und a zu bestimmen. Bei der Berechung
der Gleichung ist folgende Aussage außerst hilfreich [13]:
Satz:
Sei µ (x) = µ0 (Ax+ a) nach µ mit einer nichtsingularen Matrix A. Setze ω =(Aθ⊥)
⊥
‖Aθ⊥‖ ,
wobei ‖·‖ der Euklidischen Norm entspricht. Dann gilt:
(1
2(I + iH)Rµ
)(θ, s) =
1
‖Aθ⊥‖
(1
2(I + iH)Rµ0
)(ω, (sAθ + a) · ω)
=: φ (A, a, θ, s)
Beweis:
Wir berechnen mit ω⊥ =(Aθ⊥
) ∥∥Aθ⊥∥∥:
3 Dampfungskorrektur 26
(Rµ) (θ, s) =
∫x·θ=s
µ (x) dx =
∫θ⊥µ (sθ + y) dy
=
∫µ(sθ + tθ⊥
)dt
=
∫µ0
(A(sθ + tθ⊥
)+ a)dt
=
∫µ0
(sAθ + tAθ⊥ + a
)dt
=
∫µ0
(sAθ + t
Aθ⊥
‖Aθ⊥‖·∥∥Aθ⊥∥∥+ a
)dt
=1
‖Aθ⊥‖
∫µ0
(sAθ + tω⊥ + a
)dt
=1
‖Aθ⊥‖
∫µ0
(((sAθ + a) · ω)ω +
((sAθ + a) · ω⊥
)ω⊥ + tω⊥
)dt
=1
‖Aθ⊥‖
∫µ0
(((sAθ + a) · ω)ω + tω⊥
)dt
=1
‖Aθ⊥‖Rµ0 (ω, (sAθ + a) · ω)
Fur Funktionen g, g0 mit g (s) = g0 (αs+ β) mit Konstanten α und β gilt fur die
Hilbert-Transformation:
(Hg) (s) =1
π
∫g (t)
s− tdt
=1
π
∫g0 (αt+ β)
s− tdt
αt+β→x=
1
π
∫g0 (x)
s− x−βα
· 1
αdx
=1
π
∫g0 (x)
αs+ β − xdx
= (Hg0) (αs+ β)
Damit folgt die Behauptung. 2
Desweiteren definieren wir:
3 Dampfungskorrektur 27
gm (A, a, ϕ) :=
∫ ∞−∞
smeφ(A,a,θ,s)g (θ, s) ds
Setzen wir dies in die Konsistenzbedingungsgleichung ein, erhalten wir die vereinfachte
Schreibweise
∫ 2π
0
eikϕgm (A, a, ϕ) dϕ = 0
fur ganze Zahlen 0 ≤ m < k. Aus dieser Gleichung wollen wir nun mit numerischen
Methoden die Unbekannten A und a bestimmen. Dazu mussen wir zunachst von der
bisher kontinuierlichen Formulierung des Problems auf eine diskrete Version ubergehen.
In der Praxis haben wir die Messdaten g nicht fur jede Richtung θ, also demnach nicht
fur jeden Blickwinkel ϕ, gegeben, sondern haben nur Zugriff auf eine endliche Anzahl
p an Winkeln ϕj = 2πjp
mit j = 0, . . . , p − 1. Die gegebenen Messrichtungen sind also
θj =( cosϕjsinϕj
). Wir haben also Auswertungen des Integrals in p Punkten. Damit konnen
wir das Integral mit Hilfe der Trapezregel diskretisieren und erhalten:
gm,k (A, a) =2π
p
p−1∑j=0
eikϕjgm (A, a, ϕj) = 0
fur 0 ≤ m < k ≤ P .
Wir haben nun ein nichtlineares Gleichungssystem mit P (P+1)2
komplexen Gleichungen
fur die sechs Unbekannten aus der 2× 2-Matrix A und dem zweidimensionalen Vektor
a. Daraus ergeben sich doppelt so viele, also P (P + 1) reelle Gleichungen. Das System
ist demnach eindeutig bestimmt fur P = 2 und uberbestimmt fur P > 2. Zur Losung
konnen verschiedene Minimierungsalgorithmen verwendet werden, fur die numerischen
Tests in [13] wurde der Levenberg-Marquardt-Algorithmus eingesetzt.
3.2.2 Beweis der Konsistenzbedingung
Wir erinnern an die noch zu beweisende Aussage der Konsistenzbedingung:
Satz: Konsistenzbedingung fur die gedampfte Radontransformation
Fur f, µ ∈ S (Rn) und ganze Zahlen 0 ≤ m < k gilt:
3 Dampfungskorrektur 28
∫R
∫ 2π
0
sme12(I±iH)Rµ(θ,s)±ikϕRµf (θ, s) dϕds = 0
Grundlage des folgenden Beweises ist [12].
Beweis:
Der Beweis lasst sich ubersichtlich in mehrere Schritte unterteilen. Die ersten vier
Schritte bilden die Vorarbeit zum eigentlichen Beweis und beschaftigen sich mit der
Berechnung einiger Fouriertransformierten. Wir beginnen mit
1. Schritt: Bestimmung der Fourierkoeffizienten von Dµ (x, θ)
Wir betrachten zunachst die sogenannte Fan-Beam-Transformation
Dµ (x, θ) :=
∫ ∞0
µ (x+ tθ) dt
und bilden ihre Fourierreihe in ϕ bei festem x:
Dµ (x, θ) =∑l
pl (x) e−ilϕ
mit θ =( cosϕsinϕ
). Die Koeffizienten sind gegeben durch
pl (x) =1
2π
∫ 2π
0
eilϕDµ (x, θ) dϕ
=1
2π
∫ 2π
0
eilϕ∫ ∞0
µ (x+ tθ) dtdϕ
Mit Hilfe einer Polarkoordinatentransformation erhalten wir
3 Dampfungskorrektur 29
pl (x)ϕ→ϕ+π
=1
2π(−1)l
∫ 2π
0
∫ ∞0
µ (x− tθ) eilϕdtdϕ
=1
2π(−1)l
∫ 2π
0
∫ ∞0
µ (x− tθ) 1
−teilϕ (−t) dtdϕ
=1
2π(−1)l
∫ 2π
0
∫ ∞0
µ (x− tθ) vl (−tθ) (−t) dtdϕ
y=−tθ=
1
2π(−1)l
∫R2
µ (x− y) vl (y) dy
=1
2π(−1)l (µ ∗ vl) (x)
mit vl (rθ) := 1reilϕ fur r > 0.
2. Schritt: Berechnung der Fouriertransformierten von pl
Wir berechnen die Fouriertransformierte von vl an der Stelle ξ = ρ(cosψsinψ
), wobei wir
verwenden, dass
θ · ξ = ρ (cosϕcosψ + sinϕsinψ) = ρcos (ϕ− ψ)
vl (ξ) =1
2π
∫R2
e−ix·ξvl (x) dx
=1
2π
∫ ∞0
∫ 2π
0
e−irθ·ξvl (rθ) rdϕdr
=1
2π
∫ ∞0
∫ 2π
0
e−irρcos(ϕ−ψ)+ilϕdϕdr
ϕ→ϕ+ψ=
1
2π
∫ ∞0
∫ 2π
0
e−irρcosϕ+il(ϕ+ψ)dϕdr
=1
2π
∫ ∞0
eilψ∫ 2π
0
e−irρcosϕ+ilϕdϕdr
Die Besselfunktion erster Gattung der Ordnung l ist gegeben durch [1]
Jk (x) =i−k
2π
∫ 2π
0
eixcosϕ−ikϕdϕ
3 Dampfungskorrektur 30
Indem wir diese Funktion in die letzte Gleichung einsetzen erhalten wir:
vl (ξ) =
∫ ∞0
eilψi−l1
ρJ−l (−rρ) dr
=1
ρi−l∫ ∞0
J−l (−rρ) dr eilψ
=1
ρi−l∫ ∞0
Jl (rρ) dr eilψ
wobei wir die Eigenschaft der Besselfunktion
J−l (x) = Jl (−x) , l ∈ Z
verwenden.
Eine weitere Eigenschaft der Besselfunktion ist die folgende [1]:
∫ ∞0
Jl (t) dt = εl =
1, falls l ≥ 0
(−1)l , falls l < 0
Damit folgt:
vl (ξ) = εli−lρ−1eilψ = εli
−lvl (ξ)
Die Fouriertransformierte einer Faltung zweier Funktionen f und g berechnet sich mit-
tels [12]
(f ∗ g) = (2π)n2 f g
Fur die Fouriertransformation von pl gilt dann:
pl (ξ) =1
2π(−1)l (µ ∗ vl) (ξ) = (−1)l µ (ξ) vl (ξ) = (−1)l µ (ξ) εli
−lvl (ξ) = εlilµ (ξ) vl (ξ)
3 Dampfungskorrektur 31
3. Schritt: Bestimmung der Fourierkoeffizienten von h (θ, x · θ)
Sei nun h eine Funktion auf dem Einheitszylinder Z = S1 × R und ihre Fourierreihe
gegeben durch
h (θ, x · θ) =∑l
ql (x) e−ilϕ
ql (x) =1
2π
∫ 2π
0
eilϕh (θ, x · θ) dϕ =1
2π
∫S1
eilϕh (θ, x · θ) dx =1
2πR∗(eilϕh
)(x)
mit der adjungierten Radontransformation R∗.
4.Schritt: Berechnung der Fouriertransformierten von ql
Fur die Fouriertransformierten von ql erhalten wir unter Zuhilfenahme von [12], Theo-
rem 1.4:
ql (ξ) =1
2π
(R∗(eilϕh
))(ξ) (3.1)
= (2π)−12 |ξ|−1
[eilψh
(ξ
|ξ|, |ξ|)
+ eil(ψ+π)h
(− ξ
|ξ|,− |ξ|
)](3.2)
= (2π)−12 vl (ξ)
[h
(ξ
|ξ|, |ξ|)
+ (−1)l h
(− ξ
|ξ|,− |ξ|
)](3.3)
5. Schritt: Kern des Beweises
Wir definieren nun Funktionen h+, h− auf dem Einheitszylinder Z durch
h± (θ, σ) := (2π)12
1
2(1± sgn (σ)) µ (σθ)
Durch Einsetzen von h+ in ql nach obiger Berechnung erhalten wir:
3 Dampfungskorrektur 32
ql (ξ) = (2π)−12 vl (ξ)
[(2π)
12
1
2(1 + sgn (|ξ|)) µ
(ξ
|ξ||ξ|)
+ (−1)l (2π)12
1
2(1 + sgn (− |ξ|)) µ
(− ξ
|ξ|(− |ξ|)
)]= (2π)−
12 vl (ξ)
[(2π)
12
1
22µ
(ξ
|ξ||ξ|)]
= vl (ξ) µ (ξ)
Nun definieren wir eine Funktion u+ auf dem Einheitszylinder durch
u+ (x, θ) := h+ (θ, x · θ)−Dµ(x, θ⊥
)Wir schreiben u+ als Fourierreihe in ϕ mit Koeffizienten ul. Dazu benotigen wir
zunachst die Koeffizienten der Fourierreihe von Dµ(x, θ⊥
). Es gilt:
Dµ (x, θ) =∑l
pl (x) e−ilϕ
Mit θ =( cosϕsinϕ
)ist θ⊥ =
( cos(ϕ+π2 )
sin(ϕ+π2 )
)und damit:
Dµ(x, θ⊥
)=∑l
pl (x) e−ilϕe−ilπ2
=∑l
i−lpl (x) e−ilϕ
Somit hat die Fourierreihe von Dµ(x, θ⊥
)die Koeffizienten i−lpl. Fur die Fouriertrans-
formierten der Koeffizienten ul der Fourierreihe von u+ erhalten wir also:
ul (ξ) = ql (ξ)− i−lpl (ξ) = vl (ξ) µ (ξ)− εlµ (ξ) vl (ξ) = (1− εl) µ (ξ) vl (ξ)
Betrachtet man die Definition von εl, so sieht man, dass ul = 0 ist fur l ≥ 0. Demnach
enthalt die Fourierreihe von u+ nur Terme der Ordnung l ≤ 0. Das heißt:
3 Dampfungskorrektur 33
u+ (x, θ) =∑l<0
ul (x) e−ilϕ =∑l≥0
u−l (x) eilϕ
Die Fourierreihe enthalt also nur Exponentialterme mit Exponenten ≥ 0. Fur die Fou-
rierreihe von eu+ gilt dann:
eu+(x,θ) =∞∑j=0
(u+ (x, θ))j
j!=∞∑j=0
(∑∞l=0 u−l (x) eilϕ
)jj!
Also enthalt die Fourierreihe von eu+ nur Terme mit Exponenten ≥ 0, also Terme der
Ordnung l ≤ 0. Weiter gilt: (x · θ)m ist ein trigonometrisches Polynom vom Grad ≤ m.
Es lasst sich also darstellen als
(x · θ)m =m∑
n=−m
cneinϕ
mit Koeffizienten cn. Fur k > m ≥ 0 gilt dann:
(x · θ)m eikϕ = eikϕm∑
n=−m
cneinϕ = c−me
i(k−m)ϕ + . . .+ cmei(k+m)ϕ
Wegen k > m ist der kleinste Exponent k−m > 0, daher enthalt auch die Fourierreihe
von (x · θ)m eikϕ nur Terme mit echt positivem Exponenten. Damit enthalt die Fourier-
reihe von (x · θ)m eikϕeu+(x,θ) nur Terme echt positiver Ordnung.
Durch Integration uber 0 bis 2π erhalten wir damit:
∫ 2π
0
(x · θ)m eikϕeu+(x,θ)dϕ =
∫ 2π
0
(x · θ)m eikϕ+h+(θ,x·θ)e−Dµ(x,θ⊥)dϕ = 0
Definieren wir nun Funktionen
w± (θ, s) := sme±ikϕ+h±(θ,s)
so sehen wir, dass obiges Integral genau der adjungierten gedampften Radontransfor-
mation der Funktion w+ entspricht, d.h.
3 Dampfungskorrektur 34
R∗µw+ = 0
Auf die gleiche Weise konnen wir zeigen, dass auch
R∗µw− = 0
fur w− analog fur k > m. Damit folgt auch, da R∗µ die zu Rµ adjungierte Radontrans-
formation bezeichnet, aus der Definition zweier adjungierter Operatoren, dass
(w±, Rµf)L2(Z)=(R∗µw±, f
)L2(R2)
= 0
bezuglich des L2-Skalarproduktes. Formulieren wir die Gleichung aus, so sehen wir:
∫R
∫ 2π
0
sme±ikϕ+h±(θ,s)Rµf (θ, s) dϕds = 0
Es bleibt also noch zu zeigen, dass h± = 12
(I ± iH)Rµ ist. Dazu benotigen wir die
folgenden beiden Eigenschaften der Fouriertransformierten der Radon- und der Hil-
berttransformation:
(1) (Rµ) (θ, σ) = (2π)12 µ (σθ) (Fourier-Slice-Theorem)
(2) (Hh) (σ) =1
isgn (σ) h (σ) ([12] VII.1.11)
Dann folgt mit (1)
h± (θ, σ) = (2π)12
1
2(1± sgn (σ)) µ (σθ)
=1
2(1± sgn (σ)) (Rµ) (θ, σ)
und mit (2)
3 Dampfungskorrektur 35
Abbildung 3.1: Affin-lineare Transformation bei kunstlichen Daten, aus [13]Oben links/Mitte: Vom Computer generierte Verteilung/Abschwachung.
Unten links: Abschwachungsprototyp.Unten Mitte: Mit Hilfe des Algorithmus berechnete Abschwachung.
Rechts: Sinogramm der verrauschten Daten
h± (θ, σ) =1
2(I ± iH)Rµ (θ, σ)
und wir haben die Behauptung gezeigt. 2
3.2.3 Numerische Ergebnisse
Mit dem zuvor erlauterten Algorithmus zur Bestimmung der Abschwachungsfunktion
wurden verschiedene numerische Experimente durchgefuhrt [13]. In einem ersten Test
wurden Daten vom Computer generiert und eine Verteilungs- und Dampfungsfunktion
kunstlich erstellt. Diese Daten wurden zusatzlich mit einem Rauschen belastet, das
einer relativen mittleren quadratischen Abweichung von 12% entspricht. Fur die Be-
rechnung der gedampften Radontransformation wurden 64 Messrichtungen sowie 65
Messungen pro Richtung verwendet. Als Prototyp der Dampfungsfunktion diente ei-
ne vereinfachte Version der zuvor festgelegten tatsachlichen Abschwachung, der mit
wenigen Kenntnissen uber die zu messende Region modelliert werden kann. Das Er-
gebnis des Verfahrens ist in Bild 3.1 zu sehen. Nach nur wenigen Iterationsschritten
ergab sich selbst bei verrauschten Daten eine Abschwachung, die optisch kaum von der
tatsachlichen zu unterscheiden ist.
Fur einen weiteren Test, der einen realistischeren Fall untersucht, wurden eine Verteilungs-
3 Dampfungskorrektur 36
Abbildung 3.2: Affin-lineare Transformation bei echten Daten, aus [13]Oben links/Mitte: Verteilung/Abschwachung eines menschlichen Gehirns.
Unten links: Abschwachungsprototyp.Unten Mitte: Mit Hilfe des Algorithmus berechnete Abschwachung.
Rechts: Sinogramm der verrauschten Daten
und Abschwachungsfunktion verwendet, wie sie typischerweise bei einer Messung des
menschlichen Gehirns auftreten. Hier wurden die erhaltenen Daten zusatzlich mit einem
Rauschen von 17% Abweichung erganzt. Als Prototyp diente eine ahnliche Abbildung
wie zuvor. Auch hier ergaben sich trotz Rauschen und Vereinfachungen bei der Model-
lierung brauchbare Ergebnisse, die nahe an der tatsachlichen Abschwachung lagen. Die
Ergebnisse des realistischeren Modells werden in Bild 3.2 dargestellt.
3.3 Bestimmung mittels nichtlinearer Transformation
3.3.1 Grundidee und Rekonstruktion
In diesem Abschnitt wollen wir auf die Idee, die Dampfungsfunktion mit Hilfe ei-
ner nichtlinearen Transformation eines Prototyps zu bestimmen, naher eingehen. Die
Grundlage bilden [4], [2] und [3]. Wir gehen hier davon aus, dass f, µ ∈ S (R2,R),
mit kompaktem Trager auf Ω ⊂ R2. Uberdies sei f ≥ 0. Fur einen gegebenen Ab-
schwachungsprototyp µ0 ist dann µ gegeben durch
µ (x) = µ0 (φ (c;x))
3 Dampfungskorrektur 37
mit einer nichtlinearen Transformation φ (c;x) mit Transformationskonstanten c. Die
Idee der nichtlinearen Modellierung ruhrt daher, dass sich bestimmte typische Bewe-
gungen nicht mit linearen Verschiebungen darstellen lassen. Halt der Patient beispiels-
weise ein Bein in gleicher Position, wahrend er das andere zur Seite bewegt, entsteht
ein nichtlinearer Zusammenhang zur Ausgangsposition. Wir modellieren wie zuvor den
Zusammenhang zwischen Daten, Verteilungsfunktion und Abschwachung mittels der
gedampften Radontransformation, d.h. Rµf = g. Dieses Problem ist im Allgemeinen
nicht analytisch losbar. Um dennoch eine naherungsweise Losung zu ermitteln, macht
man einen Variationsansatz: Man sucht zunachst einen Minimierer des Funktionals
J (f, µ) :=1
2‖Rµf − g‖2L2
fur Funktionen f, µ ∈ S (R2). Fur die Wahl von φ sind verschiedene Ansatze denkbar.
Um uns auf eine engere Auswahl zu beschranken, wahlen wir hier eine Transformation
mittels B-Spline-Funktionen der Form
φ (c;x) = x+∑i
ciBi (x)
Je nach dem wie viele gewichtete Splines wir verwenden, haben wir es mit einer großeren
Anzahl an unbekannten Konstanten ci zu tun. Bei der Minimierung des obigen Funk-
tionals stoßen wir demnach im Allgemeinen auf eine hohe Zahl an moglichen Losungen.
Um dies zu verhindern, benotigen wir, wie auch schon im Falle der affin-linearen Trans-
formation, eine zusatzliche Einschrankung. Anstelle der zuvor verwendeten Konsistenz-
bedingungen fugen wir hier an das zu minimierende Funktional einen Regularisierungs-
term an. Das fuhrt zu
Jspline (f, c) := J (f, µ0 (φ (c;x))) + α ‖c‖22
mit einem positiven Regularisierungsparameter α. Der Regularisierungsterm α ‖c‖22sorgt dafur, dass Werte mit großer L2-Norm zusatzlich bestraft werden. Dieses Pro-
blem ist linear in f und nichtlinear in c. Wir versuchen also, mit Hilfe eines Algorithmus
die Verteilung f und die Parameter c, welche die Dampfungsfunktion eindeutig festle-
gen, simultan zu bestimmen.
Um mit numerischen Methoden arbeiten zu konnen, benotigen wir zunachst wieder
3 Dampfungskorrektur 38
einen Ubergang vom kontinuierlichen zum diskreten Modell. Im kontinuierlichen ist
die gedampfte Radontransformation gegeben durch
(Rµf) (θ, s) =
∫x·θ=s
f (x) e−∫L(x) µ(t)dtdx
Wir haben nun zur Auswertung p Blickwinkel θj fur j = 1, . . . , p zur Verfugung. Zu
jedem Blickwinkel haben wir q Kollimatorkanale x1, . . . , xq und damit Linien L1, . . . , Lq
zu jedem Winkel ϕj. Jede Linie wird nun noch einmal unterteilt in R Linienabschnitte
der Lange h. Damit erhalten wir: xik ist der Mittelpunkt des k-ten Segments der Linie
Li. Unser zweidimensionales Objekt ist demnach fur jeden Winkel j aufgeteilt in Punkte
xik. Als Approximation der gedampften Radontransformation erhalten wir damit:
Pµf (xi, θj) := hR∑k=1
f (xik) e−h
∑Rn=k µ(xin)
Nun sollen unsere Funktionen f, µ nicht mehr kontinuierlich, sondern ebenfalls diskret
sein. Wir schreiben f und µ also als Vektoren f ∈ Rl und µ ∈ Rm und fassen f und µ
in obiger Gleichung als Interpolation der Vektoren auf. Mit Interpolationsoperatoren
Im fur µ und I l fur f setzen wir dann
(P (µ) f)ij := P(Imµ)
(I lf)
(xi, θj)
In Bezug auf die Messdaten g soll also gelten:
P (µ) f = g
Damit erhalten wir, dass g ein Vektor der Lange qp und P (µ) eine Matrix der Dimen-
sion qp× l ist.
Die folgenden Berechnungen werden mit den diskreten Versionen der jeweiligen Funk-
tionen durchgefuhrt. Der Einfachheit und Ubersichtlichkeit halber schreiben wir aber
weiterhin Rµf statt P (µ) f .
3 Dampfungskorrektur 39
3.3.2 Der Algorithmus
Der Algorithmus zur Minimierung von Jspline verwendet die Coordinate-Search-Methode.
Grundlage ist dabei die folgende:
Setze Startwerte c0 = 0, f 0 = konstant > 0
1. Berechne Minimierer fk+1 von J(f, ck
), ckfest
2. Berechne Minimierer ck+1 von J(fk+1, c
), fk+1fest
Fur das erste, lineare Unterproblem der Minimierung fur festes ck verwenden wir
die Modified-Residual-Norm-Steepest-Descent-Methode [11], kurz MRNSD, das zwei-
te, nichtlineare Problem wird mit Hilfe von Gauss-Newton gelost [4].
Die MRNSD-Methode ist ein Verfahren zur Minimierung des Funktionals
J (f) = J(f, ck
)=
1
2‖Rµf − g‖2L2
fur festes ck. Der Regularisierungsterm ist hierbei vernachlassigbar, da er nicht von f
abhangt. Zu beachten ist hierbei, dass f nur Werte großer gleich Null annehmen soll.
Um die Nichtnegativitat zu gewahrleisten, parametrisieren wir f durch f = ez. Dadurch
erhalten wir das”modifizierte“ Problem
J (z) =1
2‖Rµ (ez)− g‖2
Hier sollte erwahnt werden, dass f nun nicht mehr den Wert 0 annehmen kann, also nur
noch echt positive Werte berucksichtigt werden. Das Modell wird dadurch allerdings
nicht wesentlich eingeschrankt und die Parametrisierung vereinfacht das Problem deut-
lich: Ohne sie ware die Beibehaltung der Bedingung f ≥ 0 im Algorithmus um einiges
komplizierter. Zur Minimierung verwendet die MRNSD-Methode ein Iterationsverfah-
ren der folgenden Gestalt:
fk+1 = fk + akpk
Hierbei ist pk die Richtung des lokalen Minimums von J, die sogenannte”Descent
Direction“, und ak die Schrittweite. Wir gehen in einem Iterationsschritt also vom alten
3 Dampfungskorrektur 40
Wert aus in Richtung des lokalen Minimums mit Schrittweite ak. Zur Bestimmung von
pk berechnen wir den Gradienten von J mit Hilfe der Kettenregel:
∇zJ (f) = diag (f)∇fJ (f) = diag (f)RTµ (Rµf − g)
Wir setzen nun
pk = diag (f)k RTµ
(Rµf
k − g)
Um nun die optimale Schrittweite zu ermitteln, eignen sich die Armijo-Goldstein-
Bedingungen [9]. Die erste Bedingung stellt sicher, dass wir uns dem Minimum des
Funktionals J mit dem gewahlten ak hinreichend annahern. Dazu sei c1 ∈ (0, 1) eine
Konstante, sodass fur ak gilt:
J(fk + akp
k)≤ J
(fk)
+ c1ak(pk)T ∇J (fk) (Armijo−Bedingung)
Die Goldstein-Bedingung garantiert im Gegenzug, dass die Schrittweite nicht zu klein
gewahlt wird, d.h. fur c2 ∈ (0, 1):
J(fk + akp
k)≥ J
(fk)
+ c2ak(pk)T ∇J (fk) (Goldstein−Bedingung)
Der zweite Teil des Minimierungsproblems wird mit Hilfe des Verfahrens von Gauss-
Newton gelost. Fur ein festes f = fk+1 minimieren wir das Funktional
J (c) = J (f, c) =1
2‖Rµf − g‖2L2
+ α ‖c‖22
nach c. Die Idee des Gauss-Newton-Verfahrens ist es nun, das Funktional mittels Tay-
lorentwicklung zu approximieren und die Naherung zu minimieren. Dazu bestimmen
wir die Taylorentwicklung von J um c:
J (c+ s) ≈ J (c) +∇J (c) s+1
2sT∇2J (c) s := J (s)
Ausgehend von einem fest gewahlten c bestimmen wir nun das Minimum von J, indem
wir J nach s minimieren. Dazu leiten wir J nach s ab und erhalten
3 Dampfungskorrektur 41
Abbildung 3.3: Armijo-Goldstein-Bedingung zur Ermittlung der optimalen Schrittwei-te, aus [9]
∇J (s) = ∇J (c) +∇2J (c) s
Aus der notwendigen Bedingung ∇J = 0 folgt dann das Gleichungssystem
−∇J (c) = ∇2J (c) s
Um dies noch expliziter aufzuschreiben, berechnen wir die erste und zweite Ableitung
von J nach c mit Hilfe der Kettenregel:
∇J (c) =
(dRµf
dc
)T(Rµf − g) + 2αc
∇2J (c) ≈(dRµf
dc
)TdRµf
dc+ Terme zweiter Ordnung + αI
Die Terme zweiter Ordnung sind hierbei vernachlassigbar, sodass sich eine vereinfachte
Darstellung der zweiten Ableitung ergibt.
Bei der numerischen Berechnung von Teilen sowohl der MRNSD- als auch der Gauss-
3 Dampfungskorrektur 42
Newton-Methode ergeben sich einige Schwierigkeiten. Der numerisch aufwandige Teil
bei MRNSD ist die Implementation von Rµ sowie R⊥µ . Im diskreten Fall ist Rµ ≈ P (µ)
eine qp×l-Matrix. Bei der Ableitung nach c entstehen zu jedem Eintrag c neue Eintrage.
Dadurch kann die Ableitung von Rµ nach c in der Gauss-Newton-Methode im Diskre-
ten aufgefasst werden als eine dreidimensionale Matrix der Dimension qp × l × c. In
der Praxis ergeben sich durch die Transformation mittels B-Spline-Funktionen 128 Frei-
heitsgrade c, wodurch typischerweise eine Dimension der Ableitung von 2562×1282×128
entsteht [4], [2].
3.3.3 Numerische Ergebnisse
Um die neue Verfahrensweise zu uberprufen, wurden Tests mit Hilfe des XCAT-Phantoms
durchgefuhrt [3]. Bei XCAT (4D Extended Cardiac Torso) handelt es sich um ein rea-
listisches und detailliertes Computermodell eines menschlichen Torso, fur das anato-
mische Strukturen und physiologische Funktionen anhand von CT-Bildern modelliert
wurden [16]. Unter Verwendung dieses Modells ist es moglich, die Ergebnisse des hier
vorgestellten Rekonstruktionsverfahrens mit den tatsachlichen Werten fur Verteilung
und Abschwachung zu vergleichen. Fur die Tests wurde die normale Atmungsaktivitat
wahrend der SPECT-Messungen berucksichtigt, wobei man von den folgenden Annah-
men ausging: Der durchschnittliche Zeitraum einer Atmung, d.h. beginnend bei 0%
eingeatmet bis dorthin zuruck, betragt 5 Sekunden. Nach 2,5 Sekunden ist mit 96%
Einatmung in etwa das Maximum erreicht. Die Bewegung wahrend einer solchen Atem-
periode liegt bei 1,2 cm.
Mit Hilfe des XCAT-Phantoms wurden die durchschnittliche Verteilung f0 innerhalb
eines Zeitraums von 5 Sekunden sowie zwei Abschwachungsfunktionen µ0 und υ be-
stimmt, sodass sich diese analog zur Atmungsverschiebung unterscheiden. Mittels des
Rekonstruktionsverfahrens durch nichtlineare Transformation wurden nun f0 und µ0
bestimmt, gegeben Messdaten g und einen Abschwachungsprototyp υ. In Bild 3.3 wer-
den diese Daten anhand des XCAT-Phantoms dargestellt.
3.4 Vergleich und Auswertung
Ein wesentlicher Unterschied der beiden Methoden besteht darin, dass im Fall der affin
linearen Transformation die Abschwachung und die Verteilung mit Hilfe der Konsis-
3 Dampfungskorrektur 43
Abbildung 3.4: Von links nach rechts: Durchschnittliche (tatsachliche) Verteilungsfunk-tion, Abschwachungsprototyp, tatsachliche Abschwachung, Sinogrammder Messdaten, aus [3]
tenzbedingung separat ermittelt werden. Wahrenddessen geschieht die Bestimmung im
nichtlinearen Fall simultan fur beide Funktionen.
In [3] wurden mit Hilfe des XCAT-Phantoms verschiedene Methoden der Bewegungs-
korrektur mit den tatsachlichen Werten fur Verteilung und Abschwachung verglichen.
Hierbei wurde der relative Fehler jeweils bei einer Rekonstruktion ohne Bewegungs-
korrektur, mittels affin-linearer sowie mittels nichtlinearer Transformation berechnet.
Die Bewegung des Patienten wurde dabei einmal als normale Atmungsaktivitat, mit 1,2
cm, sowie zum Vergleich mit 2,4 cm simuliert. Das Ergebnis der relativen Fehlerbestim-
mung ist in Tabelle 3.1 zusammengefasst. In Bild 3.4 ist der absolute Fehler dargestellt.
Bisher wurde die relativ neue Methode der nichtlinearen Transformation ausschließlich
mit Hilfe des XCAT-Phantoms getestet. Anhand der so gelieferten Sinogramme wurden
die Algorithmen unter moglichst realitatsnahen Bedingungen gepruft. Die Ergebnisse
lassen auf eine deutliche Verbesserung des Fehlers im Vergleich zum Modell der linea-
ren Transformation schließen. Eine mogliche Erklarung ware, wie bereits genannt, die
Tatsache, dass sich bestimmte typische Bewegungen des Patienten nicht mittels linea-
rer Transformation darstellen lassen. Das nichtlineare Modell ermoglicht eine deutlich
breitere Variation an Transformationen durch die hohere Zahl an Freiheitsgraden im
Algorithmus - statt sechs im affin-linearen Modell treten hier typischerweise 128 zu be-
stimmende Werte auf. Fur die Zukunft stehen Tests mit”echten“ Daten unter wirklich
realistischen Bedingungen noch aus.
Mogliche Ansatzpunkte fur Verbesserungen waren im nichtlinearen Fall unter anderem
die Art der Transformation sowie der Regularisierung. Um die Auswahl der moglichen
Transformationen mehr einzuschranken, wurden statt einer”allgemeinen“ (in der Pra-
xis nicht umsetzbaren) nichtlinearen Transformation B-Spline-Funktionen verwendet.
Aufgrund der scheinbaren Verbesserung zum affin-linearen Modell stellt sich die Fra-
3 Dampfungskorrektur 44
err (f) =‖f−f0‖2‖f0‖2
err (µ) =‖µ−µ0‖2‖µ0‖2
Bewegung keine Korrektur linear nichtlinear keine Korrektur linear nichtlinear
1,2cm 48,43% 31,16% 18,84% 29,13% 32,36% 17,12%
2,4cm 81,48% 36,50% 18,19% 39,65% 47,06% 23,31%
Tabelle 3.1: Tabelle der relativen Fehler von Verteilung und Abschwachung, aus [3]
Abbildung 3.5: Graphische Darstellung des absoluten Fehlers bei der Bewegungskor-rektur, aus [3]
Die erste Zeile bezieht sich auf die Verteilung der radioaktiven Substanz, die zweiteauf die Dampfungsfunktion.
ge, ob eine andere Art der Transformation weitere Vorteile mit sich brachte. Ebenso
konnte die Art der Regularisierung moglicherweise noch besser angepasst werden.
45
Abbildungsverzeichnis
1.1 Photonenemission aus einem dreidimensionalen Objekt mit umgebenden
Kollimatoren, aus [18]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Zweidimensionales Objekt im Koordinatensystem, aus [18] . . . . . . . 3
2.1 Parallel- und Ruckprojektion, aus [18] . . . . . . . . . . . . . . . . . . . 9
2.2 (a) Ursprungliches zu rekonstruierendes Bild, (b) Aus Projektionen ent-
standenes Sinogramm, aus [18] . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Rekonstruktion eines Punktes mit tatsachlicher Verteilung 0, aus [18] . 11
2.4 Ergebnis der gefilterten Ruckprojektion mit verschiedenen Anzahlen von
Blickwinkeln, aus [18] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.5 Diskretisiertes Modell der Objektregion aus [18] . . . . . . . . . . . . . 15
2.6 Maxima konkaver Funktionen . . . . . . . . . . . . . . . . . . . . . . . 20
3.1 Affin-lineare Transformation bei kunstlichen Daten, aus [13] . . . . . . 35
3.2 Affin-lineare Transformation bei echten Daten, aus [13] . . . . . . . . . 36
3.3 Armijo-Goldstein-Bedingung zur Ermittlung der optimalen Schrittweite,
aus [9] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.4 Von links nach rechts: Durchschnittliche (tatsachliche) Verteilungsfunk-
tion, Abschwachungsprototyp, tatsachliche Abschwachung, Sinogramm
der Messdaten, aus [3] . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.5 Graphische Darstellung des absoluten Fehlers bei der Bewegungskorrek-
tur, aus [3] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
46
Literaturverzeichnis
[1] M. Abramowitz and I. A. Stegun. Handbook of Mathematical Functions With
Formulas, Graphs, and Mathematical Tables. Dover Publications, 1972. 29, 30
[2] S. Barendt. Nonlinearly transformed attenuation prototype for spect reconstruc-
tion. Presentation, 2011. Munster. 23, 36, 42
[3] S. Barendt and J. Modersitzki. Spect reconstruction with a nonlinear transformed
attenuation prototype. Institut of Mathematics and Image Computing, University
of Lubeck. 22, 23, 36, 42, 43, 44, 45
[4] S. Barendt and J. Modersitzki. Nonlinearly transformed attenuation prototype
for spect reconstruction. Presentation, 2011. Edinburgh. 23, 36, 39, 42
[5] A. V. Bronnikov. Approximate reconstruction of attenuation map in spect ima-
ging. IEEE Transactions on Nuclear Science, 42:1483–1488, 1995. 22
[6] A. V. Bronnikov. Reconstruction of attenuation map using discrete consistency
conditions. IEEE Transactions on Medical Imaging, 19:451–462, 2000. 22
[7] L. Chang. A method for attenuation correction in radionuclide computed tomo-
graphy. IEEE Transactions on Nuclear Science, NS-25:638–643, 1978. 22
[8] V. Dicken. A new approach towards simultaneous activity and attenuation recon-
struction in emission tomography. Inverse Problems, 15:931960, 1999. 22
[9] H. Dirks. Medical Image Segmentation Including Particle Methods. PhD thesis,
Westfalische Wilhelms-Universitat Munster, 2011. 7, 40, 41, 45
[10] C. Mennessier, F. Noo, R. Clackdoyle, and L. Desbat G. Bal. Attenuation correc-
tion in spect using consistency conditions for the exponential ray transform. Phys.
Med. Biol., 44:24832510, 1999. 22
[11] J. Nagy and Z. Strakos. Enforcing nonnegativity in image reconstruction algo-
rithms. Department of Mathematics and Computer Science, Emory University,
Literaturverzeichnis 47
Atlanta; Institute of Computer Science, Academy of Sciences of the Czech Repu-
blic, Prague. 39
[12] F. Natterer. The Mathematics of Computerized Tomography. B. G. Teubner, John
Wiley and Sons, 1986. 2, 6, 9, 12, 24, 28, 30, 31, 34
[13] F. Natterer. Determination of tissue attenuation in emission tomography of op-
tically dense media. Inverse Problems, 9:731–736, 1993. 23, 24, 25, 27, 35, 36,
45
[14] F. Natterer. Inversion of the attenuated radon transform. Inverse Problems,
17:113–119, 2000. 9
[15] F. Natterer and F. Wubbeling. Mathematical Methods in Image Reconstruction.
SIAM, 2001. 2, 4, 5, 9, 13, 19, 21
[16] W. P. Segars, M. Mahesh, T. J. Beck, E. C. Frey, and B. M. W. Tsui. Realistic ct
simulation using the 4d xcat phantom, 2008. Department of Radiology, The Duke
University Medical Center, Durham; The Russell H. Morgan Department of Ra-
diology and Radiological Science, Johns Hopkins Medical Institutions, Baltimore.
42
[17] B. Setzepfandt. ESNM: Ein rauschunterdruckendes EM-Verfahren fur die Emis-
sionstomographie. PhD thesis, Westfalische Wilhelms-Universitat Munster, 1992.
13
[18] M. N. Wernick and J. N. Aarsvold. Emission Tomography - The Fundamentals of
PET and SPECT. Elsevier, 2004. 2, 3, 4, 6, 8, 9, 10, 11, 12, 13, 15, 17, 45