-
Die
Erfindung betrifft ein Verfahren zur Korrektur von Streustrahlungsfehlern
bei der Radiographie und der Computertomographie mit Flächendetektoren.
-
Aus
der Druckschrift [ZSR2005] ist bekannt, dass die Technologie der
Flachbilddetektoren (= FPD = flat Panel detectors) eine Niedrigkontrastauflösung bei
der 3-dimensionalen Rekonstruktion der Gewebedichteverteilung mit
Hilfe der C-Bogen-Computertomographie
ermöglicht.
Wegen des großen
durchstrahlten Patientenvolumens bei jedem Projektionsbild stellt
die nur sehr wenig Bildinformation enthaltende Streustrahlung hierbei
ein schwerwiegendes Problem dar. Die Auswirkungen der Streustrahlung
sind: Kontrastverlust, Verzerrungen der quantitativen Ergebnisse,
wie Cupping und Artefakte, und Erhöhung des Rauschens. Die Verwendung
von Streustrahlenrastern ist im Allgemeinen nicht ausreichend. Daher
sind zusätzliche
Korrekturverfahren erforderlich, die auf der Basis von Messungen
oder rechnerischen Schätzungen
der Streustrahlung deren negative Auswirkungen auf die Bildqualität kompensieren
sollen.
-
Mit
solchen Korrekturverfahren werden im Idealfall die quantitative
Genauigkeit Hounsfield-Units (=HU) und die ursprünglichen Kontraste zwischen
Weichteilgeweben wieder hergestellt. Zwar bleibt das Qualitätsmaß Kontrast-zu-Rauschen
erhalten, jedoch geht in der Regel mit den Korrekturen eine Zunahme
des Bildrauschens einher, so dass das Rauschen im gleichen Maße wie der
Kontrast zunimmt. Dieser Effekt ist umso stärker, je größer der relative Streustrahlungsanteil
S/P (= Streustrahlungs-zu-Primärstrahlungs-Intensität) in den
Daten ist. Dieser wirkt sich auch in den rekonstruierten Bildern
durch starke Inhomogenitäten
des Bildrauschens aus. Diese Inhomogenitäten äußern sich sowohl durch erhöhte Rauschamplituden
(=Standardabweichung des Pixelrauschens) als auch durch ihre langgezogene
Form, also unregelmäßige feine
Streifen, und Ausrichtung besonders in Richtungen, in denen die
Röntgenstrahlen
einer starken Schwächung
unterliegen, zum Beispiel Längsachsen
bei elliptischen Querschnitten, Weichteilgewebe zwischen Knochen,
und wo in der Regel auch der relative Streustrahlungsanteil S/P
relativ groß ist.
Auf diese Inhomogenität
des Bildrauschens durch Streustrahlungskorrekturen wurde bereits
in der Druckschrift [JS1982] hingewiesen.
-
Zur
Korrektur dieses Streustrahlungsproblems wurde beispielsweise in
der Druckschrift [SDB2006] eine einfache Messmethode zur Schätzung der
Streustrahlungsintensität
innerhalb des eigentlichen Messfeldes durch Interpolation von Messungen
der Streustrahlung im Kollimatorschatten außerhalb des Messfeldes vorgeschlagen.
Da aber innerhalb des Messfeldes die Streustrahlung nicht gemessen
wird, ist diese Methode nur sehr grob. Die Problematik des Bildrauschens
im Zusammenhang mit der Korrektur wird in dieser Schrift jedoch
nicht diskutiert. Die Thematik reiner Rechenverfahren zur Streustrahlungskorrektur
in der CT mit Flächendetektoren
(Cone-beam-CT) steht bisher in der Fachliteratur noch in den Anfängen. In
der Druckschrift [ZSR2005] wird ein recht allgemeiner rechnerischer
Korrekturalgorithmus beschrieben und es werden auch andere verwandte
Ansätze
zitiert. Zwar werden die wesentlichen Bildqualitätsaspekte angesprochen, jedoch fehlt
eine quantitative Bildqualitäts-Analyse.
-
Ein
Problem dieser bekannten Streustrahlungskorrekturverfahren liegt
darin, dass durch die Streustrahlungskorrektur sich das „Signal-zu-Rausch" – Verhältnis beziehungsweise das Verhältnis von
Streustrahlungsintensität
zu Primärstrahlungsintensität durch
die Streustrahlungskorrektur verschlechtert.
-
Es
ist daher Aufgabe der Erfindung, ein verbessertes Streustrahlkorrekturverfahren
zu finden, welches das „Signal-zu-Rausch" – Verhältnis möglichst wenig bis nicht verschlechtert.
-
Diese
Aufgabe wird durch die Merkmale der unabhängigen Patentansprüche gelöst. Vorteilhafte
Weiterbildungen der Erfindung sind Gegenstand untergeordneter Ansprüche.
-
Im
Folgenden wird die Erkenntnis des Erfinders beschrieben und der
mathematische und sachliche Zusammenhang von Bildqualitätsbegriffen
erläutert.
Hierbei werden die folgenden Bezeichnungen verwendet:
(y, z)
Koordinaten auf dem Detektor;
Io =
Io (y, z) gemessene ungeschwächte Intensitätsverteilung
ohne Objekt im Strahlengang;
It = It (y, z) gemessene totale Intensitätsverteilung
mit Objekt (einschließlich
Streustrahlung);
Ip = Ip (y,
z) Primärintensitätsverteilung
mit Objekt im Strahlengang (ohne Streustrahlung);
Is = Is (y, z) Intensitätsverteilung
nur der vom Objekt erzeugten Streustrahlung.
-
Im
Folgenden wird vorwiegend mit normierten Intensitäten gerechnet.
Normierung bedeutet: Division durch Io.
Für die
normierten Intensitäten
werden die folgenden Bezeichnungen eingeführt :
P = Ip/Io normierte Primärstrahlungsverteilung (ohne
Streustrahlung);
PS = It/Io normierte totale Intensitätsverteilung
(einschließlich
Streustrahlung);
S = Is/Io normierte
Streustrahlungsverteilung.
-
In
der CT werden die Projektionsbilddaten logarithmiert. Diese logarithmierten
Werte werden "CT-Projektionswerte" genannt. Hierfür werden
kleine Buchstaben verwendet. Es seien: fehlerfreie CT-Projektionswerte: p = – ln(P) (#1a);die durch
Streustrahlung gestörten
CT-Projektionswerte sind PS = – ln(PS) = – ln
(P + S) = p – ln(1
+ S/P) (#1b); der
durch Streustrahlung bedingte Störungsanteil
ist: δpS = – ln(1 +
S/P) (#1c).
-
Nachfolgend
werden die Begriffe Kontrast-, Rausch-Übertragung und differentielles
CNR bei Streustrahlung erklärt:
Für die folgenden
Betrachtungen benötigen
wir die Taylor-Formel
beziehungsweise Kettenregel der Differenziation und die Gaußsche Fehlerfortpflanzungsregel.
Für die „deterministische" Fortpflanzung von
differenziellen Änderungen ΔX von Erwartungswerten
X über
eine differenzierbare funktionale
f gilt
auf Grund der Taylor-Formel:
-
Die
Fehlerfortpflanzungsregel von Gauß für Funktionen f von unkorrelierten
Zufallsvariablen X
1, X
2,
... lautet:
-
Nachfolgend
wird der Begriff des „differentiellen
Kontrastes" beschrieben.
Hierzu wird eine kleine, infinitesimale Signaldifferenz ΔP angenommen
und gefragt, wie diese durch die Logarithmierung in Gleichung (#1a,
b) übersetzt
wird. Hierfür
wird ΔP
beziehungsweise Δp
als "differenzieller
Kontrast" bei normierten
Intensitätswerten
beziehungsweise logarithmierten Projektionswerten bezeichnet. Die
Differenzierbarkeit von P vorausgesetzt, kann man ΔP auch interpretieren
als
mit einer kleinen örtlichen
Verschiebung Δξ. Mit Streustrahlung
gilt:
ΔPS = Δ (P + S) ≈ ΔP (#4a). -
Dabei
wurde vorausgesetzt, dass die Streustrahlungsverteilung S(y, z)
im Erwartungswert räumlich sehr
niederfrequent ist, das heißt,
dass S(y, z) sich im Vergleich zur Primärstrahlungsverteilung P(y,
z) nur vernachlässigbar ändert, also
|ΔS/ΔP| ≈ 0 gilt.
-
Aus
Gleichung (#2) folgt speziell für
f = –ln
und X = P:
-
Die
Signaldifferenz bei den logarithmierten Daten wird durch einen Streustrahlungshintergrund
S wie folgt beeinflußt:
-
Dabei
wird |ΔS/ΔP| ≈ 0 vorausgesetzt.
Die Gleichung (#4c) besagt für
logarithmische CT-Projektionsdaten, dass die Streustrahlung eine
Reduzierung des differenziellen Kontrasts um den Faktor
bewirkt.
-
Unter
der Annahme der Poisson-Verteilung als Wahrscheinlichkeitsverteilung
ist σ(I)
= √I.
-
Mit
der Gauß-Fehlerfortpflanzungsregel
(#3) erhält
man
wobei
I << I
o sowie σ(I) = √
I und
berücksichtigt wurde. D.h. es gilt:
und
σ(PS) = σ(P
+ S) = σ(P) √1 + S/P (#5b). -
Somit
wird durch die Streustrahlung das Rauschen der nichtlogarithmierten
normierten Intensitäten um
den Faktor √1 + S/P erhöht.
-
Als
sinnvolles Qualitätsmaß wird das
differenzielle Kontrast-zu-Rausch-Verhältnis dCNR
eingeführt, welches
einen Quotienten aus einer differenziellen Signaländerung
und der Standardabweichung des Signalrauschens darstellt.
-
Per
definitionem und mit Gleichung (#5b) gilt:
-
Allgemein
gilt somit, dass das dCNR betragsmäßig invariant unter differenzierbaren
(monotonen) Signaltransformationen f ist, mit: dCNR(f(P)) = dCNR(P) (#7)
-
Der
Beweis ergibt sich durch die Beziehungen (#2) und (#3): Bei der
Quotientenbildung kürzen
sich die Ableitungen betragsmäßig weg,
da wegen der Voraussetzung der Monotonie von f die Ableitung nirgendwo Null
wird. Dies trifft insbesondere für
die Logarithmierung f = – ln
zu, wobei sich bei dCNR nur das Vorzeichen umkehrt. Also gilt unabhängig, ob
die Projektionsdaten logarithmiert sind oder nicht, dass die Streustrahlung eine
Reduzierung des differentiellen CNR um den Faktor
bewirkt.
-
Bei
den bisherigen Überlegungen
wurde vorausgesetzt, dass die normierte Primärintensitätsverteilung P(y, z) und die
Streuintensitätsverteilung
S(y, z) bekannt seien. In Wirklichkeit ist durch den Detektor nur
deren Überlagerung,
die „totale
Intensitätsverteilung" PS(y,
z) = P(y, z) + S(y, z), gegeben. Die Streustrahlverteilung S(y,
z) stellt dabei eine Störung
dar. Gesucht ist jedoch die ungestörte Primärstrahlverteilung P(y, z).
-
Da
die Streustrahlverteilung nicht exakt bekannt ist, muss sie durch
eine Schätzung
korrigiert werden. Die Schätzung
kann mittels Messungen, zum Beispiel durch die Beamstop-Methode gemäß [NTC2002]
oder durch Messung im Kollimatorschatten und Interpolation gemäß [SDB2006])
oder durch rechnerische Methoden gemäß [ZSR2005] erfolgen. Für die geschätzte und
für die
Korrektur verwendete Streuverteilung wird hier die Bezeichnung Scor(y, z) verwendet, für die mittels dieser Korrektur
geschätzte
Primärintensitätsverteilung schreiben
wir Pcor(y, z). Damit gilt Pcor(y, Z) = PS(y, z) – Scor(y, Z) (#9a) Pcor(y, z) = P(y, z) + εcor(y,
z) (#9b)mit
dem Streustrahlungskorrekturfehler εcor(y,
z) = S(y, z) – Scor(y, z) (#9c).
-
Unter
der Annahme, dass der Korrekturfehler relativ zu S und P klein ist,
wird im Folgenden die Auswirkung der Streustrahlungskorrektur auf
den differenziellen Kontrast, das Rauschen und das dCNR untersucht.
Dazu werden auch die logarithmierten Projektionsdaten benötigt: pcor = – ln Pcor = – ln(PS – Scor) (#10a) pcor = – ln (P
+ εcor) (#10b)
-
Es
wird nun der differenzielle Kontrast durch Streustrahlungskorrektur
wieder hergestellt. Die Störungsrechnung,
angewendet auf Gleichung (#10b), liefert:
wobei
ausgenutzt wurde, dass ε
cor sehr glatt verläuft und die Ableitung lokal
~ 0 gesetzt werden kann. Unter der Voraussetzung
eines vernachlässigbaren
Korrekturfehlers |ε
cor| << |P| ergibt sich
somit, dass sich durch eine genaue Streustrahlungskorrektur der
ursprüngliche
differenzielle Kontrast wiedergewinnen lässt, d.h. der streustrahlungsbedingte
differentielle Kontrastverlust wird durch die Korrektur durch eine
differentielle Kontrastverstärkung
um den Faktor (1 + S/P) kompensiert.
-
Die
dCNR-Erhaltung mit Streustrahlungskorrektur folgt aus der Invarianzbeziehung
Gleichung (# 7), mit: dCNR(pcor) = dCNR(pS) (#12).
-
In
Worten bedeutet dies, dass das differenzielle Kontrast-zu-Rausch-Verhältnis dCNR
durch eine genaue Streustrahlungskorrektur unverändert bleibt.
-
Trotzdem
bleibt die Erhöhung
des Rauschens durch die Streustrahlungskorrektur. Die Störungsrechnung,
angewendet auf Gleichung (#10a), liefert nämlich:
-
Unter
der Voraussetzung eines vernachlässigbaren
Korrekturfehlers |εcor| << |P| gilt demnach,
dass durch eine genaue Streustrahlungskorrektur das Rauschen der
durch Streustrahlung gestörten
logarithmierten Projektionsdaten um den Faktor (1 + S/P) erhöht wird.
-
Dieses
Phänomen
führt zu
einer sehr inhomogenen Verstärkung
des Bildrauschens, vor allem in Richtungen mit sehr kleinen Werten
von P, d.h. starker Strahlungsschwächung. Beispielsweise wird
das Rauschen bei einem Wert von S/P = 1, ein auch trotz Streustrahlenraster
nicht selten vorkommender Fall, auf das Doppelte erhöht.
-
Eine „Standard-Streustrahlungskorrektur" – basierend auf Messverfahren
oder auf Rechenmodellen –, stützt sich
im Wesentlichen auf die Gleichung (#10a) und umfasst die folgenden
Schritte:
- – Schätzung der
Streustrahlungsverteilung durch Messung oder Rechnung Scor(y,
z)
- – Subtraktion
von der gemessenen „totalen" Verteilung PS(y, z) und
- – Logarithmierung.
-
Hierbei
erhöht
sich das Rauschen gemäß Gleichung
(#13), auch wenn das dCNR sich nicht verschlechtert.
-
Erfindungsgemäß lässt sich
dieses Problem wie folgt umgehen: Die Gleichung (#10a) kann durch
Faktorisierung der Korrektur umgeschrieben werden wie folgt: Pcor = PS – Scor
Pcor = PS (1 – Scor/PS) (#14).
-
Logarithmierung
ergibt dann: pcor = – ln Pcor = – ln(PS) – ln(1 – Scor/PS)
pcor = pS – ln(1 – Scor/PS)
pcor = pS + δpcor (#15).
-
Dabei
wird der Standard-Streustrahlungskorrekturterm allgemein durch δpcor = pcor – pS (#16)definiert.
Gleichung (#16) gilt allgemein, nicht nur für subtraktive Korrekturen wie
in Gleichung (#15), sondern auch für multiplikative, wie z.B.
in der Druckschrift [ZSR2005] erwähnt.
-
Die
Funktion Scor(y, z) in Gleichung (#14) ist
sehr glatt und praktisch ohne Rauschen. Hingegen ist der Standardkorrekturterm δpcor(y, z) für die logarithmischen Daten
nicht mehr glatt, sondern hoch moduliert wegen der Division durch
PS (y, z) und ist deswegen auch verrauscht.
Das Rauschen ist streng korreliert mit dem von PS.
Tatsächlich
beruht auf der Abhängigkeit
von PS die Möglichkeit der Kontrastrestauration,
leider allerdings mit dem Nebeneffekt der Erhöhung des Rauschens.
-
Die
erfindungsgemäße Idee
besteht nun darin, vor der Korrektur gemäß Gleichung (#15) zuerst allein auf
den Streustrahlungskorrekturterm δpcor eine möglichst kontrasterhaltende
Rauschfilterung anzuwenden und dann erst den rauschgefilterten Korrekturterm
auf die unkorrigierten Daten pS zu addieren.
-
Symbolisiert
man die Rauschfilterung mit einem Operator F, dann wird Gleichung
(#15) ersetzt durch: Pcor,F = pT +
F(δpcor) (#17a).
-
Im
Allgemeinen ist die Rauschfilterung nicht-linear und 2-dimensional, kann
aber auch z.B. nur 1- dimensional in Bildzeilenrichtung gewählt werden.
Beispiele sind 1- dimensionale oder 2- dimensionale Medianfilterung.
Medianfilterung ist das einfachste Beispiel für kantenerhaltende und rauschenreduzierende
Filterung. Im Gegensatz zur „Standard- Streustrahlungskorrektur" umfasst dieses Grundverfahren
die folgenden Schritte:
- – Schätzung der Streustrahlungsverteilung
(durch Messung oder Rechnung) Scor(y, z);
- – Berechnung
des Standardkorrekturterms δpcor gemäß Gleichung
(#15) beziehungsweise (#16);
- – Rauschfilterung
des Standardkorrekturterms: F(δpcor)
- – Addition
des rauschgefilterten Standardkorrekturterms F(δpcor)
zu den (logarithmierten) gemessenen „totalen" Projektionsdaten pS.
-
Dies
ist die einfache Variante der erfindungsgemäßen Strahlungskorrektur. Gemäß einer
Weiterführung
und Verbesserung des Verfahrens kann auch eine adaptive Rauschfilterung
des Korrekturterms durchgeführt
werden, wobei nun eine rauschadaptive und verbessert kontrasterhaltende
Rauschfilterung angewendet wird. Symbolisiert man die adaptive Rauschfilterung
mit einem Operator Φ,
dann wird Gleichung (#15) ersetzt durch: Pcor, Φ =
pS + Φ (δpcor) (#17b).
-
Rauschadaptiv
bedeutet hier, dass das Rauschen umso stärker gefiltert werden soll,
je größer die Standardabweichung
ist, so dass nach der Filterung die Standardabweichung des Rauschens „homogenisiert" ist. Diese Homogenisierung
des Rauschens im Streustrahlungskorrekturterm bei den Projektionsdaten
wird sich ebenfalls in einer Homogenisierung des Rauschens in den
rekonstruierten CT-Bildern bemerkbar machen.
-
Einfache
Beispiele für
adaptive Filter sind in der Standardliteratur über Bildverarbeitung, z.B. [Jäh1997],
zu finden: rauschabhängige
normalisierte Faltung, steuerbare Filter mit von der lokalen Bildstruktur abhängigen Filtermasken;
Medianfilter und allgemeine Glättungsfilter
mit rauschabhängiger
Maskengröße.
-
Um
durch die Filterung die lokale Bildstruktur nicht zu verändern, insbesondere
z.B. Kanten möglichst nicht
zu verschleifen, ist eine vorherige Bildanalyse mittels Kantendetektion
und/oder Segmentierung zweckmäßig. Dadurch
können
unterschiedliche Bildbereiche unterschiedlich, d.h. adaptiv, gefiltert
werden. Zum Beispiel können
Kanten nur in Längs- und nicht in Querrichtung
gefiltert werden. Medianfilter haben bereits die günstige Eigenschaft,
dass sie Kanten weitgehend beibehalten, jedoch ist ihre Rauschminderungswirkung
bei Maskenvergrößerung geringer
als bei linearen Mittelungsfiltern.
-
Erfindungsgemäß ist es
besonders günstig,
die Bildanalyse, also die Kantendetektion und/oder Segmentierung,
im unkorrigierten Bild pS = pS (x,
y) durchzuführen,
da die Bildstrukturen dort im Allgemeinen deutlicher erkennbar sind
als im Korrekturterm.
-
Aus
dem unkorrigierten Projektionsbild, gegebenenfalls in den segmentierten
Bildbereichen, falls eine Segmentierung durchgeführt wurde, sind dann die Standardabweichungen
abzuschätzen,
um davon abhängig adaptiv
die Größe der lokalen
Rauschfiltermaske zur Rauschfilterung in einer Pixelumgebung zu
bestimmen.
-
Aus
den Gleichungen (# 5a-b) und P
S = e
– p S folgt:
Diese
Größe kann
für jedes
Pixel aus dem unkorrigierten Projektionsbild p
S(x,
y) berechnet werden.
-
Das
Bildrauschen wird durch die Streustrahlungskorrektur gemäß Gleichung
(#13) erhöht.
Dies kann in unserer Abschätzung
durch Multiplikation mit dem Faktor (S/P) in Gleichung (#18) noch
berücksichtigt
werden. Eine Schätzung
von (S/P) ergibt sich mit Gleichung (#9a) zu Scor/Pcor = Scor/(PS – Scor)
-
Die
folgende Größe
wird
als „Rauschzuwachsterm" bezeichnet.
-
Wenn
man annimmt, dass das Pixelrauschen statistisch unkorreliert ist,
dann lässt
sich die Standardabweichung durch Mittelung über N·N Pixel proportional zum
Faktor 1/N reduzieren. Auf diese Weise ergibt sich ein Maß zur adaptiven
Bestimmung der Maskengröße von Mittelungsfiltern
zur Rauschreduktion.
-
Die
absoluten Werte der Standardabweichung in Gleichung (#18) sind nicht
nötig,
d.h. auch der ungeschwächte
Intensitätswert
I0 ist nicht nötig. Da nur der Inhomogenitätszuwachs
des Rauschens durch die Streustrahlkorrektur, der vor allem durch
große
Werte des „Rauschzuwachsterms" bestimmt wird, reduziert werden
soll, genügt
es, die adaptive Rauschfilterung auf Bereiche mit relativ großen Werten
des „Rauschzuwachsterms" Gleichung (#19)
zu beschränken.
Angenommen für
mittlere Werte des „Rauschzuwachsterms" werde nur eine Rauschfilterung über eine
normale Maskengröße N0·N0 gemacht,
es treten aber in einem Bildbereich z.B. um einen Faktor 3 größere Werte
des Terms (#19) auf, dann muss die normale Mittelungsmaske in diesem
Bereich auf N·N
mit N = 3·N0
vergrößert werden.
Je nach Fall – z.B.
entlang von Kanten, abhängig von
der Segmentierung – können auch
rechteckige statt quadratische Mittelungsmasken vorkommen.
-
Die
erfindungsgemäße Verarbeitungskette
umfasst demnach die folgenden Schritte:
- – Schätzung der
Streustrahlungsverteilung (durch Messung oder Rechnung) Scor(y, z),
- – Berechnung
des Streustrahlungskorrekturterms δpcor gemäß Gleichung
(#15),
- – Bildanalyse
durch Kantendetektion und/oder Segmentierung des unkorrigierten
Projektionsbilds pS (x, y),
- – Bestimmung
der lokalen Rauschvarianz,
- – adaptive
Bestimmung der lokalen Rauschfiltermaskengröße,
- – adaptive
Rauschfilterung des separierten Korrekturterms: Φ(δpcor),
- – Addition
des rauschgefilterten Korrekturterms Φ(δpcor)
zu den (logarithmierten) gemessenen „totalen" Projektionsdaten pS.
-
Die
beschriebene Bildqualitätsanalyse
und der erfindungsgemäße Lösungsvorschlag
haben bezüglich der
Abspaltung und adaptiven Rauschfilterung des Streustahlungskorrekturterms
die folgenden Vorteile:
- – durch die adaptive Rauschfilterung
wird die Inhomogenität
des Bildrauschens, bestehend aus unregelmäßigen Strichen und langgezogenen „Rauschwürmern" statt isotropem „Salz und
Pfeffer", homogenisiert,
- – die
spezifische Rauschfilterung des Streustrahlungskorrekturterms reduziert
die Zunahme des Rauschens durch die Streustrahlungskorrektur beträchtlich,
- – die
kontrastrestaurierende Wirkung der Streustrahlungskorrektur kann
weitgehend erhalten bleiben,
- – dadurch
wird im Endeffekt das differenzielle Kontrast-zu-Rausch-Verhältnis dCNR gegenüber der
Standardkorrektur ohne Abspaltung und Rauschfilterung deutlich verbessert
und dem Idealfall ohne Streustrahlung angenähert,
- – der
unkorrigierte Anteil pS der Daten wurde
zunächst
unverändert
gelassen, um die Qualitätsverbesserung
unabhängig
von nachgeschalteten Bildnachverarbeitungsalgorithmen zu bewerten;
selbstverständlich können solche
bekannten Bildnachverarbeitungsalgorithmen, die gewöhnlich oder
bei Bedarf zum Einsatz kommen, auch weiterhin nach der erfindungsgemäßen Korrektur – dann auf
pcor,φ – angewendet
werden, und
- – der
Lösungsvorschlag
ist sowohl bei auf Messungen als auch auf Rechenmodellen basierenden
Streustrahlungskorrekturverfahren anwendbar.
-
Entsprechend
dem oben beschriebenen Grundgedanken des Erfinders, schlägt dieser
die Ausführung des
erfindungsgemäßen Verfahrens
in unterschiedlichen Varianten wie folgt vor:
- a)
Verfahren zur Korrektur von Streustrahlungsfehlern bei der Radiographie
und der Computertomographie, jeweils mit Flächendetektoren, mit den folgenden
Verfahrensschritten:
- – Schätzung einer
Streustrahlungsverteilung Scor(y, z),
- – Berechnung
eines Standardkorrekturterms δpcor,
- – Rauschfilterung
des Standardkorrekturterms δpcor mit F(δpcor),
- – Addition
des rauschgefilterten Standardkorrekturterms F(δpcor)
zu den logarithmierten gemessenen totalen Projektionsdaten pS.
- b) Verfahren zur Korrektur von Streustrahlungsfehlern bei der
Radiographie und der Computertomographie mit Flächendetektoren mit den folgenden
Verfahrensschritten:
- – Schätzung einer
Streustrahlungsverteilung Scor(y, z),
- – Berechnung
eines Standardkorrekturterms δpcor,
- – Bestimmung
einer lokalen Rauschvarianz,
- – adaptive
Bestimmung einer lokalen Rauschfiltermaskengröße,
- – adaptive
Rauschfilterung des Standardkorrekturterms δpcor mit Φ(δpcor),
- – Addition
des rauschgefilterten Korrekturterms Φ(δpcor)
zu den logarithmierten gemessenen totalen Projektionsdaten PS.
-
Hier
wird also in Abwandlung zum Verfahren a) die Rauschfilterung adaptiv
abhängig
von der zuvor bestimmten lokalen Rauschvarianz durchgeführt.
- c) Verfahren zur Korrektur von Streustrahlungsfehlern
bei der Radiographie und der Computertomographie mit Flächendetektoren
mit den folgenden Verfahrensschritten:
- – Schätzung einer
Streustrahlungsverteilung Scor(y, z),
- – Berechnung
eines Standardkorrekturterms δpcor,
- – Bestimmung
einer lokalen Rauschvarianz,
- – Bestimmung
eines lokalen Rauschzuwachsterms, wobei
- – eine
adaptive Rauschkorrektur ab einer vorbestimmten Größe des Rauschzuwachsterms
durchgeführt wird,
mit:
- – einer
adaptiven Bestimmung einer lokalen Rauschfiltermaskengröße,
- – einer
adaptive Rauschfilterung des Standardkorrekturterms δpcor mit Φ(δpcor), und
- – einer
Addition des rauschgefilterten Korrekturterms Φ(δpcor)
zu den logarithmierten gemessenen totalen Projektionsdaten pS,
- – und
eine nicht-adaptive Rauschkorrektur erfolgt, wenn der Rauschzuwachsterms
die vorbestimmte Größe des Rauschzuwachsterms
nicht erreicht, mit:
- – einer
Rauschfilterung des Standardkorrekturterms δpcor mit
F(δpcor), und
- – Addition
des rauschgefilterten Standardkorrekturterms F(δpcor)
zu den logarithmierten gemessenen totalen Projektionsdaten pS.
-
Dieser
Variante c) stellt eine Kombination der Verfahrensvarianten a) und
b) dar, wobei die Anwendung der jeweiligen Verfahrensvariante durch
die Größe des Rauschzuwachsterms
getriggert wird.
- d) Verfahren zur Korrektur
von Streustrahlungsfehlern bei der Radiographie und der Computertomographie mit
Flächendetektoren
mit den folgenden Verfahrensschritten:
- – Schätzung einer
Streustrahlungsverteilung Scor(y, z),
- – Berechnung
eines Standardkorrekturterms δpcor,
- – Segmentierung
eines urikorrigierten Projektionsbilds pS(x,
y),
- – Bestimmung
einer lokalen Rauschvarianz,
- – segmentabhängige Bestimmung
einer lokalen Rauschfiltermaskengröße,
- – segmentabhängige Rauschfilterung
des Standardkorrekturterms δpcor mit Φ(δpcor),
- – Addition
des rauschgefilterten Korrekturterms Φ(δpcor)
zu den logarithmierten gemessenen totalen Projektionsdaten PS.
-
In
der Verfahrensvariante d) wird schließlich das Verfahren gemäß Variante
b) dahingehend modifiziert, dass die adaptive Rauschfilterung auf
ein zuvor segmentiertes Projektionsbild angewendet wird. Hierdurch
kann nun eine segmentabhängig
Filterung durchgeführt
werden, indem beispielsweise an durch die Segmentierung erkannten
Kanten eine richtungsorientierte Filterung stattfindet und/oder
in Bereichen von Knochen und Gewebe unterschiedlich starke Rauschfilterungen
durchgeführt
werden. Ergänzend
kann auch die segmentabhängige
Rauschfilterung zumindest teilweise adaptiv in Abhängigkeit
lokaler Eigenschaften ausgeführt
werden.
-
Die
adaptive Rauschfilterung wird bei den oben genannten Verfahren bevorzugt
in Abhängigkeit
von der zuvor bestimmten lokalen Rauschvarianz durchgeführt.
-
Bei
allen oben genannten Varianten kann die Berechnung des Standardkorrekturterms δpcor mit der folgenden Gleichung durchgeführt werden: δpcor = – ln(1 – Scor/PS)wobei
T der normierten totalen Intensitätsverteilung einschließlich Streustrahlung
entspricht.
-
Alternativ
kann die Berechnung des Standardkorrekturterms δpcor auch
gemäß der folgenden
Gleichung durchgeführt
werden: δpcor = pcor – pS wobei pcor = – ln(PS) – ln(1 – Scor/PS) gilt, pS den logarithmierten totalen Projektionsdaten
und PS der normierten totalen Intensitätsverteilung
einschließlich
Streustrahlung entspricht.
-
Die
Schätzung
der Streustrahlungsverteilung Scor(y, z)
kann einerseits durch Messung, andererseits durch Berechnung ausgeführt werden.
-
Im
Rahmen der Erfindung liegt nicht nur das oben geschilderte Verfahren,
sondern auch projektive und tomographische Radiographiesysteme,
welche eine Recheneinheit mit gespeicherten Programmen aufweisen,
wobei mindestens ein Programm die oben geschilderten Verfahren ausführt.
-
Außerdem liegt
auch im Rahmen der Erfindung dedizierte Recheneinheiten, in denen
eine Bildverarbeitung von projektiven und tomographischen Aufnahmen
stattfindet, wobei die Recheneinheit einen Speicher mit Programmcode
aufweist, welcher im Betrieb die oben geschilderten Verfahren ausführen.
-
Im
Folgenden wird die Erfindung mit Hilfe der Figuren näher beschrieben,
wobei nur die zum Verständnis
der Erfindung notwendigen Merkmale dargestellt sind und folgende
Bezugszeichen verwendet werden: 1: C-Bogen-System; 2:
Röntgenröhre; 3:
Detektor; 4: C-Bogen; 5: Patient; 6:
Gehäuse; 7:
Patientenliege; 8: Steuer- und Recheneinheit; 9: Speicher;
Prg1-Prgn: Computerprogramme.
-
Es
zeigen im Einzelnen:
-
1:
C-Bogen-System mit einem Flächendetektor
und einem Steuer- und Rechensystem zur Bildaufbereitung;
-
2:
Intensitätsbild
I(y, z);
-
3:
logarithmisches Projektionsbild
-
4:
Intensitätsbild
ohne Streustrahlung;
-
5:
Intensitätsbild
mit aufaddierter Streustrahlung;
-
6:
Intensitätsbild
mit Streustrahlungskorrektur;
-
7:
Profilplot des Intensitätsbildes
aus 4;
-
8:
Profilplot des Intensitätsbildes
aus 5;
-
9:
Profilplot des Intensitätsbildes
aus 6;
-
10:
logarithmisches Projektionsbild ohne Streustrahlung;
-
11:
logarithmisches Projektionsbild mit aufaddierter Streustrahlung;
-
12:
logarithmisches Projektionsbild mit rauschgefilterter Streustrahlungskorrektur
mit einer 2-dimensionalen
Medianfilterung;
-
13:
logarithmisches Projektionsbild mit Streustrahlungskorrektur gemäß Stand
der Technik;
-
14:
Profilplot des logarithmischen Projektionsbildes aus 10;
-
15:
Profilplot des logarithmischen Projektionsbildes aus 11;
-
16:
Profilplot des logarithmischen Projektionsbildes aus 12;
-
17:
Profilplot des logarithmischen Projektionsbildes aus 13;
-
18:
logarithmisches Projektionsbild ohne Streustrahlung;
-
19:
logarithmisches Projektionsbild mit aufaddierter Streustrahlung;
-
20:
logarithmisches Projektionsbild mit rauschgefilterter Streustrahlungskorrektur
mit einer 1-dimensionalen
Medianfilterung;
-
21:
logarithmisches Projektionsbild mit Streustrahlungskorrektur gemäß Stand
der Technik.
-
In
der 1 ist ein C-Bogen-System 1, mit dem projektive
und tomographische Röntgenaufnahmen erstellbar
sind, dargestellt. Es weist einen, an einem Gehäuse 6 befestigten,
C-Bogen 4 auf, der um einen Patienten 5 bewegt
werden kann, wobei an einem Ende des C-Bogens eine Röntgenröhre 2 und
am anderen Ende des C-Bogens ein Flächendetektor (Flat Panel Detektor
FPD) 3 angeordnet sind. Der Patient 5 befindet sich
auf einer verfahrbaren Patientenliege 7, so dass der Patient 5 mit
dem zu betrachtenden Teilbereich jeweils in das Messfeld zwischen
Röntgenröhre 2 und
Detektor 3 des C-Bogensystems verfahren werden kann. Die
Steuerung des System übernimmt
eine Steuer- und Recheneinheit 8 mit einem Speicher 9,
in dem entsprechende Steuer- und Bildverarbeitungsprogramme Prg1-Prgn bei Bedarf ablaufen können, wobei sich hierin auch
Programmcode zur Durchführung
des erfindungsgemäßen Verfahrens
befindet, so dass die in diesem System auftretende Streustrahlung
entsprechend korrigiert werden kann.
-
Die
folgenden Abbildungen der 2-17 basieren
auf Simulationsrechnungen, wobei ein Simulationsobjekt (Phantom),
eine homogene Platte aus wasserähnlichem
Material, z.B. Wasser oder Plexiglas, mit einem linearen Schwächungskoeffizient μ = 0,25/cm,
einer Schichtdicke D = 20 cm, entsprechend p =μD = 5 verwendet wird. In der
Mitte der Platte ist ein kleines Quadrat aus demselben Material
mit nur 0,2 cm Dicke angeordnet, um einem schwachen Kontrastsprung
um relativ 1/1000 zu simulieren. Es ergibt sich damit: p + ∆p =
5,005.
-
Das
Projektionsbild ist in der 3 wiedergegeben.
-
Die
Bildgröße beträgt 1000 × 1000 Pixel
und das Niedrigkontrastquadrat umfasst die inneren 200 × 200 Pixel.
-
Vereinfachend
wird angenommen, dass je Detektorpixel bei ungeschwächter Strahlung
ohne Phantom im Mittelwert effektiv 400000 Röntgen-Quanten vollständig absorbiert
und in ein Messsignal umgewandelt werden, wobei der physikalische
Detektionsprozess im Einzelnen ohne Bedeutung ist.
-
Die
Strahlungsschwächung
in der Platte reduziert die Anzahl der registrierten Quanten um
den Faktor exp(–5)
außerhalb
des Niedrigkontrastquadrats auf 2695 beziehungsweise exp(–5.005)
innerhalb des Quadrats auf 2682 Quanten. Das Intensitätsbild ohne
Rauschen ist in der 5 wiedergegeben.
-
Das
Rauschen genügt
der Poisson-Statistik, wobei die Standardabweichung der Quadratwurzel
des Mittelwertes entspricht. Bei der Simulation wurde ein Zufallsgenerator
für unkorrelierte
normalverteilte Zufallszahlen benutzt. Bei den relativ großen Mittelwerten
kann die Poissonverteilung durch Normalverteilung approximiert werden.
-
Zur
Demonstration der Auswirkung einer Streustrahlungsintensität von gleicher
Stärke
wie die Primärintensität wird eine
Streustrahlungshintergrundintensität entsprechend 2682 Quanten
je Pixel angenommen, die sich im Mittelwert additiv überlagert.
Die Standardabweichung des Rauschens nimmt dabei um den Faktor √2 zu.
-
Die
gezoomten Intensitätsbilder
ohne und mit Streustrahlung sind mit Rauschen in den 5 und 6 dargestellt.
-
Die 2 zeigt
damit ein Intensitätsbild
mit einem Erwartungswert, ohne Rauschen und ohne Streustrahlung
nach der Transmission durch das simulierte Phantom. Die 6 zeigt
ein Projektionsbild nach Normierung und Logarithmierung.
-
In
den folgenden Bildern ist der Bereich mit dem schwachen Kontrastquadrat
gezoomt dargestellt.
-
Bei
den 4 bis 5 und 7 bis 9 handelt
es sich noch um Intensitätsbilder.
Die 4 zeigt das Intensitätsbild ohne Streustrahlung,
in der 5 liegt eine Überlagerung
eines konstanten Streustrahlungshintergrunds vor, wobei das Grauwertfenster
entsprechend so verschoben wurde, dass der mittlere Grauwert gleich
bleibt. Die Zunahme des Rauschens durch die Streustrahlung und die
dadurch schlechtere Erkennbarkeit des Quadrats ist hier auffällig.
-
Die 9 zeigt
eine simulierte Streustrahlungskorrektur durch digitale Subtraktion
des konstanten Erwartungswerts S des Streustrahlungshintergrunds.
Das Grauwertfenster wurde jetzt identisch wie in 7 gewählt. Der
Bildeindruck ist identisch mit dem von 8.
-
Die 7 bis 9 sind
Profilplots durch die mittleren Bildzeilen in den entsprechenden
Abbildungen in den 4 bis 6.
-
Ab
der 10 werden logarithmische Projektionsdaten und
die Wirkung des von uns vorgeschlagenen Verfahrens dargestellt.
Die 12 bis 14 sind
die zu den Intensitätsbildern
aus den 4 bis 6 korrespondierenden
Bilder der normierten und logarithmierten Projektionsdaten, wobei
die 10 die logarithmischen Projektionsdaten ohne Streustrahlung
zeigen, in der 12 eine Überlagerung des Streustrahlungshintergrunds
gezeigt wird. Hier ist das Grauwertfenster unter Beibehaltung der
Fensterbreite entsprechend so verschoben, dass der mittlere Grauwert
gleich bleibt. Die Verminderung des Kontrasts des inneren Quadrats ist
auffällig.
In der 12 ist die 11 nach
einer Standardstreustrahlungskorrektur gemäß dem Stand der Technik gezeigt,
wobei das Grauwertfenster identisch wie in 10 ist.
Die Zunahme des Rauschens gegenüber
der 10 ist auffällig.
-
Abschließend folgt
das Ergebnis einer rauschgefilterten Streustrahlungskorrektur gemäß dem erfindungsgemäßen Vor schlag
in der 13. Als Rauschfilterung wurde
hier eine 2-dimensionale Medianfilterung mit 3 × 3-Maskengröße verwendet.
Die Rauschreduktion gegenüber
der 12 ist deutlich erkennbar, jedoch erscheint das
verbleibende Bildrauschen etwas grobkörniger als das ursprüngliche
Bildrauschen in der 10.
-
Die 14 bis 17 stellen
nochmals Profilplots durch die mittleren Bildzeilen in den entsprechenden
Abbildungen in den 10 bis 13. Aus
ihnen lässt
sich bereits qualitativ erkennen, dass die Signaldifferenz durch
die Streustrahlung etwa auf die Hälfte „gestaucht" wird – siehe 15 versus 14 –, das Rauschen
durch die Standardkorrektur in 16 zunimmt
und, wie aus 17 erkennbar ist, durch die
erfindungsgemäße rauschgefilterte
Korrektur wieder etwa auf die ursprüngliche Stärke „gedämpft" wird.
-
Die 18 bis 20 entsprechen
den 10 bis 12 und
wurden zur besseren Vergleichsmöglichkeit
nochmals dargestellt. In der 21 wurde
der Korrekturterm einer zeilenweisen 1-dimensionalen Medianfilterung
mit einer Maskenlänge
von 3 Pixel unterzogen. Gegenüber
der 13 ist der Rauscheindruck feinkörniger und
wieder ähnlich
wie in der 10 beziehungsweise in der 18.
-
Zum
Zweck einer statistischen Auswertung wurden 10 Simulationsdurchläufe durchgeführt und
jeweils Bilder entsprechend den 18 bis 20 erstellt.
Die Bilder der verschiedenen Simulationsdurchläufe (Stichproben) sind statistisch
unabhängig.
Für die
Auswertung wurden jeweils Mittelwerte über die 10 Simulationen gebildet:
-
Es
bezeichnen in der nachfolgenden Tabelle:
- „Signalhintergrund" = Mittelwert über 200 × 200 Pixel
außerhalb
des Quadrats und gemittelt über
10 Bilder;
- „Signaldifferenz” = dito
innerhalb des Niedrigkontrastquadrats minus Signalhintergrund;
- „Rauschen" = Standardabweichung
des Pixelrauschens
- dCNR = Signaldifferenz/Rauschen.
-
| |
Signalhintergrund
abs | (rel) |
Signaldifferenz
abs | (rel) |
Rauschen
abs | (rel) |
dCNR
abs | (rel) |
| a)
ohne Streustrahlung |
5.00 (1.0) |
0.00499
(1.0) |
0.0194
(1.0) |
0.258
(1.0) |
| b)
mit Streustrahlung |
4.31 (0.86) |
0.00243 (0.49) |
0.0137
(0.71) |
0.178
(0.69) |
| c)
Standardkorrektur |
5.00
(1.0) |
0.00485 (0.97) |
0.0274
(1.41) |
0.177
(0.69) |
| d)
rauschgefilterte Korr. mit idim. Medianfilterung |
5.00 (1.0) |
0.00485 (0.97) |
0.0199 (1.41) |
0.244 (0.95) |
-
Die
statistische Genauigkeit der Auswertung der Stichprobe liegt bei
etwa 3 %. Man erkennt die theoretisch zu erwartenden Relationen
bei der Signaldifferenz, nämlich
einen Verlust auf die Hälfte
durch Streustrahlung, beziehungsweise beim Rauschen (Faktor √2) beziehungsweise beim dCNR
(Faktor 1/√2) deutlich wieder.
-
Die
rechte Spalte zeigt, dass man sich mit dem vorgeschlagenen Verfahren
dem Idealfall ohne Streustrahlung auf etwa 95 % in Bezug auf das
dCNR = differentielle Kontrast-zu-Rausch-Verhältnis
nähert.
-
Insgesamt
wird also in dieser Erfindung ein Verfahren zur Korrektur von Streustrahlungsfehlern
bei der Radiographie und der Computertomographie, jeweils mit Flächendetektoren,
vorgestellt, wobei zunächst
eine rechnerische oder gemessene Schätzung einer Streustrahlungsverteilung
stattfindet, danach ein Standardkorrekturterms berechnet wird, eine
Rauschfilterung des Standardkorrekturterms durchgeführt wird
und anschließend
eine Addition des rauschgefilterten Standardkorrekturterms zu den
logarithmierten gemessenen totalen Projektionsdaten stattfindet.
Die Rauschfilterung kann auch adaptiv in Abhängigkeit von einer zuvor bestimmten
lokalen Rauschvarianz durchgeführt
werden. Weiterhin kann über
die Detektordaten auch eine Kombination von einfacher und adaptiver
Rauschfilterung, vorzugsweise mit einer Entscheidung für die lineare
oder adaptive Filterung in Abhängigkeit
eines zuvor bestimmten Rauschzuwachsterms, durchgeführt werden.
Alternativ besteht auch die Möglichkeit
zunächst
eine Segmentierung der, vorzugsweise unbearbeiteten, Projektionsdaten
durchzuführen
und die Rauschfilterung segmentabhängig adaptiv auszuführen.
-
Es
versteht sich, dass die vorstehend genannten Merkmale der Erfindung
nicht nur in der jeweils angegebenen Kombination, sondern auch in
anderen Kombinationen oder in Alleinstellung verwendbar sind, ohne
den Rahmen der Erfindung zu verlassen.
-
Literaturangaben:
-
[JS1982]
-
- P. M. Joseph, R. D. Spital: "The effects of scatter in x-ray computed
tomography", Med.
Phys.9(4), July/August 1982, 464-472;
-
[Jäh1997]
-
- Bernd Jähne: "Digitale Bildverarbeitung", Springer-Verlag,
1997, 4. Auflage;
-
[SDB2006]
-
- J. H. Siewerdsen, M. J. Daly, B. Bakhtiar, D. H. Moseley,
S. Richard, H. Keller, D. A. Jaffray: "A simple, direct method for x-ray scatter
estimation and correction in digital radiography and cone-beam CT", Med. Phys.33(1), January
2006, 187-197;
-
[ZSR2005]
-
- M. Zellerhoff, B. Scholz, E.-P. Rührnschopf, T. Brunner: "Low contrast 3D reconstruction
from C-arm data",
Proceedings of SPIE. Medical Imaging 2005, Vol. 5745, 646-655.
und
