DE102007020065A1

DE102007020065A1 - Verfahren für die Erstellung von Massenbelegungsbildern anhand von in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbildern

Info

Publication number: DE102007020065A1
Application number: DE102007020065A
Authority: DE
Inventors: Philipp Dr. Bernhardt; Ernst-Peter Rührnschopf
Original assignee: Siemens AG
Current assignee: Siemens Healthcare GmbH
Priority date: 2007-04-27
Filing date: 2007-04-27
Publication date: 2008-10-30
Also published as: JP5264268B2; US20100027867A1; US8041096B2; JP2008272476A

Abstract

Für die Erstellung von streustrahlungskorrigierten Massenbelegungsbildern in der dualen Röntgenabsorptiometrie wird vorgeschlagen, die von den Schwächungsbildern in unterschiedlichen Energiebereichen gelieferte zusätzliche Information in einen inhomogenen Korrekturbildbereich (24) dazu zu verwenden, durch Inversion von Primärstrahlungsfunktion eine für eine Vielzahl von Schwächungsbildern übereinstimmende mehrdimensionale Massenbelegung zu finden und anhand der mehrdimensionalen Massenbelegung den Streustrahlungsanteil zu ermitteln.

Description

Die Erfindung betrifft ein Verfahren für die Erstellung von Massenbelegungsbildern anhand von in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbildern mit den Verfahrensschritten:

– Erzeugen von Strahlung mit einer Strahlungsquelle und Durchleuchten eines zu untersuchenden Objekts mit Hilfe der Strahlung;
– Beaufschlagen einer Detektorvorrichtung mit der Strahlung und Erfassen von Schwächungsbildern in unterschiedlichen Energiebereichen durch die Detektorvorrichtung; und
– Bestimmen von Massenbelegungsbildern durch eine der Detektorvorrichtung nachgeschaltete Auswerteeinheit.

Ein derartiges Verfahren ist aus WARP, R. J.; DOBBINS, J. T.: Quantitative evaluation of noise reduction strategies in dual-energy imaging In: Med. Phys. 30 (2), Seiten 190 bis 198, Februar 2003, bekannt. Im Rahmen der dualen Röntgenabsorptiometrie wird das zu untersuchende Objekt, bei dem es sich in der Regel um einen Patienten handelt, mit Röntgenstrahlung in unterschiedlichen Energiebereichen durchleuchtet. Die duale Röntgenabsorptiometrie kann dabei mit einer einzelnen Aufnahme oder mit einer Reihe von hintereinander durchgeführten Aufnahmen bewerkstelligt werden.
Im ersten Fall wird ein Doppeldetektor mit zwei unterschiedlichen Szintillationsmaterialien verwendet, deren Ansprechverhalten möglichst weit auseinander liegende Energieschwerpunkte hat. Im zweiten Fall werden aufeinander folgende Aufnahmen mit möglichst verschiedenen Röntgenspektren gemacht, die bei Verwendung von Röntgenröhren durch Veränderung der Röhrenspannung, durch die die Elektronen beschleunigt werden, oder durch die Auswahl von Vorfiltern erzeugt werden können.
Für jeden Bildpunkt der aufgenommenen Projektionsbilder kann aus dem Schwächungsverhalten in den unterschiedlichen Energiebereichen auf die Materialkomposition im Strahlengang zwischen der punktförmigen Röntgenquelle und den Bildpunkten geschlossen werden. Die Projektionsbilder werden nachfolgend auch als Schwächungsbilder bezeichnet. Ferner ist unter Materialkomposition die Massenbelegung der verschiedenen Materialien entlang des Strahls durch das zu untersuchende Objekt zu verstehen.
Durch eine Kombination der in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbilder lassen sich Massenbelegungsbilder erstellen, die zumindest näherungsweise die Massenbelegung der im untersuchten Objekt enthaltenen verschiedenen Materialien wiedergeben. Üblicherweise werden die Schwächungsbilder linear kombiniert, wobei die Gewichtungsfaktoren empirisch ermittelt werden. Eine mathematisch exakte Bestimmung der Massenbelegungsdichte ist jedoch mit den bekannten Verfahren kaum möglich.
Die mathematisch exakte Bestimmung der Massenbelegung wird im Übrigen auch durch die vorhandene Streustrahlung erschwert. Bereits in der Projektionsradiographie mit Flächendetektoren spielt Streustrahlung wegen des großen erfassten Raumwinkels eine erhebliche Rolle. Um die Streustrahlung zu reduzieren, werden häufig Streustrahlenraster (= anti-scatter grids) unmittelbar über der Detektoreingangsfläche eingesetzt. Die duale Röntgenabsorptiometrie stellt als quantitatives Verfahren nochmals höhere Anforderungen an die Genauigkeit der Messdaten als die einfache Projektionsbildgebung im Rahmen der Projektionsradiographie. Trotz Streustrahlenraster kann der die Daten verfälschende Streustrahlungsanteil erheblich sein. Zum Beispiel wird im Thoraxbereich gewöhnlich mit sehr kleinem Luftspalt gearbeitet, was zur Folge hat, dass trotz Streustrahlenraster die Streustrahlungsintensität noch die Primärintensität überwiegen kann, vor allem in Bildregionen mit starker Schwächung und bei höheren Photonenenergien, entspre chend Röntgenröhren-Spannungen oberhalb von 100 kV. Außerdem ist es eine empirische Tatsache, dass die Streustrahlungsanteile bei den höher- und niederenergetischen Bilddaten sehr unterschiedlich sind. Dies gilt insbesondere bei kleinem Luftspalt, also in den Fällen, in denen der Abstand zwischen dem Streuobjekt und dem Detektor gering ist. Insgesamt kann das Vorhandensein von Streustrahlung, trotz Streustrahlenraster, bei der dualen Röntgenabsorptiometrie zu unzuverlässigen und teilweise unbrauchbaren Ergebnissen, zum Beispiel zu negativen Materialdicken, führen. Die Korrektur der Streustrahlung ist daher bei der dualen Röntgenabsorptiometrie von großer Bedeutung.
Daher sind bei der dualen Röntgenabsorptiometrie zusätzlich zur Verwendung von Streustrahlenrastern rechnerische Streustrahlungskorrekturverfahren erforderlich.
Es sei an dieser Stelle angemerkt, dass die Streustrahlung nachfolgend auch als Sekundärstrahlung bezeichnet wird. Die Summe von Primärstrahlung und Sekundärstrahlung, die die gemessenen Bildwerte ergibt, wird als Gesamtstrahlung bezeichnet.
Aus HINSHAW, D. A.; DOBBINS III, J. T.: Recent Progress in noise reduction and scatter correction in dual-energy imaging. In: Proc. SPIE, 1995, Vol. 2432, Seiten 134 bis 142, ist ein Verfahren zur Streustrahlungskorrektur im Rahme der dualen Röntgenabsorptiometrie bekannt. Bei dem bekannten Verfahren wird für jedes der in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbilder eine Streustrahlungskorrektur durchgeführt, indem für einen gegebenen Bildpunkt ein empirisch ermittelter Streustrahlungsanteil in Abhängigkeit vom Bildwert bestimmt wird. Der Streustrahlungsanteil bestimmt die Form und die Breite einer Verteilungsfunktion für die Streustrahlung. Anhand der Streustrahlungsfunktion werden die Streubeiträge in benachbarten Pixeln berechnet. Das Verfahren wird anschließend für weitere Bildwerte wiederholt und die Streubeiträge in den einzelnen Pixeln aufsummiert. Somit fin det eine Faltung des mit der Detektorvorrichtung aufgezeichneten Bildes mit einer Verteilungsfunktion statt, deren Breite und Form von den Bildwerten des von der Detektorvorrichtung aufgenommenen Schwächungsbildes abhängen.
Aus FLOYD, C. B.; BAKER, J. A.; LO, J. Y.; RAVIN, C. E: Posterior Beam-Stop Method for Scatter Fraction Measurement in Digital Radiography. In: Investigative Radiology Feb. 1992, Vol. 27, Seiten 119 bis 123, ist ein messtechnisches Verfahren zur Bestimmung der Streustrahlung gemäß der Beamstop-Methode bekannt. Dieses Verfahren eignet sich für Anwendungen im Labor mit Phantomen, kaum aber für den klinischen Betrieb.
Ferner sind aus ZELLERHOFF, M.; SCHOLZ, B.; RÜHRNSCHOPF, E.-P.; BRUNNER, T.: Low contrast 3D-reconstruction from C-arm data. In: Proceedings of SPIE. Medical Imaging, 2005, Vol. 5745, Seiten 646 bis 655 verschiedene rechnerische Verfahren für die Streustrahlungskorrektur im Rahmen der Computertomographie bekannt.
Die bekannten rechnerischen Verfahren sind in der Regel jedoch ziemlich komplex und aufwendig.
Es besteht daher weiterhin ein Bedarf an vergleichsweise einfachen Korrekturverfahren, mit denen sich deutliche Bildqualitätsverbesserungen erzielen lassen.
Ausgehend von diesem Stand der Technik liegt der Erfindung daher die Aufgabe zugrunde, ein Verfahren anzugeben, mit dem sich Massenbelegungsbilder mit im Vergleich zum Stand der Technik verbesserter Qualität erstellen lassen.
Diese Aufgabe wird durch ein Verfahren mit den Merkmalen des unabhängigen Anspruchs gelöst. In davon abhängigen Ansprüchen sind vorteilhafte Ausgestaltungen und Weiterbildungen angegeben.
Bei dem Verfahren werden die Massenbelegungswerte für die Massenbelegung im Objekt bestimmt, indem eine mehrdimensionale Schwächungsfunktion invertiert wird. Da die Schwächungsfunktion eine mehrdimensionale Massenbelegung im Objekt mit Schwächungswerten der in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbilder verknüpft, kann durch eine Inversion der mehrdimensionalen Schwächungsfunktion aus den Schwächungsbildern das zugehörige Massenbelegungsbild ermittelt werden. Da die Schwächungswerte monoton von den Massenbelegungswerten abhängen und die Schwächungswerte nicht linear mit der Energie der verwendeten Strahlung skalieren, existiert immer eine eindeutige Lösung für das Inversionsproblem. Insofern kann die Massenbelegung der im Objekt vorliegenden Komponenten exakt bestimmt werden. Letztlich liegt dem Verfahren eine Konsistenzbedingung zugrunde, nämlich dass die Inversion der mehrdimensionalen Schwächungsfunktion notwendigerweise übereinstimmende Massenbelegungen ergeben muss, da die in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbilder jeweils das gleiche Objekt abbilden.
Bei einer bevorzugten Ausführungsform des Verfahrens wird ein durch Streuung hervorgerufener Sekundärstrahlungsanteil bestimmt und die Schwächungsbilder werden hinsichtlich des Sekundärstrahlungsanteils auf einen durch Schwächung erzeugten Primärstrahlungsanteil korrigiert. Dadurch kann die Qualität der durch Inversion der mehrdimensionalen Schwächungsfunktion erzeugten Massenbelegungsbilder wesentlich verbessert werden, da die Ergebnisse nicht durch die Sekundärstrahlungsanteile verfälscht werden.
Bei einer weiteren bevorzugten Ausführungsform des Verfahrens werden zur Bestimmung der Sekundärstrahlungsanteile der in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbilder diejenigen Sekundärstrahlungsanteile gesucht, die mit Primärstrahlungsanteilen verknüpft sind, aus denen sich beim Auswerten der inversen mehrdimensionalen Schwächungsfunktion jeweils die gleiche Massenbelegung für den Korrekturbildbereich ergibt. Da die Sekundärstrahlungsanteile von der Mas senbelegung abhängen und die Massenbelegung wiederum von den Primärstrahlungsanteilen abhängen, die sich erst dann ergeben, wenn die gemessenen Bildwerte hinsichtlich der Sekundärstrahlungsanteile korrigiert werden, ergeben sich implizite Gleichungen, deren Lösung durch einen Suchvorgang ermittelt werden muss. Üblicherweise wird der Suchvorgang mithilfe einer Iteration ausgeführt, durch die sich die impliziten Gleichungen lösen lassen.
Bei einer weiteren bevorzugten Ausführungsform des Verfahrens werden die Sekundärstrahlungsanteile in einem Korrekturbereich bestimmt, der einen inhomogenen Objektbereich abbildet. Dabei wird eine mehrdimensionale Schwächungsfunktion invertiert, die von mehrdimensionalen Massenbelegungen abhängt. Das Ergebnis der Inversion ist eine mehrdimensionale Massenbelegung, mit der die Schwächungsstruktur des zu untersuchenden Objekts zumindest näherungsweise bestimmt werden kann. Anhand der mehrdimensionalen Massenbelegung kann dann der Sekundärstrahlungsanteil zumindest annähernd bestimmt werden.
Vorzugsweise wird der Sekundärstrahlungsanteil ermittelt, indem eine von der mehrdimensionalen Massenbelegung abhängige Faltungsfunktion mit dem ermittelten Primärstrahlungsanteil in den Schwächungsbildern gefaltet wird.
Daneben ist es möglich, mit Hilfe von Monte-Carlo-Simulation den jeweiligen Sekundärstrahlungsanteil zu bestimmen. Die Monte-Carlo-Simulation können vorab durchgeführt und die Ergebnisse in Abhängigkeit von mehrdimensionalen Massenbelegungen tabelliert werden.
Die Streustrahlungskorrektur kann jeweils an einzelnen Bildpunkten des Schwächungsbildes durchgeführt oder anhand jeweils über einen vorbestimmten Bereich gemittelter Bildwerte bestimmt werden. Vorzugsweise wird der Sekundärstrahlungsanteil jeweils im Bereich von Stützpunkten eines über das Schwächungsbild gelegten Rasters bestimmt und der Sekundär strahlungsanteil für Bildpunkte zwischen den Stützpunkten interpoliert.
Weitere Eigenschaften und Vorteile der Erfindung gehen aus der nachfolgenden Beschreibung hervor, in der Ausführungsbeispiele der Erfindung anhand der beigefügten Zeichnung im Einzelnen erläutert werden. Es zeigen:
1 eine Vorrichtung für die duale Röntgenabsorptiometrie;
2 zwei bei verschiedenen Röhrenspannungen aufgenommene Photonenspektren einer Röntgenröhre mit einer Anode aus Wolfram;
3 den Verlauf des Massenschwächungskoeffizienten in Abhängigkeit von der Photonenenergie für verschiedene Körperbestandteile;
4 ein Schwächungsbild, mit einem inhomogenen Bildbereich und einem homogenen Bildbereich;
5 ein Diagramm, das die Abhängigkeit der Streustrahlungsintensität von verschiedenen Materialkonstellationen veranschaulicht;
6 ein Diagramm, in das das Streu-zu-Primär-Verhältnis in Abhängigkeit von der Wasserdicke für verschiedene Röhrenspannungen eingetragen ist; und
7 ein Ablaufdiagramm des Verfahrens zur Streustrahlungskorrektur;
8 einen Querschnitt durch ein virtuelles Phantom zur Veranschaulichung der Rekonstruktion von Massenbelegungsbilder;
9 eine Aufsicht auf das virtuelle Phantom aus 8;
10 eine Simulation einer niederenergetischen Röntgenaufnahme des Phantoms aus den 8 und 9;
11 eine Simulation einer hochenergetischen Röntgenaufnahme des Phantoms aus den 8 und 9;
12 ein anhand der Röntgenbilder aus den 9 und 10 rekonstruiertes Weichteilbild des virtuellen Phantoms aus den 8 und 9; und
13 ein anhand der Röntgenbilder aus den 9 und 10 rekonstruiertes Knochenbild des virtuellen Phantoms aus den 8 und 9.
1 zeigt eine Röntgenanlage 1, mit der Röntgenaufnahmen für die duale Röntgenabsorptiometrie durchgeführt werden können. Die Röntgenanlage 1 umfasst eine Röntgenröhre 2, die eine von einem Glühfaden gebildete Kathode 3 aufweist. Der Glühfaden 3 kann mit Hilfe eines Heizstromes I beheizt werden. Dabei werden von der Kathode 3 Elektronen emittiert, die mit Hilfe einer Röhrenspannung U in Richtung einer Anode 4 beschleunigt werden. Dadurch entsteht ein Elektronenstrahl 5, der in einem Brennfleck auf die Anode 4 trifft. Die in der Anode 4 abgebremsten Elektronen erzeugen Röntgenstrahlung 6, die zur Unterdrückung des niederenergetischen Teils zunächst einen Vorfilter 7 durchläuft, der die Wirkung eines spektralen Filters hat. Bei den Vorfiltern 7 handelt es sich in der Regel um dünne Kupferplatten, die in unterschiedlicher Dicke in den Strahlengang der Röntgenstrahlung 6 eingebracht werden können. Anschließend durchdringt die Röntgenstrahlung 6 einen zu untersuchenden Patienten 8.
Die durch den Patienten 8 hindurchgetretene Röntgenstrahlung 6 wird von einem Röntgendetektor 9 erfasst, der ein Schwächungsbild des Patienten 8 aufzeichnet. Dabei wird die Struk tur des die Röntgenstrahlung 6 schwächenden Materials im Patienten 8 auf den Röntgendetektor 9 projiziert. Die Schwächungsbilder enthaltenden Röntgenaufnahmen werden daher auch als Projektionsbilder bezeichnet.
Bei dem Röntgendetektor 9 handelt es sich vorzugsweise um einen Flachbilddetektor oder Flächendetektor auf Halbleiterbasis, der eine Vielzahl von Detektorelementen aufweist, mit denen sich ein digitales Röntgenbild erzeugen lässt. Die Detektorelemente nehmen jeweils einen Bildpunkt auf und werden auch als Pixel bezeichnet.
Dem Röntgendetektor 9 ist eine Auswerteeinheit 10 nachgeschaltet, die in der Regel eine Linearkombination der durch Variation der Röhrenspannung U und der Vorfilter 7 in unterschiedlichen Energiebereichen der Röntgenstrahlung 6 aufgenommenen Schwächungsbilder bildet. Das durch die Linearkombination der in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbilder erzeugte Kombinationsbild wird auf einer Anzeigeeinheit 11 dargestellt.
Bei der Linearkombination der Schwächungsbilder kann es sich beispielsweise um die Bildung einer Differenz handeln, durch die die Knochenstruktur des Patienten 8 aus dem Kombinationsbild eliminiert wird. Das auf diese Weise erzeugte Kombinationsbild enthält die Schwächungsstruktur des Weichteilgewebes, was insbesondere bei Untersuchungen der Lunge von Vorteil ist.
Für die Durchführung von Schwächungsbildern in unterschiedlichen Energiebereichen werden insbesondere die Röhrenspannung U und die Vorfilter 7 variiert. Für das Schwächungsbild im Niederenergiebereich kann zum Beispiel eine niedrigere Röhrenspannung U verwendet werden. Außerdem können die Vorfilter 7 eine geringe Materialstärke aufweisen, so dass der niederenergetische Teil des von der Röntgenröhre 2 erzeugten Spektrums nur geringfügig unterdrückt wird. Für die Schwächungsbilder im Hochenergiebereich kann dagegen eine hohe Röhren spannung U verwendet werden. Außerdem können Vorfilter 7 mit größerer Materialstärke verwendet werden, die vorzugsweise den hochenergetischen Teil des von der Röntgenröhre 2 erzeugten Röntgenspektrums passieren lassen.
Im Folgenden werden Verfahren beschrieben, mit denen sich die Massenbelegungen anhand der Schwächungsbilder weitgehend exakt bestimmen lassen.
Da die Streustrahlungskorrektur für die Brauchbarkeit der Kombinationsbilder eine wesentliche Rolle spielt, werden im Folgenden insbesondere auch Verfahren beschrieben, durch die sich die Streustrahlungseffekte reduzieren lassen.
1. Problemstellung
Die Signalbildung der Röntgenstrahlung 6, die den Patient 8 durchdringt, wird wesentlich durch das Emissionsspektrum Q_U(E), nämlich das von der angelegten Röhrenspannung U abhängige Energiespektrum der als Bremsstrahlung an der Anode emittierten Photonen, sowie durch die Transparenz T_F(E) der verwendeten spektralen Filter und die spektrale Ansprechempfindlichkeit η_D(E) des Röntgendetektors 9 bestimmt.
Die resultierenden effektiven normierten Spektralverteilungen W(E; U) sind definiert durch: W(E; U) = QU(E)TF(E)ηD(E)/cU wobei der Faktor c_U die integrierte effektive normierte Spektralverteilung auf den Wert = 1 normiert.
Beispiele für zwei effektive Spektralverteilungen, die den Röhrenspannungen 60 kV und 120 kV entsprechen, sind in 2 dargestellt. In 2 ist dabei die relative Photonenhäufigkeit N_rel/N_T je 1-keV-Intervall gegen die Photonenenergie E in keV aufgetragen, wobei N_T die Gesamtzahl der Photonen ist. Ein Röntgenspektrum 12 ist dabei einer Röhrenspannung von 70 kV zugeordnet und ein Röntgenspektrum 13 einer Röhrenspannung von 120 kV.
Die Abhängigkeit des Massenschwächungskoeffizienten für Wasser (μ/ρ)(E) von der Photonenenergie E ist in 3 dargestellt. Eine Massenschwächungskurve 14 für Wasser ist dabei in etwa deckungsgleich mit einer Massenschwächungskurve für Blut, obwohl Blut eine größere Dichte ρ als Wasser hat. Fettgewebe weist dagegen eine von der Massenschwächungskurve 14 geringfügig abweichende Massenschwächungskurve 15 auf. Eine weitere Massenschwächungskurve 16 gibt den Verlauf des Massenschwächungskoeffizienten von Knochengewebe an. Weitere Massenschwächungskurven 17 und 18 beschreiben den Verlauf des Massenschwächungskoeffizienten von Calcium und von Jod, das bei einer Photonenenergie von 33,2 keV eine K-Kante 19 aufweist. Jod wird häufig als Kontrastmittel verwendet.
Anhand von 3 wird deutlich, dass Knochenmaterial Röntgenstrahlung stärker absorbiert als Weichteilgewebe. Der Schwächungskoeffizient der Röntgenstrahlung fällt jedoch bei Knochengewebe zu höheren Energien stärker ab als die Absorption von Weichteilgewebe. Daraus ergibt sich, dass die Primärstrahlungsfunktion auch bei der Verwendung monochromatischer Röntgenstrahlung nicht mit der Energie der Röntgenstrahlung skaliert. Anhand der in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbilder können daher Massenbelegungsbilder für Knochengewebe und Weichteilgewebe ermittelt werden. Ebenso lassen sich aufgrund der unterschiedlichen Energieabhängigkeit der Massenschwächungskurve 14 für Wasser und der Massenschwächungskurve 15 für Fettgewebe Wasser und Fettgewebe separieren.
In 4 ist beispielsweise ein Schwächungsbild 20 mit einem homogenen Bildbereich 21 und einem inhomogenen Bildbereich 22 gezeigt. Im homogenen Bildbereich 21 wird ausschließlich Gewebe mit einheitlicher Energieabhängigkeit des Schwächungskoeffizienten, zum Beispiel Weichteilgewebe, abge bildet, während im inhomogenen Bildbereich 22 unterschiedliche Arten von Gewebe mit unterschiedlicher Energieabhängigkeit des Schwächungskoeffizienten, beispielsweise sowohl Weichteilgewebe als auch Knochengewebe, abgebildet werden. Anhand eines Schwächungsbildes, das sowohl den homogenen Bildbereich 21 als auch den inhomogenen Bildbereich 22 abdeckt, können nun Massenbelegungsbilder für das Weichteilgewebe und das Knochengewebe erstellt werden. Dabei ist es nicht unbedingt erforderlich, dass die Streustrahlungskorrektur über das gesamte Schwächungsbild hinweg ausgeführt wird. Die Streustrahlungskorrektur kann beispielsweise in einem Korrekturbereich 23 ausgeführt werden, der sich im homogenen Bildbereich 21 befindet, oder in einem Korrekturbereich 24 durchgeführt werden, der im inhomogenen Bildbereich 22 angeordnet ist.
In dem Korrekturbereich 24, der einem inhomogenen Objektbereich zugeordnet ist, ist die Korrektur der Streustrahlung jedoch problematisch. Denn selbst bei gleicher Schwächung der Röntgenstrahlung, das heißt, bei gleicher normierter Primärintensität P₁ oder P₂, hängt die Streustrahlung von der Materialzusammensetzung im Strahlengang ab, also von den Massenbelegungen b₁ und b₂ der beiden Materialien. Ein Beispiel hierfür ist in 5 dargestellt.
In 5 sind die Ergebnisse von Monte-Carlo-Simulationen dargestellt, die die Abhängigkeit der Streustrahlung von der Materialzusammensetzung bei gleicher Primärstrahlung zeigen. Als Röhrenspannung wurde eine Spannung von U = 70 kV angenommen. Die Breite des Luftspalts wurde auf 20 cm gesetzt. Die Fläche des Röntgendetektors 9 wurde mit 40 × 30 cm² angesetzt. Die Simulation wurde so lange durchgeführt, bis sich jeweils die gleiche Primärstrahlungskurve 25 ergab. Wenn die ungeschwächte Strahlung zunächst eine Scheibe aus Knochenmaterial mit einer Dicke von 5 cm und anschließend Wasser einer Dicke von 20 cm durchläuft, ergibt sich eine Streustrahlungsverteilung 26. Wenn dagegen eine knochenfreie homogene Wasserscheibe mit 25 cm Dicke durchlaufen wird, ergibt sich die Streustrahlungsverteilung 27. Die Primärstrahlungskurve 25 sowie die Streustrahlungsverteilung 26 und 27 sind dabei in willkürlichen Einheiten gegen die Koordinaten x_d der zentralen Detektorzeile aufgetragen. Die Streustrahlungsverteilungen 26 und 27 sind jeweils mit der gleichen Skala dargestellt.
Anhand 5 wird deutlich, dass die Streustrahlungsverteilung von der vorhandenen Materialkomposition abhängt. Um die gemessenen Projektionsbilder hinsichtlich der Streustrahlung korrigieren zu können, ist es daher notwendig, die Massenbelegung der einzelnen Materialien zu kennen.
Im Folgenden werden zwei Messwerte für dasselbe Objekt betrachtet, aufgenommen mit zwei verschiedenen Röhrenspannungen U₁ und U₂, eventuell auch verschiedenen röhrenseitigen oder detektorseitigen Filtern. Die effektiven Spektralverteilungen bezeichnen wir mit: W1(E) = W(E; U1), W2(E) = W(E; U2)
Es wird vorausgesetzt, dass sämtliche bei digitalen Röntgendetektoren 9 übliche Kalibrierungskorrekturen durchgeführt sind. Im Einzelnen sind das eine Dunkelbildsubtraktion und eine Korrektur der unterschiedlichen Empfindlichkeiten der Pixel. Außerdem wird vorausgesetzt, dass eine korrekte I₀-Normierung vorgenommen worden ist. Unter einer I₀-Normierung ist dabei Folgendes zu verstehen: Die volle ungeschwächte Intensität ohne schwächendes Objekt wird für beide Spektren bestimmt und jeder Intensitätsmesswert an jedem Pixel des Röntgendetektors 9 wird mittels Division durch den entsprechenden I₀-Wert normiert. Die folgenden Ausführungen beziehen sich immer auf normierte Intensitätswerte, seien es nun Primärstrahlungs-, Streustrahlungs- oder gemessene Gesamtstrahlungsintensitätswerte, die aus der Summe von Primärstrahlung und Streustrahlung bestehen.
In den nachfolgenden Ausführungen werden folgende Bezeichnungen verwendet:

I₁, I₂: normierte gemessene Gesamtintensitäten beim Spektrum W₁(E) oder W₂(E),
P₁, P₂: normierte, zunächst unbekannte Primärintensitäten,
S₁, S₂: normierte, zunächst unbekannte Streustrahlungsintensitäten.

Theoretisch gilt für die normierten Primärintensitäten (= Primärschwächungen) bei den Spektren W₁(E) und W₂(E) im Falle von homogenem Gewebe, wenn der jeweilige Röntgenmessstrahl dieselbe effektive Weglänge X [cm] oder Massenbelegung X [g/cm²] durchläuft,
Es sei darauf hingewiesen, dass in diesen Gleichungen μ im Exponenten als Massenschwächungskoeffizient zu verstehen ist, wenn X als Massenbelegung interpretiert wird. Wenn X dagegen als Wegstrecke aufgefasst wird, ist unter μ der lineare Schwächungskoeffizient zu verstehen.
Im Rahmen der dualen Röntgenabsorptiometrie wird die qualitative Trennung von zwei Materialien üblicherweise nach der folgenden empirischen Formel vorgenommen: b ~1 = –log(I1) + w1log(I2) b ~2 = log(I1) – w2log(I2)
Hier sind die Größen auf der linken Seite Grauwertbilder, wobei die empirischen Faktoren w₁, w₂ so zu wählen sind, dass visuell eine möglichst gute Separierung der beiden Materialien erzielt wird. I₁ und I₂ auf der rechten Seite sind die gemessenen Intensitäten im höher- und niederenergetischen Intensitätsbild. Im Allgemeinen enthalten I₁ und I₂ unterschiedliche Anteile von Streustrahlung. Somit ist die Wahl der empirischen Faktoren w₁, w₂ auch abhängig von den unterschiedlichen Streustrahlungsanteilen. Die Grauwertbilder auf der linken Seite stellen zunächst keine physikalischen Größen dar. Die durch die Linearkombination erhaltenen Grauwertbilder können erst mittels geeigneter Skalierungsfaktoren und additiver Konstanten näherungsweise als Pseudo-Materialdicken der beiden zu selektierenden Materialien, zum Beispiel als Materialdicken von Weichteilgewebe und Knochen interpretiert werden.
Im Folgenden wird daher ein Verfahren beschrieben, mit dem aus den Schwächungsbildern weitestgehend exakt nach den einzelnen Komponenten separierte Massenbelegungsbilder erstellt werden können, wobei Streustrahlungseffekte eliminiert werden.
2. Allgemeine Prozeduren einer Streustrahlungskorrektur in inhomogenen Bildbereichen
Bei den nachfolgend im Einzelnen geschilderten Verfahren werden allgemein die folgenden drei Prozeduren ausgeführt:
2.1 Inversion
Mit einer ersten Prozedur werden aus gegebenen Paaren von normierten Primärintensitätswerten P₁ und P₂ für die beiden Spektren die entsprechenden Belegungsdichten b₁ und b₂ der beiden zu separierenden Materialen, zum Beispiel Weichteilgewebe und Knochen, berechnet. Die Normierung erfolgt durch Division der gemessenen Intensitätswerte I₁ und I₂ durch die entsprechenden ungeschwächten Intensitätswerte I_1,0 und I_2,0. Diese Prozedur bezeichnen wir mit: M–1(P) = b (#1) wobei die Vektoren P = (P₁, P₂), b = (b₁, b₂) bedeuten. Diese Prozedur ist eine Punktoperation, die für jedes Pixel (x, y) des Röntgendetektors 9 einzeln durchgeführt werden kann. Der physikalische Hintergrund wird nachfolgend im Zusammenhang mit der Erläuterung zu den Ausführungsbeispielen noch näher erläutert werden.
2.2 Schätzung der Streustrahlung
Des Weiteren wird eine zweite Prozedur ausgeführt, mit der aus vorgegebenen räumlichen Primärintensitätsverteilungen P₁(x, y) und P₂(x, y) zu beiden Spektren und/oder vorgegebenen räumlichen Belegungsdichteverteilungen b₁(x, y) und b₂(x, y) der beiden zu separierenden Materialen Streuintensitätsverteilungen S₁(x, y) und S₂(x, y) für die beiden Spektren berechnet werden können. Diese Prozedur bezeichnen wir mit Sestim(P; b) = S (#2)wobei die Vektoren P = (P₁, P₂), b = (b₁, b₂), S = (S₁, S₂) bedeuten und jeweils Funktionen von (x, y) sind. Das Semikolon zwischen den Variablen auf der linken Seite soll andeuten, dass gegebenenfalls die explizite Abhängigkeit entweder von der einen oder von der anderen Variablen entfallen kann. Der Index „estim" soll darauf hinweisen, dass es sich um Rechenmodelle handelt, mit denen die Streustrahlung geschätzt werden soll.
Für die zweite Prozedur kommt eine Erweiterung eines Faltungsmodells von der in ZELLERHOFF, M.; SCHOLZ, B.; RÜHRNSCHOPF, E.-P.; BRUNNER, T.: Low contrast 3D reconstruction from C-arm data. In: Proceedings of SPIE. Medical Imaging, 2005, Vol. 5745, Seiten 646 bis 655 skizzierten Art in Frage. Ein derartiges Faltungsmodell lässt sich als relativ schneller Algorithmus implementieren. Daneben kommen auch mit heutigen Standardrechnern noch zeitaufwendige Monte-Carlo-Simulationen in Frage. Beide Alternativen werden im Zusammenhang mit konkreten Ausführungsbeispielen in den Abschnitten 4.1 und 4.2 näher erläutert.
2.3 Streustrahlungskorrektur
Als dritte Prozedur wird ein Streustrahlkorrekturalgorithmus ausgeführt, mit den vorgegebenen Intensitätsverteilungen P ~1(x, y), P ~2(x, y),welche Überlagerungen aus Primär- und Streustrahlung sind, mit gegebenen Streuintensitätsverteilungen S₁(x, y) und S₂(x, y) so korrigiert werden, dass man daraus – zumindest näherungsweise – Primärintensitätsverteilungen P₁(x, y) und P₂(x, y) erhält. Eine solche Korrekturprozedur wird als
bezeichnet. Die wichtigsten Korrekturalgorithmen werden im Abschnitt 4.3 erläutert.
3. Iterationsverfahren
Die Gleichungen (#1), (#2) und (#3) können iterativ gelöst werden. Iterationsalgorithmen bestehen allgemein aus drei Teilen: einem Iterationsbeginn mit vorbestimmten Anfangswerte, einer Iterationsvorschrift, durch die die Anfangswerte verändert werden, und einer Abbruchbedingung.
Mit den beschriebenen drei Prozeduren kann nun im Rahmen der dualen Röntgenabsorptiometrie das folgende auf zwei Komponenten angewandete Iterationsverfahren durchgeführt werden:
3.1 Iterationsbeginn
Zunächst sind nur die gemessenen Intensitätsverteilungen I₁(x, y) und I₂(x, y) für die beiden Spektren und die zugehörigen Intensitäten bei ungeschwächter Strahlung I_1,0 und I_2,0 als gegeben anzunehmen. Da wir über die Streustrahlung noch nichts wissen, setzen wir als Startwerte die unkorrigierten normierten Intensitätsverteilungen an:
Da P ⁽⁰⁾ noch nicht streustrahlungskorrigiert ist, kann es vorkommen, dass eine der beiden Materialbelegungen negativ herauskommt. Diese ist dann = 0 zu setzen. Bei der Weichteil/Knochen-Separation wird es in der Regel die Knochenkomponente sein, der zunächst negative Werte zugewiesen werden.
3.2 Iterationsvorschrift
Die Iterationsvorschrift bestimmt den Übergang von der Iteration mit dem Index n(n ≥ 0) zur nächsten Iteration mit dem Index n + 1. Im Folgenden lassen wir die Pixelkoordinaten (x, y) weg. Dabei ist zu beachten, dass wir es mit zweidimensionalen Funktionen zu tun haben. S (n) = Sestim(P (n); B (n)) (#5a) P (n+1) = Skorr(P (n), S (n)) (#5b) b (n+1) = M–1(P (n+1)) (#5c)
Wieder kann es vorkommen, dass eine der beiden Materialbelegungen negativ herauskommt. Diese ist dann = 0 zu setzen.
3.3 Iterationsabbruch
Die Iterationsfolge kann abgebrochen werden, wenn die Änderung zwischen aufeinander folgenden Iterationen kleiner als eine vorzugebende kleine Schwelle ist. In den meisten Fällen, besonders bei relativ kleinem bis moderatem Streuanteil (S <= P), genügen in der Regel ein bis drei Iterationen. Man kann auch die maximale Anzahl durchzuführender Iterationen fest oder abhängig vom Verhältnis S/P vorgeben.
4. Einzelheiten zu den auszuführenden Operationen
Im Abschnitt 2 wurden die wesentlichen mathematischen Operatoren aufgeführt, aber nur allgemein formal, noch nicht konkret beschrieben. Hier werden die einzelnen Operationen und Varianten davon beschrieben.
4.1 M^–1: Berechnung von Materialbelegungen aus Dualen-Energie-Projektionsdaten
Grundsätzlich gilt für die normierten Primärintensitäten bei den Spektren W₁(E) und W₂(E) im Falle von zwei radiologisch unterschiedlichen Materialien mit Massenschwächungskoeffizienten α₁(E) = (μ₁/ρ₁)(E) und α₂(E) = (μ₂/ρ₂)(E) in Abhängigkeit von den Massenbelegungsdichten b₁ und b₂[g/cm²], die der Röntgenstrahl durchdringt,
In Form einer Vektorfunktion lässt sich auch schreiben:
Die effektiven Spektren W₁(E) und W₂(E) können als bekannt angenommen werden, wobei die effektiven Spektren W₁(E) und W₂(E) das Emissionsspektrum der Röntgenröhre 2, die Wirkung der spektralen Filter und die energieabhängige Detektoransprechempfindlichkeit beinhalten. In 2 sind zwei typische Spektren 12 und 13 und in 3 die energieabhängigen Massenschwächungskoeffizienten dargestellt.
Die die Primärintensitäten beschreibenden Funktionen M₁(b), M₂(b) können prinzipiell als Funktionen der Massenbelegungsdicken b = (b₁, b₂) vorausberechnet oder auch experimentell bestimmt werden. Da sie bezüglich der beiden Variablen b₁ und b₂ partiell streng monoton fallen, zumindest in den Energiebereichen und Materialien, die in der medizinischen Diagnostik von Interesse sind, lässt sich die Vektoroperation M invertieren. Die Inversion kann mittels Standardverfahren der numerischen Mathematik, zum Beispiel dem Newton-Raphson-Verfahren für Vektorfunktionen durchgeführt werden. Derartige Verfahren sind beispielsweise in PRESS, FLANNERY, TEUKOLSKY, VETTERLING: Numerical Recipes. The Art of Scientific Programming, Cambridge University Press, 1989, Seiten 268 bis 273, beschrieben.
Es sei angemerkt, dass zum Zweck der Inversion vorteilhaft ist, in den Gleichung (#6a, b) auf der rechten Seite zu logarithmieren.
Ferner sei angemerkt, dass zumindest am Anfang der Iterationen der Fall auftreten kann, dass die Inversion nur mit einer der Belegungsdicken (oder beiden) außerhalb des physikalischen Definitionsbereichs erfüllt werden kann, da die korrigierten Schätzungen der Primärintensitäten P₁ und P₂ wegen der noch ungenauen Streustrahlungskorrektur noch fehlerhaft sind. Wenn die Inverse M^–1 im gesamten vorkommenden Wertebereich der normierten Messwerte P₁ und P₂ berechnet werden soll, dann muss die Vektorfunktion gemäß Gleichung (#7) über den physikalischen Definitionsbereich von b₁ und b₂ entsprechend hinaus – insbesondere auch für fiktive negative Belegungsdicken – berechnet werden.
4.2 S_estim: Schätzmodelle für Streustrahlungsverteilungen
Für die Bestimmung der Streustrahlungsverteilung kann beispielsweise eine Verallgemeinerung des bereits erwähnten Faltungsmodells oder ein Monte-Carlo-Verfahren verwendet werden.
Bei beiden Methoden kann die in den gemessenen Gesamtintensitäten enthaltene Information dazu verwendet werden, die Streustrahlung und die Massenbelegungsdichte der beiden Materialkomponenten, meist Wasser und Knochen, zu bestimmen.
4.2.1 Streustrahlungsbestimmung mit verallgemeinertem Faltungsmodell
Bei dem Faltungsmodell nach ZELLERHOFF, M.; SCHOLZ, B.; RÜHRNSCHOPF, E.-P.; BRUNNER, T.: Low contrast 3D reconstruction from C-arm data. In: Proceedings of SPIE. Medical Imaging, 2005, Vol. 5745, Seiten 646 bis 655, wird die Streustrahlungsverteilung mit folgendem Ansatz geschätzt: S(x, y) = [pe–pC(p)]**G(x, y) (#8)dabei ist p(x, y) = –ln(P(x, y)) und P(x, y) eine Schätzung für die normierte Primärverteilung.
G(x, y) ist ein zweidimensionaler Faltungskern, der die durch die Streustrahlungspropagation bewirkte Tiefpassfilterung (Verschmierung) empirisch beschreibt.
C(p) ist eine kalibrierende Gewichtsfunktion: C(p) = C(b0(p); U, Fyz, ,α, ...) (#9a)
Diese Gewichtsfunktion beschreibt das Streu-zu-Primär-Verhältnis S/P für Standardkörper, zum Beispiel für homogene Platten oder elliptische Zylinder, in Abhängigkeit von der Schichtdicke b₀(p) von wasseräquivalentem Material. Diese Schichtdicke kann aus der logarithmierten Primärschwächung p = –ln(P) auf der Basis der bekannten Aufhärtungskorrektur eindeutig bestimmt werden.
C hängt außerdem von den physikalischen Akquisitionsparametern ab: Spannung U, Strahlungsfilter, kollimierte Feldgröße F_yz des Detektors, Streustrahlenraster, Luftspalt und ähnliche Parameter.
C kann mittels Monte-Carlo-Rechnungen im Voraus berechnet und als Tabelle abgespeichert werden. Bei festen Akquisitionsparametern ist C nur von einer Variablen, nämlich der wasser-äquivalenten Schichtdicke b₀(p), abhängig. Da zwei Aufnahmen mit zwei unterschiedlichen Energiespektren bei unterschiedlichen Röhrenspannungen gemacht werden, sind zwei Tabellen C₁ und C₂ erforderlich.
Zwei typische Tabellen sind in 6 graphisch dargestellt.
In 6 gibt eine Verhältniskurve 28 das Verhältnis C₁ der Streu- zu Primärintensität bei einer Röhrenspannung U = 60 kV und einem Vorfilter 7 mit 0,1 mm Kupfer und eine Verhältniskurve 29 veranschaulicht C₂ entsprechend eine Röhrenspannung U = 120 kV und einem Vorfilter 7 mit 0,3 mm Kupfer.
Mit der Gewichtsfunktion in Gleichung (#9a) erhält man eine vereinfachte Streustrahlungskorrektur innerhalb der Iterationsschleife in den Gleichungen (#5a–c), wobei eine aus der geschätzten Primärstrahlungen P₁ und P₂ abgeleitete wasseräquivalente Weglänge und aus dieser wiederum ein Streuanteil S/P über die Tabelle C verwendet und iterativ verbessert wird.
Die hier beschriebene Verallgemeinerung des bekannten Ansatzes besteht darin, dass die Tabellen C₁ und C₂ für die beiden Spektren, entsprechend unterschiedlichen Spannungen U₁ und U₂ und eventuell unterschiedlichen Vorfilterungen nicht nur von einer Variablen, nämlich der Wasserschichtdicke, abhängen, sondern in Abhängigkeit von zwei Variablen b₁ und b₂, der Wasser- und Knochendicke, erstellt werden: C1 = C(b1, b2; U1, Fyz, , α, ...) (#9b) C2 = C(b1, b2; U2, Fyz, , α, ...) (#9c)
Die aktuellen Werte der beiden Materialdicken b₁ und b₂ ergeben sich bei jedem Iterationsschritt neu. Die neuen Werte für die Materialdicken b₁ und b₂ sind dann als neue Eingangswerte in die jetzt von den beiden Parametern b₁ und b₂ abhängigen Tabellen C₁ und C₂ zu verwenden, um eine iterativ verbesserte Schätzung der Streustrahlung zu erhalten.
Eine explizite Abhängigkeit von p oder P tritt nicht mehr in Erscheinung. Diese wird durch die Abhängigkeit von b₁ und b₂ ersetzt.
4.2.2 Streustrahlungsbestimmung mit Monte-Carlo-Methode
Grundsätzlich ist auch denkbar, zu den aktuell bei jedem Iterationsschritt neu geschätzten beiden Materialdickenverteilungen b₁(x, y) und b₂(x, y) die Streustrahlungsverteilung direkt per Monte-Carlo-Methode zu berechnen. Allerdings ist die Mon te-Carlo-Simulation sehr rechenaufwendig, so dass die Anwendung des Monte-Carlo-Verfahrens eine hohe Rechenleistung voraussetzt.
4.3 S_korr: Streustrahlungskorrekturalgorithmen
Die Streustrahlungskorrektur kann sowohl subtraktiv als auch multiplikativ durchgeführt werden. Dabei sind jeweils zwei Spektren zu korrigieren.
Die subtraktive Streustrahlungskorrektur ist besonders einfach und besteht in der Subtraktion der geschätzten Streustrahlungsintensitäten von den jeweiligen unkorrigierten normierten Intensitätsverteilungen: P( n +1)1 (x, y) = I'1 '(x, y) – S(n)1 (x, y) P(n+1)2 (x, y) = I'2 '(x, y) – S(n)2 (x, y) (#10a)(x, y) bezeichnen Pixelkoordinaten auf dem Detektor und für die Intensitäten I₁(x, y) und I₂(x, y) gilt I'1 (x, y) = I1(x, y)/I1,0 und I'2 (x, y) = I2(x, y)/I2,0
Da ein konstanter Wert nur eine sehr grobe Näherung ist, kann es vorkommen, dass in Gleichung (#5) bei der Subtraktion negative Werte auftreten. Solche physikalisch unsinnige Werte sollten nach Möglichkeit verhindert werden. Eine Gegenmaßnahme ist, anstelle der subtraktiven Korrektur in (#10a) eine multiplikative Streustrahlungskorrektur zu verwenden:
Für den Fall S(n)1 (x, y) << I1(x, y) und S(n)2 (x, y) << I2(x, y) geht die multiplikative Korrektur nach Gleichungen (#10b) in die subtraktive Korrektur gemäß Gleichungen (#10a) über.
5. Ortsabhängigkeit
Die räumliche Verteilung der Streustrahlung ist, abgesehen von Rauschen, im Allgemeinen glatt und daher niederfrequent. Dies bedeutet, dass es ausreicht, die Streustrahlung an sehr wenigen Stellen auf der Detektorfläche punktuell oder in interessieren Bereichen zu bestimmen. Die einfachste Näherung ist daher ein passender konstanter Mittelwert für die Sekundärstrahlungsintensität.
Zur pauschalen Schätzung der mittleren Streustrahlungshintergrundintensität genügt es, gemäß 4 einen geeigneten interessierenden Korrekturbereich 23 oder 24 auszuwählen und jeweils einen Mittelwert über diesen interessierenden Bereich für das nieder- und das höherenergetische Projektionsbild zu bilden. Für dieses Wertepaar I 1, I 2 wird dann das entsprechende Wertepaar S 1, S 2 der Streustrahlungsintensitäten durch Iteration bestimmt. Die hier verwendeten Überstriche sollen zum Ausdruck bringen, dass es sich um geschätzte, gemittelte oder sogar konstante Werte handelt.
Die Streustrahlungskorrektur außerhalb des Korrekturbereiches 23 oder 24 besteht dann in der Subtraktion der geschätzten Streustrahlungsintensitäten von den jeweiligen unkorrigierten normierten Intensitätsverteilungen: P ~1(x, y) = I1(x, y) – S 1 P ~2(x, y) = I2(x, y) – S 2 (#11)
Auf der linken Seite stehen die korrigierten Primärverteilungen, während S 1 und S 2 die nach dem Ende der Iteration erhalten, Streustrahlungsintensitäten sind. (x, y) bezeichnen Pixelkoordinaten auf dem Detektor. Die Tilde soll darauf hinweisen, dass es sich um korrigierte Daten handelt, also Schätzungen auf Grund einer Korrektur.
Da ein konstanter Wert nur eine sehr grobe Näherung ist, kann es vorkommen, dass in Gleichung (#11) bei der Subtraktion negative Werte auftreten. Solche physikalisch unsinnige Werte müssen auf jeden Fall verhindert werden. Eine Maßnahme ist die Wahl eines geeigneten Korrekturbereiches 23 oder 24 im Bereich starker Schwächung zur pauschalen Bestimmung der Streustrahlung. Ein Bereich mit starker Schwächung ist ein Bereich mit kleinen I_1,2-Werten.
Eine andere Maßnahme ist, anstelle der subtraktiven Korrektur in Gleichung (#11) eine multiplikative Streustrahlungskorrektur vorzunehmen:
Für den Fall S 1 << J1(x, y) und S 2 << J2(x, y) gehen Gleichungen (#12) in (#11) über.
Die Ortsabhängigkeit des Streustrahlungshintergrunds kann ferner dadurch behandelt werden, dass man das weiter oben im Zusammenhang für einen einzelnen Bildpunkt beschriebene Korrekturverfahren auf einem gleichmäßigen groben Raster von interessierenden Bereichen oder Abtastpunkten auf dem Detektor anwendet und das Ergebnis mittels Interpolation vom groben Raster auf das ursprüngliche feine Pixelraster expandiert. Die Korrekturen gemäß Gleichungen (#10) und (#11) sind dann sinngemäß zu erweitern: S 1, S 2 sind dann nicht mehr Konstan ten, sondern von den Pixelkoordinaten (x, y) abhängige Funktionen, selbst wenn diese Abhängigkeit in der Regel nur schwach ist.
Die Ortsabhängigkeit kann ferner dadurch erfasst werden, dass das für einen einzelnen Bildpunkt oder einen interessierenden Bereich beschriebene Verfahren auf einem gleichmäßigen groben Raster von interessierenden Bereichen oder von Abtastpunkten auf dem Detektor angewendet wird. Die Interpolation auf das ursprüngliche feine Raster kann dann durchgeführt werden, indem die auf dem groben Raster ermittelten Werte für die Streustrahlungsanteile mit einer breiten Verteilungsfunktion gefaltet werden.
Den Umstand, dass die Streustrahlung räumlich sehr niederfrequent ist, berücksichtigen auch herkömmliche Korrekturverfahren durch glättende Faltungsoperationen. Entsprechende Verfahren sind in ZELLERHOFF, M.; SCHOLZ, B.; RÜHRNSCHOPF, E.-P.; BRUNNER, T.: Low contrast 3D reconstruction from C-arm data. In: Proceedings of SPIE. Medical Imaging, 2005, Vol. 5745, Seiten 646 bis 655 beschrieben. Die Kombination des hier beschriebenen Ansatzes mit solchen aufwendigeren Faltungsmodellen ist grundsätzlich möglich. Beispielsweise können die Streustrahlungsanteile durch eines der hier beschriebenen Verfahren und parallel dazu mit einem herkömmlichen Verfahren bestimmt und die Ergebnisse anschließend gemittelt werden.
6. Verfahrensablauf
Die hier beschriebenen Verfahren zur Streustrahlungskorrektur sind jeweils prärekonstruktive Verfahren, die den Streustrahlungsanteil anhand der Schwächungsbilder bestimmen, ohne auf zumindest näherungsweise bestimmte dreidimensionale Volumenbilder des zu untersuchenden Objekts zurückzugreifen. Bei den Verfahren wird vielmehr die zusätzliche Information, die durch die im unterschiedlichen Energiebereichen aufgenommenen Schwächungsbilder zur Verfügung steht, dazu verwendet, den Streustrahlungsanteil abzuschätzen.
6 zeigt ein Ablaufdiagramm des Verfahrens.
Dieses Verfahren beginnt mit einer Wertzuweisung 30, der eine Inversion 31 gemäß Gleichung (#1) folgt. Dabei ergeben sich für eine mehrdimensionale Massenbelegung charakteristische Größen, anhand denen eine Schätzung 32 gemäß Gleichung (#2) durchgeführt werden kann. Anschließend erfolgt eine Korrektur 33 der Startwerte für die Primärstrahlungsanteile gemäß Gleichung (#3). Mit Hilfe einer Abfrage 34 kann schließlich entschieden werden, ob die Inversion 31 mit Hilfe der korrigierten Primärstrahlungsanteile wiederholt wird.
7. Vereinfachtes Verfahren
Neben den bisher beschriebenen Verfahren ist es auch möglich, für die einzelnen Schwächungsbilder jeweils eine separate Streustrahlungskorrektur durchzuführen und dann anhand der erhaltenen Primärstrahlungsintensitäten die Inversion gemäß Gleichung (#1) durchzuführen. Beispielsweise kann für das niederenergetische Schwächungsbild und das hochenergetische Schwächungsbild jeweils eine räumlich niederfrequente Streustrahlungsverteilung bestimmt und von den Schwächungswerten des niederenergetischen oder hochenergetischen Schwächungsbilds abgezogen werden.
Die Bestimmung der Sekundärstrahlungsanteile kann messtechnisch mithilfe der Beamstop-Methode oder mithilfe eines Rechenverfahrens erfolgen. Die Beamstop-Methode hat den Nachteil, dass sie zwei zusätzliche Messungen erfordert. Diese Messungen können allerdings wegen der abgesehen vom Rauschen niederfrequenten Natur der Streustrahlungsverteilung mit sehr geringer Dosis und nur geringer zusätzlicher Strahlenbelastung aufgenommen werden. Wenn die Sekundärstrahlungsanteile mithilfe eines der bekannten Rechenverfahren durchgeführt wird, brauchen keine zusätzlichen Messungen am Patienten durchgeführt werden. Die bekannten Verfahren zur Streustrahlungskorrektur können auch in Echtzeit durchgeführt werden, benötigen aber unter Umständen den Zugriff auf eine Datenbasis, die im Voraus im Rahmen der Systementwicklung durch Messungen oder Monte-Carlo-Simulation erstellt werden muss.
Die anschließende Invertierung der Gleichung (#1) kann im Allgemeinen nicht analytisch durchgeführt werden. Eine numerische Invertierung mit dem erwähnten Newton-Raphson-Verfahren ist dagegen in der Regel möglich. Die Invertierung kann auch vorab durchgeführt werden und in Tabellenform abgespeichert werden, so dass im konkreten Anwendungsfall die Lösung der Gleichung (#1) durch ein Nachschlagen in einer Lösungstabelle ersetzt werden kann.
8. Demonstrationsbeispiele
Anhand der 8 bis 13 wird ein Beispiel aus der medizinischen Bildgebung beschrieben, das die Rekonstruktion von Massenbelegungsbildern aus Schwächungsbildern veranschaulicht.
Ein typischer Anwendungsfall für die duale Röntgenabsorptiometrie ist die Trennung zwischen Weichteilgewebe und Knochengewebe bei einer Thorax- oder Lungenaufnahme, um Lungenknoten detektieren zu können.
Das Verfahren wird dabei anhand eines Phantoms 35 getestet. 8 zeigt einen Querschnitt durch das Phantom 35, während 9 eine Aufsicht des Phantoms 35 zeigt.
Das Phantom 35 umfasst quaderförmiges Weichteilgewebe 36, auf dem sich mittig angeordnet quaderförmiges Tumorgewebe 37 befindet, das ebenfalls aus Weichteilgewebe besteht. Über die eine Hälfte des Weichteilgewebes 36 erstreckt sich nun Knochengewebe 38, das auch das Tumorgewebe 37 zur Hälfte bedeckt.
In 10 ist nunmehr eine Röntgenaufnahme 39 im Niederenergiebereich gezeigt. Für die Röntgenaufnahme 39 im Niederenergiebereich wurde eine Aufnahme mit U = 60 kV und 0,1 mm Kupfervorfilterung durch ein Monte-Carlo-Verfahren simuliert. 11 zeigt dagegen eine Röntgenaufnahme 40, die im Hochenergiebereich aufgenommen worden ist. Insbesondere wurde eine Röhrenspannung von U = 120 kV und ein Kupfervorfilter mit 0,3 mm Dicke simuliert. In den Röntgenaufnahmen 39 und 40 wurden der Darstellbarkeit halber die unterschiedlichen Grauwerte durch Schraffuren unterschiedlicher Dichte wiedergegeben. Dabei gilt, dass die Dichte der Schraffur umso größer ist, je dunkler der entsprechende Bildbereich ist.
Während in den 10 und 11 die erhaltenen Schwächungsbilder dargestellt sind, zeigen die 12 und 13 die Massenbelegungsbilder, die mithilfe des im vorhergehenden Abschnitt 8 beschriebenen vereinfachten Verfahren erstellt worden sind. Insbesondere zeigt 12 ein Weichteilbild 41 und 13 ein Knochenbild 42. Auch beim Weichteilbild 41 und beim Knochenbild 42 wurde der Grauwert durch die Dichte der Schraffur wiedergegeben. Bei dem Weichteilbild 41 und dem Knochenbild 42 gilt dabei, dass sich umso mehr Material im Strahlengang befindet, je heller das Weichteilbild 41 oder das Knochenbild 42 ist. Anhand des Weichteilbildes 41 ist deutlich erkennbar, dass im Bereich des Tumorgewebes 37 die Helligkeit im Wesentlichen konstant ist, unabhängig davon, ob die Bildhälfte vom Knochengewebe 38 abgedeckt worden ist oder nicht. Allerdings bleibt der Bereich des Knochengewebes 38 auch im Weichteilbild 41 anhand des stärkeren Bildrauschens erkennbar. Dies ist insbesondere auch auf das schlechte Signal-zu-Rausch-Verhältnis im Bereich des Knochengewebes 38 zurückzuführen, da die Röntgenstrahlung dort stärker absorbiert wird. Eine Abhilfe hierfür ist eine Rausch-adaptive Bildfilterung, die in Bereichen mit höherem Rauschen lokal stärker glättet.
9. Vorteile
Durch den Iterationszyklus einschließlich Streustrahlungskorrektur lässt sich die Genauigkeit der Material-Separation deutlich steigern.
Ferner wird durch die iterative Einbeziehung eine Verbesserung der Streustrahlungskorrektur dadurch ermöglicht, dass das Streustrahlungsschätzmodell mit der zusätzlichen Information über die Belegungsdicken von zwei Materialkomponenten – die ohne Iteration nicht zur Verfügung stünden – verfeinert werden kann. Insbesondere kann die Streustrahlung auch bei der Anwesenheit von Weichteilgewebe und Knochen durchgeführt werden.
Wegen des zugrunde liegenden physikalisch-mathematischen Modells, das mit formalen Variablen oder Parametern beschrieben wird, lässt dich der Berechnungsvorgang systematisch an veränderte Rahmenbedingungen auf einfache Weise durch Änderung der Parameter oder Parameterfunktionen anpassen. Zum Beispiel kann eine Änderung der Spektren durch Wahl anderer Spannungen, durch Modifikation der Strahlungsfilter oder durch eine Änderung der Aufnahmegeometrie berücksichtigt werden.
Grundsätzlich eröffnet das hier beschriebene Verfahren die Möglichkeit, die Parameter der Mehr-Spektren-Bildgebung bezüglich geeigneter Bewertungsmaße bei Vorgabe der zu separierenden Materialien zu optimieren. Insbesondere kann die verwendete Spannung, die Art der Filter, und die Detektordosis für eine Verbesserung des Signal-zu-Rausch-Verhältnis im Verhältnis zur Patientendosis eingesetzt werden.
ZITATE ENTHALTEN IN DER BESCHREIBUNG
Diese Liste der vom Anmelder aufgeführten Dokumente wurde automatisiert erzeugt und ist ausschließlich zur besseren Information des Lesers aufgenommen. Die Liste ist nicht Bestandteil der deutschen Patent- bzw. Gebrauchsmusteranmeldung. Das DPMA übernimmt keinerlei Haftung für etwaige Fehler oder Auslassungen.
Zitierte Nicht-Patentliteratur

- WARP, R. J.; DOBBINS, J. T.: Quantitative evaluation of noise reduction strategies in dual-energy imaging In: Med. Phys. 30 (2), Seiten 190 bis 198, Februar 2003 [0002]
- HINSHAW, D. A.; DOBBINS III, J. T.: Recent Progress in noise reduction and scatter correction in dual-energy imaging. In: Proc. SPIE, 1995, Vol. 2432, Seiten 134 bis 142 [0009]
- FLOYD, C. B.; BAKER, J. A.; LO, J. Y.; RAVIN, C. E: Posterior Beam-Stop Method for Scatter Fraction Measurement in Digital Radiography. In: Investigative Radiology Feb. 1992, Vol. 27, Seiten 119 bis 123 [0010]
- ZELLERHOFF, M.; SCHOLZ, B.; RÜHRNSCHOPF, E.-P.; BRUNNER, T.: Low contrast 3D-reconstruction from C-arm data. In: Proceedings of SPIE. Medical Imaging, 2005, Vol. 5745, Seiten 646 bis 655 [0011]
- ZELLERHOFF, M.; SCHOLZ, B.; RÜHRNSCHOPF, E.-P.; BRUNNER, T.: Low contrast 3D reconstruction from C-arm data. In: Proceedings of SPIE. Medical Imaging, 2005, Vol. 5745, Seiten 646 bis 655 [0065]
- PRESS, FLANNERY, TEUKOLSKY, VETTERLING: Numerical Recipes. The Art of Scientific Programming, Cambridge University Press, 1989, Seiten 268 bis 273 [0078]
- ZELLERHOFF, M.; SCHOLZ, B.; RÜHRNSCHOPF, E.-P.; BRUNNER, T.: Low contrast 3D reconstruction from C-arm data. In: Proceedings of SPIE. Medical Imaging, 2005, Vol. 5745, Seiten 646 bis 655 [0083]
- ZELLERHOFF, M.; SCHOLZ, B.; RÜHRNSCHOPF, E.-P.; BRUNNER, T.: Low contrast 3D reconstruction from C-arm data. In: Proceedings of SPIE. Medical Imaging, 2005, Vol. 5745, Seiten 646 bis 655 [0109]

Claims

Verfahren für die Erstellung von Massenbelegungsbildern anhand von in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbildern mit den Verfahrensschritten: – Erzeugen von Strahlung (6) mit einer Strahlungsquelle (2) und Durchleuchten eines zu untersuchenden Objekts (8) mit Hilfe der Strahlung (6); – Beaufschlagen einer Detektorvorrichtung (9) mit der Strahlung (6) und Erfassen von Schwächungsbildern in unterschiedlichen Energiebereichen durch die Detektorvorrichtung (9); und – Bestimmen von Massenbelegungsbildern durch eine der Detektorvorrichtung (9) nachgeschaltete Auswerteeinheit (10), dadurch gekennzeichnet, dass von der Auswerteeinheit (10) durch Inversion einer mehrdimensionalen Schwächungsfunktion, durch die eine mehrdimensionale Massenbelegung im Objekt (8) mit den Schwächungswerten der in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbildern verknüpft wird, Massenbelegungswerte für die Massebelegung im Objekt (8) bestimmt werden.
Verfahren nach Anspruch 1, dadurch gekennzeichnet, dass durch die Auswerteinheit (10) ein durch Streuung hervorgerufener Sekundärstrahlungsanteil bestimmt und die Schwächungsbilder hinsichtlich des Sekundärstrahlungsanteils (26, 27) auf einen durch Schwächung erzeugten Primärstrahlungsanteil (25) korrigiert werden.
Verfahren nach Anspruch 2, dadurch gekennzeichnet, dass zur Bestimmung des Sekundärstrahlungsanteils (26, 27) der in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbilder diejenigen Sekundärstrahlungsanteile (26, 27) in einem Korrekturbildbereich (23, 24) gesucht werden, die mit Primärstrahlungsanteilen (25) verknüpft sind, aus denen sich beim Auswerten inverser Schwächungsfunktionen jeweils die gleiche Massenbelegung für den Korrekturbildbereich (23, 24) ergibt.
Verfahren nach Anspruch 3, dadurch gekennzeichnet, dass der Korrekturbildbereich (24) einen hinsichtlich der Schwächungskoeffizienten inhomogenen Objektbereich abbildet und die Schwächungsfunktion von mehrdimensionalen Massenbelegungen abhängen.
Verfahren nach Anspruch 4, dadurch gekennzeichnet, dass die Sekundärstrahlungsanteile (26, 27) in den in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbildern anhand der ermittelten mehrdimensionalen Massenbelegung näherungsweise bestimmt werden.
Verfahren nach Anspruch 5, dadurch gekennzeichnet, dass die Sekundärstrahlungsanteile (26, 27) jeweils durch Faltung der jeweiligen Primärstrahlungsanteile (25) mit einer von der mehrdimensionalen Massenbelegung abhängigen Verteilungsfunktion ermittelt werden.
Verfahren nach Anspruch 6, dadurch gekennzeichnet, dass die Verteilungsfunktion vorab mit Hilfe einer Monte-Carlo-Simulation ermittelt werden.
Verfahren nach Anspruch 5, dadurch gekennzeichnet, dass die Sekundärstrahlungsanteile (26, 27) in den jeweiligen Schwächungsbildern durch eine die mehrdimensionale Massenbelegung berücksichtigende Monte-Carlo-Simulation der Streuung am zu untersuchenden Objekt (8) ermittelt werden.
Verfahren nach einem der Ansprüche 4 bis 8, dadurch gekennzeichnet, dass nach der Bestimmung der mehrdimensionalen Massenbelegung die zugehörigen Sekundärstrahlungsanteile (26, 27) gesucht, die Primärstrahlungsanteile (25) korrigiert und die mehrdimensionale Massenbelegung ausgehend von den korrigierten Primärstrahlungsanteilen (25) erneut bestimmt wird.
Verfahren nach einem der Ansprüche 4 bis 9, dadurch gekennzeichnet, dass der Bildbereich (24) mit inhomogenen Schwächungskoeffizienten einem Körperbereich eines Patienten (8) zugeordnet ist, der Knochen oder nicht aus auf Kohlenwasserstoff basierenden Kunststoffen hergestellte Implantate umfasst.
Verfahren nach einem der Ansprüche 3 bis 10, dadurch gekennzeichnet, dass der Korrekturbildbereich (23, 24) einen einzelnen Bildpunkt der in unterschiedlichen Energiebereichen aufgenommenen Schwächungsbilder umfasst.
Verfahren nach einem der Ansprüche 3 bis 11, dadurch gekennzeichnet, dass im Korrekturbildbereich (23, 24) Bildwerte der in verschiedenen Energiebereichen aufgenommenen Schwächungsbilder gemittelt und der Sekundärstrahlungsanteil anhand der gemittelten Bildwerte bestimmt wird.
Verfahren nach einem der Ansprüche 2 bis 12, dadurch gekennzeichnet, dass aus dem Sekundärstrahlungsanteil im Korrekturbildbereich (23, 24) eine Sekundärstrahlungsintensität ermittelt wird, die von Bildwerten außerhalb des Korrekturbildbereichs (23, 24) subtrahiert wird.
Verfahren nach einem der Ansprüche 2 bis 13, dadurch gekennzeichnet, dass anhand des Sekundärstrahlungsanteils im Korrekturbildbereich (23, 24) eine Sekundärstrahlungsintensität berechnet wird, die zur Bildung eines Korrekturfaktors für die Korrektur von Bildwerten außerhalb des Korrekturbildbereichs (23, 24) herangezogen wird.
Verfahren nach einem der Ansprüche 3 bis 14, dadurch gekennzeichnet, dass Sekundärstrahlungsanteile in den Korrekturbildbereichen (23, 24) eines sich über die Schwächungsbilder erstreckenden Rasters bestimmt werden.
Verfahren nach Anspruch 15, dadurch gekennzeichnet, dass Bildwerte außerhalb der Korrekturbildbereiche (23, 24) in Abhängigkeit von dem im benachbarten Korrekturbildbereichen (23, 24) bestimmten Sekundärstrahlungsanteilen korrigiert werden.
Vorrichtung für die Aufnahme von Schwächungsbildern, dadurch gekennzeichnet, dass die Vorrichtung zur Durchführung eines Verfahrens nach einem der Ansprüche 1 bis 16 eingerichtet ist.