DE102006034132A1

DE102006034132A1 - Verfahren zur parameterfreien Bildbearbeitung

Info

Publication number: DE102006034132A1
Application number: DE102006034132A
Authority: DE
Inventors: Kurt Dr. Osterloh; Uwe Dr. Ewert; Uwe Dr. Zscherpel; Oleksandr Alekseychuk
Original assignee: Bundesanstalt fuer Materialforschung und Pruefung BAM
Current assignee: Bundesanstalt fuer Materialforschung und Pruefung BAM
Priority date: 2006-07-20
Filing date: 2006-07-20
Publication date: 2008-01-24
Also published as: EP2044572A2; WO2008009685A2; US20100040301A1; WO2008009685A3

Abstract

Um ein Verfahren anzugeben, welches einerseits die Erkennbarkeit feiner Strukturen in einem Eingangsbild erhöht und andererseits das Rauschen unterdrückt, sind folgende Verfahrensschritte vorgesehen:
- Berechnung einer Fouriertransformierten

des Eingangsbildes;
- Berechnung eines skalierten Powerspektrums P_sk(k) aus der Fouriertransformierten;
- Berechnung einer Rauschfunktion R(k), wobei für deren hohe Frequenzen k oberhalb einer vorbestimmten Grenzfrequenz k_s gilt:

- Berechnung einer Differenzfunktion D(k) = P_sk(k) - R(k);
- Berechnung einer Filterfunktion gemäß

- Multiplikation der Fouriertransformierten

mit der Filterfunktion F(k) und anschließende Durchführung einer Rück-Fouriertransformation (30) zur Berechnung eines rekonstruierten Bildes

Description

Aus dem Stand der Technik ist eine Vielzahl von Bildbearbeitungsalgorithmen und fertigen Computerprogrammen bekannt, welche zum Ziel haben, Rauschen in digitalisierten Eingangsbildern zu unterdrücken und die Erkennbarkeit feiner Details zu erhöhen.
Allen gemeinsam ist, dass die auszuwählenden Bereiche über Parameter eingestellt werden müssen, beispielsweise durch die Eingabe von Zahlen oder durch das Verschieben von virtuellen Regelknöpfen in einem grafischen User Interface. Da der Benutzer die optimale Parameterwahl nicht im Vorhinein kennt und durch suboptimale Einstellungen unerwünschte Ergebnisse entstehen, werden häufig Standardeinstellungen vorgegeben. Bei recht unterschiedlichen Bildvorgaben ist eine derartig pauschalisierte Vorgabe in den meisten Fällen eine unbefriedigende Lösung.
Zu der einfachsten Parameterwahl gehören Helligkeits- und Kontrasteinstellungen, zumal diese bereits zum visuellen Abgleichen eines Bildschirmes gehören. Um weitere Informationen, die hierbei noch nicht deutlich hervortreten, aus dem Bild herauszuholen, sind für digitale Bilder Filtertechniken entwickelt worden, die Einfluss auf das Signal/Rauschverhältnis nehmen und benachbarte Bildpunkte zueinander in Beziehung setzen. Dazu gehören Glättungsverfahren für erkennbar verrauschte Bilder sowie Kantenanhebungen mittels Differenzenbildung (z. B. Laplace-Filterung, Embossed beziehungsweise Pseudoplastdarstellungen). Der Grad, zu welchem eine Glättung oder Hervorhebung kantiger Strukturen erreicht wird, wird durch die Einstellung zahlenmäßiger Parameter gesteuert, welche der Benutzer vornimmt. Oft ist auch durch den Benutzer eine zusätzliche Entscheidung über die Form der jeweiligen Filter zu treffen. All diese Einstellungen sind vom Benutzer schwerer beherrschbar als die anschaulicheren Helligkeits- oder Kontrastregelungen. Optimale Ergebnisse werden oft nur erzielt, wenn der Benutzer eine angemessene Erfahrung gesammelt hat. Auch die aufeinander folgende Anwendung mehrerer unterschiedlicher Verfahren führt nur selten zu einer optimalen Lösung. Wiederholte Bearbeitungen ein und desselben Bildes führen erfahrungsgemäß zu recht unterschiedlichen Ergebnissen.
In der digitalen Radiographie stehen oftmals feine Details im Vordergrund, weniger großflächige Ausleuchtungen. Letztere können sogar stören, wenn in einem Bereich gut aus geleuchtete feine Strukturen in helleren Bildteilen überstrahlt werden und in dunkleren verdeckt bleiben. Aufgrund ihrer hohen Helligkeitsdynamik, die sich in digitalen Bildern über mehrere Größenordnungen erstrecken kann, sind gleichartige Feinstrukturen gleichermaßen in hellen und dunklen Bildbereichen vorhanden. Die Bildbearbeitung hat in einem solchen Fall die Aufgabe, großflächige dynamische Schwankungen herauszunehmen, ohne Details zu verwischen. Eine mögliche Technik besteht darin, vom ursprünglichen Bild ein stark geglättetes zu subtrahieren. Dazu wird aber ein Glättungsparameter gefordert. Komplexere Hochpassfilter brauchen zweierlei Angaben: die Filterkerngröße und die Filterform. Insgesamt gesehen gibt es eine breite Auswahl von Methoden und Vorgehensweisen zur digitalen Filterung.
Die Detailerkennbarkeit im Bild wird durch das Rauschen begrenzt. In ungünstigen Fällen, z. B. bei der Radiographie durch dicke Schichten mit hohem Anteil Streustrahlung, können selbst gröbere Details nicht mehr zweifelsfrei erkannt werden. Konkret tritt dieses Problem bei der Durchstrahlung von Gussteilen auf. Zur besseren Erkennbarkeit von Fehlern werden die Bilder in der Regel hochpassgefiltert, wodurch die Wirkungen des Rauschens jedoch verschärft werden. Um dies zu verhindern, können verschiedene Verfahren angewendet werden, unter denen z. B. die Median- oder Gaußfilterung zu nennen ist. Die Wahl der Filterparameter ist dabei nicht frei von Willkür und die durchgeführte Filterung ist objektiv selten optimal oder reproduzierbar.
Die DE 103 25 632 offenbart ein Verfahren zur Verbesserung der Erkennbarkeit von Strukturen auf Durchstrahlungsbildern, bei welchem die Intensitätsverteilung eines Durchstrahlungsbildes Fourier-transformiert wird und die Frequenzraum-Intensitätsverteilung gefiltert wird, indem die Gewichtung zwischen hochfrequenten und niederfrequenten Bildsignalanteilen verändert wird, wobei die Festlegung der stärker zu gewichtenden Bildsignalanteile unter Berücksichtigung einer mittleren Strukturgröße der Strukturen erfolgt, deren Erkennbarkeit verbessert werden soll. Dadurch können beispielsweise Weichteilgewebestrukturen in Röntgenbildern erkannt werden.
Nachteil eines solchen Vorgehens ist, dass ein Parameter festgelegt werden muss, welcher von der mittleren Strukturgröße der aufzulösenden Strukturen abhängt, und dass die Rauschproblematik im Allgemeinen verschlechtert wird.
Aufgabe der Erfindung ist es, ein Verfahren anzugeben, welches einerseits die Erkennbarkeit feiner Strukturen in einem Eingangsbild erhöht und andererseits das Rauschen unterdrückt. Insbesondere ist es die Aufgabe der Erfindung, ein Verfahren zur Verfügung zu stellen, welches keine Parametereinstellung oder Menuauswahl vom Benutzer verlangt. Die Bearbeitung des Bildes soll in einem einzigen Durchgang erfolgen und nicht mit Iterationen arbeiten, für die Abbruch- oder Endkriterien festzulegen wären.
Erfindungsgemäß wird diese Aufgabe mittels eines Verfahrens mit den im Anspruch 1 genannten Merkmalen gelöst.
Das erfindungsgemäße Verfahren zur Erarbeitung eines Eingangsbildes I(r) umfasst folgende Verfahrensschritte: Berechnung einer Fouriertransformierten I(k) des Eingangsbildes; Berechnung eines skalierten Powerspektrums P_sk(k) aus der Fouriertransformierten; Berechnung einer Rauschfunktion R(k), wobei für deren hohe Frequenzen k oberhalb einer vorbestimmten Grenzfrequenz k_s gilt: R(k) = P_sk(k) für ||k||>k_s; Berechnung einer Differenzfunktion D(k) = P_sk(k) – R(k); Berechnung einer Filterfunktion gemäß F(k) = α||k|| D(k); Multiplikation der Fouriertransformierten I(k) mit der Filterfunktion F(k) und anschließende Durchführung einer Rück-Fouriertransformation zur Berechnung eines rekonstruierten Bildes Ir(r).
Dabei haben von den genannten Relationen geringfügig abweichende Funktionenwerte, also solche, die innerhalb eines Bereiches von 10 % von den beanspruchten Funktionenwerten differieren, gleiche Wirkung und sind im Rahmen der Erfindung als äquivalent anzusehen. Durch die Subtraktion der Rauschverteilung vom skalierten Powerspektrum, welche für hohe Frequenzen direkt aus den Werten des skalierten Powerspektrum gewonnen wird, einerseits und der Multiplikation der Differenzfunktion mit einer im Frequenzraum linear ansteigenden Funktion andererseits wird ein Kompromiss erreicht: das zumeist hochfrequente Rauschen im Eingangsbild wird unterdrückt, während gleichzeitig feine Strukturen auf Kosten von großflächigen dynamischen Unterschieden hervorgehoben werden. Die Dynamik des gesamten Bildes wird bei dem erfindungsgemäßen Verfahren berücksichtigt und nicht nur die von Teilbereichen. Die Bearbeitung des Bildes erfolgt in einem einzigen Durchgang, wobei der Benutzer keine Bearbeitungsparameter bestimmen muss. Anschließend dürften höchstens Feinkorrekturen der Helligkeits- oder Kontrasteinstellungen erforderlich sein. Die Filterfunktion wird auf Grundlage des im Bild vorhandenen Rauschens ermittelt, so dass das Verfahren exakt reproduzierbar ist.
Im Rahmen der Erfindung bezieht sich der Begriff skaliertes Powerspektrum auf eine Funktion, welche aus dem einfachen Powerspektrum der Fouriertransformierten durch einfache Skalierung oder Transformation hervorgeht, z. B. durch Logarithmierung oder durch Wurzelbildung, wobei das skalierte Powerspektrum im letzteren Falle das Amplitudenspektrum der Fouriertransformierten darstellen würde. Oberhalb der vorbestimmten Grenzfrequenz k_s wird der Funktionsverlauf des skalierten Powerspektrums als reines Rauschen interpretiert. Die Rauschfunktion R(k), welche eine Abschätzung der Rauschverteilung im Eingangsbild darstellt, übernimmt für hohe Frequenzen oberhalb der vorbestimmten Grenzfrequenz somit die gleichen Funktionswerte wie das skalierte Powerspektrum.
Vorzugsweise gilt für die Grenzfrequenz k_s: k_s = k_max/2, wobei k_max die durch das Sampling des Eingangsbildes I(r) bestimmte Maximalfrequenz in einer Koordinatenrichtung k_x oder k_y ist. Diese spezielle Wahl kann im Ortsbereich vor der Fouriertransformation dahin gehend interpretiert werden, dass alle Intensitätsschwankungen zwischen zwei benachbarten Bildpunkten als rauschbedingt angesehen werden, während Schwankungen über größere Distanzen hinweg als Information betrachtet werden.
In bevorzugter Ausgestaltung der Erfindung ist vorgesehen, dass bei der Berechnung der Rauschfunktion R(k) deren Funktionswerte für niedrige Frequenzen ||k||<k_s aus denen der Rauschfunktion für hohe Frequenzen ||k||>k_s extrapoliert werden. Somit werden alle Frequenzen, welche größer als die Grenzfrequenz sind, dem Rauschen zugeordnet und entsprechend in die niederfrequenten Bereiche fortgesetzt.
Bei Vorliegen von gleichmäßigem, isotropem Rauschen und einer ungerichteten Musterverteilung ist bevorzugt, dass eine der Funktionen P_sk(k), R(k), D(k) oder F(k) über alle Richtungen gemittelt wird. Dies bedeutet, dass für alle Frequenzen mit k = const. eine der genannten Funktionen die gleichen Werte aufweist. Nach Durchführung eines Mittelungsverfahrens an einer der genannten Funktionen, welche dadurch rotationssymmetrisch wird, sind auch alle anderen im weiteren Verfahrensverlauf berechneten Funktionen rotationssymmetrisch.
Bevorzugt wird bei einer der Funktionen P_sk(k), R(k), D(k) oder F(k) eine gleitende Mittelwertbildung durchgeführt. Die Glättung einer der genannten Funktionen bewirkt gleichzeitig auch die Glättung einer weiteren, im weiteren Verfahrensverlauf, berechneten Funktion.
Durch die gleitende Mittelwertbildung werden feine Oszillationen der zu berechnenden Filterfunktion vermieden, welche unausweichlich zu Artefakten führen würden.
Im Falle eines richtungsabhängigen (anisotropen) Signal-/Rauschverhältnisses ist in bevorzugter Ausführungsform vorgesehen, dass eine der Funktionen P_sk(k), R(k), D(k) oder F(k) zur Glättung mit einer Gaußfunktion gefaltet wird. Die Faltung mit einer Gaußfunktion bewirkt eine Art Mittelwertbildung sowohl in radialer als auch in axialer Richtung, wenn die Art der Rauschverteilung oder die Intensitätsverteilung im Eingangsbild keine rotationssymmetrische Filterfunktion nahe legt.
Weitere bevorzugte Ausgestaltungen der Erfindung ergeben sich aus den übrigen, in den Unteransprüchen genannten Merkmalen.
Die Erfindung wird nachfolgend in einem Ausführungsbeispiel anhand der zugehörigen Zeichnungen näher erläutert. Es zeigen:
1a, b Flussdiagramme zur Veranschaulichung des erfindungsgemäßen Verfahrens (A: isotroper Fall, B: anisotroper Fall) und
2a, b, c Funktionsverläufe der Funktionen P_sk(k), R(k), D(k) und F(k).
1a zeigt ein Flussdiagramm des erfindungsgemäßen Verfahrens, welches angewendet wird, wenn sowohl die Rauschverteilung als auch die interessierenden Strukturen im Eingangsbild hinreichend isotrop, also richtungsunabhängig sind. Im Verfahrensschritt 10 wird das erfindungsgemäße Computerprogramm vom Benutzer aufgerufen. Ein Eingangsbild, welches die Intensitätsverteilung I(r) aufweist und in digitaler Form vorliegt und n·n Bildpunkte umfasst, wird eingelesen (Verfahrensschritt 12). Im Verfahrensschritt 14 wird die Fouriertransformierte I(k) des Eingangsbildes berechnet, welche n·n komplexwertige Bildpunkte aufweist. Im Verfahrensschritt 16 wird die Fouriertransformierte I(k) dupliziert und in ein separates Array von Bildpunkten kopiert, aus der durch Anwendung des erfindungsgemäßen Verfahrens eine geeignete Filterfunktion berechnet werden kann (Verfahrensschritte 18, 20, 22, 24, 26). Dabei wird im Verfahrensschritt 18 ein skaliertes Powerspektrum P_sk(k) aus der Fouriertransformierten berechnet. Eine geeignete Wahl des skalierten Powerspektrums ist der Logarithmus des Betragsquadrates der Fouriertransformierten des Eingangsbildes. Da ein isotroper Filter berechnet werden soll, wird im Verfahrensschritt 20 das skalierte Powerspektrum über alle Richtungen gemittelt, wodurch eine rotationssymmetrische Funktion entsteht. Im Verfahrensschritt 22 wird eine Rauschkorrektur vorgenommen. Hierbei wird eine Rauschfunktion R(k) berechnet, die für hohe Frequenzen k>n/2 gleich dem skalierten Powerspektrum ist. Die niedrigen Frequenzen der Rauschfunktion werden geschätzt und aus dem Funktionsverlauf für hohe Frequenzen extrapoliert. Danach wird die Differenz zwischen dem skalierten Powerspektrum P_sk(k) und der Rauschverteilung R(k) gebildet. Die entstehende Funktion weist immer noch feine Oszillationen auf, welche unerwünscht sind und zu Artefakten in den Bildrekonstruktionen führen. Daher wird im Verfahrensschritt 24 ein gleitender Mittelwert benachbarter Bildpunkte gebildet. Die derart geglättete Funktion wird in einem Verfahrensschritt 26 mit einer linear in ||k|| steigenden Funktion multipliziert. Dies entspricht einer Hochpassfilterung. Zweidimensional dargestellt besteht die Funktion, mit welcher die Differenzfunktion D(k) multipliziert wird, in einer Ursprungsgerade, dreidimensional dargestellt kommt man durch Rotation dieser Funktion um die Achse durch den Ursprung zu der Form eines umgekehrten Kegels. In der Tomographie wird dieser Filter auch als Ramachandran-Lakshminarayanan-Filter bezeichnet. Die so berechnete Filterfunktion F(k) wird im Verfahrensschritt 28 mit der Fouriertransformierten I(k) multipliziert, rücktransformiert (Verfahrensschritt 30) und die Intensitätsverteilung I_r(r) des rekonstruieren Bildes in einem Verfahrensschritt 32 dem Benutzer angezeigt. Damit endet das erfindungsgemäße Verfahren (Verfahrensschritt 34).
1b zeigt ein Flussdiagramm zur Ausführung des erfindungsgemäßen Verfahrens für den Fall, dass keine rotationssymmetrische Filterfunktion F(k) berechnet werden soll, es sich vielmehr aufgrund einer anisotropen Rauschverteilung oder Intensitätsverteilung im Eingangsbild ohne Vorzugsrichtung nicht anbietet, eine Mittelung über alle Richtungen auszuführen. Die Verfahrensschritte 20 und 24 der 1a entfallen. Um dennoch eine geeignete Glättung der Filterfunktion F(k) zu erreichen, wird in zusätzlichen Verfahrensschritten 36A, B, C eine Faltung mit einer geeigneten Gaußfunktion durchgeführt. Im erfindungsgemäßen Computerprogramm wird eine Fouriertransformation durchgeführt (36A), die entstehende Funktion mit einer Gaußfunktion multipliziert (36B) und die so entstandene Funktion rück-Fourier-transformiert (36C).
Das im Verfahrensschritt 32 angezeigte rekonstruierte Bild hat in vielen Fällen extrem hohe und niedrige Einzelwerte. Eine brauchbare Helligkeits- und Kontrastanpassung kann dadurch vorgenommen werden, dass Extremwerte eliminiert werden und nur ein mittlerer Bereich (beispielsweise ein 95 %-Bereich der Intensitätsverteilung) ausgewählt und auf das Bildausgabeformat abgebildet wird.
2a zeigt einen beispielhaften Funktionsverlauf des skalierten Powerspektrums eines verrauschten Eingangsbildes entlang einer Frequenzkoordinate k_x, welche n Rasterpunkte umfasst. Das Beispiel bezieht sich auf ein eindimensionales Eingangsbild, kann aber leicht auf 2 oder mehr Dimensionen verallgemeinert werden. Eine geeignete Wahl des skalierten Powerspektrums ist der Logarithmus des Betragsquadrates der Fouriertransformierten des Eingangsbildes. Zur Eliminierung der schnellen Oszillationen des Funktionsverlaufs wird dieser geglättet und erhält eine Form, welche in 2b dargestellt ist (P_sk(k)). In einem weiteren Verfahrensschritt wird eine Rauschfunktion R(k) bestimmt, welche für k_x>n/2 mit P_sk(k) zusammenfällt. Für die niederfrequenten Funktionswerte wird der Funktionsverlauf aus den hochfrequenten Funtionswerten geeignet extrapoliert. Da die Rauschverteilung in einer logarithmischen Darstellung für hohe Frequenzen näherungsweise linear ist, gilt dies auch für die hohen Frequenzen des skalierten Powerspektrums, so dass die Rauschfunktion R(k) in guter Näherung durch eine Gerade repräsentiert werden kann (vgl. 2b). In einem weiteren Verfahrensschritt wird die Differenzfunktion D(k) = P_sk(k) – R(k) bestimmt (2c), welche für k_x>n/2 verschwindet. Die Filterfunktion F(k) entsteht schließlich durch Multiplikation der Differenzfunktion D(k) mit einer linearen Funktion αk_x.

10: Start
12: Ausgabe des digitalisierten Eingangsbildes
14: Fouriertransformation
16: Duplikation
18: Berechnung des skalierten Powerspektrums
20: Mittelung über alle Richtungen
22: Rauschkorrektur
24: Berechnung des gleitenden Mittelwertes
26: Ram-Lac-Filter
28: Multiplikation
30: Rück-Fouriertransformation
32: Bildanzeige
34: Ende
36, B, C: Faltung mit Gaußkurve

Claims

Verfahren zur Verarbeitung eines Eingangsbildes I(r), umfassend folgende Verfahrensschritte: – Berechnung einer Fouriertransformierten I(k) des Eingangsbildes; – Berechnung eines skalierten Powerspektrums P_sk(k) aus der Fouriertransformierten; – Berechnung einer Rauschfunktion R(k), wobei für deren hohe Frequenzen k oberhalb einer vorbestimmten Grenzfrequenz k_s gilt: R(k) = P_sk(k) für ||k||k_s; – Berechnung einer Differenzfunktion D(k) = P_sk(k) – R(k); – Berechnung einer Filterfunktion gemäß F(k) = α||k|| D(k); – Multiplikation der Fouriertransformierten I(k) mit der Filterfunktion F(k) und anschließende Durchführung einer Rück-Fouriertransformation (30) zur Berechnung eines rekonstruierten Bildes I_r(r).
Verfahren nach Anspruch 1, dadurch gekennzeichnet, dass für das skalierte Powerspektrum P_sk(k) gilt: P_sk(k) = In |I(k)|².
Verfahren nach Anspruch 1 oder 2, dadurch gekennzeichnet, dass für die Grenzfrequenz k_s gilt: k_s = k_max/2, wobei k_max die durch das Sampling des Eingangsbildes I(r) bestimmte Maximalfrequenz in einer Koordinatenrichtung k_x oder k_y ist.
Verfahren nach einem der vorangehenden Ansprüche, dadurch gekennzeichnet, dass bei der Berechnung der Rauschfunktion R(k) deren Funktionswerte für niedrige Frequenzen ||k||<k_s aus denen der Rauschfunktion für hohe Frequenzen ||k||>k_s extrapoliert werden.
Verfahren nach einem der vorangehenden Ansprüche, dadurch gekennzeichnet, dass eine der Funktionen P_sk(k), R(k), D(k) oder F(k) über alle Richtungen gemittelt wird.
Verfahren nach einem der vorangehenden Ansprüche, dadurch gekennzeichnet, dass bei einer der Funktionen P_sk(k), R(k), D(k) oder F(k) eine gleitende Mittelwertbildung durchgeführt wird.
Verfahren nach einem der Ansprüche 1 bis 4, dadurch gekennzeichnet, dass eine der Funktionen P_sk(k), R(k), D(k) oder F(k) zur Glättung mit einer Gaußfunktion gefaltet wird.
Rechnerlesbares Medium mit einer oder mehreren darin gespeicherten Befehlsfolgen, um einen oder mehrere Prozessoren zum Durchführen der Verfahrensschritte, welche in einem der vorangehenden Ansprüche angegeben sind, zu veranlassen.
Computerprogramm, das es einem Computer ermöglicht, nachdem es in den Speicher des Computers geladen worden ist, ein Verfahren nach einem der Ansprüche 1 bis 7 durchzuführen.