Physikalische Grundlagen der Nuklearmedizin/ Fourier-Methoden

Einleitung

Bearbeiten

Fourier-Methoden werden unten in recht allgemeiner Form beschrieben. Wir werden sehen, dass jede beliebige Wellenform mathematisch in mehrere Sinuswellen unterschiedlicher Frequenz und Amplitude zerlegt werden kann. Hierzu führt man häufig eine sogenannte Fouriertransformation durch. Wir werden sehen, dass das resultierende Fourier-Spektrum modifiziert werden kann um bestimmte Frequenzen zu verstärken oder zu unterdrücken. Der Zweck der Übungsaufgaben ist es, die Modulationsübertragungsfunktion (MTF) und die Filterstufe eines tomographischen Rekonstruktionsprozesses im Detail zu behandeln.

Man beachte dass hier eine Tabelle heruntergeladen werden kann, mit der man die verschiedenen unten angegeben Prozesse simultan untersuchen kann.

Periodische Funktionen

Bearbeiten

Eine periodische Funktion ist eine Funktion, die sich (nach einem bestimmten Raum- oder Zeitintervall regelmäßig) wiederholt. Einfache Beispiele sind die Sinusfunktion oder die Rechteckfunktion. Die Geschwindigkeit der zeitlichen Wiederholung einer periodischen Funktion wird durch ihre Frequenz beschrieben. Die Frequenz ist die Anzahl der Perioden, die einen festen Punkt pro Zeiteinheit passieren. Die Frequenz wird üblicherweise in Perioden pro Sekunde angegeben. Ihre Si-Einheit ist Hertz mit dem Formelzeichen Hz. Es gilt:

 

Ein medizinisches Beispiel einer periodischen Funktion ist der Herzschlag, wie er im Elektrokardiogramm dargestellt wird.

Die folgende Abbildung zeigt den Graphen einer Sinusfunktion. Eine Periode ist beispielhaft eingezeichnet.

 
zeitlicher Verlauf einer Sinuswelle

Wenn eine periodische Funktion nicht zeitlich, sondern räumlich veränderlich ist, spricht man auch von einer Raumfrequenz. Diese ist entsprechend als die Anzahl der Perioden pro Raumeinheit definiert. Also zum Beispiel als Perioden pro mm, mit der Einheit:

 


 
räumlicher Verlauf einer Sinuswelle


Ein Beispiel aus der medizinisches Bildgebung für eine solchen Fluktuation ist das Bleibalkenmuster, welches häufig verwendet wird, um die räumliche Auflösung eines Abbildungssystems zu bestimmen (siehe folgende Abbildung). In diesem Fall wird die räumliche Auflösung in Linienpaare pro mm angegeben, wobei jeder Bleibalken mit seine angrenzenden Leerstelle als Linienpaar bezeichnet wird.


 
(a) Testbild zur Bestimmung der räumlichen Auflösung einer Gammakamera


 
(c) Plot der Zählrate eines Detektors aus a) oder b) in einem gut aufgelösten Bereich

Räumliche Auflösung

Bearbeiten

Eine Frage, die häufig im Rahmen der medizinischen Bildgebung auftritt, lautet: "Was ist die kleinste Struktur, die man mit einem gegebenen abbildenden System auflösen kann?" Die technische Größe, die diese Frage beantwortet, heißt räumliche Auflösung des abbildenden Systems. Man definiert sie je nach Zweck etwas unterschiedlich. Hier beschränken wir uns auf folgende einfache Definition: Die räumliche Auflösung ist der kürzeste Abstand, den zwei Linien haben dürfen, um noch als zwei getrennte und nicht als eine große verschmierte Linie gesehen werden zu können. Die Einheit der räumlichen Auflösung ist mm. Schauen wir uns zunächst ein Beispiel an. Die folgenden beiden Graphiken zeigen das gleiche Bild. Es besteht aus vertikalen Linien. Ihr Abstand und ihre Dicke nehmen nach rechts ab. Das Bild auf der linken Seite zeigt das Bild wie es von einem nahezu perfekten abbildenden System (mit sehr hoher räumlicher Auflösung) aufgenommen wurde. Man erkennt klar die einzelnen Linien, insbesondere auch am rechten Rand des Bildes. Das rechte Bild wurde mit einem qualitativ schlechteren abbildenden Verfahren aufgenommen. Man kann offenbar zum rechten Rand hin keine einzelnen Linien mehr ausmachen. In der Realität hat jedoch jedes abbildende System eine begrenzte räumliche Auflösung. Somit zeigt das rechte Bild einen realistischeren Fall.

 
(Nahezu) ideales abbildendes System
 
reales abbildendes System

Man erkennt weiterhin, dass es ziemlich schwierig ist, einen exakten Wert für die räumliche Auflösung eines abbildenden Systems anzugeben. Man muss hierzu den Bereich in der rechten Abbildung betrachten, in dem zwei benachbarte Linien sich gerade so weit überlappen, dass man sie nicht mehr als zwei getrennte Linien erkennen kann. Man erkennt sofort, dass dieser Effekt nicht bei einer bestimmten Linie auftritt, sondern über mehrere Linien allmählich zunimmt. Aber man sieht auch, dass man am linken Rand klar getrennte Linien und auf am rechten Rand überhaupt keine Linien mehr erkennen kann.

Schauen wir uns diese beiden Abbildungen nun in einer anderen Darstellung an. Die folgende Abbildung zeigt zwei Graphen. Der grüne Graph entspricht dem rechten und der rote dem linken Bild. Der obere Rand entspricht einem weißen Pixel und der untere Rand einem schwarzen Pixel. Zwischenwerte stehen entsprechend für graue Pixel. Der rote Graph würde normalerweise vertikale Linien enthalten, welche die horizontalen mit einander verbinden. Diese wurden jedoch, der Übersichtlichkeit wegen, nicht dargestellt. Die x Achse in dieser Abbildung entspricht den x Achsen in den beiden oberen Abbildungen. Somit enthalten beide Darstellungen dieselbe Information. Sie ist lediglich etwas anders aufgetragen.

 
Vergleich zwischen einem idealen und einem realen abbildenden System

Wir erkennen hier zwei Tatsachen wieder, die wir oben bereits bemerkt haben. Wir sehen jedoch auch, dass die Amplitude (der Abstand zwischen Minimum und Maximum in y Richtung) des Signals linear nach rechts abnimmt, bis die Funktion schließlich nicht mehr periodisch ist. Somit müssen wir zugeben, dass offensichtlich im rechten Teil fast alle Information verloren gegangen ist. Jedoch können wir die Information im periodischen Teil zurückgewinnen. Wir können die Amplituden entsprechend skalieren um die volle Amplitude wiederherzustellen. Leider skalieren wir jedoch auch das Rauschen. Somit funktioniert dieses Verfahren nur so lange, wie das Signal deutlich stärker ist als das Rauschen. Unten sieht man dieselbe Funktion mit einer in blau eingetragenen Einhüllenden.

 
reales abbildendes System mit blauer Einhüllender

Nun skalieren wir die Funktion und erhalten ein Bild mit höherer Auflösung. Wir halten den Kontrast auf der rechten Seite des Bildes konstant, aber wir erhöhen ihn um den Faktor 10 für die Pixel am rechten Rand und verwenden eine monotone Funktion für die dazwischenliegenden Pixel. Um die Berechnung technisch durchzuführen, benutzen wir die blaue einhüllende Funktion   und verwenden den Skalierungsfaktor  , welcher von der (Pixel-) Position abhängt. Unten sieht man zwei Bilder. Das rechte wurde mit der hier beschriebenen Methode verbessert.


 
reales abbildendes System
 
reales abbildendes System nach Kontrasterhöhung

Wir können auf dem rechten Bild deutlich mehr Linien erkennen als auf dem linken. Dennoch sehen wir am rechten Rand keine Linien mehr. Aber wir erkennen das wir durch die Skalierung des Kontrastes eine Verbesserung der räumlichen Auflösung erreichen konnten. Wir haben also die räumliche Auflösung durch ein rein rechnerisches Verfahren erhöht und dabei die Gerätschaften zur Bildaufnahme nicht verändert. Wir haben den Kontrast im rechten Teil des Bildes erhöht, jedoch war das Bild so beschaffen, dass die Raumfrequenz nach rechts zunimmt. Die bedeutet, dass wir tatsächlich die hochfrequenten Anteile des Bildes verstärkt haben. Und genau das ist die Idee bei der Bildoptimierung mit Hilfe von Fouriermethoden. Fouriermethoden erlauben es uns, selektiv bestimmte Raumfrequenzen zu verstärken und somit die räumliche Auflösung zu verbessern. Man muss jedoch aufpassen, nicht zu stark zu verstärken, weil man immer auch das Rauschen mit verstärkt und somit Gefahr läuft, ein stark verrauschtes Bild zu erhalten. Es ist wichtig, sich zu vergegenwärtigen, dass wir in diesem Beispiel noch keine Fouriermethoden verwendet haben. Wir haben lediglich ein Bild verwendet, was geschickt für unsere Zwecke konstruiert war, so das wir den Effekt der Verstärkung hoher Raumfrequenzen studieren konnten, ohne Fouriermethoden anwenden zu müssen. Das Schöne an wirklichen Fouriermethoden ist, dass man sie auf jedes Bild anwenden kann und sie meist viel besser funktionieren als der hier vorgeführte mathematische Taschenspielertrick.

Fourier-Reihen

Bearbeiten

Fast alle periodischen Funktionen, die man in der medizinischen Bildgebung braucht, können als Fourierreihe dargestellt werden. Dieser Ansatz basiert auf der Tatsache, dass jede Wellenform (mathematisch: jede periodische Funktion, die abschnittsweise stetig und monoton ist) als eine Reihe aus Sinus- und Kosinusfunktionen angesehen werden kann. Man kann sie dann wie folgt aufschreiben.

 

Wobei wir hier die Reihe nur bis zu dem Glied   bzw.   aufgeschrieben haben. Häufig hat man damit aber noch nicht genau die Funktion erreicht, die man erreichen wollte. Man hat lediglich eine sehr ähnliche Funktion. Selbst wenn man es nicht schafft, die exakte Funktion durch eine endliche Reihe darzustellen, so kann man immerhin eine Regel angeben, nach der man beliebig viele Glieder einer unendlichen Reihe berechnen kann. Mit jedem weiteren Glied nähert man sich immer besser an die gesuchte Funktion an. Ferner kann man sich beliebig genau an die gesuchte Funktion annähern, wenn man nur genügend viele Glieder berechnet. Man sagt dann auch, die unendliche Reihe ist die gesuchte Funktion. Diese Aussage ist mathematisch korrekt, allerdings kann man eine unendliche Reihe numerisch nicht genau berechnen, denn dazu müsste man ja unendliche lange rechnen. Wichtig ist nur, dass man eine Funktion durch ihre Fourierreihe beliebig genau berechnen kann, wenn man die Regel kennt, nach der die einzelnen Glieder der Reihe zu berechnen sind, und lange genug rechnet.


Eine Rechteckfunktion besitzt folgende Fourierreihe.

 

             (1)

Wobei   die Amplitude der Rechteckfunktion bezeichnet.

Auch die Gegenrichtung ist mathematisch möglich. z. B. kann eine Rechteckfunktion konstruiert werden indem, man eine große Zahl von Sinusschwingungen unterschiedlicher Amplitude aufaddiert. Die Summation der ersten 4 Terme ist in der folgenden Abbildung dargestellt. Der erste Term   ist in Abbildung (a) gezeigt die Addition des zweiten Terms in Abbildung (b) und so weiter. Man beachte, dass die Fundamentale, (auch erste Harmonische bezeichnet, Abbildung (a)) dieselbe Frequenz wie die Rechteckfunktion hat, und höhere Frequenzen die Form der Rechteckwelle nacheinander aufbauen (Abbildungen (b) bis (d)). Wir schließen daraus, dass höhere Frequenzen zur Schärfe der Flanken der Rechteckwelle beitragen.


 
(a) Erste Komponente
 
(b) Erste zwei (nicht null) Komponenten
 
(c) Erste drei (nicht null) Komponenten
 
(d) Erste vier (nicht null) Komponenten

Fourier-Spektrum

Bearbeiten

Die Fourierreihe kann auch als Frequenzspektrum dargestellt werden. Als Beispiel sind in der folgenden Abbildung die Amplituden der Frequenzkomponenten aus der Gleichung der Rechteckwelle (1) gegen die Raumfrequenz dargestellt. Man beachte, dass das Fourierspektrum dazu verwendet werden kann, diejenigen Frequenzen und Amplituden zu identifizieren, die zu einer gegebenen Wellenform beitragen. Man beachte weiterhin, dass Graphen, in denen die Amplitude gegen die Strecke aufgetragen wird, im allgemeinen als Darstellung im Ortsraum und Graphen, in denen die Amplitude gegen die Raumfrequenz aufgetragen ist, als Darstellung im Frequenzraum bezeichnet werden.

 
Fourierspektrum einer Rechteckwelle.

Fourier-Transformation

Bearbeiten

Die Fourier-Transformation ist eine elegante mathematische Methode, um Daten vom Ortsraum in den Frequenzraum zu konvertieren (siehe folgende Abbildung). Anders ausgedrückt, kann man die Frequenzen und Amplituden, aus denen eine beliebige Wellenform aufgebaut ist, leicht bestimmen, indem man die Fouriertransformierte dieser Wellenform berechnet.

 
Die Fourier-Transformation kann verwendet werden, um das Fourier-Spektrum jeder beliebigen Wellenform zu berechnen, wie hier am Beispiel der Rechteckfunktion gezeigt wird.

Die Fouriertransformation wird im Bereich der medizinischen Bildgebung weithin verwendet. Einige ihre Anwendungen sind:

  • Bestimmung der räumlichen Auflösung von abbildenden Systemen.
  • räumliche Lokalisierung in NMR (Kernspin) Bildgebung.
  • Analyse von Doppler-Ultraschall-Signalen.
  • Filterung von Bildern in der Transmissions- und Emissions-Tomographie.


Die Inverse Fourier-Transformation ist eine mathematische Methode um Daten in der umgekehrten Richtung zu konvertieren. Also vom Frequenzraum in den Ortsraum, wie in der folgenden Abbildung zu sehen ist:

 
Die inverse Fourier-Transformation kann verwendet werden um beliebige periodische Wellenformen beliebig genau anzunähern, wie hier am Beispiel der Rechteckfunktion gezeigt wird.

Zusammenfassend stellen wir fest, dass die Fouriertransformation FT uns erlaubt, die Sinuswellen zu bestimmen, aus denen eine Wellenform aufgebaut ist und dass wir mit der inversen Fouriertransformation IFT eine Wellenform aus einzelnen Sinuswellen zusammen bauen können.

Schließlich sollte noch erwähnt werden, dass die Berechnung von Fouriertransformationen mittels digitaler Computer meist mit einem speziellen Algorithmus namens Fast Fourier Transformation (FFT) durchgeführt wird.

Die Dirac-Delta Funktion

Bearbeiten

Ein für die medizinische Bildgebung wichtiger und interessanter Fall ist die Diracsche Delta-Funktion  . (Mathematisch korrekt müsste man sie eigentlich Delta- Distribution nennen, weil sie die mathematischen Eigenschaften, die es braucht, um sich Funktion nennen zu dürfen, nicht ganz erfüllt, aber außerhalb der Mathematik nimmt man es damit meist nicht so genau.) Sie hat überall den Wert 0, außer an der Stelle  . Dort kennt man ihren Wert nicht wirklich. Ferner gilt die Normierungsbedingung:

 .

Ihre Fourier-Transformierte ist in der folgenden Abbildung gezeigt:

 
Die Delta-Funktion (links) und ihre Fourier Transformierte (rechts)

Die Fouriertransformation der Delta-Funktion ist eine konstante Funktion. Die Fourierreihe der Delta-Funktion ist also die Summe aus unendlich vielen Sinus-Funktionen, die alle die gleiche Amplitude haben. Wenn wir nun anfangen, die Delta Funktion zu verbreitern (also etwa durch eine schmale Gauß-Funktion, die der Normalbedingung genügt, ersetzen), sehen wir, dass die Sinusfunktionen mit niedrigeren Frequenzen hohe Amplituden haben und dass die Amplituden mit zunehmender Raumfrequenz abnehmen.

 
Die Fourier Transformierte einer verbreiterten Delta Funktion.

Modulationsübertragungsfunktion

Bearbeiten

Diese Verbreiterung ist ein Phänomen sehr ähnlich denen, die in der medizinischen Bildgebung auftreten. Wir können annehmen, dass die vorhergehende Grafik in der Amplitude gegen Strecke folgenden Effekten ähnlich ist:

  • Dichtprofil einer Abbildung eines kleinen Lochs in einer Bleiplatte.
  • Zählratenprofil durch ein Bild einer radioaktiven Punktquelle.

Die in dieser Grafik dargestellte Funktion wird als Point Spread Function (PSF, auf Deutsch Punktantwortfunktion, jedoch meist zugunsten des englischen Begriffs nicht verwendet) bezeichnet. Ihre Fouriertransformierte wird Modulationsübertragungsfunktion MTF genannt.

 
Darstellung der MTF eines idealen und eines realen Abbildenden Systems.

Einen Vergleich von allgemeinen Begriffen und Begriffen der Bildgebung, die in diesem Themenkreis verwendet werden, ist in der folgenden Tabelle gegeben.

Gebiet Allgemeiner Begriff Begriff in der Bildgebung
Ortsraum Ausgangsfunktion Punktantwortfunktion (PSF)
Frequenzraum Fourierspektrum Modulationsübertragungsfunktion (MTF)


Weiterhin kann man die MTF bestimmen aus:

  • Linienantwortfunktion (Line-Spread-Function LSF)
  • der differenzierten Eckenantwortfunktion (Edge-Response-Function ERF).

Symbolisch können wir schreiben:

 

wobei   die mathematische Faltungsoperation bezeichnet. Anders ausgedrückt erhält man das reale Bild indem man das ideale Bild mit der PSF des Abbildungssystems faltet.

Um das ideale Bild zurückzuerhalten, muss man den Effekt der PSF auf rechnerische Weise rückgängig machen. Dies kann man leicht erreichen, indem man Fouriertransformation verwendet, weil der Faltungsprozess zu einer Multiplikation im Frequenzraum äquivalent ist. Also:

 

wobei   die Fouriertransformierte der Ortsraumfunktion   bezeichnet. Damit hat man:

 

.

Der vollständige Rekonstruktionsprozess, der auch als Entfaltung bezeichnet wird, ist durch die folgende Gleichung gegeben:

 

wobei   die inverse Fouriertransformation der Frequenzraumfunktion   bezeichnet.

Filterung des Fourierspektrums

Bearbeiten

Die oben vorgestellte Methode zur Rückgewinnung des idealen Bildes aus dem realen Bild ist ein Beispiel für die Filterung des Fourierspektrums. Anders ausgedrückt kann das Fourierspektrum, sobald es einmal für ein Bild erzeugt worden ist, gefiltert werden, so dass bestimmte Raumfrequenzen modifiziert, d.h. verstärkt oder unterdrückt werden. Das gefilterte Spektrum kann dann zurücktransformiert werden um daraus das gefilterte Bild zu gewinnen, und es so zu schärfen oder zu glätten. Dieser Prozess ist in der folgenden Abbildung schematisch dargestellt.

 
Die Filterung des Fourierspektrums

Eine Sache müssen wir noch im Detail klären, bevor wir uns näher mit den Bilddaten im Frequenzraum beschäftigen. Wir erinnern uns, dass Bilder im allgemeinen digital in einer rechteckigen Punktmatrix aufgenommen werden. Wobei die Größe dieser Punkte bestimmt, wie genau das digitale Bild seinem analogen Gegenstück entspricht. Die sich hieraus ergebende digitale räumliche Auflösung führt zu einer Grenze für die maximale Raumfrequenz die noch aufgenommen werden kann. Das üblicherweise in der digitalen Bildgebung verwendete Kriterium basiert auf dem Nyquist-Shannon-Abtasttheorem welches besagt, dass:

Wenn ein Bild Raumfrequenzkomponenten bis zu einer maximalen Raumfrequenz f enthält, dann sollten die Bilddaten mit einer Abtastfrequenz die mindestens doppelt so hoch wie f ist abgetastet werden um eine treue Wiedergabe zu erreichen.

Diese Abtastfrequenz wird üblicherweise als Nyquist-Frequenz bezeichnet. Bei geringeren Abtastfrequenzen, können die entstehenden digitalen Bilder Artefakte enthalten, die als Moiré-Effekt bezeichnet werden. Man bezeichnet dies gelegentlich auch als Alias-Effekt.

Gefilterte Rückprojektion

Bearbeiten

Der einfacheste Rückprojektionsprozess führt zu einer Unschärfe die das Bild so aussehen lässt als wäre es mit einer   Funktion verschmiert worden. In der gefilterten Rückprojektion werden Fourierfilter benutzt um den Effekt der   Verwischung zu entfernen.

In Formeln, kann man die gemessene Projektion als einer Faltung mit der Schmierfunktion wir folgt aufschreiben:


 .

Die erste Stufe des Filterungsprozesses besteht in der Fouriertransformation der gemessenen projizierten Daten. D.h.

 

.

Die korrigierte Projektion,  , erhält man dann indem man die gemessene Projektion

durch   teilt und die inverse Fouriertransformation berechnet:

 

.

Der Ausdruck   ist einfach eine steigende Ursprungsgerade (Diese wird im Englisch als Ramp bezeichnet).

Weiterhin kann man, wenn man die Funktion   variiert, gleichzeitig die Unschärfe korrigieren und bestimmte Eigenschaften des rückprojizierten Bildes unterdrücken oder verstärken. Zum Beispiel kann man das Bild schärfen und gleichzeitig:

  • feine Details hervorheben (wie im sogenannten Knochenalgorithmus in der Röntgencomputertomographie)
  • das Bildrauschen unterdrücken (wie im sogenannten Weichgewebealgorithmus in der Röntgencomputertomographie)

Man kann die Ramp-Funktion variieren in dem man sie mit einer zweiten Funktion multipliziert, z.B. mit einer Butterworth Funktion wie sie beim SPECT verwendet wird:

 


wobei:

  •   die Cutoff Frequenz bezeichnet, welche die Frequenz definiert, bei welcher der Filter nur noch 50% des Signals durchlässt.
  •  : die Ordnung des Filters definiert (je höher die Ordnung, um so steiler fällt die Durchlasskurve des Filters im Bereich der Cutoff-Frequenz ab).
 
Ramp- und Butterworthfilter (Fünfte Ordnung 60% Cutoff)
 
Ramp- und Butterworthfilter (Fünfte Ordnung 40% Cutoff)
 
Ramp- und Butterworthfilter (Dritte Ordnung 40% Cutoff)

Die Eigenschaften einer Anzahl von Filtern, die in der SPECT-Bildgebung verwendet werden, sind in der folgenden Tabelle angegeben:


Filter
Gleichung
Kommentar
Ram-Lak
Amplitude = 1 für f < fc;

Amplitude = 0 für f > fc
Entfernung von Sternartefakten; Rauschempfindlichkeit
Butterworth
 
Rauschunterdrückung
Metz
 
Rauschunterdrückung, Kontrasterhöhung
Wiener
 
Rauschunterdrückung, Kontrasterhöhung
Scramp
  für  ;

Amplitude = 0 für f > fc
Rauschunterdrückung, Kontrasterhöhung
Inverse MTF
  für f < fc;

Amplitude = 0 für f > fc
Rauschunterdrückung, Kontrastunterdrückung
Hamming
  für  ;

Amplitude = 0 für f > fc
Rauschunterdrückung
Parzen
 für f < fc;

Amplitude = 0 für f > fc
Rauschunterdrückung
Shepp-Logan
 für f < fc;

Amplitude = 0 für f > fc
Rauschunterdrückung
Hann
  für  ;

Amplitude = 0 für f > fc
Rauschunterdrückung
Man beachte, dass das x in der Gleichung für den Metz-Filter Leistungsfaktor und das S in der Wiener-Filter-Gleichung Formfaktor genannt wird.

Bei der Auswahl des Filters für eine gegeben Rekonstruktionsaufgabe muss man im allgemeinen einen Kompromiss zwischen der Rauschunterdrückung und der Verstärkung feiner Details finden. Wobei man das Raumfrequenzspektrum der interessierenden Bilddaten berücksichtigen muss. In einigen Fällen ist ferner eine Erhöhung des Kontrastes notwendig und muss somit zusätzlich zu den oben genannten Anforderungen berücksichtigt werden.