DE19528261C1 - Verfahren zur Bildrekonstruktion bei einem Computertomographen - Google Patents

Description

In der Computertomographie können durch die Linogramm-Methode CT-Bilder rekonstruiert werden. Das Linogramm-Verfahren wurde 1987 von P. Edholm vorgestellt, siehe P. R. Edholm, G. T. Herman "Linograms in Image Reconstruction from Projections", IEEE Trans. on Med. Imag., Vol. MI-6, No. 4, 1987, Seiten 301-307. Die Implementierung ist beschrieben in P. Edholm, G. T. Herman, D. A. Roberts "Image Reconstruction from Linograms: Implementation and Evaluation", IEEE Trans. on Med. Imag., Vol. 7, No. 3, 1988, Seiten 239-246. Die Linogramm-Methode ist ein spezielles Fourier-Rekonstruktionsverfahren, das ganz ohne Interpolation im Frequenzraum auskommt. Ein Nachteil üblicher Fourierrekonstruktionsverfahren besteht nämlich darin, daß nach der eindimensionalen Fast-Fouriertransformation der gemessenen Projektionen f_lk diskrete Spektrumswerte auf einem Polargitter definiert sind. Der Übergang zum kartesischen Gitter im Frequenzraum geschieht durch eine zweidimensionale Interpolation, die bei ungenügender Sorgfalt zu fehlerbehafteten Bildern führt. Diese Interpolationsfehler können bei Anwendung des Gridding-Verfahrens beliebig klein gehalten werden, siehe J. D. O'Sullivan, "A Fast Sinc Function Gridding Algorithm for Fourier Inversion in Computer Tomography", IEEE Trans. on Med. Imag., Vol. MI-4, No 4, 1985, pp 200-207. Die wesentliche Idee der Linogramm-Methode dagegen ist, durch geeignete Maßnahmen von vornherein die diskreten Spektrumswerte auf ein "quasi-kartesisches" Raster zu legen, so daß keine Interpolation im Frequenzraum mehr nötig ist.

Bei der bekannten Linogramm-Methode erzeugt man immer ein dem Bild im gesamten Meßfeld mit dem Durchmesser D_M entsprechendes Spektrum. Das hat - wie bei den direkten Fourierrekonstruktionsverfahren - zur Folge, daß im Frequenzraum mit Matrizen sehr großer Punktzahl zu rechnen ist. Bei der originalen Linogramm-Methode entstehen im Winkelbereich [-45°, 45°] im Frequenzraum spektrale Werte mit festem Raster

in ρ_x-Richtung, im Winkelbereich [45°, 135°] dagegen spektrale Werte mit festem Raster

in ρ_y-Richtung (ρ_x bzw. ρ_y sind die kartesischen Frequenzvariablen). Hat man N_P=2048 Projektionen im Winkelbereich [-45°, 135°] mit N=1024 Meßwerten je Projektion im Abtastraster

so bedeutet dies, daß beide Teilmatrizen für [-45°, 45°] bzw. [45°, 135°] im Frequenzraum je

haben.

Soll das Verfahren praktische Bedeutung haben, müssen beliebige exzentrische Bildausschnitte rekonstruiert werden können. Beim bekannten Linogramm-Verfahren sind dazu ausgehend von den beiden Teilmkatrizen sowohl in ρ_x- als auch in ρ_y-Richtung chirp-z-Transformationen auf den gewünschten Bildausschnitt D_B · D_B mit N_Pix · N_Pix Bildelemeneten im Raster

durchzuführen. Dieser Schritt ist so aufwendig, daß er den Vorteil der fehlenden Interpolation im Frequenzraum weit überwiegt.

Der Erfindung liegt die Aufgabe zugrunde, ein Fourierrekonstruktionsverfahren zur Bildrekonstruktion für die Computertomographie so auszubilden, daß beliebige exzentrische Bildausschnitte direkt rekonstruierbar sind.

Diese Aufgabe ist erfindungsgemäß gelöst durch die Merkmale des Patentanspruches. Bei dem erfindungsgemäßen Verfahren handelt es sich um ein modifiziertes Linogramm-Rekonstruktionsverfahren, mit dem das Problem der Darstellung eines exzentrischen Ausschnittes eines CT-Bildes mit stark reduziertem Aufwand gelöst ist. Es handelt sich um eine Kombination aus Linogramm- und Gridding-Methode, mit der es möglich ist, beliebige exzentrische Bildausschnitte direkt zu rekonstruieren. Die dabei entstehende Matrix im Frequenzraum hat statt

wie beim originalen Linogramm-Verfahren das Frequenzraster

(D_B ist die Seitenlänge des gewünschten Bildausschnittes). Die Punkte liegen außerdem auf einem kartesischen Raster, nicht auf dem "quasi-kartesischen" Raster der Linogramm-Methode. Besteht das Bild aus N_Pix · N_Pix = 512 · 512 Bildelementen, so hat die kartesische Frequenzmatrix 2 N_Pix · 2 N_Pix = 1024 · 1024 Punkte - unabhängig von der Wahl von D_B. Das Gridding selbst ist - im Gegensatz zum originalen Gridding-Verfahren - nur eindimensional. Die Rücktransformation des Spektrums in den Ortsraum kann mit Fast-Fourier-Transformationen durchgeführt werden, chirp-z-Transformationen wie beim originalen Linogramm-Verfahren sind nicht nötig.

Die Erfindung ist nachfolgend anhand der Zeichnung näher erläutert.

Der in der Zeichnung dargestellte Computertomograph weist eine Meßeinheit aus einer Röntgenstrahlquelle 1, die ein fächerförmiges Röntgenstrahlbündel 2 aussendet, und einen Strahlempfänger 3 auf, welcher aus einer Reihe von Einzeldetektoren, z. B. aus 1024 Einzeldetektoren, besteht. Der Fokus ist mit 11 bezeichnet. Der zu untersuchende Patient 4 liegt auf einer Patientenliege 5. Zur Abtastung des Patienten 4 wird die Meßeinheit 1,3 um ein Meßfeld 9, in dem der Patient 4 liegt, um 360° gedreht. Die Drehachse ist mit 10 bezeichnet. Dabei wird die Röntgenstrahlenquelle 1, die von einem Röntgengenerator 6 gespeist wird, gepulst oder mit Dauerstrahlung betrieben. Bei vorbestimmten Winkelpositionen der Meßeinheit 1,3 werden Sätze von Daten erzeugt, die vom Strahlenempfänger 3 einem Rechner 7 zugeführt werden, welcher aus den erzeugten Datensätzen die Schwächungskoeffizienten vorbestimmter Bildpunkte berechnet und auf einem Sichtgerät 8 bildlich wiedergibt. Auf dem Sichtgerät 8 erscheint demgemäß ein Bild der durchstrahlten Schicht des Patienten.

Ausgangspunkt für die Rekonstruktion ist ein Paralleldatensatz mit N_P Projektionen im Winkelbereich [-45°, 135°] und N Kanälen je Projektion im Abtastraster a. Mißt ein CT-Scanner Daten in Fächergeometrie, seien die Paralleldaten durch Rebinning entstanden. Sei f_lk der Meßwert bei der Projektion l im Kanal k, dann gelte

ϑ_l = -45°+lΔϑ_l l = 0, 1, 2, . . . N_P-1 (1)

p_k = (k-0,5+a_m)a k = . . ., -1, 0, 1, . . . (2)

a_m ist das Alignment.

Weiter ist Δtan ϑ bzw. Δcot ϑ konstant

Das einzige Bild, das sich aus diesen Meßwerten ohne Zusatznahmen berechnen läßt, ist wie folgt definiert:

Dabei ist

und

M_A(ρ) ist die Modulationsübertragungsfunktion der Rekonstruktion. Gleichung (5) beschreibt das bekannte Faltungs-Rückprojektions-Verfahren mit Paralleldaten.

Weiterhin kann man setzen

L_C und L_S unterliegen nicht den Einschränkungen wie bei einer Faltungs-Rückprojektions-Rekonstruktion, wo L_C der CT-Faltungskern und L_S die bei der Rückprojektion notwendige Interpolationsfunktion ist, jedoch sei L_S eine im Ortsraum wenig ausgedehnte Glättungsfunktion.

Die Grundidee des vorgeschlagenen Verfahrens besteht darin, die notwendige Faltung der Meßwerte f_lk mit dem CT-Faltungskern im Frequenzraum als Multiplikation durchzuführen und dabei durch Multiplikation mit einem Phasenfaktor das Bild im gewünschten Rekonstruktionszentrum -_z- als neuem Ursprung zu zentrieren. Nach der eindimensionalen Rücktransformation der Meßwerte in den Ortsraum können die nunmehr gefalteten, zentrierten Projektionen auf das gewünschte Bildfeld D_B begrenzt werden. Daran schließen sich eindimensionale chirp-z-Transformationen auf das gewünschte Frequenzraster

für den Winkelbereich [-45°, 45°] bzw.

für den Winkelbereich [45°, 135°] an. Das wird im folgenden theoretisch begründet.

Wählt man das Rekonstruktionszentrum bei -_z, bildet man zunächst B₀ (-_z)

Betrachte die gefalteten, "zentrierten" Projektionen

f_c(p_q, ϑ_l) kann ohne Schaden für das Bild in einem Bereich mit dem Durchmesser D_M - nun allerdings um -_z zentriert - periodisch wiederholt werden, wenn für die Periodenlänge w gilt w2 D_M=2 Na, und wenn L_C(p) für |p|<D_M auf 0 gesetzt wird

Soll das Bild in einem um --_z zentrierten Quadrat der Seitenlänge D_B · D_B, mit D_BD_M, rekonstruiert werden, so können die f′_c(p_q, ϑ_l) durch Multiplikation mit einem Rechteckfenster R_l(p_q) auf dieses Quadrat begrenzt werden, denn Werte außerhalb tragen nicht zum Bild bei. Es wird ein Rechteckfenster R_l(p_q) mit vom Projektionswinkel ϑ_l abhängiger Ausdehnung ϑ_l gewählt. Das ist möglich, solange gilt

Insbesondere kann man setzen

ϑ_l = 2 D_Bcos ϑ_l -45°ϑ_l45°

ϑ_l = 2 D_Bsin ϑ_l 45°ϑ_l135° (13)

Die mit dem Rechteckfenster R_l(p_q) auf ϑ_l begrenzten und mit L_S gefalteten f′_c(p_q, ϑ_l) können mit der Periodenlänge ϑ_l periodisch wiederholt werden, ohne daß das Bild im Quadrat D_B · D_B beeinflußt wird. Dabei entsteht

Die Fouriertransformierte

lautet

Wie beim Linogramm-Verfahren zerlegt man B_LG(-_z) in ein Teilbild B_LG,1(-_z) mit -45°ϑ_l45° und in ein Teilbild B_LG,2 (-_-z) mit 45°<ϑ_l<135°.

Alle folgenden Betrachtungen beziehen sich auf B_LG,1(-_z), für B_LG,2 (-_z)sind die Ergebnisse ähnlich.

Mit _z · _l = x_z cos ϑ_l+y_z sin ϑ_l und

gilt

_C(mΔρ) ist die an diskreten Frequenzstützstellen vorliegende Fouriertransformierte von L_C.

Ebenso kann man mit · _l = x cos ϑ_l+y sin ϑ_l und

schreiben

_S(sΔρ^l _F) ist die an diskreten Frequenzstützstellen vorliegende Fouriertransformierte von L_U.

Damit folgt für

mit

kann wie folgt veranschaulicht werden:

ist die eindimensionale Fouriertransformierte der f_lk auf das Raster

Im Frequenzraum wird

mit der Foruiertransformierten des Faltungskerns

und dem Phasenfaktor

der die Lage des Rekonstruktionszentrums regelt, multipliziert. Dann schließt sich eine FFT zurück in den Ortsraum an, dabei entstehen die gefalteten, um -_z zentrierten Projektionen f′_c(p_q,ϑ_l)

Im Ortsraum werden die f′_c(p_q,ϑ_l) durch Multiplikation mit R_l(p_q) auf den gewünschten Bildfeldbereich ϑ_l=2 D_B cos ϑ_l, d. h. auf

Werte begrenzt und mit einer chirp-z-Transformation erneut in den Frequenzraum auf das Raster

transformiert

Das Spektrum

ist in ρ_x-Richtung an diskreten Rasterpunkten sΔρ_x ^LG definiert, mit

gehört zu einem in x-Richtung mit einem idealen Rechteckfenster der Breite 2 D_B begrenzten Bild. Die Rasterwerte sΔρ_x ^LG können also direkt als Ausgangspunkte für die FFT in x-Richtung zurück in den Ortsraum verwendet werden, ohne daß eine zusätzliche Interpolation im Frequenzraum nötig ist. Dabei wird - wie beim Gridding-Verfahren - ein Gebiet der Seitenlänge D_R=2 D_B rekonstruiert. Wird das Bild D_B · D_B mit 512 · 512 Pixeln dargestellt, sind in ρ_x-Richtung 1024 Frequenzwerte notwendig. Sind nach der chirp-z- Transformation der Projektionen mehr Frequenzstützpunkte in ρ_x-Richtung vorhanden, weil die in den gefalteten Projektionen enthaltene Maximalfrequenz größer ist als die mit der Bildmatrix darstellbare Frequenz, müssen die zusätzlichen Frequenzpunkte nach dem Gridding in ρ_y-Richtung (siehe unten) als Aliasingkomponenten umgeordnet werden.

In ρ_y-Richtung sind die Rasterwerte ebenfalls äquidistant, das von sΔ_x ^LG abhängige Raster ist

Wie beim originalen Linogramm-Verfahren könnte sich nun in ρ_y-Richtung für jedes s eine chirp-z-Transformation auf das gewünschte Pixelraster Δ_y anschließen.

Um die aufwendige chirp-z-Transformation mit von sΔρ_x ^LG abhängigem Raster Δρ_y ^LG(s) zu umgehen, ist es wünschenswert, in jeder Spalte s der Frequenzmatrix eine eindimensionale Interpolation so durchzuführen, daß ein für alle s gleiches Frequenzraster Δρ_y entsteht, also eine kartesische Matrix. Soll durch diese Interpolation keine Bildverschlechterung eintreten, muß sie als eindimensionales "Gridding" realisiert werden. Die Grundidee für die Wahl der Interpolationsfunktion (ρ_y) besteht darin, das Bild B_LG,1(x, y) gedanklich durch Multiplikation mit einer Fensterfunktion T(y) in y-Richtung auf die gewünschte Breite D_R zu begrenzen, so daß bei der durch die Spektrumsabtastung in kartesischen Koordinaten hervorgerufenen periodischen Wiederholung des Bildes im Ortsraum keine Überdeckungsfehler entstehen. Durchgeführt wird die Begrenzung nicht im Ortsraum, sondern im Frequenzraum als eindimensionale Faltung des an diskreten Punkten

definierten Bildspektrums in ρ_y-Richtung mit der Fouriertransformierten (ρ_y) von T(y). Die Adreßrechnungen dürften sich dabei wegen der für jedes s äquidistanten Ausgangspunkte lΔρ_y ^LG(s) vereinfachen. Das so entstehende, in ρ_y-Richtung kontinuierliche Spektrum wird dann an den Rasterpunkten jΔρ_y, mit Δρ_y=1/D_R abgetastet. Um eine möglichst geringe Ausdehnung von (ρ) im Frequenzraum und damit eine möglichst kurze Faltung zu erhalten, wählt man die Seitenlänge des Rekonstruktionsfeldes D_R doppelt so groß wie die Seitenlänge des gewünschten Bildfeldes D_B. T(y) kann dann im Bereich D_B/2|x|D_R-D_B/2 sanft auf einen hinreichend kleinen Wert abfallen.

In ρ_x-Richtung ist kein Gridding nötig, denn hier liegt das kartesische Raster ja bereits vor. Das aus diesem Spektrum durch kartesische Abtastung und zweidimensionale Fourierrücktransformation entstehende Bild ist mit dem Faltungs-Rückprojektionsbild im gewünschten Bildbereich identisch bis auf Aliasingfehler im Ortsraum, die durch die endliche Ausdehnung von im Frequenzraum hervorgerufen werden. Durch geeignete Wahl von kann man erreichen, daß diese Aliasingfehler vernachlässigbar klein sind.

Durch die Faltung mit (ρ_y) im Frequenzraum ist das Bild im Ortsraum mit T(y) multipliziert, es hat also eine "Beule". Diese Beule kann leicht korrigiert werden, indem man das rekonstruierte Bild B_LG,1(x, y) im Bereich |y|D_B/2 durch T(y) teilt.

Nach dem Gridding und einer eventuellen Umordnung der Aliasingkomponenten hat man sowohl in ρ_x als auch in ρ_y-Richtung 1024 kartesische Werte, so daß man das Teilbild B_LG,1(-_z) durch zweidimensionale Rücktransformation von 1024 · 1024 Frequenzpunkten gewinnt. Die inneren 512 · 512 Pixel sind das Bild, das in y-Richtung durch Multiplikation mit 1/T(y) entzerrt werden muß.

Auf die gleiche Art - nur mit vertauschtem x und y - erhält man das Teilbild B_LG,2(-_z) aus

mit

Da dieses Teilbild nach dem Gridding und der Rücktransformation in x-Richtung entzerrt werden muß, können nicht einfach die Spektren von B_LG,1(-_z) und B_LG,2(-_z) addiert und gemeinsam zurücktransformiert werden. Vielmehr muß jedes Spektrum für sich in den Ortsraum transformiert werden. Das endgültige Bild entsteht aus der Summe der beiden entzerrten Teilbilder.

Zusammenfassung der Bearbeitungsschritte

Im folgenden werden die notwendigen Bearbeitungsschritte dargestellt für ein CT-Gerät der 3. Generation mit 1024 Meßwerten je Projektion und ≈1200 Fächerprojektionen im Bereich [0,2 π] (z. B. Siemens SOMATOM AR/HP).

• 1. Rebinning auf Paralleldaten f_lk:
  N = 1024 Kanäle je Projektion, N_P≈600 Projektionen mit Δ tanϑ=const im Winkelbereich [-45, 45°] bzw. Δ cot ϑ=const im Winkelbereich [45, 135°].
• 2. Zeropadding der f_lk und eindimensionale FFT der f_lk mit 2048 Werten. Durch das zeropadding entsteht im Frequenzraum das Raster
• 3. Multiplikation der _l(mΔρ) mit der Fouriertransformierten des Faltungskerns _C(mΔρ) und mit dem Phasenfaktor exp (-2 π im Δρ(x_z cos ϑ_l+y_z sin ϑ_l)).
• 4. FFT der Länge 2048 zurück in den Ortsraum, Begrenzung jeder gefalteten, zentrierten Projektion f′_c(p_q, ϑ_l) auf Werte für -45°ϑ_l45° bzw. auf Werte für 45°<ϑ_l135°.
• 5. Eindimensionale chrip-z-Transformation der begrenzten f′_c(p_q, ϑ_l) auf maximal 1024 Frequenzpunkte im Raster Die Zahl der Punkte ist abhängig von der Größe des dargestellten Bildausschnittes und von der Maximalfrequenz der Rekonstruktion.
• 6. Multiplikation der _l(sΔρ^l _F) mit der Fouriertransformierten des "Glättungskerns" _S(sΔρ^l _F).
• 7. Berechnung zweier Teilbilder für die Winkelbereiche [-45°, 45°] und [45°, 135°]: Teilbild 1: für jeden Frequenzwert sΔρ_x in ρ_x-Richtung führt man ein Gridding in ρ_y-Richtung durch, ausgehend von N_P Frequenzstützstellen im s-abhängigen Raster Δρ_y ^LG(s) auf 1024 Werte im Raster Da die Matrix komplex konjugiert ist, müssen nur die positiven s berücksichtigt werden, man hat also höchstens ≈1024 · 600 Ausgangspunkte. Nach der Umordnung der Aliasingkomponenten schließt sich eine Fouriertransformation der 1024 · 1024 Frequenzwerte in den Ortsraum an, die inneren 512 · 512 Pixel sind das Bild.Teilbild 2: Vorgehensweise wie für Teilbild 1, mit vertauschtem x und y.
• 8. Entzerrung und Addition der beiden Teilbilder.

Claims (1)

1. Verfahren zur Bildrekonstruktion für einen Computertomographen, bei dem die vom Detektor gelieferten Daten f(p_k, ϑ_l), aus denen ein Rechner ein Schnittbild berechnet, in Parallelstrahlgeometrie so vorliegen, daß Δ tan ϑ=const im Winkelbereich [-45°, 45°] bzw. Δ cot ϑ=const im Winkelbereich [45°, 135°], und p_k=(k-0,5+a_m)a (a ist das Abtastraster), wobei folgendes Rekonstruktionsverfahren angewendet wird:
   Eindimensionale FFT der f(p_k, ϑ_l) bezüglich p_k mit N_FFT2 N Werten (N ist die Kanalzahl), das Ergebnis ist _l(mΔρ), Multiplikation der _l(mΔρ) mit der Fouriertransformierten des Faltungskerns _C(mΔρ) und mit dem Phasenfaktor exp(-2 π im Δρ(x_z cos ϑ_l+ y_z sin ϑ_l)), der die Lage des Rekonstruktionszentrums regelt, FFT der Länge N_FFT zurück in den Ortsraum, Begrenzung jeder gefalteten, zentrierten Projektion f′_c(p_q, ϑ_l) auf Werte für -45°ϑ_l45° bzw. auf Werte für 45°<ϑ_l135°, eindimensionale chirp-z-Transformation der begrenzten f′_c(p_q,ϑ_l) auf das Raster für -45°ϑ_l45° bzw. für 45°<ϑ_l135°, wobei die Zahl der Punkte abhängig ist von der Größe des dargestellten Bildausschnittes und von der Maximalfrequenz der Rekonstruktion, Multiplikation der _l(sΔρ^l _F) mit der Fouriertransformierten des "Glättungskerns" _S(sΔρ^l _F), Berechnung zweier Teilbilder für die Winkelbereiche [-45°, 45°] und [45°, 135°]: Teilbild 1: für jeden Frequenzwert sΔρ_x in ρ_x-Richtung führt man ein eindimensionales Gridding, eine Interpolation mit (ρ_y) in ρ_y-Richtung durch, ausgehend von N_P Frequenzstützstellen im s-abhängigen Raster Δρ_y ^LG(s) auf 2 N_Pix Werte im Raster (N_Pix ist die Pixelzahl des CT-Bildes in einer Spalte), nach der Umordnung der Aliasingkomponenten schließt sich eine Fouriertransformation der 2 N_Pix · 2 N_Pix Frequenzwerte in den Ortsraum an, die inneren N_Pix · N_Pix Pixel sind das Bild, das in y-Richtung durch Multiplikation mit 1/T(y) entzerrt werden muß, sowie Teilbild 2: Vorgehensweise wie für Teilbild 1, mit vertauschtem x und y, und Addition der beiden Teilbilder.