JP2013127773A - デジタルx線フレームシリーズにおけるノイズ低減方法 - Google Patents
デジタルx線フレームシリーズにおけるノイズ低減方法 Download PDFInfo
- Publication number
- JP2013127773A JP2013127773A JP2012216792A JP2012216792A JP2013127773A JP 2013127773 A JP2013127773 A JP 2013127773A JP 2012216792 A JP2012216792 A JP 2012216792A JP 2012216792 A JP2012216792 A JP 2012216792A JP 2013127773 A JP2013127773 A JP 2013127773A
- Authority
- JP
- Japan
- Prior art keywords
- frame
- noise
- estimation
- motion
- compensation
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 59
- 230000009467 reduction Effects 0.000 title claims abstract description 17
- 230000033001 locomotion Effects 0.000 claims abstract description 109
- 230000000694 effects Effects 0.000 claims abstract description 29
- 238000012935 Averaging Methods 0.000 claims abstract description 19
- 238000004364 calculation method Methods 0.000 claims abstract description 13
- 238000012937 correction Methods 0.000 claims abstract description 5
- 238000001914 filtration Methods 0.000 claims description 27
- 230000001419 dependent effect Effects 0.000 claims description 14
- 238000002156 mixing Methods 0.000 claims description 5
- 230000000877 morphologic effect Effects 0.000 claims description 5
- 230000001939 inductive effect Effects 0.000 claims description 3
- 238000012545 processing Methods 0.000 abstract description 25
- 239000000203 mixture Substances 0.000 abstract description 4
- 239000013598 vector Substances 0.000 description 18
- 238000004422 calculation algorithm Methods 0.000 description 12
- 230000002123 temporal effect Effects 0.000 description 9
- 238000013459 approach Methods 0.000 description 5
- 238000001514 detection method Methods 0.000 description 5
- 230000009466 transformation Effects 0.000 description 5
- 238000000354 decomposition reaction Methods 0.000 description 4
- 230000002829 reductive effect Effects 0.000 description 4
- 230000007704 transition Effects 0.000 description 4
- 230000000903 blocking effect Effects 0.000 description 3
- 230000006872 improvement Effects 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 230000005855 radiation Effects 0.000 description 3
- 238000002583 angiography Methods 0.000 description 2
- 230000002146 bilateral effect Effects 0.000 description 2
- 238000012512 characterization method Methods 0.000 description 2
- 238000002059 diagnostic imaging Methods 0.000 description 2
- 230000003628 erosive effect Effects 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 238000002594 fluoroscopy Methods 0.000 description 2
- 238000009499 grossing Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 238000012552 review Methods 0.000 description 2
- 230000000087 stabilizing effect Effects 0.000 description 2
- 238000007619 statistical method Methods 0.000 description 2
- 231100000987 absorbed dose Toxicity 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 210000004204 blood vessel Anatomy 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 239000003795 chemical substances by application Substances 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000002586 coronary angiography Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 230000008030 elimination Effects 0.000 description 1
- 238000003379 elimination reaction Methods 0.000 description 1
- 238000011835 investigation Methods 0.000 description 1
- 238000012886 linear function Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000035790 physiological processes and functions Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000013139 quantization Methods 0.000 description 1
- 238000002601 radiography Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000011524 similarity measure Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000000844 transformation Methods 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
- 210000005166 vasculature Anatomy 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- 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/10016—Video; Image sequence
-
- 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
-
- 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/20021—Dividing image into blocks, subimages or windows
-
- 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/20036—Morphological image processing
-
- 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/20172—Image enhancement details
- G06T2207/20182—Noise reduction or smoothing in the temporal domain; Spatio-temporal filtering
-
- 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/20172—Image enhancement details
- G06T2207/20201—Motion blur correction
-
- 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/20212—Image combination
- G06T2207/20216—Image averaging
-
- 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/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30101—Blood vessel; Artery; Vein; Vascular
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
- Picture Signal Circuits (AREA)
Abstract
【課題】デジタルX線透視画像のノイズ低減方法を提供する。
【解決手段】動き推定ブロックは重複によって考慮され、動き補償の間、ブロックの重複は平坦な重み付けウィンドウで平均化される。フリッカ効果補償は、基準フレームおよび現在のフレームを、所与のフレームのLFフィルタ処理によって得られる平均値の適切なフレームで除算することによって実施される。フリッカ効果を補償するために、基準フレームの局所的な平均強度が現在のフレームの局所的な平均強度となる。基準フレームはその平均値の動き補償されたフレームで除算され、現在のフレームの平均値のフレームで乗算される。動き補償アーティファクト補正をもたらす帰納的平均化は、所与のフレームの画素および構造的類似性の計算に基づいて現在のフィルタ処理されたフレームをその当初のバージョンと混合することを伴う。
【選択図】図1
【解決手段】動き推定ブロックは重複によって考慮され、動き補償の間、ブロックの重複は平坦な重み付けウィンドウで平均化される。フリッカ効果補償は、基準フレームおよび現在のフレームを、所与のフレームのLFフィルタ処理によって得られる平均値の適切なフレームで除算することによって実施される。フリッカ効果を補償するために、基準フレームの局所的な平均強度が現在のフレームの局所的な平均強度となる。基準フレームはその平均値の動き補償されたフレームで除算され、現在のフレームの平均値のフレームで乗算される。動き補償アーティファクト補正をもたらす帰納的平均化は、所与のフレームの画素および構造的類似性の計算に基づいて現在のフィルタ処理されたフレームをその当初のバージョンと混合することを伴う。
【選択図】図1
Description
発明の詳細な説明
発明の分野
本発明はデジタルX線画像処理の領域に関し、X線を含む高エネルギ放射線の使用により得られたデジタルX線フレームの処理に関連する問題を解消するために使用することができる。具体的には、本発明は、デジタルX線フレームシリーズの画像ノイズ低減のために設計される。
発明の分野
本発明はデジタルX線画像処理の領域に関し、X線を含む高エネルギ放射線の使用により得られたデジタルX線フレームの処理に関連する問題を解消するために使用することができる。具体的には、本発明は、デジタルX線フレームシリーズの画像ノイズ低減のために設計される。
従来の技術水準
現在、臨床実務は、シャープネスの強調、解剖学的構造の強調、減法等の、デジタルX線フレームシリーズ処理のさまざまな方法を含む。デジタルX線フレームシリーズ品質改良の最も重要な方法は、画像ノイズ低減方法と考えられる。血管造影法プロシージャ中の医療機器の制御はリアルタイムに実行され、カテーテルを通して造影剤を注入することによって脈管系が検査される。患者/医療スタッフの放射線吸収線量を減少させるために、低用量が与えられると、低い信号対ノイズ(S/N比)比を有する医療画像の取得がもたらされる[Chan他、Image Sequence Filtering in Quantum-Limited Noise with Applications to Low-Dose Fluoroscopy, IEEE Trans. Medical Imaging, vol. 12, issue 3, 610-621, 1993]。そのため、デジタルX線撮影では、ノイズが診断情報を守ることができる状態でのフィルタ処理方法を開発し、高い解像度値、フレーム速度および相当なレベルの信号強度依存ノイズにてリアルタイムに動作するための課題について検討することが急務である。
現在、臨床実務は、シャープネスの強調、解剖学的構造の強調、減法等の、デジタルX線フレームシリーズ処理のさまざまな方法を含む。デジタルX線フレームシリーズ品質改良の最も重要な方法は、画像ノイズ低減方法と考えられる。血管造影法プロシージャ中の医療機器の制御はリアルタイムに実行され、カテーテルを通して造影剤を注入することによって脈管系が検査される。患者/医療スタッフの放射線吸収線量を減少させるために、低用量が与えられると、低い信号対ノイズ(S/N比)比を有する医療画像の取得がもたらされる[Chan他、Image Sequence Filtering in Quantum-Limited Noise with Applications to Low-Dose Fluoroscopy, IEEE Trans. Medical Imaging, vol. 12, issue 3, 610-621, 1993]。そのため、デジタルX線撮影では、ノイズが診断情報を守ることができる状態でのフィルタ処理方法を開発し、高い解像度値、フレーム速度および相当なレベルの信号強度依存ノイズにてリアルタイムに動作するための課題について検討することが急務である。
医療用デジタルフレームシリーズのフィルタ処理のためのアルゴリズムを開発する際に研究者が遭遇する問題がいくつかある。第1の問題は、そのレベルが信号強度に本質的に依存する画像ノイズの推定に関連する。第2の問題は、自動強度調節のシステムが使用されるときに起こる、またはX線管の不安定さによる、信号強度の考えられる急激な変化を補償する必要性である。第3の問題は、異なる種類の動作、つまり表またはX線管/検出器の動き、患者の移動または患者の体内における生理学的プロセス、血管造影剤の流れによるフレームシリーズ可変性の検討に関連する。最後に、第4の問題は、リアルタイムにフィルタ処理を実施するために、フィルタ品質と要件との均衡を観察することが必要であるという点である。この問題は、フレームの積層内の時空間平均、変換空間内の時間的フィルタ処理(ピラミッド形変換、ウェーブレット変換等)、組合わせ平均法等のフィルタ処理方法の選択も含む。これらの問題の各々について詳細に説明する。
信号依存ノイズ推定の問題について検討する。実際には、あらゆる技術水準のノイズフィルタ処理方法は、単一のデジタルフレームについては、したがってそれらのシリーズについては、ノイズスペクトル特性についての情報を使用する。比較的多数のノイズ低減方法は、ノイズ分散値をパラメータとして利用する。多くの著者の出版物は、ノイズ分散値がより正確に分かればノイズ低減アルゴリズムの品質がより高まることを実証している[Bosco他、Noise Reduction for Cfa Image Sensors Exploiting Hvs Behavior, Sensors 9(3), 1692-1713, 2009; Foi, Practical Denoising of Clipped or Overexposed Noisy Images, Proc. 16th European Signal Process. Conf., EUSIPCO, 2008]。
デジタルX線撮影におけるノイズの性質を検討する。検出器は、減衰させられた(つまり検査中の対象物を通過した)X線放射の強度を測定する。
吸収された光子のランダムな変動は光子ノイズと呼ばれる。現代の検出器では、光子ノイズが主なノイズ源である。
付加的なノイズ源は、読出しノイズ、熱的ノイズ、アンプノイズ、量子化ひずみ等の検出器ノイズを含む[Gino M., Noise, Noise, Noise]。前述のノイズ源の蓄積効果は、ガウス分布ランダム変数によってモデル化される[Bosco他、Noise Reduction for Cfa Image Sensors Exploiting Hvs Behavior, Sensors 9(3), 1692-1713, 2009; Foi他、Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754, 2008]。広く用いられているモデルによれば、線形電子回路の場合、光子ノイズと付加的なノイズ源との分散は真の信号値に直線的に依存する[Bernd Jaehne, Digital Image Processing, Springer 2005]。
σ2(I(p))=aI(p)+b (1)
ここでI(p)は、画像の画素p=(x,y)における信号強度レベルである。
ここでI(p)は、画像の画素p=(x,y)における信号強度レベルである。
したがって、デジタル画像ノイズは信号強度に直線的に依存する。多くの出版物が信号依存ノイズ推定の問題を扱っている[Bosco他、Noise Reduction for Cfa Image Sensors Exploiting Hvs Behavior, Sensors 9(3), 1692-1713, 2009; Foi他、Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754, 2008; Forstner, Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range images, In Springer Lecture Notes on Earth Sciences, 1998; Hensel他、Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006; Liu他、Automatic estimation and removal of noise from a single image, IEEE Transactions on Pattern Analysis and Machine Intelligence, 30(2), 299-314, 2008; Salmeri他、Signal-dependent Noise Characterization for Mammographic Images Denoising. IMEKO TC4 Symposium (IMEKOTC4 '08), Firenze, Italy, 2008; Waegli, Investigations into the Noise Characteristics of Digitized Aerial Images, In: Int. Arch. For Photogr. And Remote Sensing, Vol. 32-2, 1998]。[Hensel他、Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006]の著作には、デジタル医療フレームシリーズのフレームごとのノイズ推定の非パラメータ方法が記載されており、リアルタイムに動作するアルゴリズムに特に重点が置かれている。デジタルフレームノイズ推定への2パラメータ手法が[Foi他、Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754, 2008]の論文に記載されている。検出器から得られ、ガンマ補正などといった非線形変換を経ていない当初のデジタル画像のノイズは、信号ランダム変数への加算性によってモデル化される:
In(p)=I(p)+σ2(I(p))n(p) (2)
ここで、In(p)は観察中のノイズ画像の画素p中の信号強度レベルであり、σ2(I(p))はタイプ(1)の信号強度へのノイズ分散の依存度であり、n(p)∈N(0,1)は正常なランダム変数である。露出不足および露出過度、つまりダイナミックレンジの境界上の線形関数(1)の妨害を引起すセンサの非線形性を考慮に入れたノイズ分散のモデル曲線を作成する方法が示唆されている。
In(p)=I(p)+σ2(I(p))n(p) (2)
ここで、In(p)は観察中のノイズ画像の画素p中の信号強度レベルであり、σ2(I(p))はタイプ(1)の信号強度へのノイズ分散の依存度であり、n(p)∈N(0,1)は正常なランダム変数である。露出不足および露出過度、つまりダイナミックレンジの境界上の線形関数(1)の妨害を引起すセンサの非線形性を考慮に入れたノイズ分散のモデル曲線を作成する方法が示唆されている。
ダイナミックフィルタを開発する際の第2の問題は、フィルタ処理中のフレームシリーズの強度の変化、いわゆるフリッカを考慮に入れ、補償することが必要な点である。この問題は、ビデオフレーム処理の一連の異なる方法の第一段階であるフリッカ効果補償タスクに類似している。[Delon, Movie and Video Scale-Time Equalization. IEEE Transactions on Image Processing, 15(1):241-248, 2006; Pitie他、A New Robust Technique for Stabilizing Brightness Fluctuations in Image Sequence, In 2nd Workshop on Statistical Methods in Video Processing In conjunction with ECCV 2004. Springer, 2004; Wong, Improved Flicker Removal Through Motion Vectors Compensation, In Third International Conference on Image and Graphics, pp. 552-555, 2004]。自動強度調節システムまたはX線管の不安定さにより、デジタルX線フレームシリーズにフリッカが生じることがある。時間的画像シリーズ平均化においてフリッカ補償が欠如すると、フィルタ処理品質が低下する。したがって、フリッカ補償のない動きフレームにおける静止の単純な時間的平均化により、現在のフレームの推定値の相当な歪みがもたらされる。動きとブロック整合技術に基づくその補償とがある場合、フレームシーケンス強度における相当な変化により、動き推定における概算誤差と、したがってその不十分な補償とが生じる。
フリッカ推定および補償には、2つの主な局面がある。一方で、フリッカ効果は1シリーズ中の1フレームの視野に対する局所的性質を有していない。他方、それは動的過程のみにおいて、つまりフレームのシーケンスを考慮するときに生じる。これらの特徴は、フレームコントラスト/強度追従モーションにおけるいずれかの局所的な変化を識別することを伴う、フリッカ効果補償における主な問題を含むか、またはX線管強度変化によって引起こされる。[Delon, Movie and Video Scale-Time Equalization. IEEE Transactions on Image Processing, 15(1):241-248, 2006; Pitie他、A New Robust Technique for Stabilizing Brightness Fluctuations in Image Sequence, In 2nd Workshop on Statistical Methods in Video Processing In conjunction with ECCV 2004. Springer, 2004; Wong他、Improved Flicker Removal Through Motion Vectors Compensation, In Third International Conference on Image and Graphics, pp. 552-555, 2004]に広範に記載されているフリッカ効果モデリングの自然な方法は、以下の式に基づく。
I(p)=A(p)I0(p)+B(p) (3)
ここで、I0(p)は画素pの信号強度レベルであり、A(p),B(p)はフレームの視野において特定される滑らかで緩やかに変動する関数であり、I(p)は観察中のフレームの強度である。このモデルは、当初画像強度が観察中の画像強度に線形従属することを示唆するとともに、上記のフリッカ効果空間不均一性を伴う。相当数のフリッカ推定および補償方法はグリッドへのフレームセグメンテーションに基づき、グリッドは、重複するブロックと、考えられる線外値フィルタ処理によるA(p),B(p)関数のパラメータのブロックごとの推定とを含み、より正確な推定のために、モデルパラメータが動き推定と同時に算出される。
ここで、I0(p)は画素pの信号強度レベルであり、A(p),B(p)はフレームの視野において特定される滑らかで緩やかに変動する関数であり、I(p)は観察中のフレームの強度である。このモデルは、当初画像強度が観察中の画像強度に線形従属することを示唆するとともに、上記のフリッカ効果空間不均一性を伴う。相当数のフリッカ推定および補償方法はグリッドへのフレームセグメンテーションに基づき、グリッドは、重複するブロックと、考えられる線外値フィルタ処理によるA(p),B(p)関数のパラメータのブロックごとの推定とを含み、より正確な推定のために、モデルパラメータが動き推定と同時に算出される。
正確なノイズフィルタ処理方法を開発する際の第3の問題は、動きによって引起こされるフレームシーケンスの変化を考慮に入れる必要性を伴う。動きは、診断上重要ないずれかのデジタル医療用フレームシリーズの必須の属性である。非固定フレームを平均化すると、得られる画像において動きブレアーティファクトが生じる。実際には、フレームシリーズ時空間平均における動きの問題に対処するには2つの解決策がある。第1の種類の手法は、動きに関連した画像領域の検出と、これらの領域における時間的ノイズ低減の力の低下とに基づく。[Aufrichtig他、X-Ray Fluoroscopy Spatio-Temporal Filtering with Object Detection. IEEE Trans. Medical Imaging 14 (1995) 733-746; Hensel他、Motion and Noise Detection for Spatio-Temporal Filtering of Medical X-Ray Image Sequences. Biomed, Tech./Biomed. Eng. 50-1 (2005) 1106-1107; Konrad, Motion Detection and Estimation. In Bovik, A.C., ed.: Handbook of Image and Video Processing. 2nd ed. Elsevier Academic Press (2005) 253-274]。第2の種類の手法は、連続するフレーム間の画素の動き、したがって動き推定および補償を推定する[Brailean他、Noise Reduction Filters for Dynamic Image Sequences: A review. Proc. IEEE, vol. 83, pp. 1272-1292, 1995]。動きフレームにおける空電は、動き補償プロシージャの使用によって動き推定後に明示的に生成することができる。
第1の種類の最も単純なアルゴリズムは、リアルタイムでの実行のために実現可能な合理的な設備による。そのような方法では、動き検出は、たとえば、現在の平均化されたフレームと基準フレームとの差分をしきい値化することであり得る。ある現在のフレームの画素について、上記の差分が有限のしきい値より高い場合、時間的平均化は、たとえばその現在のフレームの空間において平均化することで置換される。しきい値は、当然ながらノイズ推定段階で得られたノイズレベルに相当することになる。そのような動き推定方法は、いくつかの欠点、インパルスノイズの視認性、ぎざぎざのあるエッジを有している。他方で、(たとえばブロック整合に基づいて)動き推定および補償に使用される第2の種類のアルゴリズムは、より高品質であると考えられる[Brailean他、Noise Reduction Filters for Dynamic Image Sequences: A review. Proc. IEEE, vol. 83, pp. 1272-1292, 1995]。しかし、それらの使用には障害がある。すべての種類の動きが有効に補償できるわけではなく[Bernd Jaehne, Digital Image Processing, Springer 2005]、また、高い計算要求により、これらの方法はリアルタイムで実現するのが困難である。
第4の問題は、フィルタ品質とリアルタイム実行要求との間に均衡をもたらす必要性を意味する。高周波および大きなフレームサイズ(たとえば30フレーム/秒、解像度512×512画素、15フレーム/秒、解像度1024×1024画素/フレーム)でリアルタイムに医療用フレームシリーズ処理を実施する必要性は、現代のX線診断システムにおいて高品質フィルタ処理結果をもたらす要求と矛盾する。ノイズ推定または動き推定であるそのようなプロシージャそのものが、現代のコンピュータ容量に左右される大きな計算要求を有する。
発明の開示
クレームに記載の発明は、医療用途のデジタルフレームシリーズの品質を高めることを目的とする。
クレームに記載の発明は、医療用途のデジタルフレームシリーズの品質を高めることを目的とする。
クレームに記載の発明の技術的成果は、ノイズ低減しやすい画質向上方法を呈示することである。
デジタルX線フレームの取得と、フレームごとの信号依存ノイズ分散推定と、動き推定および補償と、フリッカ効果推定および補償と、帰納的フレーム平均化とを含む、クレームに記載のデジタルフレームシリーズにおけるノイズ低減方法の技術的成果は、基準フレームにおける急激な変化と相関するノイズ画素値の形態的消去によるノイズ分散推定中に実現される。次いで、ノイズ分散の区間推定値のロバストな区分的線形近似を使用することによって、信号強度に対するノイズの依存度を表形式にし、基準デジタルフレームの画素ごとのノイズ分散推定の形態のノイズマップを決定し、動き推定および補償の際、ピラミッド方式のブロック整合が使用され、動き補償ブロック重複が平均化される際に推定ブロックが重複によって考慮され、フリッカ効果は、基準フレームおよび現在のフレームを、所与のフレームの線形LFフィルタ処理によって得られる平均値の対応フレームで除算することによって得られ、フリッカ効果を補償する際、局所的な平均基準フレーム強度は局所的な平均現在フレーム強度となり、したがって基準フレームはその平均値の動き補償されたフレームで除算され、現在のフレームの平均値のフレームで乗算され、帰納的平均化の際、動き補償アーティファクト補正は、それらの画素および構造的類似性の算出に基づいて、現在のフィルタ処理されたフレームをその当初のバージョンと混合することによって実施される。
クレームに記載の方法の考えられる実施例バージョンでは、フレームのエッジに対応する画素の算出された推定値に対する影響を減少させる(エッジの除去)ためのノイズ分散推定において、現在のフレームと動き補償された基準フレームとの差分が使用される。
クレームに記載の、デジタル医療用X線画像のノイズ推定の方法に最も近いものは、[Hensel他、Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006]の論文に記載されている変形例であり、それによると、当初の画像に基づいて信号依存ノイズ推定を作成するために、
・当初の画像のLFフィルタ処理によって真の信号を推定し、当初の画像と、ノイズ画像取得をもたらすその推定値との差分を算出し、
・当初の画像における急激な変化(エッジ、単一の「ホット」画素)に対応するノイズ画像画素値をいずれかの方法で除去し、
・真の画像推定値の強度範囲を区間に分割し、そのような各区間について、真の画像推定値の画素に対応するノイズフレーム画素値を蓄積し、
・この間隔で蓄積されたノイズ画像画素値を用いて、あらゆる区間のノイズ分散を計算する段階が必要である。
・当初の画像のLFフィルタ処理によって真の信号を推定し、当初の画像と、ノイズ画像取得をもたらすその推定値との差分を算出し、
・当初の画像における急激な変化(エッジ、単一の「ホット」画素)に対応するノイズ画像画素値をいずれかの方法で除去し、
・真の画像推定値の強度範囲を区間に分割し、そのような各区間について、真の画像推定値の画素に対応するノイズフレーム画素値を蓄積し、
・この間隔で蓄積されたノイズ画像画素値を用いて、あらゆる区間のノイズ分散を計算する段階が必要である。
本発明は、上記の論文において言及されているノイズ推定原理を用いる。主な相違は、シリーズの各フレームについて以下の動作が行なわれる点にある。
・当初の画像における急激な変化に関連したノイズフレームの画素値は、当初の画像中のエッジに対応するノイズフレームの画素値の形態的抽出によって除外される。
・信号強度に対するノイズの依存度を表す表関数をもたらすノイズ分散の区間推定のロバストな局所的線形近似が実行される。
・真の画像推定および計算された表関数に基づいて、当初の画像ノイズ分散の画素ごとの推定である画像としてノイズマップを算出する。
・当初の画像における急激な変化に関連したノイズフレームの画素値は、当初の画像中のエッジに対応するノイズフレームの画素値の形態的抽出によって除外される。
・信号強度に対するノイズの依存度を表す表関数をもたらすノイズ分散の区間推定のロバストな局所的線形近似が実行される。
・真の画像推定および計算された表関数に基づいて、当初の画像ノイズ分散の画素ごとの推定である画像としてノイズマップを算出する。
フリッカ効果補償を考慮に入れるクレームに記載の方法に最も近いものは、[Wong他、Improved Flicker Removal Through Motion Vectors Compensation. In Third International Conference on Image and Graphics, pp. 552-555, 2004]の論文に記載されている変形例である。この方法では、フレームシリーズ強度値を乗法で歪ませる平滑関数によってフリッカがモデル化される。この関数推定は、ブロック整合に基づいた動き推定と相互に実施するよう提案され、そのため区分的定数歪み関数が使用される。現在のフレームのすべての動き推定ブロックにおけるフリッカ効果値は、動き補償の結果として得られるこのブロックの平均値と基準フレームの適切なブロックとの関係に等しいと考えられる。得られたフリッカ効果推定値の区分的定常場は、平滑要求を満たすためにしきい値化および平滑化を受ける。
クレームに記載の発明では、フリッカ効果は平滑乗法関数としてモデル化もされる。しかし、フリッカ効果補償は、前フレームの局所的手段を現在のフレームの局所的手段に変形することにより実施され、そのために、前に補償されたフレーム(つまり基準フレーム)がその平均値の動き補償フレームで除算され、次いで現在のフレーム平均値のフレームで乗算される。フレーム平均値はウィンドウを用いて移動平均させることにより得られ、そのアパーチャは動き推定ブロックサイズと等しい。基準画像のノイズマップは、適切な方法で修正される。
クレームに記載の発明では、動きの推定および補償はフレームシリーズにおけるブロック整合に基づく[Wang他、Fast Algorithms for the Estimation of Motion Vectors. IEEE Trans. Image Processing, vol. 8, no. 3, pp. 435-438, March 1999]。動きベクトルの推定品質に対して大幅に影響するフリッカ効果の考察は、動き推定段階において現在のフレームおよび基準フレームがそれらの平均値フレームで除算される最も単純な手法、つまり一種の輝度正規化によって実施される。フリッカフレーム強度の考察に加え、動きの推定は、ピラミッド方式を用いることによってなされる。その理由は、医療用フレームシリーズのリアルタイムのフィルタ処理中は、速い動きの推定および補償(たとえば、検査対象物に密接に接近する際の心臓検査−冠動脈造影法中)が特に困難であるという事実にある。その結果、特に低いフレーム速度での動き補償に成功するためには、適切なブロック整合領域が大きくなりすぎる。したがって、速い動き補償中の計算の複雑度が高まるという問題を解決するために、本発明では動き推定にピラミッド形アルゴリズムが使用され、
・現在のフレームおよび動き補償されるべき基準フレームは、いくつかのレベルのピラミッド状のフレーム(検査種類に応じてせいぜい3レベル)に分解され、
・ブロック動きベクトル推定は、ピラミッドレベルの数を考慮してその半径が選択される小領域上を走査することによってLFピラミッド構成要素において行なわれ、
・ピラミッドレベルを最低から最高に変換する際、動き補償ブロックのサイズが変倍され、小半径の領域におけるサーチによってその動きベクトルが特定される。
・現在のフレームおよび動き補償されるべき基準フレームは、いくつかのレベルのピラミッド状のフレーム(検査種類に応じてせいぜい3レベル)に分解され、
・ブロック動きベクトル推定は、ピラミッドレベルの数を考慮してその半径が選択される小領域上を走査することによってLFピラミッド構成要素において行なわれ、
・ピラミッドレベルを最低から最高に変換する際、動き補償ブロックのサイズが変倍され、小半径の領域におけるサーチによってその動きベクトルが特定される。
ピラミッド方式の利用は、基準フレームにおける好適なブロック整合領域を大幅に拡張する。矩形ブロックに基づく動き補償に共通の過度なブロッキングを減少させるために、クレームに記載の発明は、重複ブロック補償を提案する[Orchard他、Overlapped Block Compensation: An Estimation-Theoretic Approach. IEEE Transactions On Image Processing, Vol. 3, No. 5, September 1994, pp. 693-699]。典型的な重複サイズは、ブロックの半径の4分の1に等しい。動き補償されたフレームを作成する際、重複したブロック領域に属する画素が平滑重み付け関数によって平均化される。
クレームに記載の発明の示差的特徴は、動き補償アルゴリズムのアーティファクトの改善の特別な方法の利用である。サーチ(整合)領域の半径が小さすぎるために概算動き補償誤差が発生する可能性があり、重複する対象物またはフレームコンテンツの急激な変化は、フィルタ処理されたフレームと、その空間で平均化されたバージョンとの混合によって減少する。混合重みは、画素および、現在のフレームの、動き補償されたフレームとの構造的類似性に基づいて計算される。画素類似性の尺度として、相互の距離が用いられる[Tomasi他、Bilateral Filtering for Gray and Color Images, In Proc. 6th Int. Conf. Computer Vision, New Delhi, India, 1998, pp. 839-846]。現在のフレームおよび補償されたフレームの画素ブロック間のユークリッド距離が、構造的類似性の計算に使用される[Buades他、Nonlocal image and movie denoising. Int. J. comput. Vision, 76(2):123-139, 2008]。画素および構造的類似性は、線形的に混合される。
効果的なノイズ低減をリアルタイムに提供することを目的としたクレームに記載の発明では、信号依存ノイズの自動推定、フレームシリーズの考えられる急激な強度変化の推定、動き推定および補償、動き補償アーティファクトの改善、ならびに時間的過剰平滑化によるカルマンフィルタが使用される。与えられたタスクを、共有PCのマルチコアプロセッサとnvidia(登録商標)のcuda(登録商標)対応のグラフィックスカードとに振り分けると、リアルタイムで高品質なフィルタ処理要求に適した結果を得ることに成功した。
発明の詳細な説明
添付の図1〜図6は、クレームに記載の技術的解決法、その考えられる実施例、および技術的成果の実現を例示する。
添付の図1〜図6は、クレームに記載の技術的解決法、その考えられる実施例、および技術的成果の実現を例示する。
X線フレームの取得は、たとえば、図1に表される装置によって実施される。X線ビーム2を発するX線管1を備える。X線ビーム2は検出器3に入射する。検出器3は、シンチレーションスクリーン(図示せず)と感光性アレイ(図示せず)とから成る。シンチレーションスクリーンは、感光性アレイの能動面と光学的に結合される。検出器3に入射したX線ビーム2は、シンチレーションスクリーンによって可視光線に変換され、検出器のセンサによってデジタル形式に変換される。
クレームに記載の発明において提案されているように、フィルタはいくつかの段階においてノイズを低減する。前述の段階の各々を詳細に検討する。
段階1。信号強度に依存するノイズ分散のフレームごとの推定。この段階では、現在のフレームのLF線形フィルタ処理によってフレーム推定が実施される[Hensel他、Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006]。速度の厳密な要件をもたらすために、リアルタイム用途では、最も単純な線形フィルタ処理(たとえば二項式フィルタ)を使用することが合理的である。得られた平滑フレームに基づいて、ノイズフレームが構築される。当初のフレームとフィルタ処理されたフレームとの差分。
最も単純なフィルタによる画像推定は非理想的である。エッジが過剰に平滑化される。当初のフレームと平滑化されたフレームとの差分を取ると、平滑領域におけるノイズ画素に加えて急激な変化(たとえば解剖学的構造エッジおよび非ノイズ画素)に対応する有限数の画素を含むノイズフレームがもたらされる。これらの画素は、推定中の分散値を大幅に歪ませることがあり、計算から除外されるものとする。そのため、原則として当初のフレームの平滑化された導関数をしきい値化することに存在し、しきい値がS/N局所的推定によって決定される様々な技術がある[Foi他、Practical poissonian-gaussian noise modeling and fitting for single-image raw-data. Image Processing, IEEE Transactions on 17, 1737-1754 (October 2008), Salmeri他、Signal-dependent Noise Characterization for Mammographic Images Denoising. IMEKO TC4 Symposium (IMEKOTC4 '08), Firenze, Italy, September 2008]。多くの細部を備えるフレーム領域では、そのような推定は不十分である。したがって、クレームに記載の発明では、エッジを除去する際、標準偏差の導関数および局所的な推定を計算することを必要としない、エッジを選択するためのより簡単な手法が提案される。非ノイズ画素を除去するこの形態的手法の本質は下記にある。
1. ノイズフレームは2つの値−正負の変化の二値フレームに分割される。
2. エッジに対応する領域を選択するために、得られたフレームに対して、その後の膨張による形態的侵食を行う。
3. 処理されたフレームを1つの二値フレーム−当初のフレームのエッジマップに一体化させる。
1. ノイズフレームは2つの値−正負の変化の二値フレームに分割される。
2. エッジに対応する領域を選択するために、得られたフレームに対して、その後の膨張による形態的侵食を行う。
3. 処理されたフレームを1つの二値フレーム−当初のフレームのエッジマップに一体化させる。
微細構造をよりよく守るため、侵食および膨張の実施は2×2の小さいマスク(2×2ウィンドウ)に基づく。
区間ノイズ分散推定値を計算する際、推定されたフレーム強度(強度範囲のエッジ)の最小値および最大値が決定され、サブ区間(たとえば32グレースケール度)が選択される。次いで、推定されたフレームの各画素について、所与の画素値を含む区間が求められ、適切なノイズフレーム画素値がこの間隔(エッジ画素は除く)においてノイズ分散推定値を計算するのに使用される。ノイズ分散の区間推定値の計算は、様々な式、たとえば中央絶対偏差を使用する一般的な不偏推定値またはロバストな推定値を使用することができる[Hensel他、Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006; Foi他、Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754, 2008]。このプロシージャは、信号強度へのノイズ分散の依存度を表す表関数をもたらす。
エッジマップを構築する際に生じる不正確さにより、区間分散推定値の計算中に概算誤差が生じる可能性がある。したがって、これらの推定値は、反復線外値除去を当該技術を用いてより正確に規定する[Hensel他、Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006]。各間隔について、この区間におけるその後のノイズ分散推定値の再計算による3つの標準的ノイズ分散に等しいしきい値の大きさを超えるノイズ画素値は反復的に除外される。
区間ノイズ分散推定値が計算されると、強度に対するノイズ分散の依存度が推定される。パラメータ的手法では、ノイズモデル(1)のパラメータ推定値をある方法(たとえば最小二乗方法、尤度関数最小化、指向的な最適化)で構築することができる。センサ線形性を考慮することもできる[Foi他、Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754, 2008]。しかし[Hensel他、Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006]に述べられているように、信号強度に対するノイズ分散の依存度を適切に示すパラメータモデルの構築は、いくつかの要因により障害を伴う可能性がある。これらの要因は、センサ非線形性、当初のフレームの非線形前処理(たとえば対数の発見)を伴う可能性がある。このため、区間分散推定値に基づいて強度に対するノイズの依存度の非パラメータ推定値に存在する、より実現可能で普遍的な手法が構築される。
クレームに記載の発明では、求められる依存度の構築への非パラメータ的手法が使用され、得られたノイズ分散の区間推定値に基づいて、補間表関数が作成される。与えられた表関数は、区間分散推定値のロバストな局所的線形近似に基づいて形成される。ロバストな方法の使用は付加的な減少する線外値の影響(区間分散推定値における概算誤差)をもたらす一方、近似局所性は、強度に対するノイズの依存度を示す曲線の複雑な傾向の反復をもたらす。したがって、当初のフレームの各強度値について得られた表関数はノイズ分散推定値を関連付ける。表の入口点は、たとえば、推定中のフレームの強度として機能することができる。
実際には、パラメータノイズ推定値の場合、当初のフレームのノイズ分散を安定させる変換が適用される手法を使用することができる[Starck他、Image Processing and Data Analysis: The Multiscale Approach. Cambridge University Press, 1998; Bernd Jaehne, Digital Image Processing, Springer 2005]。したがって、信号依存ノイズフィルタ処理の問題は、付加的な信号非依存ノイズ分散を減少させるという課題にある。クレームに記載の発明では、非パラメータ的ノイズ推定を実施しノイズマップを構築しつつ、以下の手法が提案される。推定された画像および補間表に基づいて、画像であるノイズマップが構築され、その各画素は、適切な当初の画像画素における平均平方ノイズ偏差を推定する。ノイズマップは、実用化に十分な精度の画素ごとのノイズ推定値をもたらす。ノイズ分散推定中に、フレームエッジに対応するフレーム画素の推定値の算出に対する影響を低下させる(エッジの除去)ために、現在の動き補償されたフレームと基準の動き補償されたフレームとの差分が使用される、クレームに記載の発明の実施例バージョンが可能である。
段階2。動きおよびフリッカ効果の推定および補償。クレームに記載の発明では、フレームシリーズ強度の考えられる可変性を考慮したブロック対応のサーチに基づいて動き推定が実施される。この目的のため、動きベクトルサーチプロシージャを使用する前に、現在のフレームおよび基準フレームがそれらの平均値のフレームで除算される。対応するシリーズフレームの平均値のフレームは、動き推定ブロックのサイズ、たとえば16×16画素に匹敵するフィルタ孔を備えた低域通過フィルタ処理によって得られる。重複がブロックサイズの1/4に等しいとき、動き推定ブロックが重複して与えられる。ブロック整合の際、動き推定ピラミッド形アルゴリズムが使用され、その間、
・現在のフレームおよび現在のフレームに対して動き補償されるべき基準フレームがいくつかのピラミッドレベルに分割され、
・ブロック動きベクトル推定がLFフレームピラミッドで実施され、
・下方から上方へのピラミッドレベルの移行中、ブロックサイズおよびそれらの重複が変倍され、小半径領域に対するサーチによって、所与のブロックの動きベクトルが特定され、
・動き補償されたフレームを構築する際、重複領域が平均化される。
・現在のフレームおよび現在のフレームに対して動き補償されるべき基準フレームがいくつかのピラミッドレベルに分割され、
・ブロック動きベクトル推定がLFフレームピラミッドで実施され、
・下方から上方へのピラミッドレベルの移行中、ブロックサイズおよびそれらの重複が変倍され、小半径領域に対するサーチによって、所与のブロックの動きベクトルが特定され、
・動き補償されたフレームを構築する際、重複領域が平均化される。
基準フレームから現在のフレームへの推定は、LFピラミッドのフレームにおいて実施される。せいぜい3つの分解レベルが使用される。原則として2つである。なぜなら、より多くの分解レベルを微細な物体(薄い容器)の動き推定に使用すると、動きベクトル推定値において概算誤差が存在し得るからである。下方から上方へのピラミッドレベルの移行中、ブロックサイズおよびそれらの重複が倍増され、小領域、たとえば3×3画素に対するサーチによって動きベクトルが特定される。ブロック類似性の尺度であるユークリッド平方メートル法を計算しつつ、早く算出を停止させる技術が使用される[Wang他、Fast Algorithms for the Estimation of Motion Vectors. IEEE Trans. Image Processing, vol. 8, no. 3, pp. 435-438, March 1999]。
動き推定後、補償された基準フレームが形成され、過度のブロッキングを減少させるために、ブロックの重複が平均化される。ブロックの重複を平均化するために、平滑重み付けウィンドウが使用される。さらに、フリッカ効果を補償するために、基準フレームの局所的手段が現在のフレームの局所的手段となる。
段階3。動き補償アーティファクト補正をもたらす帰納的な平均化。整合領域の半径が小さすぎることにより生じる動き補償概算誤差、重複する物体、またはフレームコンテンツの急激な変化は、クレームに記載の発明では、フィルタ処理されたフレームとその空間平均化されたバージョンとの混合によって減少される。混合重みは、現在のフレームおよび動き補償されたフレームの画素および構造的類似性に基づいて算出される。画素類似性の尺度として、相互距離が使用される[Tomasi他、Bilateral Filtering for Gray and Color Images, In Proc. 6th Int. Conf. Computer Vision, New Delhi, India, 1998, pp. 839-846]。構造的類似性の計算については、現在のフレームおよび補償されたフレームの画素ブロック間のユークリッド距離[Buades他、Nonlocal image and movie denoising. Int. J. comput. Vision, 76(2):123-139, 2008]。
Ir(p),Ic(p)を、画素p中の基準フレームおよび現在のフレームの強度値とする。
画素と構造的類似性とを組み合わせるフレームの、結果として得られる類似性の尺度の計算については、以下の式:
W(p)=λ・Ws(p)+(1-λ)・Wp(p) (5)
が使用される。
W(p)=λ・Ws(p)+(1-λ)・Wp(p) (5)
が使用される。
λは、適切な画素強度に基づいた重みに対して優勢な構造的重みを特定する。フレーム混合中、λ≒1(たとえば0.9)での混合された重み計算は、画素強度の単一インパルス値の影響を減少させることと、同時に(当初の現在の補償されていないフレームのコンテンツと同様の)正確な構造的フレームコンテンツを再構築することとをもたらす。
時間的平均化は、以下の式にしたがって、動き補償された基準フレーム(フレーム予測子)および現在のフレーム(観察されたフレーム)について一次元カルマンフィルタ処理の使用によって反復的に実施される。
発明の最良の実施例
区間ノイズ分散推定値を計算するために、フレーム強度Ie(x,y)のIminおよびImaxが求められ、強度範囲をピッチhM(hM=32)で区間Miに分割する。フレームIe(x,y)の各画素について、この画素を含む区間が求められ、Ne(x,y)の適切な画素値を用いて所与の区間Miにおけるノイズ分散推定値σ2(i)が算出され、エッジマップ値E(x,y)=1を有する画素値は除外される。
区間ノイズ分散推定値を計算する際、不偏分散推定値の式が適用される。
取得された区間分散推定値σ2(i)は、各区間Miについて、この区間におけるノイズ分散推定値の再算出による3つの標準ノイズ分散と等しいしきい値の尺度を超えるノイズ画素値が反復的に除外されるように特定される。
取得されたノイズ分散の区間推定値を使用した、強度に対するノイズ分散の依存度の推定について、補間表関数が構築される(図2の円)。この表関数は、区間分散推定値のロバストな線形近似に基づく。
各画素が当初のフレーム(図3)の適切な画素におけるノイズ分散を推定するフレームであるノイズマップは、推定されたフレームおよび補間表に基づいて構築される。
動き推定および補償は、ブロック整合およびフレームシリーズ強度の考えられる可変性を伴う。そのため、動き推定の前に、現在のフレームおよび基準フレームがそれらの平均値のフレームで除算される。平均値のフレームは、動き推定ブロックの寸法、たとえば16×16画素に匹敵するフィルタ孔を備えた低域通過フィルタ処理により取得される。動き推定ブロックは、重複がたとえば4に等しい場合には、重複によって与えられる。ブロック整合の際、動き推定ピラミッド形アルゴリズムが使用され、その間、
・現在のフレームおよび現在のフレームに対して動き補償されるべき基準フレームがいくつかのピラミッドレベルに分割され、
・ブロック動きベクトル推定がLF画像ピラミッドで実施され、
・下方から上方へのピラミッドレベルの移行中、ブロックサイズおよびそれらの重複が変倍され、小半径領域に対するサーチによって所与のブロックの動きベクトルが特定され、
・動き補償されたフレームを構築する際、重複領域が平均化される。
・現在のフレームおよび現在のフレームに対して動き補償されるべき基準フレームがいくつかのピラミッドレベルに分割され、
・ブロック動きベクトル推定がLF画像ピラミッドで実施され、
・下方から上方へのピラミッドレベルの移行中、ブロックサイズおよびそれらの重複が変倍され、小半径領域に対するサーチによって所与のブロックの動きベクトルが特定され、
・動き補償されたフレームを構築する際、重複領域が平均化される。
基準フレームから現在のフレームへの推定は、LFピラミッドのフレームにおいて実施される。せいぜい3つの分解レベルが使用される。原則として2つである。なぜなら、微細な物体(薄い容器)の動き推定のためにより多くの分解レベルを使用すると、動きベクトル推定において概算誤差が存在し得るからである。下方から上方へのピラミッドレベルの移行中、ブロックサイズおよびそれらの重複が倍増され、小領域たとえば3×3画素に対するサーチによって動きベクトルが特定される。ブロック類似性の尺度であるユークリッド平方メートル法を計算しつつ、早く算出を停止させる技術が使用される[Wang他、Fast Algorithms for the Estimation of Motion Vectors. IEEE Trans. Image Processing, vol. 8, no. 3, pp. 435-438, March 1999]。
動き推定後、補償された基準フレームが形成され、過度のブロッキングを減少させるために、ブロックの重複が平均化される。基準フレームノイズマップが同様の方法で補償される。ブロックの重複を平均化するために、平滑重み付けウィンドウ、たとえばコサインが使用される。
さらに、フリッカ効果を補償するために、基準フレームの局所的手段が現在のフレームの局所的手段となる。そのため、補償された基準フレームは、その動き補償された平均フレームで除算され、式4に従って現在のフレームの平均値のフレームで乗算される。基準フレームのノイズマップに対して、同様の平均化を行う。
現在のフレームの帰納的平均化。時間的平均は、式(6)〜(12)に従って動きおよび強度が補償されたフレームについて一次元(画素ごとの)カルマンフィルタ処理の使用によって帰納的に実施される。
動き補償誤差によって引起こされる時間的平均化アーティファクトを修正するための、フィルタ処理された正規化されたフレームと、その当初のバージョンとの混合(式(5)、(9)および(10))。現在のフレームと動きおよび強度が補償されたフレームとに基づいて混合重み(5)が算出され、構造的類似性の尺度として、現在のフレームおよび補償されたフレームの画素ブロック間のユークリッド距離が使用される。
構造的類似性の算出には、その寸法が選択されるたとえば11×11画素のスライディングブロックが与えられる。画素類似性の算出は、現在のフレームおよび基準の動き補償されたフレーム間の強度空間の距離として計算される。実際のシリーズにおける式(11)、(12)中のノイズ低減度rは、たとえばr=0.8が選択される。これは5分の1のノイズ低減に相当する。図4〜図6は、実際のシリーズ(血管造影法検査)へのフィルタの適用の結果を示す。図4は、シリーズの当初のフレームを示し、図5はフィルタ処理されたフレームを示し、図6は、図4と図5とのフレームの差分を示す。よりよく理解するため、当初のフレームおよびフィルタ処理されたフレームに対してシャープネス改善プロシージャを行う。
1 X線管、2 X線ビーム、3 検出器。
Claims (2)
- デジタルX線フレームシリーズのノイズ低減方法であって、デジタルX線フレームシリーズの取得と、信号依存ノイズ分散のフレームごとの推定と、動き推定および補償と、フリッカ効果推定および補償と、帰納的フレーム平均化とを含み、
信号強度依存ノイズ分散のフレームごとの推定は、当初のフレームのエッジに対応するノイズフレームの画素値の形態的消去を伴い、
信号強度に対するノイズの表形式依存度は、ノイズ分散の区間推定値のロバストな区分的線形近似の使用によって得られ、
ノイズマップは、当初のデジタルフレームのノイズ分散の画素ごとの推定として計算され、
動きおよびフリッカ効果の推定および補償の間、ブロック整合のためのピラミッド方式が利用され、
動き推定ブロックは重複によって考慮され、動き補償中に、ブロック重複が平滑重み付けウィンドウを適用することによって平均化され、
フリッカ効果補償は、基準フレームおよび現在のフレームを、所与のフレームのLFフィルタ処理によって得られる平均値の適切なフレームで除算することによって実施され、
フリッカ効果を補償するために、基準フレームの局所的な平均強度は現在のフレームの局所的な平均強度となることを特徴とし、
基準フレームは、その平均値の動き補償されたフレームで除算され、現在のフレームの平均値のフレームで乗算され、
動き補償アーティファクト補正による帰納的平均化は、所与のフレームの画素および構造的類似性の計算に基づいて、現在のフィルタ処理されたフレームをその当初のバージョンと混合することを伴う、方法。 - フレームのエッジに対応する計算された推定画像画素への影響を減少させることを目的とするノイズ分散推定は、現在のフレームと基準の動き補償されたフレームとの差分を要する、請求項1に記載の方法。
Applications Claiming Priority (2)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
EA201101502A EA017302B1 (ru) | 2011-10-07 | 2011-10-07 | Способ подавления шума серий цифровых рентгенограмм |
EA201101502 | 2011-10-07 |
Publications (1)
Publication Number | Publication Date |
---|---|
JP2013127773A true JP2013127773A (ja) | 2013-06-27 |
Family
ID=47136908
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2012216792A Pending JP2013127773A (ja) | 2011-10-07 | 2012-09-28 | デジタルx線フレームシリーズにおけるノイズ低減方法 |
Country Status (6)
Country | Link |
---|---|
US (1) | US20130089247A1 (ja) |
EP (1) | EP2579207A3 (ja) |
JP (1) | JP2013127773A (ja) |
KR (1) | KR20130038794A (ja) |
CN (1) | CN103177423A (ja) |
EA (1) | EA017302B1 (ja) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2015052890A1 (ja) * | 2013-10-10 | 2015-04-16 | キヤノン株式会社 | X線画像処理装置および方法、並びにx線撮像装置、プログラム |
GB2567668A (en) * | 2017-10-20 | 2019-04-24 | Nokia Technologies Oy | Deflickering of a series of images |
Families Citing this family (27)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
JP6071444B2 (ja) * | 2012-11-07 | 2017-02-01 | キヤノン株式会社 | 画像処理装置及びその作動方法、プログラム |
CN105144675B (zh) | 2013-03-15 | 2018-09-14 | 耶鲁大学 | 用于处理具有传感器相关噪声的成像数据的技术 |
KR102104403B1 (ko) * | 2013-05-28 | 2020-04-28 | 한화테크윈 주식회사 | 단일영상 내의 안개 제거 방법 및 장치 |
CN104240184B (zh) * | 2013-06-08 | 2017-09-26 | 通用电气公司 | 噪声标准差的估算方法和*** |
WO2015049103A1 (en) * | 2013-10-01 | 2015-04-09 | Agfa Healthcare | Method for noise reduction in an image sequence |
WO2015052351A1 (en) * | 2013-10-11 | 2015-04-16 | Mauna Kea Technologies | Method for characterizing images acquired through a video medical device |
CN104182948B (zh) * | 2013-12-23 | 2015-07-22 | 上海联影医疗科技有限公司 | 一种相关性噪声的估计方法 |
US10043243B2 (en) * | 2016-01-22 | 2018-08-07 | Siemens Healthcare Gmbh | Deep unfolding algorithm for efficient image denoising under varying noise conditions |
KR20240017978A (ko) * | 2016-01-29 | 2024-02-08 | 인튜어티브 서지컬 오퍼레이션즈 인코포레이티드 | 광 레벨 적응 필터 및 방법 |
CN106355559B (zh) * | 2016-08-29 | 2019-05-03 | 厦门美图之家科技有限公司 | 一种图像序列的去噪方法及装置 |
EP3438922B1 (en) * | 2017-07-27 | 2019-11-20 | Siemens Healthcare GmbH | Method for processing at least one x-ray image |
TWI680675B (zh) * | 2017-12-29 | 2019-12-21 | 聯發科技股份有限公司 | 影像處理裝置與相關的影像處理方法 |
US10891720B2 (en) | 2018-04-04 | 2021-01-12 | AlgoMedica, Inc. | Cross directional bilateral filter for CT radiation dose reduction |
US11080898B2 (en) | 2018-04-06 | 2021-08-03 | AlgoMedica, Inc. | Adaptive processing of medical images to reduce noise magnitude |
FR3083415B1 (fr) * | 2018-06-29 | 2020-09-04 | Electricite De France | Traitement d'un bruit impulsionnel dans une sequence video |
CN110730280B (zh) * | 2018-07-17 | 2021-08-31 | 瑞昱半导体股份有限公司 | 噪声等化方法与噪声去除方法 |
CN110008958B (zh) * | 2018-08-01 | 2021-04-16 | 金华市灵龙电器有限公司 | 直流无刷电机控制机构 |
IT201900001225A1 (it) * | 2019-01-28 | 2020-07-28 | General Medical Merate S P A | Metodo predittivo per controllare un'apparecchiatura radiologica e apparecchiatura radiologica che lo implementa |
KR102304124B1 (ko) * | 2019-02-21 | 2021-09-24 | 한국전자통신연구원 | 학습기반 3d 모델 생성 장치 및 방법 |
JP7403994B2 (ja) * | 2019-08-21 | 2023-12-25 | 富士フイルムヘルスケア株式会社 | 医用画像処理装置および医用画像処理方法 |
DE102019122667A1 (de) * | 2019-08-22 | 2021-02-25 | Schölly Fiberoptic GmbH | Verfahren zur Unterdrückung von Bildrauschen in einem Videobildstrom, sowie zugehöriges medizinisches Bildaufnahmesystem und Computerprogrammprodukt |
CN111242831B (zh) * | 2020-01-20 | 2022-11-08 | 暨南大学 | 一种基于Zernike矩抗几何攻击的可逆鲁棒水印方法 |
CN113709324A (zh) * | 2020-05-21 | 2021-11-26 | 武汉Tcl集团工业研究院有限公司 | 一种视频降噪方法、视频降噪装置及视频降噪终端 |
CN113160072B (zh) * | 2021-03-19 | 2023-04-07 | 聚融医疗科技(杭州)有限公司 | 一种基于图像金字塔的鲁棒自适应帧相关方法及*** |
CN115988995A (zh) * | 2021-06-18 | 2023-04-18 | 深透医疗公司 | 用于实时视频去噪的***和方法 |
CN115661006B (zh) * | 2022-12-29 | 2023-03-10 | 湖南国天电子科技有限公司 | 一种海底地貌图像去噪方法 |
CN116249018B (zh) * | 2023-05-11 | 2023-09-08 | 深圳比特微电子科技有限公司 | 图像的动态范围压缩方法、装置、电子设备和存储介质 |
Family Cites Families (15)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5923775A (en) * | 1996-04-04 | 1999-07-13 | Eastman Kodak Company | Apparatus and method for signal dependent noise estimation and reduction in digital images |
EP1294194B8 (en) * | 2001-09-10 | 2010-08-04 | Texas Instruments Incorporated | Apparatus and method for motion vector estimation |
US7187794B2 (en) * | 2001-10-18 | 2007-03-06 | Research Foundation Of State University Of New York | Noise treatment of low-dose computed tomography projections and images |
EP1392047B1 (en) * | 2002-08-22 | 2013-10-30 | Samsung Electronics Co., Ltd. | Digital document processing for image enhancement |
US8107535B2 (en) * | 2003-06-10 | 2012-01-31 | Rensselaer Polytechnic Institute (Rpi) | Method and apparatus for scalable motion vector coding |
US20070092007A1 (en) * | 2005-10-24 | 2007-04-26 | Mediatek Inc. | Methods and systems for video data processing employing frame/field region predictions in motion estimation |
RU2322694C2 (ru) * | 2006-03-09 | 2008-04-20 | Общество с ограниченной ответственностью "Комэксп" | Способ обработки изображений |
KR20080044737A (ko) * | 2006-11-16 | 2008-05-21 | 주식회사 메디슨 | 초음파 영상 처리 방법 |
RU2342701C1 (ru) * | 2007-08-15 | 2008-12-27 | Российская Федерация, от имени которой выступает Министерство обороны Российской Федерации | Способ комплексирования цифровых многоспектральных полутоновых изображений |
US8731062B2 (en) * | 2008-02-05 | 2014-05-20 | Ntt Docomo, Inc. | Noise and/or flicker reduction in video sequences using spatial and temporal processing |
WO2010096701A1 (en) * | 2009-02-20 | 2010-08-26 | Mayo Foundation For Medical Education And Research | Projection-space denoising with bilateral filtering in computed tomography |
EP2302588B1 (en) * | 2009-08-21 | 2012-02-29 | Telefonaktiebolaget L M Ericsson (Publ) | Method and apparatus for estimation of interframe motion fields |
US20110134315A1 (en) * | 2009-12-08 | 2011-06-09 | Avi Levy | Bi-Directional, Local and Global Motion Estimation Based Frame Rate Conversion |
CN102014240B (zh) * | 2010-12-01 | 2013-07-31 | 深圳市蓝韵实业有限公司 | 一种实时医学视频图像去噪方法 |
US8538114B2 (en) * | 2011-06-06 | 2013-09-17 | Kabushiki Kaisha Toshiba | Method and system utilizing parameter-less filter for substantially reducing streak and or noise in computer tomography (CT) images |
-
2011
- 2011-10-07 EA EA201101502A patent/EA017302B1/ru not_active IP Right Cessation
-
2012
- 2012-09-27 EP EP12186390.6A patent/EP2579207A3/en not_active Withdrawn
- 2012-09-28 CN CN2012103719525A patent/CN103177423A/zh active Pending
- 2012-09-28 JP JP2012216792A patent/JP2013127773A/ja active Pending
- 2012-10-05 US US13/645,518 patent/US20130089247A1/en not_active Abandoned
- 2012-10-05 KR KR1020120110554A patent/KR20130038794A/ko not_active Application Discontinuation
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2015052890A1 (ja) * | 2013-10-10 | 2015-04-16 | キヤノン株式会社 | X線画像処理装置および方法、並びにx線撮像装置、プログラム |
JP2015073800A (ja) * | 2013-10-10 | 2015-04-20 | キヤノン株式会社 | X線画像処理装置および方法、並びにx線撮像装置、プログラム |
US10016176B2 (en) | 2013-10-10 | 2018-07-10 | Canon Kabushiki Kaisha | X-ray image processing apparatus and method, and X-ray imaging apparatus |
GB2567668A (en) * | 2017-10-20 | 2019-04-24 | Nokia Technologies Oy | Deflickering of a series of images |
GB2567668B (en) * | 2017-10-20 | 2022-03-02 | Vivo Mobile Communication Co Ltd | Deflickering of a series of images |
Also Published As
Publication number | Publication date |
---|---|
US20130089247A1 (en) | 2013-04-11 |
EA201101502A1 (ru) | 2012-10-30 |
EP2579207A3 (en) | 2013-07-10 |
KR20130038794A (ko) | 2013-04-18 |
EP2579207A2 (en) | 2013-04-10 |
CN103177423A (zh) | 2013-06-26 |
EA017302B1 (ru) | 2012-11-30 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP2013127773A (ja) | デジタルx線フレームシリーズにおけるノイズ低減方法 | |
US8873825B2 (en) | Method of noise reduction in digital X-rayograms | |
JP6145889B2 (ja) | 放射線画像処理装置および方法並びにプログラム | |
JP6156847B2 (ja) | 放射線画像処理装置および方法並びにプログラム | |
EP3435326B1 (en) | Method for processing at least one x-ray image | |
Cerciello et al. | A comparison of denoising methods for X-ray fluoroscopic images | |
EP2477154B1 (en) | Noise assessment method for digital X-ray films | |
US9922409B2 (en) | Edge emphasis in processing images based on radiation images | |
CN111353958B (zh) | 图像处理方法、装置及*** | |
US9619893B2 (en) | Body motion detection device and method | |
EP2002648A1 (en) | Temperature artifact correction | |
JP2019126524A (ja) | 放射線画像処理装置、散乱線補正方法及びプログラム | |
Hosseinian et al. | Assessment of restoration methods of X-ray images with emphasis on medical photogrammetric usage | |
Horkaew et al. | Structural adaptive anisotropic NAS-RIF for biomedical image restoration | |
Nguyen | Digital Radiography with a Consumer Camera: Image Denoising and Deblurring | |
Castillo-Amor et al. | Reduction of blooming artifacts in cardiac CT images by blind deconvolution and anisotropic diffusion filtering | |
EP1001370A1 (en) | Method for correcting artefacts in a digital image | |
Gu et al. | Scatter correction by non-local techniques | |
Jiang et al. | A contrast enhancement filter for improving stent visibility in interventional X-ray fluoroscopy: an initial study | |
JP2020089599A (ja) | 画像処理装置、画像処理方法およびプログラム | |
Ruikar et al. | A Comparison of Filtering Techniques for Image Quality Improvement in Computed Tomography |