DE102010040041B3

DE102010040041B3 - Verfahren zur Korrektur von durch zeitliche Veränderungen von Schwächungswerten auftretenden Artefakten

Info

Publication number: DE102010040041B3
Application number: DE102010040041A
Authority: DE
Inventors: Andreas Fieselmann; Frank Dennerlein
Original assignee: Siemens AG
Current assignee: Siemens Healthcare GmbH
Priority date: 2010-08-31
Filing date: 2010-08-31
Publication date: 2012-01-26
Anticipated expiration: 2030-09-01
Also published as: US8588497B2; US20120051619A1

Abstract

Die Erfindung betrifft ein Verfahren zur Korrektur von durch zeitliche Veränderungen von Schwächungswerten während der Rotationsaufnahme ein- oder zweidimensionaler Röntgen-Projektionsbilder mit einer Winkelgeschwindigkeit in einem zwei- oder dreidimensionalen, insbesondere mittels einer gefilterten Rückprojektion rekonstruierten Bilddatensatz eines Zielbereichs auftretenden Artefakten, wobei ein linearer, analytisch hergeleiteter, filterartiger Zusammenhang zwischen den Schwächungswerten des Bilddatensatzes zu einer Referenzzeit und den realen Schwächungswerten verwendet wird, in dem sich die Schwächungswerte des Bilddatensatzes aus einer Summe jeweils einer Anwendung eines winkelgeschwindigkeitsabhängigen Wichtungsfaktors, einer Punktbildfunktion und einer zeitlichen Ableitung, ausgewertet zu der Referenzzeit, auf die realen Schwächungswerte ergibt, wobei der lineare Zusammenhang invertiert wird und zur Korrektur der invertierte lineare Zusammenhang auf die Schwächungswerte des Bilddatensatzes angewendet wird.

Description

Die Erfindung betrifft ein Verfahren zur Korrektur von durch zeitliche Veränderungen von Schwächungswerten während der Rotationsaufnahme ein- oder zweidimensionaler Röntgen-Projektionsbilder mit einer Winkelgeschwindigkeit in einem zwei- oder dreidimensionalen, insbesondere mittels einer gefilterten Rückprojektion rekonstruierten Bilddatensatz eines Zielbereichs auftretenden Artefakten.
Es ist in der Röntgenbildgebung seit langem bekannt, zweidimensionale Bilddatensätze aus eindimensionalen Projektionsbildern oder dreidimensionale Bilddatensätze aus zweidimensionalen Projektionsbildern zu rekonstruieren, wofür verschiedene Rekonstruktionsverfahren bekannt sein, beispielsweise die gefilterte Rückprojektion und die iterative Rekonstruktion. Die Projektionsbilder werden dabei zeitlich aufeinander folgend unter unterschiedlichen Projektionswinkeln genommen, wobei Röntgenstrahler und Detektor meist eine Kreisbahntrajektorie beschreiben.
Ein Problem tritt immer dann auf, wenn sich während des Aufnahmevorgangs, also von der Aufnahme des ersten Projektionsbildes bis zur Aufnahme des letzten Projektionsbildes, die Schwächungswerte an bestimmten Punkten des Zielbereichs zeitlich verändern. Dies sei anhand eines Beispieles näher erläutert.
Der Schlaganfall ist die dritthäufigste Todesursache in westlichen Ländern. zur Behandlung kann die intra-arterielle Thrombolyse verwendet werden, bei der durch einen Katheter lokal in einer Gehirnarterie ein Medikament verabreicht wird, mit Hilfe dessen der Thrombus entfernt werden kann bzw. das diesen auflöst. Dabei ist es hilfreich, kurz vor und/oder während und/oder nach dem Eingriff eine tomographische Blutflussmessung durch Perfusionsbildgebung durchzuführen.
Dabei sind heute die Perfusions-CT-Bildgebung und die Perfusions-MR-Bildgebung bekannt, bei denen ein Kontrastmittel injiziert wird und eine Serie von Bilddatensätzen aufgenommen wird, um die Ausbreitung des Kontrastmittels im Gefäßsystem des Patienten zu beobachten. Es können aus den Bilddatensätzen dann Zeit-Schwächungskurven bestimmt werden, aus denen der Blutfluss, im Beispiel der zerebrale Blutfluss, ermittelt wird.
Da derartige Eingriffe häufig mit Unterstützung einer C-Bogen-Röntgeneinrichtung, beispielsweise eines Angiographiesystems, vorgenommen werden, wäre es für den Workflow bei den Eingriffen von großem Vorteil, wenn die Perfusionsbildgebung mit dem C-Bogen durchgeführt werden könnte.
Für Standard-Rekonstruktionsverfahren ist es jedoch eine Voraussetzung, dass die Schwächungswerte während der Aufnahme der Projektionsbilder wenigstens im Wesentlichen zeitlich konstant bleiben. Da die Rotationszeit des C-Bogens mit dem Röntgenstrahler und dem Detektor, die beispielsweise drei bis fünf Sekunden betragen kann, jedoch wesentlich länger ist als die einer CT-Einrichtung, welche beispielsweise 0,5 Sekunden betragen kann, kommt es bei den langsameren C-Bogen-Aufnahmen zu Artefakten, welche durch die Veränderung der Schwächungswerte durch das Einströmen beziehungsweise Ausströmen von Kontrastmittel begründet sind. Diese verfälschen aber die Ergebnisse.
Zur Lösung dieser Problematik sind im Stand der Technik bereits zwei Verfahren bekannt.
So wurde vorgeschlagen, einen dynamischen Rekonstruktionsansatz zu wählen, um einen Perfusionsbilddatensatz aus mit einer C-Bogen-Röntgeneinrichtung aufgenommenen Projektionsbildern zu rekonstruieren. Dieser Ansatz beruht auf einem ähnlichen Verfahren wie die algebraischen Rekonstruktionstechniken und ist in den Artikeln von Serowy, S. et al., „A Jacobi-like Solution to the Model Based Tomographic X-Ray Perfusion Imaging”, Proc. IEEE NSS/MIC, 2007, 4, Seiten 3085 bis 3088, sowie C. Neukirchen und S. Hohmann, „An Iterative Approach for Model-Based Tomographic Perfusion Estimation”. Proc. Fully-3D, 2007, Seiten 104 bis 107, beschrieben. Dieses Verfahren ist jedoch äußerst langsam und damit ineffizient. Zudem können die vorgeschlagenen Optimierungsverfahren Lösungen erzeugen, die nur lokalen Minima entsprechen und damit eine nicht zufrieden stellende Bildqualität liefern.
Weiter wurde vorgeschlagen, einen dynamischen Rekonstruktionsalgorithmus zu wählen, der aufbauend auf der gefilterten Rückprojektion die Inkonsistenz der Datensätze durch ein zeitliches Interpolationsverfahren reduziert. Dieses Verfahren ist in den Artikeln von P. Montes und G. Lauritsch, „A Temporal Interpolation Approach for Dynamic Reconstruction in Perfusion CT”, Med. Phys. 2007, 34, Seiten 3077 bis 3092, sowie A. Fieselmann et al., „A Dynamic Reconstruction Approach for Cerebral Blood Flow Quantification with an Interventional C-arm CT”, Proc. IEEE ISBI, 2010, Seiten 53 bis 56, beschrieben. Auch dieses Verfahren ist eher langsam und ist nicht optimal für typische Aufnahmegeometrien von C-Bogen-Angiographie-Anlagen geeignet, bei der der C-Bogen abwechselnd in unterschiedlichen Richtungen rotiert.
Die gattungsbildende Druckschrift CA 2 360 656 A1 schlägt vor, ein neues Rekonstruktionsverfahren auf der Theorie der pseudo-inversen Systemmatrizen zu basieren. Damit ist es möglich, mit jeder Projektion ein Update des zu rekonstruierenden tomographischen Bildes zu ermöglichen, genauso wie bei jedem Ereignis (PET). Modelle der physikalischen Prozesse der Bildaufnahme können durch die Systemmatrix genauso eingebracht werden wie Korrekturen für andere physikalische und instrumentelle Faktoren, beispielsweise im Falle der PET die Schwächung im Körper.
Ein Verfahren zur Kalibrierung von Computertomographen bei Zwei- oder Mehrspektrentomographie ist aus DE 10 2006 049 664 A1 bekannt. Im Rahmen der Materialzerlegung wird vorgeschlagen, die Inversionsformel durch eine Entwicklung nach Basisfunktionen zu bestimmen, wobei konkret einer Taylorentwicklung vorgeschlagen wird. Die dabei auftretenden Koeffizienten sollen durch Scannen eines kalibrierten Phantoms ermittelt werden.
Verfahren und Vorrichtungen zur Analyse von Ultraschallbildern sind aus US 2006/0078182 A1 bekannt. Dabei soll ein interessierendes Gebiet eines Bildes in einen Satz von Minisegmenten unterteilt werden, wobei für jedes Minisegment repräsentative Intensitätswerte für wenigstens zwei Arten von Substanzen ermittelt werden sollen. Auf diese Intensitätswerte soll eine Schwächungskorrektur angewandt werden, um eine Karte des interessierenden Bereichs bezüglich der Substanzen zu ermitteln.
Die US 2005/0135664 A1 betrifft ein Verfahren und eine Vorrichtungen zur Rekonstruktion von Volumendaten aus Projektionsdaten. Dabei soll ein Tomosynthese-Projektionsdatensatz gemeinsam mit zusätzlicher Information über das Objekt genutzt werden, um eine bestimmte Energiefunktion zu minimieren, so dass ein bestimmter Satz von Einschränkungen erhalten wird. Diese sollen dann zur Aktualisierung eines volumetrischen Bildes auf das rekonstruierte volumetrische Bild angewendet werden. Bei der Tomosynthese werden nur sehr wenige Projektionen verwendet. Endergebnis ist ein Datensatz, in dem jedes Voxel einer einzigen Materialklasse mit einer Komponente zugeordnet wird.
Aus der US 2004/0136608 A1 ist die algebraische Rekonstruktion von Bildern aus nicht den gleichen Abstand aufweisenden Daten bekannt. Dabei soll letztlich eine erste räumliche Verteilung von Datenpunkten auf eine zweite räumliche Verteilung von Datenpunkten abgebildet werden, wobei eine Matrix von Interpolationskoeffizenten, die die zweite räumliche Verteilung auf die erste räumliche Verteilung abbildet, invertiert wird.
Der Erfindung liegt daher die Aufgabe zugrunde, ein Verfahren anzugeben, mit dem eine verbesserte und insbesondere schnellere Korrektur von durch zeitliche Veränderungen der Schwächungswerte während der Aufnahme der Projektionsbilder hervorgerufenen Artefakten, insbesondere in der Perfusionsbildgebung, möglich ist.
Zur Lösung dieser Aufgabe ist bei einem Verfahren der eingangs genannten Art erfindungsgemäß vorgesehen, dass ein linearer, analytisch hergeleiteter, filterartiger Zusammenhang zwischen den Schwächungswerten des Bilddatensatzes zu einer Referenzzeit und den realen Schwächungswerten verwendet wird, in dem sich die Schwächungswerte des Bilddatensatzes aus einer Summe jeweils einer Anwendung eines winkelgeschwindigkeitsabhängigen Wichtungsfaktors, einer Punktbildfunktion und einer zeitlichen Ableitung, ausgewertet zu der Referenzzeit, auf die realen Schwächungswerte ergibt, wobei der lineare Zusammenhang invertiert wird und zur Korrektur der Artefakte der invertierte lineare Zusammenhang auf die Schwächungswerte des Bilddatensatzes angewendet wird.
Dabei sei vorab angemerkt, dass eine nullte Ableitung zur Referenzzeit wie bekannt die Auswertung der Funktion an der Referenzzeit zur Folge hat.
Der Erfindung liegt also die grundlegende Erkenntnis zugrunde, dass der Zusammenhang zwischen den Schwächungswerten des Bilddatensatzes und den realen Schwächungswerten analytisch hergeleitet und in einen invertierbaren Zusammenhang überführt werden kann. Dies ist für verschiedene Aufnahmegeometrien, beispielsweise die Fächerstrahlgeometrie, die Parallelstrahlgeometrie oder die Kegelstrahlgeometrie, möglich, indem beispielsweise die Zeitabhängigkeit über den Projektionswinkel der einzelnen, ja zeitlich hintereinander aufgenommenen Projektionsbilder verschlüsselt wird, dessen Zeitabhängigkeit ja wiederum bekannt ist, beispielsweise aus der Winkelgeschwindigkeit gefolgert werden kann. Schreibt man dann mithin die bekannte Formel für die Projektionen unter den verschiedenen Winkeln zeitabhängig und kombiniert diese Formel mit der ebenso bekannten, zeitabhängig formulierten Formel für die gefilterte Rückprojektion, so ergibt sich daraus eine Formel, die als ein Filtermodell aufgefasst werden kann, in der sich die rekonstruierten Schwächungswerte des Bilddatensatzes als räumliches Integral über eine zur Referenzzeit ausgewertete Beitragsfunktion ergeben, in der auch die realen Schwächungswerte enthalten sind.
Setzt man in der Beitragsfunktion eine Taylorentwicklung um die Referenzzeit, die beispielsweise als der mittlere Zeitpunkt des Aufnahmezeitraums der Projektionsbilder angesetzt werden kann, für die realen Schwächungswerte an, ergibt sich letztlich, da auch zeitliche Ableitungen eine lineare Operation sind, ein linearer Zusammenhang mit analytisch bestimmbaren Operatoren, hier der zeitlichen Ableitung, ausgewertet zur Referenzzeit, einem Wichtungsfaktor und einer Punktbildfunktion, wobei aufgrund der Taylorentwicklung immer eine Punktbildfunktion und ein Wichtungsfaktor einer zeitlichen Ableitung zugeordnet sind.
Es kann also der lineare Zusammenhang durch die Formel a = (Σ iW_iP_iD_i)b beschrieben werden, worin a ein die gemessenen Schwächungswerte beschreibender Vektor, b ein die realen Schwächungswerte beschreibender Vektor, W_i eine den Wichtungsfaktor angebende Diagonalmatrix, P_i eine die Punktbildfunktion beschreibende Matrix und D_i eine eine zeitliche Ableitung i-ter Ordnung beschreibende Matrix ist. Der Term nullter Ordnung kann sich dabei, wenn der zugehörige Gewichtungsfaktor eins ist, einfach als die Punktbildfunktion P₀ ergeben, was dann ja auch dem statischen Fall entspricht: Die Verschmierung wird durch die (statische) Punktbildfunktion bestimmt. Damit ergibt sich in dieser Schreibweise ein tatsächliches lineares Filtermodell.
Insgesamt können also die einzelnen Beiträge in ihrer Summe als eine einzige Matrix A beschrieben werden, so dass dann a = Ab gilt. Die Matrix A kann automatisch analytisch ermittelt und entsprechend auch, beispielsweise über einen bekannten Invertierungsalgorithmus, invertiert werden, b = A^–1a, so dass die realen Schwächungswerte zur Referenzzeit erhalten werden können. Dabei kann vorgesehen sein, dass der invertierte lineare Zusammenhang durch eine Singulärwertzerlegung ermittelt wird, wobei selbstverständlich auch andere Algorithmen eingesetzt werden können.
Wie viele Terme, insbesondere bezüglich der zeitlichen Ableitungen, tatsächlich berücksichtigt werden müssen, ergibt sich aus der Art der Zeitabhängigkeit, wobei üblicherweise die erste und ggf. die zweite Ableitung ausreichend sein können; dies wird für das Beispiel der Perfusion im Folgenden noch näher erläutert werden.
Erfindungsgemäß wird also zusammenfassend vorgeschlagen, die gefilterte Rückprojektion eines zeitlich dynamischen Objekts als ein Filter zu beschreiben, das sich mathematisch exakt herleiten lässt. Das heißt, wenn das Eingangssignal (beispielsweise die zeitliche Änderung der Kontrastmittelkonzentration in einer Arterie) genau bekannt ist, kann man mit Hilfe des durch den Zusammenhang beschriebenen Filters den rekonstruierten Bilddatensatz bestimmen. Das Filter lässt sich mathematisch exakt herleiten, was im Rahmen der Entwicklung der vorliegenden Erfindung gezeigt wurde.
Wie soeben gezeigt wurde, kann das Filter als lineares Gleichungssystem (Matrix A) beschrieben werden. Somit kann das Filter invertiert werden, um Rekonstruktionsartefakte zu eliminieren und/oder zu reduzieren.
Dabei kann mit besonderem Vorteil also die gefilterte Rückprojektion als Rekonstruktionsverfahren eingesetzt werden, welche bekanntlich sehr schnell arbeitet. Somit ist das Verfahren bereits schneller als die algebraischen Ansätze. Ein weiterer Vorteil der vorliegenden Erfindung ist, dass das Verfahren leicht an unterschiedliche Rotationsrichtungen, beispielsweise eines C-Bogens, angepasst werden kann.
In besonders vorteilhafter Ausgestaltung der vorliegenden Erfindung kann jedoch vorgesehen sein, dass die Korrektur lokal angewendet wird, insbesondere im Bereich eines Artefakts und/oder einer artefaktgefährdeten Struktur, insbesondere eines Blutgefäßes. Dies ist möglich, da der Faltungskern lokal begrenzt ist: Sowohl die Punktbildfunktion als auch die Ableitung sind lokale Operationen, so dass das erfindungsgemäße Verfahren in einem kleinen Unterbereich des Bilddatensatzes eingesetzt werden kann, was die Größe der zu invertierenden Matrizen reduziert und somit das Verfahren weiter stark beschleunigt und einen Vorteil gegenüber dem eingangs beschriebenen Stand der Technik bildet. Korrekturen können beispielsweise, indem Artefakte im Bilddatensatz identifiziert werden oder indem Vorwissen über artefaktgefährdete Bereiche genutzt wird, im Falle der Perfusionsbildgebung beispielsweise die Lage großer Blutgefäße, insbesondere von Arterien, speziell nur dort vorgenommen werden, wo sie benötigt sind, so dass die Gesamtrechenzeit stark reduziert wird. Solche lokalen Bereiche können in einem aus eindimensionalen Projektionsbildern rekonstruierten zweidimensionalen Bilddatensatz beispielsweise eine Größe von 40 mal 40 Pixeln bei 40 Zeitschritten, also 40 Projektionsbildern, haben.
Selbstverständlich ist es aber auch denkbar, dass der gesamte Bilddatensatz und/oder ein ausgewählter Bereich des Bilddatensatzes in lokal zu korrigierende Einheiten zerlegt werden. Dabei können Größen gewählt werden, die dem Wirkungsbereich insbesondere der Punktbildfunktion entsprechen und/oder die mit den verwendeten Algorithmen schnell und stabil berechnet werden können. Im dreidimensionalen Fall, worin also zweidimensionale Projektionsbilder zu einem dreidimensionalen Bilddatensatz rekonstruiert werden, hat sich gezeigt, dass die Korrektur beispielsweise schichtweise erfolgen kann, da beispielsweise bei zerebralen Blutflussmessungen in z-Richtung nur eine geringe Streuung vorliegt.
Ferner kann vorgesehen sein, dass die Lage der Struktur automatisch, insbesondere mittels eines anatomischen Atlasses, ermittelt wird. Wenn, wie beispielsweise bei der Perfusionsbildgebung, von Anfang an bekannt ist, dass Probleme beispielsweise bei der Anflutung in Arterien zu erwarten sind, so kann die Lage dieser Arterien als artefaktgefährdete Struktur aus einem anatomischen Atlas vorab bestimmt werden, so dass dann automatisch lokal um diese Strukturen korrigiert werden kann, ohne dass ein Eingriff eines Benutzers zur Auswahl der Bereiche notwendig ist.
Wie bereits erwähnt, liegt ein besonders vorteilhaftes Anwendungsgebiet des erfindungsgemäßen Verfahrens in der Perfusionsbildgebung, so dass vorgesehen sein kann, dass die Rotationsaufnahme im Rahmen einer Perfusionsbildgebung erfolgt, insbesondere mit einem C-Bogen. So können Artefakte reduziert werden und es ist beispielsweise eine genauere Bestimmung des Blutflusses möglich. Das erfindungsgemäße Verfahren erlaubt einen Einsatz der Perfusionsbildgebung auch an Röntgeneinrichtungen mit einem C-Bogen, die sich häufig ohnehin in einem Raum für Eingriffe, bei denen solche Informationen benötigt werden, befinden. Perfusionszeitkurven lassen sich häufig zeitlich lokal linearisieren oder wenigstens in zweiter Ordnung annähern (beispielsweise, wenn die Aufnahme das Erreichen der maximalen Anflutung und das langsame Abklingen einschließt), so dass zweckmäßigerweise vorgesehen sein kann, dass die Terme der Summe nur bis zur ersten Ableitung oder bis zur zweiten Ableitung betrachtet werden. So wird der Rechenaufwand weiter reduziert und das Verfahren kann schneller ausgeführt werden.
In weiterer vorteilhafter Weiterbildung der vorliegenden Erfindung kann vorgesehen sein, dass zur Ermittlung der einer n-ten Ableitung zugeordneten Punktbildfunktion ein durch ein Polynom n-ten Grades in der Zeit beschriebener realer Schwächungswertpunkt vorwärtsprojiziert wird und hieraus durch Rückprojektion die Punktbildfunktion bestimmt wird. Es ist also möglich, die Punktbildfunktion auch im Rahmen einer Simulationsmessung zu bestimmen, indem ein Punkt mit einem zeitlich veränderlichen Schwächungswert betrachtet wird. Dieser kann im Rahmen einer Vorwärtsprojektion genutzt werden, um virtuelle Projektionsbilder zu erhalten, welche zum Erhalten der Punktbildfunktion dann wie üblich rekonstruiert werden. Dabei ist selbstverständlich auf die Normalisierung zu achten. Der Vorteil einer solchen Vorgehensweise ist, dass der analytisch ermittelte Zusammenhang von einer perfekten Kreisbahn ausgeht, die im Realfall nicht immer vorliegt. Die Abweichungen sind aber bekannt, so dass sie im Rahmen einer derartigen Simulation mitberücksichtigt werden können, so dass das Verfahren noch exakter wird. Konkret kann zur Ermittlung der Punktbildfunktion nullter Ordnung beispielsweise ein konstanter zeitlicher Schwächungswertverlauf angesetzt werden, für die erste Ordnung ein linearer etc.
Es sei an dieser Stelle noch allgemein angemerkt, dass das erfindungsgemäße Verfahren auch im Rahmen der Perfusions-Tomosynthese vorteilhaft angewandt werden kann, wo eine Rekonstruktion, vorzugsweise einer nicht-axialen Schicht, aus einem begrenzten Projektionswinkelintervall durchgeführt wird. Auch ist es anwendbar bei Rekonstruktionen, bei denen die zur Verfügung stehenden Projektionen mit einem größeren Winkelunterschied, beispielsweise von 5° statt 1°, aufgenommen werden. Auch in solchen Fällen können die Punktbildfunktionen nullter und höherer Ordnung, die in dem Filtermodell verwendet werden, berechnet werden.
Weitere Vorteile und Einzelheiten der vorliegenden Erfindung ergeben sich aus den im Folgenden beschriebenen Ausführungsbeispielen sowie anhand der Zeichnung. Dabei zeigen:
1 eine schematische Darstellung des Filtermodells,
2 eine Darstellung zur Ermittlung der Punktbildfunktion aus einer Simulation und
3 ein mögliches Artefakt und seine Reduzierung.
Das erfindungsgemäße Verfahren wird im Folgenden der einfacheren Darstellung halber für den Fall einer Rekonstruktion eines zweidimensionalen Bilddatensatzes aus eindimensionalen Projektionsbildern, die unter verschiedenen, auf einer Kreisbahn liegenden Projektionsgeometrien, also Projektionswinkeln, aufgenommen werden, erläutert. Selbstverständlich kann das hier Beschriebene auch auf die Rekonstruktion dreidimensionaler Bilddatensätze aus zweidimensionalen Projektionsbildern angewandt werden, wobei beispielsweise der Feldkamp-Algorithmus der gefilterten Rückprojektion angesetzt werden kann.
Die Herleitung des genutzten Filtermodells sei zunächst für den Fall der Fächerstrahlgeometrie skizziert.
Es wird davon ausgegangen, dass der Röntgenstrahler mit einer konstanten Winkelgeschwindigkeit ω_s auf einer kreisförmigen Bahn des Radius R um den Ursprung des Koordinatensystems rotiert. Die Position a(λ(t)) der Quelle zum Zeitpunkt t kann dann geschrieben werden als: a(λ(t)) = (Rcos(λ(t)), Rsin(λ(t)))^T (1) λ(t) = ω_s·t + λ₀ (2) worin λ(t) die Winkelposition des Röntgenstrahlers zum Zeitpunkt t ist und λ₀ der Winkel zum Zeitpunkt t = 0. Für zeitabhängige Schwächungswerte μ(x, λ(t)) an Positionen x = (x, y)^T lassen sich die Projektionen p(λ(t), γ) unter dem Fächerwinkel γ mit der üblichen Definition der δ-Funktion schreiben: p(λ(t), γ) = ∫∫ +∞ / –∞μ(x, λ(t))·δ(γ*(x, λ(t)) – γ)dxdy (3), worin γ*(x, λ) den Fächerwinkel bei der Röntgenstrahlungsquellenposition λ definiert, der sich mit dem Pixel an der Stelle x schneidet. Dieser Winkel kann berechnet werden zu
worin die Einheitsvektoren e_u = (–sin(λ), cos(λ))^T und e_w = (cos(λ), sin(λ))^T verwendet werden. Der Schwächungswert μ^r(r, t_r) an einem Ort r bei der Referenzzeit t_r kann durch gefilterte Rückprojektion geschrieben werden:
worin h_ramp(γ) der Rampenfilter ist. Die Funktion w(λ, γ) ist eine Gleitfensterfunktion, die auch redundante Daten in der Fächerstrahlgeometrie kompensiert. Die Funktion w(λ, γ) ist Null außerhalb eines Winkelintervalls der Länge λ_D und nimmt Funktionswerte m(λ, γ) innerhalb dieses Intervalls an. Die Funktion m(λ, γ) kann beispielsweise wie in dem Artikel von F. Noo et al., „Image reconstruction from fan-beam projections an less than a short scan”, Phys. Med. Biol. 47, 2002, 2525–2546, gewählt werden.
Kombiniert man nun die Formeln (3) und (5), ändert die Reihenfolge der Integration und wertet die Faltung mit der Delta-Funktion aus, erhält man die Gleichungen
Die Gleichungen (7) und (8) können als ein Filter interpretiert werden, das die realen Schwächungswerte μ in die rekonstruierten Schwächungswerte μ^r des Bilddatensatzes transformiert. Die Funktion χ(x, r, t_r) kann als die Beiträge der Orte x zur Rekonstruktion bei r verstanden werden.
Nun wird die Taylor-Entwicklung von μ(x, λ(t)) um λ(t_r) berechnet, um die Eigenschaften von χ(x, r, t_r) in Bezug auf die Ableitungen von μ zu berechnen.
Dabei wird angenommen, dass μ(x, λ(t)) in dem Intervall als eine Taylorreihe ausgedrückt werden kann, das dem aktuellen Gleitfenster w entspricht. Da die zeitlichen Veränderungen, insbesondere im Rahmen der Perfusion, in der Praxis kontinuierliche, niederfrequente Funktionen sind, ist dies möglich. Betrachtet man die Gleichung (2), so erkennt man, dass Ableitungen zweiter und höherer Ordnung von λ(t) verschwinden und man erhält:
Kombiniert man nun (8), (9) und (10) und ändert man die Reihenfolge von Summe und Integration, erhält man schließlich
Dabei kann P_n(x, r, λ(t_r)) als die Punktbildfunktion interpretiert werden, die mit der n-ten zeitlichen Ableitung von μ(r, t_r) und der inversen Winkelgeschwindigkeit ω_s des Röntgenstrahlers gewichtet wird. Als Winkel λ(t_r) wird dabei in diesem Ausführungsbeispiel zweckmäßigerweise das Zentrum des Gleitfensters gewählt.
Eine Darstellung des so erarbeiteten Filtermodells zeigt 1, wobei ein mögliches Aussehen der Punktbildfunktion P_n mit angedeutet wurde.
Eine ähnliche Herleitung ist auch für andere Geometrien möglich, wobei das Filtermodell grundsätzlich, setzt man als a und b jeweils einen die Schwächungswerte des Bilddatensatzes bzw. die realen Schwächungswerte enthaltenden Vektor an, als lineares Gleichungssystem geschrieben werden kann: a = (Σ iW_iP_iD_i)b = Ab (13), worin W_i eine die Wichtung enthaltende Diagonalmatrix ist, P_i die die jeweilige Punktbildfunktion beschreibende Matrix und D_i eine die Differentiation n-ten Grades nach der Zeit beschreibende Matrix (beispielsweise eine Differenzenbildung zwischen zeitlich benachbarten Punkten).
Zur Korrektur von Artefakten wird nun der Zusammenhang (13) genutzt, indem die Matrix A, beispielsweise durch eine Singulärwertzerlegung, invertiert wird.
Die Punktbildfunktionen können dabei analytisch berechnet werden, es kann jedoch auch eine Art „Vermessung” in einer Simulation stattfinden, um Abweichungen von einer idealen Kreisbahn berücksichtigen zu können. Dies ist schematisch in 2 für P₀ und P₁ dargestellt.
Ausgegangen wird von einem Punkt im Raum, insbesondere im Ursprung eines Koordinatensystems, dessen Schwächungswert einem bestimmten zeitlichen Verlauf folgt. Für die nullte Ordnung wird ein konstanter zeitlicher Verlauf 1 angenommen, beispielsweise bei Eins, für die erste Ordnung ein linearer Anstieg 2, beispielsweise mit Steigung Eins.
Für diese Punkte mit ihren Verläufen wird nun in einem Schritt 3 vorwärtsprojiziert, das heißt, es werden virtuelle Projektionsbilder ermittelt, die dann in einem Schritt 4 mittels gefilterter Rückprojektion zu einem Bilddatensatz rekonstruiert werden, der den Punktbildfunktionen entspricht. In den Schritten 3 können dabei Abweichungen von einer idealen Kreisbahngeometrie berücksichtigt werden.
Selbstverständlich lassen sich auch Punktbildfunktionen höherer Ordnung so bestimmen, wobei aber grundsätzlich in allen Fällen die Normalisierung beachtet werden muss.
Ein konkreter Ablauf des erfindungsgemäßen Verfahrens kann sich dann beispielsweise folgendermaßen gestalten. Zunächst werden die die Punktbildfunktionen beschreibenden Matrizen bestimmt, sei es rechnerisch oder über eine Simulation nach 2, und beispielsweise für mehrere Korrekturvorgänge vorgehalten.
Nach Aufnahme und Rekonstruktion der Bilddatensätze erfolgt dann die Korrektur, wobei mehrere Fälle denkbar sind.
Wie auch beispielsweise aus 1 bereits ersichtlich war, sind alle linearen Operationen, insbesondere auch die Punktbildfunktion, auf den Schwächungswerten lokal, so dass mit besonderem Vorteil die Korrektur lokal angewandt werden kann, insbesondere in Bereichen um Artefakte herum bzw. um artefaktgefährdete Strukturen herum. Die Bereiche können automatisch und/oder manuell definiert werden, wobei eine automatische Ermittlung, beispielsweise durch Artefaktdetektionsalgorithmen, bevorzugt aber durch Ermittlung der Lage artefaktgefährdeter Strukturen, beispielsweise im Falle der Perfusionsmessung größerer Blutgefäße, insbesondere anhand eines anatomischen Atlasses, erfolgt.
Beispielhaft sei dazu auf 3 verwiesen. Dort zeigt das linke Teilbild 5 angedeutet einen Schnitt durch den Kopf 6, wie er in einem anatomischen Atlas vorkommen kann. Eine Arterie 7 ist beispielhaft gezeigt. Es kann nun angenommen werden, dass es in einem Bereich 8, der beispielsweise 40 mal 40 Pixel groß sein kann, um die Arterie 7 zu Artefakten kommen kann, wenn eine Perfusionsmessung, beispielsweise zur Bestimmung des Blutflusses im Gehirn, durchgeführt werden soll.
Einen rekonstruierten zweidimensionalen Bilddatensatz 9 in einer dem Schnitt des Teilbildes 5 entsprechenden Schicht, der noch nicht korrigiert ist, zeigt das mittlere Teilbild. Ersichtlich ist die Arterie 7 nicht zu erkennen, da sie von Artefakten 10 aufgrund der Änderung der Kontrastmittelkonzentration während der Aufnahme der eindimensionalen Projektionsbilder, die beispielsweise drei Sekunden andauern kann, umgeben ist.
Wird nun aber die Matrix A für den Bereich 8 bestimmt, wobei es vorliegend ausreichend sein kann, die Ableitungen nur bis zur ersten oder zur zweiten Ordnung zu betrachten, da sich Perfusionszeitverläufe oft linear oder quadratisch bereits hinreichend genau annähern lassen, und wird dann die invertierte Matrix A^–1 auf den Bereich 8 in dem Bilddatensatz 9 angewandt, Pfeil 11, so erhält man den als rechtes Teilbild dargestellten korrigierten Bilddatensatz 12, in dem die Signale 13 der Arterie 7 klar zu erkennen sind, da die umgebenden Artefakte 10 reduziert oder gar gänzlich beseitigt sind.
Sind die Bereiche, in denen Artefakte vorliegen oder vorliegen können, nicht genau definiert, so kann alternativ vorgesehen sein, dass der gesamte Bilddatensatz 9 in Bereiche aufgeteilt wird, die dann jeder für sich lokal korrigiert werden, um eine schnellere Berechnung zu ermöglichen.
Selbstverständlich ist es grundsätzlich auch möglich, das gesamte Bild auf einmal zu korrigieren.
Wird grundsätzlich mit den selben Bildaufnahmeparametern gearbeitet und werden grundsätzlich Ableitungen bis zur selben Ordnung betrachtet oder lassen sich die Anwendungsfälle in entsprechende Gruppen einteilen, können die invertierte Matrix A^–1 bzw. die mehreren invertierten Matrizen A^–1 dann selbstverständlich auch dauerhaft vorgehalten werden, beispielsweise in einer Speichereinheit der Röntgeneinrichtung, um dann direkt angewendet zu werden.

Claims

Verfahren zur Korrektur von durch zeitliche Veränderungen von Schwächungswerten während der Rotationsaufnahme ein- oder zweidimensionaler Röntgen-Projektionsbilder mit einer Winkelgeschwindigkeit in einem zwei- oder dreidimensionalen, insbesondere mittels einer gefilterten Rückprojektion rekonstruierten Bilddatensatz eines Zielbereichs auftretenden Artefakten, dadurch gekennzeichnet, dass ein linearer, analytisch hergeleiteter, filterartiger Zusammenhang zwischen den Schwächungswerten des Bilddatensatzes zu einer Referenzzeit und den realen Schwächungswerten verwendet wird, in dem sich die Schwächungswerte des Bilddatensatzes aus einer Summe jeweils einer Anwendung eines winkelgeschwindigkeitsabhängigen Wichtungsfaktors, einer Punktbildfunktion und einer zeitlichen Ableitung, ausgewertet zu der Referenzzeit, auf die realen Schwächungswerte ergibt, wobei der lineare Zusammenhang invertiert wird und zur Korrektur der Artefakte der invertierte lineare Zusammenhang auf die Schwächungswerte des Bilddatensatzes angewendet wird.
Verfahren nach Anspruch 1, dadurch gekennzeichnet, dass der lineare Zusammenhang durch die Formel a = (Σ iW_iP_iD_i)b beschrieben wird, worin a ein die gemessenen Schwächungswerte beschreibender Vektor, b ein die realen Schwächungswerte beschreibender Vektor, W_i eine den Wichtungsfaktor angebende Diagonalmatrix, P_i eine die Punktbildfunktion beschreibende Matrix und D_i eine eine zeitliche Ableitung i-ter Ordnung beschreibende Matrix ist.
Verfahren nach Anspruch 1 oder 2, dadurch gekennzeichnet, dass der invertierte lineare Zusammenhang durch eine Singulärwertzerlegung ermittelt wird.
Verfahren nach einem der vorangehenden Ansprüche, dadurch gekennzeichnet, dass die Korrektur lokal angewendet wird, insbesondere im Bereich eines Artefakts und/oder einer artefaktgefährdeten Struktur, insbesondere eines Blutgefäßes.
Verfahren nach Anspruch 4, dadurch gekennzeichnet, dass der gesamte Bilddatensatz und/oder ein ausgewählter Bereich des Bilddatensatzes in lokal zu korrigierende Einheiten zerlegt werden.
Verfahren nach Anspruch 5, dadurch gekennzeichnet, dass bei einem dreidimensionalen Bilddatensatz die Korrektur schichtweise erfolgt.
Verfahren nach einem der Ansprüche 4 bis 6, dadurch gekennzeichnet, dass die Lage der Struktur automatisch, insbesondere mittels eines anatomischen Atlasses, ermittelt wird.
Verfahren nach einem der vorangehenden Ansprüche, dadurch gekennzeichnet, dass die Rotationsaufnahme im Rahmen einer Perfusionsbildgebung erfolgt, insbesondere mit einem C-Bogen.
Verfahren nach Anspruch 8, dadurch gekennzeichnet, dass die Terme der Summe nur bis zur ersten Ableitung oder bis zur zweiten Ableitung betrachtet werden.
Verfahren nach einem der vorangehenden Ansprüche, dadurch gekennzeichnet, dass zur Ermittlung der einer n-ten Ableitung zugeordneten Punktbildfunktion ein durch ein Polynom n-ten Grades in der zeit beschriebener realer Schwächungswertpunkt vorwärtsprojiziert wird und hieraus durch Rückprojektion die Punktbildfunktion bestimmt wird.