5. Diskretisierung

Neben der Diskretisierung der Umkehrformeln der Radontransformation (4) ergibt sich eine andere Möglichkeit durch folgenden algebraischen Ansatz: Die Funktion aus (3) beschreibt auf dem quadratischen Bildgebiet ², welches das zu rekonstruierende Objekt enthält, die Dichteverteilung der Materie. Auf ²\ ² soll verschwinden. Mit Hilfe der Methode der finiten Elemente wollen wir eine Näherungslösung des Systems (4) berechnen. Dazu zerlegen wir wie in Abb. 5 skizziert das Quadrat ² in n² gleich große quadratische Elemente, im Folgenden Pixel genannt. Wir suchen eine Bestapproximation an die Lösung von (4) im n²-dimensionalen Raum

S = [₁ ,...,_n²]

mit den Basisfunktionen _i : ² --> , i = 1,...,n², die durch

(k-1)n+(x) =		1,	falls x im Pixel in der k-ten Zeile und	(5)
			-ten Spalte der Zerlegung liegt;
		0,	sonst.

(k, = 1,...,n) definiert sind. Die Funktionen _i, i = 1,...,n², sind also gerade die Indikatorfunktionen für die n² Pixel. Wir stellen die gesuchte Näherungslösung als Linearkombination der Basisfunktionen ₁ ,...,_n² dar:

	(x) =	n2	ii	(x)
		i = 1

Abbildung 5. Approximation des Objekts durch endlichdimensionalen Teilraum.

Bei dieser Wahl des Raumes S ist in jedem Pixel der Zerlegung konstant und |_{²\ ²} 0. Das Bild, welches der Mediziner später für seine Diagonse verwendet, entsteht durch Darstellung des Wertes im i-ten Pixel (also von _i) als Grauwert. Die Diskretisierung ist also dem Problem angepasst. In der Praxis sind zur Zeit Auflösungen von 512×512 Pixel gängig.

   Eingesetzt in (4) ergibt sich für i = 1,...,m

	Li		(x)dx =	n2	j		Li	j	(x)dx = bi.    (6)
				j = 1

Mit A = ((a_ij )), A ^m,n2, wobei

aij	=		Li	j	(x)dx  i = 1,...,m,  j = 1,...,n2,
	=	(1,..., n2)T n2,
b	=	(b1,...,bm)T m,

kann man (6) als lineares Gleichungssystem

Au = b    (7)

mit dem unbekannten Koeffizientenvektor und den gemessenen Werten der rechten Seite b schreiben. Die Formel für a_ij besagt, dass das (i, j)-te Element von A gerade die Länge des Schnitts des i-ten Strahls mit dem j-ten Pixel ist.

Die Zahl der Messungen (d. h. die Zahl der Gleichungen im linearen Gleichungssystem (7)) wird gewöhnlich wesentlich größer sein als die Anzahl der Pixel, um die Qualität der Rekonstruktion zu erhöhen. Somit ergibt sich ein überbestimmtes lineares Gleichungssystem sehr großer Dimension m×n², dessen Koeffizientenmatrix A aber schwach besetzt ist, d.h. nur wenige Einträge von A sind von Null verschieden, vgl. Abb. 6 und Tab. 2.

Abbildung 6. Struktur der Matrix A für n = 32, n = 36, nP = 42. Von Null verschiedene Elemente sind schwarz markiert.

n n nP m n2 nnz(A) MB 512 580 672 389760 262144 0.2% 1683.3 256 290 336 97440 65536 0.3% 210.8 128 290 336 97440 16383 0.6% 105.4 128 145 168 24360 16383 0.6% 26.4 64 290 336 97440 4096 1.2% 52.7 64 145 168 24360 4096 1.2% 13.2 32 290 336 97440 1024 2.3% 26.4 32 145 168 24360 1024 2.3% 6.6 32 72 84 6048 1024 2.3% 1.7

Tabelle 2. Werte für die Rekonstruktion aus Abschnitt 7. In der Spalte nnz(A) ist der Anteil der von Null verschiedenen Elemente in A, in der letzten Spalte der Speicheraufwand in Megabyte angegeben.

Es gilt die folgende Abschätzung:

Satz 1 In einer Zeile der Matrix A sind höchstens 2n-1 Elemente von Null verschieden.

Diese Aussage folgt durch einfache geometrische Überlegungen, die am Besten an Hand von Skizzen nachzuvollziehen sind. Wir geben dennoch auch einen formalen Beweis.

Beweis. Die Elemente der i-ten Zeile von A sind die Längen des Schnitts der Geraden L_i mit den Pixeln der Zerlegung von ². Ist L_i parallel zu einer der beiden Achsen, so schneidet L_i genau n Pixel, da wir eine Kante immer nur einem Pixel zuordnen. Es bleibt also nur der Fall L_i streng monoton zu untersuchen.
   Ohne Einschränkung sei L_i streng monoton fallend (den Fall L_i streng monoton wachsend behandelt man analog). Hat L_i eine Steigung größer als -1, so steht die Annahme L_i schneide drei oder mehr Pixel der j-ten Spalte der Zerlegung (1 j n) sofort im Widerspruch dazu, dass die Pixel quadratisch sind. Da j beliebig war, schneidet L_i in jeder Spalte der Zerlegung höchstens zwei Pixel. Es bleibt noch zu zeigen, dass es nicht möglich ist, dass in jeder Spalte zwei Pixel geschnitten werden. Dies beweisen wir durch Widerspruch: Angenommen, dies wäre der Fall. Nehmen wir also an, in der k-ten Spalte, 1 k n, wären dies Pixel und +1, 1 < n, dann müssen es in der (k+1)-ten Spalte Pixel +1 und +2 sein, usw. Offensichtlich bekommen wir die maximale Anzahl von Pixeln mit nicht leerem Schnitt mit L_i für k = 1 und = 1. In der n-ten Spalte würden dann die Pixel n und (n+1) geschnitten, im Widerspruch dazu, dass die n-te Spalte nur n Pixel hat. Es ist also nicht möglich, dass in der ersten und der letzten Spalte zwei Pixel von L_i geschnitten werden.
   Analog schneidet ein Strahl mit einer Steigung kleiner oder gleich -1 in jeder Zeile der Zerlegung höchstens zwei Pixel; jedoch niemals zwei Pixel in der ersten und in der letzten Zeile.
   Insgesamt haben also höchstens 2n-1 Pixel einen nichtleeren Schnitt mit L_i.

Das lineare Gleichungssystem (7) ist daher überbestimmt und i.A. nicht lösbar oder nicht eindeutig lösbar. Wir suchen deshalb eine Lösung im Sinne der kleinsten Quadrate (Carl F. Gauß, 1809), d.h.

Q() := ||b - A||² Q()  für alle ⁿ²     (8)

   Bekanntlich ist Lösung der Normalengleichungen

A^TAu= A^Tb.

Die Lösung ist genau dann eindeutig bestimmt, wenn A vollen Spaltenrang hat.

WeiterZurück

Auszüge aus dem Artikel "Mathematik fürs Leben am Beispiel der Computertomographie"
Autoren: Marlis Hochbruck · Jörg-M. Sautter

Kompletter Artikel (pdf)