DE102009043213A1

DE102009043213A1 - Effiziente Korrektur von Polychromieeffekten bei der Bildrekonstruktion

Info

Publication number: DE102009043213A1
Application number: DE102009043213A
Authority: DE
Inventors: Holger Dr. Kunze; Holger Scherl
Original assignee: Siemens AG
Current assignee: Siemens Healthcare GmbH
Priority date: 2009-09-28
Filing date: 2009-09-28
Publication date: 2011-03-31
Also published as: US20110075899A1; US8259896B2

Abstract

Die Erfindung betrifft die Ermittlung von bzgl. Polychromieartefakten korrigierten Absorptionskoeffizienten für ein aus einer Mehrzahl bzgl. der Absorptionseigenschaften sich unterscheidenden Materialtypen zusammengesetztes Objekt. Dabei wird eine Vielzahl von Röntgenstrahl-Projektionen aus unterschiedlichen Positionen für das Objekt aufgenommen und durch ein Rekonstruktionsverfahren mittels der aufgenommenen Projektionen ein erster Satz von Absorptionskoeffizienten ermittelt. Mittels Vorwärtsprojektion werden Projektionen berechnet und die aufgenommen Projektionen mit Hilfe der berechneten Projektionen im Sinne von mit monochromen Röntgenstrahlen aufgenommenen Projektionen korrigiert. Durch Rekonstruktion der korrigierten Projektionen wird schließlich ein zweiter, bzgl. Polychromieartefakten korrigierter Satz von Absorptionskoeffizienten ermittelt. Dabei wird eine formelmäßige Beschreibung einer Polychromie berücksichtigende Vorschrift verwendet, welche zu bestimmende Parameter umfasst und die Vorschrift für die Berechnung von Projektionen mittels Vorwärtsprojektion verwendet. Die Bestimmung der Parameter erfolgt im Zuge der Berechnung von Projektionen mittels Vorwärtsprojektion. Das Verfahren vereinigt Schritte herkömmlicher Verfahren und ist damit effizienter.

Description

Die Erfindung betrifft ein Verfahren, eine Vorrichtung und ein Computerprogramm zur Ermittlung von bzgl. Polychromieartefakten korrigierten Absorptionskoeffizienten für ein aus einer Mehrzahl von bzgl. der Absorptionseigenschaften sich unterscheidenden Materialtypen zusammengesetztes Objekt.
Die Computer-Tomographie (CT) stellt ein Diagnose- und Messverfahren für Medizin und Prüftechnik zur Verfügung, mit dessen Hilfe innere Strukturen eines Patienten und Prüfobjektes untersucht werden können, ohne dabei operative Eingriffe an dem Patienten durchführen oder das Prüfobjekt beschädigen zu müssen. Dabei wird von dem zu untersuchendem Objekt eine Anzahl Projektionen aus verschiedenen Winkeln aufgenommen. Aus diesen Projektionen läst sich eine 3D Beschreibung des Objekts berechnen.
1 zeigt eine typische CT-Anordnung mit einer Röntgenquelle 1 in einer ersten Position, die ein Röntgenstrahlenbündel 2 für eine erste Projektion aussendet, das in einem Detektor 3 an einer ersten Position nach gewiesen wird, nachdem es das zu untersuchende Objekt 4 durchdrungen hat. Die Daten des Detektors gelangen in einen Auswertungsrechner 5, der die Rekonstruktion vornimmt, und werden anschließend auf einer Anzeigeeinheit 6 dargestellt. Die Röntgenquelle bewegt sich idealerweise auf einer Kreisbahn, wobei zahlreiche Projektionen aufgenommen werden. zeigt beispielhaft die Röntgenquelle in einer anderen Position 11, wobei ein Röntgenstrahlenbündel 12 für eine andere Projektion ausgesandt wird, das dann im Detektor an einer anderen Position 13 nachgewiesen wird.
Bei der anschließenden Verarbeitung der Messdaten geht man von der Annahme aus, dass für die gemessenen Daten das Lambert-Beer'sche Gesetz
gilt. Dabei stellt I die mittels des Detektors 3 gemessene Intensität, I₀ die unabgeschwächte Intensität und μ(s) den Abschwächungskoeffizienten am Ort s dar. Diese Beziehung gilt jedoch nur für monochromatische Strahlung und nicht für die polychromatische Strahlung einer Röntgenröhre. Für diese gilt vielmehr die von der Photonenenergie E abhängige Beziehung
Beachtet man diese Abhängigkeit nicht, weist das rekonstruierte Bild Artefakte auf, die den rekonstruierten Schwächungswert um bis zu einige Prozent (verglichen mit einer Rekonstruktion aus monochromatisch gewonnenen Messdaten) verfälschen [Buz04].
In der Literatur sind verschiedene Verfahren bekannt, wie aus den so gemessenen Werten näherungsweise von der Photonenenergie unabhängige Werte berechnet werden können. Das für diese Anmeldung relevante Verfahren ist iterativer Natur und wird im Folgenden skizziert.
Zunächst wird angenommenen, dass die Daten monochromatisch aufgenommen wurden. Es wird eine Rekonstruktion des Objekts berechnet. Anschließend wird das Objekt in unterschiedliche Bereiche segmentiert (zum Beispiel Weichgewebe und Knochen). Hierfür sind in der Literatur mehrere Ansätze beschrieben. Im einfachsten Fall kann dies mittels eines Schwellwertverfahrens erfolgen.
Von den so segmentierten Bildern werden nun für jede Gewebeklasse künstliche Projektionen berechnet. Relevante Informationen sind jeweils der Abschwächungsanteil und die Strahllänge jedes Röntgenstrahls durch den segmentierten Gewebeanteil.
Aus den errechneten Projektionswerten und/oder Materialdicken der einzelnen Gewebeklassen lässt sich ein Korrekturfaktor bestimmen, mit dem der ursprünglich gemessene Wert korrigiert werden kann, bevor eine zweite, endgültige Rekonstruktion des Objekts berechnet wird [Cas04]. Das Korrekturverfahren ist in 2 schematisch für drei verschiedene Gewebeklassen dargestellt.
Metallähnliche Elemente wie Zahnimplantate oder Hüftprothesen erzeugen so starke Artefakte, dass eventuell das beschriebene Verfahren nicht mehr angewendet werden kann. In [Bal06] wird ein Metallartefaktkorrekturverfahren beschrieben, das ebenfalls die Berechnung verschiedener Projektionsbilder durch verschiedene Gewebeklassen erfordert. In einem ersten Schritt werden diejenigen Bereiche auf dem Projektionsbild bestimmt, dessen zugehörige Strahlen durch Metall verlaufen. In einem zweiten Schritt wird ein Projektionsbild eines Modellvolumens bestimmt. Im Modellvolumen wurden Metallbereiche durch benachbarte Gewebeklassen ersetzt. Auf Basis dieser beiden Projektionen lassen sich nun die originalen Projektionswerte korrigieren, so dass ein artefaktfreies Bild rekonstruiert werden kann. In 3 ist dieses Verfahren schematisch zusammengefasst.
Die Projektionsberechnung kann folgendermaßen erfolgen: Das Bild wird in rechteckförmige Pixel unterteilt, die jeweils über die gesamte Pixelfläche hinweg den Wert des Abschwächungskoeffizienten an dem dazugehörigen Abtastpunkt aufweisen. Das Linienintegral entlang der Strahlen kann dann als Summe der Abtastwerte gewichtet mit der durchstrahlten Länge durch die dazugehörigen Pixel berechnet werden [Sid85]. Diese Vorgehensweise ist in 4 veranschaulicht.
Aus der Literatur sind weitere Projektionsalgorithmen bekannt, die sich von dem beschriebenen Verfahren durch die Berechnung des Gewichts, mit dem ein Pixel in die Projektionsberechnung eingeht, unterscheiden [Mue98].
Bei den beschriebenen Verfahren sind bei den geschilderten Verfahren mindestens zwei Projektionsberechnungen nötig. Zudem muss zusätzlich zum Speicher für das Objekt noch weiterer Speicher für die segmentierten Objekte bereit gehalten werden.
Die Erfindung hat zur Aufgabe, die Effizienz von herkömmlichen Verfahren zur Korrektur von Polychromieartefakten zu verbessern.
Die Aufgabe wird durch die Ansprüche gelöst.
Erfindungsgemäß wird ein Verfahren zur Ermittlung von bzgl. Polychromieartefakten korrigierten Absorptionskoeffizienten für ein aus einer Mehrzahl von bzgl. der Absorptionseigenschaften sich unterscheidenden Materialtypen zusammengesetztes Objekt vorgeschlagen. Bei dem Verfahren wird eine Vielzahl von Röntgenstrahl-Projektionen aus unterschiedlichen Positionen für das Objekt aufgenommen. Durch Rekonstruktion der aufgenommenen Projektionen wird ein erster Satz von Absorptionskoeffizienten ermittelt, wobei die Rekonstruktion eine für monochrome Röntgenstrahlung geltende Absorptionsvorschrift benutzt. Für die Rekonstruktion wird ein Rekonstruktionsverfahren verwendet, z. B. gefilterte Rückprojektion, Feldkampalgorithmus, algebraische Rekonstruktionstechniken (ART/SART), Grangeat-Verfahren oder Fourierverfahren.
Mittels Vorwärtsprojektion werden dann Projektionen berechnet. Die aufgenommenen bzw. gemessenen Projektionen werden mit Hilfe der berechneten Projektionen im Sinne von mit monochromen Röntgenstrahlen aufgenommenen Projektionen korrigiert. Schließlich wird durch Rekonstruktion auf Basis der korrigierten Projektionen ein zweiter, bzgl. Polychromieartefakten korrigierter Satz von Absorptionskoeffizienten ermittelt. Dabei wird für die Korrektur eine formelmäßige Beschreibung einer Polychromie berücksichtigende Vorschrift verwendet, welche zu bestimmende Parameter umfasst. Gemäß der Erfindung erfolgt die Bestimmung der Parameter im Zuge bzw. der Berechnung von Projektionen mittels Vorwärtsprojektion bzw. parallel, zeitgleich oder simultan dazu.
Die Erfindung erlaubt eine effizientere Korrektur von Polychromieartefakten. Im Gegensatz zu herkömmlichen Verfahren, die zunächst die Formeln für die Schwächung der polychromen Strahlung bestimmen und dann die Projektionsberechnung (Vorwärtsprojektion) durchführen werden erfindungsgemäß im Rahmen der Projektionsberechnung die Formeln aufgestellt. Insbesondere brauchen keine zusätzlichen Projektionen für die Bestimmung von Formelparametern durchgeführt zu werden.
Gemäß einer bevorzugten Ausgestaltung der Erfindung wird eine Segmentierung eines Volumens des Objektes im Zuge der Berechnung von Projektionen mittels Vorwärtsprojektion vorgenommen. Diese Ausgestaltung ermöglicht eine noch aufwandsärmere Implementierung, da nur ein Projektionsschritt für die Berechnung der Projektionswerte erfolgen muss und zudem kein zusätzlicher Speicher für die Segmentierung angelegt werden muss. Der erforderliche Mehraufwand für die Segmentierung während der Projektionsberechnung bleibt durch die Verwendung einer Pipeline-Struktur vernachlässigbar. Sowohl die Rechenzeit, als auch der Ressourcenaufwand vermindert sich deutlich durch die vorgeschlagene Lösung.
Gemäß einer Weiterbildung erfolgt eine Korrektur für die Verhinderung von Metallartefakten im Zuge der Berechnung von Projektionen mittels Vorwärtsprojektion.
Die Erfindung umfasst auch eine Vorrichtung zur Durchführung eines erfindungsgemäßen Verfahrens. Die vorgeschlagenen abgeänderten Verfahren eignen sich besonders für eine Implementierung auf Spezialhardware basierend auf FPGA- und/oder ASIC-Technologie. Andere Spezielhardware (wie DSP-Lösungen oder Spezialprozessoren wie die Cell Broadband Engine oder Grafikbeschleuniger) können bei dem vorgeschlagenen Ansatz ebenfalls ausgenutzt werden.
Die Erfindung wird im Folgenden im Rahmen von Ausführungsbeispielen anhand von Figuren näher beschrieben. Es zeigen:
1: typische Anordnung für CT-Aufnahmen
2: Korrekturverfahren für Polychromieartefakte
3: Korrekturverfahren für Polychromieartefakte mit zusätzlicher Korrektur von Metallartefakten
4: Illustration eines erfindungsgemäßen Vorgehens
5: Berechnung des Linienintegrals durch ein Voxel
6: Schematische Darstellung des bisherigen Ansatzes
7: Schematische Darstellung des neuen Ansatzes
2 zeigt ein herkömmliches Verfahren zur Ermittlung von bezüglich Polychromie-Artefakten korrigierten Absorptionskoeffizienten für Objekte, welche aus einer Mehrzahl von unterschiedlichen Materialtypen zusammengesetzt sind. In der folgenden detaillierteren Beschreibung dieses Verfahrens orientiert sich die Darstellung an dem in Cas04 aufgestellten bimodalen Energiemodell. Der Ablauf eines derartigen Verfahrens ist insbesondere auch in 5.2 der genannten Quelle dargestellt. Als Ausgangspunkt für das Verfahren werden eine Vielzahl von Röntgenstrahlprojektionen 3 aus unterschiedlichen Positionen bzw. Richtungen aufgenommen (Schritt S1). Die so erhaltenen Projektionen bzw. Sinogramme 2 bilden die Grundlage für die Rekonstruktion des Objektes 1, genauer gesagt der Absorptionseigenschaften des untersuchten Objekts 1. Die Projektionen werden zunächst für eine Rückprojektion S2 verwendet. Diese geht von monochromer Röntgenstrahlung aus, d. h. sie setzt die Gültigkeit des Lambert-Beer'schen Gesetzes voraus. Als Resultat dieser Rückprojektion (ein Beispiel für ein Verfahren zur Rückprojektion siehe US 6,771,733 B2 ) erhält man einen ersten Satz von Schwächungskoeffizienten 3 für das Objekt 1, welche bezüglich der Polychromie des verwendeten Strahles nicht korrigiert sind.
Dieser Satz von Schwächungskoeffizienten bzw. diese Objektrekonstruktion 3 wird in einem nächsten Schritt S4 in verschiedene Materialklassen bzw. Materialtypen segmentiert. Diese Segmentierung entspricht einer Klassifizierung, d. h. jeder der rekonstruierten Absorptions-Koeffizientenwerte ist einer Materialklasse zugeordnet. Die einzelnen Materialklassen unterscheiden sich – wie oben gesagt – durch ihre Absorptionseigenschaften, d. h. Materialien, die sich nur in anderen physikalischen Eigenschaften als dem Absorptionsvermögen unterscheiden, werden hier als eine Materialklasse behandelt. In 2 sind exemplarisch 3 verschiedene Materialklassen gezeigt, nämlich Weichgewebe 41, normales Gewebe 42 und Knochen 43. Diese Klassifizierung kann beispielsweise mit einem sogenannten Schwellenwertverfahren durchgeführt werden. Dieses basiert auf den signifikanten Änderungen der Absorptionsrate für verschiedene Materialtypen. Danach wird, wie in [Cas04] beschrieben, für eine Richtung die Dicke des Objektes 1 bestimmt, differenziert nach Abschnitten mit verschiedenen Gewebetypen. Nach Segmentierung und Dickenbestimmung wird (für die entsprechende Richtung) eine Vorwärtsprojektion S5 durchgeführt. Diese besteht in der simulierten Bestrahlung mit Röntgenstrahlen. Dabei findet üblicherweise eine Berechnung für eine Vielzahl von parallelen Strahlen statt, wobei sich die Position der Strahlen in der Ebene senkrecht zur Strahlungsrichtung und deren Dichte an den Detektorgegebenheiten orientieren. Diese Vorwärtsprojektion S5 simuliert mit Hilfe eines mathematischen Ansatzes für die Schwächung polychromatischer Röntgenstrahlung die Strahlschwächung bzw. Projektionen beim Durchstrahlen der einzelnen Gewebetypen 51, 52, 53. Der mathematische Ansatz ist von Parametern abhängig, deren Bestimmung das Ziel dieser Simulation ist. Als Modell für die entsprechende Schwächung bzw. Absorption polychromer Röntgenstrahlung kann beispielsweise ein in [Cas04] beschriebenes bimodales Energiemodell verwendet werden, bei dem aus physikalischen Gründen von zwei im Strahl dominierenden Energien ausgegangen wird.
Ein mathematischer Ausdruck für die Abschwächung der Intensität gemäß diesem Modell ist beispielsweise in der Gleichung 3.6 in [Cas04] wiedergegeben. Möglich ist aber beispielsweise auch ein nicht physikalischer motivierter Ansatz, der geeignete Ansatzfunktionen für die Beschreibung der entsprechenden Abschwächungskurven verwendet. Wie in [Cas04] beschrieben, wird nun entlang der einzelnen Strahlen abschnittsweise für die verschiedenen Materialtypen die Intensitätsabschwächung ermittelt. Dabei ist zu berücksichtigen, dass im Inneren des Objektes liegende Materialbereiche durch weiter außen liegende Bereiche eines anderen Materialtyps abgeschirmt werden können. Diese Abschirmung ist bei der Berechnung der Strahlintensitätsabschwächung zu berücksichtigen. Dies geschieht z. B., indem für den Abschwächungsbeitrag eines innen liegenden Bereichs berücksichtigt wird, dass der Strahl bereits mit abgeschwächter Intensität diesen Bereich erreicht, wobei sich diese abgeschwächte Intensität mathematisch durch den verwendeten Ansatz ausdrücken lässt. Man erhält so für die Intensität eines Strahles nach Durchdringen des gesamten Objektes einen Ausdruck, der von dem verwendeten mathematischen Ansatz und dessen Parametern für die einzelnen Materialschichten und deren Dicke abhängt. Unbekannte Größen sind hier die Parameter für die einzelnen Materialschichten. Dieser Ausdruck für die Intensitätsabschwächung wird dann gleich der vom Detektor gemessenen Intensitätsabschwächung für eine Vielzahl von Strahlen, eine Vielzahl Gleichungssystem mit den Parametern der Gewebeklassen als Unbekannte. Dieses Gleichungssystem ist in der Regel überbestimmt und lässt sich als Minimierungsproblem durch Minimierung der Unterschiede und berechneter Intensitätsabschwächung lösen. Auf diese Weise werden die Werte für die Parameter ermittelt (Schritt S6), mit denen ein formelmäßiger Zusammenhang für die Intensitätsabschwächung der polychromen Röntgenstrahlung vorliegt.
Mittels dieser mathematischen Beschreibung für die Abschwächung polychromer Strahlung werden die gemessenen Projektionsdaten hinsichtlich der Polychromie der Röntgenstrahlung korrigiert, d. h. derart korrigiert, dass sie bei polychromer Röntgenstrahlung aufgenommenen Daten entsprechen (Schritt S7). Man erhält korrigierte Projektionen 6. Anschließend wird eine Rückprojektion (Schritt S8) durchgeführt, für welche nun das Beer'sche Gesetzt korrekt angewendet ist. Es ergibt sich ein zweiter, nun bezüglich der Monochromie korrigierter Satz von Abschwächungskoeffizienten 7, d. h. die endgültige Rekonstruktion des Objekts.
In 3 ist eine Erweiterung konventioneller Verfahren zur Korrektur gezeigt, bei dem metallische Bereiche des Objekts behandelt werden (zusätzlich bei Klassifizierung: normales Gewebe 44 und Metall 45; die Projektionen der zur Berechnung der Vorschriften für polychromatische Strahlung für die verschiednen Gewebetypen außer Metall ist mit Bezugszeichen 5 bezeichnet). Ein solches Verfahren ist genauer in [Bal06] beschrieben. Die Motivation für diese Weiterentwicklung ist vornehmlich, dass bei Gewebeuntersuchungen mitunter Metall durchstrahlt wird, beispielsweise in der Form von Metallprothesen, Chirurgieklammern oder Zahnfüllungen. Diese führen zu erheblichen Störungen. Vereinfacht dargestellt besteht die Weiterentwicklung darin, dass bei der Segmentierung als metallisch erkannte Bereiche für die Simulierung der Abschwächung entfernt und durch in der Nachbarschaft liegende Materialtypen modelliert werden (Schritt S9 in 3). Es liegt auf der Hand, dass der berechnete Abschwächungswert für einen Strahl, bei dem in der Simulation metallisches Material durch umliegendes Material ersetzt wurde, nicht mehr dem gemessenen korrespondieren kann, da das ersetzte Material eine andere Absorptionseigenschaft hat. Der errechnete Abschwächungswert für derartige Strahlen wird daher nicht für das Optimierungsproblem herangezogen, was insofern kein Problem darstellt, als dass das Optimierungsproblem in einem überbestimmten Gleichungssystem besteht. Dagegen wird im Zuge der Korrektur für metallische Bereiche durchdringende Strahlen der gemessene Wert durch den berechneten, bei dem der Einfluss des Metalls entfernt wurde, ersetzt (Schritt S10 in 3). Somit erhält man (einen durch die Substitution der berechneten Strahlen vollständigen) Satz von Projektionsdaten, die teilweise aus gemessenen und teilweise aus simulierten Daten (nämlich für die metallischen Bereiche) bestehen. Wieder wird eine Korrektur der Projektionsdaten mit Rücksicht auf die Polychromie durchgeführt und danach in einer Rückprojektion der korrigierte Datensatz gewonnen.
Der Erfindung liegt die Beobachtung zugrunde, dass bei derartigen Verfahren sich die Effizienz steigern lässt, indem Schritte parallel durchgeführt werden. Herkömmlich werden folgende Schritte sequentiell durchgeführt:

a. Segmentierung des Volumens (Schritt S4 aus 2)
b. Berechnung der Intensitätsabschwächung der polychromen Strahlen durch die verwendeten Gewebeklassen und Bestimmung von Vorschriften für die Abschwächung der polychromatischen Strahlung (Schritte S5 und S6 aus 2)
c. Berechnung der Korrektur (Schritt S7 aus 2)

Erfindungsgemäß werden die Schritte b. und c. parallel durchgeführt. In einer unten detaillierter beschrieben Weiterbildung werden alle drei Schritte a.–c. (zumindest teilweise) parallel durchgeführt.
Anhand 4 werden Prinzipien eines möglichen Vorgehens erläutert. Dabei handelt es sich um eine prinzipielle Darstellung, die zwecks Verständlichkeit darauf verzichtet, komplexere Details von realen Anwendungen wiederzugeben.
Es ist ein Objekt O gezeigt, welches mittels Röntgenstrahlung untersucht wird. Das Objekt O ist gerastert dargestellt, wobei Schnittpunkte der Rasterlinien sog. Voxel oder Pixel darstellen sollen, d. h. Positionen, für die der Schwächungskoeffizient μ ermittelt wird. Exemplarisch sind die Voxel v₁, v₂, v₃, v_n und v_n+1 benannt.
Unten und rechts in der Zeichnung sind zwei Detektorpositionen D₁ und D₂ gezeigt, für welche Aufnahmen (Durchstrahlung des Objektes O mit Röntgenstrahlung) durchgeführt werden. In der Realität werden erheblich mehr Aufnahmen durchgeführt. Zudem entspricht der Abstand zwischen den Detektorpositionen D₁ und D₂ nicht den realen Gegebenheiten.
Der Detektor (Strahlungsdetektor) besteht üblicherweise aus einzelnen Elementen (Photozelle, Photomultiplier, CCD-Sensor, ...), welche Strahlungsintensitäten registrieren. Die Größe und Anordnung der Elemente bestimmt die Granularität der Aufnahme, d. h. jedem Element und damit der Position jedes Elements ist ein gemessener Intensitätswert zuzuordnen. Exemplarisch sind für die Detektorposition D₁ mittels vertikalen Strichen die Positionen angedeutet, wo Intensitätsmessungen vorgenommen werden. Exemplarisch sind die gemessene Intensitäten I₁, I₂, I₃, I_n, I_n+1 und für die Detektorposition D₂ die Intensität I_y eingezeichnet. Diese Messergebnisse bilden die Grundlage für die Rekonstruktion des Objekts O, d. h. für die Ermittlung der Schwächungskoeffizient μ (auch als Grauwerte im Hinblick auf eine Visualisierung des Objektes O bezeichnet). Das Objekt O setzt sich aus drei verschiedenen Materialtypen A, B und C zusammen.
Die Erfindung zielt auch auf eine Rekonstruktion der Schwächungskoeffizient μ ab, welche die Polychromie der bei der Messung verwendeten Strahlung berücksichtigt. Es wird wie folgt vorgegangen:
Zunächst rekonstruiert man das Objekt ohne Berücksichtigung von Polychromieeffekten. Dies geschieht mittels einer sog. gefilterten Rückprojektion, welche von der Gültigkeit des Lambert-Beer'sche Gesetz ausgeht. Verfahren dafür sind z. B. in Bu04 beschrieben. Man erhält einen ersten Satz (nicht korrigierter) Schwächungskoeffizienten μ₁(v_i) für alle Voxel v_i. Dieser erste Satz von Schwächungskoeffizienten wird für die Segmentierung gebraucht.
Die Segmentierung entspricht einer Klassifizierung bzgl. der Materialklasse, d. h. in dem vorliegenden Fall einer Zuordnung zu einer der Materialtypen A, B oder C. Für die Segmentierung wird z. B. ein Schwellenwertverfahren bzw. ein Clustering-Verfahren verwendet. Dieses unterteilt die Skala für mögliche Schwächungswerte in Bereiche. Diejenigen Schwächungswerte, die in denselben Bereich liegen, werden dann auch derselben Materialklasse zugeordnet.
Die Segmentierung kann als separater Schritt nach der Rückprojektion erfolgen. Im Rahmen dieses Ausführungsbeispiels wird sie jedoch im Zuge der Vorwärtsprojektion durchgeführt.
Nach der Rekonstruktion des ersten Satzes von Schwächungskoeffizienten erfolgt eine Vorwärtsprojektion. Diese Vorwärtsprojektion simuliert die reale Strahlintensitätsabschwächung. D. h. für die Ermittlung der Abschwächung wird nicht mehr das Lambert-Beer'sche Gesetz verwendet. Stattdessen wird ein mathematischer Ansatz gemacht, mittels dem die Polychromie der Strahlung berücksichtig wird. Dieser Ansatz kann, aber muss nicht nach physikalischen Prinzipien aufgestellt werden. Beispiele für eine derartige Beschreibung sind in Cas04 zu finden. Eine entsprechende Formel für die Intensitätsabschwächung hat folgendes Aussehen:
I(d) = f(I_o, d, α_M,i), wobei I_o die Anfangsintensität des Strahls, d die Dicke des durchdrungenen Materials und α_M,i Parameter darstellen.
Dabei unterscheidet der Index i die verwendeten Parameter, z. B. i aus {1, 2, 3}. Es wird hier der Einfachheit angenommen, dass alle Parameter von dem Materialtyp M abhängen (d. h. M ist im vorliegenden Fall aus der Menge {A, B, C}). Die Parameter α_M,i sind natürlich auch abhängig von den Strahleigenschaften (Energiespektrum) des für die Messung verwendeten Strahls. Dies braucht aber im Folgenden als explizite Abhängigkeit nicht berücksichtigt zu werden.
Die Vorwärtsprojektion erfolgt durch ein simuliertes Durchstrahlen des Objekts O. Dieses simulierte Durchstrahlen emuliert die zu Beginn erfolgte Messung. Entsprechend orientiert sich Position und Anzahl der Strahlen an den Messergebnissen, d. h. für jeden registrierten Messwert wird die Abschwächung eines Strahls simuliert. Diese Vielzahl von Strahlen ist auf der linken Seite der Fig. durch Pfeile angedeutet. Exemplarisch sind für die Strahlen S₁, S₂ und S_v Bezugszeichen angegeben. Der Einfachheit halber ist für die vertikale Richtung nur ein Strahl S_v angegeben, obwohl auch in dieser Richtung eine Vielzahl von Strahlen simuliert wird.
Es wird nun schrittweise bzw. voxelweise vorgegangen. Das Voxel v₁ wird mittels Schwellewertverfahren der Materialklasse A zugeordnet. Ab dem zweiten Schritt bzw. dem zweiten Voxel v₂ wird zusätzlich überprüft, ob eine Materialklassenänderung stattgefunden hat. Eine Materialklassenänderung wird bei Übergang von Voxel v_n auf Voxel v_n+1 detektiert. Dies triggert die Materialdickenberechnung. Die Dicke der durchdrungenen Schicht des Materialtyps A wird aus dem Produkt der Anzahl der Schritte und der Schrittweite berechnet. Auf diese Weise werden die Voxel klassifiziert und die durchdrungenen Dicken d_A, d_B und d_C berechnet. Wenn das Ende des Objekts erreicht ist, wird eine mathematische Beschreibung der Intensitätsabschwächung des Strahls aufgestellt. Mit der Formel I(d) = f(I_o, d, α_M,i) und der Beobachtung, dass innere Schichten abgeschirmt sind, d. h. eine reduzierte Intensität sehen, ergibt sich für die berechnete Intensität I_d nach Durchdringen des Objekts: I_d = f(f(f(I_o, d_A, α_A,i), d_B, α_B,i), d_C, α_C,i) (1)
Mit der Gleichsetzung von gemessener und berechneter Intensität (I_d = I_y) für eine Vielzahl von Strahlen ergibt sich ein überbestimmtes Gleichungssystem, aus dem durch Minimierung der Abweichungen zwischen gemessenen und berechneten Werten die Parameter α_M,i bestimmt werden.
Die so erhaltene Beschreibungen für den Effekt der Polychromie werden verwendet, um die Messdaten (z. B. gemessene Intensitäten I₁, I₂, I₃, I_n, I_n+1) so zu korrigieren, dass eine erneute gefilterte Rückprojektion basierend auf dem Lambert-Beer'sche Gesetz korrigierte Abschwächungskoeffizienten ergibt.
Es ist sinnvoll, das Gleichungssystem wegen der Anzahl der Parameter schrittweise durch entsprechende Auswahl von Gleichungen zu lösen.
Zunächst werden für die einzelnen Strahlen S die Dicken der durchdrungen Schichten bestimmt. Dies wird im Folgenden anhand des Strahls S_y beschrieben. Es wird schrittweise (in der Figur ist die Schrittweite entsprechend dem Voxelabstand gewählt) entlang dem Strahl S_y entlang gegangen. Dabei wird für jedes Voxel (v₁, v₂, v₃, ... v_n, v_n+1 ...) eine Klassifizierung des Voxels durchgeführt. Gleichzeitig wird überprüft, ob sich die Gewebeklasse ändert. Beim Übergang von Voxel v_n, auf Voxel v_n+1 wird ein Gewebeübergang festgestellt. Die durchstrahlte Länge des Gewebes A wird bei diesem Punkt berechnet und abgespeichert. Entsprechend wird weiter verfahren, bis das Objekt O durchdrungen ist.
Als Resultat dieser Simulation des Strahls S_y erhält man

– die durchstrahlten Gewebedicken d_A, d_B, d_C
– Klassifizierungen, d. h. eine Gewebezuordnung für die auf dem Strahl befindlichen Voxel
– die simulierte Projektion bzw. berechnete Intensität nach Durchdringen des Objekts

Dies wird für alle Strahlen aller Projektionen durchgeführt. Dabei kann die Klassifizierung beendet werden, wenn alle Voxel klassifiziert sind. Im vorliegenden Beispiel ist die Klassifizierung bereits mit der Simulation der Strahlen für eine der beiden dargestellten Projektionen vollständig.
Die (berechneten) Strahlintensitäten werden mit Hilfe der Dicken d_A, d_B, d_C ermittelt. Dabei könnte man das Gleichungssystem (1) in aller Allgemeinheit lösen. Ein besser durchführbares Vorgehen besteht in einer schrittweisen Lösung, die die Gewebetypen nacheinander behandelt. Dabei werden zunächst Strahlen selektiert, die nur einen Gewebetyp (im vorliegenden Fall A) durchstrahlen.
D. h. es wird nach Strahlen S gesucht, bei denen für die bestimmten Dicken gilt d_A > 0, d_B = 0, d_C = 0.
Die Formel (1) reduziert sich dann auf I_d = f(I_o, d_A, α_A,i) (2)
Mit der Gleichsetzung von gemessener und berechneter Intensität (I_d = I_y) erhält man die α_A,i als Ergebnis eines Optimierungsproblems für eine hinreichende Anzahl von Strahlen, für die (2) gilt. Als nächster Schritt wird nach Strahlen S gesucht, bei denen für die bestimmten Dicken gilt d_A ≥ 0, d_B > 0, d_C = 0.
Die Formel (1) reduziert sich dann auf I_d = f(f(I_o, d_A, α_A,i), d_B, α_B,i) (3)
Da die α_A,i schon bestimmt sind, erhält man mit der Gleichsetzung von gemessener und berechneter Intensität eine Gleichungssystem mit lediglich den α_B,i als Unbekannten. Auf diese Weise wird die Intensitätsabschwächung nacheinender für die einzelnen Gewebetypen errechnet. Damit können alle Intensitäten I_d erhalten und für die Korrektur der gemessenen Werte (Schritt ... in 3) verwendet werden.
Das Vorgehen lässt sich erweitern auf den Fall von 4, d. h. bei Vorliegen von Metallartefakten, indem bei der Strahlsimulation Voxel, die einem Metall zugeordnet würden, durch Voxel eines benachbarten Gewebetyps ersetzt werden.
In 4 sind die Strahlen so gewählt, dass die Voxel auf den Strahlen liegen. Dies ist nicht notwendig, wie anhand von 5 verdeutlicht wird.
In 5 ist das relevante Volumen in Würfel aufgeteilt. Exemplarisch ist einer der Würfel W hervorgehoben. Diese Würfel sind dabei so gewählt, dass jedem der Würfel ein Voxel und damit ein Schwächungskoeffizient μ zugeordnet ist. Über den Schwächungskoeffizienten kann jeder Würfel einer Gewebeklasse zugeordnet werden (Segmentierung). Die Dicke einer durchstrahlten Gewebeklasse ergibt sich dann als Summe der in den passierten Würfeln durchdrungenen Längen d, die wiederum aus Eintrittspunkt E und Austrittspunkt A erhalten werden.
6 und 7 kontrastieren in schematischer Weise die herkömmliche (6) und die erfindungsgemäße (7) Vorgehensweise.
Bei der Anwendung eines einfachen Segmentierungsverfahrens, wie zum Beispiel eines Schwellwert-basierten Segmentierungsverfahrens, wird demzufolge die Segmentierung erst während der Projektionsberechnung durchgeführt, anstelle das Objekt vor der Projektionsberechnung zu segmentieren. Bei dem geänderten Projektor wird für den betrachteten Voxel und Strahl das benötigte Gewicht bestimmt. Anschließend wird mittels Schwellwertverfahren anhand des Werts des Pixels entschieden, auf welchen Projektionswert einer Gewebeklasse der gewichtete Wert des Pixels addiert wird. Außerdem kann die durchstrahlte Materiallänge entsprechend aktualisiert werden.
Sollte es zudem, wie zum Beispiel im Metallartefaktkorrekturverfahren nach [Bal06], nötig werden, eine zusätzliche Modellprojektion durch einen veränderten Volumeninhalt zu berechnen, können immer noch die erforderlichen Projektionswerte mittels eines einzigen Projektionsschritts vorgenommen werden. Dafür ist es allerdings notwendig, dass die einzelnen Klassen oder Werte im Volumen in einer speziellen Form codiert werden, die dann im Projektionsschritt decodiert werden kann. Soll zum Beispiel wie dargestellt, die Modellprojektion, in der die Volumenwerte in Metallregionen durch Werte von benachbarten Volumenwerten ersetzt wurden, und die Metallprojektion selbst berechnet werden, können die Werte der Metallregionen durch abgeänderte Werte der Modellklassen, ersetzt werden. Dies kann zum Beispiel durch das Setzen eines Bits oder durch die Addition eines Offsets erreicht werden. Dadurch ist es dann möglich, im Projektionsschritt selbst zu entscheiden, bei welchen Projektionsbildern die Aktualisierungen durchzuführen sind. Die entsprechenden Werte der einzelnen Gewebeklassen (Metall, zu Knochen korrigiertes Metallgewebe, etc.) sind aus der gegebenen Codierung zu erhalten.
Sollte ein komplexeres Segmentierungsverfahren erforderlich sein, das nicht während der Projektionsberechnung durchführbar ist, kann trotzdem die Projektionsberechnung der einzelnen Gewebeklassen und/oder Modelle in einem Schritt durchgeführt werden, indem erst – wie oben beschrieben – während der Projektionsberechnung anhand der Codierung des Voxelwerts entschieden wird, auf welchem Projektionsbild einer Gewebeklasse der gewichtete Wert addiert wird.
Durch die Verwendung einer Pipeline-Struktur kann die Entscheidung, auf welche der beiden Projektionen der gewichtete Pixelwert addiert wird, und die anschließende Addition zeitgleich mit der Berechnung des nächsten Gewichts erfolgen.
Literatur:

Bu04 Buzug: Einführung in die Computertomographie: 1. Auflage 2004. Springer. ISBN 3-540-20808-9
Cas04 Van de Casteele: Model-based approached for Beam Hardening Correction and Resolution Measurements in Microtomography, Dissertation, Universiteit Antwerpen, 2004
Bal06 Bal and Spies: Metal artefact reduction in CT using tissue-class Modelling and adaptive prefiltering, Med. Phys. 33(8), 2006
Sid85 Siddon: Fast calculation of the exact radiological path for a threedimensional CT array, Medical Physics, 12(2): 252–255, March 1985
Mue98 Mueller: Fast and accurate three-dimensional reconstruction from Cone-Beam projection data using Algebraic Methods, Dissertation, Ohio State Univ., 1998

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 Patentliteratur

US 6771733 B2 [0030]

Zitierte Nicht-Patentliteratur

Buz04 [0005]
Cas04 [0009]
Bal06 [0010]
Sid85 [0011]
Mue98 [0012]
Cas04 [0030]
Cas04 [0031]
Cas04 [0031]
Cas04 [0032]
Cas04 [0032]
Bal06 [0034]
Bu04 [0041]
Cas04 [0044]
Bal06 [0065]

Claims

Verfahren zur Ermittlung von bzgl. Polychromieartefakten korrigierten Absorptionskoeffizienten für ein aus einer Mehrzahl von bzgl. der Absorptionseigenschaften sich unterscheidenden Materialtypen zusammengesetztes Objekt, bei dem – eine Vielzahl von Röntgenstrahl-Projektionen aus unterschiedlichen Positionen für das Objekt aufgenommen wird, – durch ein Rekonstruktionsverfahren mittels der aufgenommenen Projektionen ein erster Satz von Absorptionskoeffizienten ermittelt wird, – mittels Vorwärtsprojektion Projektionen berechnet werden, – die aufgenommenen Projektionen mit Hilfe der berechneten Projektionen im Sinne von mit monochromen Röntgenstrahlen aufgenommenen Projektionen korrigiert werden, und – durch das Rekonstruktionsverfahren mittels der korrigierten Projektionen ein zweiter, bzgl. Polychromieartefakten korrigierter Satz von Absorptionskoeffizienten ermittelt wird, wobei – eine formelmäßige Beschreibung einer Polychromie berücksichtigende Vorschrift verwendet wird, welche zu bestimmende Parameter umfasst, – die Vorschrift für die Berechnung von Projektionen mittels Vorwärtsprojektion verwendet wird, und – die Bestimmung der Parameter im Zuge der Berechnung von Projektionen mittels Vorwärtsprojektion erfolgt.
Verfahren nach Anspruch 1, dadurch gekennzeichnet, dass eine Segmentierung eines Volumens des Objektes im Zuge der Berechnung von Projektionen mittels Vorwärtsprojektion erfolgt.
Verfahren nach einem der Ansprüche 1 oder 2, dadurch gekennzeichnet, dass eine Korrektur für die Verhinderung von Metallartefakten im Zuge der Berechnung von Projektionen mittels Vorwärtsprojektion erfolgt.
Vorrichtung, welche zur Durchführung eines Verfahrens nach einem der Ansprüche 1 bis 3 ausgestaltet ist.
Computerprogramm, welches ein Verfahren nach einem der Ansprüche 1 bis 3 implementierenden Code umfasst.