DE19711963C2 - Verfahren zur Bildrekonstruktion bei einem Computertomographen - Google Patents

Description

Es sind Computertomographen bekannt, bei denen ein Flächende­ tektor vorgesehen ist, der aus einer Matrix von Detektorele­ menten besteht. Er ist demgemäß aus einer Reihe von Zeilen aufgebaut, wobei jede Zeile eine Reihe von Detektorelementen enthält. Der Detektor wird von einem pyramidenförmigen Rönt­ genstrahlenbündel getroffen und erlaubt je nach Ausdehnung des Flächendetektors in z-Richtung, d. h. in Richtung der System­ achse, um die das Meßsystem aus Röntgenstrahler und Detektor gedreht wird, die simultane Abtastung mehrerer Schichten. Zur Abtastung eines Volumens des Untersuchungsobjektes rotiert das Meßsystem kontinuierlich um die Systemachse, während eine Re­ lativbewegung des Untersuchungsobjektes gegenüber dem Meßsy­ stem in Richtung der Systemachse erfolgt.

Für die Bildrekonstruktion aus Daten eines Spiralscanners mit einer einzigen Detektorzeile sind in der Literatur verschiede­ ne Verfahren bekannt. Beispielsweise wird in der EP 0 430 549 A2 ein Verfahren angegeben, bei dem eine Interpolation in Rich­ tung der z-Achse durchgeführt wird, indem man die gewichtete Summe der beitragenden Strahlen berechnet, wobei die Gewichte vom Abstand der Strahlen zur Bildebene abhängen. Bei einem Mehrzeilendetektor ist der Abstand eines Strahls zur Bildebene wegen der Neigung der Strahlen um ihren Cone-Winkel nicht mehr über den ganzen Strahl konstant, sondern wird zu einer Funktion der betrachteten Position längs des Strahls. Dadurch wird es erforderlich, die Spiralinterpolation bei der Rückprojektion für jedes Voxel gestreut durchzuführen. Die dazu erforderliche Kegelstrahl-Rückprojektion hat aber einen sehr hohen Rechen­ aufwand. In der DE 44 38 988 A1 wird anstelle einer Kegel­ strahl-Rückprojektion eine Gewichtung der einzelnen Detektorzeile mit anschließender gewöhnlicher Rückprojektion vorgeschlagen. Dabei wird jedoch der Cone-Winkel der Strahlen vernachlässigt und es ergibt sich ein Verlust an z-Schärfe und eine Zunahme von Cone-Beam-Artefakten.

In der US 5 406 479 A ist ein Verfahren zur Vermeidung des durch den Cone-Winkel hervorgerufenen Fehlers bei der Bildrekonstruk­ tion für einen (Einzeilen-)Computertomographen im Spiralbe­ trieb beschrieben, bei dem die vom Detektor gelieferten Daten zunächst von Fächer- auf Parallelgeometrie umsortiert und um­ interpoliert werden und anschließend ein Bild aus der Menge der erzeugten Parallelprojektionen berechnet wird, wobei das Bild so definiert ist, daß jeder Strahl mit einem Gewicht zum Bild beiträgt. Dabei wird die Berechnung des Bildes durch Faltung der transformierten Parallelprojektionen und Gewichtungsfunk­ tion im Fourierraum und anschließender Rückprojektion durchge­ führt.

Der Erfindung liegt die Aufgabe zugrunde, ein Verfahren zur Bildrekonstruktion bei einem Mehrzeilendetektor-Computertomo­ graphen im Spiralbetrieb anzugeben, das eine genaue Bildrekon­ struktion innerhalb kurzer Zeit ermöglicht.

Diese Aufgabe ist erfindungsgemäß gelöst durch die Merkmale des Patentanspruches. Das erfindungsgemäße Verfahren ist ein ap­ proximatives Verfahren und stellt eine sehr effiziente Imple­ mentierung als Fourier-Rekonstruktion dar. Dieses Vorgehen vermeidet die nicht korrekte Vertauschung von z-Interpolation und Faltung in Gleichung 21 und erlaubt die Berücksichtigung der Neigungswinkel der Strahlen (Cone-Winkel), ohne daß eine aufwendige Kegelstrahl-Rückprojektion erforderlich wäre.

Die Erfindung ist nachfolgend anhand der Zeichnung näher er­ läutert. Es zeigen

Fig. 1 die wesentlichen Teile eines Computertomographen zur Erläuterung des Erfindungsgedankens und Fig. 2 bis 4 geometrische Darstellungen zur Erläuterung des Bildrekonstruktionsverfahrens bei dem Computer­ tomographen gemäß Fig. 1.

In Fig. 1 ist mit 1 der Fokus eines Röntgenstrahlers 9 und mit 2 der Detektor eines Computertomographen bezeichnet. Der De­ tektor 2 besteht aus einer Reihe paralleler Zeilen, von denen jede eine Reihe von Detektorelementen enthält. Eine Primä­ strahlenblende blendet ein pyramidenförmiges Röntgenstrahlen­ bündel ein, das genau auf dem Detektor 2 auftrifft. Von diesem Röntgenstrahlenbündel sind in der Fig. 1 beispielhaft nur der Mittelstrahl 4 und ein weiterer Strahl 3 dargestellt. Zur Ab­ tastung eines in einem Meßfeld 5 liegenden Untersuchungsobjek­ tes rotiert das Meßsystem aus Röntgenstrahler 9 und Detektor 2 um eine Systemachse 6. Die dabei gebildeten Detektorsignale werden einem Rechner 7 zugeführt, der ein Bild des untersuchten Volumens des Patienten berechnet, das auf einem Monitor 8 wie­ dergegeben wird. Für die Volumenabtastung erfolgt eine Rela­ tivbewegung in Richtung der Systemachse 6 (z-Richtung) zwischen dem Meßsystem 1, 2, 9 und dem im Meßfeld 5 liegenden Untersuchungsobjekt.

Definition der Geometrie:

Fig. 1 erläutert schematisch die Geometrie des zugrundegeleg­ ten Mehrzeilendetektor-Computertomographen. Der Detektor 2 ist als Zylinderoberfläche ausgebildet. Der Radius dieses Zylin­ ders ist R_f + R_d. Der Fokus 1 liegt auf der Zylinderachse. Der Detektor 2 und der Strahler 9 sind starr miteinander verbunden. Während der Aufnahme bewegt sich der Strahler 9 relativ zum zu untersuchenden Objekt auf einer Spiralbahn mit Radius R_f. Der Abstand des Detektors 2 vom Drehzentrum sei R_d.

Der Projektionswinkel ist mit α bezeichnet. Innerhalb einer Projektion ist ein Strahl durch seine Parameter β und ζ_Det ein­ deutig bestimmt. Dabei ist β der Fächerwinkel des Strahls und ζ_Det die z-Position des zu ihm gehörenden Detektorelements re­ lativ zur z-Position des Fokus.

ζ_Det = z_Det - z_F(α). (1)

Die spiralförmige Fokusbahn hat die Steigung s und ist gegeben durch

z_F(α) = z_F,0 + α . s. (2)

z_F,0 ist die z-Position der Spirale für α = 0.

Da natürlich vom System nur diskrete Meßwerte aufgenommen wer­ den können, sind α, β und ζ_Det diskret:

Δα ist das Projektionswinkelinkrement zwischen aufeinander­ folgenden Projektionen. N_P,2π ist die Zahl der während eines Vollumlaufs aufgenommenen Projektionen, N_P ist die Zahl der insgesamt vorliegenden Projektionen. N ist die Zähl der Ka­ näle pro Zeile. Der Einfachheit halber beschränken wir die fol­ genden Betrachtungen auf den Fall geradzahliger N. AM ist das sog. Alignment im Fächerwinkel.

Δζ_Det ist der Abstand zweier Detektorzeilen in z-Richtung, N_rows ist die Zahl der Detektorzeilen und A Q ist ein mögliches Ali­ gnment in der z-Richtung. Die in das Drehzentrum skalierte Ko­ ordinate ζ_Det wird mit ζ bezeichnet. Mithin gilt

Der bei α_n, β_m und ζ_Det,q gemessene Strahl trägt den Meßwert p(α_n, β_m, ζ_Det,q) oder kurz p(n, m, q).

Fig. 2 zeigt die Projektion der Anordnung aus Fig. 1 in die x-y-Ebene. Die beiden Variablen α und β kennzeichnen eindeu­ tig die Projektion eines Strahls in diese Ebene. Den gleichen Zweck erfüllen auch die beiden Variablen θ und p, die mit α und β über die beiden folgenden Gleichungen verknüpft sind:

p = -R_fsinβ. (7)

Während α und β zur Beschreibung der Strahlen in Fächergeome­ trie verwendet werden, dienen die Parameter θ und p der Be­ schreibung der Strahlen in Parallelgeometrie, nach der Durchführung des im folgenden erklärten Rebinning, in dem von Fächer- auf eine spezielle Parallelgeometrie uminterpoliert wird.

Das Rebinning:

In diesem ersten Verarbeitungsschritt des Algorithmus werden die Mehrzeilen-Fächerdaten in Paralleldaten umsortiert und um­ interpoliert. Dieser Schritt wird für alle Detektorzeilen un­ abängig voneinander durchgeführt. Dabei wird zunächst die Schräglage der Strahlen sowie ihre jeweils verschiedene z-Po­ sition vernachlässigt, d. h. die Strahlen werden zunächst so be­ handelt als seien sie in einem zweidimensionalen CT-Gerät in einer Ebene z = const. gemessen worden. Erst später wird anhand der genauen geometrischen Lage und Orientierung der Original­ strahlen die den erzeugten Parallelstrahlen zuzuordnende Lage und Orientierung ermittelt.

Das Rebinning ist nachfolgend anhand der Fig. 2 genauer er­ klärt. Ein Strahl in Parallelgeometrie ist gegeben durch θ und p. In den Parallelprojektionen liegen Abtastwerte an den dis­ kreten Stellen

vor, wobei üblicherweise N^par = N und N par|p = N_p gesetzt wird, d. h. man berechnet die gleiche Zahl an Parallelprojektionen mit der gleichen Anzahl Kanäle wie in Fächergeometrie.

Der Begriff Parallelprojektion bezeichnet hier einen Satz par­ alleler Strahlen mit gleichem θ_l und q.

Das Rebinning ist im wesentlichen beschrieben durch (6) und (7). Auflösen dieser Gleichungen liefert eine Vorschrift für die Auswahl der Fächerstrahlen, die gesuchten Parallelstrahlen entsprechen. Man erhält:

Außerdem wird gesetzt:

q = const. (12)

Die gewünschten Abtastwerte in den Parallelprojektionen können anhand der Gleichungen (10), (11) und (12) errechnet werden.

Die Berechnung eines Paralleldatums p(Θ_l, p_k, p)oder kurz p(l, k, q) aus den Fächerdaten p(n, m, q) wird wie folgt voll­ zogen:

Zunächst werden aus Gleichung 10 und Gleichung 11 die Indizes , und des Fächerstrahls ermittelt, der dem Strahl l, k, q in Parallelgeometrie entspricht. Es ergibt sich:

Da und im allgemeinen nicht ganzzahlig sind, muß zwischen den nächsten Nachbarn interpoliert werden. Für den Spezialfall bilineare Interpolation erhält man mit

n_lo = floor(), n_hi = ceil()
m_lo = floor(), m_hi = ceil()

p(l, k, q) = (n_hi - )[(m_hi - )p(n_lo, m_lo, q) + ( - m_lo)p(n_lo, m_hi, q) + ( - n_lo)[(m_hi - ) p(n_hi, m_lo, q) + ( - m_lo)p(n_hi, m_hi, q)]

Im folgenden wird die Lage und Orientierung der erzeugten Par­ allelstrahlen bestimmt.

Das p-t-Koordinatensystem in Fig. 2 sei gegenüber dem x-y- System um den Winkel 0 - π um die z-Achse verdreht. p bezeichnet die Position eines Strahls innerhalb einer Parallelprojektion. t bezeichnet eine Position auf einem Strahl, d. h. t ist eine Koordinate in Strahllängsrichtung. Die Stelle t = 0 entspricht dem Schnittpunkt des Strahls mit einer Geraden, die sowohl auf dem Strahl als auch auf der z-Achse senkrecht steht, d. h. der Punkt t = 0 ist derjenige Punkt auf dem Strahl mit dem kleinsten Abstand zur z-Achse. Die z-Position eines Strahls bei t = 0 ist ein weiterer die Lage des Strahls beschreibender Parameter. Wir bezeichnen ihn mit z_t=0. Er ist gegeben durch

Dabei ist z_F(θ_l, p_k) die Fokus z-Position der Strahlen (θ_l, p_k). Man kann sie errechnen, indem man zunächst (10) in (11) einsetzt:

und dann (2) anwendet. Das Ergebnis ist:

Wegen des asin(.)-Terms ist die z-Position der Strahlen inner­ halb einer Parallelprojektion kanalabhängig. Allerdings haben alle Strahlen einer Parallelprojektion den gleichen Cone-Win­ kel (das ist der Winkel der Strahlen gegenüber der x-y-Ebene). Aus der Vielzahl der erzeugten Parallelprojektionen wird nun im folgenden ein Bild rekonstruiert.

Definition der Rekonstruktionsvorschrift:

Die Projektionen aller Strahlen der Parallelprojektionen zu den Winkeln θ = θ_lo + vπ (v beliebig ganzzahlig) in die Bildebene sind parallel, wobei die Strahlrichtung für geradzahlige v en­ gegengesetzt zu der für ungeradzahlige v ist. Bei der Berech­ nung des Bildes summieren wir daher nur über einen Halbumlauf (Projektionswinkelbereich θ_lo = θ...π d. h. l₀ = 0...(N par|P,π - 1)), be­ rücksichtigen aber für alle l₀ jeweils alle Parallelprojektio­ nen mit θ = θ_lo + vπ für alle v und q. Abb. 3 veranschaulicht ein Abtastmuster eines Mehrzeilendetektors als Schnitt entlang der z-Achse, d. h. die Zeichenebene ist die z- t-Ebene. Die dargestellten Strahlen sind die bei p_k = 0.

Ähnlich wie bei vielen anderen approximativen Verfahren soll in unserer Definition des zu rekonstruierenden Bildes das Ge­ wicht, mit dem ein Strahl zum Bild beiträgt, vom Abstand in z- Richtung des Strahls zur Bildebene abhängen. Dieser Abstand verändert sich jedoch im allgemeinen entlang eines Strahls, so­ daß das Gewicht, mit dem ein Strahl zur Bildebene beiträgt ebenso in Strahllängsrichtung variiert. Indem man für ein fe­ stes l₀, alle beitragenden Parallelprojektionen p(θ_lo + vπ, p_k, q) (für alle (v, q)) mit ihren zugehörigen Gewichtsfunktionen ver­ sieht und aufaddiert erhält man ein gewichtetes Mittel aus die­ sen Parallelprojektionen mit sich sowohl in Strahllängsrichtung als auch in Kanalrichtung verändernden Ge­ wichten. Diese gemittelte Parallelprojektion wird als Approxi­ mation einer unter dem Projektionswinkel θ_lo, vollständig in der Bildebene gemessenen Parallelprojektion verwendet. Die Gewich­ tungsfunktionen müssen aus Normierungsgründen so beschaffen sein, daß sich die Gewichte aller aufzusummierenden Strahlen jeweils zu 1 ergänzen. Die auf eine bestimmte Parallelprojek­ tion anzuwendende normierte Gewichtsfunktion h_n(.) ergibt sich, wenn man die abstandsabhängige Gewichtsfunktion h(.) durch die Summe aller beitragenden Gewichte teilt:

Weil die Orientierung der Parallelprojektionen für gerade bzw. ungerade v gerade entgegengesetzt ist, müssen für eine der bei­ den Fälle die beiden Koordinatenrichtungen p und t umgedreht werden. Dies wird durch r(v) bewerkstelligt. Es gilt:

Um die Notation zu vereinfachen, führen wir modifizierte Pro­ jektionen und Gewichtungsfunktionen ein, in denen für ungerad­ zahlige v die Koordinatenrichtungen umgekehrt sind.

Dabei ist der Abstand des Strahls in z-Richtung von der Bild­ ebene bei z_img gegeben durch

Unter Benutzung dieser gemittelten Parallelprojektion kann nun ein gewöhnliches Parallelrückprojektionsbild definiert werden:

Der Term in eckigen Klammern ist die gemittelte Parallelpro­ jektion.

Die Abhängigkeit von t in hn(.) berücksichtigt die Schräglage der Strahlen. Die Auswertung der vorigen Gleichung würde fol­ gende Schritte erfordern:

• - Gewichtete Summation über v und q für alle t
• - Faltung der gewichteten Summe mit L(.)
• - Zuschlag zum betrachteten Voxel x, y, z.

Da sich die t- und die p-Abhängigkeit in hn(.) nicht separie­ ren lassen, kann die Reihenfolge von Summation über (v, q) und Faltung über k nicht vertauscht werden. Es muß also strengge­ nommen für jedes t und damit für jedes Pixel eine eigene Fal­ tung durchgeführt werden. Der Rechenaufwand dafür ist aber für eine praktische Anwendung zu hoch. Führt man die Faltung den­ noch vor der gewichteten Summation aus, so kann die Interpola­ tion in z-Richtung bei der Rückprojektion erfolgen. Man begeht dazu allerdings einen Fehler. Das erfindungsgemäße Verfahren vermeidet das unzulässige Vertauschen der Reihenfolge von ge­ wichteter Summation und Faltung ohne einen hohen Rechenaufwand zur Folge zu haben. Das Verfahren kommt ohne die aufwendige Kegelstrahl-Rückprojektion aus und berücksichtigt dennoch die Cone-Winkel der Strahlen.

Durch periodische Wiederholung der gefalteten gemittelten Pro­ jektion erreicht man, daß das Spektrum in Kanalrichtung diskret wird. Das Bild ist im interessierenden Bereich von dieser pe­ riodischen Wiederholung nicht betroffen, wenn man davon aus­ geht, daß die gefalteten Projektionen endliche Ausdehnung haben und die Periode w entsprechend gewählt wird. Das Ergeb­ nis ist:

Innerhalb einer Parallelprojektion ist hn(.) eine für jeden Kanal p_k verschiedene eindimensionale Funktion von t. Für jede Bildebene z_img und für jeden Wert v, q und p_k erhalten wir also eine unterschiedliche Funktion hn(t,...), die wir im voraus be­ rechnen und abspeichern können. Statt eine Tabelle der Funkti­ onswerte abzuspeichern, approximieren wir die Gewichtsfunktion mit einer Fourierreihe und speichern die Fourierkoeffizienten. Wegen der speziellen Symmetrieeigenschaften der spiralförmigen Fokusbahn muß diese Tabelle nicht für verschiedene Projekti­ onswinkel abgespeichert werden. Eine Änderung des Projektions­ winkels entspricht nämlich einer Änderung von Z_img. Im folgenden ist erläutert, wie (22) im Frequenzbereich implemen­ tiert werden kann, was zu einer sog. Fourierrekonstruktion führt.

Implementierung im Frequenzbereich:

Zunächst approximieren wir die Gewichtungsfunktionen wie oben erwähnt durch Fourierreihen:

Dadurch werden die Gewichtsfunktionen periodisch. Die Periode Ω muß geeignet gewählt werden.

Setzt man dies in (22) ein, so erhält man zunächst:

In dieser Form läßt sich leicht das Spektrum des Bildes durch zweidimensionale Fouriertransformation berechnen. Man erhält:

Die Auswertung dieser Gleichung beinhaltet folgende Schritte:

Multiplikation der Projektionen p(θ_lo + vπ, p_k, q) mit den Gewich­ tungskoeffizienten c_µ(p_k, z_img, l₀, v, q) und Summation über alle bei­ tragenden (v, q).

FFT der aufsummierten, mit den Koeffizienten multiplizierten Projektionen entlang der p_k-Richtung, und zwar für alle µ.

Multiplikation des so erzeugten Spektrums mit (.), wiederum für alle µ.

Die obengenannten Schritte sind für alle l_o durchzuführen.

Abb. 4 zeigt das Spektrum einer einzelnen gemittelten Pro­ jektion unter einem bestimmten Projektionswinkel. Das Spektrum des Bildes ergibt sich durch Überlagerung aller dieser Beiträge über einen Halbumlauf.

Das resultierende Spektrum kann durch verschiedene Techniken in ein kartesisches Raster uminterpoliert (z. B. H. Schomberg, J. Timmer: "The Gridding Method for Image Reconstruction by Fourier Transformation,", IEEE Transactions on Medical Ima­ ging, Vol. 14, No. 3, September 1995) und dann mittels zweidi­ mensionaler inverser FFT in den Ortsraum zurücktransformiert werden.

Claims (1)

1. Approximatives Verfahren zur Bildrekonstruktion bei einem Mehrzeilendetektor-Computertomographen im Spiralbetrieb, bei dem die vom Detektor (2) gelieferten Daten zunächst für jede Detektorzeile getrennt von Fächer- auf Parallelgeometrie um­ sortiert und uminterpoliert werden, anschließend ein Bild aus der Menge der erzeugten Parallelprojektionen berechnet wird, wobei das Bild so definiert ist, daß jeder Strahl (3, 4) mit einem Gewicht zum Bild beiträgt, das von seinem jeweiligen, sich in Strahllängsrichtung ändernden Abstand zur Bildebene abhängt, der durch diese Abstandsabhängigkeit entstehende Ge­ wichtsverlauf jedes Strahls (3, 4) in Strahllängsrichtung schließlich durch eine Fourierreihe mit den Koeffizienten c_µ, µ = -N_µ...N_µ approximiert wird, indem zunächst jeder Kanal einer jeden beitragenden Parallelprojektion mit den zugehörigen c_µ multipliziert wird, sodann über alle beitragenden Parallel­ projektionen summiert und anschließend für alle µ eine FFT in Kanalrichtung durchgeführt und das Ergebnis für alle µ mit dem Spektrum eines Faltungskerns (.) multipliziert wird, wo­ bei das sich ergebende Spektrum der zu einem Projektionswin­ kel gehörenden Bildkomponente damit nur Werte auf diskreten Stützstellen hat, die auf 2N_µ + 1 parallelen Linien liegen, und das Gesamtspektrum sich durch Summation über alle Projek­ tionswinkel ergibt und schließlich durch bekannte Gridding- Techniken auf ein kartesisches Frequenzraster uminterpoliert und anschließend in den Ortsbereich (Bild) zurücktransfor­ miert wird.