JP2013175182A

JP2013175182A - デジタルｘ線フレームのノイズ低減方法

Info

Publication number: JP2013175182A
Application number: JP2013031039A
Authority: JP
Inventors: Vasil Evich Merkur Ev Sergey; セルゲイ・バシリエビッチ・メルクリエフ
Original assignee: ZAKRYTOE AKCIONERNOE OBSHCHESTVO IMPUL'S
Current assignee: ZAKRYTOE AKCIONERNOE OBSHCHESTVO IMPUL'S
Priority date: 2012-02-22
Filing date: 2013-02-20
Publication date: 2013-09-05
Also published as: WO2013125984A3; US20130216117A1; EA201200255A1; CN103295190A; KR20130096658A; WO2013125984A2; EP2631870A1; US8873825B2

Abstract

【課題】デジタルＸ線フレームのノイズ低減方法を提供する。
【解決手段】本発明はデジタルＸ線画像処理の領域に関し、Ｘ線を含む高エネルギ放射線の使用により得られたデジタルＸ線フレームのノイズ除去に関連する問題を解消するために使用することができる。具体的には、本発明は、デジタルＸ線フレームシリーズの画像ノイズ低減のために設計される。本発明の課題は、デジタルＸ線フレームシリーズのノイズ除去のための多機能な方法を開発すること、特に、画像のテクスチャ領域内の局所的に方向付けられた対象物エッジに沿って向けられた、不連続性の形式の残留ノイズおよびアーティファクトのレベルを低減すること、および、残留ＬＦノイズのレベルを低減すること、および、小さな細部の過度のフィルタ処理（過剰平滑化）を排除することである。上記方法は、デジタルＸ線画像処理品質の向上をもたらす。
【選択図】図１

Description

発明の詳細な説明
発明の分野
本発明はデジタルＸ線画像処理の領域に関し、Ｘ線を含む高エネルギ放射線の使用により得られたデジタルＸ線フレームのノイズ除去に関連する問題を解消するために使用することができる。具体的には、本発明は、デジタルＸ線フレームシリーズの画像ノイズ低減のために設計される。

従来の技術水準
現在、臨床実務は、シャープネスの強化、解剖学的構造の強調、関心領域（region of interest：ＲＯＩ）自動検出等の、デジタルＸ線フレームシリーズ処理のさまざまな方法を含む。デジタルＸ線フレームシリーズ品質向上の最も重要な方法は、画像ノイズ除去方法と考えられる。画像ノイズ除去は、視覚認識および低コントラストの細部の視認性を向上させ、他の処理アルゴリズム演算の品質を向上させる。そのため、デジタルＸ線撮影では、ノイズが診断情報を守ることができる状態でのフィルタ処理方法を開発し、高い解像度値、フレーム速度および相当なレベルの信号強度依存ノイズにてリアルタイムに動作するための課題について検討することが急務である。

ボスコ他、Ｈｖｓ挙動を利用するＣｆａ画像センサのためのノイズ低減、センサーズ９（３）、ｐｐ．１６９２−１７１３、２００９

デジタルフレームシリーズのフィルタ処理のためのアルゴリズムを開発する際に研究者が遭遇する問題が２つある。第１の問題は、そのレベルが信号強度に本質的に依存する画像ノイズの評価に関連する。第２の問題は、フィルタ動作品質とフィルタ処理速度との一定のバランスをもたらすフィルタ処理方法の開発である。各々の問題について詳細に検討する。

信号依存ノイズ推定の問題について検討する。各々の優れたデジタル画像フィルタ処理方法において、ノイズスペクトルおよび振幅分布データが実際に適用される。比較的多数のノイズ除去方法は、ノイズ分散値をパラメータとして利用する。多くの著者の出版物は、デジタル画像のノイズが、ホワイトノイズの正規振幅およびスペクトル分布を有する定常ランダムプロセスと考えられるモデルを実証している。一時的なノイズ低減の問題、すなわち分散値が信号強度に依存するノイズについても調査中であり、ノイズ分散値がより正確に分かればノイズ除去アルゴリズムの結果がよりよくなることが明らかになった［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］。

デジタルＸ線画像ノイズ源について検討する。検出器は、減衰させられた（つまり検査中の対象物を通過した）Ｘ線放射の強度を測定する。

吸収された光子のランダムな変動は光子ノイズと呼ばれる。現代の検出器では、光子ノイズが主なノイズ源である。付加的なノイズ源は、読出しノイズ、熱的ノイズ、アンプノイズ、量子化ひずみ等の検出器ノイズを含む。前述のノイズ源の蓄積効果は、ガウス分布ランダム変数によってモデル化される［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, TIP-03364-2007-FINAL］。広く用いられているモデルによれば、線形電子回路の場合、光子ノイズと付加的なノイズ源との分散は真の信号値に直線的に依存する［Bernd Jaehne, Digital Image Processing, Springer 2005］。

ここで、ｕ（ｐ）は、画像の画素ｐ＝（ｉ，ｊ）における信号強度レベルであり、σ^２（ｕ（ｐ））は、信号強度依存によるノイズ分散値である。

したがって、デジタル画像ノイズは信号強度に直線的に依存する。多くの出版物が信号依存ノイズ推定の問題を扱っている［Bosco他、Noise Reduction for Cfa Image Sensors
Exploiting Hvs Behavior, Sensors 9(3), p. 1692-1713, 2009; Foi他、Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Image Processing, IEEE TRANSACTIONS, TIP-03364-2007-FINAL; Foerstner, 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, Proceedings BVM 2006, Hamburg, Germany, March 19-21, p.46-50, 2006; Liu他、Automatic estimation and removal of noise from a single image, IEEE Transactions on Pattern Analysis and Machine Intelligence, p.1-35, October, 2006; Salmeri他、Signal-dependent Noise Characterization for Mammographic Images Denoising. 16th IMEKO TC4 Symposium (IMEKOTC4 '08), Florence, 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］の著作には、デジタル医療フレームシリーズのフレームごとのノイズ推定の非パラメータ方法が記載されており、リアルタイムに動作するアルゴリズムに特に重点が置かれている。デジタルフレームノイズ推定への２パラメータ手法が［Foi他、Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754, 2008］の論文に記載されている。検出器から得られ、ガンマ補正などといった非線形変換を経ていないデジタル原画像のノイズは、信号ランダム変数への加算性によってモデル化される：

ここで、ｕ_ｎ（ｐ）は観察中のノイズ画像の画素中の信号強度レベルであり、σ^２（ｕ（ｐ））はタイプ（１）の信号強度へのノイズ分散の依存度であり、ｎ（ｐ）∈Ｎ（０，１）は通常のランダム変数である。露出不足および露出過度、つまりダイナミックレンジの境界上の線形関数の妨害を引起すセンサの非線形性を考慮に入れたノイズ分散のモデル曲線を作成する方法が示唆されている。

第２の問題は、フィルタ動作品質と速い動作速度とのバランスをもたらすことである。現代のＸ線診断複合体が高い解像度の画像を素早く処理するための要件は、高品質なフィルタ処理結果を提供したいという願望と矛盾する。現在でも、デジタル画像ノイズ除去の課題をカバーする非常に多様な数学的ツールがある。一般に、以下の２つのフィルタ処理方法が各々のノイズ除去方法をほぼカバーすることができる：
・変換領域内でのフィルタ処理
・原画像領域内でのフィルタ処理。

フィルタ処理方法の第１のグループ（変換領域内でのフィルタ処理）は、特定の変換を使用することによりソース画像を変換する。続いて、特定の規則に基づいて、適切な非線形演算子を用いて、たとえば係数の処理が行なわれる。このような処理の目的は、ノイズエネルギを低減し、恐らくは画像のシャープネスを向上させることである。再構築中の画像の逆変換がこのフィルタ処理プロセスを終了させる。このグループは、フーリエおよびウェーブレット変換：[Strang他、Wavelets and Filter Banks, Wellesley-Cambridge Press.; Vetterli他、Wavelets and Subband Coding, Prentice Hall, Englewood Cliffs, NJ, 1995]、複素ウェーブレット変換[Kingsbury, The dual-tree complex wavelet transform: A new technique for shift invariance and directional filters, In Proc. 8th IEEE DSP Workshop, Utah, Aug. 9-12, 1998, paper no. 86; Selesnick, The double-density dual-tree discrete wavelet transform. IEEE Trans. Signal Processing, vol. 52, no. 5, pp. 1304-1314, May 2004]、カーブレット変換[Candes他、Curvelets, multiresolution representation, and scaling laws, In SPIE conference on Signal and Image Processing: Wavelet Applications in Signal and Image Processing VIII, San Diego, 2000]等に基づく多くの方法を含む。所与のクラスフィルタの動作速度および動作品質は、使用される変換および係数処理技術によって規定される。

フィルタ処理方法の第２のグループ（原画像領域内でのフィルタ処理）も、異方性拡散[Weickert, Anisotropic Diffusion in Image Processing. Tuebner Stuttgart, 1998]、変分アルゴリズム[Rudin他、Nonlinear total variation based noise removal algorithms, Physica D, 60:259-268, 1992]、バイラテラルフィルタ処理[Tomasi他、Bilateral Filtering for Gray and Color Images, In Proc. 6th Int. Conf. Computer Vision, New Delhi, India, 1998, pp. 839-846]、非局所平均（ＮＬ平均）フィルタ[Buades他、On image denoising methods, SIAM Multiscale Modeling and Simulation, 4(2):490-530, 2005]、カーネル回帰[Takeda他、Kernel regression for image processing and reconstruction, IEEE Transactions on Image Processing, vol. 16, no. 2, pp. 349-366, February 2007]、最適空間適応フィルタ[Kervrann他、Optimal spatial adaptation for patch-based image denoising. IEEE Transactions on Image Processing, vol. 15, no. 10, October 2006]等の多数の手法を含む。原則として、これらの技術は、恐らくは（画像構造の局所的な方向などの）幾何学的データを考慮に入れた、画像自体の中の画素または画素のグループ（パッチ）間の類似性検索と、これらの画素強度の任意の修正（平均化）アルゴリズムの適用とを伴い、修正範囲は画素（パッチ）の類似性の度合いに依存する。

発明の開示
クレームに記載の発明では、２つのステップを含むデジタルＸ線画像ノイズ除去技術が提供される。第１のステップは、ノイズの信号強度依存推定を含む。第２のステップは、ノイズ除去と、結合バイラテラル非局所フィルタ処理、残留ＬＦノイズのピラミッド型除去および過剰平滑化補正等のいくつかの補助的なステップとを含む。

クレームに記載の、デジタル医療用Ｘ線画像のノイズ推定の方法に最も近いものは、［Hensel他、Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences, Springer, 2006］の論文に記載されている変形例であり、それによる
と、原画像に基づいて信号依存ノイズ推定を作成するために、
・原画像のＬＦフィルタ処理によって真の信号を推定し、原画像と、ノイズ画像取得をもたらすその推定値との差分を計算し、
・原画像における急激な変化（エッジ、単一の「ホット」画素）に対応するノイズ画像画素値をいずれかの方法で除去し、
・真の画像推定値の強度範囲を区間に分割し、そのような各区間について、真の画像推定値の画素に対応するノイズフレーム画素値を蓄積し、
・この間隔で蓄積されたノイズ画像画素値を用いて、あらゆる区間のノイズ分散を計算する段階が必要である。

本発明は、上記の論文において言及されているノイズ推定原理を用いる。主な相違は以下の通りである：１）原画像における急激な変化に関連したノイズフレームの画素値は、原画像中のエッジに対応するノイズフレームの画素値の形態的抽出によって除外される；２）信号強度に対するノイズの依存度を表わす表関数をもたらすノイズ分散の区間推定値のロバストな局所的線形近似が実行される；３）真の画像推定および計算された表関数に基づいて、原画像ノイズ分散の画素ごとの推定である画像としてノイズマップを演算する。

ノイズは、定数分散の加算であると想定される。すなわち、ノイズ画像のモデルは、各画素ｐについてσ^２（ｕ（ｐ））＝σ^２＝ｃｏｎｓｔである式（２）によって表わされる。ここで引用されるのはＮＬ平均の基本スキームであり、これは、ＮＬ平均の改良、促進および実際の適用のさまざまな局面に関して調査の大きな流れを引起した[Bilcu他、Fast non-local means for image de-noising, In Digital Photography III Proc. of SPIE IS&T, Russel A. Martin, Jeffrey M. DiCarlo, and Nitin Sampat, Eds., February 2007, vol. 6502; Dauwe他、A Fast Non-Local Image Denoising Algorithm, In Proc. SPIE Electronic Imaging, vol. 6812, San Jose, USA, Jan 2008; Mahmoudi他、Fast image and video denoising via nonlocal means of similar neighborhoods, IEEE Signal Processing Letters, 12(12):839-842, 2005; Peter他、Robust estimation approach for nonlocal-means denoising based on structurally similar patches, Int. J. Open Problems Compt. Math., vol. 2, no. 2, June 2009; Wang他、Fast non-local algorithm for image denoising, In Proc. of IEEE International Conference on Image Processing (ICIP), 2006, pp. 1429-1432]。

ＮＬ平均フィルタは演算することが非常に困難である。ＮＬ平均フィルタ動作を加速させる可能な方法の１つは、[Buades他、A review of image denoising algorithms, with a new one, Multiscale Modeling and Simulation (SIAM interdisciplinary journal), Vol. 4 (2), pp. 490-530, 2005]の論文に記載されており、そこで著者はベクトルＮＬ平均の使用を提案している。

上記のフィルタ処理基本スキーム（３）、（４）には以下の多くの欠点がある：
・画像のテクスチャ領域内の局所的に方向付けられた対象物エッジに沿って向けられた、不連続性の形式の高レベルの残留ノイズおよび特徴的なアーティファクト;
・残留ＬＦノイズの存在；
・小さな細部の過度のフィルタ処理（過剰平滑化）。

クレームに記載の発明の狙いは、特に、画像のテクスチャ領域内の局所的に方向付けられた対象物エッジに沿って向けられた、不連続性の形式の残留ノイズおよび特徴的なアーティファクトのレベルを低減する可能性があり、残留ＬＦノイズのレベルを低減する可能性があり、小さな細部の過度のフィルタ処理（過剰平滑化）を排除する可能性がある、機能性が拡大したＸ線画像ノイズ除去方法の開発である。

原画像の取得、信号強度に依存したノイズ分散の推定、結合バイラテラル非局所フィルタ処理、残留ＬＦノイズ低減、過剰平滑化補正を含む、クレームに記載の発明における技術的成果、すなわちデジタルＸ線画像処理品質の向上は、信号強度に依存したノイズ分散の推定の段階において、原画像中のエッジに対応するノイズ画像の形態的画素値抽出が行なわれ、区間分散推定値のロバストな局所的線形近似に基づいて信号強度に対するノイズ分散の依存度を表わす表関数が得られ、原画像ノイズ分散の画素ごとの推定である画像としてノイズマップが計算され、結合バイラテラル非局所フィルタ処理の段階において、ＮＬ平均フィルタおよびバイラテラルフィルタの重みを組合わせることによって平均化が行なわれ、残留ＬＦノイズ低減の段階において、バイラテラルフィルタによるＬＦフィルタ処理のピラミッドスキームが用いられ、最終的な画像が得られ、過剰平滑化補正は、フィルタ処理された画像および原画像の双方の構造的類似性および画素類似性に基づいて最終的な画像を原画像と混合することによって行なわれる、ことによって達成される。

クレームに記載の発明の考えられる実施例バージョンは、結合バイラテラル非局所フィルタ処理の段階において、線形ＬＦフィルタによって平滑化された画像について平均重みが計算される、というものである。

クレームに記載の発明の考えられる実施例バージョンは、結合バイラテラル非局所フィルタ処理の段階において、ＮＬ平均フィルタおよびバイラテラルフィルタの固定重み部の形式の結合重み関数が用いられ、これらのフィルタ処理方法は結合をもたらす、というものである。

クレームに記載の発明の考えられる実施例バージョンは、ＮＬ平均フィルタおよびバイラテラルフィルタの分離を伴う結合バイラテラル非局所フィルタ処理の段階において、バイラテラルフィルタの分離可能なバージョンが用いられる、というものである。

発明の詳細な説明
添付の図１〜図５は、クレームに記載の技術的解決法、その考えられる実施例、および技術的成果の実現を例示する。

Ｘ線画像取得のための装置を示す図である。頭部のＸ線画像を示す図である。図２に表わされる画像について、信号強度に対する通常のノイズ偏差の依存度を表わす表関数（５‐実線）を示す図である。図２に表わされる画像についての通常のノイズ偏差マップを示す図である。クレームに記載の方法の効果を示す図である。クレームに記載の方法の効果を示す図である。クレームに記載の方法の効果を示す図である。クレームに記載の方法の効果を示す図である。クレームに記載の方法の効果を示す図である。

Ｘ線フレームの取得は、たとえば、図１に表わされる装置によって実施される。装置は、Ｘ線ビーム２を照射するＸ線管１を備える。Ｘ線ビーム２は検出器３に入射する。検出器３は、シンチレーションスクリーン（図示せず）と感光性アレイ（図示せず）とから成る。シンチレーションスクリーンは、感光性アレイの活性面と光学的に結合される。検出器３に入射したＸ線ビーム２は、シンチレーションスクリーンによって可視光線に変換され、検出器のセンサによってデジタル形式に変換される。

クレームに記載の発明において提案されているように、フィルタはいくつかの段階においてノイズを低減する。前述の段階の各々を詳細に検討する。

段階１。信号強度に依存するノイズ分散の推定。この段階では、現在の画像のＬＦ線形フィルタ処理によって原画像推定が実施される［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 (October 2008), Salmeri他、Signal-dependent Noise Characterization for Mammographic Images Denoising. IMEKO TC4 Symposium (IMEKOTC4 '08), Firenze, Italy, September 2008］。多くの細部を備える画像領域では、そのような推定は不十分である。したがって、クレームに記載の発明では、エッジを除去する際、標準偏差の導関数および局所的な推定を計算することを必要としない、エッジを選択するためのより簡単な手法が提案される。非ノイズ画素を除去するこの形態的手法の本質は下記にある。
１．ノイズ画像は２つの値−正負の変化の二値画像に分割される。
２．領域に対応するエッジを選択するために、得られた画像に対して、その後の拡張による形態的侵食を行う。
３．処理された画像を１つの二値画像−原画像のエッジマップにまとめる。

微細構造をよりよく守るため、侵食および拡張の実施は２×２の小さいマスク（２×２ウィンドウ）に基づく。

区間ノイズ分散推定値を演算する際、推定された画像強度（強度範囲のエッジ）の最小値および最大値が決定され、サブ区間（たとえば３２グレースケール度）が選択される。次いで、推定された画像の各画素について、所与の画素値を含む区間が求められ、適切なノイズ画像画素値がこの間隔（エッジ画素は除く）においてノイズ分散推定値を演算するのに使用される。ノイズ分散の区間推定値の演算は、様々な式、たとえば中央絶対偏差を使用する一般的な不偏推定値またはロバストな推定値を使用することができる［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］。各間隔について、この区間におけるその後のノイズ分散推定値の再計算による３つの標準的ノイズ分散に等しいしきい値の大きさを超えるノイズ画素値は反復的に除外される。

区間ノイズ分散推定値が演算されると、強度に対するノイズ分散の依存度が推定される。パラメータ的手法では、ノイズモデル（１）のパラメータ推定値をある方法（たとえば最小二乗方法、尤度関数最小化、指向的な最適化）で構築することができる。センサ線形性を考慮することもできる［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］。したがって、信号依存ノイズフィルタ処理の問題は、付加的な信号非依存ノイズ分散を減少させるという課題にある。クレームに記載の発明では、非パラメータ的ノイズ推定を実施しかつノイズマップを構築しつつ、以下の手法が提案される。推定された画像および補間表に基づいて、画像であるノイズマップが構築され、その各画素は、適切な原画像画素における平均平方ノイズ偏差を推定する。

ノイズマップは、実用化に十分な精度の画素ごとのノイズ推定値をもたらす。図４は、図２に表わされる画像について構築されたノイズマップの例を示す。

この例は、露出過度（センサ飽和／フレア）の際のノイズ推定アルゴリズム演算、頭部の周りの光のハロー（図２の位置４において十分に見られる）を示す。

段階２。結合フィルタ処理。画素ブロックサイズが大きくなるにつれてフィルタ処理の標準的なスキーム（３）を使用すると、Ｓ／Ｎのピーク値において測定されるフィルタ処理品質は向上する。品質の向上に加えて、画素パッチサイズが大きくなるにつれてベクトルＮＬ平均（４）を使用すると、フィルタの速度性能も上がる。しかし、非常にテクスチャ化された領域において急な強度変化（エッジ）に沿った画素パッチサイズが大きくなると、局所的に方向付けられた対象物エッジに沿って向けられた、不連続性の形式の十分に高いレベルの残留ノイズおよび特徴的なアーティファクトが観察され得る。テクスチャ化された領域に残留ノイズおよび特徴的なアーティファクトが存在することは、類似のパッチ数の減少によって引起されるノイズ除去レベルの減少によって説明できる。高レベルの残留ノイズは、処理された画像の平滑な領域に見られる。クレームに記載の発明では、標準的な方法の上記の限界を軽減することを目的としたバイラテラル統合およびＮＬ平均統合に基づく結合ノイズ除去方法が開発された。この方法を、結合フィルタ処理アルゴリズムまたは結合フィルタと呼ぶ。結合非局所バイラテラルフィルタ処理は以下の原理に基づく。

結合非局所バイラテラルフィルタ処理の原理。平均化の際、バイラテラル（輝度）および非局所（構造的）の２つのマトリックスに基づいて、現在の画素への近さを用いて、すべての画素の重みを計算する。ＮＬ平均フィルタの場合、現在の画素と比較画素とが非局所平均フィルタの意味で類似していれば、すなわちそれらのブロックが構造的に類似していれば、構造的重みが優勢になる。一方、パッチが互いに類似しているのではなく現在の画素および比較画素の強度に類似していれば、バイラテラル重みの影響が増大する。残留ノイズを低減させるために、最も単純なローパスフィルタによる原画像のフィルタ処理によって得ることができる平滑化された画像にわたって平均重みを計算する。

ここで与えられる結合フィルタ処理スキーム（５）、（６）では、さまざまな方法によって決定され得る重み関数Ｖの選択が重要な役割を果たす。結合フィルタ処理の速度効果の高いアルゴリズムを開発するために、輝度および構造的重みの比率が所与である形式で重み関数を形成することにより、単純であるが品質の点で成功的な結合フィルタのバージョンを生成する可能性がある。

結合フィルタのベクトルバージョンでは、

である。

分離可能なバイラテラルフィルタを適用した場合の結合フィルタ処理の分割可能な形式（９）により、フィルタ処理速度がさらに向上して、ＮＬ平均フィルタのベクトル化されたバージョンの速度にフィルタ処理速度を本質的に近付けることができる。

重み（８）を備えた結合フィルタ処理方法（６）の効果が図５（Ａ，Ｂ，Ｃ，Ｄ，Ｅ）に示されている。図５ＡはＸ線試験パターンの断片を示し、図５Ｂは標準的なＮＬ平均フィルタのベクトル化されたバージョン（４）によって処理されたＸ線試験パターンの断片を示し、図５Ｃは結合フィルタによって処理されたＸ線試験パターンの断片を示し、図５Ｄは図５Ａと図５Ｂとの差を示し、図５Ｅは図５Ａと図５Ｃとの差を示す。結合フィルタを適用した結果、標準的なＮＬ平均フィルタを使用した時に現れる数字のエッジに沿って、目に見えるアーティファクトのノイズ除去が行なわれる。

段階３。ＬＦノイズ低減。定数分散のノイズ実現と同一強度の平滑な画像の領域内のバイラテラルフィルタおよびＮＬ平均フィルタは、最も単純な線形近隣平均化フィルタの特徴を有する（すなわち、移動平均の役割を果たす）。平均化の有限領域の適用により、ＬＦに対応するスペクトルノイズ成分の一部は、低減されることはなく、平滑な画像の領域において十分に見ることができる。

クレームに記載の発明では、線形回帰を用いたバイラテラルフィルタを使用して、画像のピラミッドにおいて残留ＬＦノイズの低減が行なわれる[Burt他、The Laplacian pyramid as a compact image code, IEEE Trans. Commun., vol. COM-31, pp. 532-540, Apr. 1983; Tomasi他、Bilateral Filtering for Gray and Color Images, In Proc. 6th Int. Conf. Computer Vision, New Delhi, India, 1998, pp. 839-846]。画像のピラミッドは以下のように構築される。

ここで、Ｉ（ｐ）は原画像であり、Ａ_ｋ（ｐ），Ｈ_ｋ（ｐ）は、画像ピラミッドのｋレベルのローパス（近似）画像および高周波数（詳細）画像であり、Ｒｅｄｕｃｅは、画像サンプリング周波数をフィルタ処理して２度低減させる（スケーリング）演算子であり、すなわちＲｅｄｕｃｅ（Ａ）＝（Ａ＊Ｌ_ｒ）↓２（Ｌ_ｒは線形ローパスフィルタであり、↓２はスケーリングである）であり、Ｅｘｐａｎｄは、クラウディング（crowding）サンプルおよび画像フィルタ処理演算子であり、すなわちＥｘｐａｎｄ（Ａ）＝（Ａ↑２）＊Ｌ_ｅ（↑２はクラウディングサンプルであり、Ｌ_ｅはローパスフィルタである）である。

ノイズ分散分布は、ローパスピラミッド成分（１０）にわたって推定される。そのため、標準偏差σを有するホワイトガウスノイズが結合フィルタの入力に付加される。近隣平均化フィルタとの平滑な画像の領域内のフィルタの上記の類似性のために、ピラミッドレベルの数字に対するローパスノイズ二乗平均平方根誤差依存度が得られる。

ここで、σ_０は、結合フィルタによって処理された画像ノイズの二乗平均平方根誤差レベルであり、σ_ｋは、ピラミッド（１０）の画像ノイズＡ_ｋ（ｐ）の二乗平均平方根誤差レベルであり、ｑ（ｋ，Ｌ_ｒ）は、画像Ａ_ｋ（ｐ）における二乗平均平方根誤差レベルの減少の係数である。

残留ノイズを低減させるために、ピラミッド（１０）の画像Ａ_ｋ（ｐ）はバイラテラルフィルタによってフィルタ処理される。ぎざぎざのアーティファクトを低減させるために[Buades他、The staircasing effect in neighborhood filters and its solution, IEEE Transactions on Image Processing, Vol. 15 (6), pp. 1499-1505, 2006]、線形回帰を用いたバイラテラルフィルタが使用される。

段階４。過剰平滑化補正。クレームに記載の発明では、有用なＨＦ情報の冗長な平滑化に対しては、フィルタ処理された画像とその元のバージョンとを混合することによって補正が行なわれる。フィルタ処理された画像およびその元のバージョンの画素類似性および構造的類似性について混合重みが演算される。構造的類似性を演算するために、画像画素パッチ間のユークリッド距離が演算される。画素類似性（すなわち、強度類似性）を演算するために、適切な画像の画素強度間の距離が計算される。

最終的な画素、および画像類似性が計算された画素の構造的尺度を演算するために、構造的類似性が以下の形式で組合わせられる。

そして、過剰平滑化補正は、以下の規則に従ってフィルタ処理された画像と元のノイズ画像との一体化を伴う。

補正（１７）による（１６）のような類似性尺度の定義により、画像の単一の急な変化（いわゆるバースト）の除去が可能になる。

発明の最良の実施例

ノイズ画像Ｎ_ｅ（ｐ）に基づく分散演算からエッジを除去する際、正および負の変化を有する画像が形成される。

エッジに対応する大きなエリアを選択するために、所与の二値画像に対して、２×２マスクを使用することによって、その後の拡張を伴う形態的侵食が行なわれる。

そして、これらの画像は、原画像エッジのマップを得るために一体化される。

区間ノイズ分散推定値を演算するために、画像Ｉ_ｅ（ｐ）の最小強度Ｉ_ｍｉｎおよび最大強度Ｉ_ｍａｘが得られ、ピッチｈ_Ｍが選択され、次いで強度範囲がピッチｈ_Ｍ（たとえば、ｈ_Ｍ＝３２）を有する区間Ｍ_ｉに分割される。各画像画素Ｉ_ｅ（ｐ）について、この画素値を含む区間が求められ、画像Ｎ_ｅ（ｐ）の適切な画素値を用いて所与の区間Ｍ_ｉにおけるノイズ分散推定値σ^２（ｉ）が計算され、エッジマップの値がＥ（ｐ）＝１である画素値は除外される。区間ノイズ分散推定値を計算する際、不偏推定値の式が用いられる。

得られた区間ノイズ分散推定値σ^２（ｉ）は、各区間Ｍ_ｉについて、この区間におけるノイズ分散推定値のその後の再計算によるノイズの３標準偏差と等しいしきい値の大きさを超えるノイズ画素値を反復的に除外するように特定される。

強度に対するノイズ分散の依存度を推定するために、得られたノイズ分散の区間推定値に基づいて補間表関数が作成される。この表関数は、ノイズ分散区間推定値のロバストな区分的線形近似を使用することによって作成される。そのため、区間Ｍ_ｉのグリッド上では、（区間Ｍ_ｉにおける多数の点ｎ_ｉに依存可能にされ得る）近似化ステップｈ_Ｉおよび半径ｒ_Ｉが選択される。表関数の以前のステップ中に得られた値の近似化のために以下の規則が使用される。

ここで、ｋは近似クノーテ（knote）の数であり、ｍ_ｋは区間Ｍ_ｉ（ｋ・ｈ_Ｉ）の中心であり、パラメータａ，ｂは大きさが最小の偏差に基づいて計算される[Press他、Numerical Recipes in C: The Art of Scientific Computing, Second Edition. New York: Cambridge University Press, 1999]。

評価画像および補間表に基づいてノイズマップが生成され、画像の各画素は原画像の適切な画素におけるノイズ分散を評価する。

フィルタ処理品質とフィルタ速度との優れた対応関係をもたらすために、式（９）を使用して結合バイラテラル非局所フィルタ処理が行なわれ、ここで、λ＝０．８５であり、平均化領域の半径はｔ＝５であり、画素ブロックの半径はｆ＝４であるように選択され、１つの画素の重複を有する状態でブロックが選択される。強度および構造について画素類似性を計算するために、連続的な遷移を有するしきい値関数Ｆ（ｒ，σ，ｋ，ｈ）＝ｅｘｐ（−ｍａｘ（ｒ^２−ｋ・σ^２，０）／ｈ^２）（ｋおよびｈはパラメータである）が使用される。

構造的類似性を計算するための過剰平滑化補正では、半径５のブロックが使用される。重み関数の形態を規定するパラメータは、結合フィルタ重み関数のパラメータと一致する。強度について画素類似性の重み関数を演算する際、幾何学的距離は考慮に入れられない。

１Ｘ線管、２Ｘ線ビーム、３検出器。
[参照文献]

Claims

原画像の取得と、信号依存ノイズ分散の推定と、結合バイラテラル非局所フィルタ処理と、残留ＬＦノイズの低減と、過剰平滑化補正とを含む、デジタルＸ線フレームのノイズ低減方法であって、
前記信号依存ノイズ分散の推定の段階は、原画像中のエッジに対応するノイズ画像の画素値の形態的抽出と、ノイズ分散の区間推定値のロバストな局所的線形近似を用いることによる、信号強度に対するノイズ依存度を表わす表関数の取得と、原画像ノイズ分散の画素ごとの推定である画像としてのノイズマップの演算とを伴い、
前記結合バイラテラル非局所フィルタ処理の段階は、ＮＬ平均フィルタ重みとバイラテラルフィルタ重みとを組合わせることによる平均化を伴い、
前記残留ＬＦノイズの低減の段階は、バイラテラルフィルタによるローパスフィルタ処理のピラミッドスキームを使用し、
最終的な画像が得られ、
フィルタ処理された画像と原画像との構造的類似性および画素類似性に基づいて最終的な画像と原画像とを組合わせることによって、過剰平滑化補正が行なわれることを特徴とする、方法。
前記結合バイラテラル非局所フィルタ処理の段階において、線形ローパスフィルタによって平滑化された原画像に基づいて平均重みが計算される、請求項１に記載の方法。
前記結合バイラテラル非局所フィルタ処理の段階において、ＮＬ平均フィルタ重みおよびバイラテラルフィルタ重みの比率が固定された形式の結合重み関数が用いられ、所与のフィルタ処理は最終的な統合を用いて分類する、請求項１に記載の方法。
バイラテラルフィルタからのＮＬ平均フィルタの分離を伴う前記結合バイラテラル非局所フィルタ処理の段階において、バイラテラルフィルタの分離可能なバージョンが用いられる、請求項１または２に記載の方法。