Iteratives Bildrekonstruktionsverfahren für die Computertomographie

Inhaltsverzeichnis

Vorwort

1. Einleitung

2. Computertomographie
2.1. Physikalische und mathematische Grundlagen
2.2. Systemkomponenten
2.3. Prinzip des Messverfahrens und Scanner-Typen
2.4. Bildrekonstruktionsverfahren
2.4.1. Direkte Rückprojektion
2.4.2. Gefilterte Rückprojektion
2.4.3. Direkte Fourier-Rekonstruktion
2.4.4. Algebraische Rekonstruktion

3. Iteratives Bildrekonstruktionsverfahren
3.1. Terminologie
3.2. Mathematische Formulierungen
3.2.1. Gewichtung der einzelnen Matrixelemente
3.2.1.1. Ohne Berücksichtigung des Strahlprofils
3.2.1.2. Mit Berücksichtigung des Strahlprofils
3.2.2. Ermittlung der Schwächungskoeffizientenänderung (Abbildung in dieser Leseprobe nicht enthalten)
3.2.2.1. „lineare Gewichte“
3.2.2.2. „quadratische Gewichte“
3.2.2.3. Endgültige Schwächungskoeffizientenänderung
3.3. Berechnungen
3.3.1. Vorgegebene Werte
3.3.2. Zu berechnende Zwischenwerte
3.3.2.1. Matrixelementgröße (am)
3.3.2.2. Einzelstrahlabstand Dh
3.3.2.3. Bestimmung von (Dxs)
3.3.2.4. Beteiligte Matrixelemente
3.3.2.5. Weglängen li,j,k(z´,s´) durch die einzelnen Matrixelemente
3.3.3. Berechnung der Formeln aus Kapitel 3.2.

4. Beschreibung des Programmes ibrek
4.1.Beschreibung der Funktionen
4.1.1. Hauptprogramm (main)
4.1.2. Einlesen des Headers (lese_header)
4.1.3. Einlesen der Daten (lese_daten)
4.1.4. Sortieren der Daten (sort_struktur)
4.1.5. Einlesen des Speichernamens (speichername)
4.1.6. Einlesen der Matrixgröße (m_groesse)
4.1.7. Berechnung der Pixelgröße (pixel_groesse)
4.1.8. Festlegung der Strahlanzahl (dim_strahl)
4.1.9. Einlesen der Strahlart (strahl_art)
4.1.9.1. Rechteckprofil
4.1.9.2. Parabelprofil
4.1.9.3. Realprofil
4.1.9.4. Flächennormierung
4.1.10. Art der Schwächungswertänderung (aend_auswahl)
4.1.11. Initialisierung der „Felder“ (ini_feld)
4.1.12. Initialisierung der Schwächungswertmatrix (ini_matrix)
4.1.13. Gesamtanzahl der Iterationen / „update“ - Kriterium (ini_it)
4.1.14. Berücksichtigung von Winkel ³ 180° (tausch)
4.1.15. Strahlabstand vom Koordinatenursprung (x=0) (abstand)
4.1.16. Anzahl Elemente pro Zeile (anzahl_elemente)
4.1.17. Strahlwege durch die einzelnen Matrixelemente (gewmael)
4.1.17.1. Weglänge für a ¹ 0° und a ¹ 90°
4.1.17.2. Weglänge für a = 0° oder a = 90°
4.1.18. Berechnung der Gewichte der einzelnen Matrixelemente (gewicht)
4.1.19. Schwächungswertänderung nach Gl. (3.16) (lineare_aenderung)
4.1.20. Schwächungswertänderung nach Gl. (3.17) (quadrat_aenderung)
4.1.21. Aktualisierung der Schwächungswertmatrix (update)
4.1.22. Darstellen der rekonstruierten Matrix (zeige)
4.2. Installation des Bildrekonstruktionsprogrammes ibrek
4.2.1. ibrek
4.2.2. rek_bild.dat
4.2.3. rek_bild_sav
4.2.4. rek_bild.name
4.2.5. IDL
4.2.6. my_idl_setup
4.2.7. my_idl_startup
4.2.8. man_bild.pro
4.2.9. rek_bild_idl
4.3. Bedienung des Bildrekonstruktionsprogrammes ibrek
4.3.1. Beschreibung der zu übergebenden Parameter

5. Erste Ergebnisse
5.1. Verwendete Prüfkörper
5.1.1. „Frau“
5.1.2. „Stern“
5.1.3. „FH-Heilbronn“
5.2. Variation der Rekonstruktionsparameter
5.2.1. „Update“ Kriterium bei einem Iterationsschritt
5.2.2. Anzahl der Iterationen unter Berücksichtigung der Auswertemethode (Dµ)
5.2.3. „Update“ Kriterium bei fünf Iterationsschritten
5.2.4. Strahlprofil
5.2.5. Strahlbreite
5.3. Variation der Aufnahmeparameter
5.3.1. Anzahl der Winkelstellungen

6. Zusammenfassung

7. Ausblick

Literaturverzeichnis

Anhang

Messfileformat

Vorwort

Ich danke all denen, die mich während meines Studiums unterstützt haben und die dazu beigetragen haben, dass ich mein Studium mit Erfolg abschließen konnte

Zu besonderem Dank bin ich meinem betreuenden Dozenten Herrn Prof. Dr. K. Rauschnabel verpflichtet. Sowohl seine hervorragende Betreuung, als auch seine innovativen Impulse haben maßgeblich zum Gelingen dieser Arbeit beigetragen.

Danken möchte ich auch der Systemadministratorin, Frau Dipl. Phys. A. Szasz, die mir bei auftretenden Netzwerkproblemen stets hilfreich zur Seite stand.

Mein Dank gilt auch Herrn Dipl. Ing. (FH) Dieter Scheinecker, der mir bei der Visualisierung der rekonstruierten Bilder behilflich war.

Hiermit erkläre ich, dass die vorliegende Arbeit selbstständig von mir erstellt wurde. Sämtliche verwendeten Hilfsmittel sind im Literaturverzeichnis aufgeführt.

Die vorliegende Arbeit darf im Rahmen der Forschung und Lehre weiterverwendet werden.

Eine kommerzielle Nutzung bedarf einer gesonderten Erlaubnis des Autors.

Heilbronn, im März 1995

( Arno Holzwarth )

1. Einleitung

Im Jahre 1967 waren alle grundlegenden Voraussetzungen für die Erfindung der Computertomographie (CT) gegeben: Rekonstruktionsmathematik, Computertechnologie und Röntgendetektoren. Zu dieser Zeit, als Godfrey Newbold Hounsfield begann, sich mit der Computertomographie zu beschäftigen, war noch nicht absehbar, welche zentrale Stellung dieses Verfahren schon bald in der gesamten Technik einnehmen sollte.

In den letzten Jahren hat sich die Computertomographie in sehr vielen Bereichen als unverzichtbares Hilfsmittel etabliert. Nicht nur bei medizinischen Anwendungen, wie z.B. der Früherkennung von Tumoren, sondern auch bei der zerstörungsfreien Materialprüfung oder auch bei Verschleißmessungen findet die CT ihre Anwendung.

Um dieser Entwicklung folgen zu können, und um an der Weiterentwicklung von zur Bildrekonstruktion benötigten Verfahren mitzuarbeiten, wurde in einem Gemeinschaftsprojekt der Fachhochschule Heilbronn, der Fachhochschule Gießen und der Phillips-Universität Marburg ein voll funktionsfähiger Translations-Rotations-Scanner gebaut.

Sowohl der Computertomograph als auch die dazu benötigte Bildrekonstruktionssoftware wurde im Rahmen von mehreren Diplomarbeiten entwickelt.

Die CT ist ein spezielles Röntgen-Schichtaufnahmeverfahren, mit dem sich ein zweidimensionales Bild einer Objektschicht berechnen lässt. In der klassischen Röntgentechnik wird der vom untersuchten Objekt erzeugte Schatten auf photographischem Film fixiert, oder auf einem Bildschirm dargestellt. Dabei wird immer ein dreidimensionaler Körper auf eine zweidimensionale Ebene abgebildet. Verfahrensbedingt überlagern sich die Abbildungen aller Objektstrukturen, die im Strahlengang hintereinander liegen. Das Bild entsteht durch die Überlagerung von Informationen aus allen Objektschichten und wird deshalb auch Superpositionsbild genannt.

Bei der CT tritt dieser Überlagerungseffekt prinzipiell nicht auf, da nur eine einzige dünne Objektschicht ausgewählt und untersucht wird und nur Informationen dieser einen Schicht zur Entstehung des Bildes beitragen. Störende Einflüsse anderer Schichten werden vermieden. Ein solches von Überlagerungen freies Bild nennt man Substitutionsbild oder Tomogramm.

Auf Grund der beschränkten Leistungsfähigkeit der zur Bildrekonstruktion verwendeten Personalcomputer konnten zunächst nur die weniger rechenintensiven Rekonstruktionsmethoden, wie zum Beispiel die gefaltete Rückprojektion verwirklicht werden.

Durch die Verfügbarkeit von HP-UX Workstations mit sehr großen Rechenkapazitäten sind in der FH-Heilbronn die Voraussetzung gegeben, auch sehr rechenintensive Problemstellungen zu bearbeiten.

Im Rahmen der vorliegenden Diplomarbeit soll ein iteratives Bildrekonstruktionsverfahren für die Computertomographie in der Programmiersprache C entwickelt werden. Abweichend von den bisher vorhandenen iterativen Verfahren ist bei der angestrebten Lösung sowohl die Strahlbreite, als auch der Intensitätsverlauf innerhalb des verwendeten monoenergetischen Röntgenstrahls zu berücksichtigen. Die Implementierung des Programmes ist auf einer HP-UX Workstation unter Zuhilfenahme der Programmentwicklungsumgebung SoftBench zu bewerkstelligen.

2. Computertomographie

Ziel der Computertomographie ist die zerstörungsfreie Prüfung von Materialien, bzw. die Untersuchung von Weichteilen im menschlichen Körper ohne chirurgischen Eingriff.

Zu diesem Zweck wird das Untersuchungsobjekt mit einem senkrecht zu der untersuchenden Schicht stehenden Röntgenstrahl durchstrahlt. Anhand der auf der anderen Seite des Objektes austretenden abgeschwächten Röntgenstrahlung werden die sogenannten Absorptionskoeffizienten (µ) der durchstrahlten Schicht ermittelt.

2.1. Physikalische und mathematische Grundlagen

Die Röntgenstrahlen, die durch das Untersuchungsobjekt hindurchgehen, werden beim Durchdringen des Körpers entsprechend der Dichte des Gewebes abgeschwächt.

Die Eigenschaften von Gewebe Röntgenstrahlung zu schwächen beruht auf unterschiedlichen Wechselwirkungsprozessen beim Durchgang von Röntgenstrahlen durch Materie, wie der photoelektrischen Absorption und der Compton-Streuung. Jeder dieser Prozesse hat eine bestimmte Wahrscheinlichkeit, die außerdem eine Funktion der Strahlenergie ist. Die in einer Röntgenröhre erzeugte polychromatische Röntgenstrahlung enthält ein Spektrum von unterschiedlich energiereichen Quanten.

Für monoenergetische Strahlung wird die Röntgenstrahlungsintensität I0 beim Durchgang durch eine Schicht von Material der Dicke d auf eine Intensität I geschwächt Dies wird durch folgende Beziehung beschrieben:

Abbildung in dieser Leseprobe nicht enthalten (2.1)

Der lineare Schwächungskoeffizient hängt auch von der Strahlungsenergie ab. Für niedrige Energien ist er groß und verringert sich bei höheren Energien.

In einem Röntgenspektrum, d.h. bei polychromatischer Strahlung werden beim Durchgang durch Materie die niedrigen Energien schneller ausgefiltert als die hohen Energien. Die Folge ist, dass sich das effektive des Gewebes verringert, wenn der Röntgenstrahl das zu untersuchende Objekt durchsetzt.

Dieser Effekt wird Strahlaufhärtung [5] genannt. Um zu vermeiden, dass im CT-Bild eine Inhomogenität durch diesen Effekt entsteht, muss die Strahlaufhärtung sorgfältig auskalibriert werden.

Die Strahlungsquelle, welche in dem Computertomographen der FH-Heilbronn verwendet wird, ist keine Röntgenröhre die polychromatische Strahlung emittiert, sondern das radioaktives Präparat 241Am, welches eine monochromatische Strahlung aussendet.

Da bei diesem Gerät also nur monochromatische Strahlung auftritt, muss der Effekt der Strahlaufhärtung im weiteren Verlauf nicht weiter berücksichtigt werden.

Durchdringt nun ein Röntgenstrahl ein Objekt, so kann man sich dieses Objekt als eine Aneinanderreihung vieler würfelförmiger Körper mit unterschiedlichen Schwächungskoeffizienten µi vorstellen.

Abbildung in dieser Leseprobe nicht enthalten

Unter der Voraussetzung, dass alle Würfel dieselbe Dicke d besitzen, ergibt sich ausgehend von Gleichung (2.1) folgende Beziehung.

Abbildung in dieser Leseprobe nicht enthalten

Im Allgemeinen liegt jedoch kein Medium vor, das sich aus würfelförmigen Objekten zusammensetzt. Auch die Vereinfachung, dass innerhalb der Quader die Schwächungskoeffizienten konstant sind, kann nicht als realistisch betrachtet werden. Werden diese Vereinfachungen nicht gemacht, so ergibt sich aus Gleichung (2.2) der allgemeine Fall für die Abschwächung eines Röntgenstrahls beim Durchgang durch Materie.

Abbildung in dieser Leseprobe nicht enthalten

Durchdringt ein paralleles Strahlenbündel einen inhomogenen Körper, so ist der kontinuierlich verteilte Absorptionskoeffizient µ eine Funktion der Ortskoordinaten: µ = µ(x,y)

Aus Gleichung (2.3) folgt:

Abbildung in dieser Leseprobe nicht enthalten

Die in (2.4) ermittelte Integralgleichung läßt sich nicht ohne weiteres nach der gesuchten Verteilung µ(x,y) auflösen. Zunächst ist es nur möglich, die Gesamtabsorption zu ermitteln:

Abbildung in dieser Leseprobe nicht enthalten

Aufgabe der Rekonstruktionssoftware ist es, mittels aufwendiger mathematischer Methoden aus dieser Gesamtabsorption ein detailliertes Bild der Verteilung µ(x,y) der Schwächungskoeffizienten zu berechnen.

2.2. Systemkomponenten

Bei der Computertomographie handelt es sich um ein - im Vergleich zur photographischen und Röntgenaufnahme - sehr komplexes Verfahren. Zum Erhalt der Abbildung müssen mehrere Phasen eines Gesamtsystems durchlaufen werden, die sich verschiedener Techniken bedienen. Abb. 2.1 zeigt schematisch die zur Durchführung erforderlichen Verfahrensschritte und Systemkomponenten.

Abbildung in dieser Leseprobe nicht enthalten

Abb. 2.1: Blockschaltbild der Systemkomponenten in der Computertomographie

Zunächst werden mit Hilfe eines Messsystems, dem eigentlichen Tomographen, an der zu untersuchenden Objektschicht Absorptionsmessungen durchgeführt. Die erhaltenen Profile werden aufbereitet und im Computer gespeichert. Mit Hilfe einer entsprechenden Software wird aus den Profilen eine zweidimensionale bildliche Darstellung der gesuchten Verteilung der Absorptionskoeffizienten Abbildung in dieser Leseprobe nicht enthalten ermittelt. Auf diese Weise gewinnt man den Schwächungskoeffizienten für jeden Ort innerhalb der Objektschicht. Das Ergebnis dieser Rechnungen ist zunächst eine Koeffizientenmatrix, bei der jeder Wert den Absorptionskoeffizienten des Materials an der betreffenden Stelle der Schicht angibt. Die Position des Wertes in der Matrix entspricht der Lage des Bildpunktes in der untersuchten Schicht.

Zur Darstellung und Beurteilung kann die Matrix direkt in dieser Form von Zahlen auf einem Drucker ausgegeben werden. Sie kann aber auch weiterverarbeitet und durch andere Peripheriegeräte in geeigneter Weise visualisiert werden. Dies kann zum Beispiel durch Darstellung der Koeffizienten als Grauwerte auf einem Monitor geschehen. Es ist möglich, die Bildmatrix zusätzlich zu manipulieren und durch eine Farbdarstellung die Beurteilungsmöglichkeit zu verbessern.

Man kann das Gesamtverfahren der Computertomographie im Wesentlichen in vier Schritte gliedern:

1) Das Messsystem zur Gewinnung der Profile:

Neben der Lösung mechanischer Probleme ist dies hauptsächlich eine physikalisch- technische Aufgabe zur effektiven Messung ionisierender Strahlung.

2) Das Rekonstruktionsverfahren zur Ermittlung der zweidimensionalen Verteilung der Absorptionskoeffizienten aus den eindimensionalen Profilen:

Hierbei handelt es sich im Wesentlichen um ein mathematisches Problem.

3) Die Präsentation der Ergebnisse:

Ziel ist eine für den jeweiligen Zweck geeignete und für den Anwender anschauliche Darstellung des Schichtaufbaus.

4) Die Bildanalyse:

Bestimmung von Objektstrukturgrößen oder Ermittlung von Dichteverteilungen innerhalb des Objektes.

2.3. Prinzip des Messverfahrens und Scanner-Typen

Bei der Computertomographie werden zunächst die Absorptionsprofile aufgenommen. Hierzu wird das zu untersuchende Objekt mit einem Röntgenstrahl durchstrahlt. Auf der gegenüberliegenden Seite der Strahlungsquelle befindet sich ein Detektorsystem, das die aus dem Objekt austretende, abgeschwächte Strahlung registriert. Die Detektoren wandeln die auftreffende Strahlung in elektrische Impuls um. Diese werden über eine A/D-Wandlerkarte in digitale Signale umgeformt und somit einer Speicherung durch den Computer zugänglich gemacht. Nachdem das Objekt aus mehreren, radialsymmetrischen Richtungen durchstrahlt wurde, wird mittels der so gewonnenen Datensätze mit einer Bildrekonstruktionssoftware das ursprüngliche Bild der untersuchten Schicht errechnet. Bei der Berechnung wird jedem Schwächungswert ein Zahlenwert zugewiesen. Durch die Darstellung dieser Zahlen als Grauwerte kann eine Visualisierung auf einem Sichtgerät (z.B. Monitor) erfolgen.

Einfacher Translations-Rotations-Scanner

Da es sich bei dem in der FH-Heilbronn verwendeten Scanner um einen einfachen Translations-Rotations-Scanner handelt, soll dieser im weiteren Verlauf etwas näher erläutert werden.

Das Prinzip der Datenerfassung lässt sich anschaulich an dem messtechnischen Aufbau gemäß Abb. 2.2 erklären:

Strahlenquelle und Detektor sind auf einem gemeinsamen Gestell mechanisch fest miteinander verbunden, so dass sich Quelle und Detektor immer genau gegenüberstehen. Bei einer bestimmten Winkelstellung, zum Beispiel bei Abbildung in dieser Leseprobe nicht enthalten wird die Schwächung des Strahls längs des gesamten Weges durch das Objekt gemessen. Durch Parallelverschiebung des Systems Strahlenquelle-Detektor in Richtung der Translations-Koordinate Abbildung in dieser Leseprobe nicht enthalten erfasst man das Absorptionsprofil für diese Winkelstellung. Ein Absorptionsprofil, auch Projektion genannt, wird demnach durch lineare translatorische Bewegungen bei einem bestimmten Projektionswinkel gewonnen.

Nach Beendigung der Messung einer Projektion durch mehrere Translationsbewegungen wird das Messsystem um einen bestimmten Winkel gedreht, es erfolgt die Aufnahme einer weiteren Projektion. Auf diese Weise erhält man eine Vielzahl von Projektionen. Wird zum Beispiel mit einem Winkelintervall Abbildung in dieser Leseprobe nicht enthalten gemessen, so erhält man insgesamt 180/2 = 90 Projektionen.

Bei Verwendung von monochromatischer Strahlung genügt die Beschränkung auf einen gesamten Winkelbereich von 180°. Das Nichtauftreten der Strahlaufhärtung [8] bei monochromatischer Röntgenstrahlung ist der Grund dafür, dass eine Projektion bei einem Winkel Abbildung in dieser Leseprobe nicht enthalten gegenüber der Projektion bei einem Winkel Abbildung in dieser Leseprobe nicht enthalten keine weitere Information enthält.

Vor dem Detektor befindet sich der Kollimator (eine Blende mit dünnem Spalt), der die Aufgabe hat, aus dem räumlich ausgedehnten Strahlenbündel einen kleinen Bereich auszublenden, der schließlich zum Detektor gelangt. Der Kollimator ist für ein gutes räumliches Auflösungsvermögen der Apparatur notwendig. Die Spaltbreite muss so groß sein, dass die Anzahl der am Detektor ankommenden Quanten noch ausreicht, um den statistischen Fehler der Messung ausreichend klein zu halten.

Abbildung in dieser Leseprobe nicht enthalten

Abb. 2.2: Einfacher Translations-Rotations-Scanner mit einem Detektor[1]

Das Prinzip des einfachen Translations-Rotations-Scanners mit einem Detektor wurde bei den ersten Tomographen eingesetzt, hat aber neben dem Vorteil der geringen Streustrahlung zwei schwerwiegende Nachteile:

- Die von der Strahlenquelle erzeugte Gesamtstrahlung wird nur zu einem geringen Grad ausgenutzt, da durch die erforderliche scharfe Ausblendung eines kleinen Strahlenbündels auch nur ein kleiner Teil der Strahlung zur Informationsgewinnung genutzt wird.
- Es dauert sehr lange, bis ein vollständiger Satz von Projektionen vorliegt, da nicht gleichzeitig Translations- und Rotationbewegungen möglich sind. Im günstigsten Fall liegt die Dauer der Messung im Minutenbereich. In der medizinischen Anwendung bewirken lange Scanzeiten eine große Strahlenbelastung des Patienten und eine Erhöhung der Gefahr von Bildverfälschungen durch Patientenbewegungen (Bewegungsartefakte).

Aufgrund der schwerwiegenden Nachteile dieses Systems wurden seit der Einführung der Computertomographie erhebliche Anstrengungen unternommen, die Geräte insbesondere durch Verkürzung der Aufnahmezeiten leistungsfähiger zu machen. Dies führte zur Entwicklung verschiedener Abtastsysteme, die sich im Wesentlichen durch Anzahl und Anordnung der Detektoren unterscheiden. Die nachfolgende Gegenüberstellung zeigt die Vor- und Nachteile der einzelnen Systeme.

Mehrdetektor-Translations-Rotations-Scanner

Das Funktionsprinzip ist gleich dem des Translations-Rotations-Scanners. Die Rotationsabstände D j sind allerdings größer, da mehrere Detektoren verwendet werden, die bei einer Translationsbewegung gleichzeitig die Projektionen verschiedener Winkelstellungen aufnehmen. Dies wird durch einen fächerförmigen Strahlenverlauf der Quelle und einer nicht parallelen Anordnung der Detektoren ermöglicht, so dass bei einer Winkelstellung des Systems mehrere Messrichtungen vorliegen.

Vorteile dieses Prinzips:

- Es wird eine erheblich kürzere Gesamtmessdauer benötigt, da mehrere Messrichtungen bzw. Projektionen pro Winkelstellung eine kleinere Anzahl von Winkelstellungen erfordert. Bei Verwendung von n Detektoren ist nur der n-te Teil der Rotationsstellungen notwendig. Es werden Scanzeiten bis zu 20 Sekunden erreicht.
- Durch kürzere Scanzeiten ergeben sich in der medizinischen Anwendung geringere Strahlenbelastungen der Patienten.
- Die Strahlenquelle wird besser ausgenutzt, da mehrere Bereiche des Strahlenbündels zur Informationsgewinnung verwendet werden.
- Unproblematische technische Realisierung aufgrund des einfachen Prinzips.

Fächerstrahl-Verfahren (Fan-Beam-Verfahren)

Es wird ein Rotationsscanner verwendet, bei dem die Strahlenquelle fest mit den Detektoren verbunden ist. Die Ausnutzung der Strahlenquelle ist optimal, weil der Strahlenfächer den gesamten Objektquerschnitt erfasst und ein vielzelliges Detektorsystem (bis zu 1024 Detektoren pro Schicht [5]) viele Messpunkte simultan aufnehmen kann. Die translatorische Bewegung entfällt, die Rotation des Systems kann kontinuierlich verlaufen. Ein Besonderheit des Verfahrens ist die etwas aufwendigere - zur Weiterverarbeitung notwendige - Datenaufbereitung, da aus den gemessenen Zentralprojektionen Parallelprojektionen gewonnen werden müssen.

Ringdetektor-Gerät

Dieser Scannertyp verkürzt die Messzeiten und trägt aufgrund der Vermeidung von Ringartefakten infolge der Röntgenstrahlgeometrie zu einer besseren Bildqualität bei.

Das Ringdetektor-Gerät ist eine Weiterentwicklung des Fächerstrahl-Gerätes und hat stationäre Detektoren, die auf einem Vollkreis um das Untersuchungsobjekt angeordnet sind. Das einzige bewegte Element ist die Strahlungsquelle innerhalb der feststehenden Detektoren. Das System ist technisch aufwendig. Dies wirkt sich nachteilig auf den Preis und die Störungsanfälligkeit des Gerätes aus.

Die höhere Geschwindigkeit der Fächerstrahlscanner wird durch größere und kompliziertere Detektorsysteme erreicht. Der auf das Detektorsystem treffende Streustrahlenanteil ist wegen der Strahlgeometrie groß.

2.4. Bildrekonstruktionsverfahren

Die Prinzipien der Bildrekonstruktion lassen sich gut anhand der Ausführung eines Parallelstrahlgerätes, dem einfachen Translations-Rotations-Scanner (Abb. 2.2), erklären.

Beim Parallelstrahlprinzip liefert ein Strahl einen Messwert, eine Anzahl von Strahlen liefert ein Messwertprofil. Nach Logarithmierung des Quotienten aus gemessener Intensität Io in Luft ohne Objekt und der gemessenen Intensität I( j , h ) nach Durchstrahlung des Objekts erhält man eine sogenannte Projektion p:

Abbildung in dieser Leseprobe nicht enthalten (2.6)

Mit Gl. (2.5) ergibt sich daraus:

Abbildung in dieser Leseprobe nicht enthalten (2.7)

In Wirklichkeit ist µ(x,y) eine kontinuierliche zweidimensionale Funktion, die unendlich viel Information enthält. Zur ihrer Bestimmung wäre ebenfalls unendlich viel Information notwendig. Das Problem der Rekonstruktion ist es, aus einer endlichen Anzahl von Daten Rückschlüsse auf µ(x,y) zu ziehen. Dies ist nur möglich, indem die Objektschicht gedanklich diskretisiert wird (Abb. 2.3). Unter Diskretisierung versteht man die Zerlegung der Objektschicht in Bereiche (z.B. Quadrate) in denen der Schwächungskoeffizient einen konstanten Wert besitzt. Dieser Wert stellt somit den Mittelwert der tatsächlichen Verteilung innerhalb dieses Bereiches dar. Das Resultat dieser Vereinfachung ist die Reduzierung der zu ermittelnden Information auf ein endliches Maß.

Die Objektschicht baut sich gedanklich aus vielen kleinen diskreten Teilbildflächen (Pixel) auf, die aber aufgrund der endlichen Dicke der betrachteten Objektschicht (endliche Höhe des Kollimatorspalts) eigentlich Volumenelemente (Voxel) darstellen.

In Abbildung (2.3) ist das Abtastprinzip des Translations-Rotations-Scanners dargestellt.

Abbildung in dieser Leseprobe nicht enthalten

Abb. 2.3: Abtastung des diskretisierten Objekts nach dem Parallelstrahlprinzip

Zunächst erfolgt eine Messung unter dem Winkel j =0. Die schrittweise lineare Abtastbewegung wird in Richtung der momentanen Translationskoordinate h vorgenommen. Die Werte der diskreten Projektion p( j =0°, h ) ergeben sich, wenn man nach Gl. (2.7) entlang jedes Messstrahls die Produkte aus den Schwächungskoeffizienten (angegeben in relativen Einheiten 1/LE ) und den durchstrahlten Wegen (in Längeneinheiten LE) der zugehörigen Volumenelemente aufsummiert. Dies sei hier beispielsweise für den Projektionswert getan, der in der Schichtmatrix zur zweiten Zeile von unten gehört. Es ergibt sich ein Projektionswert Abbildung in dieser Leseprobe nicht enthalten.

Nach vollständiger Abtastung des Objektes, d.h. nachdem die Projektion p( j =0°, h ) aufgenommen wurde, wird das Objekt relativ zum Messsystem Quelle-Detektor gedreht und es können weitere Projektionen aufgenommen werden (beispielsweise bei Abbildung in dieser Leseprobe nicht enthalten).

Somit erhält man schließlich einen vollständigen Satz von Projektionen, dessen Größe unter anderem von der Größe des Objektes und der gewünschten Ortsauflösung abhängt.

Für die Auflösung der Projektionen nach der gesuchten Verteilung der Schwächungskoeffizienten gibt es unterschiedliche Vorgehensweisen:

- Wenn es sich um einen unendlich dünnen, monochromatischen Röntgenstrahl („pencil beam“) handelt, können die Gleichungen unmittelbar diskretisiert werden, d.h. man kann sie in ein lineares algebraisches Gleichungssystem mit bis zu einigen 100000 Unbekannten überführen, welches dann gelöst werden muss.
- Es ist möglich, die Diskretisierung zugunsten einer analytischen Behandlung erst bei der numerischen Auswertung der Lösungsformeln vorzunehmen.

Das zentrale Problem der Bildrekonstruktion ist die Umkehrung der Gleichung (2.5), um die Werte der Dichtefunktion aus den gemessenen Projektionswerten (Linienintegrale) zu erhalten. Im Jahre 1917 löste als erster Radon dieses Problem analytisch (inverse Radontransformation)

Als Vereinfachung wird die zu rekonstruierende Ebene in Quadrate, sogenannte „Pixel“ unterteilt. Innerhalb jedes dieser Pixel ist der zugehörige Abschwächungsbeitrag, also der Betrag der Dichtefunktion, als konstant anzusehen. Die Matrixelemente stellen die in der Rekonstruktion zugeordneten Abschwächungswerte eines jeden Pixels dar.

Die in der Computertomographie verwendeten Rekonstruktionsverfahren lassen sich in

- die algebraische Rekonstruktion,
- die direkte und gefilterte Rückprojektion und
- die direkte Fourier-Rekonstruktion aufgliedern.

Die beiden letzten Verfahren werden häufig unter dem Oberbegriff "Analytische Verfahren" zusammengefasst.

2.4.1. Direkte Rückprojektion

Die Methode der direkten Rückprojektion wurde von den Pionieren der Computertomographie benutzt, ergab jedoch nur unbrauchbare Bilder, weil hierbei eine verfahrensmäßig bedingte Unschärfe erzeugt wird, die an den Effekt der Verwischung aus der konventionellen Röntgentechnik erinnert.

Bei der direkten Rückprojektion werden die Profile über das gesamte Bild in der der Aufnahmerichtung der Projektion entgegengesetzten Richtung verschmiert. Hierbei liefert jeder Profilwert einen gleichen Beitrag zu allen Punkten entlang der Projektionslinie. Diese Unschärfe tritt besonders deutlich bei der Rekonstruktion eines einzelnen Punktes (Delta-Funktion) hervor. Das Profil aus jeder Richtung besteht aus einer einzigen Spitze. Bei der Rekonstruktion entsteht ein Stern mit seinem Kreuzungspunkt am Ort des ursprünglichen Punktes. Er besteht aus so vielen Strahlen wie die Anzahl der Projektionen.

2.4.2. Gefilterte Rückprojektion

Das am häufigsten verwendete Verfahren zur Bildrekonstruktion bei kommerziellen Computertomographen ist die gefilterte Rückprojektion. Aufgrund der sehr geringen Rechenzeit wird dieses Verfahren am häufigsten in der Medizin eingesetzt.

Bei der gefilterten Rückprojektion werden die Profile p( j , h ) vor der Rückprojektion gefiltert und dadurch gerade so modifiziert, dass die verfahrensbedingte, unerwünschte Faltung der direkten Rückprojektion zwar nicht vermieden, jedoch weitgehend kompensiert wird. Erst diese gefilterten Profile werden rückprojiziert, d.h. gemäß der Systemgeometrie auf die einzelnen Bildpunkte der Bildmatrix verteilt. Die Filterung ist nichts anderes, als eine Faltung der Daten mit einem geeigneten Faltungskern, welcher der ersten unerwünschten Faltung entgegenwirkt. Der Faltungsprozess entspricht einer Hochpassfilterung, die hohen Signalfrequenzen werden verstärkt.

2.4.3. Direkte Fourier-Rekonstruktion

Grundlage der direkten Fourier-Rekonstruktion ist der definierte Zusammenhang zwischen den eindimensionalen Fourier-Transformierten der gemessenen Projektionen p( j , h ) und der zweidimensionalen Transformierten der Verteilung der Schwächungskoeffizienten µ(x,y). Aufgrund dieser Beziehung ist es möglich die Schwächungskoeffizienten zu ermitteln [1].

Den genauen Zusammenhang beschreibt das sogenannte Projektionstheorem [3].

2.4.4. Algebraische Rekonstruktion

Die algebraischen Bildrekonstruktionsverfahren sind die historisch ältesten in der Computertomographie. Aufgrund des sehr hohen Rechenaufwandes fanden sie in der letzten Zeit in der Medizin kaum noch Anwendung, da ihnen die neueren Integraltransformationsverfahren in dieser Beziehung weit überlegen sind. Mit der Entwicklung von immer leistungsfähigeren Rechnern (Workstations, Parallelrechner) können heutzutage jedoch diese Verfahren wieder ernsthaft in Betracht gezogen werden, da diese schnellen Rechner relativ geringe Rekonstruktionszeiten ermöglichen.

Ein großer Vorteil der algebraischen Rekonstruktionsverfahren liegt darin, dass bei diesen Methoden im Gegensatz zu den analytischen Verfahren fehlende Datensätze (Winkelstellungen) kein Problem darstellen. Auch brauchen die Abstände der einzelnen Projektionen innerhalb eines Schwächungsprofils bei den algebraischen Methoden nicht äquidistant sein, um ein artefaktfreies Bild zu liefern. Dadurch, dass keine äquidistante Abtastung erfolgen muss, können bestimmte Regionen mit hohem Informationsgehalt durch kleinere Abtastintervalle sehr genau untersucht werden. Bei der vollständigen Absorption der Röntgenstrahlung bei zum Beispiel einer Winkelstellung, tritt bei der gefalteten Rückprojektion ein durch das ganze Bild gehendes, fächerförmiges Artefakt auf. Die algebraischen Methoden hingegen weisen diese - das ganze Bild störende - Form der Artefakte nicht auf.

Man stellt sich die Objektschicht in kleine kubische Zellen zerlegt vor (siehe Abb. 2.3). In diesem Modell entspricht einem gemessenen Projektionswert p( j i , h i ) die Summe der Produkte aus den Schwächungskoeffizienten µk aller vom Strahl getroffenen Zellen und den innerhalb der einzelnen Zellen zurückgelegten Wegen D z k,j (s.a. Gl. (2.7)).

Abbildung in dieser Leseprobe nicht enthalten (2.8)

Schreibt man diese Summengleichung für jeden Projektionswert gemäß Gl. (2.7) explizit auf, so ergibt sich folgendes Gleichungssystem für eine bestimmte Projektion, beispielsweise der Projektion unter dem Winkel j i (Abb. 2.4):

Abbildung in dieser Leseprobe nicht enthalten

Abb. 2.4: Gleichungssystem einer Projektion [1]

Verfährt man für alle anderen Projektionen auf die gleiche Weise, so ergibt sich schließlich ein Gleichungssystem mit Abbildung in dieser Leseprobe nicht enthalten Unbekannten. Die Größe des Gleichungssystems reduziert sich jedoch drastisch, da jeder Projektionsstrahl die meisten Zellen der Bildmatrix überhaupt nicht durchquert und die zugehörigen durchstrahlten Wege deshalb Null sind.

Vor der Auflösung des Gleichungssystems nach den gesuchten Schwächungskoeffizienten müssen die durchstrahlten Wege der Elemente, die von einem Projektionsstrahl getroffen werden, erst noch ermittelt werden. Abb. 2.5 zeigt die hierbei auftretende Problematik. Der Strahl trifft viele Pixel nur teilweise und hat außerdem eine endliche Ausdehnung. Der durchstrahlte Weg durch ein Pixel, beispielsweise durch den schraffierten Bereich, ist nicht eindeutig.

Bei einer anderen, häufiger verwendeten Methode (in Abb.2.5 oben), wird einfach der Mittelpunkt des Pixels auf das Profil projiziert und dann zwischen zwei abgetasteten Projektionswerten interpoliert.

Abbildung in dieser Leseprobe nicht enthalten

Abb. 2.5: Bestimmung der durchstrahlten Wege [1]

Letztendlich ergibt sich das eigentliche Gleichungssystem nach Abb. (2.6), in dem nur noch diejenigen Absorptionskoeffizienten enthalten sind, deren Zellen von den jeweiligen Projektionsstrahlen auch tatsächlich durchquert wurden. Außer den gesuchten Absorptionskoeffizienten sind im Gleichungssystem keine Unbekannten mehr enthalten. Es kann prinzipiell gelöst werden (Abb. 2.6).

Bei Aufnahme von 4 Projektionen könnte man die 6 angegebenen Gleichungen gewinnen, aus denen die 4 Unbekannten eindeutig zu bestimmen wären. Die Schwierigkeiten bei der praktischen Realisierung dieses vermeintlich einfachen, exakten algebraischen Verfahrens werden durch ein anderes Beispiel deutlich. Bei der Aufteilung einer Objektschicht in beispielsweise 256 mal 256 Elemente sind 65356 unbekannte Werte zu bestimmen und man benötigt zur eindeutigen Bestimmung aller Unbekannten 65536 linear unabhängige Gleichungen. Die Lösung eines solchen Gleichungssystems führt selbst beim Einsatz größerer Rechner zu erheblichen Schwierigkeiten:

- Inakzeptable Rechenzeit
- Große Ungenauigkeiten durch Rundungsfehler
- Unterbestimmtheit des Gleichungssystems, da die Anzahl der zu bestimmenden Größen die Anzahl der vorliegenden Messdaten übersteigt

Abbildung in dieser Leseprobe nicht enthalten

Abb. 2.6: Algebraische Bildrekonstruktion durch exakte Lösung [1]

Aufgrund der erheblichen Mängel des Verfahrens, kann die Lösung des Gleichungssystems praktisch nur iterativ approximiert werden. In Abb.2.7 wird die prinzipielle Vorgehensweise bei einer solchen iterativen Approximation bildlich dargestellt.

Abb.2.7 stellt gleichzeitig das grobe Prinzip des hier entwickelten iterativen Bildrekonstruktions-verfahrens dar.

Vereinfachung bei dem in Abb.2.7 dargestellten Beispiel:

- Die Weglängen des Röntgenstrahls innerhalb der an einer Richtung beteiligten Matrixelemente sind gleich groß.
- Der Röntgenstrahl ist unendlich dünn („pencil beam“) und trifft jeweils nur eine Matrix- elementreihe.

Aufgrund dieser Vereinfachungen kann der Projektionswert für einen Strahl durch ein einfaches aufsummieren aller Schwächungskoeffizienten erfolgen. Die Differenz zwischen dem gemessenen und dem so berechneten Projektionswert, kann aufgrund der gleichen Weglängen auch zu gleichen Teilen auf die am behandelten Strahl beteiligten Pixel verteilt werden, d.h. alle Gewichte sind in diesem Fall gleich groß (s.a. 3.2.1.1.).

Da in der Wirklichkeit der Röntgenstrahl eine endliche Ausdehnung besitzt, und darüber hinaus die Weglängen innerhalb der an dem Strahl beteiligten Matrixelemente nicht gleich groß sind, wird eine Möglichkeit gesucht, diese realen Faktoren in den mathematischen Gleichungen zu berücksichtigen. In dem nun folgenden Kapitel 3 wird erläutert, wie diese Einflüsse berücksichtigt werden müssen, um eine Berechnung der Projektionswerte (Pi,j) und der Gewichte (wi,j(z,s)) zu ermöglichen. Die berechneten Gewichte legen dann fest um welchen Wert (Abbildung in dieser Leseprobe nicht enthalten) die Matrixelemente geändert werden müssen, um eine Annäherung an die gesuchte Matrix zu erzielen.

Abbildung in dieser Leseprobe nicht enthalten

Abb. 2.7: Einfaches Beispiel einer algebraischen Bildrekonstruktionstechnik [5]

3. Iteratives Bildrekonstruktionsverfahren

Wie schon erwähnt besitzt ein realer Röntgenstrahl zum einen eine endliche Strahlbreite und zum anderen einen bestimmten Intensitätsverlauf innerhalb des Röntgenstrahls.

Im folgenden Kapitel wird erläutert, wie diese realen Einflüsse berücksichtigt werden müssen, um eine Berechnung der Projektionswerte (Pi,j) und der Gewichte (wi,j(z,s)) zu ermöglichen. Die berechneten Gewichte legen dann fest, um welchen Wert (Abbildung in dieser Leseprobe nicht enthalten) die Pixelwerte geändert werden müssen, um eine Annäherung an die zu den gemessenen Schwächungsprofilen gehörende Matrix zu erzielen.

Es gibt unterschiedliche Ansätze die Gewichtsfaktoren der am Strahl beteiligten Matrixelemente zu bestimmen.

- Drehung der Bildmatrix [7]
- Drehung des Röntgenstrahls

Bei der Drehung der Bildmatrix sind die Gewichte immer gleich 1. Die Probleme beim Drehen der Bildmatrix liegen in der Interpolation der Matrixelemente.

Das hier entwickelte Verfahren geht von einer feststehenden Matrix aus, auf die der Röntgenstrahl aus unterschiedlichen Winkelstellungen einfällt. In diesem Fall ergeben sich für die Gewichtung der einzelnen Matrixelemente individuelle Werte. Dieses Vorgehen beinhaltet ebenfalls gewisse Schwierigkeiten, die im weiteren Verlauf näher erläutert werden.

Ein weiteres Problem bei einem ausgedehnten Röntgenstrahl ist, dass ein und dasselbe Pixel von mehreren Messstrahlen (j) aus unterschiedlichen Richtungen (i) getroffen werden kann. Aus diesem Grund muss eine sinnvolle Mittelung aller Schwächungskoeffizientenänderungen (Dµi,j(z,s)) durchgeführt werden, so dass sich für jedes Matrixelement eine endgültige Schwächungskoeffizientenänderung (Abbildung in dieser Leseprobe nicht enthalten) ergibt.

3.1. Terminologie

Im folgenden Abschnitt werden die im weiteren Ablauf benutzten Bezeichnungen der Einzelgrößen definiert.

Abbildung in dieser Leseprobe nicht enthalten

[...]