Gefilterte Rückprojektion

Die gefilterte Rückprojektion (auch FBP für filtered back projection) ist ein auf der Radon-Transformation beruhendes Verfahren zur Bildrekonstruktion, das in erster Linie in der Computertomographie verwendet wird. Die FBP rekonstruiert aus einem Satz eindimensionaler Projektionen verschiedener Richtung ein ursprüngliches 2D-Bild. Hierzu werden die Projektionen zuerst gefiltert und dann in der jeweiligen Richtung über die Bildfläche gewischt („rückprojiziert“). Das Verfahren hat den großen Vorteil, dass es schnell ist, da es wenig Rechenleistung benötigt. In der SPECT sowie in der PET wurde sie mittlerweile von den iterativen Rekonstruktionsverfahren verdrängt.

Überblick

Aufnahme von Projektionen mittels Durchleuchtung mit Röntgenstrahlen.

In zahlreichen tomografischen Anwendungen kann die Projektion eines Objekts in verschiedene Richtungen gemessen werden, nicht aber direkt das Innere des Objektes. Beispiele hierfür sind die Durchleuchtung von Objekten mittels Röntgenstrahlung, Neutronen oder Ultraschall, oder auch projektive Messungen in der Quantenmechanik. Die Menge der Projektionen in alle Richtungen enthält die vollständige Information über das Innere des Objektes. Sie ist allerdings nur schwierig direkt interpretierbar, da jeder Punkt aus einer Überlagerung von Eigenschaften besteht, die beim Projizieren aufsummiert wurden. Originalbild und Projektionen sind durch eine bijektive Transformation miteinander verknüpft, die mathematisch durch die Radon-Transformation beschrieben wird. Die gefilterte Rückprojektion ist eine Implementierung der Transformation, die es erlaubt, aus einer Menge von Projektionen das Originalbild zu errechnen.

Algorithmus

Schematische Darstellung der gefilterten Rückprojektion einer Kreisscheibe.

Die gefilterte Rückprojektion wird auf diskreten Daten durchgeführt, also Pixelbildern der Projektionen für eine endliche Anzahl an gleichmäßig verteilten (typischerweise einigen hundert) Projektionswinkeln. Die eindimensionalen Projektionen können so ebenfalls zu einem zweidimensionalen Bild, dem sogenannten Sinogramm, zusammengefasst werden.

Die gefilterte Rückprojektion basiert auf der einfachen Idee, dass die Projektionen ähnlich, wie sie aus dem Bild herausprojiziert wurden, auch wieder zurückprojiziert werden können. Die naive Implementierung der ungefilterten Rückprojektion funktioniert aber nicht, weil dabei jeder Punkt der Projektion über die gesamte Bildfläche verschmiert wird anstatt nur an seine ursprüngliche 2D-Position. Man würde ein Bild erhalten, in dem jeder Originalpunkt mit einer Punktspreizfunktion der Form 1/|r| verbreitert wird und das somit einer sehr unscharfen Version des Originalbildes entspricht.

Den funktionierenden Algorithmus erhält man aus der inversen Radontransformation. Sei $f(x,y)$ das Originalbild und $p_{\phi }(z)$ die Projektion in Richtung des Winkels $\phi$ . Dann ist^[1]

f(x,y)=\int _{0}^{\pi }\left(\int _{-\infty }^{\infty }p_{\phi }(z')\cdot g(z-z')\,\mathrm {d} z'\right)\,\mathrm {d} \phi

.

Hierbei wurden die Koordinaten $x=z\sin(\phi )$ und $y=z\cos(\phi )$ in Polarkoordinaten ausgedrückt. Das Integral über $\phi$ drückt aus, dass über die Beiträge aller Winkel summiert werden muss. Das Integral über $z'$ ist eine Faltung der Projektion $p_{\phi }(z)$ mit einem geeigneten Hochpassfilter $g(z)$ , die dem Verfahren ihren Namen verleiht. Tatsächlich ist die inverse Radontransfomation ein schlecht gestelltes Problem, die einen irregulären Filterkern $g(z)$ verlangt. Um einen diskreten Filterkern zu erhalten, muss eine leichte Weichzeichnung als Regularisierung in Kauf genommen werden, die je nach Parameter zu unterschiedlich breitem $g(z)$ führt.

In der Praxis wird die Filterung als schnelle Faltung im Fourierraum implementiert, wo sie eine besonders elegante Form annimmt. Dort nämlich wird die (diskrete) Fouriertransformierte ${\tilde {p}}_{\phi }(k)={\mathcal {F}}p_{\phi }(z)$ mit dem Filter im Fourierraum ${\tilde {g}}(k)={\mathcal {F}}g(z)$ multipliziert, der dort die besonders einfache Form

{\tilde {g}}(k)=|k|

annimmt. Alle Frequenzkomponenten werden also proportional zum Betrag ihrer Frequenz gewichtet. Zur Vermeidung von starkem Rauschen im Ergebnis wird auch hier ein zusätzlicher Tiefpassfilter in Form einer Fensterfunktion $w(k)$ angewendet, der die ganz hohen Frequenzen unterdrückt. Für die Auswahl geeigneter Fenster hat die Literatur dutzende verschiedene Beispiele hervorgebracht, zum Beispiel Hamming, Hanning oder Blackman. Außerdem muss man bei der DC-Komponente $k=0$ aufpassen, die den Mittelwert von $g$ bestimmt. Für diskrete Frequenzbins der Breite $\Delta k$ muss ${\tilde {g}}(0)$ dem Bin-Mittelwert von $|k|$ entsprechen, also ${\tilde {g}}(0)=\Delta k/4$ anstatt 0. Unter Anwendung der Fouriertransformation ${\mathcal {F}}$ und ihrer Inversen ${\mathcal {F}}^{-1}$ erhält man

f(x,y)=\int _{0}^{\pi }{\mathcal {F}}^{-1}{\bigl (}{\mathcal {F}}(p_{\phi })(k)\cdot |k|\cdot w(|k|){\bigr )}\,\mathrm {d} \phi

.

Die Implementierung der Fouriertransformation und ihrer Inversen erfolgt effizient als FFT.

Die Schritte der gefilterten Rückprojektion sind also:

Transformiere jede Projektion $p_{\phi }(z)$ in den Fourierraum.
Multipliziere die Transformierte ${\tilde {p}}_{\phi }(k)$ mit ${\tilde {g}}(k)=|k|$ und einer Fensterfunktion $w(k)$ . Setze ${\tilde {g}}(0)=\Delta k/4$ .
Rücktransformation in den Ortsraum ergibt die gefilterte Projektion $p'_{\phi }(z)$ .
Rückprojektion von $p'_{\phi }(z)$ in die Bildebene von $f(x,y)$ , indem die 1D-Projektion in eine zweite Dimension gestreckt, dann um den Winkel $-\phi$ gedreht und schließlich zum Bild addiert wird.

Die Summe der Rückprojektionen aller Winkel ergibt das Bild $f(x,y)$ .

Alternativmethode Gefiltertes Schichtgramm

Eine alternative Methode der gefilterten Rückprojektion ist das Gefilterte Schichtgramm („filtered layergram“),^[1] wo die Reihenfolge der beiden Schritte Filterung und Rückprojektion vertauscht wird. Dabei werden also die ungefilterten Bilder zuerst rückprojiziert und daraufhin das 2D-Bild gefiltert. Das ist deshalb möglich, weil beide Schritte lineare Funktionen auf den Messdaten darstellen und daher mathematisch vertauscht werden können. Aufgrund der diskreten Rasterung der Daten ergeben sich aber leichte Unterschiede in den Ergebnissen, zum Beispiel im Einfluss von Bildrauschen und Diskretisierungsfehler. In der Praxis wird die gefilterte Rückprojektion immer der Schichtgramm-Methode vorgezogen, weil dabei der rechenintensivste Schritt, die Fouriertransformation, auf jedem 1D-Bild einzeln ausgeführt werden kann, sobald dieses gemessen wurde, und somit nicht abgewartet werden muss, bis alle Projektionen aufgenommen sind.

Beispiel

Ein Beispiel soll den Algorithmus verdeutlichen:

Man stelle sich ein von vorne beleuchtetes Aquarium vor, bei dem sich die Konturen der Fische an einer dahinter befindlichen Leinwand abbilden. Führt man Lichtquelle und Leinwand in gleichen Winkelschritten um das Aquarium herum, erhält man viele Projektionen der Fische auf der Leinwand. Um aus diesen Daten hinterher die Position der Fische im Aquarium zu bestimmen, nimmt man jede einzelne der aufgenommenen Projektionen und projiziert sie auf das Volumen des Aquariums zurück (daher der Name des Verfahrens). Es ist klar, dass hierbei die Tiefeninformation nicht berücksichtigt wird, d. h., das endlich tiefe Bild des Fisches auf der Leinwand wird in Projektionsrichtung über das Bild verschmiert. Dieser Fehler lässt sich jedoch durch Anwendung eines geeigneten Bildfilters bei der Reprojektion wirkungsvoll unterdrücken.

Die Punktspreizfunktion der ungefilterten Rückprojektion ist $1/r$ , wobei $r$ der Betrag im Ortsraum ist; das bedeutet, wenn das abzubildende Objekt nur aus einem Punkt mit den Koordinaten $\left(x',y',z'\right)$ besteht (Delta-Distribution), so erhält die ungefilterte Rückprojektion als Bild ein Signal am Ort $\left(x,y,z\right)$ , das proportional zu ${\frac {1}{\sqrt {(x-x')^{2}+(y-y')^{2}+(z-z')^{2}}}}$ ist. Die „Filterung“ entspricht mathematisch einer Faltung. Mit Hilfe des Faltungssatzes kann entfaltet werden, indem das rückprojizierte Bild in den Fourierraum transformiert, mit $|k|={\sqrt {k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}}$ , dem Betrag im Fourierraum, multipliziert und anschließend wieder zurück in den Ortsraum transformiert wird. Gemäß dem Abtasttheorem kann der Filter bei einer bestimmten Raumfrequenz abgeschnitten werden. Außerdem ist zu beachten, dass Daten einer Computertomographie diskret vorliegen und nicht, wie mathematisch eigentlich nötig, kontinuierlich. Deshalb ist das CT nicht „exakt“ und es gibt nicht den idealen Filter, da beispielsweise zwischen den einzelnen diskreten Punkten gemittelt (Shepp-Logan-Filter) oder nicht gemittelt (Ram-Lak-Filter) werden kann. Je nachdem, welchen Filter man benutzt, wird das Bild entweder kontrastreicher, aber verrauschter, oder kontrastärmer, aber rauschreduzierter.