Die Computertomographie erm¨oglicht es, Querschnitte eines lebenden menschlichen K¨or-pers graphisch darzustellen. Im Gegensatz zu einer herk¨ommlichen R¨ontgenuntersuchung, bei der man die Projektion des Objekts auf eine Ebene betrachtet, werden bei der Com-putertomographie die R¨ontgenstrahlen in der Querschnittsebene aus vielen verschiedenen Richtungen durch den K¨orper geschickt. Die Daten liefern um den Patienten rotierende R¨ontgenquellen und Detektoren (Abbildung 3.10). Der K¨orper absorbiert die
R¨ontgen-D
D Q
Q
Abbildung 3.10: Prinzipieller Aufbau eines Computertomographen mit R¨ontgenquellenQ und Detektoren D.
strahlung unterschiedlich. Aus der Abschw¨achung der Strahlung k¨onnen dann die Ab-schw¨achungskoeffizienten in der Querschnittsebene bestimmt und damit zwischen dich-ter und weniger dichdich-ter Madich-terie (z.B. Knochen versus Gewebe) undich-terschieden werden.
Ubrigens beschrieb G.N. Hounsfield (zusammen mit J. Ambrose) 1972 das Prinzip der¨ Computertomographie als erster und erhielt daf¨ur 1979 den Nobelpreis der Medizin.
Wir wollen in diesem Abschnitt die Frage beantworten:
Frage: Wie kann das Schnittbild aus den Daten der R¨ontgendetektoren rekonstruiert werden?
Die Schw¨achung der R¨ontgenstrahlung beim Durchgang durch Materie geschieht nach dem Lambert-Beer-Gesetz
I(d) =I(0)e−µd, (3.18)
wobei I(d) die Intensit¨at des Strahles am Ort d und d die Dicke des homogenen Medi-ums seien. Die Konstante µ > 0 ist der Abschw¨achungskoeffizient und ist ein Maß f¨ur die St¨arke der Abschw¨achung. Genauer gesagt ist µ eine Funktion des Materials und der Energie der R¨ontgenstrahlung. Das Gesetz (3.18) kann mit der Tunnelwahrscheinlichkeit aus Abschnitt 3.2 verglichen werden, die ebenfalls exponentiell von der Dicke des Hin-dernisses abh¨angt. Wegen der unterschiedlichen Materialien im menschlichen K¨orper ist der Koeffizientµeine Funktion des Ortes. Die Abschw¨achung der Strahlungsintensit¨atdI entlang eines kleinen Wegst¨uckes dx ist dann proportional zu der Intensit¨at:
dI =−µ(x)I dx
oder dI
dx =−µ(x)I.
Integration dieser Differentialgleichung
−µ(x) = 1 I
dI dx = d
dxlnI
¨uber einen Weg Γ von einer Quelle am Ort Q zu einem Detektor am Ort D liefert dann Z
Γ
µ(x)dx=−lnI(D) + lnI(Q) =−lnI(D)
I(Q). (3.19)
Unser Ziel lautet, den ortsabh¨angigen Schw¨achungskoeffizienten µ eines Querschnitts zu bestimmen. Die Apparatur rotiere in nα gleichverteilten Schritten zwischen αmin und αmax, und jedes Strahlenb¨undel bestehe aus np gleichverteilten Strahlen zwischen pmin
undpmax(Abbildung 3.10). Die Messungen mit dieser Geometrie liefernm :=nαnp Werte der nαnp Linienintegrale (3.19) ¨uber µ(x) entlang der Geradenst¨ucke Γ1, . . . ,Γm von den Quellen zu den Detektoren. Bezeichnen wir mitQiundDi die durch dasi-te Geradenst¨uck Γi verbundenen Quellen und Detektoren und mit I(Qi) und I(Di) (i = 1, . . . , m) die zugeh¨origen Intensit¨aten, so gilt nach (3.19)
Z
Γi
µ(x)dx=bi, i= 1, . . . , m, (3.20) wobei
bi =−lnI(Di) I(Qi)
die experimentell ermittelten Werte sind. Die Gleichung (3.20) ist die Radontransforma-tion der Funktion µ. Bei (geeignet gew¨ahlten) unendlich vielen Strahlen ist µ eindeutig durch die rechten Seiten bestimmt. Bei endlich vielen Strahlen l¨aßt sich µ nur approxi-mativ bestimmen.
Eine M¨oglichkeit, µ zu bestimmen, ist die Diskretisierung der Umkehrformel der Ra-dontransformation. Wir w¨ahlen im folgenden einen sehr einfachen algebraischen Ansatz aus [5]: Wir suchen eine N¨aherungsl¨osung µ des Systems (3.20). Dazu zerlegen wir die Querschnittsebene in n2 gleich große quadratische Elemente, im folgenden Pixel genannt (Abbildung 3.11). Wir suchen eine Bestapproximation an die L¨osung von (3.20) im n2 -dimensionalen Raum
S = span{φ1, . . . , φn2} mit den Basisfunktionen φi :R2 →R, definiert durch
φ(k−1)n+`(x) =
( 1 : falls x im Pixel der k-ten Zeile und `-ten Spalte liegt 0 : sonst.
n 3 2 1
2
D P
Pmax
min Q
α
Abbildung 3.11: Zur Approximation des Abschw¨achungskoeffizienten µ.
Die Funktionenφi sind also die Indikatorfunktionen dern2 Pixel. Die N¨aherungsl¨osungµ ist als Linearkombination aller φi dargestellt:
µ(x) =
n2
X
i=1
µiφi(x). (3.21)
Dies bedeutet, daß µ auf jedem Pixel der Zerlegung konstant ist. Das medizinische Bild entsteht dann durch Darstellung der Werteµials Grauwerte. Das Bild wird umso genauer, je mehr Pixel verwendet werden.
Setzen wir (3.21) in (3.20) ein, so folgt n¨aherungsweise bi =
Z
Γi
µ(x)dx=
n2
X
j=1
µj
Z
Γi
φj(x)dx, i= 1, . . . , m. (3.22) Setzen wir
aij = Z
Γi
φj(x)dx, u= (µ1, . . . , µn2)>, b = (b1, . . . , bm)>
und A= (aij), k¨onnen wir (3.22) als lineares Gleichungssystem
Au=b (3.23)
mit der Unbekanntenuschreiben. Die Zahl der Messungen (d.h. die Zahl der Gleichungen in (3.22)) wird im allgemeinen wesentlich gr¨oßer als die Anzahl der Pixel sein, um die Qualit¨at der Rekonstruktion zu erh¨ohen. In diesem Fall ist (3.23) ein ¨uberbestimmtes
Gleichungssystem und daher im allgemeinen nicht l¨osbar oder nicht eindeutig l¨osbar. Wir suchen deshalb eine L¨osung im Sinne der kleinsten Quadrate, d.h., u ist die L¨osung von
Q(u) = min
v∈Rn2Q(v) mit Q(v) := kAv−bk22, (3.24) und kxk22 :=x>x ist die euklidische Norm auf Rn2.
Wie k¨onnen wir nun (3.24) l¨osen? Wir behaupten, daß das Minimierungsproblem ¨aqui-valent ist zur L¨osung derGaußschen Normalengleichungen
A>Au=A>b. (3.25)
Satz 3.3 Seien A ∈ Rm×n2 und b ∈ Rm. Dann sind die Probleme (3.24) und (3.25)
¨aquivalent.
Beweis: Wir berechnen diei-te Richtungsabteilung von Q:
1
s(Q(v+sei)−Q(v)) = 1
s(Av−b+sAei)>(Av−b+sAei)− 1
s(Av−b)>(Av−b)
= 2(Av−b)>Aei+s(Aei)>Aei,
wobei ei der i-te Einheitsvektor desRm sei. Der Grenzwert s→0 ergibt
∂Q(v)
∂ui
= 2(Av−b)>Aei und damit
∇uQ(v) = (2(Av−b)>A)> = 2A>(Av−b).
Ist nun u eine L¨osung von (3.24), gilt ∇uQ(u) = 0, und u l¨ost (3.25). Ist umgekehrt u eine L¨osung von (3.25), so ist u ein kritischer Punkt von Q. Die Funktion u ist ein Minimum vonQ, weilA>Awegenv>(A>A)v =kAvk22 ≥0 f¨urv ∈Rn2 positiv semidefinit
ist. ¤
Wir bemerken, daß das lineare Gleichungssystem (3.25) im allgemeinen nur iterativ gel¨ost werden kann. ¨Ublicherweise werden n2 = 5122 Pixel zur Bestimmung des Quer-schnittsbildes verwendet. Verwendet man nun m > n2 = 262144 R¨ontgenstrahlen, ist ein riesiges Gleichungssystem mit einer Matrix mit mindestens 262144×262144 Elementen zu l¨osen. Die Koeffizientenmatrix enth¨alt sehr viele Nullen, da jeder R¨ontgenstrahl nur eine kleine Anzahl von Pixeln schneidet. Doch selbst wenn nur die von Null verschiedenen Eintr¨age abgespeichert werden, erfordert die Matrix einen Speicherbedarf von etwa zwei Gigabyte.
4 Str¨ omungen
Eine große Anzahl von sich bewegenden Teilchen kann durch makroskopische Variablen wie Teilchendichte und mittlere Teilchengeschwindigkeit beschrieben werden. R¨aumliche und zeitliche ¨Anderungen dieser Gr¨oßen gen¨ugen h¨aufig Erhaltungsgleichungen, die formal zu partiellen Differentialgleichungen ¨aquivalent sind. In diesem Kapitel betrachten wir verschiedene Str¨omungen, die auf Erhaltungsgleichungen f¨uhren.