DE102007040519B4 - Verfahren zur Reduktion des Bildrauschens im Rahmen der Aufnahme eines Bildes mit zwei unterschiedlichen Röntgenspektren - Google Patents

Abstract

Verfahren zur Reduktion des Bildrauschens im Rahmen der Aufnahme wenigstens eines strahlungsbasierten Bildes eines Bildaufnahmebereichs mit zwei unterschiedlichen Röntgenspektren mit den Schritten – Aufnahme von Rohbildern des Bildaufnahmebereichs mit den beiden unterschiedlichen Röntgenspektren mit jeweils paarweise einander zugeordneten Messwerten und – zur Separation unterschiedlicher Materialien im Bildaufnahmebereich Anwendung wenigstens eines einen Übergang von einem Messwertepaar in ein zugeordnetes Rekonstruktionswertepaar beschreibenden Inversionsoperators mit einer integrierten Rauschfilterung auf die aufgenommenen Rohbilder, wobei die Rauschfilterung in den Inversionsoperator integriert wird, indem dieser für jeweils mit den unterschiedlichen Röntgenspektren aufgenommene geglättete Messwertepaare der Rohbilder bestimmt wird, wobei als Inversionsoperator wenigstens eine Inversionsmatrix bestimmt wird zu jedem geglätteten Messwertepaar der Rohbilder, wobei eine Koeffizientenmatrix zur Bildung der Inversionsmatrix in Abhängigkeit von einem zugehörigen geglätteten Messwertepaar invertiert wird.

Description

• Die Erfindung betrifft ein Verfahren zur Reduktion des Bildrauschens im Rahmen der Aufnahme wenigstens eines strahlungsbasierten Bildes eines Bildaufnahmebereichs mit zwei unterschiedlichen Röntgenspektren.
• Bei der sogenannten Zweispektren-Projektionsradiographie bzw. der Dual-Energie-Projektionsbildgebung, oft auch abgekürzt als DE-Bildgebung bezeichnet, wird ein zu untersuchendes Objekt, beispielsweise ein Patient, mit zwei unterschiedlichen Röntgenspektren aufgenommen, um so zwei Projektionsbilder (Rohbilder) des Aufnahmebereichs zu erzeugen.
• Durch eine geeignete Kombination der beiden Rohbilder können so zwei radiologisch unterschiedliche Materialien wie beispielsweise Weichteilgewebe und Knochen voneinander getrennt werden.
• Im Rahmen einer weit verbreiteten eher qualitativen Methode der Dual-Energie-Projektionsbildgebung werden nur Grauwertbilder erzeugt, während bei der quantitativen Dual-Energie-Bildgebung physikalische Größen rekonstruiert werden, also Rekonstruktionswerte wie Materialdicken (in cm) oder Massenbelegungsflächendichten (in g/cm²), im Folgenden auch kurz Massenbelegungen oder Massenbelegungsdichten, aus den Rohbilddaten bestimmt werden.
• Ein bekannter Nachteil der DE-Bildgebung ist die starke Zunahme des Bildrauschens im Vergleich zu den Rohbildern. Für eine exakte quantitative Rekonstruktion im Hinblick auf die Bestimmung von Rekonstruktionswerten ist ein mathematisch im Allgemeinen schlecht konditioniertes nichtlineares Gleichungssystem zu lösen. Damit ist aber eine Erhöhung des Bildrauschens verbunden.
• Deshalb werden bereits derzeit verschiedene Rauschfilterverfahren eingesetzt. Diese basieren jedoch entweder nicht auf dem physikalisch korrekten nichtlinearen Modell bzw. stellen lediglich eine Bildnachverarbeitung dar, bei der auf geschickte Weise negative Korrelationen zwischen separierten Materialbildern zum Zweck der Rauschreduktion ausgenutzt werden. Derartige Verfahren sind beispielsweise aus den Artikeln „Quantitative evaluation of noise reduction strategies in dual-energy imaging” von R.J. Warp und J.T. Dobbins aus Med. Phys. 30 (2), Feb 2003, und „An Algorithm for noise suppression in Dual Energy CT Material Density Images” von W.A. Kalender, E. Klotz und L. Kostaridou aus IEEE Trans. Med Imaging, Vol. 7, No. 3, September 1988, 218–224, sowie „A correlated noise reduction algorithm for dual-energy digital subtraction angiography” von C.H. McCollough, M.S. VanLysel, W.W. Peppler und C.A. Mistretta aus Med. Phys. 16 (6), Nov/Dec 1989, 873–880, bekannt.
• Diese Rauschfilterungen bzw. Verfahren zur Reduktion des Bildrauschens sind somit ansatzbedingt nicht optimal.
• Aus der Patentschrift US 4 792 900 A ist es bekannt, materialspezifische Filterfunktionen auf Röntgenbilddaten anzuwenden, wobei die Röntgenbilddaten zwei bei verschiedenen Röntgenenergien aufgenommene Bilder umfassen, um anschließend aus den gefilterten Röntgenbilddaten materialspezifische Bilddaten zu erzeugen.
• Aus der Patentschrift US 5 115 394 A ist es bekannt, eine tabellierte Funktion zu erzeugen und zu verwenden, um ausgehend von Röntgenintensitäten bei zwei verschiedenen Röntgenenergien eine photoelektrische Komponente und eine Compton-Komponente zu berechnen. Weiterhin ist es aus dieser Patentschrift bekannt, diese einzelnen Komponenten anschließend zu filtern, um das Signal-Rausch-Verhältnis zu erhöhen.
• Aus den Seiten 375–377 der Druckschrift „Imaging Systems for Medical Diagnostics” (Hrsg. Erich Krestel, Berlin, 1990) sind die allgemeinen Grundlagen der Bildgebung mittels zweier verschiedener Röntgenenergien bekannt.
• Aus der Patentschrift US 7 010 092 B2 ist es bekannt, Röntgenbildgebung bei zwei verschiedenen Röntgenenergien durchzuführen, indem eine Röntgenstrahlung nach dem Durchgang durch das Untersuchungsobjekt aufgespalten wird, und mit einem unterschiedlichen Filter energieabhängig gefiltert wird.
• Aus dem Zeitschriftenartikel „Generalized image combinations in dual KVP digital radiography” (Lehmann et al., Med. Phys. 8(5), 659 (1981)) ist eine nicht-lineare Transformation zur Bestimmung einer photoelektrischen Komponente und einer Compton-Komponente aus Röntgenintensitäten bei zwei verschiedenen Röntgenenergien bekannt. Weiterhin ist aus diesem Zeitschriftenartikel die Verwendung von Inversionsoperatoren bekannt.
• Aus dem Zeitschriftenartikel „Density and atomic number measurements with spectral x-ray attenuation method” (Heismann et al., J. Appl. Phys. 94, 2073 (2003)) ist ein Algorithmus zur Extraktion der Dichte und der mittleren Kernladungszahl mittels einer Röntgenbildgebung bei zwei verschiedenen Röntgenenergien bekannt.
• Der Erfindung liegt die Aufgabe zu Grunde, ein Verfahren zur Reduktion des Bildrauschens im Rahmen der Aufnahme wenigstens eines strahlungsbasierten Bildes eines Bildaufnahmebereichs mit zwei unterschiedlichen Röntgenspektren anzugeben, das diesbezüglich verbessert ist.
• Zur Lösung dieser Aufgabe ist ein Verfahren dieser Art vorgesehen, das die folgenden Schritte aufweist:
• – Aufnahme von Rohbildern des Bildaufnahmebereichs mit den beiden unterschiedlichen Röntgenspektren mit jeweils paarweise einander zugeordneten Messwerten und
• – zur Separation unterschiedlicher Materialien im Bildaufnahmebereich Anwendung wenigstens eines einen Übergang von einem Messwertepaar in ein zugeordnetes Rekonstruktionswertepaar beschreibenden Inversionsoperators mit einer integrierten Rauschfilterung auf die aufgenommenen Rohbilder, wobei die Rauschfilterung in den Inversionsoperator integriert wird, indem dieser für jeweils mit den unterschiedlichen Röntgenspektren aufgenommene geglättete Messwertepaare der Rohbilder bestimmt wird, wobei als Inversionsoperator wenigstens eine Inversionsmatrix bestimmt wird zu jedem geglätteten Messwertepaar der Rohbilder, wobei eine Koeffizientenmatrix zur Bildung der Inversionsmatrix in Abhängigkeit von einem zugehörigen geglätteten Messwertepaar invertiert wird.
• Es werden also zunächst im Rahmen der Durchführung einer Zweispektren-Projektionsradiographie bzw. einer Dual-Energie-Projektionsbildgebung mit zwei unterschiedlichen Röntgenspektren bzw. Strahlungsspektren zwei Projektionsbilder als Rohbilder erzeugt. Diese Rohbilder zeigen den gleichen Bildaufnahmebereich, so dass sich die jeweils zu identischen Pixeln der Bildaufnahmen gehörenden Messwerte als Paare einander zuordnen lassen.
• Um nun eine Trennung zweier radiologisch unterschiedlicher Materialien wie beispielsweise von Weichteilgewebe und Knochen im Bildaufnahmebereich bzw. in den Bildaufnahmen zu erhalten, wird ein Inversionsoperator auf diese aufgenommenen Rohbilder angewandt, der den Übergang von einem solchen Messwertpaar für identische Pixel der Rohbilder in ein zugeordnetes Rekonstruktionswertepaar, also in ein Paar zu rekonstruierender physikalischer Größen wie einer Materialdicke bzw. einer Massenbelegungsdichte und dergleichen, beschreibt.
• Dieser Inversionsoperator, der somit der Bestimmung der Rekonstruktionswerte dient, weist erfindungsgemäß eine integrierte Rauschfilterung auf. Die Filterung ist also in den Operator einbezogen bzw. an einen Inversionsoperator, mit dem die Rekonstruktionswerte bestimmt werden können, gekoppelt.
• Zur Vereinfachung ist hier jeweils von Messwertepaaren bzw. Rekonstruktionswertepaaren die Rede. Denkbar ist es aber ebenso, weitere Projektionsbilder mit jeweils nochmals anderen Energien aufzunehmen, bei denen sich dann wiederum für identische Pixel Messwerte einander zuordnen lassen. In diesem Fall handelt es sich nicht mehr um Messwertepaare, sondern allgemeiner um Messwertetupel, wobei ein entsprechender Inversionsoperator in diesem Fall den Übergang eines solchen Messwertetupels in ein Rekonstruktionswertetupel beschreiben würde. Auch eine Reduktion des Bildrauschens für solche Mehr-Energie-Projektionsbilder ist selbstverständlich von der Erfindung erfasst, auch wenn im Folgenden zur Vereinfachung stets auf den Fall von zwei Energien bzw. zwei unterschiedlichen Spektren für die Aufnahme eingegangen wird.
• Selbstverständlich muss die Aufnahme der Rohbilder nicht direkt im Vorfeld der Anwendung des Inversionsoperators mit der integrierten Rauschfilterung erfolgen, sondern die Rohbilder können bereits vorab aufgenommen sein bzw. als Aufnahmen vorliegen, die später, um abschließend das korrekte rauschgefilterte Bild zu erhalten, in der erfindungsgemäßen Art und Weise verarbeitet werden.
• Die Messwerte, die im Folgenden als p₁ und p₂ bezeichnet werden und die mit den beiden unterschiedlichen Spannungen der Röntgenröhre aufgenommen wurden, ergeben sich aus dem Quotienten der gemessenen Intensitäten vor und hinter einem aufzunehmenden Objekt im Bildaufnahmebereich im selben Detektorpixel durch die Formeln p₁ = –ln(I₁/I₁₀), p₂ = –ln(I₂/I₂₀), wobei I_i die geschwächten, I_i0 die ungeschwächten Intensitäten sind.
• Die logarithmierten Messwerte werden im Folgenden einfach als Messwerte bezeichnet.
• Um von einem Messwertepaar zu einem zugehörigen Rekonstruktionswertepaar zu gelangen, kann nun ein Inversionsoperator angegeben werden, der diesen Übergang beschreibt.
• Erfindungsgemäß ist die Rauschfilterung in den Inversionsoperator integriert, indem dieser für jeweils mit den unterschiedlichen Strahlungsspektren aufgenommene geglättete Messwertepaare der Rohbilder bestimmt wird. Hintergrund dieses Vorgehens ist, dass der Inversionsoperator wesentlich durch Materialeigenschaften und die Röntgenspektren bestimmt wird. Daher ist es nicht zweckmäßig bzw. in einem gewissen Sinne übertrieben, zufällige Messwertfehler, die beispielsweise lediglich durch das Quantenrauschen entstehen, als Schwankungen der Materialeigenschaften im Bildaufnahmebereich falsch zu interpretieren, indem der Inversionsoperator selbst auf Grund dieser Messwertfehler verändert wird. Daher geht die Erfindung davon aus, dass der Inversionsoperator nicht für die aufgenommenen Messwertepaare direkt, sondern für geglättete Messwertepaare der Rohbilder bestimmt wird. Die beiden Rohbilder für die niedrigere bzw. höhere Spannung am Röntgenstrahlungserzeuger werden also jeweils geglättet und diese geglätteten Daten für die Bestimmung des Inversionsoperators verwendet.
• Insbesondere können die Messwertepaare der Rohbilder adaptiv geglättet werden. Dabei handelt es sich um eine Kanten ausschließende Glättung. Selbstverständlich können ergänzend weitere hier nicht genannte Glättungsverfahren bzw. andere zweckmäßige Glättungsverfahren zur Glättung der Rohbilder herangezogen werden.
• Anschließend kann jeweils für jedes identische Pixelpaar der geglätteten Rohbilder mit einem zugeordneten geglätteten Messwertepaar ein Inversionsoperator bestimmt werden. Es wird also für jedes Messwertepaar, das durch die Messwerte zu jeweils identischen Pixeln der Rohbilder gebildet wird, ein Inversionsoperator bestimmt, der den Übergang zwischen den Messwerten und gewünschten Rekonstruktionswerten bestimmt. Da der Inversionsoperator für die geglätteten Messwertepaare angegeben wird, enthält dieser bereits eine in den Operator selbst einbezogene Rauschfilterung.
• Anschließend kann wenigstens ein für jeweils mit den beiden unterschiedlichen Strahlungsspektren aufgenommene geglättete Rohbilder bestimmter Inversionsoperator zur Bestimmung eines Rekonstruktionswertepaars auf ein entsprechendes tatsächliches Messwertepaar der ungeglätteten Rohbilder angewandt werden.
• Der Inversionsoperator, der für geglättete Rohdaten bestimmt wurde, wird also in Form einer Operatorglättung nunmehr zur Bestimmung der Rekonstruktionswertepaare, also beispielsweise einer Materialdicke und dergleichen, auf die ursprünglichen Rohbilddaten vor der Glättung angewandt. Auf diese Art und Weise werden rauschreduzierte Rekonstruktionswerte erhalten.
• Das Verfahren kann weitgehend automatisiert mit Hilfe eines Programmmittels beispielsweise auf einer Steuerungseinrichtung einer entsprechenden Röntgenaufnahmeeinrichtung ablaufen.
• Erfindungsgemäß wird als Inversionsoperator wenigstens eine Inversionsmatrix bestimmt, insbesondere eine 2×2-Matrix, zu jedem geglätteten Messwertepaar der Rohbilder. Die Bestimmung einer 2×2-Matrix zur Beschreibung des Übergangs von den Messwerten zu Rekonstruktionswerten ist dann zweckmäßig, wenn wie im Regelfall Paare von Messwerten vorliegen, also Projektionsbilder mit zwei unterschiedlichen Spektren bzw. im Rahmen einer Dual-Energie-Aufnahme aufgenommen wurden.
• Jeweils einem mit den beiden unterschiedlichen Strahlungsspektren aufgenommenen Paar von Messwerten kann als wenigstens ein Rekonstruktionswertepaar ein Paar von Massenbelegungsdichten und/oder Materialdicken zugeordnet werden. Dabei werden die Materialdicken in der Regel in cm und die Massenbelegungsdichten in g/cm² angegeben. Selbstverständlich können ebenso weitere oder andere Rekonstruktionswerte bzw. -paare aus den Messwerten bestimmt werden.
• Im Rahmen des erfindungsgemäßen Verfahrens kann jeweils einem mit den unterschiedlichen Strahlungsspektren aufgenommenen Paar von geglätteten Messwerten ein theoretisches Wertepaar zugeordnet werden. Dies bedeutet, dass jedem Paar von echt gemessenen Messwerten (p₁, p₂) ein theoretisches Wertepaar (M₁(b), M₂(b)) zugeordnet wird, mit dem das Messwertepaar zu identifizieren ist. Dabei bezeichnet der Vektor b = (b₁, b₂) ein Paar von Rekonstruktionswerten b_l und b₂.
• Für die theoretischen Werte M₁(b) sowie M₂(b), die die logarithmische Primärschwächung beschreiben, gilt

  

  wobei W₁(E) bzw. W₂(E) die Aufnahmespektren bezeichnen, und α₁(E) = (μ₁/ρ₁)(E) sowie α₂(E) = (μ₂/ρ₂)(E) die Massenschwächungskoeffizienten angeben, wobei mit μ_i die Schwächungen und mit ρ_i die Dichten angegeben werden. Es wird jeweils in den Grenzen von 0 bis zur jeweiligen Röhrenspannung, multipliziert mit der Elementarladung, integriert.
• Die effektiven Spektren W₁(E) und W₂(E), die die Emissionsspektren der Röntgenröhre bei den beiden verschiedenen Spannungen sowie Spektralfilter und eine energieabhängige Detektoransprechempfindlichkeit beinhalten, sind bekannt und so normiert, dass sich das Integral zu 1 ergibt. Grundsätzlich kann die Zuordnung zwischen einem Messwertepaar und einem Rekonstruktionswertepaar vorausberechnet und in einer Tabellenform abgespeichert werden. Bei einer Diskretisierung von 1000 Werten für die logarithmierten Messwerte umfasst die Tabelle dann 2·10⁶ Einträge.
• Für das erfindungsgemäße Verfahren ist jedoch eine andere Vorgehensweise basierend auf einer Linearisierung angebracht.
• Hierzu wird erfindungsgemäß jedem Rekonstruktionswertepaar ein Paar von diskreten Energien zugeordnet.
• Dies bedeutet, dass zunächst (äquivalente) diskrete Energien E₁ = E₁(b₁, b₂) sowie E₂ = E₂(b₁, b₂) angegeben werden.
• Zu jedem theoretischen Wertepaar kann dann eine Koeffizientenmatrix mit von den diskreten Energien abhängigen Koeffizienten bestimmt werden, die den Übergang von einem Rekonstruktionswertepaar zum zugeordneten theoretischen Wertepaar angibt. Es lassen sich also lineare Beziehungen M₁(b₁, b₂) = α₁(E₁)b₁ + α₂(E₁)b₂ sowie M₂(b₁, b₂) = α₁(E₂)b₁ + α₂(E₂)b₂ angeben, die den Zusammenhang zwischen den einzelnen Komponenten M₁ sowie M₂ der theoretischen Wertepaare und den entsprechenden Komponenten b₁ und b₂ der Rekonstruktionswertepaare angeben. Die Koeffizienten dieser linearen Beziehungen, also die Koeffizienten α_i(E_k), bilden eine Koeffizientenmatrix, deren Einträge jeweils Massenschwächungskoeffizienten α_i darstellen. Diese Matrix wird im Folgenden als Matrix A bezeichnet. Um die Abhängigkeit der Matrix A von den Rekonstruktionswertepaaren darzustellen, kann diese mit einem entsprechenden Index (b) versehen werden,

  

  wobei für die Messwertepaare

  

  gilt.
• Daraufhin wird erfindungsgemäß die Koeffizientenmatrix zur Bildung eines Inversionsoperators und zur anschließenden Bestimmung eines einem tatsächlichen Wertepaar zugeordneten Rekonstruktionswertepaars in Abhängigkeit von einem zugehörigen geglätteten Messwertepaar invertiert. Die Inversion, d. h. die Rückrechnung von einem Paar von Messwerten p = (p₁, p₂) auf das zugehörige Paar von Rekonstruktionswerten b = (b₁, b₂), beispielsweise von Materialdicken, kann durch Inversion der Matrix A erfolgen.
• Es ist also erfindungsgemäß im Unterschied zum bisherigen Vorgehen nicht unbedingt erforderlich, eine Tabelle anzugeben, in der jedem Messwertepaar ein Paar von Rekonstruktionswerten zugeordnet wird, sondern die Tabelle kann auch so aufgebaut werden, dass jedem Messwertepaar zunächst die zugehörige Inversionsmatrix bzw. deren vier Koeffizienten zugeordnet werden. Bei der Inversion der Matrix geht ihre Abhängigkeit vom Rekonstruktionswertepaar in eine Abhängigkeit vom Messwertepaar p über, so dass

  

  gilt, wobei

  

  die vom Messwertepaar p abhängige invertierte Matrix bezeichnet. Dabei ist in diese invertierte Matrix, z. B. wie im Folgenden geschildert, eine Rauschfilterung integriert.
• Wird der nichtlineare Zusammenhang, der vorstehend beschrieben wurde, z. B. bei einer Verwendung einer konstanten Inversionsmatrix vernachlässigt, so können sich beispielsweise bei einer Materialdickenbestimmung beachtliche Fehler in der Größenordnung von ein bis zwei Zentimetern ergeben. Dem wird erfindungsgemäß entgegengetreten, indem wie beschrieben eine nichtlineare Inversion mit einer messwerteabhängigen Matrix durchgeführt wird.
• Da davon auszugehen ist, dass die Messwerte, also die Rohdaten, durch Rauschen gestört sind, tritt zu dem exakten Messwertepaar noch ein Fehler, nämlich das Rauschen, additiv hinzu. Das Rauschen der Rohdaten würde ohne eine in die Matrix integrierte Rauschfilterung auf die mit Hilfe der Inversionsmatrix berechneten Rekonstruktionswerte übertragen und bei einer schlechten Kondition der Matrix sogar verstärkt. Dies wäre selbst dann der Fall, wenn die Matrix messwertunabhängig wäre, erst recht aber bei einer messwertabhängigen Inversionsmatrix bzw. einem messwertabhängigen Inversionsoperator. Das Rauschen wird in einem solchen Fall noch durch die Inversionsmatrix selbst mit dem Resultat einer Rauschverstärkung verändert. In Simulationen lässt sich nachweisen, dass beispielsweise bei einer Materialkombination von zwanzig Zentimetern Wasser und einigen Zentimetern Knochen das Pixelrauschen um dreißig Prozent zunimmt, wenn die Inversion messwertabhängig durchgeführt wird, also in direkter Abhängigkeit vom aufgenommenen Rohwertepaar.
• Anders ist die Situation im erfindungsgemäßen Fall, wenn also der Operator in Abhängigkeit von einem zugehörigen geglätteten Messwertepaar bestimmt wird, nicht für das tatsächliche Messwertepaar. Selbstverständlich kann die Rauschfilterung nicht erfindungsgemäß auch auf andere Art und Weise in den Inversionsoperator integriert werden.
• Im Vordergrund steht jedoch der erfindungsgemäße Gedanke, die Inversionsmatrix bzw. den Inversionsoperator auf Grund der Messwertfehler nicht mehr zu ändern, also diesen für bereits geglättete Rohbilder zu bestimmen. Anschließend kann mit diesem geglätteten Operator unter Anwendung auf die ungeglätteten Rohbilddaten mit einem deutlich geringeren Fehler ein Rekonstruktionswertepaar bestimmt werden.
• Diese Vorgehensweise lässt sich mit der Formel

  

  beschreiben, bei der der Index (Sp) der inversen Matrix A^–1 angibt, dass die Matrix zu adaptiv vorgeglätteten Rohdaten gehört. Der Operator der adaptiven Glättung wird durch das Symbol S bezeichnet. Das sich damit ergebende modifizierte Rekonstruktionswertepaar wird mit b ^~ bezeichnet.
• Erfindungsgemäß wird also vorgeschlagen, zu jedem Messwertepaar die vier Koeffizienten einer 2×2-Matrix A^–1 zu tabellieren, anstatt jedem Messwertepaar nur das zugehörige Paar von Rekonstruktionswerten zuzuordnen. Bei Mehr-Energie-Aufnahmen sind entsprechend größere Matrizen zu tabellieren. Durch die vorgeschlagene Form der Tabellierung wird eine größere Flexibilität hinsichtlich der Anwendbarkeit von Verfahren zur Rauschfilterung erreicht.
• Durch die vorgeschlagene Operatorglättung lässt sich das Rauschen als Standardabweichung des Pixelrauschens um typischerweise etwa fünfundzwanzig Prozent reduzieren, wodurch das Signal-zu-Rausch-Verhältnis entsprechend verbessert wird.
• Dies bietet den Vorteil, dass die Patientendosis, also die Belastung des Patienten durch die Röntgenstrahlung, deutlich verringert werden kann, typischerweise auf 0,75² = 56% des bisherigen Wertes, ohne dass damit eine Verschlechterung der Bildqualität einherginge.
• Das Verfahren lässt sich selbstverständlich auch auf Projektionsbildserien anwenden, die zur Bildrekonstruktion in der Computertomographie weiterverarbeitet werden.
• Außerdem kann nach der Anwendung des Inversionsoperators wenigstens eine weitere Rauschfilterung durchgeführt werden. Das Filterverfahren gemäß der Erfindung, das der eigentlichen Materialseparation vorgeschaltet ist, indem die Filterung bereits z. B. dadurch in den Operator eingebunden ist, dass dieser für geglättete Daten bestimmt wird, kann also vorteilhafterweise mit verschiedenen Nachverarbeitungs-Rauschfilterverfahren und -algorithmen kombiniert werden, um so eine weitere Verbesserung der (Qualität der rekonstruierten) Daten zu erhalten.
• Weiterhin wird zweckmäßigerweise auf die aufgenommenen Rohbilder zur Eliminierung systematischer Ungenauigkeiten vor der Glättung wenigstens ein Kalibrierungs- und/oder Korrekturverfahren angewandt und/oder es wird Streustrahlung aus den Rohbildern eliminiert. Für optimale Ergebnisse ist dementsprechend vorauszusetzen, dass die Messdaten abgesehen vom Rauschen keine systematischen Ungenauigkeiten mehr enthalten und dass auch die Streustrahlung eliminiert worden ist.
• Das Verfahren kann seitens einer Recheneinrichtung und/oder automatisch und/oder bedienergestützt durchgeführt werden, insbesondere unter Verwendung wenigstens eines Programmmittels.
• Es kann also z. B. eine Röntgeneinrichtung bzw. eine Einrichtung zur Erstellung strahlenbasierter Aufnahmen vorgesehen sein, die über eine Steuerungseinrichtung in Form einer Recheneinrichtung verfügt, die wiederum Zugriff auf wenigstens ein Programmmittel hat bzw. auf der Programmmittel abgelegt sind, mit deren Hilfe, gegebenenfalls automatisch nach der Aufnahme der Rohbilder, die erfindungsgemäße Bestimmung von Rekonstruktionswerten auf Basis der Operatorglättung durchgeführt wird. Die Rekonstruktionswertbestimmung ist vollautomatisch möglich, also z. B. derart, dass ein Bediener lediglich die Aufnahme der Messdaten startet, woraufhin das weitere Verfahren der Bildrekonstruktion vollkommen autark abläuft. Es ist aber ebenso denkbar, dass die erfindungsgemäße Rauschfilterung nach der Aufnahme der Rohbilddaten durch einen Bediener separat initiiert wird bzw. unter Berücksichtigung weiterer Bedienereingaben, beispielsweise hinsichtlich eines bestimmten auszuwählenden Glättungsverfahrens oder einer Auswahl geeigneter Nachverarbeitungs-Rauschfilteralgorithmen, durchgeführt wird.
• Beim erfindungsgemäßen Verfahren mit Operatorglättung sind zu jedem Messwertepaar je vier reelle Zahlen, nämlich die vier Koeffizienten der Inversionsmatrix, zu tabellieren. Dies bedeutet, dass die Tabelle bei einer Diskretisierung der Messwerte in jeweils tausend Werte 4·10⁶ Einträge umfassen würde. Des Weiteren ist es erforderlich, dass, zumindest temporär, seitens der Recheneinrichtung neben den Rohbilddaten auch auf die geglätteten Rohbilddaten zugegriffen werden kann.
• Fällt ein konkretes Messwertepaar nicht in das Diskretisierungsraster bzw. nicht exakt in dieses Raster, so können die vier Inversionsmatrixkoeffizienten durch Interpolationsverfahren z. B. bilinear aus der Tabelle gewonnen werden.
• Dies bedeutet z. B. im Fall bilinearer Interpolation für jedes Pixel einen Aufwand von etwa zwanzig Additionen und Multiplikationen, nämlich von 4·4 Operationen für die bilineare Interpolation und zusätzlich vier Operationen für die Matrix-Vektor-Multiplikation. Die zusätzliche Bestimmung der bilinearen Interpolationsgewichte umfasst etwa sechs Multiplikationen. Damit ergibt sich ein insgesamt durchaus vertretbarer Speicher- bzw. Rechenaufwand.
• Weitere Vorteile, Merkmale und Einzelheiten der Erfindung ergeben sich anhand der folgenden Ausführungsbeispiele sowie aus den Zeichnungen. Dabei zeigen:
• 1 energieabhängige Massenschwächungskoeffizienten für verschiedene Materialien,
• 2 typische Dual-Energie-Spektren für die Dual-Energie-Projektionsbildgebung,
• 3 und 4 Darstellungen der quantitativen Ungenauigkeit bei Verwendung einer konstanten Inversionsmatrix im Vergleich zur exakten nichtlinearen Inversion und
• 5 und 6 simulierte Ergebnisse eines erfindungsgemäßen Verfahrens zur Reduktion des Bildrauschens.
• In der 1 sind energieabhängige Massenschwächungskoeffizienten (μ/ρ)(E) in cm²/g für verschiedene Materialien aufgetragen, wobei μ die Schwächung, ρ die Dichte und E die Energie bezeichnet. Dabei sind auf der x-Achse 1 die Werte der Photonenenergie in keV (Kiloelektronenvolt) angegeben, während auf der y-Achse 2 in logarithmischer Form die zugehörigen Massenschwächungskoeffizienten in cm² / g aufgetragen sind. Die Kurve 3 gibt den energieabhängigen Verlauf des Massenschwächungskoeffizienten für Jod an, die Kurve 4 denjenigen für Kalzium, die Kurve 5 den Verlauf für Knochengewebe und die Kurven 6 sowie 7 die energieabhängigen Massenschwächungskoeffizienten für Wasser und Fettgewebe.
• Dabei weist Jod gemäß der Kurve 3 im gesamten Energiebereich den höchsten Massenschwächungskoeffizienten auf, während dieser für Fettgewebe gemäß der Kurve 7 am niedrigsten liegt. Am ehesten vergleichbar mit dem Verlauf der Kurve 7 für Fettgewebe ist der Verlauf des Massenschwächungskoeffizienten in Abhängigkeit von der Energie für Wasser gemäß der Kurve 6.
• Knochengewebe sowie Kalzium weisen im Vergleich zu den Kurven 6 und 7 bei niedrigeren Energien wesentlich höhere Massenschwächungskoeffizienten auf. Bei höheren Energien nähert sich die Kurve 5 mehr und mehr an die Kurven 6 bzw. 7 für Wasser und Fettgewebe an. Kalzium weist demgegenüber gemäß der Kurve 4 auch im Bereich von 140 keV noch einen erkennbar größeren Massenschwächungskoeffizienten auf.
• Die 2 zeigt vier typische Dual-Energie-Spektren gemäß den Kurven 8, 9, 10 und 11. Bei diesen Spektren handelt es sich um normierte Spektren für verschiedene Spannungen und Vorfilterungen. Dabei ist auf der x-Achse 12 jeweils die Quantenenergie in keV aufgetragen, während auf der y-Achse 13 die normierte Intensität aufgetragen ist.
• Die Kurve 8 gibt die energieabhängige Intensität für eine Spannung von 60 kV bei einer Vorfilterung durch 0,1 mm Kupfer wieder, während bei der Kurve 9 bei gleicher Spannung eine Vorfilterung durch 0,3 mm Kupfer zu Grunde liegt.
• Die Kurven 10 und 11 sind jeweils Spannungen von 150 kV zugeordnet, wobei die Kurve 10 das Spektrum bei einer Vorfilterung mit 0,1 mm Kupfer angibt, die Kurve 11 dasjenige bei einer Vorfilterung mit 0,3 mm Kupfer.
• Während die Kurven 8 und 9 für die niedrigeren Spannungen entsprechend bei niedrigeren Quantenenergien einen deutlichen Peak zeigen und der jeweilige Intensitätswert bei 60 keV bereits auf 0 abgefallen ist, sind die Spektren bei den höheren Spannungen breiter und weisen stattdessen zwei kleinere Peaks auf. Die unterschiedlichen Vorfilterungen führen zu weiteren Abweichungen bei den Kurvenformen, z. B. bei der niedrigen Spannung zu einer Verschiebung des Peaks in den Bereich höherer Energien. Da sich die jeweiligen Spektren bei niedriger bzw. höherer Spannung möglichst wenig überlappen sollten, ist es zweckmäßig z. B. bei 60 kV mit nur 0,1 mm Kupfer, bei 150 kV aber mit 0,3 mm Kupfer vorzufiltern.
• In den 3 und 4 sind Darstellungen der quantitativen Ungenauigkeit bei Verwendung einer konstanten Inversionsmatrix im Unterschied zu einer messwertabhängigen Matrix dargestellt.
• Dabei ist in der 3 auf der x-Achse 14 eine tatsächliche Knochendicke in Gramm pro cm² (g/cm²) aufgetragen, während auf der y-Achse 15 die rekonstruierte Dicke, ebenfalls in Gramm pro cm², dargestellt ist. Die exakte, nichtlineare Inversion entspricht für Knochengewebe der Kurve 16. Die Kurve 17 gibt den deutlich abweichenden Verlauf bei Verwendung einer konstanten Matrix, also einer Näherung statt der exakten Matrix, bei 20 cm Wasser und 0 cm Knochengewebe wieder, die Kurve 18 die Ergebnisse für eine konstante Matrix bei 20 cm Wasser und 5 cm Knochengewebe. Die exakten Verhältnisse, die einer nichtlinearen Inversion entsprechen, sind für Wasser durch die Kurve 19, die Näherungswerte für 20 cm Wasser bei 5 cm Knochengewebe durch die Kurve 20 und für 20 cm Wasser bei 0 cm Knochengewebe durch die Kurve 21 gegeben.
• Die zugehörigen Rekonstruktionsfehler in Gramm pro cm² (g/cm²) sind auf der y-Achse 22 der 4 gegen die wahre Knochendicke in Gramm pro cm² auf der x-Achse 23 aufgetragen.
• Die Kurven 24 für Knochen bzw. 25 für Wasser geben jeweils die Fehler bei 20 cm Wasser und 5 cm Knochen wieder, während die Kurven 26 und 27 für Knochen bzw. Wasser den Rekonstruktionsfehler bei 20 cm Wasser und 0 cm Knochengewebe zeigen.
• Diesen Kurven 24–27 ist zu entnehmen, dass die Rekonstruktionsfehler durchaus im Bereich mehrerer Gramm pro cm² liegen können. Selbst im Fall von 20 cm Wasser und 5 cm Knochengewebe gemäß den Kurven 24 und 25 sind bei hohen tatsächlichen Knochendichten Abweichungen von fast 2 Gramm pro cm² möglich. Auch bei sehr geringen tatsächlichen Knochendichten liegen die Abweichungen über einem Gramm pro cm². Bei 20 cm Wasser und 0 cm Knochen steigt der Fehler auf etwa +/–5 g/cm² bei einer tatsächlichen Knochendicke von 10 g/cm² an.
• Dementsprechend ist die Vernachlässigung der Messwertabhängigkeit der Matrix für eine aussagekräftige quantitative Bildgebung nicht akzeptabel.
• Die 5 und 6 zeigen schließlich simulierte Ergebnisse eines erfindungsgemäßen Verfahrens zur Reduktion des Bildrauschens. Die Darstellung bezieht sich auf ein Simulationsobjekt in Form einer homogenen Schicht aus Wasser mit einer Dicke von 20 cm, über der kreisförmige Scheiben aus Knochen mit Durchmessern von je 20 Pixeln angeordnet sind, deren Belegungsdicken zwischen 50 und 250 mg/cm² in Stufen von 50 mg/cm² variieren. Als Strahlungsquelle wird eine Röntgenröhre mit einer Wolframanode und einer internen Filterung, die 2,5 mm Aluminium entspricht, verwendet. Der Detektor ist ein Flat-Panel-Detektor mit Caesiumiodid mit einer Dichte von 100 mg/cm² als Szintillatormaterial. Als Strahlungsspektren werden Spektren mit einer Spannung von 70 kV mit einer Vorfilterung von 0,1 mm Kupfer bzw. mit einer Spannung von 150 kV mit einer Vorfilterung von 0,3 mm Kupfer verwendet. Für das Rauschen des Messsignals am Detektor wird eine rauschäquivalente Quantenzahl von 1000 zu Grunde gelegt, die relative Standardabweichung des Rauschens bei jedem Detektorpixel entspricht

  
• Die erste Zeile 28 zeigt das ungefilterte Knochenbild gemäß der Gleichung (1)

  
• In der zweiten Zeile 29 ist das Ergebnis der Rekonstruktion der erfindungsgemäßen Operatorglättungsmethode nach der Gleichung (2)

  

  gezeigt. Damit wird eine verbesserte Objekterkennung ermöglicht. Dementsprechend weist das Ergebnis weniger Fehler auf.
• Die Zeile 30 schließlich zeigt das Rekonstruktionsergebnis bei der Kombination der Operatorglättungsmethode mit einer nachfolgenden Bildglättung gemäß der Formel (3)

  

  In der 6 sind auf der waagerechten Achse 31 jeweils Pixelzahlen aufgetragen, während auf der senkrechten Achse 32 der Fehler in Gramm pro cm² (g/cm²) bei einer Rauschfilterung gemäß den Matrixformeln (1), (2) und (3) aufgetragen ist.
• Dem ist zu entnehmen, dass bei der exakten Inversion ohne Filterung gemäß der Spalte 33 (zur Matrixformel (1)) noch größere Fehler auftreten, die bei der Operatorglättungsmethode der Erfindung gemäß der Spalte 34 (Formel (2)) bereits deutlich verringert sind. Eine weitere Verringerung der Fehler wird durch eine Kombination mit einer nachträglichen Bildfilterung mit einer 5×5-Maske gemäß der Spalte 35 (Formel (3)) erreicht. Bei Kombination der Operatorglättungsmethode mit einer nachträglichen Bildglättung kann der Fehler durchgängig unter 0,2 g/cm² gehalten werden. Auch bei einfacher Verwendung der Operatorglättungsmethode gemäß der Spalte 34 ohne nachträgliche Glättung treten kaum Fehler über 0,4 g/cm² auf.
• Somit ermöglicht die erfindungsgemäße Rauschfilterung, die in einen Inversionsoperator integriert ist, eine deutliche Verbesserung der quantitativen Bildgebung im Rahmen der Zweispektren-Projektionsbildgebung.

Claims (8)

1. Verfahren zur Reduktion des Bildrauschens im Rahmen der Aufnahme wenigstens eines strahlungsbasierten Bildes eines Bildaufnahmebereichs mit zwei unterschiedlichen Röntgenspektren mit den Schritten – Aufnahme von Rohbildern des Bildaufnahmebereichs mit den beiden unterschiedlichen Röntgenspektren mit jeweils paarweise einander zugeordneten Messwerten und – zur Separation unterschiedlicher Materialien im Bildaufnahmebereich Anwendung wenigstens eines einen Übergang von einem Messwertepaar in ein zugeordnetes Rekonstruktionswertepaar beschreibenden Inversionsoperators mit einer integrierten Rauschfilterung auf die aufgenommenen Rohbilder, wobei die Rauschfilterung in den Inversionsoperator integriert wird, indem dieser für jeweils mit den unterschiedlichen Röntgenspektren aufgenommene geglättete Messwertepaare der Rohbilder bestimmt wird, wobei als Inversionsoperator wenigstens eine Inversionsmatrix bestimmt wird zu jedem geglätteten Messwertepaar der Rohbilder, wobei eine Koeffizientenmatrix zur Bildung der Inversionsmatrix in Abhängigkeit von einem zugehörigen geglätteten Messwertepaar invertiert wird.
2. Verfahren nach Anspruch 1, wobei die für die jeweils mit den beiden unterschiedlichen Röntgenspektren aufgenommenen und geglätteten Rohbilder bestimmte Inversionsmatrix zur Bestimmung eines Rekonstruktionswertepaars auf ein entsprechendes tatsächliches Messwertepaar der ungeglätteten Rohbilder angewandt wird.
3. Verfahren nach einem der vorangehenden Ansprüche, wobei die Inversionsmatrix eine 2×2-Matrix ist.
4. Verfahren nach einem der vorangehenden Ansprüche, wobei jeweils einem mit den beiden unterschiedlichen Röntgenspektren aufgenommenen Paar von Messwerten als wenigstens ein Rekonstruktionswertepaar ein Paar von Massenbelegungsdichten und/oder Materialdicken zugeordnet wird.
5. Verfahren nach einem der vorangehenden Ansprüche, wobei zu jedem theoretischen Wertepaar eine Koeffizientenmatrix mit von diskreten Energien abhängigen Koeffizienten bestimmt wird, die den Übergang von einem Rekonstruktionswertepaar zum zugeordneten theoretischen Wertepaar angibt.
6. Verfahren nach einem der vorangehenden Ansprüche, wobei nach der Anwendung der Inversionsmatrix wenigstens eine weitere Rauschfilterung durchgeführt wird.
7. Verfahren nach einem der vorangehenden Ansprüche, wobei auf die aufgenommen Rohbilder zur Eliminierung systematischer Ungenauigkeiten vor der Glättung wenigstens ein Kalibrierungs- und/oder Korrekturverfahren angewandt wird und/oder Streustrahlung aus den Rohbildern eliminiert wird.
8. Verfahren nach einem der vorangehenden Ansprüche, wobei das Verfahren seitens einer Recheneinrichtung und/oder automatisch und/oder bedienergestützt durchgeführt wird, insbesondere unter Verwendung wenigstens eines Programmmittels.

Images