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 relativ schnell ist und wenig Rechenleistung benötigt. In der SPECT sowie in der PET wurde sie mittlerweile von den iterativen Rekonstruktionsverfahren verdrängt, bildet aber nach wie vor das Rückgrat vieler Rekonstruktionsanwendungen.

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.

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.^[1]

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

${\displaystyle 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 ${\displaystyle x=z\sin(\phi )}$ und ${\displaystyle y=z\cos(\phi )}$ in Polarkoordinaten ausgedrückt. Das Integral über ${\displaystyle \phi }$ drückt aus, dass über die Beiträge aller Winkel summiert werden muss.

Das Integral über ${\displaystyle z'}$ ist eine Faltung der Projektion ${\displaystyle p_{\phi }(z)}$ mit einem geeigneten Hochpassfilter ${\displaystyle g(z)}$, die dem Verfahren ihren Namen verleiht. Tatsächlich ist die inverse Radontransfomation ein schlecht gestelltes Problem, die einen irregulären Filterkern ${\displaystyle 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 ${\displaystyle g(z)}$ führt.

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

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

annimmt. Alle Frequenzkomponenten werden also proportional zum Betrag ihrer Frequenz gewichtet. In der Praxis weicht der Filterkern von diesem idealen Rampenfilter ab, weil bei einer endlichen Faltungslänge der ideale Filterkern ${\displaystyle g(z)}$ zwangsweise abgeschnitten wird, was alle Frequenzkomponenten der Diskreten Fourier-Transformation verändert. Am stärksten ist die DC-Komponente ${\displaystyle {\tilde {g}}(0)}$ betroffen, die dann anders als bei der idealen Rampe nicht exakt Null wird (Siehe Grafik rechts). Der so angepasste Filterkern mit ansonsten linearem Frequenzverlauf wird als Ram-Lak Filter (benannt nach Ramachandran und Lakshminarayanan) bezeichnet und nimmt auf diskreten Sampeln die Form

${\displaystyle g(z)_{\text{Ram-Lak}}={\begin{cases}1/4&z=0\\-1/(\pi ^{2}z^{2})&z\ \mathrm {ungerade} \\0&\mathrm {ansonsten} \end{cases}}}$

an.^[3]

Zur Vermeidung von starkem Rauschen im Ergebnis kann ein zusätzlicher Tiefpassfilter in Form einer Fensterfunktion ${\displaystyle w(k)}$ angewendet werden, der die ganz hohen Frequenzen unterdrückt. Eine populäre Wahl ist der Shepp-Logan Filter, der die hohen Frequenzen mit einer Sinc-Funktion gerade stark genug heruntergewichtet, um Oszillationen im Filterkern zu verhindern. Fensterfunktionen können auch individuell angepasst werden, um einen beliebigen gewünschten Schärfeeindruck des Bildes zu erhalten. Unter Anwendung der Fouriertransformation ${\displaystyle {\mathcal {F}}}$ und ihrer Inversen ${\displaystyle {\mathcal {F}}^{-1}}$ erhält man

${\displaystyle 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:

1. Transformiere jede Projektion ${\displaystyle p_{\phi }(z)}$ in den Fourierraum.
2. Multipliziere die Transformierte ${\displaystyle {\tilde {p}}_{\phi }(k)}$ mit ${\displaystyle {\tilde {g}}(k)}$, die aus der Fouriertransformierten des Filterkerns und ggf. einer Fensterfunktion ${\displaystyle w(k)}$ vorberechnet wurde. (Filterung im Frequenzraum)
3. Rücktransformation in den Ortsraum ergibt die gefilterte Projektion ${\displaystyle p'_{\phi }(z)}$.
4. Rückprojektion von ${\displaystyle p'_{\phi }(z)}$ in die Bildebene von ${\displaystyle f(x,y)}$, indem die 1D-Projektion in eine zweite Dimension gestreckt, dann um den Winkel ${\displaystyle -\phi }$ gedreht und schließlich zum Bild addiert wird.

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

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 ${\displaystyle 1/r}$, wobei ${\displaystyle r}$ der Betrag im Ortsraum ist; das bedeutet, wenn das abzubildende Objekt nur aus einem Punkt mit den Koordinaten ${\displaystyle \left(x',y',z'\right)}$ besteht (Delta-Distribution), so erhält die ungefilterte Rückprojektion als Bild ein Signal am Ort ${\displaystyle \left(x,y,z\right)}$, das proportional zu ${\displaystyle {\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 ${\displaystyle |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.

Die gefilterte Rückprojektion ist ein relativ effizientes Verfahren, das auf moderner Hardware in Echtzeit ausgeführt werden kann. Die Laufzeit skaliert sowohl direkt mit der Größe der Ein- und Ausgangsdaten, als auch aufgrund der verwendeten FFT mit dem Logarithmus der Projektionslänge. Direkte Fourier-Methoden basierend auf dem Zentralschnitt-Theorem können theoretisch eine bessere Laufzeit als die FBP erreichen, da nicht jeder Pixel unabhängig aus allen Projektionen berechnet werden muss, sondern das gesamte Bild in einer kombinierten Berechnung entsteht. Bei den häufig niedrigen Bildauflösungen ist der Vorteil aber gering. Iterative und algebraische Methoden wie z. B. die „Algebraische Rekonstruktions-Technik“ ART sind in der Regel deutlich rechenintensiver als die FBP und eignen sich daher oft nur für Spezialanwendungen.

In praktischen Rekonstruktions-Anwendungen ist die Geometrie oft nicht exakt in zweidimensionale Schichten aufteilbar, wo der Algorithmus direkt angewendet werden kann. In der Computertomographie beispielsweise wurden in den 1990er Jahren Mehrzeilen-CT-Geräte eingeführt, wo die äußeren Zeilen schräg unter einem Kegelwinkel aufgenommen werden, und die Projektionen verschiedener Richtungen einander überschneiden. Für solche Fälle wurden zahlreiche Varianten entwickelt, welche die komplexere Geometrie exakt oder näherungsweise berücksichtigen.

Eine alternative Methode der gefilterten Rückprojektion ist das Gefilterte Schichtgramm („filtered layergram“),^[2] 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 in der Regel der Schichtgramm-Methode vorgezogen, weil dabei ein rechenintensiver 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.

Der häufig eingesetzte FDK-Algorithmus (benannt nach den Autoren Feldkamp, Davis und Kress) erweitert die gefilterte Rückprojektionen auf in axialer Richtung ausgedehnte 3D-Volumina, wenn das Objekt von einer punktförmigen Strahlungsquelle einmal kreisförmig umfahren wird. Genau wie bei der zweidimensionalen gefilterte Rückprojektion werden auch hier die Projektionen in azimutaler Richtung, also entlang der Drehrichtung, gefiltert. Die Rückprojektion erfolgt dann allerdings nicht zeilenweise in einzelne Schichten. Stattdessen werden für jeden Voxel im 3D-Volumen aus jeder Projektion die (bereits gefilterten) Strahlen herausinterpoliert, welche den jeweiligen Voxel schneiden würden. Das Verfahren erzeugt hierdurch ein angenähertes Rekonstruktionsbild, das im Fall von kleinen Kegelwinkeln nahe an eine exakte Rekonstruktion herankommt, und für größere Kegelwinkel stärker werdende sogenannte Kegelstrahl-Artefakte aufweist.^[4][5]

Basierend auf der Erkenntnis, dass der Rampenfilter in eine Abfolge aus einer Ableitung und einer Hilbert-Transformation zerlegt werden kann, wurde eine neue Klasse von gefilterten Rückprojektionsalgorithmen entwickelt. Die Ableitung hat den Vorteil, dass sie lokal ist, und deshalb im Gegensatz zur normalen gefilterten Rückprojektionen in bestimmten Situationen auch ohne vollständig gemessenen Daten korrekte Rekonstruktionen erzeugen kann. Da die Ableitung und die Hilbert-Transformation während einer Rekonstruktion auch in verschiedene Richtungen ausgeführt werden können, wurde es hiermit möglich, in bestimmten Geometrien auch Daten, die mit großen Kegelwinkeln gemessen wurden, exakt zu rekonstruieren, obwohl keine 2D-Ebene existiert, in der vollständige Daten aus allen Richtungen für die klassische FBP vorliegen.^[6][5]

Beispiele solcher Algorithmen sind die von Tuy^[7] (1983), Smith^[8] (1985), Grangeat^[9] (1987) und Katsevich^[10] (2002). Bis in die 2020er Jahre wurde eine Vielzahl an Varianten solcher Algorithmen entwickelt, die zunehmend effizienter wurden. Aufgrund spezieller Voraussetzungen und ihrer Komplexität konnten sie allerdings die herkömmliche FBP bislang in praktischen Anwendungen nicht ersetzen.^[11]

• Dekonvolution

• Universität Giessen, Klinik für Nuklearmedizin (Memento vom 8. Mai 2006 im Internet Archive)
• Universität Freiburg (PDF-Datei; 1,44 MB) (Die Seite existiert leider nicht)
• Schattenspiele. Uni Bonn (Memento vom 15. Januar 2013 im Internet Archive)
• W. Schlegel, J.Bille (Hrsg.): Medizinische Physik 2. Springer, Berlin Heidelberg New York 2002, ISBN 3-540-65254-X

1. ↑ Richard Gordon, Gabort Herman: Three-Dimensional Reconstruction from Projections: A Review of Algorithms. In: International Review of Cytology. 38. Jahrgang, 1974, S. 111–151, doi:10.1016/S0074-7696(08)60925-0 (englisch). 
2. ↑ ^{a b} Thorsten M. Buzug: Computed Tomography. Springer Berlin Heidelberg 2010, ISBN 978-3-642-07257-4.
3. ↑ Azriel Rosenfeld, Avinash C. Kak: Digital image processing and analysis. 1982, ISBN 978-0-12-597301-4 (englisch). 
4. ↑ L. A. Feldkamp, L. C. Davis, and J. W. Kress: Practical cone-beam algorithm. In: Journal of the Optical Society of America A. 1. Jahrgang, Nr. 6, 1984, S. 612–619, doi:10.1364/JOSAA.1.000612 (englisch, optica.org). 
5. ↑ ^{a b} Henrik Turbell: Cone-Beam Reconstruction Using Filtered Backprojection. In: Doktorarbeit Uni Linköping. Februar 2001; abgerufen im 1. Januar 1 (englisch). 
6. ↑ Rolf Clackdoyle, Michel Defrise: Tomographic reconstruction in the 21st century. In: IEEE Signal Processing Magazine. 27. Jahrgang, Nr. 4, 14. Juni 2010, S. 60–80, doi:10.1109/MSP.2010.936743 (englisch, ieee.org). 
7. ↑ Heang K. Tuy: An inversion formula for cone-beam reconstruction. In: SIAM Journal on Applied Mathematics. 43. Jahrgang, Nr. 3, 1983, S. 546–552, doi:10.1137/0143035 (englisch, purdue.edu [PDF]). 
8. ↑ Bruce D. Smith: Image reconstruction from cone-beam projections: necessary and sufficient conditions and reconstruction methods. In: IEEE transactions on medical imaging. 4. Jahrgang, Nr. 1, März 1985, S. 14–25, doi:10.1109/TMI.1985.4307689 (englisch, mit.edu [PDF]). 
9. ↑ Pierre Grangeat: Mathematical framework of cone beam 3D reconstruction via the first derivative of the radon transform. In: Medical Imaging Techniques. 1497. Jahrgang, 5. Juni 1990, S. 66–97, doi:10.1007/BFb0084509 (englisch). 
10. ↑ Alexander Katsevich: Theoretically Exact Filtered Backprojection-Type Inversion Algorithm for Spiral CT. In: SIAM Journal on Applied Mathematics. 62. Jahrgang, Nr. 6, 2002, S. 2012–2026, doi:10.1137/S0036139901387186 (englisch, ucf.edu). 
11. ↑ Xiaochuan Pan, Emil Y Sidky, Michael Vannier: Why do commercial CT scanners still employ traditional, filtered back-projection for image reconstruction? In: Inverse Problems. 25. Jahrgang, Nr. 12, 1. Dezember 2009, ISSN 0266-5611, S. 123009, doi:10.1088/0266-5611/25/12/123009 (englisch).