Bildrekonstruktionsverfahren in der Computer Tomographie

Anthoula Papadopoulou, 05.03.2015

Unter dem Begriff der Bildrekonstruktion wird in der Computertomographie der Prozess der Transformation von Sensor-Daten (Messdaten) in aussagekräftigen Bilder eines Objektes bezeichnet. Das Faszinierende dabei sind nicht nur die physikalisch-technischen Prinzipien, sondern vielmehr die mathematischen Methoden mit deren Hilfe eine räumliche Verteilung aus ihren Projektionsaufnahmen rekonstruiert werden kann. Dabei gibt es eine Vielzahl von Varianten für die Modellierung, die sich durch unterschiedlichen Schwierigkeitsstufen unterscheiden lassen. ^[1]

Grundlagen der Messung

Im Rahmen einer Computertomographie wird die Intensität $//<![CDATA[ \begin{array}{l}I_0\end{array} //]]>$ bzw. die Schwächung der Röntgenstrahlung beim Durchgang durch das untersuchte Objekt gemessen. Daraus lässt sich der Projektionswert $//<![CDATA[ \begin{array}{l}P\end{array} //]]>$ errechnen, wodurch die Struktur des Objektes mithilfe der Dichteverteilung dargestellt werden kann (Abb.1). ^[2]

Die Messung sowie die Auswertung der Daten beruhen auf folgenden Annahmen: ^[3]

• Die Stahlen $//<![CDATA[ \begin{array}{l}L\end{array} //]]>$ sind Geraden und die Strahlenergie ist monoenergetisch.

• Die Dichtefunktion $//<![CDATA[ \begin{array}{l}f\end{array} //]]>$ beschreibt die Schwächungsverteilung des Objektes. Daraus lässt sich die Objektstruktur darstellen.

• Das beschränkte räumliche Objekt wird durch eine Folge von ebenen Parallelschnitten ersetzt.

Für die Intensität gilt ein exponentielles Schwächungsgesetz. Bei konstanter Dichte $//<![CDATA[ \begin{array}{l}f\end{array} //]]>$ und einer Weglänge $//<![CDATA[ \begin{array}{l}l\end{array} //]]>$ im Schnitt des Objektes ergibt sich ein reduzierter Wert der Anfangsintensität:

$//<![CDATA[ \begin{array}{l}I = I_0 \ exp(-f * l)\end{array} //]]>$   (1)

Bei inhomogenen Objekten, die eine Variation der Dichte in einem Schnitt des Objektes aufweisen, kann die Schwächung als Summe, beziehungsweise als Integral wie folgt beschrieben werden:

$//<![CDATA[ \begin{array}{l}I = I_0 \ exp(-\int_L \! f(x,y) \, \mathrm{d}l)\end{array} //]]>$   (2)

Wird das Integral isoliert, so ergibt sich die folgende Beziehung zur Darstellung der Projektionswerte $//<![CDATA[ \begin{array}{l}P\end{array} //]]>$:

$//<![CDATA[ \begin{array}{l}\int_L \! f(x,y) \, \mathrm{d}l = ln \frac{I_0}{I} = p(L)\end{array} //]]>$   (3)

Aus diesem Wert sollte zunächst der lineare Schwächungskoeffizient $//<![CDATA[ \begin{array}{l}\mu\end{array} //]]>$ für jeden Punkt der Fläche bestimmt werden. Das CT-Bild ist nämlich eine Darstellung der räumlichen Verteilung des linearen Schwächungskoeffizienten $//<![CDATA[ \begin{array}{l}\mu (x, y)\end{array} //]]>$ bzw. $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$. Die errechneten Schwächungskoeffizienten werden dabei über sogenannten CT-Zahlen relativ zur Schwächung von Wasser angegeben. Diese Zahlen werden auch als Hounsfield Units bezeichnet. ^[2]

Die Radon-Transformation

Bei inhomogenen Objekten ist es prinzipiell unmöglich die räumliche Verteilung des linearen Schwächungskoeffizienten $//<![CDATA[ \begin{array}{l}\mu (x, y)\end{array} //]]>$ direkt zu bestimmen. Die Information über die gesuchte Verteilung ist nach der Messung nur in Form von Projektionswerten vorhanden. ^[2]

Nach Radon kann eine zweidimensionale Objekteigenschaft exakt beschrieben werden, wenn eine unendliche Anzahl von Linienintegralen vorliegt. Sei $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ die Verteilung der Schwächungskoeffizienten, so kann das Objekt anhand aller Linienintegrale bestimmt werden:

$//<![CDATA[ \begin{array}{l}\int_{-\infty}^{+\infty} \! f(x, y) \, \mathrm{d}l\end{array} //]]>$   (4)

Eine Integration über alle Punkte und Richtungen kann somit die ganze Fläche des Objektes darstellen. Dabei sind einige Linienintegrale jedoch identisch, was zu Problemen führt. Um diese Herausforderung zu bewältigen wurde ein neues Ordnungsschema gewählt, bei dem alle Linienintegrale nur einmal vorkommen: Die sogenannte Hesse-Normalform. Als Ausgangspunkt der Problemlösung gilt die Annahme, dass die Projektionen von $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ unter beliebigen Winkeln bekannt sind. Unter dem Begriff Projektionen wird ein Intergral über $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ entlang einer Schar von Projektionsgeraden wahrgenommen. Diese Projektionsgeraden werden durch den Winkel $//<![CDATA[ \begin{array}{l}θ\end{array} //]]>$ zwischen der Gerade und der x-Achse, sowie dem Abstand $//<![CDATA[ \begin{array}{l}s\end{array} //]]>$ der Gerade vom Koordinatenursprung beschrieben (Abb.2). ^[4]

Wenn eine Projektionsgerade durch die Gleichung $//<![CDATA[ \begin{array}{l}g = x \cos \theta + y \sin θ\end{array} //]]>$ (5) beschrieben werden kann, so kann eine Projektion unter einem bestimmten Winkel ($//<![CDATA[ \begin{array}{l}θ = const.\end{array} //]]>$) durch folgender Beziehung beschrieben werden:

$//<![CDATA[ \begin{array}{l}\int_g \! f(x, y) \, \mathrm{d}l = p_\theta(s)\end{array} //]]>$   (6)

Die Gesamtheit aller Projektionen unter allen Projektionswinkeln $//<![CDATA[ \begin{array}{l}\theta\end{array} //]]>$ ($//<![CDATA[ \begin{array}{l}0° ≤ \theta ≤ 180°\end{array} //]]>$ und $//<![CDATA[ \begin{array}{l}s_min ≤ s ≤ s_max\end{array} //]]>$) wird als Funktion von $//<![CDATA[ \begin{array}{l}\theta\end{array} //]]>$ und $//<![CDATA[ \begin{array}{l}s\end{array} //]]>$ betrachtet $//<![CDATA[ \begin{array}{l}(p (\theta, s))\end{array} //]]>$. Die Transformation von $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ nach $//<![CDATA[ \begin{array}{l}p (\theta, s)\end{array} //]]>$ erfolgt nach Radon. Radon hat das erste mathematische Verfahren zur Rekonstruktion von Projektionen als Radon-Transformation von $//<![CDATA[ \begin{array}{l}f\end{array} //]]>$ beschrieben. Deshalb wird diese Funktion als Radon-Transformierte der Funktion $//<![CDATA[ \begin{array}{l}f(x, y)\end{array} //]]>$ bezeichnet (Abb.3). Die Beziehung lautet:

$//<![CDATA[ \begin{array}{l}Rf(\theta, s) = \int_g \! f(x, y) \, \mathrm{d}l = p(\theta, s)\end{array} //]]>$   (7)

Von den gemessenen Projektionen $//<![CDATA[ \begin{array}{l}p (\theta, s)\end{array} //]]>$ kann jedoch die Funktion $//<![CDATA[ \begin{array}{l}f(x, y)\end{array} //]]>$ nicht direkt abgeleitet werden. Das Rekonstruktionsproblem besteht darin, wie die inverse Transformation $//<![CDATA[ \begin{array}{l}R^{-1}\end{array} //]]>$ zu bestimmen ist. ^{[5] [6]}

Mathematische Verfahren zur Bildrekonstruktion

Zur Lösung des Problems der Bildrekonstruktion gibt es unterschiedliche Verfahren. Eine Möglichkeit, die zahlreichen Methoden zu klassifizieren, ist deren Aufteilung in

• analytische Verfahren,
• algebraische Verfahren sowie
• statistische Verfahren.

In der folgenden Abbildung werden die Verfahren mit zusätzlicher Angabe der wichtigsten Algorithmen zusammengefasst: ^[2]

Gegenüberstellung der Methoden

Der wesentliche Unterschied zwischen analytischen und iterativen bzw. algebraischen Verfahren liegt in der Art der Verarbeitung der Messdaten:

Die analytische Bildrekonstruktion fasst die Funktion $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ als kontinuierliche mathematische Funktion auf und somit werden auch die Projektionsdaten $//<![CDATA[ \begin{array}{l}p (\theta, s)\end{array} //]]>$ als eine kontinuierliche Funktion angenommen. Das Rekonstruktionsproblem besteht im Wesentlichen darin, eine Lösung für die Integralgleichung zu finden. Eine Diskretisierung der Messdaten erfolgt erst nach der Lösung des Integrals und wird bei der Implementierung der Inversionsformel berücksichtigt. Analytische Rekonstruktionsverfahren werden heutzutage wegen ihrer Robustheit und ihrer hoher Leistung bevorzugt. Voraussetzung für die Anwendung von algebraischen Methoden ist die Diskretisierung der Messdaten. Ziel ist dabei das Gleichungssystem $//<![CDATA[ \begin{array}{l}p = M * \mu\end{array} //]]>$  (8) zu lösen. Dafür werden aber iterative Verfahren verwendet, welche langsamer konvergieren.

Im Rahmen dieses Artikels wurden die Projektionen mithilfe der Radon-Transformation beschrieben. Im Folgenden wird in die Methoden eingegangen, mit denen das Bild durch eine inverse Radontransformation rekonstruiert werden kann. Dabei gibt es zwei Möglichkeiten die Radon-Transformation zu invertieren. Die eine Möglichkeit stellt das Fourier-Scheiben-Theorem und die andere die gefilterte Rückprojektion dar.

Fourier-Scheiben-Theorem

Die Radon-Transformierte $//<![CDATA[ \begin{array}{l}p (\theta, s)\end{array} //]]>$ der Funktion $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ reicht alleine nicht zur Lösung des gegebenen Problems. Aus diesem Grund wird die Fourier-Transformation, bzw. das Fourier-Scheiben-Theorem zur Lösung des Problems herangezogen. Diese stellt den einfachsten Lösungsansatz dar.

Folgende Zusammenhänge liegen dabei vor:

• Sei eine Funktion $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ sowie deren 2D-Fouriertransformierte $//<![CDATA[ \begin{array}{l}F (u, v)\end{array} //]]>$:

$//<![CDATA[ \begin{array}{l}f (x, y) \underrightarrow{2D - FT} F(u, v)\end{array} //]]>$   (9)

• Sei $//<![CDATA[ \begin{array}{l}p_\theta (s)\end{array} //]]>$ eine Projektion von $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ zu einem beliebigen Winkel $//<![CDATA[ \begin{array}{l}\theta\end{array} //]]>$ und $//<![CDATA[ \begin{array}{l}P_\theta (w)\end{array} //]]>$ deren 1D-Transformierte:

$//<![CDATA[ \begin{array}{l}p_\theta (s) \underrightarrow{1D - FT} P_\theta (w)\end{array} //]]>$   (10)

Die 1D-Transformierte $//<![CDATA[ \begin{array}{l}P_\theta (w)\end{array} //]]>$ beschreibt die Werte von $//<![CDATA[ \begin{array}{l}F (u, v)\end{array} //]]>$ auf einem Radialstrahl, der mit einem Winkel $//<![CDATA[ \begin{array}{l}\theta\end{array} //]]>$ zur x-Achse durch den Ursprung geht. Um die Funktion $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ zu bestimmen, wird aus allen Projektionen deren 1D-Fouriertransformierten $//<![CDATA[ \begin{array}{l}P_\theta (w)\end{array} //]]>$ gebildet. Diese Werte sind auf dem zu Winkel $//<![CDATA[ \begin{array}{l}\theta\end{array} //]]>$ gehörenden Radialstrahl in die Funktion $//<![CDATA[ \begin{array}{l}F (u, v)\end{array} //]]>$ einzutragen. Zu guter Letzt ist die inverse 2D-Fouriertransformation von $//<![CDATA[ \begin{array}{l}F (u, v)\end{array} //]]>$ erforderlich um an die Funktion $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ zu gelangen. Laut dem Fourier-Scheiben-Theorem kann also durch die inverse Fourier-Transformation die ursprüngliche Funktion $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ aus den Projektionen $//<![CDATA[ \begin{array}{l}p (\theta, s)\end{array} //]]>$ rekonstruieren werden (Abb. 4).

Das einzige Problem liegt bei der Umstellung von polaren auf kartesischen Koordinaten. Die polaren Koordinaten werden in erster Linie in die Transformation eingeführt um die Berechnung zu vereinfachen. Das Problem erscheint bei der Umstellung von polaren auf kartesischen Koordinaten, da sie nicht genau übereinstimmen können. Im polaren Koordinatensystem liegen die Werte im Bereich des Ursprunges näher aneinander und verbreiten sich im Raum mit einer zunehmenden Entfernung vom Ursprung. Dies führt zu einer auftretenden Unschärfe, welche bei der Rückprojektion ohne Filterung zu sehen ist. Unter idealen Bedingungen kann mithilfe des Fourier-Scheiben-Theorems ein perfektes Bild erzeugt werden, jedoch ist diese Methode im Fall einer nicht sorgfältigen Durchführung als rechenintensiver und fehleranfälliger zu betrachten. ^[1]

Gefilterte Rückprojektion

Die gängigste Methode zur Bildrekonstruktion ist die gefilterte Rückprojektion. Dadurch soll die aufgrund der Fourier-Transformation entstehende Unschärfe, die im vorherigen Kapitel erklärt wurde, verschwinden. Im Folgenden wird die Herleitung der grundlegenden Formel dargestellt. Die gesuchte Funktion $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ kann durch die inverse Fourier-Transformation beschrieben werden:

$//<![CDATA[ \begin{array}{l}f(x, y) = \iint\limits_{-\infty}^{+\infty}\ F(u, v) e^{i2\pi(ux + vy)} \,du\,dv\end{array} //]]>$   (11)

Im Fourier-Raum werden ebene Polarkoordinaten zur Vereinfachung der Berechnung eingeführt. Der Übergang von den kartesischen Koordinaten  $//<![CDATA[ \begin{array}{l}u\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}v\end{array} //]]>$ zu den Polarkoordinaten $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$, $//<![CDATA[ \begin{array}{l}\theta\end{array} //]]>$ wird durch $//<![CDATA[ \begin{array}{l}u = w * cos\,{\theta}\end{array} //]]>$   (12), $//<![CDATA[ \begin{array}{l}v = w * sin\,{\theta}\end{array} //]]>$   (13) und $//<![CDATA[ \begin{array}{l}\,du\,dv = w * dw * d\theta\end{array} //]]>$   (14) gegeben. Das Substituieren der Variablen sowie der Integrationsgrenzen in Gleichung 11 liefert:

$//<![CDATA[ \begin{array}{l}f(x, y) = \iint\limits_{0}^{2 \pi \infty}\ F(\theta, w) e^{i\,2\,\pi\,w(x\,cos\,\theta + y\,sin\,\theta)} \,w\,dw\,w\,d\theta\end{array} //]]>$   (15)

$//<![CDATA[ \begin{array}{l}f(x, y) = \iint\limits_{0-\infty}^{\pi + \infty}\ F(\theta, w) e^{i\,2\,\pi\,w(x\,cos\,\theta + y\,sin\,\theta)} \,|w|\,dw\,w\,d\theta\end{array} //]]>$   (16)

Durch die letzte Umformung wurden die Integrationsgrenzen geändert, die Integrationsfläche bleibt jedoch identisch. Dabei muss der Betrag von $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$ eingesetzt werden, da jetzt über negative Werte von $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$ integriert wird. Die Formel kann nochmal umformuliert werden, indem die Abkürzung $//<![CDATA[ \begin{array}{l}s = x * cos\,{\theta} + y * sin\,{\theta}\end{array} //]]>$ eingeführt wird:

$//<![CDATA[ \begin{array}{l}f(x, y) = \iint\limits_{0 - \infty}^{\pi+\infty}\ F(\theta, w) e^{i\,2\,\pi\,w\,s} \,|w|\,dw\,w\,dθ\end{array} //]]>$   (17)

Die Fourier-Transformierte $//<![CDATA[ \begin{array}{l}F (\theta, w)\end{array} //]]>$ im Integral entspricht nach dem Fourier-Scheiben-Theorem der 1D-Fourier-Transformierten $//<![CDATA[ \begin{array}{l}P_\theta (w)\end{array} //]]>$ der gemessenen Projektion. Die Gleichung kann somit auch wie folgt formuliert werden:

$//<![CDATA[ \begin{array}{l}f(x, y) = \iint\limits_{0-\infty}^{\pi + \infty}\ P_{\theta}(w)\,e^{i\,2\,\pi\,w\,s} \,|w|\,dw\,w\,d\theta = \int\limits_{0}^{\pi}\ ( \int\limits_{-\infty}^{+\infty}\ P_{\theta}(w)\,e^{i\,2\,\pi\,w\,s} \,|w|\,dw)\,d\theta\end{array} //]]>$    (18)

Um das Verfahren der gefilterten Rückprojektion besser erläutern zu können wird die Gleichung 18 in zwei Teile separat betrachtet. Das innere Integral beinhaltet die Multiplikation der Fourier-Transformierten $//<![CDATA[ \begin{array}{l}P_\theta (w)\end{array} //]]>$ mit dem Faktor $//<![CDATA[ \begin{array}{l}|w|\end{array} //]]>$ vor der Rücktransformation, was zu einer Faltung im Fourier-Raum und eine Filterung im Ortsraum resultiert. Der innere Integral wird als $//<![CDATA[ \begin{array}{l}(p_\theta) ̃(s)\end{array} //]]>$ bezeichnet:

$//<![CDATA[ \begin{array}{l}\widetilde{p_{\theta}}(s) = \int_{-\infty}^{+\infty}P_{\theta}(w) e^{i 2 \pi w s} |w| dw\end{array} //]]>$   (19)

Das Integral $//<![CDATA[ \begin{array}{l}(p_\theta) ̃(s)\end{array} //]]>$ unterscheidet sich eigentlich von der gemessenen Projektion $//<![CDATA[ \begin{array}{l}p_\theta (s)\end{array} //]]>$ nur durch den Faktor $//<![CDATA[ \begin{array}{l}|w|\end{array} //]]>$. Ohne Berücksichtigung des Faktors $//<![CDATA[ \begin{array}{l}w\end{array} //]]>$ wurde die gemessene Projektion $//<![CDATA[ \begin{array}{l}p_\theta (s)\end{array} //]]>$ resultieren. Das äußere Integral lässt sich unter Berücksichtigung der Gleichungen 5 und 18 wie folgt darstellen:

$//<![CDATA[ \begin{array}{l}f(x, y) = \int_0^{\infty} \widetilde{p_{\theta}}(s) d\theta = \int_0^{\infty} \widetilde{p_{\theta}}(s) (x cos \theta + y cos \theta) d\theta\end{array} //]]>$   (20)

Die Betrachtung der Gleichung verdeutlicht, dass die gesuchte Verteilung $//<![CDATA[ \begin{array}{l}f (x, y)\end{array} //]]>$ aus den gefilterten Projektionen $//<![CDATA[ \begin{array}{l}(p_\theta) ̃(s)\end{array} //]]>$ sich zusammensetzt. Um den Schwächungskoeffizient $//<![CDATA[ \begin{array}{l}\mu\end{array} //]]>$ an einem beliebigen Punkt $//<![CDATA[ \begin{array}{l}(x, y)\end{array} //]]>$ zu berechnen, muss von allen gefilterten Projektionen der Wert an der Stelle $//<![CDATA[ \begin{array}{l}s = x * cos\,{\theta} + y * sin\,{\theta}\end{array} //]]>$ genommen und aufaddiert werden. Die Faltung bzw. Filterung spielt bei der Rückprojektion eine wesentliche Rolle, indem sie zur Erstellung klarer Bilder beiträgt. Der Einfluss der Filterung wird bei komplexen und inhomogenen Objekten ersichtlicher. ^{[5] [1]}

Algebraische Verfahren

Bei der algebraischen Rekonstruktion kommt ein ganz anderer Lösungsweg zum Einsatz. Ziel ist dabei das Rekonstruktionsproblem aus den gemessenen Projektionswerten mithilfe von linearen Gleichungssystemen zu lösen. Das Objekt wird in Pixeln geteilt, woraus sich die Bildmatrix zusammensetzt. Bei diesem Verfahren ist es sinnvoll, die zweidimensionale Bildmatrix in einem Vektor $//<![CDATA[ \begin{array}{l}f\end{array} //]]>$ der Dimension $//<![CDATA[ \begin{array}{l}d = (N x 1)\end{array} //]]>$ umzuwandeln: ^[1]

$//<![CDATA[ \begin{array}{l}f = \begin{pmatrix} f_1 \\ f_2 \\ \vdots \\ f_d \end{pmatrix}\end{array} //]]>$   (21)

Für jeden Strahl, der durch das Bild geht, lässt sich der entsprechende Projektionswert wie folgt berechnen:

$//<![CDATA[ \begin{array}{l}p_i = \sum_{j = 1}^d a_{ij} * f_j\end{array} //]]>$   (22)

Dabei repräsentiert $//<![CDATA[ \begin{array}{l}a_{ij}\end{array} //]]>$ den Gewichtungsfaktor. Betrachtet man einen Projektionsstrahl, der das Objekt in einer beliebigen Richtung durchläuft, dann können einem Bildpunkt bzw. Pixel $//<![CDATA[ \begin{array}{l}f\end{array} //]]>$ mehrere Gewichtungsfaktoren zugeordnet werden. Aus der Summe aller Strahlen $//<![CDATA[ \begin{array}{l}P\end{array} //]]>$ und Pixeln $//<![CDATA[ \begin{array}{l}F\end{array} //]]>$ entsteht ein lineares Gleichungssystem, das wie folgt beschrieben werden kann: ^[1]

$//<![CDATA[ \begin{array}{l}p = M * f\end{array} //]]>$   (23)

Hierbei sind $//<![CDATA[ \begin{array}{l}p\end{array} //]]>$ und $//<![CDATA[ \begin{array}{l}M\end{array} //]]>$ bekannt, während $//<![CDATA[ \begin{array}{l}f\end{array} //]]>$ die Unbekannte darstellt. Das daraus resultierende Gleichungssystem kann somit groß werden, was negativ auf die Berechnung sowie die benötigte Rechenzeiten auswirkt. Aus diesem Grund wird die Anwendung von algebraischen Verfahren in der Praxis nicht bevorzugt.

Literatur

1. Cherry, S. R.; Sorenson, J. A.; Phelps, M. E.: Physics in Nuclear Medicine. Saunders. Philadelphia, Pennsylvania, 2003. ISBN 978-1-4160-5198-5
2. Kalender, W. A.: Computertomographie. Publicis Corporate Publishing. Erlangen, 2006. ISBN 978-3895780820.
3. Schott, D.: Bildrekonstruktion am Computer als Studentenprojekt. Hochschule Wismar. Wismar, 2005.
4. Schlegel, W.; Bille, J.: Medizinische Physik 2. Springer. Berlin, Heidelberg, New York, 2002. ISBN 978-3-642-56259-4.
5. Revens, M.; Zinterhof, P.: Application of numbertheoretic methods to the numerics of computerized tomography (CT). Oldenbourg. Wien, München, 1992.
6. Natterer, F.: The mathematics of computerized tomography. John Wiley & Sons Ltd and B G Teubner. Stuttgart, 1986. ISBN 3-519-02103-X.