DE19802850A1

DE19802850A1 - Bildrekonstruktionsverfahren für die 3D-Rekonstruktion

Info

Publication number: DE19802850A1
Application number: DE19802850A
Authority: DE
Inventors: Karl Dr Ing Barth; Karl Wiesent
Original assignee: Siemens AG
Current assignee: Siemens AG
Priority date: 1998-01-26
Filing date: 1998-01-26
Publication date: 1999-07-29
Also published as: JPH11267117A; JP4350190B2; US6125163A

Description

Die Computertomographie ist ein mittlerweile weit verbreite tes Verfahren in der Medizintechnik und zur Materialprüfung. In einer zweidimensionalen Variante wird z. B. ausgehend von einer Röntgenquelle von verschiedenen Positionen aus eine Schicht des Objekts durchstrahlt. Die geschwächten Intensitä ten werden dann hinter dem Objekt von einem eindimensionalen Detektorarray gemessen und in digitaler Form einem Rechner zugeführt. Die Verteilung der Schwächungswerte in der durch strahlten Objektschicht wird dann mittels eines mathemati schen Rekonstruktionsverfahrens im Rechner ermittelt. Die am weitesten verbreitete Methode ist die gefilterte Rückprojek tion. In ihrer allgemeinen Form kann sie durch folgende Schritte charakterisiert werden:

1. Vorverarbeitung (Gewichtung) der Meßwerte
2. Faltung der Meßwerte (=Filterung)
3. Gewichtung der gefalteten Daten
4. Rückprojektion
5. Skalierung des Ergebnisses.

Rechenintensiv und damit zeitkritisch sind dabei vor allem die Faltung und die Rückprojektion.

Verzichtet man auf die Faltung, so entsteht ein sehr ver waschenes Bild, das sogenannte Layergram, das keine Details erkennen läßt.

Allgemein wird die diskrete Faltung durch folgende Formel be schrieben:

Dabei stellt x(n) die Reihe der Meßwerte dar, h(n) ist der sogenannte Faltungskern und das Ergebnis y(n) die Reihe der gefalteten (gefilterten) Daten.

Ein Wert y(n) des Faltungsergebnisses ist damit eine gewich tete Summe aller Eingangswerte x(n). Da ein Meßwertsatz x(n) in der Praxis nur aus endlich vielen, allgemein N, Werten be steht, reduziert sich obige Formel auf eine endliche Summe. Um einen einzelnen gefalteten Wert zu berechnen, sind also im allgemeinen N Multiplikationen und N Additionen nötig. Sind beim verwendeten Faltungskern h(n) nur M Komponenten von Null verschieden und ist M < N, so führt dies zu einer weiteren Reduktion des Rechenaufwandes. Für die Berechnung eines ge filterten Wertes y(n) sind dann nur maximal M Multiplika tionen und Additionen nötig.

Die in der Computertomographie üblicherweise verwendeten Kerne haben volle Länge, d. h. M = 2N-1. Die Anzahl der arithme tischen Operationen wird also durch N, die Länge der Meßwert felder x(n), bestimmt. Es handelt sich bei diesen Kernen um Varianten, die vom sog. Rampenfilter abgeleitet werden. Die beiden bekanntesten sind der Kern von Ramachandran und Lakshminarayanan sowie der Kern von Shepp und Logan.

Die verwendeten CT-Kerne haben üblicherweise folgende Eigen schaften:

1. Symmetrie, d. h. h(n) = h(-n). Deshalb genügt es z. B. je weils nur die Zentralkomponente h(0) und die rechte Hälfte h(n), mit n = 1, 2, . . ., anzugeben.
2. Fehlender Gleichanteil, d. h. die Summe über alle Kernkom ponenten ist Null. Dies entspricht dem Fakt, daß im Fou rierraum das Rampenfilter H(s) im Nullpunkt verschwindet.
3. Die Zentralkomponente h(0) ist die einzige positive Kom ponente, alle anderen sind negativ oder Null. Wegen der abschließenden Skalierung des gesamten Rekonstruktionser gebnisses kann der Faltungskern ohne Beschränkung der Allgemeinheit so normiert werden, daß die Zentralkompo nente den Wert h(0) = 1 hat. Für n → ∞ geht der Betrag der Kernkomponenten gegen Null.

Die oben beschriebene Faltung im Ortsraum ist äquivalent zu einer komponentenweisen Multiplikation im Fourier-Raum. Ab einer gewissen Kernlänge und bei Verwendung von Kernen voller Länge kann es numerisch effektiver sein, die Filterung multi plikativ im Fourier-Raum durchzuführen, d. h. zweimal die schnelle Fouriertransformation (FFT) und dazwischen eine numerisch sehr effiziente multiplikative Filterung anzuwen den:

FFT (hin): x(n) → X(s)
Y(s) = X(s).H(s)
FFT (zurück): Y(s) → y(n).

Ist die Anzahl M der von Null verschiedenen Komponenten h(n) deutlich geringer als die Länge N des Meßfeldes, so führt dies nur bei der Ortsraumvariante zu einer Reduktion des Rechenaufwandes, nicht jedoch bei der Fouriermethode.

Ein Trend in der Medizintechnik besteht in der Vermessung und Rekonstruktion ganzer Volumina in möglichst kurzer Zeit. Dies geschieht durch den Einsatz von Flächendetektoren statt eines eindimensionalen Detektorarrays. Möglich ist auch eine Volu menrekonstruktion aus den Meßdaten einer Rotationsangiogra phie. Fig. 1 zeigt den prinzipiellen Aufbau eines solchen Ge rätes. Solche Geräte nutzen die Röntgenstrahlung wesentlich besser als Schnittbildtomographen, bei denen nur ein ebener Fächer ausgeblendet wird. Als Ausführungsbeispiel für diese Erfindung diene eine solche Angiographieanlage (C-Bogen-Ge rät).

Die Röntgenquelle 1 läuft dabei auf einem Kreisbogen 2 um das Objekt 3. Durchstrahlt wird dabei jeweils ein ganzes Teil volumen des Objektes 3. Die geschwächten Intensitäten hinter dem Objekt 3 werden von einem Flächendetektor 4, meist einem Röntgenbildverstärker (RBV), gemessen. Für diese Art von Geometrie hat Feldkamp, LA, David LC and Kress JW: Practical cone-beam algorithm. J. Opt. Soc. Am., 1: 1612-619, 1984 einen effizienten approximativen 3D-Rekonstruktionsalgorithmus vor geschlagen. Dieser Algorithmus ist auch vom Typ gefilterte Rückprojektion und ähnelt stark dem oben beschriebenen 2D-Algorithmus. Insbesondere wird beim Feldkamp-Algorithmus die Filterung der Meßwerte des Flächendetektors nicht etwa zwei dimensional vorgenommen, sondern nur zeilenweise, also ein dimensional wie bei der klassischen Schnittbildtomographie mit einem 1-Zeilen Detektor. Alles was im ersten Abschnitt zur Faltung gesagt wurde, bleibt unverändert gültig. Von dem von der Röntgenquelle 1 ausgehenden, pyramidenförmigen Rönt genstrahlenbündel ist nur der Fächer 5 und auf dem Flächen detektor 4 nur die Detektorzeile 6 aus einer Reihe von Detek torelementen gezeigt.

Typische Daten für ein solches C-Bogen-Gerät sind ca. 50 Pro jektionen mit jeweils 1024×1024 Meßwerten. Beim Feldkamp-Algorithmus fallen damit über 50000 eindimensionale Faltungen über Meßwertfelder x(n) der jeweiligen Länge N = 1024 an. Bei Verwendung normaler CT-Kerne voller Länge, wie etwa dem von Shepp und Logan, ist dies ein enormer Rechenaufwand, der für manche Anwendungen, z. B. Einsatz während der Intervention, zu prohibitiv langen Rechenzeiten führt.

In vielen praktischen Fällen ist das Objekt 3 breiter als der Detektor, d. h. das zu rekonstruierende VOI (= Volume of Interest) ist zwar jeweils im Strahlengang, nicht jedoch der gesamte Objektquerschnitt. Die Fig. 2 beschreibt diesen Sach verhalt. Hier ist das Meßfeld mit 7 bezeichnet. Die Meßwert profile haben dann nicht mehr die Gestalt wie in der Fig. 3 mit einer Abnahme gegen Null am Rande, sondern sie weisen Sprungstellen auf, wie in der Fig. 4 gezeigt. Diese unnatür lichen Sprungstellen sind Fehlerstellen in den Meßdaten x(n), die sich bei Verwendung ausgedehnter Kerne auf alle Daten y(n) fortpflanzen. Theoretisch ist aus solchen abgeschnitte nen Projektionen eine exakte Rekonstruktion nicht möglich. In der Praxis werden aber dennoch recht brauchbare Ergebnisse erzielt. Im wesentlichen sind zwei Methoden bekannt, dieses Problem anzugehen:

Extrapolation der Meßwerte
Übergang zu lokalen Faltungskernen.

Extrembeispiel eines lokalen Kernes ist der sog. Laplace-Kern, der einer zweifachen Differentiation entspricht:

h(0) = 1
h(1) = h(-1) = -1/2
h(n) = 0 sonst.

Verwendet man diesen Kern, so wird nicht das ursprüngliche Objekt f(x) rekonstruiert, sondern eine modifizierte Vertei lung λf(x). Diese weist die gleichen Kanten auf, wie f(x), sie sind sogar noch überhöht, informiert aber nicht mehr über den wahren Wert in homogenen Bereichen. Ist das Objekt z. B. ein homogener Zylinder, so hängen die rekonstruierten Werte in der Mitte durch, wie in der Fig. 5a skizziert. Im Ver gleich zur normalen CT-Rekonstruktion, z. B. mit dem Shepp und Logan Kern, wird durch den kurzen Laplace-Kern verstärkt Rau schen in das Bildzentrum übertragen. Um diesen Effekt zu mil dern, verwendet man auch geglättete Varianten dieses diffe renzierenden Kernes, z. B.:

. . ., 0, -1/4, -1/4,1, -1/4, -1/4,0, . . .
oder
. . ., 0, -1/6, -1/3,1, -1/3, -1/6,0, . . .

Verzichtet man auf die Faltung, so entsteht das bereits er wähnte Layergram mit einem Profil wie in der Fig. 5b. Aus mathematischen Gründen wird das Layergram auch mit λ^-1f(x) bezeichnet. Das Weglassen der Faltung ist äquivalent zu einer Faltung mit dem sog. Einheitskern, der definiert ist durch h(0) = 1 und h(n) = 0 sonst. Der Einheitskern hat nicht die typi sche Eigenschaft der CT-Kerne, daß die Summe über alle Kompo nenten Null ist. In der Fig. 5c ist das exakt rekonstruierte Profil eines homogenen Zylinders dargestellt.

Als allgemeine λ-Rekonstruktion wird eine Linearkombination aus diesen beiden Varianten bezeichnet. Die mathematische Theorie geht zurück auf Kennan T. Smith and F. Keinert, Mathematical Foundations of Computed Tomography, Appl. Optics, Vol. 24, No. 23, 1. Dec. 1985, Seiten 3950-3957. Die Grundidee ist, daß sich in der gewichteten Überlagerung in ursprünglich homogenen Objektbereichen das konkave und das konvexe Verhalten der beiden Einzelvarianten zumindest teil weise kompensieren und die Randinformation trotzdem erhalten bleibt.

Eine weitere Möglichkeit ist, einen Standard CT-Kern voller Länge, z. B. den von Shepp und Logan, symmetrisch zum Null punkt nach jeweils L < N Komponenten abzuschneiden, und wahl weise die verbleibenden Komponenten so zu modifizieren, daß der Übergang auf Null glatt erfolgt.

Es sei hier ausdrücklich erwähnt, daß das Ziel einer Rekon struktion nicht die möglichst genaue Berechnung von Röntgen schwächungskoeffizienten ist, sondern die Darstellung signi fikanter Information für den Arzt. Deshalb kommen in der medizinischen Computertomographie je nach Fragestellung ver schiedene Faltungskerne zum Einsatz. Bei der 3D-Rekonstruk tion aus den Daten einer Rotationsangiographie geht es um die Darstellung feiner hochkontrastiger Gefäßbäume und Gefäßano malien. Die Fig. 6 verdeutlicht das Spannungsfeld für eine Kernoptimierung (Bereich der Faltungsoptimierung). Zielgrößen sind:

1. Ergebnisqualität, und nachgeordnet
2. Numerische Effizienz.

Kurze Kerne der Länge M « N haben folgende Eigenschaften:

1. Sie sind in ihrer Wirkung lokal, d. h. etwaige Störstellen wie bei abgeschnittenen Projektionen bleiben in ihrer Wir kung lokal.
2. Sie sind numerisch effizient.

Der Erfindung liegt die Aufgabe zugrunde, ein Bildrekonstruk tionsverfahren für die Computertomographie zu optimieren, so daß sich eine Beschleunigung der Faltung und damit eine schnelle Bilderzeugung bei einem C-Bogengerät mit Flächen detektor ergibt.

Diese Aufgabe ist erfindungsgemäß gelöst durch die Merkmale des Patentanspruchs 1.

Die Erfindung ist nachfolgend anhand eines in der Fig. 7 dar gestellten Ausführungsbeispiels näher erläutert.

In der Fig. 7 ist ein Lagerungstisch 8 dargestellt, auf dem ein Patient 9 liegt und der mit Hilfe eines Stativs 10 an der Decke 11 des Untersuchungsraumes aufgehängt ist. Der Lage rungstisch 8 kann auch unabhängig vom eigentlichen C-Bogenge rät montiert sein, z. B. auf dem Boden des Untersuchungsrau mes. Zur Anfertigung von Röntgenbildern sind ein Röntgen strahler 12 und der Flächendetektor 4 vorgesehen. Der Rönt genstrahler 12 und der Flächendetektor 4 sind an einem C-Bogen 14 befestigt, welcher an einem Sockel 15 verstellbar gelagert ist. Die Bildwiedergabe erfolgt auf einer Monitor ampel 16. Volumendaten für die Erzeugung von dreidimensiona len Bildern können gewonnen werden, wenn der Röntgenstrahler 12 und der Flächendetektor 4 um die Systemachse 17 gedreht werden.

Der Flächendetektor 4 besteht beispielsweise aus einer Matrix von Detektorelementen, die der jeweils empfangenen Strah lungsintensität entsprechende Ausgangssignale digitalisiert einem Rechner 18 zuführen, der ein ausgewähltes Volumen rekonstruiert.

Vorgeschlagen wird eine CT-Rekonstruktion mit einem Filter, dessen Koeffizienten außerhalb der Zentralkomponente sich dem Betrage nach wie eine Exponentialfunktion verhalten.

Beispiel ist ein Filter h(k) mit folgenden Koeffizienten:

h(0) = 1
h(1) = -(1-a)/2
h(k) = h(1).a^k-1
für k < 1, und
h(k) = h(-k) für k < 0.

Der Parameter a kann dabei frei aus dem Intervall [0,1] ge wählt werden. Er dient zur Optimierung des Ergebnisses. Die ses Filter hat für a < 1 alle die im 1. Abschnitt erwähnten Eigenschaften eines CT-Kernes. Insbesondere ist die Summe aller Koeffizienten Null.

Für a = 0 ergibt sich das Laplace Filter . . .0,-1/2,1,- 1/2,0, . . ., und damit die reine X-Rekonstruktion.

Für den Grenzfall a = 1 bekommt man den Einheitskern . . .0,1,0,. . ., und damit das Layergram.

Das Filter hat zwar volle Länge, fällt aber exponentiell und damit schließlich sehr schnell gegen Null. Das Standardfilter von Shepp und Logan (siehe z. B. Heinz Morneburg (ed), Bild gebende Systeme für die medizinische Diagnostik, Publicis MCD Verlag, 1995) fällt hingegen nur mit 1/k². Das Exponential filter hat damit einen gewissen lokalen Charakter. Durch Variation von a kann gesteuert werden, wie schnell die Werte im zentralen Teil fallen sollen, z. B. schneller als Shepp-Logan, oder erst langsamer. In der Tabelle gemäß Fig. 8a, b sind die jeweils ersten zehn Werte für verschiedene Belegun gen von a aufgelistet. Die erste Zahl einer jeden Zeile ist der laufende Index gefolgt von den Beträgen der jeweiligen Komponente des Shepp-Logan Kernes und des zum Parameter a ge hörenden jeweiligen Exponentialfilters. Die beiden letzten Spalten enthalten die Partialsummen zu den beiden Kernen, d. h. die Summe der jeweiligen Komponenten vom ersten bis zum aktuellen Index. Wegen der Normierung der Zentralkomponente auf eins und der Symmetrie ist die Konvergenz der Partialsum men gegen 0.5 äquivalent zur Nullsummeneigenschaft.

Bei Rekonstruktion von Gefäßbäumen aus mit Rotationsangiogra phie gemessenen Patientendaten kann mit solchen Exponential kernen eine Bildqualität erreicht werden, die der eines Shepp und Logan Kernes voller Länge entspricht bzw. diese sogar übertrifft.

Rekursive Filter sind in der Nachrichtentechnik wohl bekannt. Sie sind von Natur aus in ihrer Wirkung stark unsymmetrisch. Eine Symmetrisierung kann aber erreicht werden, indem man die Filterung von links nach rechts (steigend) und von rechts nach links (fallend) vornimmt und die beiden Ergebnisse mit telt.

Im folgenden soll gezeigt werden, daß das Ergebnis einer Fal tung mit dem oben definierten Exponentialfilter auch erreicht wird, wenn man von den Originaldaten das Ergebnis einer ge mittelten rekursiven Filterung der Ordnung eins abzieht.

Ein rekursives Filter der Ordnung eins ist definiert durch die Vorschrift

y(n) = a.y(n-1) + b.x(n).

Dieses kann umgeschrieben werden zu

Schreibt man u(n) für die aufsteigende Richtung und v(n) für das Ergebnis in fallender Richtung, so erhält man

Dies ist identisch zum Ergebnis einer normalen Faltung mit dem symmetrischen Kern h, definiert durch h(0) = 2.(1-b) und h(i) = -b.aⁱ.

Normiert man die Zentralkomponente auf eins und setzt b = 1-a so hat man die oben beschriebene Filterklasse.

Die rekursive Filterung ist numerisch sehr effizient. Um z. B. einen Wert u(n) zu berechnen, ist nur eine Multiplikation und eine Addition nötig. Würde man die Faltung mit dem ausgedehn ten Exponentialkern wie üblich durchführen, so wären zur Be rechnung eines Wertes y(n), wie in Abschnitt 1 beschrieben, N Multiplikationen und Additionen nötig, bei der Realisierung über die rekursive Filterung aber nur zwei Multiplikationen und drei Additionen.

Die Faltung mit dem ausgedehnten Exponentialkern kann also mit der gleichen Effizienz implementiert werden, wie die Fal tung mit dem kurzen Laplace Kern.

Bei dem Ausführungsbeispiel Rotationsangiographie ist N = 1024 und mit der neuen Methode ergibt sich größenordnungsmäßig eine hundertfache Beschleunigung der Faltung.

Das erfindungsgemäße Verfahren ist bei dem beschriebenen Aus führungsbeispiel im Rechner 18 durch Software realisiert. Dies gilt insbesondere für das rekursive Filter. Dieses kann jedoch auch hardwaremäßig realisiert werden.

Die Datengewinnung ist in Verbindung mit Fig. 7, d. h. in Ver bindung mit einem C-Bogengerät beschrieben. Auch ein üblicher Computertomograph mit einem Tragring für den Röntgenstrahler und den Flächendetektor ist anwendbar.

Claims

1. Bildrekonstruktionsverfahren für die 3D-Rekonstruktion aus den Detektordaten eines Röntgengeräts, insbesondere eines Computertomographen, mit einem Röntgenstrahler (12) und einem Flächendetektor (4), die auf einer verstellbaren Haltevor richtung (14, 15) zur Durchstrahlung eines Meßfelds (7) mit Hilfe eines pyramidenförmigen oder kegelförmigen Röntgen strahlenbündels unter verschiedenen Richtungen gelagert sind, bei dem ein Rechner (18) zur Bildrekonstruktion aus den vom Flächendetektor (4) gelieferten Signalen vorgesehen ist, wo bei die Rekonstruktion eines dreidimensionalen Teilvolumens des Objekts (3) durch gefilterte Rückprojektion erfolgt und die Komponenten des Faltungskerns außerhalb der Null-Komponente einer Exponentialfunktion entsprechen.

2. Verfahren nach Anspruch 1, bei dem die Filterung durch ein rekursives Filter realisiert ist.

3. Verfahren nach Anspruch 1 oder 2, bei dem die Filterung dadurch realisiert ist, daß die Eingangsdaten von beiden Sei ten her einer rekursiven Filterung der Ordnung 1 unterzogen werden und die Summe der Ergebnisse von den verdoppelten Ori ginaldaten abgezogen wird.

4. Verfahren nach einem der Ansprüche 1 bis 3, bei dem zur Filterung ein rekursives Filter beliebiger Ordnung verwendet wird.