DE69830494T2 - Verfahren zur Verbesserung von Artefakten in digitalen Bildern - Google Patents

Verfahren zur Verbesserung von Artefakten in digitalen Bildern Download PDF

Info

Publication number
DE69830494T2
DE69830494T2 DE69830494T DE69830494T DE69830494T2 DE 69830494 T2 DE69830494 T2 DE 69830494T2 DE 69830494 T DE69830494 T DE 69830494T DE 69830494 T DE69830494 T DE 69830494T DE 69830494 T2 DE69830494 T2 DE 69830494T2
Authority
DE
Germany
Prior art keywords
image
gradient
multiscale
representation
scale
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Lifetime
Application number
DE69830494T
Other languages
English (en)
Other versions
DE69830494D1 (de
Inventor
Koen Van De Velde
Piet Vuylsteke
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Agfa HealthCare NV
Original Assignee
Agfa Gevaert NV
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Agfa Gevaert NV filed Critical Agfa Gevaert NV
Application granted granted Critical
Publication of DE69830494D1 publication Critical patent/DE69830494D1/de
Publication of DE69830494T2 publication Critical patent/DE69830494T2/de
Anticipated expiration legal-status Critical
Expired - Lifetime legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/90—Dynamic range modification of images or parts thereof
    • G06T5/94—Dynamic range modification of images or parts thereof based on local image properties, e.g. for local contrast enhancement
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10116—X-ray image
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20016—Hierarchical, coarse-to-fine, multiscale or multiresolution image processing; Pyramid transform

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)
  • Facsimile Image Signal Circuits (AREA)

Description

  • ERFINDUNGSGEBIET
  • Die vorliegende Erfindung betrifft die digitale Bildverarbeitung. Insbesondere betrifft sie ein Verfahren zum Korrigieren von Artefakten, die aus einer nichtlinearen Behandlung von Komponenten (Koeffizienten) einer mehrskaligen Bilddarstellung herrühren.
  • Das Verfahren läßt sich zum Beispiel in einem Röntgensystem in der Medizin einsetzen, wie etwa in einem Computerröntgen-(CR)-System oder einem Computertomographie-(CT)-System.
  • ALLGEMEINER STAND DER TECHNIK
  • Bei Bildgebungssystemen, bei denen das letzte ausgegebene Bild von einem menschlichen Betrachter ausgewertet werden muß, kommt es zu einem Problem, wenn das Vorlagenbild, wie es von einer Bilderfassungseinrichtung erhalten worden ist, Detailinformationen mit unterschiedlichen Grobheitsgraden innerhalb eines breiten Amplitudenbereichs enthält. Diese Situation ist typisch, wenn der Sensor ein gutes Signal-Rausch-Verhältnis über einen großen Dynamikbereich hinweg aufweist, was bei CR und CT der Fall ist. In üblichen Arbeitsumgebungen wie etwa einem Krankenhaus ist keine Zeit, um für jedes individuelle Bild (wenn überhaupt) die optimale Fenster-/Pegeleinstellung auszuwählen, weshalb eine Notwendigkeit dafür besteht, ein einzelnes Bild auf einem Film oder auf einem Monitorschirm anzuzeigen, der alle relevanten diagnostischen Details mit einem annehmbaren Kontrast über den ganzen Dynamikbereich hinweg offenlegt. Dieses Problem wurde auf dem Gebiet der digitalen Radiographie erkannt, siehe:
    Maack I., Neitzel u., Optimized image processing for routine digital radiography, Proceedings International Symposium CAR '91, S. 109, Springer Verlag.
  • Es sind viele Verfahren zur Kontrastverstärkung entwickelt worden, wie etwa die gemeinhin bekannten Techniken der Unschärfenmaskierung (Cocklin M.L., Gourlay A.R., Jackson P.H., Kaye G., Kerr I.H. und Lams P., Digital processing of chest radiographs, Image and Vision Computing, Band 1, Nr. 2, 67–78 (1983)) und adaptive Histogrammodifikation (Pizer S.M., Amburn E.P., Austin J.D., Cromartie R., Geselowitz A., Greer T., Ter Haar Romery B., Zimmerman J.B. und Zuiderveld K., Adaptive histogram equalisation and its variations, Computer Vision, Graphics and Image Processing, Band 39, 355–368 (1987)). Diese traditionellen Verfahren sind mit zwei Problemen behaftet:
    • 1. Sie legen eine bestimmte Filtergröße fest, wohingegen diagnostische Details möglicherweise in ein und demselben Bild bei mehreren Skalen auftreten können. Es kann wünschenswert sein, Details im Bild bei jeder möglichen Skala zu verstärken.
    • 2. Artefakte werden in der Nähe von signifikanten Signalpegelübergängen erzeugt, z.B. an Grenzen zwischen Knochen und Weichgewebe innerhalb des Bilds. Diese Artefakte können die Diagnose beeinträchtigen.
  • Zur Überwindung des ersten Problems sind in den vergangenen letzten Jahren Mehrskalenverfahren für die Kontrastverstärkung entwickelt worden, siehe Vuylsteke P., Schoeters E., Multiscale Image Contrast Amplification (weiter als MUSICA-Algorithmus bezeichnet), Proceedings of SPIE, Band 2167, 551–560 und Lu J., Healy D.M., Contrast enhancement via multiscale gradient transformation, in Proceedings of SPIE: Wavelet applications, Orlando, FL, USA (1994);
    Fan J., Laine A., Contrast enhancement by multiscale and nonlinear operators, in Wavelets in medicine and biology, Akram Aldroubi und Michael Unser Hrsg., CRC Press (1996).
  • Während dieses erste Problem effektiv gelöst wird, scheinen diese Verfahren auch bezüglich des zweiten Problems eine gute Leistung zu zeigen. Der MUSICA-Algorithmus wird routinemäßig in mehreren Krankenhäusern auf der ganzen Welt eingesetzt, und es gibt keine Beschwerden über Artefakte in CR.
  • Es wurde jedoch festgestellt, daß man beim Einsatz dieser Verfahren an CT-Bildern sorgfältiger sein muß, da sie schärfere Signalpegelübergänge aufweisen und dennoch Artefakte dort entstehen könnten.
  • AUFGABEN DER ERFINDUNG
  • Eine Aufgabe der vorliegenden Erfindung besteht in der Bereitstellung eines Verfahrens, das Artefakte korrigiert, die aus einer nichtlinearen Behandlung von Komponenten (Koeffizienten) in einer mehrskaligen Darstellung herrühren.
  • Weitere Aufgaben der vorliegenden Erfindung ergeben sich aus der folgenden Beschreibung.
  • DARSTELLUNG DER ERFINDUNG
  • Die Aufgaben der vorliegenden Erfindung werden durch ein Verfahren wie in Anspruch 1 dargelegt gelöst.
  • Das Verfahren der vorliegenden Erfindung basiert auf einer mehrskaligen Gradientendarstellung, wie sie etwa beschrieben wird in: Mallat S., Zhong S., Characterization of signals from multiscale edges, IEEE Trans. on Pattern analysis and Machine Intelligence, Band 14, Nr. 7 (1992)). Eine ähnliche Technik ist aus EP-A-0 712 092 bekannt.
  • Das Verfahren weist zwei grundlegende Schritte auf. Der erste Schritt ist der nichtlineare Verarbeitungsschritt, z.B. ein Kontrastverstärkungsschritt (z.B. unter Verwendung eines Algorithmus ähnlich dem von Lu und Healy (siehe oben)).
  • Der zweite und neuartige Teil ist ein Artefaktkorrekturschritt. Mit diesem Schritt können Artefakte korrigiert werden, die von einer beliebigen nichtlinearen Behandlung der Koeffizienten in der mehrskaligen Gradientendarstellung verursacht werden, und seine Verwendung ist deshalb nicht notwendigerweise auf das Problem der Kontrastverstärkung beschränkt.
  • AUSFÜHRLICHE BESCHREIBUNG DER ERFINDUNG
  • Mehrskalige Gradientendarstellung
  • Ein digitales Bild f ist eine Abbildung von einem quadratischen Gitter von Pixeln P auf die realen Zahlen R. Die mehrskalige Gradientendarstellung ist eine Zerlegung des Bilds in eine Reihe (etwa K, K=1, 2, 3,...) von vektorwertigen Bildern auf P (d.h. Abbildungen von P auf R2) W →0ƒ ,..., W →K–1ƒ und ein Restbild SKƒ jeweils mit der Größe des Vorlagenbilds (keine Unterabtastung findet statt). Es sein S0ƒ=ƒ und für j = 1,..., K, Sjƒ wird durch Anwenden eines Glättungsfilters aus Sj–1ƒ erhalten, dann wird W →jƒ definiert als der (diskrete) Gradient von Sjƒ. Auf W →0ƒ ,..., W →K–1ƒ, SKƒ wird als die mehrskalige Gradientendarstellung (MGR) des digitalen Bilds f Bezug genommen.
  • Die Glättungs- und Gradientenfilterung werden durch eine zweidimensionale Faltung erhalten.
  • Eine effiziente Implementierung von solchen Filtern basiert auf dem nacheinander Anwenden eines horizontalen 1-D-Faltungsfilters, gefolgt von einem vertikalen 1-D-Faltungsfilter, angewendet auf das Ergebnis des ersten Filters. Wenn das horizontale Filter mit a und das vertikale Filter mit b bezeichnet wird, dann wird das resultierende 2-D-Filter üblicherweise mit a...b bezeichnet.
  • Die Filter bei mehreren Skalen werden aus einem Basisfilter mit der kleinsten Skala erzeugt. Die entsprechenden Filter bei nachfolgenden größeren Skalen j=1,...K werden erhalten durch Einsetzen von 2j–1 nullwertigen Koeffizienten zwischen zwei beliebigen Koeffizienten des Basisfilters.
  • Bei einer bevorzugten Ausführungsform wird ein Basisglättungsfilter h0 mit Koeffizienten (1/16, 4/16, 6/16, 4/16, 1/16) verwendet.
  • Der unterstrichene Filterabgriff wird als der Filterursprung definiert (d.h. der Mittenabgriff mit einem Index Null). Der Basisgradientenfilter wird gewählt als g0=(1/2, –1/2).
  • Die entsprechenden Filter bei nachfolgenden größeren Skalen werden erhalten durch Einsetzen von Nullkoeffizienten wie oben beschrieben:
    h1 = (1/16, 0, 4/16, 0, 6/16, 0, 4/16, 0, 1/16), g1 = (1/2, 0, –1/2),
    h2 = (1/16, 0, 0, 0, 4/16, 0, 0, 0, 6/16, 0, 0, 0, 4/16, 0, 0, 0, 1/16)
    g2 = (1/2, 0, 0, 0, –1/2).
  • Wenn ein Filter eine geradzahlige Anzahl von Abgriffen (z.B. g0) aufweist, dann weist das Filter bei der nächst größeren Skala eine ungeradzahlige Anzahl von Abgriffen auf und der Filterursprung wird um eine Position nach rechts verschoben.
  • Entsprechende Kombinationen dieser Filter in horizontaler und vertikaler Richtung werden ausgewählt, um die 2-D-Filter zu erzeugen, die für die MGR-Zerlegung benötigt werde: S0f = f Sjf = (hj–1 ...hj–1·Sj–1f, für 0<j≤K W →jf = ((gj ...δ0)·Sjf, (δ0...gj)·Sjf) (1)
  • Hier ist · der Faltungsoperator und δ0 ist das Kronecker-Delta (d.h. das 1-D-Filter, bei dem der Mittenabgriff einen Koeffizienten gleich 1 aufweist und alle anderen Abgriffe einen Koeffizienten gleich 0 aufweisen).
  • Somit bedeutet gj...δ0 die Filterung durch gj nur in der horizontalen Richtung und δ0...gj ist gleichwertig zu dem Anwenden des Filters gj nur in der vertikalen Richtung.
  • Sjf ist eindeutig eine geglättete Version von Sj–1f, während W →jf der (diskrete) Gradient von Sjf bei Skala j ist.
  • Wenn ein Bild f mit einem Filter gefaltet wird, wird das Bild an den Rändern auf kontinuierliche Weise durch Reflexion an den Bildrändern verlängert. Nach der Durchführung der Faltung mit diesem verlängerten Bild wird das Ergebnis wieder auf den Pixelsatz des Vorlagenbilds eingeschränkt.
  • Die Zerlegung ist in 1 schematisch dargelegt (in dieser Figur bezeichnen die hochgestellten Zahlen (1) und (2) die horizontale bzw. vertikale Komponente der Gradienten).
  • Es ist wohlbekannt, daß diese Zerlegung umkehrbar ist, d.h., das Vorlagenbild kann aus der Darstellung durch einen Rekonstruktionsvorgang, der mit der Bezeichnung inverser mehrskaliger Gradient (IMG) bezeichnet wird, ganz wiederhergestellt werden. Dazu werden zwei zusätzliche Filter P0 und q0 mit der Einschränkung benötigt, daß p0·h0+q0·g0 = δ0 (2)p0 und q0 sind 1-D-Basisfilter.
  • Die Filter pj und qj bei größeren Skalen j werden erzeugt durch Einsetzen von 2j-1 Nullkoeffizienten zwischen aufeinanderfolgenden Abgriffen, wie oben beschrieben.
  • Die IMG-Rekonstruktion wird durchgeführt, wie in Zong, S, „Edge representation from wavelet transform maxima", Doktorarbeit, New York University, 1990, wie folgt erläutert: Sjƒ=(pj⊗pj)·Sj+1ƒ+(qj⊗lj)·Wj ( 1)ƒ+(lj⊗qj)·Wj (2)ƒ, j≥0. (3)
  • Das 1-D-Filter lj ist definiert als lj=1/2(δ0+pj·hj).
  • (Hochgestellte Zahlen (1) und (2) bezeichnen horizontale und vertikale Richtung).
  • Das Bild wird aus dem niedrig aufgelösten Bild SKf und den Gradienten W →0ƒ, ..., W →K–1ƒ durch die wiederholte Verwendung dieser Gleichung für j=K–1,...,0, wiederhergestellt (siehe 2).
  • Die Filter p0 und q0 sollten ausreichend glatt sein, um beim Rekonstruieren der Bilder aus einem Satz von nichtlinear modifizierten Transformationskoeffizienten Diskontinuitäten zu vermeiden (wie das immer der Fall ist, wenn eine Kontrastverstärkung vorgenommen wird).
  • Bei einer bevorzugten Ausführungsform p0=(1/4, 0, 1/2, 0, 1/4) und q0=(–1/32,–5/32,–13/32,–25/32,25/32,13/32,5/32,1/32).
  • Diese Filter genügen der Einschränkung (2).
  • Bei einer ersten Ausführungsform bleibt die Anzahl der Pixel bei jeder Skala konstant. Eine derartige Darstellung wird üblicherweise als ein mehrskaliger Stapel bezeichnet.
  • Bei einer zweiten Ausführungsform wird die Anzahl der Pixel in jedem Schritt der Zerlegung durch Unterabtastung reduziert. Eine derartige Darstellung wird üblicherweise mit einer mehrskaligen Pyramide bezeichnet.
  • Die MGR-Zerlegung wird dann definiert als:
  • Figure 00070001
  • Hier bezeichnet ↓ das Herunterrechnen um einen Faktor von zwei, wobei nur die Pixel mit geradzahligen Koordinaten (2x, 2y) beibehalten werden, wobei (x, y) sich auf die Pixelkoordinaten im aktuellen Gitter bezieht.
  • Die Pixel des entstehenden Bilds werden von (2x, 2y) zu (x, y) umbenannt. Die Anzahl der Pixel wird in beiden Dimensionen in jeder Stufe der pyramidenförmigen Zerlegung halbiert. Die Rekonstruktion aus der MGR-Pyramide erfolgt mit der Gleichung
  • Figure 00070002
  • Der Rekonstruktionsvorgang (5) ist gleichwertig einer Unterabtastung des Ergebnisses der Rekonstruktion (3), wobei berücksichtigt wird, daß für jedes willkürliche Bild f gilt (p'0⊗p'0)↓ƒ=(p0⊗p0)ƒ, wenn p0 zwischen aufeinanderfolgenden Abgriffen Nullkoeffizienten aufweist.
  • Der letztere Zustand wird erfüllt, wenn p0 = (1/4, 0, 1/2, 0, 1/4) ist Gleichung (5) beschreibt, wie eine unterabgetastete Version ↓S'jf aus S'jf erhalten wird. Als nächstes wird das rekonstruierte Bild S'jf in vollständiger Größe aus ↓S'jf und W →jf auf die folgende Weise berechnet (W →jf ist der Gradient von S'jf). Zuerst werden die Pixel mit den Koordinaten (x, y) der unterabgetasteten Rekonstruktion ↓S'jf (im weiteren als Abtastpixel bezeichnet) neuen Koordinatenwerten (2x, 2y) zugeordnet.
  • Die Pixel des rekonstruierten Bilds mit voller Größe an den Positionen (2x+1, 2y), die zwischen den beiden Abtastpixeln (2x, 2y) und (2x+2, 2y) liegen, werden berechnet über:
  • Figure 00080001
  • Die Pixel mit den Koordinaten (2x, 2y+1), die zwischen den beiden Abtastpixeln (2x, 2y) und (2x, 2y+2) liegen, werden berechnet über:
  • Figure 00080002
  • Die Pixel mit den Koordinaten (2x+1, 2y+1) zwischen den vier Abtastpixeln (2x, 2y), (2x+2, 2y), (2x, 2y+2) und (2x+2, 2y+2) werden berechnet über:
  • Figure 00080003
  • Der (durch die Gleichungen 6–8 beschriebene) letztere Vorgang wird hier als die gradientenbasierte Interpolation (GBI) bezeichnet.
  • Das pyramidenförmige Rekonstruktionsverfahren ist in 4 dargestellt.
  • Kontrastverstärkungsschritt
  • Zuerst wird die mehrskalige Gradientendarstellung W →0ƒ,..., W →K–1ƒ, SKƒ (oder W →'0ƒ, ..., W →'K–1ƒ, S'Kƒ im Fall der Pyramide) des Bilds f berechnet.
  • Um den Kontrast im Bild zu verstärken, sollen die Transformationskoeffizienten derart modifiziert werden, daß Details mit kleiner Amplitude auf Kosten derjenigen mit größerer Amplitude vergrößert werden, und zwar dies gleichförmig bei allen Skalen. Deshalb müssen die Transformationskoeffizienten auf nichtlineare, monotone Weise modifiziert werden.
  • Bei einer bevorzugten Ausführungsform wird die folgende Kontrastverstärkungsfunktion verwendet:
    Figure 00090001
    m ist eine Obergrenze für die Transformationskoeffizienten, 0<p<1 bestimmt den Grad an Nichtlinearität und c verhindert, daß die Funktion y für kleine Werte von x zu groß wird. Dies geschieht, um eine unnötige Überverstärkung des Rauschens im Bild zu verhindern. Der Wert von c sollte deshalb zu dem Rauschgehalt des Bilds in Beziehung stehen, d.h., c wird als die geschätzte Standardabweichung des Rauschens in dem Bild genommen.
  • Die Funktion y ist kontinuierlich und nimmt monoton ab, weist einen konstanten Wert
    Figure 00090002
    für 0<x<c auf und erreicht den Wert 1 für x=m. (siehe 5).
  • Die modifizierte mehrskalige Darstellung V →0, ..., V →K–1, SKƒ (alternativ V →'0,...,V →'K–1,S'Kƒ im Fall einer pyramidenförmigen Zerlegung) ist definiert als
    Figure 00100001
    Für alle Pixel i und alle Skalen 0≤j<k
  • Dies hat den folgenden Effekt: Gradientenvektoren mit kleiner Größe werden um einen größeren Faktor multipliziert als Gradientenvektoren mit großer Größe, so daß Details mit kleiner Amplitude auf Kosten von auffallenderen bevorzugt werden.
  • Das verstärkte Ausgangsbild f wird dann berechnet, indem die IMG an V →0,..., V →K–1, SKƒ (oder V →0,...,V →K–1,S'Kƒ) durchgeführt wird.
  • Gegebenenfalls werden die Grauwerte dann wieder linear auf den Originalbereich skaliert.
  • Man beachte, daß bei dieser Modifikation nur die Gradientengrößen beeinflußt werden, nicht ihre Richtungen. Dies geschieht, um bei Rekonstruktion Artefakte zu verhindern. Dieses Ziel wird jedoch nur teilweise erreicht, weil das Verfahren die Transformationskoeffizienten so behandelt, als wenn sie voneinander unabhängig wären.
  • Dies stimmt jedoch nicht ganz, da der Gradient an der j-ten Skala tatsächlich eine geglättete (möglicherweise unterabgetastete) Version des Gradienten an der j-1-ten Skala ist. Wenn die Gradienten auf die obige Weise unabhängig modifiziert werden, geht letztere Eigenschaft verloren und sogar ist V →0,..., V →K–1, SKƒ nicht länger eine MGR.
  • Das bedeutet, daß kein Bild g derart existiert, daß V →0,..., V →K–t, SKƒ die MGR von g ist.
  • Das gleiche gilt natürlich im Fall einer pyramidenförmigen Zerlegung V →'0,..., V →'K–1, S'Kƒ.
  • Infolge dessen sind die Gradienten des rekonstruierten Bilds f nicht länger parallel zu denen des Vorlagenbilds f, wie in Betracht gezogen wurde.
  • Die Verzerrung der Randgradientenorientierung aufgrund einer unabhängigen Modifikation von MGR-Pixeln bei nachfolgenden Skalen ist ein unerwünschter Effekt und kann zu sichtbaren Artefakten führen.
  • Artefaktkorrekturschritt
  • Das obige Problem wird gelöst, wenn die Rekonstruktion auf einen Satz von Gradientenbildern X →0,X →1,...,X →K–1 anstelle von V →0,...,V →K–1 (oder V →'0,...,V →'K–1) angewendet wird, die den folgenden Einschränkungen genügen
    • D1. Es existiert ein Bild g derart, daß X →0,X →1,...,X →K–1,SKg die MGR von g ist
    • D2. für alle X →j(i)∥⁣W →jƒ(i)(∥⁣W →jƒ(i)) für alle 0≤j<k und alle Pixel i
    • D3. X →j ist mehr entzerrt als W →jƒ(W →'jƒ).
  • In D3 ist der Ausdruck „mehr entzerrt" in dem Sinne von Gleichung (10) zu verstehen: X →j wird als mehr entzerrt als W →jƒ bezeichnet, wenn |W →jƒ(i)| „klein" ist, dann |X →j(i)|>|W →j(i)|, während wenn |W →jƒ(i)| „groß" ist, dann |X →j(i)|<|W →jƒ(i)|.
  • Mit anderen Worten versucht das Artefaktkorrekturverfahren, die Gradientenvektoren des verstärkten (D3) Bilds zu den Gradientenrichtungen des Vorlagenbilds umzulenken, und zwar auf eine Weise, die mit der mehrskaligen Gradientenrepräsentationsstruktur (D1) übereinstimmt.
  • Nachfolgend werden zwei alternative Ausführungsformen von Verfahren beschrieben, die einen derartigen Satz X →0,X →1,...,X →K–1 konstruieren.
  • - Abwechselnde konvexe Projektionen
  • Es ist wohlbekannt, daß alternatives Projizieren auf zwei konvexe nichtdisjunkte Teilsätze eines linearen Raums ein Konvergenzverfahren ist, das einen Punkt am Schnittpunkt dieser beiden Sätze liefert (siehe Youla D., Webb H., Image restoration by the method of convex projections, IEEE Trans. on Medical Imaging 1, 81–101 (1982)).
  • Mit diesem Wissen werden in einer ersten Ausführungsform Gradientenbilder konstruiert, die die obigen Anforderungen D1, D2 und D3 erfüllen.
  • Das Verfahren beginnt bei der modifizierten Darstellung V →0,..., V →K–1, SKƒ (oder von V →'0,...V →'K–1,S'Kƒ), definiert in (10), und iteriert eine abwechselnde Sequenz von Projektionsoperationen wie folgt in die beiden Vektorräume:
    • 1. Projektion von mehrskaligen Gradienten in den Raum W, d.h. alle Sätze von Vektorbildern {X →0,X →1,...,X →K–1} derart, daß ein Bild g existiert, für das X →j = W →jg (oder X →j = W →'jg) für 0≤j<K
    • 2. Projektion von Sätzen von Vektorbildern {X →0,X →1,...,X →K–1} in den Raum V|, derart, daß für alle Pixel i und alle Skalen 0≤j<K, X →j(i)'∥⁣W →jƒ(i) oder (X →'j(i)∥⁣W →jƒ(i)).
  • Das Ergebnis der ersten Projektion wird in die zweite Projektion eingespeist. Die resultierenden Vektorbilder werden wieder in die erste Projektionsstufe eingespeist, usw.
  • Die erste Projektion von mehrskaligen Gradienten auf den Raum W geschieht dadurch, daß an einem Satz von Vektorbildern der mehrskalige Gradientenrekonstruktionsvorgang durchgeführt wird (Gleichungen (3) oder (4), (5), (6), (7) und nach Anwenden des mehrskaligen Gradientenzerlegungsvorgangs (Gleichung (1) oder (4) an das resultierende Bild). Diese Projektionsoperation wird mit PW bezeichnet.
  • Die zweite Projektionsoperation auf den Vektorraum V| von Bildern X →'j, deren vektorwertige Pixelwerte parallel zu entsprechenden Gradientenpixeln der mehrskaligen Gradientendarstellung W →jƒ ist, erfolgt, indem die orthogonale Projektion auf W →jƒ(i)(W →ƒ(i) für jeden Wert von i und j genommen wird. Die orthogonale Projektion eines Vektors ν → auf dem Vektor w → ist gleich
  • Figure 00130001
  • Projektion auf V| wird mit Pv| bezeichnet.
  • Das abgegebene Bild erhält man durch Anwenden des IMG-Rekonstruktionsprozesses auf das Ergebnis des Iterierens der obigen Alternierung aus Projektionsoperationen PV| ∘PW...PV| ∘PW{V →0,...,V →K–1,SKƒ} (11)und analog im Fall des Pyramidenverfahrens.
  • Für eine Konvergenz werden etwa 10 Iterationen benötigt, d.h., die Operationen PW und Pv| müssen 10 mal durchgeführt werden. Dieser Algorithmus führt zu einer mehrskaligen Gradientendarstellung, die allen drei obigen Einschränkungen D1., D2. und D3. genügt.
  • - Mehrskalige Gradientenprojektion
  • Eine zweite Ausführungsform des Artefaktkorrekturprozesses wird als nächstes beschrieben, der auf der mehrskaligen Stapelausführungsform der MGR basiert; danach wird eine dritte Ausführungsform der Artefaktkorrektur auf der Basis der pyramidenförmigen Darstellung offenbart.
  • Die zweite Ausführungsform basiert auf der Tatsache, daß für jedes j, 0≤j<K–1, ein linearer Operator Fj auf den Raum von Vektorbildern derart existiert, daß für jedes Bild f gilt W →j+1ƒ = Fj(W →jƒ) (12)
  • Dies bedeutet, daß das mehrskalige Gradientenbild bei einer spezifischen Skala aus dem mehrskaligen Gradientenbild mit der vorausgegangenen kleineren Skala durch eine lineare Filteroperation Fj () berechnet werden kann.
  • Tatsächlich
  • Figure 00140001
  • Fj wirkt somit auf ein Vektorbild V → als Fj(V →) = bj·hj, hj·bj)·V → (14)
  • Das Ergebnis dieser Operation ist ebenfalls ein Vektorbild. Anhand (13) ist das Filter bj definiert durch die Beziehung bj·gj = gj+1·hj (15)
  • Somit können die Filterkoeffizienten von bj aus h0 und g0 berechnet werden. Dies geschieht auf die folgende Weise.
  • Filter werden oftmals durch ihre z-Transformierte dargestellt. Die z-Transformierte A(z) eines Filters a ist eine Funktion der durch A(z) = Σ akzk definierten komplexen Zahl z, wobei die Summe über alle Filterabgriffe k läuft, in denen ak die Koeffizienten des Filterkerns sind.
  • Mit unserer Auswahl von h0 und g0
    H0(z)=1/16(z–2+4z–1+6+4z+z2) und G0(z)=1/2(1–z)
  • Einer der Vorteile der z-Transformierten besteht darin, daß sie eine Faltung in ein Produkt umwandelt, d.h. für zwei Filter a und b ist die z-Transformierte von a·b normal A(z)B(z).
  • Deshalb kann bj leicht mit Hilfe seiner z-Transformierten berechnet werden.
  • Aus (15) Bj(z)=Gj+1(z)Hj(z)/Gj(z).
  • Dies ergibt B0(z)=1/16(z–3+5z–2+10z–1+10+5z+z2), weshalb
    b0=(1/16,5/16,10/16,10/16,5/16,1/16),
    b1=(1/16,0,5/16,0,10/16,0,10/16,0,5/16,0,1/16), usw.
  • Das Artefaktkorrekturverfahren gemäß einer zweiten Ausführungsform beginnt mit dem Zerlegen des verstärkten Bilds f in seine mehrskalige Gradientendarstellung, wodurch man U →0,...,U →K–1,SK erhält.
  • Gemäß dem Kriterium D2 müssen diese Vektoren in Richtung der mehrskaligen Gradienten des Vorlagenbilds umgelenkt werden.
  • Für j=0,...,K–1 soll Pj den pixelmäßigen orthogonalen Projektionsoperator auf W →jƒ bezeichnen. 1–Pj ist dann die pixelmäßige orthogonale Rejektion w.r.t. W →jƒ. Die orthogonale Rejektion ν →⊥ eines Vektors ν → w.r.t. der Vektor w → ist gegeben durch die Formel
  • Figure 00150001
  • Infolge des Projizierens von U →0 auf W →0ƒ gemäß P0 gehen einige Informationen verloren, und dies kann bei Rekonstruktion zu einem Kontrastverlust führen.
  • Diese verlorenen Informationen befinden sich in (1–P0)U →0. (1–P0)U →0 weist eine Nullkomponente entlang W →0ƒ auf.
  • Jedoch weist F0{((1–P0)U →0} keine Nichtnullkomponente entlang W →0ƒ auf. Außerdem kann F0{((1–P0)U →0} obwohl er formell zur Skala 1 gehört, als ein Beitrag zur Skala 0 angesehen werden, da er auf keinerlei Weise eine geglättete Version von P0U0 ist (er wird erhalten aus ((1–P0)U →)
  • Deshalb addieren wir (1/2)P0F0{1–P0)U →0} zu P0U →0.
  • Der Faktor 1/2 kommt deshalb ins Spiel, weil das in der Operation F0 verwendete Filter b0 eine Verstärkung (=Summe aller Filterabgriffe) gleich 2 aufweist, so daß durch das Addieren lediglich von P0F0{1–P0)U →0} ein Beitrag addiert würde, der um das Doppelte zu groß ist.
  • Indem dieser Vorgang wiederholt wird, wird das Vektorbild bei der kleinsten Skala X →0 definiert durch:
  • Figure 00160001
  • Dies wird der volle Beitrag zur Skala 0 sein.
  • Das Vektorbild bei der nächsten Skala Xi wird auf ähnliche Weise beginnend bei F0{X →0} erhalten
    Figure 00160002
    und so weiter und analog
  • Figure 00160003
  • Das ausgegebene Endbild wird dann rekonstruiert durch Anwenden des IMG auf {X →0,...,X →K–1,SK}
  • Gemäß unseren Erkenntnissen genügt {X →0,...,X →K–1,SK} den Einschränkungen D2. und D3. genau und D1. ungefähr.
  • Dieser Algorithmus ist schneller als der alternierende konvexe Projektionsalgorithmus und behält den in dem Kontrastverstärkungsschritt erhaltenen Kontrast besser bei.
  • Es wird nun angegeben, welche Änderungen vorgenommen werden sollten, wenn das pyramidenförmige Zerlegungsverfahren verwendet wird, gemäß einer dritten Ausführungsform der Artefaktkorrektur.
  • Es existiert weiterhin ein linearer Operator F'0 derart, daß für jedes Bild f W →'j+1ƒ = F'0(W →'jƒ) (20)aber in diesem Fall wirkt F'0 auf ein Vektorbild V → als F'0(V →) = ↓(b0⊗h0,h0⊗b0)·V → (21)
  • Das Filter b0 ist das gleiche wie oben definiert, doch beinhaltet der Operator F'0 auch die durch ↓ bezeichnete Unterabtastung. Die verworfenen Komponenten der Projektion (1–Pj) werden auf eine Weise ähnlich der zweiten Ausführungsform zurückaddiert, doch muß der Effekt der Unterabtastung durch ordnungsgemäße Interpolation umfaßt werden.
  • Bei der bevorzugten Ausführungsform wurde eine bilineare Interpolation gewählt, die mit dem Operator I bezeichnet ist. Dieser Operator wird auf jede der Beiträge bei größeren Skalen angewendet, um die Bildabmessung der nächst kleineren Skala anzupassen.
  • Die korrigierten Gradientenbilder X →j werden dann erhalten durch:
  • Figure 00170001
  • Wenn ein Gradientenfeld mit einer bestimmten größeren Skala auf ein Gradientenfeld mit einer kleineren Skala projiziert wird, müssen natürlich die entsprechenden Pixel identifiziert werden, d.h., wenn Y → k Skalen größer ist als Z →, dann sollte das Pixel (x, y) von Y → mit dem Pixel (x2k,y2k) von Z → identifiziert werden.

Claims (7)

  1. Verfahren zum Korrigieren von Artefakten in einem verarbeiteten digitalen Bild f, dargestellt durch eine mehrskalige Gradientendarstellung V={V →0,...,V →K–1,SKƒ}, umfassend Gradientenbilder bei mehreren Auflösungsniveaus und ein Restbild, wobei die mehrskalige Gradientendarstellung v erhalten wird durch eine nichtlineare Verarbeitung einer mehrskaligen Gradientendarstellung W={W →0ƒ,...,W →K–1ƒ,SKƒ} eines Vorlagenbilds f, gekennzeichnet durch die folgenden Schritte: (i) Erzeugen einer mehrskaligen Gradientendarstellung durch Modifizieren der Gradientenbilder von V, um die Richtungen ihrer Gradientenvektoren pixelmäßig gleich denen von W zu machen, und (ii) Rekonstruieren eines ausgegebenen Endbilds aus dieser modifizierten mehrskaligen Gradientendarstellung durch Anlegen einer inversen mehrskaligen Gradiententransformierten (IMG), wobei die inverse Transformierte derart ist, daß sie bei Anwendung an W f ergibt.
  2. Verfahren nach Anspruch 1, wobei – eine mehrskalige Gradientendarstellung X={X →0,X →1,...,X →K–1,SK},K>1 erzeugt wird und wobei das abgegebene Endbild aus der mehrskaligen Gradientendarstellung X={X →0,X →1,...,X →K–1,SK},K>1 rekonstruiert wird, wobei die mehrskalige Gradientendarstellung {X →0,X →1,...,X →K–1,SK},K>1 den folgenden Einschränkungen genügt (1) ein Bild g≠ƒ existiert derart, daß X →j für j=0,...,K-1,. das Gradientenbild von g bei der j-ten Skala darstellt, (2) in jedem Pixel weist X →j die gleiche Richtung wie das mehrskalige Gradientenbild W →jƒ des Ursprungsbilds f bei der j-ten Skala auf, (3) wenn das mehrskalige Gradientenvorlagebild W →j eine geringe Größe aufweist, dann ist die Größe von X →j größer als die Größe des Gradientenvorlagebilds; wenn das Gradientenvorlagebild eine große Größe aufweist, dann ist die Größe X →j kleiner als die Größe des Gradientenvorlagebilds bei jedem Pixel und bei jeder Skala j.
  3. Verfahren nach Anspruch 1 oder 2, wobei die mehrskalige Gradientendarstellung des Vorlagenbilds erhalten wird durch (i) Anwenden von eindimensionalen Gradientenfiltern unabhängig auf die Zeilen und Spalten einer digitalen Bilddarstellung des Vorlagenbilds, wodurch man die horizontale und vertikale Komponente des mehrskaligen Gradientenbilds bei der kleinsten Skala erhält, und Anwenden eines zweidimensionalen Tiefpaßfilters auf die digitale Darstellung des Vorlagenbilds, wodurch man eine Annäherung des Vorlagenbilds bei einer nächst größeren Skala erhält, (ii) Identifizieren der obigen Operationen als den ersten Schritt einer iterativen Schleife und Durchführen von zusätzlichen Iterationen unter Verwendung des Annäherungsbilds, das sich aus einer vorausgegangenen Iteration ergibt, anstelle des Vorlagenbilds, wobei die Filter auf die aktuelle Skala angepaßt sind, wodurch man Gradienten- und Annäherungsbilder des Vorlagenbilds bei den nachfolgend größeren Skalen erhält.
  4. Verfahren nach Anspruch 3, wobei die Gradientenfilter Filterkoeffizienten (0,5, –0,5) oder ein Mehrfaches davon haben und wobei die Tiefpaßfilterung erhalten wird durch Anwenden eines horizontalen eindimensionalen Tiefpaßfilters auf die Zeilen eines Bilds, gefolgt von einem vertikalen eindimensionalen Tiefpaßfilter, angewendet auf die Spalten des resultierenden Bilds des horizontalen Filters, wobei beide Tiefpaßfilter Koeffizienten aufweisen ausgewählt aus dem folgenden Satz: (0,5, 0,5), (0,25, 0,5, 0,25), (0,0625, 0,25, 0,375, 0,25, 0,0625).
  5. Verfahren nach Anspruch 2, wobei die Einschränkungen (1) bis (3) erfüllt werden durch (A) Erhalten einer anfänglichen Abschätzung der mehrskaligen Gradientendarstellung X durch: a) Anwenden einer inversen mehrskaligen Grandietentransformierten (IMG) auf die mehrskalige Darstellung des verarbeiteten Bilds f und danach Zerlegen des resultierenden Bilds in eine mehrskalige Gradientendarstellung X', b) Umlenken aller vektorwertigen Pixel von X', in (a) erhalten, durch eine pixelmäßige orthogonale Projektion auf die Pixel der mehrskaligen Gradientendarstellung W des Vorlagenbilds, wodurch man eine mehrskalige Gradientendarstellung X' erhält, so daß die Pixel von X' und W bei allen Skalen identische Richtungen aufweisen, (B) Identifizieren der obigen Operationen als dem ersten Schritt einer iterativen Schleife und Ausführen von zusätzlichen Iterationen unter Verwendung des Ergebnisses X' von Schritt (b) aus der vorausgegangenen Iteration als der Eingabe von Schritt (a) der aktuellen Iteration anstelle von V, wodurch man nachfolgende verbesserte Abschätzungen der mehrskaligen Gradientendarstellung X' erhält, (C) Erhalten einer mehrskaligen Gradientendarstellung X als das Ergebnis X'' von Schritt (2) nach der letzten Iteration.
  6. Verfahren nach Anspruch 2, wobei die Einschränkungen (1) bis (3) erfüllt werden durch: – Anwenden einer inversen mehrskaligen gradiententransformierten (IMG) auf die mehrskalige Darstellung V des verarbeiteten Bilds f und nachfolgendes Zerlegen des resultierenden Bilds in eine mehrskalige Gradientendarstellung U = {U →0,U →1,...,U →K–1,SK} – Erhalten von X →0 als die Summe
    Figure 00220001
    wobei Y →0k erhalten wird aus U →0 durch Anwenden der folgenden Schritte: Nehmen von U →0 als die erste Eingabe von A1 und A. Wiederholen für j=0 bis k: Extrahieren der Komponente des Eingabevektorbilds parallel zu W →0ƒ und der Komponente des zu W →0ƒ orthogonalen Eingabevektorbilds durch pixelmäßige orthogonale Projektion dieses Eingabevektorbilds auf W →0ƒ, wobei die Parallelkomponente das Ausgabevektorbild von A1 ist, wenn j<k, Anwenden eines Operators Fj auf die orthogonale Komponente, wobei der Operator Fj derart ist, daß man W →j+1ƒ bei Einwirkung auf W →jƒ erhält, Dividieren dieses Ergebnisses durch einen Faktor von 2, Verwendung dieses Ergebnisses als die nächste Eingabe zu Schritt A1, B. Nehmen der letzten Ausgabe von A1 entsprechend dem Schleifenindex j=k und gegebenenfalls Interpolieren dieses Ergebnisses, um die gleiche Anzahl von Pixeln wie W →0ƒ zu erhalten, wobei die Ausgabe von B Y →0k ist, – für i=1,..., K–1, Erhalten von X →i als die Summe
    Figure 00230001
    wobei Y →ik erhalten wird aus Fi–1(X →j–1) durch den folgenden Vorgang: Nehmen von Fi–1(X →j–1) als die erste Eingabe von C1 C. Wiederholen für j=i bis K: 1. Extrahieren der Komponente des Eingabevektorbilds parallel zu W →jƒ und der Komponente des Eingabevektorbilds orthogonal zu W →jƒ durch pixelmäßige orthogonale Projektion dieses Eingabevektorbilds auf W →jƒ, wobei die parallele Komponente die Ausgabe zu C1 ist, 2. falls j<k, Anwenden eines Operators Fj auf die orthogonale Komponente, wobei der Operator Fj derart ist, daß man W →j+1ƒ bei Einwirkung auf W →jƒ erhält, und Dividieren dieses Ergebnisses durch 2 unter Verwendung dieses resultierenden Vektorbilds als die nächste Eingabe von C1. D. Nehmen der letzten Ausgabe von C1 und gegebenenfalls Interpolieren dieses Ergebnisses, um die gleiche Anzahl von Pixeln wie W →jƒ zu erhalten.
  7. Verfahren nach Anspruch 1, wobei die nichtlineare Verarbeitung eine kontrastverstärkende Verarbeitung ist.
DE69830494T 1998-11-10 1998-11-10 Verfahren zur Verbesserung von Artefakten in digitalen Bildern Expired - Lifetime DE69830494T2 (de)

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
EP98203786A EP1001370B1 (de) 1998-11-10 1998-11-10 Verfahren zur Verbesserung von Artefakten in digitalen Bildern

Publications (2)

Publication Number Publication Date
DE69830494D1 DE69830494D1 (de) 2005-07-14
DE69830494T2 true DE69830494T2 (de) 2006-03-16

Family

ID=8234317

Family Applications (1)

Application Number Title Priority Date Filing Date
DE69830494T Expired - Lifetime DE69830494T2 (de) 1998-11-10 1998-11-10 Verfahren zur Verbesserung von Artefakten in digitalen Bildern

Country Status (3)

Country Link
EP (1) EP1001370B1 (de)
JP (1) JP2000172839A (de)
DE (1) DE69830494T2 (de)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6956975B2 (en) 2001-04-02 2005-10-18 Eastman Kodak Company Method for improving breast cancer diagnosis using mountain-view and contrast-enhancement presentation of mammography
EP1349113A3 (de) * 2002-03-20 2006-01-11 Ricoh Company Bildverarbeitungsvorrichtung und -verfahren
EP2017786A1 (de) * 2007-07-20 2009-01-21 Agfa HealthCare NV Verfahren zur Generierung eines kontrastverstärkten, mehrskaligen Bildes
EP2026278A1 (de) 2007-08-06 2009-02-18 Agfa HealthCare NV Verfahren zur Erhöhung des Kontrasts in einem Bild

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP0712092A1 (de) * 1994-11-10 1996-05-15 Agfa-Gevaert N.V. Verfahren zur Bildverbesserung

Also Published As

Publication number Publication date
JP2000172839A (ja) 2000-06-23
EP1001370A1 (de) 2000-05-17
DE69830494D1 (de) 2005-07-14
EP1001370B1 (de) 2005-06-08

Similar Documents

Publication Publication Date Title
DE102004007637B4 (de) Verfahren zum Erzeugen eines Bildes mit erhöhter Auflösung unter Verwendung einer Mehrzahl von Bildern mit niedriger Auflösung
DE69917708T2 (de) Blutgefässabbildung mit adaptiver Mittelwertbildung
DE602004003845T2 (de) Bildverarbeitungsvorrichtung zur Reduktion von Pixelrauschen
DE60014338T2 (de) Verfahren und gerät zur bildverarbeitung
DE3420576C2 (de) Anordnung zum Reprojizieren von Bildern aus mehreren eindimensionalen Projektionen in der Computer-Tomographie
DE69922983T2 (de) Bildgebungssystem und-verfahren
DE602005004410T2 (de) System und verfahren zur korrektur zeitlicher artefakte in tomographischen bildern
DE102019208496B4 (de) Computerimplementierte Verfahren und Vorrichtungen zum Bereitstellen eines Differenzbilddatensatzes eines Untersuchungsvolumens und zum Bereitstellen einer trainierten Generatorfunktion
DE102006005803A1 (de) Verfahren zur Rauschreduktion in bildgebenden Verfahren
DE19827034A1 (de) Iteratives Filtersystem für medizinische Bilder
DE102007013570A1 (de) Verfahren zur Rauschverminderung in digitalen Bildern mit lokal unterschiedlichem und gerichtetem Rauschen
DE102010029281A1 (de) Verfahren und Bildrekonstruktionseinrichtung zur Rekonstruktion von Bilddaten
DE19635017A1 (de) Verfahren und Vorrichtung zur Erhöhung der Bildschärfe
DE4224568C2 (de) Vorrichtung und Verfahren zur Bildung der Anzeige eines dreidimensionalen sequentiellen tomografischen Flächenschattierungsbildes
DE102010034099B4 (de) Iterative Bildfilterung mit anisotropem Rauschmodell für ein CT-Bild
DE69627755T2 (de) Verfahren und Vorrichtung zur temporären Rauschfilterung einer Bildfolge
DE102006017097A1 (de) Verfahren und System zur Behandlung des Rauschens in einem erzeugten Bild
EP1302899A2 (de) Vorrichtung und Verfahren zur Verarbeitung von Digitalbildern
DE102015117893A1 (de) Strahlbildungsvorrichtung, Ultraschallbildgebungsvorrichtung und Strahlbildungsverfahren
DE10356174A1 (de) Verfahren und Einrichtung zur Tomosynthese-Bildverbesserung unter Verwendung von Querfilterung
US6788826B1 (en) Method for correcting artefacts in a digital image
DE10305221B4 (de) Verfahren zur Reduktion von Rauschstrukturen in Arrays von Pixelwerten
DE69830494T2 (de) Verfahren zur Verbesserung von Artefakten in digitalen Bildern
DE102005045179A1 (de) Bildverbindung auf der Grundlage unabhängiger Rauschbedingungen
DE102012217613A1 (de) Verfahren zur Reduzierung von Artefakten bei der Rekonstruktion von Bilddatensätzen aus Projektionsbildern

Legal Events

Date Code Title Description
8364 No opposition during term of opposition
8320 Willingness to grant licences declared (paragraph 23)
8327 Change in the person/name/address of the patent owner

Owner name: AGFA HEALTHCARE NV, MORTSEL, BE