DE102010011911B4

DE102010011911B4 - Tomosyntheseverfahren mit einer iterativen Maximum-A-Posteriori-Rekonstruktion

Info

Publication number: DE102010011911B4
Application number: DE102010011911.3A
Authority: DE
Inventors: Dr. Jerebko Anna
Original assignee: Siemens Healthcare GmbH
Current assignee: Siemens Healthcare GmbH
Priority date: 2010-03-18
Filing date: 2010-03-18
Publication date: 2019-09-26
Anticipated expiration: 2030-03-19
Also published as: US20110229007A1; DE102010011911A1; US8503755B2

Abstract

Verfahren zum Ermitteln eines Schätzwertes eines Schwächungskoeffizienten µeines Voxels j, das einem Abschnitt eines Gewebes (10) zugeordnet ist, wobei der Schwächungskoeffizient µ die Schwächungsrate einer Strahlung (4) in dem Gewebe (10) angibt, wobei die Strahlung (4) von einem Emitter (2) abgegeben wird und von einem Detektor (18) erfasst wird und wobei das Verfahren folgende Schritte aufweist:- Schätzen des Rauschschätzwertes Δµ̂ des Schwächungskoeffizienten µ eines Voxels (S1);- Schätzen eines mittleren Schwächungskoeffizienten µ̅ (S2) ;- Berechnen der Paramater einer Prior-Funktion U(p) aus dem Rauschschätzwert Δµ̃ (S9); und- Iteratives Berechnen des Schätzwertes eines Schwächungskoeffizienten µeines Voxels j aus mittels des Detektors (18) erfassten Projektionsdaten, aus dem Rauschschätzwert Δµ̃, aus dem mittleren Schwächungskoeffizienten µ̅ und aus der Prior-Funktion U(p) (S10, S11) mittels folgender Formel (S10):wobei β ein Optimierungsparameter ist, di die erwartete Photonenzahl ist, die aus der Strahlenquelle entlang der Projektion i austritt,die gemessenen Photonenanzahl in der Projektion i ist und ldie Schnittlänge des Strahls durch das Voxel j in der Projektion i ist undwobei β mittels folgender Formel bestimmt wird (S9):wobeider Mittelwert über alle Voxel des Volumenelementes über alle i und j ist.

Description

Die vorliegende Erfindung betrifft ein iteratives Tomosyntheseverfahren, das eine iterative Maximum-A-Posteriori-Rekonstruktion verwendet, und insbesondere die analytische Bestimmung der Parameter einer bei der Maximum-A-Posteriori-Rekonstruktion verwendeten Prior-Funktion.
Bei einem Tomosyntheseverfahren wird ein dreidimensionales Bild aus einer Mehrzahl zweidimensionaler Bilder erzeugt. Mittels einer Röntgeneinrichtung mit einer Röntgenstrahlquelle und einem Detektor wird ein zweidimensionales Bild bzw. eine Projektion des zu untersuchenden Gewebes, z.B. eines zu untersuchenden Organs, erzeugt, das der Röntgenstrahl durchläuft. Das zweidimensionale Bild stellt hierbei die Dämpfung des Gewebes in dem Volumen dar, das der Strahl durchlaufen hat. Ein zweites zweidimensionales Bild bzw. eine zweite Projektion des gleichen Gewebes bzw. Volumens wird aufgenommen, nachdem die Strahlquelle und/oder der Detektor in eine zweite Stellung bewegt wurde. Nachdem eine Mehrzahl zweidimensionaler Bilder aufgenommen wurde, kann ein dreidimensionales Tomosynthesebild mittels einer Rekonstruktion erzeugt werden. Die vorliegende Erfindung kann bei der dreidimensionalen Tomosynthese bei einer Röntgentechnik, Computertomographie oder bei der Magnetresonanztomographie angewendet werden.
Ein Anwendungsgebiet der eingangs erwähnten dreidimensionalen Bildgebungsverfahren ist die Mammographie. Eine typischerweise in der Mammographie verwendete Bilderzeugungsvorrichtung umfasst eine schwenkbare Röntgenstrahlquelle und einen stationären Röntgendetektor. Das zu untersuchende Gewebe wird über dem stationären Detektor positioniert. Anschließend wird die Röntgenquelle in mehreren Schritten geschwenkt, beispielsweise in einem Bereich von +/- 25°, und es werden eine Mehrzahl zweidimensionaler Röntgenbilder aus unterschiedlichen Schwenkstellungen der Röntgenstrahlquelle mit dem ortsfesten Detektor aufgenommen. Selbstverständlich ist es auch möglich, eine Mehrzahl ortsfester Röntgenstrahlquellen zu verwenden oder die Röntgenstrahlungsquelle lediglich translatorisch zu verschieben. Auch der Detektor kann entgegen der Bewegung der Röntgenquelle verschoben oder geschwenkt werden. Aus der Mehrzahl zweidimensionaler Röntgenbilder wird mittels der Rekonstruktion ein dreidimensionales Bild erzeugt. Bildgebende Verfahren und Vorrichtungen für die Mammographie des Standes der Technik sind beispielsweise in der DE 10 2006 046 741 A1 , DE 10 2008 004 473 A1 und der DE 10 2008 028 387 A1 beschrieben.
Im Stand der Technik werden zur Rekonstruktion eines dreidimensionalen Bildes aus einer Mehrzahl zweidimensionaler Bilder so genannte gefilterte Rückprojektionen verwendet, die beispielsweise in Imaging Systems for Medical Diagnostics, Arnulf Oppelt, Publicis Corporate Publishing, Erlangen, ISBN 3-89578-226-2, Kapitel 10.5 beschrieben sind. Diese gefilterten Rückprojektionsrekonstruktionsverfahren stellen rekonstruierte Bilder mit einem hohen Kontrast und einer hohen Detailtreue dar, aber verlieren bei der Tomosynthese mit eingeschränktem Abtastwinkel aufgrund der fehlenden Daten Information über die relative Gewebedichte. Dies wird dadurch verursacht, dass bestimmte Filterkerne niederfrequente Anteile entfernen. Im Allgemeinen wird die digitale Brusttomosynthese (DBT - Digital Breast Tomosynthesis) durch unvollständige Daten und eine schlechte Quantenstatistik beeinträchtigt, die durch die Gesamtdosis beschränkt ist, die in der Brust absorbiert wird. Die Brust besteht hauptsächlich aus Drüsengewebe, Fettgewebe, Bindegewebe und Blutgefäßen. Die Röntgenschwächungskoeffizienten dieser Gewebetypen ähneln sich stark, was die Auswertung dreidimensionaler Mammographiebilder erheblich erschwert. Der Hauptanwendungsbereich von bildgebenden Verfahren in der Mammographie ist die frühzeitige Erkennung krebshaltigen Gewebes. Erschwerend kommt hierbei hinzu, dass krebshaltiges Gewebe einen ähnlichen Röntgenschwächungskoeffizienten wie andere Gewebetypen aufweist. Mammographieverfahren sind beispielsweise in Imaging Systems for Medical Diagnostics, Arnulf Oppelt, Kapitel 12.6, Publicis Corporate Publishing, Erlangen, ISBN 3-89578-226-2 beschrieben.
Ein weiteres Problem der zuvor erwähnten gefilterten Rückprojektionsrekonstruktionsverfahren sind Artefakte außerhalb der Ebene (sogenannte Out-of-Plane-Artefakte), die durch eine Filterung zusammen mit den Bildmerkmalen verstärkt werden. Es wurden statistische iterative Rekonstruktionsverfahren vorgeschlagen, die die Ähnlichkeiten zwischen den berechneten und gemessenen Projektionen maximieren und eine Rauschunterdrückung, beispielsweise durch Prior-Funktionen, ermöglichen.
Es sind so genannte Maximum-Likelihood-Verfahren bekannt, mit denen der Schätzwert µ des Schwächungskoeffizienten, beispielsweise des Brustvolumenschwächungskoeffizienten, bestimmt werden kann, der die logarithmische Wahrscheinlichkeits-Funktion L(µ) maximiert: $μ_{max} =_{μ}^{argmax} {L (μ)};$
Die Maximum-Likelihood-Rekonstruktion führt zu vergleichsweise guten Ergebnissen in der DBT und konvergiert in 4 bis 5 Iterationen. Ein Nachteil derartiger Rekonstruktionen ist, dass ohne die Verwendung von Prior- oder Straf-Funktionen übermäßig verrauschte Bilder entstehen. Es wurde die Verwendung von Glättungs-Prior-Funktionen vorgeschlagen. Die logarithmische Wahrscheinlichkeits-Funktion L(µ) wird hierbei in eine Log-Posterior-Funktion Φ(µ) geändert: $Φ (μ) = L (μ) - β U (μ);$
wobei U(µ) eine Glättungs-Prior-Funktion ist, die die Differenzen der Werte benachbarter Pixel reduziert bzw. bestraft. Der Parameter β > 0 ist ein Steuerungsparameter.
Ein Nachteil dieser Verfahren ist jedoch, dass die optimalen Parameter der Prior-Funktionen empirisch ermittelt werden müssen, was üblicherweise eine Suche über einen großen Wertebereich umfasst, was bei hochauflösenden Bildern unerwünscht ist, insbesondere wenn mehrere Prior-Funktionen beurteilt werden müssen.
Die US 7 526 060 B2 und die US 7 308 071 B2 offenbaren die Verwendung statistischer Methoden bei der Auswertung von CT-Bildern. Die Veröffentlichungen Chan, C. [et al.], An anatomically based regionally adaptive prior for MAP reconstruction in emission tomography, IEEE Nuclear Science Symposium Reference Record (2007) 4137 -4141, die Veröffentlichung Das, M. [et al.], Dose Reduction in breast tomosynthesis using a penalized maximum likelihood reconstruction, Proceedings of SPIE, Band 7258.725855, die Veröffentlichung De Man B., Basu S., An efficient maximum-likelihood reconstruction algorithm for transmission tomography and a generalized Geman prior, Abstracts of RSNA 2005, die Veröffentlichung Lalush D. S., Tsui, B.M.W., Simulation evaluation of Gibbs prior distributions for use in maximum a posteriori SPECT reconstruction, IEEE Transactions on Medical Imaging, Band 11, Nr. 2 (1992) 267-275 und die Veröffentlichung Lange K., Fessler, J.A., Globally convergent algorithms for maximum a posteriori transmission tomography, IEEE Transaction on Image Proceedings, Band 4, Nr. 10 (1995) 1430-1438 offenbaren statistische Verfahren zum Auswerten von CT-Bildern.
Die Erfindung stellt sich daher die Aufgabe, ein verbessertes Tomosynthese-Rekonstruktionsverfahren anzugeben, das ein dreidimensionales Bild mit einem reduzierten Rauschen bereitstellt, ohne den Kontrast der Kanten des Gewebes zu beeinträchtigen.
Die Aufgabe wird erfindungsgemäß durch ein Verfahren zum Ermitteln eines Schätzwertes eines Schwächungskoeffizienten µ_j eines Voxels j gelöst, das einem Abschnitt eines Gewebes zugeordnet ist. Der Schwächungskoeffizient µ gibt die Schwächung pro Länge einer Strahlung in dem Gewebe an, wobei die Strahlung von einem Emitter, z.B. einer Röntgenstrahlröhre, abgegeben wird und von einem Detektor erfasst wird. Der Schwächungskoeffizient repräsentiert die Gewebedichte. Der Rauschschätzwert Δµ̂ des Schwächungskoeffizienten µ eines Voxels wird geschätzt. Vorzugsweise wird der Rauschschätzwert Δµ̂ des Schwächungskoeffizienten µ eines Voxels aus der Gewebedicke und Informationen über die abgegebene Strahlung geschätzt. Der Rauschschätzwert Δµ̂ muss kleiner als der kleinste zu erhaltende Gradient sein. Der Rauschschätzwert Δµ̂ ist vorzugsweise der maximale Rauschwert, der in einem Bild erwartet wird, wenn keine Glättung erfolgt.
Der Rauschschätzwert Δµ̂ kann für ein gegebenes Röntgenspektrum beispielsweise als Funktion der Gewebedicke, beispielsweise der Brustdicke, und der Röhrenladung (mAs) unter Verwendung von Regressionsverfahren bestimmt werden. Im Allgemeinen kann der Rauschpegel für eine beliebige Brustdicke durch die Auswahl der Röhrenladung und der Photonenenergie (keV) konstant gehalten werden.
Der mittlere Gewebeschwächungskoeffizient µ̅ wird geschätzt. Vorzugsweise wird der mittlere Gewebeschwächungskoeffizient µ̅ auf Grundlage der Energie der vom Emitter abgegebenen Strahlung und auf Grundlage einer Gewebedichteschätzung geschätzt. Beispielsweise kann der mittlere Gewebedämpfungsschätzwert µ̅ auf Grundlage der Photonenenergie und einer Brustdickeschätzung geschätzt werden, die durch segmentierte Originalprojektionen unter Verwendung herkömmlicher Verfahren ermittelt werden.
Der Rauschschätzwert Δµ̂ und der mittlere Gewebeschwächungskoeffizient µ̅ sowie ihre Abhängigkeit von den zuvor genannten physikalischen Größen können mittels statistischer Untersuchungen an einer Mehrzahl von Patienten oder Gewebeproben ermittelt werden. Der Rauschschätzwert Δµ̂ und der mittlere Gewebeschwächunskoeffizient µ̅ können durch eine statistische Studie an fünf Patienten oder Gewebeproben oder einer größeren Anzahl von Patienten bzw. Gewebeproben bestimmt werden. Der Zusammenhang zwischen dem Rauschschätzwert Δµ̂ sowie dem mittlere Gewebeschwächungskoeffizienten µ̅ und den zuvor genannten physikalischen Größen kann durch ein Regressionsverfahren bestimmt werden. Der Rauschschätzwert Δµ̂ ist vorzugsweise der höchste zu kompensierende Rauschwert.
Die Parameter einer Prior-Funktion U(µ) werden aus dem erwarteten Rauschwert Δµ̂ und dem mittleren Gewebedämpfungsschätzwert µ̅ berechnet. Anschließend wird der Schätzwert des Schwächungskoeffizienten µ_j eines Voxels aus mittels des Detektors erfassten Projektionsdaten, dem Rauschschätzwert Δµ̂, dem mittleren Gewebeschwächungskoeffizienten µ̅ und der Prior-Funktion U(µ) iterativ berechnet. Der erwartete Rauschwert Δµ̂ ist derjenige Rauschwert, der in einem Bild erwartet wird, wenn keine Glättung mittels der Prior-Funktion erfolgt.
Der Ausdruck Voxel bezeichnet eine Volumeneinheit in der dreidimensionalen Bildverarbeitung. Er setzt sich aus den Ausdrücken Pixel und Volumen zusammen. Die Prior-Funktion stellt sicher, dass Kanten erhalten bleiben, während das Rauschen unterdrückt wird. Die Schwächung kann aufgrund einer Absorption verursacht werden. Der Schwächungskoeffizient kann ein Absorptionskoeffizient sein. Die Schwächung kann auch eine Spin-Relaxation sein.
Die Prior-Funktion U(µ) kann eine Geman-Prior-Funktion sein. Die Geman-Prior-Funktion kann mittels folgender Gleichung bestimmt werden: $U (μ) = \sum_{j} \sum_{k} \frac{Δ μ_{j k}^{2} σ^{m}}{2 {(\sqrt{Δ μ_{j k}^{2} / 2 + σ^{2} / 2})}^{m}};$
wobei σ ein Optimierungsparameter und m eine Konstante ist und wobei µ_j und µ_k Schwächungskoeffizienten benachbarter Voxel sind, die je einem Gewebeabschnitt zugeordnet sind, und j und k lineare Voxel-Indizes über drei Dimensionen sind. Die Werte der Schwächungskoeffizienten µ_k und µ_j werden mit jeder Iteration aktualisiert. D.h., bei der ersten Iteration werden diese Werte mit Schätzwerten vorbelegt, beispielsweise mit dem mittleren Gewebeschwächungswert. Bei den anschließenden Iterationen werden für die Schwächungskoeffizienten µ_k und µ_j diejenigen Werte verwendet, die bei der vorangegangenen Iteration bestimmt wurden. Es kann beispielsweise eine Nachbarschaft von 3x3x3 Voxeln ausgewertet werden, um die Differenzen µ_k der Gewebeschwächungskoeffizienten zu bestimmen.
Mittels der Parameter σ und m kann die Prior-Funktion für das jeweilige Anwendungsgebiet optimiert werden. Um die Kanten des Bildes möglichst gut zu erhalten und die Konvergenz sicherzustellen, werden σ und m vorzugsweise so ausgewählt, dass die Prior-Funktion eine konvexe Form annimmt.
Bevor das Bild rekonstruiert wurde, sind die Differenzen Δµ_jk nicht bekannt. Wird das Rekonstruktionsverfahren auf Daten angewendet, die mittels Röntgenstrahlen erzeugt wurden, wird angenommen, dass das Rauschen in der Projektion i ortsunabhängig ist. Daher wird zur Bestimmung des optimalen Wertes für den Parameter σ in der Geman-Prior-Funktion der Term Δµ_jk durch den Rauschschätzwert, vorzugsweise den maximalen Rauschschätzwert, Δµ̃ ersetzt. Wurden die Bilddaten mittels eines Magnetresonanztomographen aufgenommen, kann das Rauschen ortsabhängig sein. In diesem Fall kann der Wert Δµ̃_j vom Ort des Voxels j abhängig sein.
Es ist jedoch auch möglich, die Werte µ_j und µ_k und somit Δµ_jk mittels herkömmlicher Rekonstruktionsverfahren zu bestimmen, beispielsweise mittels einer linearen und/oder gefilterten Rückproduktion. Hierdurch kann der Parameter σ genauer bestimmt werden.
Falls die Prior-Funktion eine konvexe Form aufweisen soll, wird σ so gewählt, dass die erste Ableitung der Prior-Funktion ihren stärksten Krümmungspunkt bei Werten von Δµ erreicht, die dem Rauschschätzwert Δµ̂ entsprechen. Vorzugsweise ist m ≤ 16/17.
Falls die Prior-Funktion eine nicht-konvexe Form aufweisen soll, wird σ so gewählt, dass die erste Ableitung der Prior-Funktion ihr Maximum bei Werten von Δµ erreicht, die dem Rauschschätzwert Δµ̂ entsprechen. Vorzugsweise ist m = 2.
Das Verfahren weist ferner den Schritt des iterativen Berechnens des Schätzwertes des Schwächungskoeffizientens µ_j eines Voxels j zur Rekonstruktion eines dreidimensionalen Bildes mittels folgender Formel auf: $μ_{j}^{n + 1} = μ_{j}^{n} + \frac{μ_{j}^{n}}{\sum_{i} Y_{i} l_{i j} + β μ_{j}^{n} \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ^{n})} \frac{\partial}{\partial μ_{j}} Φ (μ^{n});$
wobei β ein Optimierungsparameter ist, und Φ(µⁿ) = L(µ) - βU(µ) eine Log-Posterior-Funktion ist, L(µ) eine Zielfunktion ist, Yi die gemessene Photonenanzahl der Projektion i ist und lij die Schnittlänge des Strahls i durch das Voxel j ist.
Dieses Verfahren wird als Maximum-A-Posterior-Verfahren mit einem Newton-basierten Aktualisierungsverfahren bezeichnet, das eine einfache quadratische Glättungs-Prior-Funktion verwendet. Die Parameter σ und β können bei allen Iterationsschritten konstant sein.
Die Zielfunktion L(µ) wird gemäß einer ersten Ausführungsform der Erfindung wie folgt bestimmt: $L (μ) = \sum_{i} {- d_{i} e^{- 〈 l_{i}, μ 〉} - Y_{i} 〈 l_{i}, μ 〉} + c_{;}$
wobei $〈 l_{i}, μ 〉 = \sum_{j} l_{i j} μ_{j}$
ist, i ein Projektionsindex ist, j der Voxelindex ist, d_i die erwartete Photonenanzahl ist, die aus der Strahlungsquelle entlang der Projektion i austritt, Y_I die erfasste Photonenanzahl in der Projektion i ist und c eine (zu vernachlässigende) Konstante ist.
Der Steuerparameter β der Geman-Prior-Funktion wird gemäß dieser Ausführungsform mittels folgender Formel bestimmt: $β = \frac{Δ \tilde{μ} \cdot \bar{\sum_{i} Y_{i} l_{i j}}}{\bar{μ} \cdot \frac{\partial}{\partial μ_{j}} U (μ) {|_{Δ μ_{j} = Δ \tilde{μ}} - Δ \tilde{μ} \cdot \bar{μ} \cdot \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ) |}_{^{^{^{_{^{^{Δ μ_{_{j}} = Δ \tilde{μ}}}}}}}}};$
wobei $\bar{\sum_{i} Y_{i} l_{i j}}$
der Mittelwert über alle Voxel des Volumenelementes über alle i und j ist. Dieser Mittelwert kann zur Beschleunigung der Iterationsverfahren einmal vor den Iterationen berechnet werden.
Der Steuerparameter β kann vor der ersten Iteration unter Zuhilfenahme der Geman-Prior-Funktion U(µ) bestimmt werden. Zu diesem Zeitpunkt ist jedoch die Differenz benachbarter Schwächungskoeffizienen Δµ_jk noch nicht bekannt. Werden die Messdaten mittels Röntgenstrahlung ermittelt, kann angenommen werden, dass das Rauschen in einer Projektion konstant ist. In diesem Fall wird daher der Ausdruck Δµ_jk durch den maximalen Rauschschätzwert Δµ̃ ersetzt. Wurden die Messdaten mittels eines Magnetresonanztomographen erzeugt, kann das Rauschen ortsabhängig sein, wodurch der Ausdruck Δµ_jk in der Geman-Prior-Funktion U(µ) durch einen ortsabhängigen Rauschschätzwert Δµ̃_j ersetzt werden kann. Wie zuvor erwähnt wurde, kann die Differenz der Schwächungskoeffizienten Δµ_jk zuvor durch ein herkömmliches Rückprojektionsverfahren ermittelt werden. In diesem Fall kann β genauer bestimmt werden.
Durch die Bestimmung der Parameter β und σ kann der Rauschpegel effektiv kontrolliert werden, während die Kanten und Ränder der Gewebebereiche erhalten bleiben. Darüber hinaus werden mit dem vorliegenden Verfahren der Signalrauschabstand (SNR) und der Kontrast/Rausch-Abstand (CNR) erheblich verbessert. Der verbesserte Signalrauschabstand und der verbesserte Kontrast/Rausch-Abstand können bewirken, dass eine niedrigere Strahlendosis erforderlich ist, obwohl dreidimensionale Bilder besserer Qualität erzielt werden, was einerseits die Sicherheit der Behandlung verbessert und andererseits die Nebenwirkungen der Behandlung reduziert. Die Reduzierung des Wertes bzw. die Penale, d. h. die sogenannte Penalty, für Voxel-Intensitätsabweichungen kleiner oder gleich dem höchsten Rauschschätzwert ist proportional zu diesen Abweichungen. Gleichzeitig werden die Kanten mit Voxel-Intensitätsabweichungen, die höher als der höchste Rauschschätzwert Δµ̂ sind, lediglich mit einem kleinen Korrekturwert bzw. einer kleineren Penale versehen, die demjenigen bzw. derjenigen des höchsten Rauschschätzwertes entspricht. Bei nicht-konvexen Optimierungsverfahren ist es möglich, eine Geman-Prior-Funktion mit m = 2 auszuwählen, bei der die Penale bzw. der Korrekturwert für Kantenvoxel nahezu null ist.
Der Startwert µ_j ⁰ kann mit dem gleichen Wert für alle j festgelegt werden. Alle Voxel in dem zu untersuchenden Volumen könnten mit dem durchschnittlichen Schwächungswert µ des zu untersuchenden Gewebes initialisiert werden. Eine schnellere Konvergenz kann erzielt werden, wenn der Startwert µ_j ⁰ mit einem normalisierten Rückprojektionswert vorbelegt wird, wie beispielsweise in der Veröffentlichung Ludwig J., Mertelmeier T., Kunze H., Härer W., A Novel Approach for Filtered Projection in Tomosynthesis Based on Filter Kernels Determined by Iterative Reconstruction Techniques, Krupinski E., Lecture Notes in Computer Science 5116, Digital Mammography, 9th International Workshop, IWDM 2008, pp. 612-620, Springer Verlag Berlin Heidelberg (2008) beschrieben ist. Durchschnittlich sind drei Iterationen erforderlich. Es können jedoch bis zu 20 Iterationen erforderlich sein. Das iterative Bestimmungsverfahren kann abgebrochen werden, wenn ein Endekriterium erfüllt ist. Ein mögliches Endekriterium ist, dass die Differenz aus µ_j ⁿ⁺¹ und µ_j ⁿ kleiner als ein vorbestimmter Schwellenwert wird.
Das Rekonstruktions-Verfahren wird mittels einer iterativen Gradienten-Optimierungs-Aktualisierungs-Formel durchgeführt. Der Schritt des iterativen Berechnens des Schätzwertes des Schwächungskoeffizienten µ_j eines Voxels j wird gemäß der ersten Ausführungsform der Erfindung mittels folgender Formel durchgeführt: $μ_{j}^{n + 1} = μ_{j}^{n} + μ_{j}^{n} \frac{\sum_{i} l_{i j} [d_{i} e^{- 〈 l_{i}, μ 〉} - Y_{i}] - β \frac{\partial}{\partial μ_{j}} U (μ^{n})}{\sum_{i} Y_{i} l_{i j} + β μ_{j}^{n} \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ^{n})};$
wobei β ein Optimierungsparameter ist, $〈 l_{i}, μ 〉 = \sum_{j} l_{i j} μ_{j}$
ist, Yi die gemessenen Photonenanzahl in der Projektion i ist und lij die Schnittlänge des Strahls durch das Voxel j in der Projektion i ist und
wobei β mittels folgender Formel bestimmt wird: $β = \frac{Δ \tilde{μ} \cdot \bar{\sum_{i} Y_{i} l_{i j}}}{\bar{μ} \cdot \frac{\partial}{\partial μ_{j}} U (μ) {|_{Δ μ_{j} = Δ \tilde{μ}} - Δ \tilde{μ} \cdot \bar{μ} \cdot \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ) |}_{^{^{^{_{^{^{Δ μ_{_{j}} = Δ \tilde{μ}}}}}}}}};$
wobei $\bar{\sum_{i} Y_{i} l_{i j}}$
der Mittelwert über alle Voxel des Volumenelementes über alle i und j ist.
Diese Formel entspricht der zuvor genannten Formel zur Bestimmung von µ_j , wobei die Zielfunktion L(µ) und die eingangs beschriebene Gleichung für Φ(µ) in die Formel eingesetzt wurden.
Der Steuerparameter β kann vor der ersten Iteration unter Zuhilfenahme der Geman-Prior-Funktion U(µ) bestimmt werden. Zu diesem Zeitpunkt ist jedoch die Differenz benachbarter Schwächungskoeffizienten Δµ_jk noch nicht bekannt. Werden die Messdaten mittels Röntgenstrahlung ermittelt, kann angenommen werden, dass das Rauschen in einer Projektion konstant ist. In diesem Fall wird daher der Ausdruck Δµ_jk durch den, vorzugsweise maximalen, Rauschwert Δµ̃ ersetzt. Wurden die Messdaten mittels eines Magnetresonanztomographen erzeugt, kann das Rauschen ortsabhängig sein, wodurch der Ausdruck Δµ_jk in der Geman-Prior-Funktion U(µ) durch einen ortsabhängigen Rauschschätzwert Δµ̃_j ersetzt werden kann. Wie zuvor erwähnt wurde, kann die Differenz der Schwächungskoeffizienten Δµ_jk zuvor durch herkömmliche Rückprojektionsverfahren ermittelt werden. Dabei kann β genauer bestimmt werden.
Alternativ hierzu kann das Rekonstruktions-Verfahren gemäß einer zweiten Ausführungsform der Erfindung eine iterative konvexe Optimierungs-Aktualisierungs-Formel verwenden. Hierfür wird der Schritt des iterativen Berechnens des Schätzwertes des Schwächungskoeffizienten µ_j eines Voxels j mittels folgender Formel durchgeführt: $μ_{j}^{n + 1} = μ_{j}^{n} + μ_{j}^{n} \frac{\sum_{i} l_{i j} [d_{i} e^{- 〈 l_{i}, μ 〉} - Y_{i}] - β \frac{\partial}{\partial μ_{j}} U (μ^{n})}{\sum_{i} l_{i j} 〈 l_{i}, μ^{n} 〉 d, e^{- 〈 l_{i}, μ^{n} 〉} + β μ_{j}^{n} \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ^{n})};$
wobei β ein Optimierungsparameter ist, $〈 l_{i}, μ 〉 = \sum_{j} l_{i j} μ_{j}$
ist, Yi die gemessenen Photonenanzahl in der Projektion i ist und l_ij die Schnittlänge des Strahls durch das Voxel j in der Projektion i ist und
wobei β mittels folgender Formel bestimmt wird: $β = \frac{Δ \tilde{μ} \cdot \sum_{i} l_{i j} 〈 l_{i}, \bar{μ} 〉 d, e^{- 〈 l_{i}, \bar{μ} 〉}}{\bar{μ} \cdot \frac{\partial}{\partial μ_{j}} U (μ) {|_{Δ μ_{j} = Δ \tilde{μ}} - Δ \tilde{μ} \cdot \bar{μ} \cdot \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ) |}_{^{^{^{_{^{^{Δ μ_{j} = Δ \tilde{μ}}}}}}}}};$
wobei $〈 l_{i}, \bar{μ} 〉 = \sum_{i} l_{i j} \bar{μ}$
ist.
Der Steuerparameter β kann vor der ersten Iteration unter Zuhilfenahme der Geman-Prior-Funktion U(µ) bestimmt werden. Zu diesem Zeitpunkt ist jedoch die Differenz benachbarter Schwächungskoeffizienten Δµ_jk noch nicht bekannt. Werden die Messdaten mittels Röntgenstrahlung ermittelt, kann angenommen werden, dass das Rauschen in einer Projektion konstant ist. In diesem Fall wird daher der Ausdruck Δµ_jk durch Δµ̃ ersetzt. Wurden die Messdaten mittels eines Magnetresonanztomographen erzeugt, kann das Rauschen ortsabhängig sein, wodurch der Ausdruck Δµ_jk in der Geman-Prior-Funktion U(µ) durch einen ortsabhängigen Rauschschätzwert Δµ̃_j ersetzt werden kann. Wie zuvor erwähnt wurde, kann die Differenz der Schwächungskoeffizienten Δµ_jk zuvor durch herkömmliche Rückprojektionsverfahren ermittelt werden. Hierdurch kann β genauer bestimmt werden.
Die ermittelten Schwächungskoeffizienten µ_j werden auf einer Anzeigeeinrichtung, beispielsweise einem Bildschirm, als Grauwerte oder in einer Falschfarbendarstellung dargestellt. Die ermittelten Schwächungskoeffizienten können auch auf einem Ausdruck, beispielsweise als Intensitäten oder Falschfarben dargestellt werden.
Das beschriebene Verfahren kann auch zur Kontrastverbesserung verwendet werden, falls zuvor eine dreidimensionale Rekonstruktion, beispielsweise durch eine ungefilterte RückProjektion durchgeführt wurde. Die ungefilterte Rückprojektion stellt Initialisierungswerte µ_j ⁰ bereit, so dass eine schnelle Konvergenz erzielt werden kann.
Die Erfindung betrifft auch eine Gewebedichte-Ermittlungseinrichtung mit einem Emitter, der eine Strahlung in ein Gewebe abgibt, und einem Detektor, der die Strahlung erfasst, nachdem sie das Gewebe durchlaufen hat. Die Gewebedichte-Ermittlungseinrichtung umfasst eine Rauschwert-Schätzeinrichtung, die zum Schätzen des Rauschschätzwertes Δµ̃ des Schwächungskoeffizienten µ eines Voxel ausgebildet ist. Die Rauschwert-Schätzeinrichtung ist vorzugsweise so ausgebildet, dass sie den erwarteten Rauschwert Δµ̃ des Schwächungskoeffizienten µ eines Voxel aus der Gewebedichte und der Information über die abgegebene Strahlung, insbesondere aus der Brustdicke und der Röhrenladung (mAs), zu bestimmen. Der Rauschpegel kann im Allgemeinen für eine beliebige Brustdicke durch Auswahl der Röhrenladung und der Photonenenergie (keV) konstant gehalten werden. Der Rauschschätzwert Δµ̃ muss kleiner als der kleinste zu erhaltende Gradient sein. Der Rauschschätzwert Δµ̂ ist der maximale Rauschschätzwert, der in einem Bild erwartet wird, wenn keine Glättung erfolgt.
Die Gewebedichte-Ermittlungseinrichtung umfasst ferner eine Gewebedämpfungs-Mittelwert-Schätzeinrichtung, die zum Schätzen eines mittleren Gewebeschwächungskoeffizienten µ̅ ausgebildet ist, wobei die Gewebedämpfungs-Mittelwert-Schätzwerteinrichtung den mittleren Gewebeschwächungskoeffizienten µ̅ vorzugsweise auf Grundlage einer Gewebedichteschätzung und der Energie der vom Emitter emittierten Strahlung ermittelt. Der Gewebedämpfungsschätzwert µ̅ kann beispielsweise aus der Photonenenergie und der Brustdichte ermittelt werden, die aus segmentierten Originalprojektionen unter Verwendung bestehender Verfahren geschätzt wird. Der erwartete Rauschwert Δµ̃ und/oder der mittlere Gewebedämpfungsschätzwert µ̅ können mittels Regressionsverfahren ermittelt werden, wie zuvor erläutert wurde.
Die Gewebedichte-Ermittlungseinrichtung weist ferner eine Prior-Funktion-Bestimmungseinrichtung auf, die dazu ausgebildet ist, Parameter einer Prior-Funktion U(µ) aus dem Rauschschätzwert Δµ̃ und dem mittleren Gewebeschwächungskoeffizienten µ̅ zu berechnen. Die Gewebedichte-Ermittlungseinrichtung umfasst ferner eine Iterationseinrichtung, die dazu ausgebildet ist, den Erwartungswert eines Schwächungskoeffizienten µ_j eines Voxel j aus mittels des Detektors (18) erfassten Projektionsdaten, dem Rauschschätzwert Δµ̃, dem mittleren Gewebedämpfungsschätzwert µ̅ und der Prior-Funktion U(µ) zu berechnen. Der Rauschschätzwert Δµ̃ ist derjenige Rauschwert, der in einem Bild erwartet wird, wenn keine Glättung mittels der Prior-Funktion erfolgt.
Der bestimmte Schwächungskoeffizient µ_j kann als Grauwert oder als eine Falschfarbe auf einer Anzeigeeinrichtung dargestellt werden. Somit ist es möglich, dichteres Gewebe von weniger dichtem Gewebe zu unterscheiden, da der Schwächungskoeffizient mit der Gewebedichte zunimmt, wobei gleichzeitig ein höherer Signal/Rausch-Abstand und Kontrast/Rausch-Abstand als mit Verfahren des Standes der Technik ermöglicht werden. Aufgrund des verbesserten Signal/Rausch-Abstandes und Kontrast/Rausch-Abstandes kann die Dosis des Strahlungsemitters reduziert werden, wodurch die Nebenwirkungen reduziert werden.
Die Prior-Funktion-Bestimmungseinrichtung kann dazu ausgebildet sein, die Parameter einer Geman-Prior-Funktion zu berechnen. Die Geman-Prior-Funktion U(µ) kann folgende Gleichung aufweisen: $U (μ) = \sum_{j} \sum_{k} \frac{Δ μ_{j k}^{2} σ^{m}}{2 {(\sqrt{Δ μ_{j k}^{2} / 2 + σ^{2} / 2})}^{m}};$
wobei β ein Optimierungsparameter und m eine Konstante ist und wobei µ_k und µ_j Schwächungskoeffizienten benachbarter Voxel sind, die je einem Gewebeabschnitt zugeordnet sind, und j und k lineare Voxel-Indizes über drei Dimensionen sind. Die Werte der Schwächungskoeffizienten µ_k und µ_j werden mit jeder Iteration aktualisiert. D.h., bei der ersten Iteration werden diese Werte mit Schätzwerten vorbelegt, beispielsweise mit dem mittleren Gewebeschwächungswert. Bei den anschließenden Iterationen werden für die Schwächungskoeffizienten µ_k und µ_j diejenigen Werte verwendet, die bei der vorangegangenen Iteration bestimmt wurden. Es kann beispielsweise eine Nachbarschaft von 3x3x3 Voxeln ausgewertet werden.
Bevor das Bild rekonstruiert wurde, sind die Differenzen der Schwächungskoeffizienten Δµ_jk nicht bekannt. Wird das Rekonstruktionsverfahren auf Daten angewendet, die mittels Röntgenstrahlen erzeugt wurden, wird angenommen, dass das Rauschen in der Projektion ortsunabhängig ist. Daher wird zur Bestimmung des optimalen Wertes für den Wert σ in der Geman-Prior-Funktion der Term Δµ_jk durch den Rauschschätzwert, vorzugsweise den maximalen Rauschschätzwert, Δµ̃ ersetzt. Wurden die Bilddaten mittels eines Magnetresonanztomographen aufgenommen, kann das Rauschen ortsabhängig sein. In diesem Fall kann der Wert Δµ̃_j vom Ort des Voxels j abhängig sein.
Es ist jedoch auch möglich, die Werte µ_j und µ_k und somit Δµ_jk mittels herkömmlicher Rekonstruktionsverfahren zu bestimmen, beispielsweise mittels einer linearen und/oder gefilterten Rückproduktion. Hierdurch kann σ genauer bestimmt werden.
Die Gewebedichte-Ermittlungseinrichtung ist vorzugsweise dazu ausgebildet, dass, falls die Geman-Prior-Funktion eine konvexe Form aufweisen soll, σ so gewählt wird, dass die erste Ableitung der Geman-Prior-Funktion ihre stärkste Krümmung bei Werten von Δµ erreicht, die dem zu Rauschschätzwert Δµ̃ entsprechen, und vorzugsweise m ≤ 16/17 ist. Eine konvexe Form der Geman-Prior-Funktion ist bevorzugt, damit die Kanten und die Ränder der Gewebestrukturen beim Reduzieren des Rauschens erhalten bleiben. Die Geman-Prior-Funktion-Bestimmungseinrichtung ist vorzugsweise dazu ausgebildet, dass, falls die Geman-Prior-Funktion eine nicht-konvexe Form aufweisen soll, σ so gewählt wird, dass die erste Ableitung der Geman-Prior-Funktion ihr Maximum bei Werten von Δµ erreicht, die dem Rauschschätzwert Δµ̃ entsprechen, und vorzugsweise m = 2 ist.
Die Iterationseinrichtung ist dazu ausgebildet, den Schätzwert des Schwächungskoeffizienten eines Voxels zur dreidimensionalen Bildrekonstruktion mittels des folgenden Maximum-A-Posteriori-Verfahren zu bestimmen, das ein Newton-basiertes Aktualisierungsverfahren mit einer einfachen quadratischen Glättungs-Prior-Funktion U(µ) ist: $μ_{j}^{n + 1} = μ_{j}^{n} + \frac{μ_{j}^{n}}{\sum_{i} Y_{i} l_{i j} + β μ_{j}^{n} \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ^{n})} \frac{\partial}{\partial μ_{j}} Φ (μ^{n});$
wobei β ein Optimierungsparameter ist, Φ(µⁿ) = L(µ) - βU (µ) ist, L(µ) eine Zielfunktion ist, Yi die gemessene Photonenanzahl in der Ebene i ist und lij die Schnittlänge des Strahls durch das Voxel j in der Ebene i ist. Der Anfangswert µ_j ⁰ des iterativen Verfahrens kann so gewählt werden, wie zuvor hinsichtlich des Verfahrens ausgeführt wurde.
Die Zielfunktion L(µ) lautet gemäß der ersten Ausführungsform der Erfindung: $L (μ) = \sum_{i} {- d_{i} e^{- 〈 l_{i}, μ 〉} - Y_{i} 〈 l_{i}, μ 〉} + c_{;}$
wobei $〈 l_{i}, μ 〉 = \sum_{j} l_{i j} μ_{j}$
ist, i ein Projektionsindex ist, j der Voxelindex ist, d_i die erwartete Photonenanzahl ist, die aus der Strahlungsquelle entlang der Projektion i austritt, Y_I die erfasste Photonenanzahl in der Projektion i ist und c eine Konstante ist.
Die Iterationseinrichtung ist gemäß dieser Ausführungsform dazu ausgebildet, β mittels folgender Formel zu bestimmen: $β = \frac{Δ \tilde{μ} \cdot \bar{\sum_{i} Y_{i} l_{i j}}}{\bar{μ} \cdot \frac{\partial}{\partial μ_{j}} U (μ) {|_{Δ μ_{j} = Δ \tilde{μ}} - Δ \tilde{μ} \cdot \bar{μ} \cdot \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ) |}_{^{^{^{_{^{^{Δ μ_{_{j}} = Δ \tilde{μ}}}}}}}}};$
wobei $\bar{\sum_{i} Y_{i} l_{i j}}$
der Mittelwert über alle Voxel des Volumenelementes über alle i und j ist.
Der Steuerparameter β kann vor der ersten Iteration unter Zuhilfenahme der Geman-Prior-Funktion U(µ) bestimmt werden. Zu diesem Zeitpunkt ist jedoch die Differenz benachbarter Schwächungskoeffizienten Δµ_jk noch nicht bekannt. Werden die Messdaten mittels Röntgenstrahlung ermittelt, kann angenommen werden, dass das Rauschen in einer Projektion konstant ist. In diesem Fall wird daher der Ausdruck Δµ_jk durch Δµ̃ ersetzt. Wurden die Messdaten mittels eines Magnetresonanztomographen erzeugt, kann das Rauschen ortsabhängig sein, wodurch der Ausdruck Δµ_jk in der Geman-Prior-Funktion U(µ) durch einen ortsabhängigen Rauschschätzwert Δµ̃_j ersetzt werden kann. Wie zuvor erwähnt wurde, kann die Differenz der Schwächungskoeffizienten Δµ_jk zuvor durch herkömmliche Rückprojektionsverfahren ermittelt werden, wodurch β genauer bestimmt werden kann.
Die Iterationseinrichtung ist dazu ausgebildet, die dreidimensionale Rekonstruktion mittels einer iterativen Gradienten-Optimierungs-Aktualisierungs-Formel durchzuführen. Die iterative Berechnung des Schätzwertes des Schwächungskoeffizienten µ_j eines Voxels j wird gemäß der ersten Ausführungsform der Erfindung hierbei mittels folgender Formel durchgeführt: $μ_{j}^{n + 1} = μ_{j}^{n} + μ_{j}^{n} \frac{\sum_{i} l_{i j} [d_{i} e^{- 〈 l_{i}, μ 〉} - Y_{i}] - β \frac{\partial}{\partial μ_{j}} U (μ^{n})}{\sum_{i} Y_{i} l_{i j} + β μ_{j}^{n} \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ^{n})};$
wobei β ein Optimierungsparameter ist, $〈 l_{i}, μ 〉 = \sum_{i} l_{i j} μ_{j}$
ist, Yi die gemessenen Photonenanzahl in der Projektion i ist und lij die Schnittlänge des Strahls durch das Voxel j in der Projektion i ist und
wobei β vorzugsweise mittels folgender Formel bestimmt wird: $β = \frac{Δ \tilde{μ} \cdot \bar{\sum_{i} Y_{i} l_{i j}}}{\bar{μ} \cdot \frac{\partial}{\partial μ_{j}} U (μ) {|_{Δ μ_{j} = Δ \tilde{μ}} - Δ \tilde{μ} \cdot \bar{μ} \cdot \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ) |}_{^{^{^{_{^{^{Δ μ_{j} = Δ \tilde{μ}}}}}}}}};$
wobei $\bar{\sum_{i} Y_{i} l_{i j}}$
der Mittelwert über alle Voxel des Volumenelementes über alle i und j ist.
Diese Formel entspricht der zuvor genannten Formel zur Bestimmung von µ_j , wobei die Zielfunktion L(p), und die eingangs erläuterte Gleichung für Φ(µ) eingesetzt wurden.
Der Steuerparameter β kann vor der ersten Iteration unter Zuhilfenahme der Geman-Prior-Funktion U(µ) bestimmt werden. Zu diesem Zeitpunkt ist jedoch die Differenz benachbarter Schwächungskoeffizienten Δµ_jk noch nicht bekannt. Werden die Messdaten mittels Röntgenstrahlung ermittelt, kann angenommen werden, dass das Rauschen in einer Projektion konstant ist. In diesem Fall wird daher der Ausdruck Δµ_jk durch Δµ̃ ersetzt. Wurden die Messdaten mittels eines Magnetresonanztomographen erzeugt, kann das Rauschen ortsabhängig sein, wodurch der Ausdruck Δµ_jk in der Geman-Prior-Funktion U(µ) durch einen ortsabhängigen Rauschschätzwert Δµ̃_j ersetzt werden kann. Wie zuvor erwähnt wurde, kann die Differenz der Schwächungskoeffizienten Δµ_jk zuvor durch herkömmliche Rückprojektionsverfahren ermittelt werden, wodurch β genauer bestimmt werden kann.
Alternativ hierzu kann die Iterationseinrichtung gemäß einer zweiten Ausführungsform der Erfindung dazu ausgebildet sein, die dreidimensionale Rekonstruktion mittels einer iterativen konvexen Optimierungs-Aktualisierungs-Formel durchzuführen. Hierfür wird das iterative Berechnen des Schätzwertes des Schwächungskoeffizienten µ_j eines Voxels j mittels folgender Formel durchgeführt: $μ_{j}^{n + 1} = μ_{j}^{n} + μ_{j}^{n} \frac{\sum_{i} l_{i j} [d_{i} e^{- 〈 l_{i}, μ 〉} - Y_{i}] - β \frac{\partial}{\partial μ_{j}} U (μ^{n})}{\sum_{i} l_{i j} 〈 l_{i}, μ^{n} 〉 d, e^{- 〈 l_{i}, μ^{n} 〉} + β μ_{j}^{n} \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ^{n})};$
wobei β ein Optimierungsparameter ist, $〈 l_{i}, μ 〉 = \sum_{i} l_{i j} μ_{j}$
ist, Y_i die gemessenen Photonenanzahl in der Projektion i ist und lij die Schnittlänge des Strahls durch das Voxel j in der Projektion i ist und
wobei β mittels folgender Formel bestimmt wird: $β = \frac{Δ \tilde{μ} \cdot \sum_{i} l_{i j} 〈 l_{i}, \bar{μ} 〉 d_{i} e^{- 〈 l_{i}, μ^{n} 〉}}{\bar{μ} \cdot \frac{\partial}{\partial μ_{j}} U (μ) {|_{Δ μ_{j} = Δ \tilde{μ}} - Δ \tilde{μ} \cdot \bar{μ} \cdot \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ) |}_{^{^{^{_{^{^{Δ μ_{j} = Δ \tilde{μ}}}}}}}}};$
wobei $〈 l_{i}, \bar{μ} 〉 = \sum_{i} l_{i j} \bar{μ}$
ist.
Der Steuerparameter β kann vor der ersten Iteration unter Zuhilfenahme der Geman-Prior-Funktion U(µ) bestimmt werden. Zu diesem Zeitpunkt ist jedoch die Differenz benachbarter Schwächungskoeffizienten Δµ_jk noch nicht bekannt. Werden die Messdaten mittels Röntgenstrahlung ermittelt, kann angenommen werden, dass das Rauschen in einer Projektion konstant ist. In diesem Fall wird daher der Ausdruck Δµ_jk durch Δµ̃ ersetzt. Wurden die Messdaten mittels eines Magnetresonanztomographen erzeugt, kann das Rauschen ortsabhängig sein, wodurch der Ausdruck Δµ_jk in der Geman-Prior-Funktion U(µ) durch einen ortsabhängigen Rauschschätzwert Δµ̃_j ersetzt werden kann. Wie zuvor erwähnt wurde, kann die Differenz der Schwächungskoeffizienten Δµ_jk zuvor durch herkömmliche Rückprojektionsverfahren ermittelt werden, um β genauer zu parametrisieren.
Die Erfindung betrifft ferner ein Computerprogramm, das Programmcodemittel aufweist, die in einem Computer geladen werden und die zum Ausführen aller Schritte des zuvor beschriebenen Verfahrens ausgebildet sind.
Die Erfindung betrifft ferner ein medizinisches System mit der zuvor beschriebenen Gewebedichte-Ermittlungseinrichtung und Positionierungsmitteln, um zumindest einen Teil eines Organs eines Patienten räumlich gegenüber dem Emitter und dem Detektor festzulegen, und einer Anzeigeeinrichtung, die dazu ausgebildet ist, Repräsentationen der Werte µ_j anzuzeigen. Das Organ kann die Brust sein. Es ist auch möglich die gesamte Bauchregion mittels der Positionierungsmittel örtlich festzulegen und ein Bild des Bauchbereichs aufzunehmen. Darüber hinaus können mittels der Positionierungsmittel ein Körperbereich mit einer Verletzung, beispielsweise einer Fraktur, örtlich festgelegt werden, um die innere Verletzung zu untersuchen. Die Anzeigeeinrichtung kann ein Bildschirm oder ein Drucker sein. Die Repräsentationen der Gewebeschwächungskoeffizienten können Intensitäten, Grauwerte oder Falschfarben sein.
Die Erfindung wird nun unter Bezugnahme auf die beigefügten Zeichnungen detaillierter erläutert. Es zeigen:

1 ein Röntgengerät;
2 die physikalischen Grundlagen der Schwächung von Röntgenstrahlung in Gewebe;
3 ein Flussdiagramm des erfindungsgemäßen Verfahrens;
4 Details einer erfindungsgemäßen Gewebedämpfungswert-Ermittlungseinrichtung;
5 und 6 die erste und die zweite Ableitung der Geman-Prior-Funktion;
7 bis 9 einen Schnitt durch die Brust, der mit unterschiedlichen bildgebenden Verfahren rekonstruiert wurde; und
10 bis 12 Bildrekonstruktionen des kleinsten sichtbaren Gefäßes, das mittels unterschiedlicher Rekonstruktionsverfahren rekonstruiert wurde.

1 zeigt ein Mammographiegerät als eine Ausführungsform eines Röntgengeräts. Das Mammographiegerät weist eine schwenkbare Röntgenstrahlröhre 2 auf, die Röntgenstrahlen 4 erzeugt. Die Röntgenstrahlröhre 2 ist innerhalb eines Winkelbereichs von -φ₁ bis +φ₂ um einen Winkel α um den Mittelpunkt M schwenkbar. Jeder Projektion i ist ein Winkel α der Röntgenstrahlröhre 2 zur Symmetrieachse 6 zugeordnet. Das Mammographiegerät weist eine erste Kompressionsplatte 8 und eine zweite Kompressionsplatte 12 auf, zwischen denen das zu untersuchende Organ 10, d. h. die weibliche Brust, zur Untersuchung eingebettet wird. Das Röntgengerät weist einen Detektor 18 mit einer Mehrzahl von Messelementen 14, 16 auf, die auf dem Detektor 18 matrixartig angeordnet sind.
Zum Aufnehmen einer zweidimensionalen Projektion wird die Röntgenstrahlröhre 2 in die erste Stellung i = 1 geschwenkt. Anschließend wird ein Röntgenstrahl 4 erzeugt, der durch die erste und zweite Kompressionsplatte 8, 12 sowie durch das zu untersuchende Organ 10 verläuft und von den Messelementen 14, 16 des Detektors 18 erfasst wird. Die Röntgenstrahlröhre 2 wird in die zweite Stellung i = 2 geschwenkt und es wird abermals ein Röntgenstrahl 4 erzeugt, der das zu untersuchende Organ 10 unter einem anderen Winkel, und folglich in einer anderen Projektion i durchläuft. Die Röntgenstrahlen 4 werden nach dem Durchlaufen des zu untersuchenden Organs 10 von den Messelementen 14, 16 des Detektors 18 erfasst. Das Mammographiegerät kann n Projektionen des zu untersuchenden Organs 10 aufnehmen. Die Anzahl der aufgenommenen Projektionen n befindet sich vorzugsweise im Bereich von 20 bis 50 Projektionen.
Die zweidimensionalen Projektionen werden von der Gewebedichte-Ermittlungseinrichtung 20 in eine dreidimensionale Darstellung gewandelt. Mittels der Tastatur 24 kann ein Schnitt durch die dreidimensionale Darstellung ausgewählt werden, der auf der Anzeigeeinrichtung 22 dargestellt wird.
Es wird auf 2 Bezug genommen, um die physikalischen Grundlagen der Ermittlung der Gewebedichte mittels Röntgenstrahlen zu erläutern.
Die Grundlagen zur Messung von Gewebe mittels Röntgenstrahlung sind detailliert in Imaging Systems for Medical Diagnostics, Arnulf Oppelt, Publicis Corporate Publishing, Erlangen, ISBN 3-89578-226-4, Kapitel 9 und 10 erläutert.
Die Röntgenstrahlröhre 2 emittiert einen Röntgenstrahl 4, der vom Detektor 18 erfasst wird. Der Röntgenstrahl 4 wird im Gewebe 10 gedämpft. Jedes Gewebeelement 30, das als Voxel auf der Anzeigeeinrichtung 22 dargestellt wird, weist einen individuellen Dämpfungskoeffizienten µ_j auf. Die Dämpfung im Gewebe 10 erfolgt exponentiell. Die vom jeweiligen Messelement 14, 16 des Detektors 18 erfasste Intensität lässt sich wie folgt berechnen: $J (r, φ) = J_{0} \cdot e^{- \int μ d s}$
Amplitudendiskret kann die vom jeweiligen Messelement 14, 16 des Detektors 18 erfasste Intensität J wie folgt dargestellt werden: $J (r, φ) = J_{0} \cdot e^{- \sum_{j = 1}^{J} l_{ij} μ_{j}};$
wobei µ_j der Dämpfungskoeffizient eines Gewebeabschnitts bzw. Voxels ist, J₀ die Intensität des Röntgenstrahls vor dem Eintritt in das Gewebe ist, i die jeweilige Projektion ist, j der Voxelindex ist und l_ij die Schnittlänge ist. Die Schnittlänge lij bezeichnet folglich die Länge, durch die der Röntgenstrahl das Voxel j mit dem Dämpfungskoeffizient µ_j durchquert.
Die Röntgenstrahlröhre 2 wird, wie eingangs beschrieben, um einen Winkel α geschwenkt, um anschließend weitere zweidimensionale Projektionen aufzunehmen, die anschließend von der Gewebedichte-Ermittlungseinrichtung 20 zu einem dreidimensionalen Abbild weiterverarbeitet werden.
Der Schwächungskoeffizient µ_j korreliert mit der Gewebedichte. Wie eingangs erwähnt wurde, umfasst die weibliche Brust Drüsengewebe, Fettgewebe, Bindegewebe und Blutgefäße deren Röntgenstrahlungsdämpfungskoeffizienten sehr ähnlich sind. Zusätzlich kommt erschwerend hinzu, dass sich der Röntgenstrahlenschwächungskoeffizient von krebsbefallenem Gewebe zahlenmäßig in der Nähe der Röntgenstrahlenschwächungskoeffizienten der zuvor genannten Gewebearten befindet. Folglich muss für die Tomosynthese im Bereich der Mammographie ein besonders rauscharmes Rekonstruktionsverfahren bereitgestellt werden, das darüber hinaus die Kanten und Ränder der zuvor genannten Gewebearten deutlich darstellt.
Der Detektor 18 ermittelt die Intensität des Röntgenstrahls 4 mittels der matrixartig angeordneten Messelemente 14 und 16. Ein Messelement 14, 16 kann nicht den Gewebeschwächungskoeffizient µ_j darstellen, sondern lediglich einen Wert ermitteln, der angibt, wie stark der Röntgenstrahl 4 im Gewebe 10 gedämpft bzw. abgeschwächt wurde, sofern die anfängliche Intensität J₀ des Röntgenstrahls bekannt ist. Da jedoch mehrere Projektionen i unter unterschiedlichen Winkeln α des Gewebes 10 aufgenommen werden, wobei das Gewebe 10 in seiner Orientierung verbleibt, ist es möglich, ein dreidimensionales Bild in der Gewebedämpfungswert-Ermittlungseinrichtung 20 zu ermitteln.
Nachdem eine Mehrzahl zweidimensionaler Schnittbilder in unterschiedlichen Projektionen i ermittelt worden sind, wird das erfindungsgemäße Rekonstruktionsverfahren begonnen, das jetzt detaillierter unter Bezugnahme auf 3 beschrieben wird.
Am Anfang des Iterationsverfahrens wird der Iterationsindex n mit 0 initialisiert. Im Schritt S1 wird der Rauschschätzwert Δµ̃ des Schwächungskoeffizienten µ eines Voxels geschätzt. Die Schätzung erfolgt bei einem gegebenen Röntgenspektrum als eine Funktion der Brustdicke und der Röhrenladung (mAs) unter Verwendung von Regressionsverfahren und zuvor durchgeführten statistischen Untersuchungen. Im Allgemeinen kann der Rauschpegel für eine beliebige Brustdicke durch Auswahl der Röhrenladung und der Photonenenergie (keV) konstant gehalten werden. Der Rauschschätzwert Δµ̃ muss kleiner als der kleinste zu erhaltende Gradient sein. Der Rauschschätzwert Δµ̂ ist der maximale Rauschwert, der in einem Bild erwartet wird, wenn keine Glättung mittels einer Prior-Funktion erfolgt.
Das Verfahren fährt anschließend mit dem Schritt S2 fort. Im Schritt S2 wird der mittlere Gewebedämpfungsschätzwert bzw. mittlere Gewebeschwächungsschätzkoeffizient µ̅ geschätzt. Der Gewebeschwächungsschätzkoeffizient µ̅ wird unter Verwendung von Regressionsverfahren und statistischen Untersuchungen aus der Photonenenergie und der Brustdichteschätzung geschätzt. Der Gewebeschwächungsschätzkoeffizient kann auch die aus segmentierten Originalprojektionen unter Verwendung bekannter Verfahren aus der Photonenenergie und der Brustdickeschätzung erhalten werden. Solche Verfahren sind beispielsweise in O. Alonzo-Proulx, A. H. Tyson, G. E. Marwdsley, M. J. Yaffe, Effect of Tissue Thickness Variation in Volumeric Breast Density Estimation, Proc. IWDM 2008, pp. 659 - 666 beschrieben. Der Rauschschätzwert Δµ̃ und die Brustdichte können aus dem segmentierten rekonstruierten Brustvolumen µ⁰ berechnet werden, das zum Initialisieren verwendet wird.
Im Schritt S3 wird anschließend geprüft, ob die Geman-Prior-Funktion konvex sein soll. Als Geman-Prior-Funktion U(µ) wird folgende Gleichung verwendet: $U (μ) = \sum_{j} \sum_{k} \frac{Δ μ_{j k}^{2} σ^{m}}{2 {(\sqrt{Δ μ_{j k}^{2} / 2 + σ^{2} / 2})}^{m}};$
wobei σ ein Optimierungsparameter und m eine Konstante ist.
Bevor das Bild rekonstruiert wurde sind die Differenzen Δµ_jk nicht bekannt. Wird das Rekonstruktionsverfahren auf Daten angewendet, die mittels Röntgenstrahlen erzeugt wurden, wird angenommen, dass das Rauschen in der Projektion ortsunabhängig ist. Daher wird zur Bestimmung des optimalen Wertes für den Wert σ in der Geman-Prior-Funktion der Term Δµ_jk durch den Rauschschätzwert, vorzugsweise durch den maximalen Rauschschätzwert, Δµ̃ ersetzt.
Es ist jedoch auch möglich, die Werte µ_j und µ_k und somit Δµ_jk mittels herkömmlicher Rekonstruktionsverfahren zu bestimmen, beispielsweise mittels einer linearen und/oder gefilterten Rückproduktion. Hierdurch kann der Parameter σ genauer bestimmt werden.
Es ist bevorzugt, dass die Geman-Prior-Funktion eine konvexe Form aufweist, um Kanten und Gewebebereichsgrenzen in der Darstellung zu erhalten und die Konvergenz sicherzustellen, während das Rauschen unterdrückt wird.
Wurde im Schritt S3 festgelegt, dass die Geman-Prior-Funktion konvex sein soll, wird zum Schritt S4 fortgefahren. σ wird so gewählt, dass die erste Ableitung der Geman-Prior-Funktion ihre stärkste Krümmung bei Werten von Δµ̃ erreicht, die dem erwarteten Rauschpegel entsprechen. Die Konstante m wird wie folgt festgelegt: m ≤ 16/17.
Wird im Schritt S3 bestimmt, dass die Geman-Prior-Funktion nicht konvex sein soll, wird zum Schritt S5 fortgefahren. Hierbei wird σ so gewählt, dass die erste Ableitung der Geman-Prior-Funktion ihr Maximum bei Werten von Δµ̃ erreicht, die dem erwarteten Rauschpegel entsprechen. Die Konstante m wird wie folgt festgelegt: m = 2.
Nach dem Schritt S4 oder S5 fährt das Verfahren mit dem Schritt S6 fort und bestimmt den Steuerparameter β. Der Steuerparameter β ist stets größer 0. Der Steuerparameter β wird mittels folgender Formel bestimmt: $β = \frac{Δ \tilde{μ} \cdot \bar{\sum_{i} Y_{i} l_{i j}}}{\bar{μ} \cdot \frac{\partial}{\partial μ_{j}} U (μ) {|_{Δ μ_{j} = Δ \tilde{μ}} - Δ \tilde{μ} \cdot \bar{μ} \cdot \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ) |}_{^{^{^{_{^{^{Δ μ_{j} = Δ \tilde{μ}}}}}}}}};$
wobei lij die Schnittlänge des Strahls durch das Voxel j in der Ebene i ist, $\bar{\sum_{i} Y_{i} l_{i j}}$
der Mittelwert über alle Voxel des Volumenelementes über alle i und j ist.
Die Bestimmung der Schnittlänge lij wurde zuvor anhand der 2 erläutert. Yi ist die gemessene Photonenanzahl in der Projektion i.
Der Steuerparameter β kann vor der ersten Iteration unter Zuhilfenahme der Geman-Prior-Funktion U(µ) bestimmt werden. Zu diesem Zeitpunkt ist jedoch die Differenz benachbarter Schwächungskoeffizienten Δµ_jk noch nicht bekannt. Werden die Messdaten mittels Röntgenstrahlung ermittelt, kann angenommen werden, dass das Rauschen in einer Projektion i konstant ist. In diesem Fall wird daher der Ausdruck Δµ_jk durch Δµ̃ ersetzt. Wurden die Messdaten mittels eines Magnetresonanztomographen erzeugt, kann das Rauschen ortsabhängig sein, wodurch der Ausdruck Δµ_jk in der Geman-Prior-Funktion U(µ) durch einen ortsabhängigen Rauschschätzwert Δµ̃_j ersetzt werden kann. Wie zuvor erwähnt wurde, kann die Differenz der Schwächungskoeffizienten Δµ_jk zuvor durch herkömmliche Rückprojektionsverfahren ermittelt werden.
Anschließend werden im Schritt S7 die Schwächungskoeffizienten µ_j ⁿ und µkⁿ , die benachbarten Voxel zugeordnet sind, bestimmt. Die Indizes j und k sind lineare Indizes über drei Dimensionen. Die Werte der Schwächungskoeffizienten µ_j ⁿ und µ_k ⁿ können diejenigen sein, die bei einer vorangegangenen Iteration ermittelt wurden. Diese Schwächungskoeffizienten können bei der ersten Iteration mit einem, beispielsweise für alle Voxel identischen, Schätzwert initialisiert werden. Die Schwächungskoeffizienten µ_j und µ_k können auch mittels herkömmlicher dreidimensionaler Rekonstruktionsverfahren bestimmt werden, um die Konvergenz zu beschleunigen.
Anschließend werden im Schritt S8 die Differenzen der Indizes Δu_ik berechnet. Es gilt: $k \in N_{j};$
Die Differenz Δµ_ik kann beispielsweise in einer Nachbarschaft von 3x3x3 Voxeln bestimmt werden.
Der Schritt S9, der nach dem Schritt S8 ausgeführt wird, bestimmt die Geman-Prior-Funktion aus Δµ_jk , σ, m, j und k. Die Geman-Prior-Funktion U(µ) wird als Glättungs-Prior-Funktion verwendet, die Differenzen in Werten benachbarter Pixel bestraft bzw. penalisiert.
Nach dem Schritt S9 fährt das Verfahren mit dem Schritt S10 fort. Im Schritt S10 wird der Schwächungskoeffizient µ_j iterativ durch ein Maximum-A-Posterior-Verfahren bestimmt, das ein Newton-basiertes Aktualisierungsverfahren mit einer einfachen quadratischen Glättungs-Prior-Funktion U(µ) verwendet.
Die Gleichung des iterativen Maximum-A-Posterior-Verfahrens lautet wie folgt: $μ_{j}^{n + 1} = μ_{j}^{n} + μ_{j}^{n} \frac{\sum_{i} l_{i j} [d_{i} e^{- 〈 l_{i}, μ 〉} - Y_{i}] - β \frac{\partial}{\partial μ_{j}} U (μ^{n})}{\sum_{i} Y_{i} l_{i j} + β μ_{j}^{n} \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ^{n})}$
wobei $〈 l_{i}, μ 〉 = \sum_{j} l_{i j} μ_{j}$
ist, d_i die erwartete Photonenanzahl ist, die aus der Strahlungsquelle entlang der Projektion i austritt, β der Optimierungsparameter ist, Yi die gemessene Photonenanzahl in der Ebene i ist, n der Iterationsparameter ist und lij die Schnittlänge des Strahls durch das Voxel j in der Ebene i ist.
Der Anfangswert µ_j ⁰ kann mit einem einheitlichen Wert für alle Voxel initialisiert werden. Beispielsweise kann der Anfangswert auf den durchschnittlichen Dämpfungswert bzw. Schwächungskoeffizienten des Brustgewebes festgelegt werden. Eine schnellere Konvergenz kann erreicht werden, falls der Anfangswert mit einem Wert festgelegt wird, der mittels eines normalisierten Rückprojektionsverfahrens bestimmt wurde.
Im Schritt S11, der auf den Schritt S10 folgt, wird geprüft, ob sich µ_j bei der letzen Iteration um eine Differenz verändert hat, die kleiner als ein Schwellenwert ist. Falls die Veränderung von µ_j nicht kleiner als ein Schwellenwert ist, wird der Iterationsparameter n um 1 inkrementiert und zum Schritt S3 zurückgekehrt und eine weitere Iteration durchgeführt. Falls die Änderung von µ_j bei der letzten Iteration kleiner als ein vorbestimmter Schwellenwert ist, wird das Verfahren beendet und µ_j ist ausreichend genau bestimmt.
4 zeigt Details der Gewebedichte-Ermittlungseinrichtung 20. Über einen Eingangsanschluss 32 erhält die Gewebedichte-Ermittlungseinrichtung 20 Daten von der Röntgenstrahlröhre 2, beispielsweise die Röhrenladung, die erwartete Photonenanzahl di, die Intensität des Röntgenstrahls 4 beim Austritt aus der Röntgenstrahlröhre 2 etc. Über einen zweiten Anschluss 34 erhält die Gewebedichte-Ermittlungseinrichtung 20 Daten vom Detektor 18, beispielsweise die empfangene Intensität des Röntgenstrahls 4 an den jeweiligen matrixartig angeordneten Messelementen 14, 16.
Die Gewebedichte-Ermittlungseinrichtung 20 umfasst eine Schwächungskoeffizient-Speichereinrichtung 36, die den zuvor beschriebenen Schritt S7 ausführt und die Schwächungskoeffizienten µ_j und µ_k speichert. Eine Differenzbildungseinrichtung 38 der Gewebedichte-Ermittlungseinrichtung 20 führt den zuvor beschriebenen Schritt S8 aus. Die Gewebedichte-Ermittlungseinrichtung 20 umfasst ferner eine Rauschwert-Schätzeinrichtung 40, die den zuvor beschriebenen Schritt S1 ausführt. Der Schritt S2 wird von einer Gewebedämpfung-Mittelwert-Schätzeinrichtung 42 der Gewebedichte-Ermittlungseinrichtung 20 ausgeführt. Eine Prior-Funktion-Bestimmungseinrichtung 44 der Gewebedichte-Ermittlungseinrichtung 20 ist dazu ausgebildet, die Schritte S3, S4, S5 und S9 auszuführen. Eine Iterationseinrichtung 46 der Gewebedichte-Ermittlungseinrichtung 20 führt die Schritte S6, S10 und S11 aus. Die ermittelten Schwächungskoeffizienten µ_j werden über einen dritten Anschluss 48 der Gewebedichte-Ermittlungseinrichtung 20 an eine Anzeigeeinrichtung 22 ausgegeben.
Es versteht sich, dass die Gewebedichte-Ermittlungseinrichtung 20 aus diskreten Bauteilen aufgebaut sein kann. Die Gewebedichte-Ermittlungseinrichtung 20 kann auch durch einen Computer implementiert sein, der dazu ausgebildet ist, das unter Bezug auf 3 beschriebene Verfahren auszuführen, das in einem Computerprogramm implementiert ist.

Das erfindungsgemäße Verfahren reduziert das Rauschen und verbessert den Signal/Rausch-Abstand sowie den Kontrast/Rausch-Abstand.

Verfahren	Rauschpegel an einer merkmalslosen Stelle - Standardabweichung	Signal/Rausch-Abstand	Kontrast/Rausch-Abstand
Gefilterte Rückprojektion	1,59 %	59,35	2,27
Maximum-Likelihood-Verfahren ohne Prior-Funktionen	1,35 %	61,41	2,04
Maximum-A-Posteriori-Verfahren mit Geman-Prior-Funktion	0,63 %	171,46	4,80

Aufgrund des verbesserten Signal/Rausch-Abstandes und des verbesserten Kontrast/Rausch-Abstandes können niedrigere Strahlungsdosen verwendet werden, was die Nebenwirkungen der diagnostischen Bildgebung reduziert.
5 zeigt die erste Ableitung der Geman-Prior-Funktion für 0 < σ ≤ 0,0006 bei den höchsten Krümmungspunkten. 6 zeigt die zweite Ableitung der Geman-Prior-Funktion, die durchwegs positiv definit ist.
7 bis 9 zeigen einen Schnitt durch ein Brustvolumen. 7 wurde mittels einer gefilterten Rückprojektions-Rekonstruktion ermittelt, die einen hohen Kontrast und eine hohe Auflösung aufweist. Jedoch unterscheidet sich die Helligkeit der Voxel, die dem äußeren Fettgewebe entsprechen, nicht von den Voxeln, die dem dichteren Gewebe im Mittelpunkt der Brust entsprechen. 8 ist eine iterative Maximum-Likelihood-Gradienten-Rekonstruktion ohne Prior-Funktionen und Rauschreduzierung. Sie ergibt eine etwas realistischere Dichtedarstellung des Brustgewebes. 9 ist eine iterative Maximum-A-Posteriori-Rekonstruktion mit einer Geman-Prior-Funktion.
10 bis 12 zeigen das dünnste sichtbare Gefäß in einem Zentralschnitt mit einem merkmalsfreien Hintergrund. 10 wurde mit einer gefilterten Rückprojektion erzeugt. 11 wurde mittels fünf Iterationen einer Maximum-Likelihood-Rekonstruktion ohne Prior-Funktionen erzeugt. 12 zeigt eine Maximum-A-Posteriori-Rekonstruktion mit einer Geman-Prior-Funktion nach der fünften Iteration.
Die Erfindung hat den Vorteil, dass das Rauschen reduziert wird und der Kontrast erhöht wird. Dadurch kann dichteres Gewebe beispielsweise Drüsengewebe, von Fettgewebe unterschieden werden. Darüber hinaus kann krebshaltiges Gewebe besser von anderem Gewebe, beispielsweise Fettgewebe, Drüsengewebe, Gefäßen etc., unterschieden werden.
Bezugszeichenliste

M: Drehpunkt
2: Emitter, Röntgenstrahlröhre
4: Strahlung, Röntgenstrahl
6: Symmetrieachse
8: erste Kompressionsplatte
10: Gewebe, Organ
12: zweite Kompressionsplatte
18: Detektor
14: Messelement
16: Messelement
20: Gewebedichte-Ermittlungseinrichtung
22: Anzeigeeinrichtung
24: Tastatur
30: Voxel
32: erster Anschluss
34: zweiter Anschluss
36: Schwächungskoeffizient-Speichereinrichtung
38: Differenzbildungseinrichtung
40: Rauschwert-Schätzungseinrichtung
42: Gewebedämpfungs-Mittelwert-Schätzeinrichtung
44: Prior-Funktion-Bestimmungseinrichtung
46: Iterationseinrichtung
48: dritter Anschluss

Claims

Verfahren zum Ermitteln eines Schätzwertes eines Schwächungskoeffizienten µ_j eines Voxels j, das einem Abschnitt eines Gewebes (10) zugeordnet ist, wobei der Schwächungskoeffizient µ die Schwächungsrate einer Strahlung (4) in dem Gewebe (10) angibt, wobei die Strahlung (4) von einem Emitter (2) abgegeben wird und von einem Detektor (18) erfasst wird und wobei das Verfahren folgende Schritte aufweist: - Schätzen des Rauschschätzwertes Δµ̂ des Schwächungskoeffizienten µ eines Voxels (S1); - Schätzen eines mittleren Schwächungskoeffizienten µ̅ (S2) ; - Berechnen der Paramater einer Prior-Funktion U(p) aus dem Rauschschätzwert Δµ̃ (S9); und - Iteratives Berechnen des Schätzwertes eines Schwächungskoeffizienten µ_j eines Voxels j aus mittels des Detektors (18) erfassten Projektionsdaten, aus dem Rauschschätzwert Δµ̃, aus dem mittleren Schwächungskoeffizienten µ̅ und aus der Prior-Funktion U(p) (S10, S11) mittels folgender Formel (S10): $μ_{j}^{n + 1} = μ_{j}^{n} + μ_{j}^{n} \frac{\sum_{i} l_{i j} [d_{i} e^{- 〈 l_{i}, μ 〉} - Y_{i}] - β \frac{\partial}{\partial μ_{j}} U (μ^{n})}{\sum_{i} Y_{i} l_{i j} + β μ_{j}^{n} \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ^{n})};$
wobei β ein Optimierungsparameter ist, di die erwartete Photonenzahl ist, die aus der Strahlenquelle entlang der Projektion i austritt, $〈 l_{i}, μ 〉 = \sum_{i} l_{i j} μ_{j}$
die gemessenen Photonenanzahl in der Projektion i ist und l_ij die Schnittlänge des Strahls durch das Voxel j in der Projektion i ist und wobei β mittels folgender Formel bestimmt wird (S9): $β = \frac{Δ \tilde{μ} \cdot \bar{\sum_{i} Y_{i} l_{i j}}}{\bar{μ} \cdot \frac{\partial}{\partial μ_{j}} U (μ) {|_{Δ μ_{j} = Δ \tilde{μ}} - Δ \tilde{μ} \cdot \bar{μ} \cdot \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ) |}_{^{^{^{_{^{^{Δ μ_{j} = Δ \tilde{μ}}}}}}}}};$
wobei $\bar{\sum_{i} Y_{i} l_{i j}}$
der Mittelwert über alle Voxel des Volumenelementes über alle i und j ist.
Verfahren zum Ermitteln eines Schätzwertes eines Schwächungskoeffizienten µ_j eines Voxels j, das einem Abschnitt eines Gewebes (10) zugeordnet ist, wobei der Schwächungskoeffizient µ die Schwächungsrate einer Strahlung (4) in dem Gewebe (10) angibt, wobei die Strahlung (4) von einem Emitter (2) abgegeben wird und von einem Detektor (18) erfasst wird und wobei das Verfahren folgende Schritte aufweist: - Schätzen des Rauschschätzwertes Δµ̃ des Schwächungskoeffizienten µ eines Voxels (S1); - Schätzen eines mittleren Schwächungskoeffizienten µ̅ (S2) ; - Berechnen der Paramater einer Prior-Funktion U(µ) aus dem Rauschschätzwert Δµ̃ (S9); und - Iteratives Berechnen des Schätzwertes eines Schwächungskoeffizienten µ_j eines Voxels j aus mittels des Detektors (18) erfassten Projektionsdaten, aus dem Rauschschätzwert Δµ̃, aus dem mittleren Schwächungskoeffizienten µ̅ und aus der Prior-Funktion U(p) (S10, S11) mittels folgender Formel: $μ_{j}^{n + 1} = μ_{j}^{n} + μ_{j}^{n} \frac{\sum_{i} l_{i j} [d_{i} e^{- 〈 l_{i}, μ 〉} - Y_{i}] - β \frac{\partial}{\partial μ_{j}} U (μ^{n})}{\sum_{i} l_{i j} 〈 l_{i}, μ^{n} 〉 d_{i} e^{-}^{〈 l_{i}, μ^{n} 〉} + β μ_{j}^{n} \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ^{n})};$
wobei β ein Optimierungsparameter ist, $〈 l_{i}, μ 〉 = \sum_{j} l_{i j} μ_{j}$
ist, di die erwartete Photonenzahl ist, die aus der Strahlenquelle entlang der Projektion i austritt, Y_i die gemessenen Photonenanzahl in der Projektion i ist und l_ij die Schnittlänge des Strahls durch das Voxel j in der Projektion i ist und wobei β mittels folgender Formel bestimmt wird: $β = \frac{Δ \tilde{μ} \cdot \sum_{i} l_{i j} 〈 l_{i}, \bar{μ} 〉 d_{i} e^{- 〈 l_{i}, μ^{n} 〉}}{\bar{μ} \cdot \frac{\partial}{\partial μ_{j}} U (μ) {|_{Δ μ_{j} = Δ \tilde{μ}} - Δ \tilde{μ} \cdot \bar{μ} \cdot \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ) |}_{^{^{^{_{^{^{Δ μ_{j} = Δ \tilde{μ}}}}}}}}};$
wobei $〈 l_{i}, \bar{μ} 〉 = \sum_{j} l_{i j} \bar{μ}$
ist.
Verfahren nach Anspruch 1 oder 2, dadurch gekennzeichnet, dass - der Schritt des Schätzens des Rauschschätzwertes Δµ̃ des Schwächungskoeffizienten µ eines Voxels (S1)die Gewebedicke und eine Information über die abgegebene Strahlung (4) berücksichtigt; und/oder - der Schritt des Schätzens eines mittleren Schwächungskoeffizienten µ̅ (S2) auf Grundlage einer Gewebedichteschätzung und der Energie der vom Emitter (2) abgegebenen Strahlung (4) durchgeführt wird.
Verfahren nach einem der Ansprüche 1 bis 3, dadurch gekennzeichnet, dass der Schritt des Berechnens der Parameter der Prior-Funktion U(p) die Parameter einer Geman-Prior-Funktion berechnet, die folgende Gleichung aufweist: $U (μ) = \sum_{j} \sum_{k} \frac{Δ μ_{j k}^{2} σ^{m}}{2 {(\sqrt{Δ μ_{j k}^{2} / 2 + σ^{2} / 2})}^{m}};$
wobei σ ein Optimierungsparameter ist, der aus dem Rauschschätzwert Δµ̃ bestimmt wird, und m eine Konstante ist und wobei Δµ_jk die Differenz der Schwächungskoeffizienten benachbarter Voxel ist, die je einem Gewebeabschnitt zugeordnet sind, und j und k lineare Voxel-Indizes über drei Dimensionen sind.
Verfahren nach einem der Ansprüche 1 bis 4, dadurch gekennzeichnet, dass, falls die Prior-Funktion U(µ) eine konvexe Form aufweisen soll, σ so gewählt wird, dass die erste Ableitung der Prior-Funktion U(µ) ihre stärkste Krümmung bei Werten von Δµ erreicht, die dem Rauschschätzwert Δµ̃ entsprechen (S4) und, falls die Prior-Funktion U(µ) eine nicht-konvexe Form aufweisen soll, σ so gewählt wird, dass die erste Ableitung der Prior-Funktion U(µ) ihr Maximum bei Werten von Δµ erreicht, die dem Rauschschätzwert Δµ̃ entsprechen (S5).
Gewebedichte-Ermittlungseinrichtung (20), mit - einem Emitter (2), der eine Strahlung (4) in ein Gewebe (10) abgibt; - einem Detektor (18), der die Strahlung (4) erfasst, nachdem sie das Gewebe (10) durchlaufen hat; - einer Rauschwert-Schätzeinrichtung (40), die zum Schätzen des Rauschschätzwertes Δµ̃ des Schwächungskoeffizienten µ eines Voxels ausgebildet ist; - einer Gewebedämpfung-Mittelwert-Schätzeinrichtung (42), die zum Schätzen eines mittleren Schwächungskoeffizienten µ̅ ausgebildet ist; - einer Prior-Funktion-Bestimmungseinrichtung (44), die dazu ausgebildet ist, die Parameter einer Prior-Funktion U(p) aus dem Rauschschätzwert Δµ̃ zu berechnen; und - einer Iterationseinrichtung (46), die dazu ausgebildet ist, den Schätzwert eines Schwächungskoeffizienten µ_j eines Voxels j aus mittels des Detektors (18) erfassten Projektionsdaten, aus dem erwarteten Rauschwert Δµ̃, aus dem mittleren Gewebedämpfungsschätzwert µ̅ und aus der Prior-Funktion U(µ) mittels folgender Formel zu bestimmen: $μ_{j}^{n + 1} = μ_{j}^{n} + μ_{j}^{n} \frac{\sum_{i} l_{i j} [d_{i} e^{- 〈 l_{i}, μ 〉} - Y_{i}] - β \frac{\partial}{\partial μ_{j}} U (μ^{n})}{\sum_{i} Y_{i} l_{i j} + β μ_{j}^{n} \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ^{n})};$
wobei β ein Optimierungsparameter ist, di die erwartete Photonenzahl ist, die aus der Strahlenquelle entlang der Projektion i austritt, $〈 l_{i}, μ 〉 = \sum_{j} l_{i j} μ_{j} ist,^{\bar{\sum_{i} Y_{i} l_{i j}}}$
der Mittelwert über alle Voxel des Volumenelementes über alle i und j ist, Yi die gemessenen Photonenanzahl in der Projektion i ist und l_ij die Schnittlänge des Strahls durch das Voxel j in der Projektion i ist und wobei die Iterationseinrichtung dazu ausgebildet ist, β mittels folgender Formel zu bestimmen: $β = \frac{Δ \tilde{μ} \cdot \bar{\sum_{i} Y_{i} l_{i j}}}{\bar{μ} \cdot \frac{\partial}{\partial μ_{j}} U (μ) {|_{Δ μ_{j} = Δ \tilde{μ}} - Δ \tilde{μ} \cdot \bar{μ} \cdot \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ) |}_{^{^{^{_{^{^{Δ μ_{j} = Δ \tilde{μ}}}}}}}}};$
wobei $\bar{\sum_{i} Y_{i} l_{i j}}$
der Mittelwert über alle Voxel des Volumenelementes über alle i und j ist.
Gewebedichte-Ermittlungseinrichtung (20), mit - einem Emitter (2), der eine Strahlung (4) in ein Gewebe (10) abgibt; - einem Detektor (18), der die Strahlung (4) erfasst, nachdem sie das Gewebe (10) durchlaufen hat; - einer Rauschwert-Schätzeinrichtung (40), die zum Schätzen des Rauschschätzwertes Δµ̃ des Schwächungskoeffizienten µ eines Voxels ausgebildet ist; - einer Gewebedämpfung-Mittelwert-Schätzeinrichtung (42), die zum Schätzen eines mittleren Schwächungskoeffizienten µ̅ ausgebildet ist; - einer Prior-Funktion-Bestimmungseinrichtung (44), die dazu ausgebildet ist, die Parameter einer Prior-Funktion U(p) aus dem Rauschschätzwert Δµ̃ zu berechnen; und - einer Iterationseinrichtung (46), die dazu ausgebildet ist, den Schätzwert eines Schwächungskoeffizienten µ_j eines Voxels j aus mittels des Detektors (18) erfassten Projektionsdaten, aus dem erwarteten Rauschwert Δµ̃, aus dem mittleren Gewebedämpfungsschätzwert µ̅ und aus der Prior-Funktion U(µ) mittels folgender Formel zu bestimmen: $μ_{j}^{n + 1} = μ_{j}^{n} + μ_{j}^{n} \frac{\sum_{i} l_{i j} [d_{i} e^{- 〈 l_{i}, μ 〉} - Y_{i}] - β \frac{\partial}{\partial μ_{j}} U (μ^{n})}{\sum_{i} l_{i j} 〈 l_{i}, μ^{n} 〉 d_{i} e^{- 〈 l_{i}, μ^{n} 〉} + β μ_{j}^{n} \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ^{n})};$
wobei β ein Optimierungsparameter ist, di die erwartete Photonenzahl ist, die aus der Strahlenquelle entlang der Projektion i austritt, $〈 l_{i}, μ 〉 = \sum_{j} l_{i j} μ_{j}$
ist, Y_i die gemessenen Photonenanzahl in der Projektion i ist und l_ij die Schnittlänge des Strahls durch das Voxel j in der Projektion i ist und wobei die Iterationseinrichtung dazu ausgebildet ist, β mittels folgender Formel zu bestimmen: $β = \frac{Δ \tilde{μ} \cdot \sum_{i} l_{i j} 〈 l_{i}, \bar{μ} 〉 d_{i} e^{- 〈 l_{i}, μ^{n} 〉}}{\bar{μ} \cdot \frac{\partial}{\partial μ_{j}} U (μ) {|_{Δ μ_{j} = Δ \tilde{μ}} - Δ \tilde{μ} \cdot \bar{μ} \cdot \frac{\partial^{2}}{\partial μ_{j}^{2}} U (μ) |}_{^{^{^{_{^{^{Δ μ_{j} = Δ \tilde{μ}}}}}}}}};$
wobei $〈 l_{i}, \bar{μ} 〉 = \sum_{j} l_{i j} \bar{μ}$
ist.
Gewebedichte-Ermittlungseinrichtung (20) nach Anspruch 6 oder 7, dadurch gekennzeichnet, dass die Gewebedämpfung-Mittelwert-Schätzeinrichtung (42) den mittleren Schwächungskoeffizienten µ̅ auf Grundlage einer Gewebedichteschätzung und der Energie der vom Emitter (2) emittierten Strahlung (4)ermittelt; und/oder die Rauschwert-Schätzeinrichtung (40) dazu ausgebildet ist, den erwarteten Rauschwert Δµ̃ des Schwächungskoeffizienten µ eines Voxels aus der Gewebedicke und einer Information über die abgegebene Strahlung (4) zu schätzen.
Gewebedichte-Ermittlungseinrichtung (20) nach einem der Ansprüche 6 bis 8, dadurch gekennzeichnet, dass die Prior-Funktion-Bestimmungseinrichtung (44) dazu ausgebildet ist, die Prior-Funktion U(µ) als eine Geman-Prior-Funktion zu berechnen, die folgende Gleichung aufweist: $U (μ) = \sum_{j} \sum_{k} \frac{Δ μ_{j k}^{2} σ^{m}}{2 {(\sqrt{Δ μ_{j k}^{2} / 2 + σ^{2} / 2})}^{m}};$
wobei σ ein Optimierungsparameter ist, der aus dem Rauschschätzwert Δµ̃ bestimmt wird, und m eine Konstante ist und wobei Δµ_jk die Differenz der Schwächungskoeffizienten benachbarter Voxel ist, die je einem Gewebeabschnitt zugeordnet sind, und j und k lineare Voxel-Indizes über drei Dimensionen sind.
Gewebedichte-Ermittlungseinrichtung (20) nach einem der Ansprüche 6 bis 9, dadurch gekennzeichnet, dass die Prior-Funktion-Bestimmungseinrichtung (44) dazu ausgebildet ist, dass, falls die Prior-Funktion U(µ) eine konvexe Form aufweisen soll, σ so zu wählen, dass die erste Ableitung der Prior-Funktion U(µ) ihre stärkste Krümmung bei Werten von Δµ erreicht, die dem Rauschschätzwert Δµ̃ entsprechen und, falls die Prior-Funktion U(µ) eine nicht-konvexe Form aufweisen soll, σ so zu wählen, dass die erste Ableitung der Prior-Funktion U(µ) ihr Maximum bei Werten von Δµ erreicht, die dem Rauschschätzwert Δµ̃ entsprechen.
Computerprogramm geladen in einen Speicher eines Computers, umfassend Mittel, die zur Ausführung aller Schritte eines Verfahrens nach einem der Ansprüche 1 bis 5 eingerichtet sind, wenn das Copmuterprogramm auf dem Computer ausgeführt wird.
Medizinisches System, mit - der Gewebedichte-Ermittlungseinrichtung (20) nach einem der Ansprüche 6 bis 10; - Positionierungsmitteln (8, 12), um zumindest einen Teil eines Organs (10) eines Patienten räumlich gegenüber dem Emitter (2) und dem Detektor (18) festzulegen; und - einer Anzeigeeinrichtung (22), die dazu ausgebildet ist, Repräsentationen der Gewebeschwächungskoeffizienten µ_j anzuzeigen.