DE19626775A1

DE19626775A1 - Schnelle Faltung von Projektionen

Info

Publication number: DE19626775A1
Application number: DE19626775A
Authority: DE
Inventors: Karl Dr Ing Schwarz
Original assignee: Siemens AG
Current assignee: Siemens AG
Priority date: 1996-07-03
Filing date: 1996-07-03
Publication date: 1998-01-08
Also published as: US5943434A

Description

Die Erfindung betrifft ein Verfahren für die Bildrekonstruk tion in der bildgebenden Technik zur Durchführung einer schnellen Faltung mit der Transformationslänge M unter Zulas sung geringer Überfaltungsfehler, wobei gemessene Projektio nen der Länge N mit einem modifizierten Filterkern unter Ver wendung der "Fast Fourier Transform" (FFT) und der "Inverse Fast Fourier Transform" (IFFT) gefaltet werden.

In der bildgebenden Technik, z. B. der Computertomographie, werden Projektionen eines Patienten, welche beispielsweise von einem um den Patienten umlaufenden radiologischen Meßsy stem ermittelt werden, zur Bildrekonstruktion mehreren mathe matischen Operationen unterzogen, wobei diese u. a. mit einem Filterkern gefaltet werden. Durch die Faltung, welche im we sentlichen einer Hochpaßfilterung entspricht, werden weitrei chende Verschmierungen einzelner Objektdetails im rekonstru ierten Bild vermieden, die sich bei einer unmittelbaren Rück projektion gemessener Schwächungsprofile (Projektionen) er geben würden. Die Anwendung der Faltung ist dabei wesentliche Voraussetzung für die Ermittlung eines Sofortbildes, d. h. die Bereitstellung des berechneten Bildes unmittelbar nach Beendigung des Meßvorganges. Da die Rekonstruktion von Bil dern folglich ein zeitkritischer Prozeß ist, wird zur Berech nung der Faltung, d. h. zur schnellen Ausführung der Faltung in einem Rechner, die sogenannte "schnelle Faltung" einge setzt, welche in der Regel zyklisch und diskret berechnet wird.

Ist ein aus einer Projektion gebildeter mit einem Filterkern zu faltender Vektor N Elemente lang, so muß die Transforma tionslänge M nach der Theorie der schnellen Faltung M 2N- 1 sein, um fehlerfreie Ergebnisse aus der schnellen Faltung zu erhalten. Da man aber bei der Durchführung der schnellen Faltung unter Verwendung der Fast Fourier Transformation und der Inversen Fast Fourier Transformation aus Vektorbiblio theken meist nur Zweierpotenzen mit der Transformationslänge M = 2^m zur Verfügung hat, wählt man als Transformationslänge M folglich die kleinste Zweierpotenz, welche die Bedingung M = 2m 2N-1 erfüllt.

Die Wahl von M ist einfach, wenn N eine Zweierpotenz S^m-1 ist. Ist die Länge N des zu faltenden Vektors aber keine Zweierpo tenz (z. B. N = 1536), muß die Transformationslänge M auf die nächste Zweierpotenz aufgerundet werden (z. B. hier 4096). Dabei entsteht ein unnötiger Overhead von M-(2N-1) (z. B. hier 1025) Werten, welche in der schnellen Faltung verarbei tet werden müssen und zu einer unerwünschten hohen Rechenzeit der schnellen Faltung führen können.

Um dieser Rechenzeiterhöhung entgegenzuwirken, kann man unter Zulassung von Überfaltungsfehlern ein in Fig. 1 veranschau lichtes Verfahren zur schnellen "Faltung mit der nächsten Zweierpotenz" verwenden. Ein aus einer Projektion der Länge N (2^m-1 N < 2^m, z. B. N ≈ 0,6 _*2^m) gebildeter Vektor wird da bei durch Anfügen von S Nullen an die Projektion der Länge N auf die Transformationslänge M = N+S = 2^m erweitert. Der zur schnellen Faltung verwendete Filterkern h(k) im Ortsbereich (= Impulsantwort) wird für eine Transformationslänge von M berechnet, wobei die Impulsantwort zur Bildverbesserung auch um wenige Werte kürzer als M gewählt werden kann. Nach der schnellen Faltung mit der Transformationslänge M des Vektors mit dem etwas zu kurz gewählten Filterkern h(k) sind im ge filterten Vektor der Länge N die S+1 mittleren Werte richtig, während am rechten und linken Rand der Projektion insgesamt N-S-1 Werte verfälscht sind. Dabei entsteht zum einen auf grund der zu kurzen Transformationslänge M des Filterkerns h(k) ein Überfaltungsfehler und zum anderen wird ein Teil der Projektion der Länge N mit falschen Filterkoeffizienten gewichtet. Da aber einerseits diese falschen Filterkoeffi zienten am äußeren Rand des Filterkerns liegen und somit sehr klein sind, und andererseits das Meßobjekt am Rand auch ge ringere Schwächungswerte aufweist, ist bei einer hinreichend großen Anzahl von S an eine Projektion der Länge N angefügten Nullen der Fehler vernachlässigbar.

Implementiert man das Verfahren zur schnellen "Faltung mit der nächsten Zweierpotenz " auf einem Rechner, bestehen zwei Möglichkeiten der Durchführung der Faltung: Zum einen die Transformation einer Projektion mit der Transformationslänge M/2 mit Nachverarbeitung, welche die unvollständige Transfor mation korrigiert, und zum anderen die Transformation von zwei Projektionen, je eine im Real- und eine im Imaginärteil einer Transformation der Transformationslänge M.

Ist nun aber die Anzahl S (z. B. S < 0,5 _*N) von an die Pro jektion der Länge N angefügten Nullen zu gering, führt das vorstehend beschriebene Verfahren zu einem zu großen, nicht mehr tolerierbaren Fehler am Außenrand des rekonstruierten Bildes.

Der Erfindung liegt daher die Aufgabe zugrunde, ein Verfahren der eingangs genannten Art derart auszubilden, daß der auf tretende Fehler im rekonstruierten Bild verringert und die Rechenzeit für die schnelle Faltung verkürzt wird.

Nach der Erfindung wird diese Aufgabe gelöst durch ein Ver fahren für die Bildrekonstruktion in der bildgebenden Technik zur Durchführung einer schnellen Faltung mit der Transforma tionslänge M unter Zulassung geringer Überfaltungsfehler, wobei jeweils p (oder 2p) von n gemessenen Projektionen der Länge N mit einem modifizierten Filterkern (k) unter Ver wendung der "Fast Fourier Transform" (FFT) und der "Inverse Fast Fourier Transform" (IFFT) gleichzeitig in einem Schritt gefaltet werden, umfassend folgende Schritte:

a) Zusammenfassung der p Projektionen der Länge N zu einem Vektor und Anfügung von S Nullen hinter jeder der p Pro jektionen der Länge N zu einem Vektor mit der Transformationslänge M,
b) Wahl der Parameter p und S derart, daß die Transformationslänge M = p (N+S) = 2^m (2)zur Durchführung der schnellen Faltung eine Zweierpotenz ist,
c) Entwurf des Filterkerns (k) im Ortsbereich derart, daß er wenigstens im wesentlichen gleich einer Länge N+S ist, und Auffüllen des Filterkerns (k) mit Nullen auf die Transformationslänge M gemäß
d) Zyklische Durchführung der schnellen Faltung mit der Transformationslänge M gemäß mit k = 0(1) (M-1),wobei das Faltungsergebnis einer Projektion i die Formy_i(k) = y(k + (i-1) (N+S)) (5)mit k = 0(1) (N-1) und
i = 1(1)p

aufweist. Durch die gleichzeitige Transformation mehrerer (p) in einem Vektor (k) zusammengefaßter Projektionen der Länge N in einem Schritt der schnellen Faltung mit einer größeren Transformationslänge kann also der Overhead von Werten, wel cher bei der definitionsgemäßen schnellen Faltung auftritt, falls die Länge N der Projektionen keine Zweierpotenz ist, verringert werden. Zudem wird durch die Zulassung geringer Fehler bei der Durchführung der schnellen Faltung nach dem erfindungsgemäßen Verfahren die Effektivität gegenüber der definitionsgemäßen schnellen Faltung weiter erhöht, ohne daß die auftretenden gegenläufigen und sich annähernd kompensie renden Fehler das Ergebnis der Bildrekonstruktion merklich stören. Erfüllt nach dem erfindungsgemäßen Verfahren die An zahl S von an eine Projektion der Länge N angefügten Nullen dabei die Bedingung S N-1 und liegt die Filterlänge des Filterkerns (k) im Ortsbereich zwischen 2N-1 und 2S+1, er hält man sogar keine Überfaltungsfehler und es ergibt sich trotzdem der Vorteil, daß bei bestimmten Längen von M und N ein Zeitvorteil gegenüber der herkömmlichen definitions gemäßen schnellen Faltung erreicht wird (z. B. N = 682, M = 4096, p = 3). Im Falle, daß die Beziehung S < N-1 gilt, entstehen zwar Überfaltungsfehler bei der schnellen Faltung nach dem erfindungsgemäßen Verfahren, die aber bei Anwendung bei spielsweise in der Computertomographie und bei hinreichend großer Wahl von S an eine Projektion der Länge N angefügten Nullen vernachlässigt werden können. Die Filterlänge im Orts bereich wird dabei in der Größenordnung N+S gewählt und mit Nullen auf die Transformationslänge M aufgefüllt. Wird die Filterlänge im Ortsbereich geringfügig kürzer oder länger als N+S gewählt, erhält man zwar je nach Berechnungsfall gering fügig bessere oder schlechtere Ergebnisse bei der Bildrekon struktion, man hat aber immer den Vorteil deutlich schneller als die definitionsgemäße Faltung zu sein.

Nach der Faltung eines Vektors (k) mit dem modifizierten Faltungskern (k) nach dem erfindungsgemäßen Verfahren sind die S+1 mittleren Werte jeder gefilterten Projektion korrekt berechnet. Am rechten und linken Rand jeder Projektion sind dann insgesamt N-S-1 Werte verfälscht. In diesen verfälschten Werten sind zwei Fehlerbeträge enthalten, die in Fig. 2 veran schaulicht werden, welche die fehlerhafte Faltung von p = 4 Projektionen zeigt. Zum einen tritt ein Überfaltungsfehler wegen der zu kurzen Länge S der Nullsequenz auf, so daß ein Teil einer anschließenden Projektion mit Filterkoeffizienten gewichtet und addiert wird (vgl. Fig. 2, schraffierter Be reich, bei Projektion ₄(k)). Zum anderen ergibt sich ein weiterer Fehler dadurch, daß durch die zu kurz gewählte Fil terlänge des Filterkerns (k) im Ortsbereich von N+S ein Teil der Projektion nicht gefiltert wird (vgl. Fig. 2, schraf fierter Bereich, bei Projektion ₁(k)). Beide auftretenden Fehler sind aber gegenläufig, da kurz aufeinanderfolgende Projektionen einen ähnlichen Verlauf bei identischen Kanälen aufweisen. Auch die dort wirksamen Filterkoeffizienten sind nicht sehr unterschiedlich und die Steigung ist umgekehrt. Dadurch ist eine Verbesserung des Fehlerverhaltens zu erzie len. Da weiterhin die dort wirksamen Filterkoeffizienten sehr klein sind, ist eine Auswirkung auf das rekonstruierte Bild mit einem Einfluß von mehr als einer Hounsfield-Einheit (HU) bei hinreichend kleiner Wahl von N-S zu vermeiden. Wenn das Meßfeld einen Radius R hat, so ergibt sich der Radius r, bei dem im Bild kein Unterschied zur exakten Faltung besteht, aus der folgenden Gleichung 6 (vgl. auch Fig. 5):

r = ((S+1)/N) _*R (6)

Implementiert man das erfindungsgemäße Verfahren zur schnel len Faltung auf einem Rechner bestehen, wie im Falle der schnellen "Faltung mit der nächsten Zweierpotenz", die zwei bereits erwähnten Möglichkeiten der Durchführung der Faltung:
Zum einen die Transformation von p Projektionen mit der Transformationslänge M/2 mit Nachverarbeitung und zum anderen die Transformation von 2_*p Projektionen, je p im Real- und p im Imaginärteil einer Transformation der Transformationslänge M. Wenn die Anzahl der Projektionen n dabei nicht durch p teilbar ist, wird bei der letzten Faltung die letzte gültige Projektion wiederholt, um den Vektor zu füllen. An dieser Stelle einen Übergang zu der definitionsgemäßen schnellen Faltung zu machen, würde bedeuten, daß auch wieder der origi nale Filterkern eingesetzt werden müßte. Ein Auffüllen des Vektors bei der letzten Faltung mit Nullen würde dabei nur den Fehler der Überfaltung eliminieren, die absichtlich zu kurze gewählte Filterlänge des Filterkerns würde aber blei ben. Somit könnte der Kompensationseffekt beim Auftreten bei der Fehler nicht genutzt werden.

Ein Ausführungsbeispiel der Erfindung ist in den beigefügten Zeichnungen dargestellt. Es zeigen:

Fig. 1 ein Beispiel einer fehlerbehafteten Faltung einer Projektion nach dem Verfahren der schnellen "Faltung mit der nächsten Zweierpotenz",

Fig. 2 ein Beispiel einer fehlerbehafteten Faltung von p = 4 Projektion nach dem erfindungsgemäßen Verfahren der schnellen Faltung,

Fig. 3 einen Computertomographen zur Durchführung des erfin dungsgemäßen Verfahrens der gleichzeitigen schnellen Faltung gemessener Projektionen mit einem modifizier ten Faltungskern,

Fig. 4 ein Beispiel einer gleichzeitigen schnellen Faltung von 2p = 6 Projektionen mit Überfaltungsfehlern nach dem erfindungsgemäßen Verfahren,

Fig. 5 Bereiche mit geringfügig verfälschter Bildinformation bei der schnellen Faltung nach dem erfindungsgemäßen Verfahren für N = 1536, M = 8192, p = 3 und R = 25 cm,

Fig. 6 eine prinzipielle Darstellung der Fehlerkompensation bei der gleichzeitigen schnellen Faltung von p = 3 Pro jektionen,

Fig. 7 den absoluten Fehler e(k) der schnellen Faltung ge messener Projektionen eines Phantoms nach dem defini tionsgemäßen Verfahren bei der vollen Länge des Fil terkerns im Ortsbereich von 2N-1,

Fig. 8 den absoluten Fehler e(k) der schnellen Faltung ge messener Projektionen eines Phantoms nach dem erfin dungsgemäßen Verfahren bei der Länge des Filterkerns im Ortsbereich von N+S,

Fig. 9 den absoluten Fehler e(k) bei der schnellen Faltung gemessener Projektionen eines Phantoms nach dem er findungsgemäßen Verfahren bei einer verkürzten Länge des Filterkerns im Ortsbereich von (N+S)-4,

Fig. 10 ein rekonstruiertes Bild eines Phantoms anhand simu lierter Projektionen unter Anwendung des definitions gemäßen Verfahrens der schnellen Faltung,

Fig. 11 ein rekonstruiert es Bild eines Phantoms anhand simu lierter Projektionen unter Anwendung des erfindungs gemäßen Verfahrens der schnellen Faltung,

Fig. 12 das Differenzbild der beiden Bilder der Fig. 10 und 11 (maximale Differenz = ±1 HU),

Fig. 13 ein rekonstruiertes Bildes des Phantoms anhand des gleichen Datensatzes wie in Fig. 10 mit meßwertunab hängigem Rauschen unter Anwendung des definitionsge mäßen Verfahrens der schnellen Faltung,

Fig. 14 ein rekonstruiertes Bildes des Phantoms anhand des gleichen Datensatzes wie in Fig. 10 mit meßwertunab hängigem Rauschen unter Anwendung des erfindungsgemä ßen Verfahrens der schnellen Faltung,

Fig. 15 das Differenzbild der beiden Bilder der Fig. 13 und 14 (maximale Differenz = ±1 HU),

Fig. 16 ein rekonstruiertes Bild des Phantoms anhand des gleichen Datensatz wie in Fig. 10 mit meßwertabhängi gem Rauschen unter Anwendung des definitionsgemäßen Verfahrens der schnellen Faltung,

Fig. 17 ein rekonstruiertes Bild des Phantoms anhand des gleichen Datensatzes wie in Fig. 10 mit meßwertabhän gigem Rauschen unter Anwendung des erfindungsgemäßen Verfahrens der schnellen Faltung, und

Fig. 18 das Differenzbild der beiden Bilder der Fig. 16 und 17 (maximale Differenz = ±1 HU).

Die Fig. 3 zeigt einen Computertomographen 1 zur Durchführung des erfindungsgemäßen Verfahrens der gleichzeitigen schnellen Faltung gemessener Projektionen mit einem modifizierten Fal tungskern. Der Computertomograph 1 weist ein Meßsystem aus einer Röntgenstrahlenquelle 3, die ein fächerförmiges Rönt genstrahlenbündel 4 aussendet, und einem Strahlenempfänger 5 auf, welcher aus einer Reihe von Einzeldetektoren besteht. Der Fokus der Röntgenstrahlenquelle 3, von dem das Röntgen strahlenbündel 4 ausgeht, ist mit 2 bezeichnet. Der zu un tersuchende Patient P liegt auf einer Patientenliege 6.

Zur Durchführung einer radiologischen Untersuchung eines Pa tienten P wird das Meßsystem 3, 5 um ein Meßfeld 10, in dem der Patient P liegt, um 360° gedreht. Ein Motor 12 treibt hierzu den Drehtisch 11 an. Die Drehachse, welche rechtwink lig zu dem fächerförmigen Röntgenstrahlenbündel 4 steht, ist mit A bezeichnet. Dabei wird die Röntgenstrahlenquelle 3, die von einem Röntgengenerator 7 gespeist wird, gepulst oder mit Dauerstrahlung betrieben. Bei vorbestimmten Winkelpositionen α des Meßsystems 3, 5 werden Projektionen von Schichten des Patienten aufgenommen, wobei die zugehörigen Datensätze der Meßdaten vom Strahlenempfänger 5 einem Rechner 8 zugeführt werden, welcher aus den erzeugten Datensätzen die Schwä chungskoeffizienten vorbestimmter Bildpunkte berechnet und auf einem Monitor 9 bildlich wiedergibt. Auf dem Monitor 9 erscheint demgemäß ein Bild der durchstrahlten Schicht des Patienten P.

Auf dem Rechner 8 wird unter anderem auch die Berechnung der gleichzeitigen schnellen Faltung der aufgenommenen Projektio nen mit einem modifizierten Filterkern durchgeführt. Die Fig. 4 zeigt ein Beispiel einer solchen schnellen Faltung nach dem erfindungsgemäßen Verfahren mit 2p = 6 Projektionen der Länge N = 1536 Kanäle, wobei die ersten p Projektionen ₁(k) bis ₃(k) im Realteil und die zweiten p Projektionen ₄(k) bis ₆(k) im Imaginärteil des komplexen Vektors (k) der Länge M = 8192 abgelegt sind. Durch das erfindungsgemäße Verfahren können somit in einer Fast Fourier Transformation der Trans formationslänge M = 8192 = 2¹³ gleichzeitig sechs Projektio nen verarbeitet werden, wo hingegen bei der herkömmlichen definitionsgemäßen Vorgehensweise nur zwei Projektionen mit einer Fast Fourier Transformation der Transformationslänge M = 4096 verarbeitet werden könnten. Auf jede Projektion ₁(k) bis ₆(k) folgen jeweils S = M/p - N = 1194,7 Nullen (zweimal 1195 und einmal 1194). Die Länge des Filter kerns (k) im Ortsbereich beträgt nicht, wie im Falle der definitionsgemäßen schnellen Faltung, 2N-1 = 3071, sondern nur N+S = 1536 + 1195 = 2731 Werte. Die restlichen Werte des Filterkerns (k) sind mit Nullen aufgefüllt, so daß insge samt im vorliegenden Fall folgende Gleichungen gelten:

Somit ergeben sich nach der Faltung am Anfang und Ende jeder gefilterten Projektion nur (N-S-1)/2, im vorliegenden Fall also ungefähr 170 verfälschte Werte. Diese Werte liegen al lerdings im allgemeinen sehr nahe am richtigen Ergebnis der exakten Faltung, da der Anteil der bei der Überfaltung von einer anderen nahegelegenen Projektion zuviel addiert wird (vgl. Fig. 4, a), durch das Fehlen eines etwa zahlenmäßig gleichen Anteils (vgl. Fig. 4, b) aus der zu bearbeitenden Projektion wieder kompensiert wird (vgl. auch Fig. 6). Der Grund dafür liegt in der Ähnlichkeit aufeinanderfolgender Projektionen.

Fig. 5 veranschaulicht in diesem Zusammenhang die Bereiche in einem Meßfeld 10 (vgl. Fig. 3), welche bei einer Rekonstruk tion eines Bildes nach dem erfindungsgemäßen Verfahren der schnellen Faltung geringfügig verfälschte Bildinformationen aufweisen. Hat das Meßfeld 10 dabei einen Durchmesser von 2R = 50 cm, welchem im vorliegenden Fall N = 1536 Kanäle ent sprechen, dann entsprechen 170 Kanälen in etwa 5,5 cm, d. h., daß nur in dem in Fig. 5 schraffiert gezeichneten Ringbereich mit einer Dicke von 5,5 cm eine geringfügige Verfälschung der Bildinformation auftritt. Im inneren Bereich des Meßfeldes 10 mit einem Durchmesser von 39 cm (r = 19,5 cm) ergibt sich da gegen kein Unterschied im rekonstruierten Bild, wenn bei der Bildrekonstruktion anstelle des definitionsgemäßen Verfahrens der schnellen Faltung mit dem erfindungsgemäßen Verfahren der schnellen Faltung gearbeitet wird. Auch ein Zoom innerhalb dieses Bereiches liefert aufgrund der in diesem Bereich kor rekt berechneten schnellen Faltung ein fehlerfreies Bild.

Im folgenden wird in einem Vergleich der beiden Berechnungs verfahren der Faltung gemessener Projektionen einer Wasser scheibe vom Durchmesser 49 cm mit einem Filterkern nach dem definitionsgemäßen und dem erfindungsgemäßen Verfahren der schnellen Faltung die bei dem erfindungsgemäßen Verfahren der schnellen Faltung von p = 3 Parallelprojektionen der Länge N = 1536 Kanäle auftretende Fehlerkompensation betrachtet. In beiden Berechnungsfällen wird dabei mit dem an sich bekannten Shepp-Logan Kern gefaltet.

Um die auftretenden Effekte besser kenntlich machen zu kön nen, werden zunächst einige Vektoren definiert. Zuerst ein Vektor der Transformationslänge M, der nur eine Projektion _i(k) enthält.

Dann die zyklische Faltung dieses Vektors mit der nach Glei chung 3 modifizierten Impulsantwort des Filterkerns.

Addiert man alle p Signale auf, so erhält man aus Lineari tätsgründen das gleiche Faltungsergebnis wie bei der schnel len Faltung mit (k) nach Gleichung 4.

Die gefaltete einzelne Projektion y_i(k) ist dann:

y_i(k) = y(k+(i-1) (N+S)) für k = 0(1)N-1 (10)

Die ohne Fehler gefaltete Projektion ergibt sich nach folgen der Gleichung:

Repräsentativ für alle gefalteten Projektionen wird nun der Fehler bei Projektion 1 näher untersucht. Der Gesamtfehler ergibt sich dabei aus der Differenz der erfindungsgemäßen schnellen Faltung der Projektion 1 zur definitionsgemäßen schnellen Faltung der Projektion 1 nach Gleichung 11.

e(k) = y₁(k) - y_c1(k) = y(k) - y_c1(k) für k = 0(1)N-1 (12)

Dieser Fehler setzt sich aus drei Anteilen zusammen. Der er ste Anteil besteht dabei aus dem Fehler

(k) = ₁(k) - y_c1(k) für k = 0(1)N-1, (13)

welcher aus der zu kurzen Länge des Filterkerns im Ortsbe reich herrührt, während die Fehler 2 und 3 aus den Ausschwin gern von ₂(k) und ₃(k) der beiden anderen gefilterten Pro jektionen herrühren, die zu Überfaltungen führen. Fig. 6 zeigt im übrigen den prinzipiellen, nicht maßstäblichen Verlauf dieser Faltung, um die Effekte aufzuzeigen.

Die Fig. 7 bis 9 veranschaulichen den absoluten Fehler e(k), welcher bei unterschiedlichen Längen des Filterkerns im Orts bereich bei der Berechnung nach dem erfindungsgemäßen Verfah ren der schnellen Faltung gegenüber dem definitionsgemäßen Verfahren der schnellen Faltung auftritt. Fig. 7 zeigt den ab soluten Fehler e(k) für die volle Länge des Filterkerns im Ortsbereich von 2N-1, welcher bei Anwendung des erfindungs gemäßen Verfahrens der schnellen Faltung auftritt. Dagegen veranschaulicht Fig. 8 den absoluten Fehler e(k), wenn nach dem erfindungsgemäßen Verfahren der schnellen Faltung mit der verkürzten Länge des Filterkerns im Ortsbereich von N+S ge rechnet wird. Hierbei addieren sich die drei oben angegebenen Fehler zu dem deutlich kleineren Fehler e(k). Durch eine wei tere geringfügige Verkürzung der Länge des Filterkerns im Ortsbereich auf N+S-4 (vgl. Fig. 9) kann der absolute Fehler e(k) der schnellen Faltung gemäß dem erfindungsgemäßen Ver fahren gegenüber einem nicht in der Länge modifizierten Fil terkern für die Wasserscheibe noch weiter verbessert werden.

Zur Verdeutlichung, daß sich die Fehler, welche bei der Be rechnung der schnellen Faltung nach dem erfindungsgemäßen Verfahren auftreten, nur unwesentlich in der Bilddarstellung in der Computertomographie auswirken, sind in den Fig. 10 bis 18 die rekonstruierten Bilder eines simulierten Phantoms dar gestellt. Das Phantom wird von einer Wasserscheibe (0 HU) mit einem Durchmesser von 49 cm und von einem Knochen (1000 HU) mit einem Durchmesser von 2 cm bei einem Abstand des Mit telpunktes des Knochens von 23 cm vom Mittelpunkt der Wasser scheibe gebildet. Bei der Simulation wurden n=528 Parallel projektionen einer Länge N = 1536 Kanäle aufgenommen. Die Re konstruktion der in den Fig. 10 bis 18 gezeigten Bilder er folgte unter Verwendung des Shepp-Logan Filterkerns.

Während Fig. 10 ein rekonstruiertes Bild des Phantoms anhand der n = 528 gemessenen Parallelprojektionen unter Verwendung des definitionsgemäßen Verfahrens der schnellen Faltung zeigt, ist in Fig. 11 ein rekonstruiertes Bild des Phantoms anhand der n = 528 Parallelprojektionen unter Verwendung des erfindungsgemäßen Verfahrens gezeigt, bei dem, wie beschrie ben, jeweils drei in Real- und Imaginärteil aufgespaltenen Projektionen p gleichzeitig gefaltet wurden. Fig. 12 zeigt da bei das Differenzbild der beiden in Fig. 10 und 11 dargestell ten Bilder, wobei sich die beiden Bilder in den Fig. 10 und 11 maximal um ±1 HU unterscheiden. Der Fehler tritt dabei nur, wie beschrieben, im äußeren Kreisring von 5 cm auf. Des wei teren wurde der gleiche Datensatz in den Fig. 13 bis 15 mit meßwertunabhängigem und in den Fig. 16 bis 18 mit meßwertab hängigem Rauschen untersucht, wobei in allen Fällen, wie den Fig. 15 und 18 zu entnehmen ist, in den Differenzbildern des definitionsgemäßen Verfahrens zu dem erfindungsgemäßen Ver fahren die Fehler ± 1 HU sind.

Das erfindungsgemäße Verfahren der schnellen Faltung verletzt zwar die Definition der schnellen Faltung, aber solange der auftretende Fehler unterhalb einer gewissen Schwelle liegt, ist er in der Computertomographieanwendung vernachlässigbar. Das bedeutet, daß ein Fehler, der sich im rekonstruierten Bild einer gemessenen Projektion im Bereich bis zu einer HU auswirkt, durchaus toleriert werden kann.

Das erfindungsgemäße Verfahren der schnellen Faltung kann die Rechenzeit erheblich reduzieren. Im Falle des vorstehend beschriebenen Beispiels mit der Transformationslänge M = 8192 und 2p = 6 Projektionen kann man eine Reduktion der Multiplika tion um etwa 28% pro Projektion gegenüber der definitionsge mäßen Faltung mit M = 4096 und zwei Projektion erreichen (vgl. Tabelle 1).

In der Tabelle 1 wird dabei angenommen, daß die definitions gemäße schnelle Faltung der Transformationslänge M mit einem reellen Filterkern 4M log₂ M + 2M reelle Multiplikationen hat.

Tabelle 1

Reduktion der Zahl der Multiplikationen pro Pro jektion bei dem erfindungsgemäßen Verfahren der schnellen Faltung am Beispiel 2p = 6

Im übrigen ist eine Voraussetzung für die Anwendung des er findungsgemäßen Verfahrens, daß eine Prozessor- und Hardware architektur zur Verfügung steht, bei der eine Verdopplung der Transformationslänge der Fast Fourier Transformation nicht zu einem überproportionalen Anstieg der Berechnungszeit führt.

Des weiteren ist die Anwendung des erfindungsgemäßen Verfah rens der schnellen Faltung nicht auf die Computertomographie beschränkt, sondern es sind auch andere Anwendungsfälle denk bar. So besteht die Möglichkeit, das erfindungsgemäße Verfah ren einzusetzen, wenn die Impulsantwort eines Filterkerns im Zeitbereich am Rand betragsmäßig relativ klein ist oder/und aufeinanderfolgende Vektoren der zu faltenden Signale ähnlich sind.

Darüber hinaus arbeitet das erfindungsgemäße Verfahren der schnellen Faltung nicht nur mit dem Shepp-Logan-Filterkern, sondern auch mit anderen Filterkernen, z. B. dem an sich be kannten Kern von Ramachandran und Lakshminarayanan.

Claims

Verfahren für die Bildrekonstruktion in der bildgebenden Technik zur Durchführung einer schnellen Faltung mit der Transformationslänge M unter Zulassung geringer Überfaltungs fehler, wobei jeweils p (oder 2p) von n gemessenen Projektio nen der Länge N mit einem modifizierten Filterkern (k) un ter Verwendung der "Fast Fourier Transform" (FFT) und der "Inverse Fast Fourier Transform" (IFFT) gleichzeitig in einem Schritt gefaltet werden, umfassend folgende Schritte:
- a) Zusammenfassung der p Projektionen der Länge N zu einem Vektor und Anfügung von S Nullen hinter jeder der p Pro jektionen der Länge N zu einem Vektor mit der Transformationslänge M,
- b) Wahl der Parameter p und S derart, daß die Transformationslänge M = p (N+S) = 2^m (2)zur Durchführung der schnellen Faltung eine Zweierpotenz ist,
- c) Entwurf des Filterkerns (k) im Ortsbereich derart, daß er wenigstens im wesentlichen gleich einer Länge N+S ist, und Auffüllen des Filterkerns (k) mit Nullen auf die Transformationslänge M gemäß
- d) Zyklische Durchführung der schnellen Faltung mit der Transformationslänge M gemäß mit k = 0(1) (M-1),wobei das Faltungsergebnis einer Projektion i die Formy_i(k) = y(k + (i-1) (N+S)) (5)mit k = 0(1) (N-1) und
  i = 1(1)paufweist.