DE102005053498B4

DE102005053498B4 - Verfahren zur Beschleunigung der Streustrahlungskorrektur in einem Computertomographiesystem sowie Computertomographiesystem zur Durchführung des Verfahrens

Info

Publication number: DE102005053498B4
Application number: DE102005053498A
Authority: DE
Inventors: Ernst-Peter Rührnschopf; Bernhard Dr. Scholz; Alois Nöttling
Original assignee: Siemens AG
Current assignee: Siemens Healthcare GmbH
Priority date: 2005-11-09
Filing date: 2005-11-09
Publication date: 2008-07-24
Anticipated expiration: 2025-11-10
Also published as: US7801266B2; DE102005053498A1; US20070104310A1

Abstract

Verfahren zur Beschleunigung der Streustrahlungskorrektur in einem Computertomographiesystem (1) mit einer Strahlungsquelle (2) und einem großflächig ausgebildeten Detektor (3) mit einer Vielzahl von Detektorzeilen, wobei ein Objekt (7) aus einer Vielzahl von Projektionswinkeln abgetastet wird, aus den gemessenen Schwächungswerten der Strahlung Projektionsdaten erstellt werden, die zur Rekonstruktion tomographischer Darstellungen weiterverarbeitet werden, wobei eine Strahlaufhärtungskorrektur direkt auf den Projektionsdaten p_T vorgenommen wird, dadurch gekennzeichnet, dass die Streustrahlungskorrektur ebenfalls direkt auf den Projektionsdaten p_T erfolgt und hierfür ein in der Strahlaufhärtungskorrektur verwendeter Aufhärtungskorrekturterm mit einem Verstärkungsfaktor multipliziert wird.

Description

Die Erfindung betrifft ein Verfahren zur Beschleunigung der Streustrahlungskorrektur in einem Computertomographiesystem mit einer Strahlungsquelle und einem großflächig ausgebildeten Detektor mit einer Vielzahl von Detektorzeilen, wobei ein Objekt aus einer Vielzahl von Projektionswinkeln abgetastet wird, aus den gemessenen Schwächungswerten der Strahlung Projektionsdaten erstellt werden, die zur Rekonstruktion tomographischer Darstellungen weiterverarbeitet werden, wobei eine Strahlaufhärtungskorrektur direkt auf den Projektionsdaten p vorgenommen wird. Die Erfindung betrifft weiterhin ein Computertomographiesystem zur Durchführung dieses Verfahrens.
Die Computer-Tomographie (CT) ist ein Verfahren, mit dem die 3-dimensionale Rekonstruktion der Gewebedichteverteilung oder genauer definiert die 3-dimensionale Verteilung des linearen Schwächungskoeffizienten für Röntgenstrahlung des Körperinneren mit hoher räumlicher Auflösung und großer quantitativer Genauigkeit der Dichteauflösung möglich ist.
Es ist heute technisch möglich, mit kurzer Aufnahmezeit einen großen Volumenbereich abzubilden und zu rekonstruieren. Hierzu werden meist entweder Flachbilddetektoren (Flat-Panel-Detectors = FPD) verwendet, mit denen bereits mit einer einzigen "Halbrotation", genauer einer Rotation über > 180° + Fächerwinkel, durch Conebeam-Computertomographie (CBCT) aussagekräftige CT-Darstellungen begrenzter Volumina rekonstruiert werden können. Andererseits werden Vielzeilen-Detektoren mit Mehrfachrotationen um den Patienten und kontinuierlichem Vorschub verwendet, die durch die so genannte Spiral-CT den vollständigen Körper eines Patienten tomographisch darstellen können.
Ein schwerwiegendes Problem ist die im Aufnahmeobjekt, meist dem Patienten, erzeugte Streustrahlung, die mit der Fläche der Detektoren und dem gleichzeitig entsprechend größeren durchstrahlten Volumen zunimmt. Dabei kann die Streustrahlung in ihrer Intensität die Größenordnung der ungestreuten direkten Primärstrahlung erreichen und in Extremfällen sogar überwiegen. Die Folge ist eine Verzerrung der quantitativen Gewebedichte-Rekonstruktion, die sich typischerweise z. B. bei einem homogenen zylindrischen Körper durch eine kontinuierliche Abdunkelung vom Rand zur Mitte auffällig bemerkbar macht. Dieser Effekt wird "Cupping" (=Schüsseleffekt) genannt. Dieser Schüsseleffekt kann die Erkennbarkeit kleiner pathologischer Strukturen niedrigen Kontrastes sehr beeinträchtigen, zumal Fehler von mehreren 100 HU Abweichung zwischen Objektrand und -mitte auftreten können. 1 HU (= Hounsfield Unit) entspricht dem 1/1000 der Dichte von Wasser.
Bei den herkömmlichen CTs mit Einzeilen- oder Mehrzeilen-Detektoren wird die Streustrahlung durch Kollimation wirkungsvoll unterdrückt. Dies geschieht zum Beispiel durch Einblendung der Aufnahmeschichtdicke bei Einzeilendetektoren oder des Schichtvolumens beim Mehrzeilendetektor durch axiale Kollimation in Richtung der Rotationsachse. Es kann auch eine Kapselung beziehungsweise Verwendung seitlich angeordneter Kollimatorbleche zwischen benachbarten Detektorelementen, also eine laterale Kollimation, genutzt werden.
Bei Flächendetektoren, wie sie vornehmlich bei modernen C-Bogen-Geräten verwendet werden, wird die Strahlung, schon aus Gründen des Strahlenschutzes, durch Kollimatoren auf das Messfeld – in der Regel die Detektorfläche, in manchen Applikationen auch nur auf eine Teilfläche – eingeschränkt. Die im Patienten entstehende Streustrahlung kann durch diese Maßnahmen aber nicht unterdrückt werden. Im Übrigen besteht die Möglichkeit einer lateralen Kollimation beim Flächendetektor nicht.
Es können auch vor die Detektoreingangsfläche montierte Streustrahlenraster (anti-scatter grids) verwendet werden. Ihr Nutzen für die CBCT-Bildgebung wird zwar, wie aus der Schrift von J. H. Siewerdsen, D. H Moseley, B. Bakhtiar, S. Richard, D. A. Jaffray: „The influence of antiscatter grids an soft-tissue detectability in cone-beam computed tomography with flat-panel detectors", Med. Phys. 31(12), December 2004, 3506 bis 3520, hervorgeht, kontrovers diskutiert, der Einsatz ist aber zumindest bei hohem Streustrahlungsanteil zu empfehlen. In der Regel reicht jedoch die Reduktion der Streustrahlung durch Streustrahlenraster nicht aus, so dass zusätzliche Streustrahlungskorrekturverfahren erforderlich sind.
In der Schrift von R. Ning, X. Tang, D. L. Conover: „X-ray scatter suppression algorithm for cone beam volume CT", Proc. SPIE, Vol. 4682, 2002, 774 bis 781, wurde vorgeschlagen, die Streustrahlung in wenigen Projektionsrichtungen zu messen und zu korrigieren. Für die Messung der Streustrahlung muss dabei eine Vorrichtung vorhanden sein, mit der eine "Strahlstopp"-Trägerplatte patientennah in den Strahlengang zwischen Röntgenquelle und Patient eingefahren werden kann. Diese erforderliche zusätzliche Messprozedur erscheint jedoch für den normalen klinischen Arbeitsablauf nicht akzeptabel.
In der Schrift von M. Zellerhoff, B. Scholz, E.-P. Rührnschopf, T. Brunner: „Low contrast 3D reconstruction from Carm data", Proceedings of SPIE, Medical Imaging 2005, Vol. 5745, Seiten 646 bis 655, wird ein rein rechnerischer Streustrahlungskorrekturalgorithmus für CBCT vorgeschlagen. Der Algorithmus ist zwar relativ schnell, die Verlängerung der gesamten Rekonstruktionszeit ist jedoch nicht vernachlässigbar. So benötigt der Algorithmus für 500 Projektionsbilder der Größe 1024×1024 Pixel auf einem PC mit 3 GHz Taktzeit immer noch etwa 1 Minute.
Seitens der klinischen Applikationen, insbesondere durch den zunehmenden Einsatz von C-Bogen-Conebeam-CT-Anlagen mit FPD in der Interventionellen Medizin, besteht die dringende Forderung nach möglichst kurzer Rechenzeit für die Volumen-Rekonstruktion. Zusätzlich implementierte Korrekturalgorithmen zur Verbesserung der Bildqualität werden nur dann akzeptiert, wenn sie die Rechenzeit praktisch nur unmerklich verlängern.
Auf die Druckschrift EP 1 566 771 A1 wird verwiesen. Diese offenbart eine einfache Streustrahlungskorrektur auf Projektionsdaten, allerdings wird der Effekt der Strahlungsaufhärtung in dieser Schrift nicht berücksichtigt.
Es ist daher Aufgabe der Erfindung, ein Verfahren zur Streustrahlungskorrektur zu finden, welches einen geringeren Rechenaufwand bei im Wesentlichen gleich bleibenden Bildergebnissen erfordert.
Diese Aufgabe wird durch die Merkmale der unabhängigen Patentansprüche gelöst. Vorteilhafte Weiterbildungen der Erfindung sind Gegenstand untergeordneter Ansprüche.
Die Erfinder haben erkannt, dass der Aufwand des letztgenannten rechnerischen Streustrahlungskorrekturalgorithmus zum großen Teil auf die je Pixelwert erforderlichen Operationen Delogarithmierung und Logarithmierung im Korrekturschritt und auf eine 2-dimensionale Faltung beziehungsweise Bildexpansion innerhalb des verwendeten Streustrahlungs-Schätzmodells zurückzuführen ist und daher ein Korrekturverfahren gefunden werden sollte, welches auf derart aufwendige Rechenschritte bei der Korrektur jedes Pixels verzichtet.
Weiterhin haben die Erfinder gesehen, dass eine ganz ähnliche Bildauswirkung wie durch die Streustrahlung durch die spektrale Strahlaufhärtung stattfindet. Obwohl physikalisch völlig anderer Art, können die Bildwirkungen des Cuppings beider Ursachen praktisch nicht von einander unterschieden werden. Die Aufhärtungskorrektur lässt sich auf einfachen „Table-lookup" für jedes Pixel jedes Projektionsbildes zurückführen und ist daher sehr schnell.
Die Idee besteht nun darin, den Streustrahlungseffekt durch eine entsprechende Verstärkung der Aufhärtungskorrektur zu realisieren und dabei die Schnelligkeit der Aufhärtungskorrektur auszunutzen. Eine solche Streustrahlungskorrektur wäre sozusagen zum rechnerischen „Nulltarif" erreichbar. Insbesondere würden die aufwendigen Operationen Delogarithmierung und Logarithmierung und Bildexpansion entfallen.
Bezüglich bekannter Strahlaufhärtungskorrekturverfahren wird beispielhaft auf die Schriften von H. H. Barrett, W. Swindell: Radiological Imaging, Vol. 2, Academic Press, 1981, Seiten 438 bis 441, und von A. C. Kak, M. Slaney: Principles of Computerized Tomographi Imaging, IEEE Press New York, 1979, Seiten 120 bis 123, verwiesen.
Als Hauptproblem bleibt die Bestimmung des streustrahlungsabhängigen Verstärkungsfaktors. Die dazu nötigen Informationen können beispielsweise auf die gleiche Art gewonnen werden, wie es in der oben genannten Schrift von M. Zellerhoff, B. Scholz, E.-P. Rührnschopf, T. Brunner: „Low contrast 3D reconstruction from C-arm data", Proceedings of SPIE, Medical Imaging 2005, Vol. 5745, Seiten 646 bis 655, beschrieben wird. Es wird dabei eine Datenbasis benutzt, welche für typische Streukörper – Wasserzylinder mit elliptischem Querschnitt – die S/P-Verhältnisse (Intensität der Streustrahlung im Verhältnis zur Primärstrahlung) in Abhängigkeit aller die Streustrahlung beeinflussenden Parameter wie Dicke des Streukörpers, Luftspalt = Abstand des Streukörpers zum Detektor, Detektorfeldgröße, Spannung der Röntgenröhre, Dicke und Material der röhrenseitigen Vorfilterung und Streustrahlenraster enthält. Diese Datenbasis, abgelegt in einer mehrdimensionalen Tabelle, kann im Voraus zum Beispiel durch aufwendige Monte-Carlo-Simulationsrechnungen und Abgleichung mit Messungen erstellt werden.
Es ist anzumerken, dass das Korrekturverfahren, da es direkt auf den Projektionsdaten operiert, unabhängig von der Aufnahmegeometrie ist und prinzipiell jeder für die Aufnahmegeometrie geeignete Rekonstruktionsalgorithmus verwendet werden kann. Jedoch muss diese Aufnahmegeometrie bei der Erstellung der S/P-Datenbasis berücksichtigt werden.
Entsprechend der hier geschilderten Grundidee schlagen die Erfinder ein Verfahren zur Beschleunigung der Streustrahlungskorrektur in einem Computertomographiesystem mit einer Strahlungsquelle und einem großflächig ausgebildeten Detektor mit einer Vielzahl von Detektorzeilen vor, wobei hierbei ein Objekt aus einer Vielzahl von Projektionswinkeln abgetastet wird, aus den gemessenen Schwächungswerten der Strahlung Projektionsdaten erstellt und zur Rekonstruktion tomographischer Darstellungen weiterverarbeitet werden, wobei eine Strahlaufhärtungskorrektur direkt auf den Projektionsdaten p_T vorgenommen wird. Erfindungsgemäß erfolgt dabei die Streustrahlungskorrektur ebenfalls direkt auf den Projektionsdaten, wobei zur Streustrahlungskorrektur ein in der Strahlaufhärtungskorrektur ohnehin verwendeter Aufhärtungskorrekturterm mit einem Verstärkungsfaktor multipliziert wird. Hierdurch wird erreicht, dass die aufwendigen Rechenoperationen des Logarithmierens und Delogarithmierens vermieden werden und so die Rechengeschwindigkeit insgesamt beschleunigt wird.
Eine konkrete Ausführung des Verfahrens kann darin bestehen, dass für die kombinierte Aufhärtungs- und Streustrahlungskorrektur korrigierte Projektionsdaten p^* = p_T + ρ·δ_a mit dem Verstärkungsfaktor ρ = 1 + 1/δ_a log(l + S/P_a) berechnet werden, wobei δ_a den Korrekturterm für die Strahlaufhärtungskorrektur, S = I_S/I_o die normierte Streustrahlungsverteilung und P_a = I_A/I_o die normierte Primärstrahlung mit Aufhärtung ohne Streustrahlung beschreiben.
Hierbei kann das Verhältnis S/P_a von normierter Streustrahlung S zur normierten Primärstrahlung mit Aufhärtung ohne Streustrahlung P_a unter Berücksichtigung mindestens eines repräsentativen Pixels geschätzt werden. Alternativ kann das Verhältnis S/P_a auch aus Simulationsdaten eines Phantoms mit zum untersuchten Objekt ähnlicher Form, Größe und Struktur entnommen werden. Eine weitere Möglichkeit besteht darin, das Verhältnis S/P_a durch ein iteratives Verfahren anzunähern.
Vorteilhaft ist weiterhin, wenn der Verstärkungsfaktor ortsabhängig relativ zum Detektor bestimmt wird. Dabei kann der Verstärkungsfaktor auch in Abhängigkeit von der betrachteten Detektorzeile bestimmt werden.
Im Falle nichtrotationssymmetrischer Objekte zumindest im Scanbereich schlagen die Erfinder auch vor, dass der Verstärkungsfaktor zusätzlich projektionswinkelabhängig bestimmt wird.
Gemäß einer weiteren besonderen Ausführung des Verfahrens kann der Verstärkungsfaktor nochmals mit einer Funktion multipliziert werden, deren Wertebereich zwischen 0 und 1 liegt und deren Wert zum Objektrand auf 0 abfällt. Hierbei ist es auch vorteilhaft, die Funktion mit einem Exponenten zu potenzieren, so dass der Werteabfall zum Rand des Objektes hin verstärkt wird.
Entsprechend dem oben dargelegten Verfahren wird auch ein CT-System vorgeschlagen, welches mit einer Strahlungsquelle und einem großflächig ausgebildeten Detektor mit einer Vielzahl von Detektorzeilen ausgestattet ist, mit dem ein Objekt aus einer Vielzahl von Projektionswinkeln abgetastet wird, wobei aus den gemessenen Schwächungswerten der Strahlung Projektionsdaten erstellt werden, und einer Steuer- und Recheneinheit, die mindestens einen Speicher für Programme und Daten aufweist, und mit gespeicherten Programmen und Projektionsdaten eine Rekonstruktion tomographischer Darstellungen durchführt. Erfindungsgemäß sind im CT-System Computerprogramme oder Programm-Module gespeichert, die im Betrieb die Verfahrensschritte eines der voranstehenden Verfahrensansprüche durchführen.
Im Folgenden wird die Erfindung anhand bevorzugter Ausführungsbeispiele mit Hilfe der Figuren näher beschrieben, wobei nur die zum Verständnis der Erfindung notwendigen Merkmale dargestellt sind. Hierzu werden die nachfolgenden Bezugszeichen verwendet: 1: CT-System/C-Bogen-Gerät; 2: Röntgenröhre; 3: Detektor; 4: Systemachse; 5: Gantrygehäuse/C-Bogen; 6: Patientenliege; 7: Patient; 8: Öffnung im Gantrygehäuse; 9: Steuer- und Recheneinheit; 10: Speicher; 11: idealer Verlauf der Projektionsdaten ohne Streustrahlung und Strahlaufhärtung; 12: Verlauf der Projektionsdaten mit Streustrahlung und Strahlaufhärtung; 13: Verlauf der Projektionsdaten ohne Streustrahlung mit Strahlaufhärtung; 14: Verlauf des Strahlaufhärtungsanteil; 15: Verlauf des Streustrahlanteils; 16: erfindungsgemäß korrigierte Projektionsdaten.
Es zeigen im Einzelnen:
1 Schematische Darstellung eines CT-Systems mit Gantry zur Durchführung des erfindungsgemäßen Verfahrens;
2 Schematische Darstellung eines C-Bogen-Gerätes zur Durchführung des erfindungsgemäßen Verfahrens;
3 Gegenüberstellung idealer und mit dem erfindungsgemäßen Verfahren korrigierter Projektionsdaten im Projektionswinkel 0°; und
4 Gegenüberstellung idealer und mit dem erfindungsgemäßen Verfahren korrigierter Projektionsdaten im Projektionswinkel 90°.
Die 1 zeigt beispielhaft in 3D-Darstellung ein erfindungsgemäßes CT-System 1, welches zur Durchführung des erfindungsgemäßen Verfahrens geeignet ist. Es besteht in bekannter Weise aus einem Gantrygehäuse 5, in dem die hier nicht sichtbare Gantry angeordnet ist, auf der eine Röntgenröhre 2 als Strahlenquelle und ein großflächiger Vielzeilendetektor 3 befestigt sind. Zum Scan des auf der Patientenliege 6 liegenden Patienten 7 wird dieser kontinuierlich für einen Spiralscan oder diskontinuierlich für einen Kreisscan durch die Öffnung 8 im Gantrygehäuse 5 entlang der Systemachse 4 geschoben, während die auf der Gantry befestigte Röntgenröhre 2 und der gegenüberliegende Detektor 3 um die Systemachse 4 rotieren. Die Steuerung des Systems erfolgt durch eine Steuer- und Recheneinheit 9, die auch über einen Speicher 10 verfügt, in dem die Steuer- und Auswerteprogramme Prg₁-Prg_n des CT-Systems 1 gespeichert sind, die zum Betrieb des Systems aufgerufen und abgearbeitet werden. Bevorzugt befinden sich hier auch die Computerprogramme Prg_x und abgespeicherten Look-up-Tabellen, die das erfindungsgemäße Verfahren ausführen.
Die bevorzugte Ausführungsform eines Röntgentomographiesystems, für die das hier beschriebene erfindungsgemäße Verfahren verwendet wird, ist in der 2 in Querschnitt dargestellt. Dort wird ein grundsätzlich bekanntes C-Bogen-CT-Gerät 1 gezeigt, welches an seinem C-Bogen 5 eine Röntgenröhre 2 aufweist, der gegenüber ein großflächiger Detektor 3 angeordnet ist. Röntgenröhre 2 und Detektor 3 können mit dem C-Bogen 5 um den im Strahlungsfeld auf einer Patientenliege 6 angeordneten Patienten 7 und die Systemachse 4 geschwenkt werden, wobei aus den Detektorausgangsdaten ebenso wie im konventionellen CT mit Gantry tomographische Bilder des Patienten rekonstruiert werden können. Die Steuerung, Auswertung der Detektordaten und Rekonstruktion der tomographischen Darstellungen erfolgt auch hier in einer Steuer- und Recheneinheit 9, die einen Speicher 10 zur Abspeicherung bekannter und erfindungsgemäßer Computerprogramme Prg₁-Prg_n und Datentabellen aufweist. Derartige C-Bogen-Geräte werden vornehmlich in der Koronar-Angiographie verwendet.
Durch den hier sehr großflächig ausfallenden Detektor 3 und die kaum realisierbare Anwendung von Kollimatoren, wie sie in konventionellen CT-Systemen mit ein- bis mehrzeiligen Detektoren zur Streustrahlreduktion genutzt werden, erfordert zusätzliche Maßnahmen, um die Bildwirkung von Streustrahlung rechnerisch zu unterdrücken. Da gleichzeitig jedoch die Notwendigkeit eines schnellen Bildaufbaues besteht und die Re chenleistung begrenzt beziehungsweise teuer ist, wird das nachfolgend nochmals hergeleitete und begründete erfindungsgemäße Verfahren zur Streustrahlungskorrektur vorgeschlagen.
Um die Grundidee der Erfindung darzustellen, ist es zunächst notwendig, das Prinzip der Strahlaufhärtungskorrektur zu erklären. Es werden dabei die folgenden Bezeichnungen verwendet:

Q_U(E): spannungsabhängiges Emissionsspektrum,
τ_F(E) = exp(–μ_F(E)T): Transparenz des verwendeten spektralen Filters (z. B. Kupfer),
η_D(E): spektrale Ansprechempfindlichkeit des Detektors,
W_U(E) = Q_U(E)τ_F(E)η_D(E)/c_U: effektive normierte Spektralverteilung,
Normierfaktor und
eU: maximale Energie des Röntgenspektrums in eV.

Bei der Wasserkorrektur wird vereinfachend angenommen, die Schwächung eines Röntgenstrahls im Objekt (Patienten) werde allein durch wasseräquivalentes Material hervorgerufen. Dabei bedeutet "wasseräquivalent": die Energieabhängigkeit des Massenschwächungskoeffizienten (μ/ρ)(E) sei identisch zu Wasser und Unterschiede beruhten nur auf lokalen Dichteunterschieden.
Betrachtet man nun einen Messstrahl, der das Objekt durchdringt, wobei die Koordinate längs seines Weges mit x bezeichnet wird und der ortsabhängige lineare energieabhängige Schwächungskoeffizient mit μ(x, E) = ρ(x) α(x, E) angegeben ist, so errechnet sich der Massenschwächungskoeffizient α mit: α(x, E) = μ(x, E)/ρ(x).
Der Projektionswert zum betrachteten Messstrahl ist dann:
wobei der Index_a auf „Aufhärtung" hinweisen soll.
Dazu wird eine "äquivalente Wasserdicke" b_U = b_U(p_a) auf folgende Weise bestimmt:
α_W(E) sei der energieabhängige Massenschwächungskoeffizient von Wasser;
der polychromatische logarithmische Projektionswert zu einem Messstrahl, der längs einer Weglänge (Belegungsdicke) b in Wasser (ρ = 1 g/cm³) geschwächt wird, ist:
Diese Funktion kann für jede Spannung U vorausberechnet oder auch experimentell bestimmt werden. Sie steigt monoton mit b und kann invertiert werden, z. B. durch inverse Interpolation.
Zu jedem Messwert p_a kann eine äquivalente Wasserdicke b_U = b_U (p_a) = b so bestimmt werden, dass p_a = f_U(b) gilt, nämlich durch Inversion der letzten Gleichung und es ergibt sich: bU = f–1U (pα).
Mit b_U ist es dann möglich, auf den entsprechenden Projektionswert umzurechnen, der im Idealfall bei einem monochromatischen "Spektrum" mit Photonen nur einer einzigen Referenzenergie E₀ gemessen worden wäre:
mit bU = f–1U (pα) ergibt sich der aufhärtungskorrigierte Projektionswert zu * = αW(E0)bU = αW(E0)f–1U p = αW(E0)bU = αW(E0)f–1U (pα) = pα + δα mit δα = αW(E0)f–1U (pα) – pα
Der Term δa stellt die explizite Form des Aufhärtungskorrekturterms dar, wie er nachfolgend benutzt wird.
Dies vorweg geschickt soll im Folgenden die Grundidee der Erfindung dargestellt werden.
Es werden die folgenden Bezeichnungen verwendet:

(y, z):: Koordinaten auf dem Detektor;
I_o = I_o (y, z):: „Kalibrierbild” = gemessene „ungeschwächte" Intensitätsverteilung ohne Objekt im Strahlengang;
I_T = I_T (y, z):: gemessene totale Intensitätsverteilung mit Objekt im Strahlengang einschließlich Streustrahlung;
I_A = I_A (y, z):: Primär-Intensitätsverteilung mit Objekt im Strahlengang mit Aufhärtung, aber ohne Streustrahlung;
I_S = I_S (y, z):: Intensitätsverteilung nur der vom Objekt erzeugten Streustrahlung;
T = I_T/I_o:: normierte totale Intensitätsverteilung mit Objekt im Strahlengang einschließlich Streustrahlung;
P_a = I_A/I_o:: normierte Primärstrahlungsverteilung mit Aufhärtung, aber ohne Streustrahlung;
S= I_S/I_o:: normierte Streustrahlungsverteilung.

Dabei ist zu beachten, dass nur die Verteilungen I_o und I_T beziehungsweise T durch Messung gegeben sind. Für die Skizzierung der Grundidee wird von der vereinfachenden Annahme ausgegangen, dass auch S und P_a bekannt seien.
Es gilt: T = Pa + S (#1)
Die CT-Bildrekonstruktion erfolgt mit negativen logarithmierten normierten Intensitätsdaten, welche üblicherweise „Projektionsdaten" genannt und hier mit dem kleinen Buchstaben p bezeichnet werden. Aus der Gleichung (#1) folgt dann: pT = –log (Pa + S) = –log (Pa(1 + S/Pa)) = –log (Pa) – log (1 + S/Pa)oder kurz: PT = pa – δs (#2)mit δs = log(1 + S/Pa) (#3)
Die Gleichung (#3) beschreibt den Streustrahlungs-Korrekturterm, der zu den durch Streustrahlung gestörten Projektionsdaten p_T addiert werden muss, um die – nur noch durch Aufhärtung gestörten – Projektionsdaten p_a zu erhalten: pa = pT + δs (#4)
Die aufhärtungskorrigierten und für die Rekonstruktion zu verwendenden Projektionsdaten werden mit p^* bezeichnet. Die Aufhärtungskorrektur beziehungsweise Wasserkorrektur führt p_a in p^* über: p* = pa + δa (#5)dabei ist δa der oben für die Strahlaufhärtungskorrektur hergeleitete Aufhärtungskorrekturterm.
Will man die Streustrahlungskorrektur durch „Verstärkung" der Aufhärtungskorrektur formulieren, dann folgt aus den Gleichungen (#4) und (#5) p* = pT + δa + δs p* = pT + ρ·δa (#6)mit dem „Verstärkungsfaktor" ρ = 1 + δs/δa (#6a) ρ = 1 + 1/δa·log(1 + S/Pa)
Man beachte, dass alle Größen in den vorangegangenen Gleichungen ortsabhängig sind, d. h. Funktionen der Ortskoordinaten (y, z) auf dem Detektor sind. Dies gilt insbesondere für den „Verstärkungsfaktor" ρ = ρ(y, z) in der Gleichung (#6a).
Eine wesentliche Erkenntnis der Erfindung liegt darin, dass die beiden Terme in Gleichung (#6a), der Aufhärtungsterm δa(y, z) und der Streustrahlungsterm δs(y, z), einen ähnlichen räumlichen Verlauf aufweisen. Betrachtet man beispielsweise den Verlauf der beiden Korrekturterme beim Scan eines Zylinders, so verlaufen beide von der Mitte aus zunächst sehr flach und zum Rand hin stärker abfallend.
Dieser Verlauf ist in den 3 und 4 für zwei unterschiedliche Projektionswinkel dargestellt. In der 3 ist der Einfluss von Strahlaufhärtungskorrektur und Streustrahlungskorrektur bei einem Projektionswinkel von 0° für ein elliptisches Phantom dargestellt, wobei zu bemerken ist, dass der Strahlkegel, so wie es in der 2 gezeigt ist, bei einem Projektionswinkel von 0° das Phantom beziehungsweise den Patienten nicht vollständig umfasst.
Auf der Abszisse sind die Detektorpixel einer mittleren Bildzeile in willkürlichen Einheiten aufgetragen, während die Ordinate die negativen logarithmischen Projektionswerte, ebenfalls in willkürlichen Einheiten, wiedergibt. Die Kurve 11 zeigt den Verlauf fehlerfreier Projektionsdaten mit exakter Aufhärtungs- und Streustrahlungskorrektur. Mit der Kurve 12 sind die gemessenen Projektionswerte ohne Strahlaufhärtungs korrektur und ohne Streustrahlungskorrektur gezeigt, während die Kurve 13 die Projektionsdaten ohne Streustrahlanteil, jedoch mit Strahlaufhärtungskorrektur zeigt. Entsprechend sind unten die Differenzkurven dargestellt, die den Verlauf des Strahlaufhärtungseffektes 14 und den Verlauf des Streustrahlanteil 15 zeigen.
Die gleiche Darstellung ist in der 4, jedoch für einen Projektionswinkel von 90°, gezeigt, wobei hier am seitlichen Abfall der Kurven erkennbar ist, dass der Strahlenkegel das Phantom von der Seite vollständig umschließt. Die Bezugszeichen entsprechen den in der 3 verwendeten.
Da die stärkste Auswirkung auf das rekonstruierte Bild, das Cupping, durch die großen Projektionswerte, d. h. die Messstrahlen mit der stärksten Schwächung im Objekt, bestimmt. wird, kann ρ in Gleichung (#6a) durch einen geeigneten Mittelwert in der Umgebung von Pixeln in der Bildmitte beziehungsweise von Pixeln mit maximalem Projektionswert ersetzt werden.
Schreibt man für diese repräsentative Konstante ρ = ρ_o, dann wird die exakte Gleichung (#6) ersetzt durch die Näherung p* o = pT + ρo·δa (#7a)
Dies bedeutet also, dass die Streustrahlungskorrektur durch Verstärkung des Strahlaufhärtungs-Korrekturterms δa um einen Faktor ρ_o ersetzt werden kann und dadurch aufwendige Rechenschritte eingespart werden können.
Eine Genauigkeitsverbesserung ohne wesentliche Erhöhung der Rechenzeit kann erreicht werden, wenn anstelle eines globalen Verstärkungsfaktors für jede Zeile z im Projektionsbild jeweils ein spezifischer Verstärkungsfaktor ρ_o = ρ_o(z) eingeführt wird, der für jedes Pixel innerhalb der Zeile konstant bleibt.
Es hat sich gezeigt, dass dieses Korrekturverfahren nach Gleichung (#7a) bei einem typischen Objekt mit elliptischem Querschnitt zum Objektrand hin zu wenig abfällt und dort zu einer störenden Aufhellung führt, d. h. das Cupping wird zwar reduziert, aber es entsteht eine breite helle Ringzone am Objektrand.
Um diesem Effekt gegenzusteuern, kann man den Streustrahlungs-Korrekturterm δs_o = (ρ_o – 1)·δa in Gleichung (#7a) durch Multiplikation mit einer geeigneten Formfunktion h(y, z) modifizieren, deren Wertebereich zwischen 0 und 1 liegt und die zum Objektrand hin auf 0 abfällt. Dazu kann beispielsweise die normierte und eventuell auch geglättete Projektionsfunktion p_T selbst verwendet werden. In der Notation wird sie mit Überstrich dargestellt:
ρ1 = P0 – 1 ρ*1 (y, z) = pT(y, z) + (1 + p'1 h(y, z))δα(y, z) (#7b)
Dadurch erhöht sich zwar der Rechenaufwand wieder etwas und es geht Geschwindigkeit verloren, weil pixelweise im Vergleich zum Algorithmus (#7a) noch mindestens jeweils eine Addition, ein Speicherzugriff und eine Multiplikation dazukommen. Dieser modifizierte Algorithmus benötigt jedoch nach wie vor keine Delogarithmierung und Logarithmierung und ist daher immer noch schneller als der bekannte Algorithmus der Streustrahlungskorrektur, der auf den gemessenen Intensitäten und nicht den Projektionsdaten ausgeführt wird.
Eine weitere Verbesserung der Korrektur kann durch Reduktion der Dichteüberhöhung am Objektrand erreicht werden, beispielsweise durch eine Verallgemeinerung des Ansatzes in Gleichung (#7b), mit dem das Abfallen der Formfunktion dadurch verstärkt wird, dass man sie mit einem Exponenten N>=1 potenziert. Es gilt dann: hN(y, z) = [h(y, z)]N und P*1 (y, z) = PT(y, z) + (1 + ρ'1 hN(y, z))δα(y, z) (#7c)
Wie oben erwähnt wurde bei der Darstellung der Grundidee der Erfindung vereinfachend angenommen, dass zur Schätzung von ρ_o die Kenntnis von S und P_a, zumindest für repräsentative Pixel, bekannt sei. In den Gleichungen(#3) und (#6a) sieht man, dass nur das Verhältnis S/P_a eingeht. Dieses kann mit Hilfe der Datenbasis einer S/P-Tabelle geschätzt werden. Jedoch setzt die Suche in der S/P-Tabelle die Kenntnis der äquivalenten Wasserschichtdicke (Weglänge) voraus. Diese wird dann durch einen an sich bekannten Algorithmus im Rahmen der Wasserkorrektur aus dem logarithmierten Projektionswert P_a = –log(P_a) berechnet werden.
Die Gewinnung des in Gleichung (#6a) benötigten Quotienten S/P_a aus der S/P-Datenbasis unter Verwendung des genannten Algorithmus wird abkürzend mit S/Pa = g(pa) (#8)beschrieben.
Dieser Term hängt von p_a ab, p_a seinerseits ist aber erst nach der Schätzung des Streu-/Primär-Verhältnisses S/Pa bekannt. Dieses implizite Problem kann iterativ gelöst werden. Zur Vereinfachung der Bezeichnungen bedeuten nachfolgend p_a und p_T nur skalare Werte, da man sich zur Bestimmung von ρ_o auf einen Pixelmittelwert beschränken kann. Die Iterationsvorschrift für p_a lautet dann wie folgt: Start: pa (0) = pT (#9a) Iteration (n ≥ 0): pa (n+1) = PT + log(1 + g(pa (n))) (#9b)
In der Regel reichen zwei bis drei Iterationen aus.
Da das Korrekturverfahren direkt auf den Projektionsdaten operiert, ist es unabhängig von der Aufnahmegeometrie, jedoch muss diese bei der Erstellung der S/P-Datenbasis berücksichtigt werden, und prinzipiell kann jeder geeignete Rekonstruktionsalgorithmus verwendet werden.
Es wird nun nochmals auf die 3 und 4 hingewiesen, in denen mit dem Bezugszeichen 16 versehen die Ergebnisse der erfindungsgemäßen Streustrahlungskorrektur gezeigt sind. Es zeigt sich, dass die Unterschiede, insbesondere im 90° Projektionswinkel zum idealen Verlauf der Projektionsdaten sehr gering sind, so dass durch diese mit sehr wenig Rechenaufwand durchführbare Streustrahlungskorrektur kaum Qualitätseinbußen in der Bilddarstellung zu erwarten sind.
Konkret kann die Streustrahlungskorrektur als verstärkte Wasserkorrektur mit konstantem Verstärkungsfaktor je Projektionsbild somit entsprechend den folgenden Arbeitsschritten durchgeführt werden:
Am Anfang wird in Abhängigkeit von den für die Akquisition verwendeten Aufnahmeparametern, wie Feldgröße, Streustrahlenraster, Vorfilterung, etc. die entsprechende vorausberechnete mehrdimensionale, noch von der Röhrenspannung und dem Luftspalt abhängige S/P-Tabelle eingelesen.
Danach wird für jedes aus logarithmierten Daten bestehende Projektionsbild ein globaler Verstärkungsfaktor ρ wie folgt bestimmt:

1. zu Röhrenspannung und eventuell abzuschätzendem Luftspalt spezifische S/P-Teiltabelle extrahieren;
2. Projektionsbild reduzieren (z. B. durch 8×8 Downsampling oder Binning; Median- oder Tiefpassfilterung);
3. Pixel mit maximalem Projektionswert und Mittelwert aus z. B. 5×5 Pixeln einschließlich Maximum bestimmen;
4. verwende diesen Mittelwert für p_S und berechne mittels Iterationsalgorithmus gemäß den Gleichungen (#9a, b) eine Schätzung für p_a;
5. bestimme den Aufhärtungskorrekturterm δα = αW(E0)f–1U (pα) – pα
6. bestimme den Verstärkungsfaktor ρ entsprechend Gleichung (#7b), wobei der Term S/Pa im Sinne der Bedeutung von Gleichung (#8) mittels der S/P-Tabelle auszuwerten ist.

Anschließend wird für jedes Pixel je Projektionsbild die verstärkte Wasserkorrektur gemäß Gleichung (#6) mit dem vorher bestimmten konstanten Verstärkungsfaktor ρ durchgeführt.
Eine Variante dieser Vorgehensweise kann darin bestehen, dass statt dem Maximalwert im Schritt 3 der Mittelwert in einer so genannten „Dominante" verwendet wird. Danach wird im modifizierten Schritt 3 in einer vorgegebenen Dominante, das ist meist eine vorgegebene kleine ROI (= region of interest) in der Mitte des Detektors, der Mittelwert dieser reduzierten Region berechnet.
Erfindungsgemäß kann die Streustrahlungskorrektur auch als verstärkte Wasserkorrektur mit konstantem Verstärkungsfaktor je Bildzeile jedes Projektionsbildes stattfinden. Hier wird anstelle eines globalen Verstärkungsfaktors für jede Zeile z im Projektionsbild ein spezifischer Verstärkungsfaktor ρ = ρ(z) eingeführt, der aber für jedes Pixel innerhalb der Zeile konstant bleibt. Durch diese Modifikation des Verfahrens, die den Rechenaufwand zwar etwas erhöht, wird eine Genauigkeitsverbesserung erreicht. Der Programmablauf je Projektionsbild ändert sich dann wie folgt:
Zunächst wird für jedes Projektionsbild aus logarithmierten Daten zur allgemeinen Vorbereitung:

• zur Röhrenspannung und eventuell zum abzuschätzenden Luftspalt eine spezifische S/P-Teiltabelle extrahiert;
• das Projektionsbild reduziert, z. B. durch 4×4 Downsampling oder Binning und Median- oder Tiefpassfilterung;

Es folgt nun die Bestimmung des zeilenabhängigen Verstärkungsfaktors ρ(z) mit:
Bestimmung für jede Zeile z' im reduzierten Projektionsbild je eines Verstärkungsfaktors ρ(z') mit den Schritten:

• Bestimmung des Pixels mit maximalem Projektionswert und Mittelwert aus mehreren Pixeln einschließlich Maximum bestimmen; (alternativanalog zur Variante 5.2: Mittelwert über z. B. 5 cm in der Umgebung der Zeilenmitte)
• verwende diesen Mittelwert für p_S und berechne mittels Iterationsalgorithmus gemäß Gin. (# 9a, b) eine Schätzung für p_a;
• bestimme den Aufhärtungskorrekturterm δa = δa(z') gemäß δα = αW(E0)f–1U (pα) – Pα
• bestimme den Verstärkungsfaktor ρ = ρ(z') entsprechend Gl. (#7b), wobei der Term S/Pa im Sinne der Bedeutung von Gl. (# 8) mittels der S/P-Tabelle auszuwerten ist.

Es erfolgt nun die Expansion des auf dem reduzierten Zeilenraster bestimmten Verstärkungsfaktors ρ(z') auf das ursprüngliche feine Raster durch lineare Interpolation; vorher kann eine Glättung durch 1-dimensionale Faltung durchgeführt werden.
Schließlich wird für jede Projektionsbildzeile z die Wasserkorrektur mit dem Verstärkungsfaktor ρ(z) durchgeführt.
Insgesamt wird also durch die Erfindung eine Streustrahlungskorrektur zur Verfügung gestellt, die auf die Aufhärtungskorrektur aufgesetzt ist und auf den logarithmierten Projektionsdaten stattfindet. Es ist somit nicht mehr notwendig, die als Projektionen vorliegenden CT-Daten vor und zur Streustrahlungskorrektur zu delogarithmieren und nach erfolgter Korrektur wieder zu logarithmieren, so dass ein nicht unerheblicher Rechenaufwand eingespart werden kann.

Claims

Verfahren zur Beschleunigung der Streustrahlungskorrektur in einem Computertomographiesystem (1) mit einer Strahlungsquelle (2) und einem großflächig ausgebildeten Detektor (3) mit einer Vielzahl von Detektorzeilen, wobei ein Objekt (7) aus einer Vielzahl von Projektionswinkeln abgetastet wird, aus den gemessenen Schwächungswerten der Strahlung Projektionsdaten erstellt werden, die zur Rekonstruktion tomographischer Darstellungen weiterverarbeitet werden, wobei eine Strahlaufhärtungskorrektur direkt auf den Projektionsdaten p_T vorgenommen wird, dadurch gekennzeichnet, dass die Streustrahlungskorrektur ebenfalls direkt auf den Projektionsdaten p_T erfolgt und hierfür ein in der Strahlaufhärtungskorrektur verwendeter Aufhärtungskorrekturterm mit einem Verstärkungsfaktor multipliziert wird.
Verfahren gemäß dem voranstehenden Patentanspruch 1, dadurch gekennzeichnet, dass für die kombinierte Aufhärtungs- und Streustrahlungskorrektur korrigierte Projektionsdaten p^* = p_T + ρ·δ_a mit dem Verstärkungsfaktor ρ = 1 + 1/δ_a log(1 + S/P_a) berechnet werden, wobei δ_a den Korrekturterm für die Strahlaufhärtungskorrektur, S = I_S/I_o die normierte Streustrahlungsverteilung und P_a = I_A/I_o die normierte Primärstrahlung mit Aufhärtung ohne Streustrahlung beschreiben.
Verfahren gemäß dem voranstehenden Patentanspruch 2, dadurch gekennzeichnet, dass das Verhältnis S/P_a von normierter Streustrahlung S zur normierten Primärstrahlung mit Aufhärtung ohne Streustrahlung P_a unter Berücksichtigung mindestens eines repräsentativen Pixels geschätzt wird.
Verfahren gemäß dem voranstehenden Patentanspruch 2, dadurch gekennzeichnet, dass das Verhältnis S/P_a von normierter Streustrahlung S zur normierten Primärstrahlung mit Aufhärtung ohne Streustrahlung P_a aus Simulati onsdaten eines Phantoms mit zum untersuchten Objekt ähnlicher Form, Größe und Struktur entnommen wird.
Verfahren gemäß dem voranstehenden Patentanspruch 2, dadurch gekennzeichnet, dass das Verhältnis S/P_a von normierter Streustrahlung S zur normierten Primärstrahlung mit Aufhärtung ohne Streustrahlung P_a durch ein iteratives Verfahren angenähert wird.
Verfahren gemäß einem der voranstehenden Patentansprüche 1 oder 5, dadurch gekennzeichnet, dass der Verstärkungsfaktor ortsabhängig relativ zum Detektor ist.
Verfahren gemäß einem der voranstehenden Patentansprüche 1 bis 6, dadurch gekennzeichnet, dass der Verstärkungsfaktor abhängig von der betrachteten Detektorzeile ist.
Verfahren gemäß einem der voranstehenden Patentansprüche 1 oder 7, dadurch gekennzeichnet, dass der Verstärkungsfaktor zusätzlich projektionswinkelabhängig ist.
Verfahren gemäß einem der voranstehenden Patentansprüche 1 bis 8, dadurch gekennzeichnet, dass der Verstärkungsfaktor nochmals mit einer Funktion multipliziert wird, deren Wertebereich zwischen 0 und 1 liegt und deren Wert zum Objektrand auf 0 abfällt.
Verfahren gemäß dem voranstehenden Patentanspruch 9, dadurch gekennzeichnet, dass die Funktion mit einem Exponenten potenziert wird.
CT-System (1) mit: 11.1 einer Strahlungsquelle (2) und 11.2 einem großflächig ausgebildeten Detektor (3) mit einer Vielzahl von Detektorzeilen, mit dem ein Objekt (7) aus einer Vielzahl von Projektionswinkeln abgetastet wird, wobei aus den gemessenen Schwächungswerten der Strahlung Projektionsdaten erstellt werden, und 11.3 einer Steuer- und Recheneinheit (9), die mindestens einen Speicher für Programme (Prg_x) und Daten aufweist, wobei mit den gespeicherten Programmen und Projektionsdaten eine Rekonstruktion tomographischer Darstellungen durchgeführt wird, dadurch gekennzeichnet, dass 11.4 Computerprogramme (Prg_x) oder Programm-Module gespeichert sind, die im Betrieb die Verfahrensschritte eines der voranstehenden Verfahrensansprüche durchführen.