Stand: 13. Oktober 2003

Hiermit versichere ich, diese Arbeit selbst¨ andig verfaßt und keine anderen als die angegebenen Quellen und Hilfsmittel benutzt zu haben.

Karlsruhe, den

Zusammenfassung

Gegenstand der vorliegenden Arbeit sind die Simulation der Computertomographie (CT) sowie die Implementierung und der Vergleich verschiedener Verfahren zur Rekonstruktion der untersuchten Objekte.

Kapitel 1 befaßt sich mit den Grundlagen der CT unter verschiedenen Randbedingungen. Es werden die ausgenutzten physikalischen Prozesse erw¨ ahnt, gebr¨ auchliche Bauarten von CT-Scannern vorgestellt und die Simulation des Scanvorganges erl¨ autert. Die Modellierung einiger in der Realit¨ at unvermeidlicher St¨ orungen der Datengewinnung ist das Thema von Kapitel 2.

Mit den mathematischen Verfahren zur Bildrekonstruktion besch¨ aftigt sich Kapitel 3. Deren tats¨ achliche Implementierungen folgen in Kapitel 4. Den Abschluß der Betrachtungen bilden ¨ Uberlegungen zur Bewertung der erreichten

Rekonstruktionsqualit¨ at und deren praktische Anwendung. Kapitel 5 definiert geeignete Qualit¨ atskriterien, Kapitel 6 wendet diese unter verschiedenen Randbedingungen auf die zur Verf¨ ugung stehenden Rekonstruktionsverfahren an.

Der Anhang schließlich liefert denjenigen Lesern, welche die implementierte Software selbst anwenden m¨ ochten, Beschreibungen der verwendeten Datenformate, der Testumgebungen sowie der Programmierschnittstellen der wesentlichen Pakete.

Inhalt

1  Simulierte Computertomographie  1

1.1  Physikalische Grundlagen  1

1.2  Zweidimensionale Scanner  2

1.2.1  F acherstrahlscanner  3

1.2.2  Parallelstrahlscanner  3

1.3  Wahl des Rekonstruktionsbereiches  5

1.4  Phantome  5

1.5  Simulation der Datengewinnung  7

2  St orquellen  10

2.1  Photonenstatistik  10

2.2  Streuung (scatter)  12

2.3  Strahlaufh artung (beam hardening)  13

2.4  Reihenfolge der St orungen  14

3  Verfahren zur Bildrekonstruktion  15

3.1  Kontinuierliche Rekonstruktionsverfahren  15

3.1.1  Gefilterte R uckprojektion  16

3.1.1.1  Differenzierungsoperator  16

3.1.1.2  Hilberttransformation  16

3.1.1.3  R uckprojektion  18

3.1.2  Faltungsverfahren f ur parallele Strahlen  18

3.1.3  Rekonstruktion durch Fouriertransformation  20

3.1.4  Rho Filtered Layergram  20

3.2  Diskrete Rekonstruktionsverfahren  20

3.2.1  Additive algebraische Rekonstruktion  21

3.2.2  Tricks  23

i

3.3 Nichtiterative Reihenentwicklungsverfahren . . . . . . . . . . . . . . . . . 23 3.3.1 Harmonische Dekomposition . . . . . . . . . . . . . . . . . . . . . 24 3.3.2 Polynomiale Dekomposition . . . . . . . . . . . . . . . . . . . . . . 25 3.3.3 Effiziente Berechnung . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.4 Hilfsoperatoren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.4.1 Parallelisierungsverfahren f¨ ur divergente Strahlen (rebinning) . . . 27 3.4.2 Sch¨ atzen der mittleren Schw¨ achung aus den Rohdaten . . . . . . . 28 3.4.3 Normalisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4 Implementierungen 30

4.1 hilbertOperator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.2 backProjection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.3 convolutionOperator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.4 additiveART . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.5 polynomialDecomp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.6 rebinning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

5 Abstandsmaße 41

5.1 Abstandskennwerte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.2 Schnittdiagramme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

6 Rekonstruktionsqualit¨ at 44

6.1 Ideale Rohdaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 6.1.1 Faltungsverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 6.1.2 Algebraische Rekonstruktion . . . . . . . . . . . . . . . . . . . . . 44 6.1.3 Polynomiale Dekomposition . . . . . . . . . . . . . . . . . . . . . . 46 6.2 Auswirkungen der Photonenstatistik . . . . . . . . . . . . . . . . . . . . . 46 6.3 Auswirkungen der klassischen Streuung . . . . . . . . . . . . . . . . . . . 50 6.4 Realistische Rohdaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 6.5 Zeitaufwand . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

A Dateiformate 55

A.1 Phantomdefinition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 A.2 Scannerbeschreibung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 A.3 Graphikformate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 A.3.1 Scannergeometrie . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 A.3.2 Skalierungsinformation . . . . . . . . . . . . . . . . . . . . . . . . . 58 A.3.3 Transformationskennungen . . . . . . . . . . . . . . . . . . . . . . 58

ii

B Kommandosyntax 60

C Modulschnittstellen 62

Kapitel 1

Simulierte Computertomographie

1.1 Physikalische Grundlagen

R¨ ontgenstrahlen werden beim Durchgang durch ein Objekt in Abh¨ angigkeit von den Eigenschaften des durchstrahlten Materials ged¨ ampft ¹ . Das Ausmaß der Schw¨ achung des Strahls wird durch den sogenannten Schw¨ achungskoeffizienten µ (engl. attenuation coefficient) bestimmt. Dieser Koeffizient hat die Einheit [cm −1 ].

F¨ ur eine anf¨ angliche Strahlungsintensit¨ at J 0 ergibt sich nach dem Durchlaufen einer homogenen Schicht der Dicke d eine Restintensit¨ at J d von J d = J 0 e −µd (Schw¨ achungsgesetz). (1.1)

Betrachtet man anstelle des Strahls die einzelnen Photonen, so k¨ onnen diese aufgrund ihrer Quanteneigenschaften nat¨ urlich nicht abgeschw¨ acht werden. Jedes einzelne kann nur entweder das Objekt durchlaufen oder absorbiert werden. Analog zu Gleichung (1.1) sind aber Aussagen m¨ oglich ¨ uber die Wahrscheinlichkeit P , mit der ein einzelnes Photon ungest¨ ort das Objekt durchl¨ auft. Diese ist gerade J d

= e −µd . P = (1.2)

J 0

F¨ ur eine Schicht der Dicke 1cm mit dem Schw¨ achungskoeffizienten µ erf¨ ullt P damit die Gleichung µ = − ln P. (1.3)

Sind N solcher Schichten hintereinander angeordnet, ist die Wahrscheinlichkeit P ^{(N )} , daß ein Photon sie alle ungest¨ ort durchquert, gerade das Produkt der einzelnen Wahrscheinlichkeiten P i f¨ ur i = 1, . . . , N − 1. Der Schw¨ achungskoeffizient aller Schichten zusammen berechnet sich daher als N −1 N −1

µ N = − ln P ^{(N )} = − ln P i = (− ln P i ) . (1.4) i=0 i=0

¹ Die physikalischen Effekte, die der CT zugrunde liegen, k¨ onnen hier nat¨ urlich nur oberfl¨ achlich behandelt werden. F¨ ur eine genauere Darstellung sei auf die einschl¨ agige Physikliteratur,

beispielsweise [Str90] S. 412ff, verwiesen.

1

2 Kapitel 1: Simulierte Computertomographie

Faßt man die Gleichungen (1.3) und (1.4) zusammen, so erh¨ alt man die Formel

in der die µ i die Schw¨ achungskoeffizienten der einzelnen Schichten bezeichnen. Verallgemeinern wir dies nun auf eine beliebige Aufeinanderfolge von beliebig dicken Materialschichten. Es sei weiterhin eine Gerade L so parametrisiert, daß L(t) f¨ ur jedes t ∈ (−∞, ∞) der Punkt auf L ist, der t L¨ angeneinheiten von einem beliebigen aber festen Referenzpunkt - ebenfalls auf L - entfernt ist. Dann ist, wenn µ(L(t)) gerade die Schw¨ achung des Materials an der Stelle L(t) angibt, die Schw¨ achung des gesamten Strahls entlang L zu berechnen ¨ uber das Integral

In der Realit¨ at mißt man diese Gr¨ oße nicht direkt, sondern unter Ausnutzung von Gleichung (1.3) ¨ uber die Durchgangswahrscheinlichkeit P von Photonen. Allerdings kann auch P nur indirekt gemessen werden, indem etwa ein Referenzdetektor die von der R¨ ontgenquelle in eine andere Richtung emittierten Photonen unged¨ ampft z¨ ahlt. Das Verh¨ altnis der von beiden Detektoren gemessenen Photonenraten ist dann ein Maß f¨ ur die Schw¨ achung des betrachteten Strahls.

Diese Art der Messung birgt in sich aber bereits eine wesentliche Fehlerquelle, da die Emission und Absorption von Photonen stochastische Prozesse sind. Dies wird in Abschnitt 2.1 n¨ aher beleuchtet werden. Dort findet sich auch eine detailliertere Beschreibung des Meßprinzips.

1.2 Zweidimensionale Scanner

Man beschr¨ ankt sich h¨ aufig auf Strahlen in einer zweidimensionalen Ebeneder Rekonstruktionsebene - und definiert in dieser einen geeigneten Punkt als Rekonstruktionszentrum. In der Rekonstruktionsebene sind alle Strahlen eindeutig identifiziert durch ihren Abstand vom Rekonstruktionszentrum und ihren Winkel θ zu einer festen Referenzgeraden, etwa der positiven y-Achse.

Betrachten wir die Verteilung der Schw¨ achung des Materials in der Rekonstruktionsebene als eine Funktion f (r, φ) in den Polarkoordinaten (r, φ). Dann ” berechnet“ jeder durch

(, θ) beschriebene Strahl ² bei seinem Durchgang durch das betrachtete Objekt das folgende Integral:  √

² Nat¨ urlich nur im Idealfall. Es m¨ ussen zum Beispiel hinreichend viele Photonen betrachtet und St¨ orungen vernachl¨ assigt werden.

1.2. Zweidimensionale Scanner 3

Gleichung (1.7) definiert die sogenannte Radontransformierte der Funktion f . Tats¨ achlich k¨ onnen nat¨ urlich nur endlich viele konkrete Werte dieser Funktion gemessen bzw. in der Simulation berechnet werden. Diese werden nachfolgend als Rohdaten bezeichnet.

In der Praxis erfolgt die Sammlung der Rohdaten mit Hilfe einer Anordnung von einer oder mehreren R¨ ontgenquellen auf der einen und R¨ ontgendetektoren auf der anderen Seite. Zwischen Strahlquelle(n) und Detektoren wird das zu untersuchende Objekt plaziert. W¨ ahrend des Meßvorganges wird die Meßanordnung schrittweise um das Objekt gedreht. F¨ ur jede Drehstellung mißt jeder Detektor die Schw¨ achung des ihm zugeordneten Strahls.

Die so f¨ ur jeden Drehwinkel getrennt gesammelten Daten werden als Sicht (engl. view ) bezeichnet. Eine Sicht ist also der R¨ ontgenschatten des durchstrahlten Objekts unter einem festen Betrachtungswinkel.

Wie die Strahlen einer Sicht geometrisch zueinander in Beziehung stehen ist dabei durch den Aufbau des verwendeten Scanners bestimmt. Haupts¨ achlich kommen die zwei nachfolgend beschriebenen Grundtypen zum Einsatz.

1.2.1 F¨ acherstrahlscanner

Ein Scanner mit einer R¨ ontgenquelle und einem Satz von Sensoren auf einem Kreisbogen sammelt bei einem vollen Umlauf der Strahlenquelle um das gescannte Objekt f¨ ur jeden Drehwinkel der Anordnung die Schw¨ achung eines divergierenden Strahlenb¨ undels. Aufgrund der Form des Strahlenb¨ undels wird diese Bauart als F¨ acherstrahlscanner (engl. fanbeam) bezeichnet.

Die Strahlen einer Sicht schneiden sich also - wie nicht anders zu erwarten - alle an der Position der Strahlenquelle. Jeder Meßwert wird ¨ uber den Drehwinkel θ des mittleren

Strahls bzw. der gesamten Scanneranordnung zur positiven y-Achse und den Winkel α des jeweiligen Detektors identifiziert. Der Winkel zwischen benachbarten Strahlen ist bei dieser Art Scanner im allgemeinen konstant.

Die Geometrie eines F¨ acherstrahlscanners ist charakterisiert durch den Abstand der Strahlenquelle vom Drehmittelpunkt, die Anzahl der Detektoren, den ¨ Offnungswinkel

des Strahlf¨ achers sowie den Rotationswinkel der Meßanordnung.

1.2.2 Parallelstrahlscanner

Ein Parallelstrahlscanner arbeitet mit einem Satz paralleler Strahlen, wobei benachbarte Strahlen im allgemeinen konstanten Abstand voneinander haben. Hierzu wird entweder eine R¨ ontgenquelle in jeder Drehstellung zus¨ atzlich zur Drehung linear verschoben oder f¨ ur jeden Strahl eine getrennte Strahlenquelle verwendet. Letzteres beschleunigt den Vorgang der Datengewinnung bei realen Scannern nat¨ urlich immens. Bei dieser Bauart ist nur eine Drehung der Meßanordnung um 180 Grad notwendig, da bei einer weiteren Drehung nur bereits gemessene Strahlen in umgekehrter Richtung erneut gemessen w¨ urden.

1.3. Wahl des Rekonstruktionsbereiches 5

Die Meßwerte bzw. die zugeh¨ origen Strahlen werden gekennzeichnet durch den Drehwinkel θ der Meßanordnung sowie den Abstand des jeweiligen Strahls vom Drehmittelpunkt. Offensichtlich ist das genau die gleiche Parametrisierung, wie sie von Gleichung 1.7 verwendet wird.

Ein Parallelstrahlscanner wird zum einen charakterisiert durch den Rotationswinkel der Meßanordnung. Dieser betr¨ agt zumeist 180 Grad. Des weiteren sind entscheidend der Abstand der beiden ¨ außersten Strahlen, der bei jeder Drehstellung der gleiche ist, sowie der Abstand zweier benachbarter Strahlen bzw. deren Anzahl.

1.3 Wahl des Rekonstruktionsbereiches

Der Rekonstruktionsbereich ist der Teil der Rekonstruktionsebene, in dem sich ein Objekt befinden darf, um seine Dichteverteilung rekonstruieren zu k¨ onnen. In den hier betrachteten Simulationen ist dieser Bereich grunds¨ atzlich der Kreis, der durch die ¨ außersten Strahlen w¨ ahrend einer Messung tangiert wird. Das Zentrum dieses Kreises ist das Rekonstruktionszentrum. Dieses ist offensichtlich mit dem Mittelpunkt des Kreises identisch, der von der Meßanordnung w¨ ahrend eines vollen Umlaufs beschrieben wird. Zweckm¨ aßigerweise verwenden wir nachfolgend das Rekonstruktionszentrum zugleich als Koordinatenursprung. Der Radius E des Rekonstruktionsbereiches h¨ angt ab von Typ und Abmessungen des verwendeten Scanners. F¨ ur einen Parallelstrahlscanner mit maximalem Strahlabstand max ist offensichtlich gerade E = max , (1.8)

w¨ ahrend im Falle eines F¨ acherstrahlscanners mit F¨ acher¨ offnungswinkel α und Umlaufradius D der Strahlenquelle gilt α

E = D sin . (1.9)

2

1.4 Phantome

Um die Datengewinnung simulieren zu k¨ onnen, ist die Kenntnis der bei der Messung vorhandenen Dichteverteilung in der Rekonstruktionsebene erforderlich. Es ist hierbei zu beachten, daß die Funktionsf¨ ahigkeit der sp¨ ater betrachteten Rekonstruktionsverfahren m¨ oglicherweise von der genauen Form dieser Verteilung abh¨ angt. Ein Verfahren, das etwa einen Sch¨ adelschnitt sehr exakt rekonstruiert, kann bei Betrachtung eines anderen K¨ orperbereiches durchaus v¨ ollig versagen. Dies ist umso wahrscheinlicher, da die Rohdaten - wie bereits erw¨ ahnt - nur eine endliche Teilmenge aller Strahlsummen repr¨ asentieren, also notwendigerweise unvollst¨ andig sind. Man erwartet allerdings, daß geringe Unterschiede in der Dichteverteilung auch nur geringen Einfluß auf die erzielte Rekonstruktionsqualit¨ at haben. Daher werden in der Simulation k¨ unstliche Verteilungen verwendet, die den praktisch relevanten Objekten m¨ oglichst ¨ ahnlich sind. Solche Verteilungen heißen Phantome.

6

Als Objekt verwenden alle hier beschriebenen Tests das in Abb. 1.3 gezeigte Phantom eines Sch¨ adelschnittes mit eingebetteten Strukturen. Die Grauwerte des Phantoms stellen die Schw¨ achung von R¨ ontgenquanten einer bestimmten aber frei gew¨ ahlten Wellenl¨ ange dar. Je h¨ oher der Grauwert, um so h¨ oher ist dabei die Wahrscheinlichkeit, daß ein Quant beim Durchgang durch den betreffenden Pixel absorbiert oder gestreut wird. Ein Grauwert von 0 (schwarz) bedeutet ”

der Photonen. Dies entspricht - etwas idealisiert - der Schw¨ achung durch die das Objekt umgebende Luft.

Das verwendete Phantom wird aus einer Phantombeschreibungsdatei erzeugt. Deren Format ist in Anhang A.1 beschrieben. Das oben abgebildete Sch¨ adelphantom ist definiert durch die folgende Beschreibung: # Phantombeschreibungsdatei #

# Verwendete Basisobjekte und Parameter: # # ellipse # triangle # segment # # Beginn Kopfphantom ellipse ellipse ellipse ellipse ellipse # Zwei Sicheln segment segment segment segment # Haematom segment 1371 328 1102 434 94 0.263 0 segment 1395 328 1125 445 352 0.391 0

1.5. Simulation der Datengewinnung 7

# segment segment triangle

triangle # # Ende Kopfphantom Mit # beginnende Zeilen sind Kommentare. Alle Koordinaten und Maße verstehen sich in 1 10 mm. Der Radius des Rekonstruktionsbereiches ist fest auf 30cm eingestellt. Die Einheit der Schw¨ achungen (Parameter value) ist cm −1

. 1.5 Simulation der Datengewinnung

Der Vorgang der Datengewinnung vollzieht die Berechnung“ der

”

Radontransformierten (1.7) durch R¨ ontgenstrahlen idealisiert nach. Man betrachte hierzu einen beliebigen Strahl, der durch die Parameter (, θ) gem¨ aß Abschnitt 1.2 beschrieben ist. Weiterhin seien die Koordinatenfunktionen r ,θ (t) und φ ,θ (t) so definiert, daß, wenn t von −∞ bis ∞ l¨ auft, (r ,θ (t), φ ,θ (t)) gerade den Strahl (, θ) uberstreicht. Dann heißt das Integral ¨

die Strahlsumme des Strahls (, θ).

Da in der Praxis - wie auch in der Simulation - nur endlich viele Strahlen betrachtet werden k¨ onnen, identifizieren wir diese durch ( i , θ i ) f¨ ur 0 ≤ i < I und arbeiten im folgenden nur noch mit der Koordinatendiskretisierung p(nd, m∆).

Es ist also = nd und θ = m∆. Hierbei bezeichnet d den Abstand zweier benachbarter Strahlen ³ , ∆ den Drehwinkel zwischen der Aufnahme zweier aufeinanderfolgender Sichten. Als Wertebereiche dieser Indizes w¨ ahlt man n ∈ [−N..N ] sowie m ∈ [0..M − 1].

Daraus erh¨ alt man bei Beschr¨ ankung auf den Rekonstruktionsbereich mit dem Radius E f¨ ur den Definitionsbereich von p die Gleichungen

Der maximale Strahlabstand vom Rekonstruktionszentrum ist also gerade E, was nicht uberrascht, da E in Abschnitt 1.3 gerade so festgelegt wurde. Sichtenwinkel gr¨ oßer als π ¨

brauchen deshalb nicht ber¨ ucksichtigt zu werden, weil die Orientierung eines Strahls f¨ ur die Berechnung seiner Strahlsumme irrelevant ist.

³ Diese Notation setzt implizit einen Parallelstrahlscanner voraus. Die Darstellung der Strahlen beim F¨ acherstrahlscanner erfolgt analog.

8

Numeriert man nun die Pixel im abgetasteten Bild mit j im Bereich 0 ≤ j < J, so wird aus dem Integral in Gleichung (1.10) die Gleichung f¨ ur die diskrete Strahlsumme

Hierbei sind (r j , φ j ) die polaren Koordinaten des j-ten Pixels und δ j (nd, m∆) die Schnittl¨ ange des Strahls (nd, m∆) mit dem j-ten Pixel. Es ist also insbesondere δ j (nd, m∆) = 0 f¨ ur alle Pixel j, die vom betrachteten Strahl nicht geschnitten werden.

Beispielrohdaten eines Parallelstrahlscanner zeigt Abbildung 1.4. Jede Zeile enth¨ alt die Daten einer Sicht, also den R¨ ontgenschatten des Sch¨ adelphantoms unter einem Drehwinkel der Meßanordnung. Die erste Bildzeile stellt die von den einzelnen Detektoren gemessenen Schw¨ achungen bei senkrecht von oben nach unten durch das Phantom laufenden Strahlen dar. Nach unten wird die Meßanordnung gegen den Uhrzeigersinn bis auf 180 Grad gedreht. Die Punkte einer Zeile entsprechen den gemessenen Schw¨ achungen der einzelnen Strahlen der jeweiligen Sicht. Die senkrechte Koordinatenachse des Rohdatenbildes ist also mit θ ∈ [0, π] beschriftet, die waagerechte mit ∈ [−E, E]. Schreibt man die Pixel aus Originalbild und Rohdaten aus Abb. 1.4 jeweils zeilenweise hintereinander, so kann man die Dichteverteilung f (r j , φ j ) und die gemessenen

Strahlsumme (1.13) - im Vorgriff auf sp¨ atere ¨

Dichteverteilung betreffend - in der Form

f (1.14)

schreiben. Die Elemente r ij der sogenannten Projektionsmatrix R beschreiben hierbei die L¨ ange des Schnittes des i-ten Strahls mit dem j-ten Pixel. R ist also eine I × J-Matrix, die abh¨ angig ist von der Geometrie des verwendeten Scanners.

1.5. Simulation der Datengewinnung 9

Die Schreibweise in Gleichung (1.14) hat dabei einen entscheidenden Vorteil im Vergleich zu (1.7): Sie ist nicht festgelegt auf lineare Strahlen.

Betrachtet man statt R¨ ontgenstrahlen etwa Elektronenstrahlen in einem magnetischen Feld, so sind deren Bahnen keine Geraden, sondern Kreisb¨ ogen. Ihre Schw¨ achungen stellen deshalb nicht mehr die Radontransformierte der Dichteverteilung in der Rekonstruktionsebene dar. Die Projektionsmatrix l¨ aßt sich aber auch hier angeben. Rekonstruktionsverfahren, die von der Darstellung des Problems als lineares Gleichungssystem ausgehen, sind also flexibler als solche zur Berechnung der zu (1.7) inversen Transformation.