(57)【要約】【課題】デジタルＸ線フレームのノイズ低減方法を提供する。【解決手段】本発明はデジタルＸ線画像処理の領域に関し、Ｘ線を含む高エネルギ放射線の使用により得られたデジタルＸ線フレームのノイズ除去に関連する問題を解消するために使用することができる。具体的には、本発明は、デジタルＸ線フレームシリーズの画像ノイズ低減のために設計される。本発明の課題は、デジタルＸ線フレームシリーズのノイズ除去のための多機能な方法を開発すること、特に、画像のテクスチャ領域内の局所的に方向付けられた対象物エッジに沿って向けられた、不連続性の形式の残留ノイズおよびアーティファクトのレベルを低減すること、および、残留ＬＦノイズのレベルを低減すること、および、小さな細部の過度のフィルタ処理（過剰平滑化）を排除することである。上記方法は、デジタルＸ線画像処理品質の向上をもたらす。【選択図】図１
2013175182.tif 000002
【特許請求の範囲】
【請求項１】
原画像の取得と、信号依存ノイズ分散の推定と、結合バイラテラル非局所フィルタ処理と、残留ＬＦノイズの低減と、過剰平滑化補正とを含む、デジタルＸ線フレームのノイズ低減方法であって、
前記信号依存ノイズ分散の推定の段階は、原画像中のエッジに対応するノイズ画像の画素値の形態的抽出と、ノイズ分散の区間推定値のロバストな局所的線形近似を用いることによる、信号強度に対するノイズ依存度を表わす表関数の取得と、原画像ノイズ分散の画素ごとの推定である画像としてのノイズマップの演算とを伴い、
前記結合バイラテラル非局所フィルタ処理の段階は、ＮＬ平均フィルタ重みとバイラテラルフィルタ重みとを組合わせることによる平均化を伴い、
前記残留ＬＦノイズの低減の段階は、バイラテラルフィルタによるローパスフィルタ処理のピラミッドスキームを使用し、
最終的な画像が得られ、
フィルタ処理された画像と原画像との構造的類似性および画素類似性に基づいて最終的な画像と原画像とを組合わせることによって、過剰平滑化補正が行なわれることを特徴とする、方法。
【請求項２】
前記結合バイラテラル非局所フィルタ処理の段階において、線形ローパスフィルタによって平滑化された原画像に基づいて平均重みが計算される、請求項１に記載の方法。
【請求項３】
前記結合バイラテラル非局所フィルタ処理の段階において、ＮＬ平均フィルタ重みおよびバイラテラルフィルタ重みの比率が固定された形式の結合重み関数が用いられ、所与のフィルタ処理は最終的な統合を用いて分類する、請求項１に記載の方法。
【請求項４】
バイラテラルフィルタからのＮＬ平均フィルタの分離を伴う前記結合バイラテラル非局所フィルタ処理の段階において、バイラテラルフィルタの分離可能なバージョンが用いられる、請求項１または２に記載の方法。
【発明の詳細な説明】
【技術分野】
【０００１】
発明の詳細な説明
発明の分野
本発明はデジタルＸ線画像処理の領域に関し、Ｘ線を含む高エネルギ放射線の使用により得られたデジタルＸ線フレームのノイズ除去に関連する問題を解消するために使用することができる。具体的には、本発明は、デジタルＸ線フレームシリーズの画像ノイズ低減のために設計される。
【背景技術】
【０００２】
従来の技術水準
現在、臨床実務は、シャープネスの強化、解剖学的構造の強調、関心領域（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］。
【０００６】
デジタルＸ線画像ノイズ源について検討する。検出器は、減衰させられた（つまり検査中の対象物を通過した）Ｘ線放射の強度を測定する。
【０００７】
【数１】
2013175182.tif 000003
【０００８】
吸収された光子のランダムな変動は光子ノイズと呼ばれる。現代の検出器では、光子ノイズが主なノイズ源である。付加的なノイズ源は、読出しノイズ、熱的ノイズ、アンプノイズ、量子化ひずみ等の検出器ノイズを含む。前述のノイズ源の蓄積効果は、ガウス分布ランダム変数によってモデル化される［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］。
【０００９】
【数２】
2013175182.tif 000004
【００１０】
ここで、ｕ（ｐ）は、画像の画素ｐ＝（ｉ，ｊ）における信号強度レベルであり、σ２（ｕ（ｐ））は、信号強度依存によるノイズ分散値である。
【００１１】
したがって、デジタル画像ノイズは信号強度に直線的に依存する。多くの出版物が信号依存ノイズ推定の問題を扱っている［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］の論文に記載されている。検出器から得られ、ガンマ補正などといった非線形変換を経ていないデジタル原画像のノイズは、信号ランダム変数への加算性によってモデル化される：
【００１２】
【数３】
2013175182.tif 000005
【００１３】
ここで、ｕｎ（ｐ）は観察中のノイズ画像の画素中の信号強度レベルであり、σ２（ｕ（ｐ））はタイプ（１）の信号強度へのノイズ分散の依存度であり、ｎ（ｐ）∈Ｎ（０，１）は通常のランダム変数である。露出不足および露出過度、つまりダイナミックレンジの境界上の線形関数の妨害を引起すセンサの非線形性を考慮に入れたノイズ分散のモデル曲線を作成する方法が示唆されている。
【００１４】
第２の問題は、フィルタ動作品質と速い動作速度とのバランスをもたらすことである。現代のＸ線診断複合体が高い解像度の画像を素早く処理するための要件は、高品質なフィルタ処理結果を提供したいという願望と矛盾する。現在でも、デジタル画像ノイズ除去の課題をカバーする非常に多様な数学的ツールがある。一般に、以下の２つのフィルタ処理方法が各々のノイズ除去方法をほぼカバーすることができる：
・変換領域内でのフィルタ処理
・原画像領域内でのフィルタ処理。
【００１５】
フィルタ処理方法の第１のグループ（変換領域内でのフィルタ処理）は、特定の変換を使用することによりソース画像を変換する。続いて、特定の規則に基づいて、適切な非線形演算子を用いて、たとえば係数の処理が行なわれる。このような処理の目的は、ノイズエネルギを低減し、恐らくは画像のシャープネスを向上させることである。再構築中の画像の逆変換がこのフィルタ処理プロセスを終了させる。このグループは、フーリエおよびウェーブレット変換：[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］の論文に記載されている変形例であり、それによる
と、原画像に基づいて信号依存ノイズ推定を作成するために、
・原画像のＬＦフィルタ処理によって真の信号を推定し、原画像と、ノイズ画像取得をもたらすその推定値との差分を計算し、
・原画像における急激な変化（エッジ、単一の「ホット」画素）に対応するノイズ画像画素値をいずれかの方法で除去し、
・真の画像推定値の強度範囲を区間に分割し、そのような各区間について、真の画像推定値の画素に対応するノイズフレーム画素値を蓄積し、
・この間隔で蓄積されたノイズ画像画素値を用いて、あらゆる区間のノイズ分散を計算する段階が必要である。
【００１９】
本発明は、上記の論文において言及されているノイズ推定原理を用いる。主な相違は以下の通りである：１）原画像における急激な変化に関連したノイズフレームの画素値は、原画像中のエッジに対応するノイズフレームの画素値の形態的抽出によって除外される；２）信号強度に対するノイズの依存度を表わす表関数をもたらすノイズ分散の区間推定値のロバストな局所的線形近似が実行される；３）真の画像推定および計算された表関数に基づいて、原画像ノイズ分散の画素ごとの推定である画像としてノイズマップを演算する。
【００２０】
【数４】
2013175182.tif 000006
【００２１】
ノイズは、定数分散の加算であると想定される。すなわち、ノイズ画像のモデルは、各画素ｐについてσ２（ｕ（ｐ））＝σ２＝ｃｏｎｓｔである式（２）によって表わされる。ここで引用されるのはＮＬ平均の基本スキームであり、これは、ＮＬ平均の改良、促進および実際の適用のさまざまな局面に関して調査の大きな流れを引起した[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]の論文に記載されており、そこで著者はベクトルＮＬ平均の使用を提案している。
【００２３】
【数５】
2013175182.tif 000007
【００２４】
上記のフィルタ処理基本スキーム（３）、（４）には以下の多くの欠点がある：
・画像のテクスチャ領域内の局所的に方向付けられた対象物エッジに沿って向けられた、不連続性の形式の高レベルの残留ノイズおよび特徴的なアーティファクト;
・残留ＬＦノイズの存在；
・小さな細部の過度のフィルタ処理（過剰平滑化）。
【００２５】
クレームに記載の発明の狙いは、特に、画像のテクスチャ領域内の局所的に方向付けられた対象物エッジに沿って向けられた、不連続性の形式の残留ノイズおよび特徴的なアーティファクトのレベルを低減する可能性があり、残留ＬＦノイズのレベルを低減する可能性があり、小さな細部の過度のフィルタ処理（過剰平滑化）を排除する可能性がある、機能性が拡大したＸ線画像ノイズ除去方法の開発である。
【００２６】
原画像の取得、信号強度に依存したノイズ分散の推定、結合バイラテラル非局所フィルタ処理、残留ＬＦノイズ低減、過剰平滑化補正を含む、クレームに記載の発明における技術的成果、すなわちデジタルＸ線画像処理品質の向上は、信号強度に依存したノイズ分散の推定の段階において、原画像中のエッジに対応するノイズ画像の形態的画素値抽出が行なわれ、区間分散推定値のロバストな局所的線形近似に基づいて信号強度に対するノイズ分散の依存度を表わす表関数が得られ、原画像ノイズ分散の画素ごとの推定である画像としてノイズマップが計算され、結合バイラテラル非局所フィルタ処理の段階において、ＮＬ平均フィルタおよびバイラテラルフィルタの重みを組合わせることによって平均化が行なわれ、残留ＬＦノイズ低減の段階において、バイラテラルフィルタによるＬＦフィルタ処理のピラミッドスキームが用いられ、最終的な画像が得られ、過剰平滑化補正は、フィルタ処理された画像および原画像の双方の構造的類似性および画素類似性に基づいて最終的な画像を原画像と混合することによって行なわれる、ことによって達成される。
【００２７】
クレームに記載の発明の考えられる実施例バージョンは、結合バイラテラル非局所フィルタ処理の段階において、線形ＬＦフィルタによって平滑化された画像について平均重みが計算される、というものである。
【００２８】
クレームに記載の発明の考えられる実施例バージョンは、結合バイラテラル非局所フィルタ処理の段階において、ＮＬ平均フィルタおよびバイラテラルフィルタの固定重み部の形式の結合重み関数が用いられ、これらのフィルタ処理方法は結合をもたらす、というものである。
【００２９】
クレームに記載の発明の考えられる実施例バージョンは、ＮＬ平均フィルタおよびバイラテラルフィルタの分離を伴う結合バイラテラル非局所フィルタ処理の段階において、バイラテラルフィルタの分離可能なバージョンが用いられる、というものである。
【００３０】
発明の詳細な説明
添付の図１〜図５は、クレームに記載の技術的解決法、その考えられる実施例、および技術的成果の実現を例示する。
【図面の簡単な説明】
【００３１】
【図１】Ｘ線画像取得のための装置を示す図である。
【図２】頭部のＸ線画像を示す図である。
【図３】図２に表わされる画像について、信号強度に対する通常のノイズ偏差の依存度を表わす表関数（５‐実線）を示す図である。
【図４】図２に表わされる画像についての通常のノイズ偏差マップを示す図である。
【図５Ａ】クレームに記載の方法の効果を示す図である。
【図５Ｂ】クレームに記載の方法の効果を示す図である。
【図５Ｃ】クレームに記載の方法の効果を示す図である。
【図５Ｄ】クレームに記載の方法の効果を示す図である。
【図５Ｅ】クレームに記載の方法の効果を示す図である。
【発明を実施するための形態】
【００３２】
Ｘ線フレームの取得は、たとえば、図１に表わされる装置によって実施される。装置は、Ｘ線ビーム２を照射するＸ線管１を備える。Ｘ線ビーム２は検出器３に入射する。検出器３は、シンチレーションスクリーン（図示せず）と感光性アレイ（図示せず）とから成る。シンチレーションスクリーンは、感光性アレイの活性面と光学的に結合される。検出器３に入射したＸ線ビーム２は、シンチレーションスクリーンによって可視光線に変換され、検出器のセンサによってデジタル形式に変換される。
【００３３】
クレームに記載の発明において提案されているように、フィルタはいくつかの段階においてノイズを低減する。前述の段階の各々を詳細に検討する。
【００３４】
段階１。信号強度に依存するノイズ分散の推定。この段階では、現在の画像のＬＦ線形フィルタ処理によって原画像推定が実施される［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］。したがって、信号依存ノイズフィルタ処理の問題は、付加的な信号非依存ノイズ分散を減少させるという課題にある。クレームに記載の発明では、非パラメータ的ノイズ推定を実施しかつノイズマップを構築しつつ、以下の手法が提案される。推定された画像および補間表に基づいて、画像であるノイズマップが構築され、その各画素は、適切な原画像画素における平均平方ノイズ偏差を推定する。
【００４３】
ノイズマップは、実用化に十分な精度の画素ごとのノイズ推定値をもたらす。図４は、図２に表わされる画像について構築されたノイズマップの例を示す。
【００４４】
この例は、露出過度（センサ飽和／フレア）の際のノイズ推定アルゴリズム演算、頭部の周りの光のハロー（図２の位置４において十分に見られる）を示す。
【００４５】
段階２。結合フィルタ処理。画素ブロックサイズが大きくなるにつれてフィルタ処理の標準的なスキーム（３）を使用すると、Ｓ／Ｎのピーク値において測定されるフィルタ処理品質は向上する。品質の向上に加えて、画素パッチサイズが大きくなるにつれてベクトルＮＬ平均（４）を使用すると、フィルタの速度性能も上がる。しかし、非常にテクスチャ化された領域において急な強度変化（エッジ）に沿った画素パッチサイズが大きくなると、局所的に方向付けられた対象物エッジに沿って向けられた、不連続性の形式の十分に高いレベルの残留ノイズおよび特徴的なアーティファクトが観察され得る。テクスチャ化された領域に残留ノイズおよび特徴的なアーティファクトが存在することは、類似のパッチ数の減少によって引起されるノイズ除去レベルの減少によって説明できる。高レベルの残留ノイズは、処理された画像の平滑な領域に見られる。クレームに記載の発明では、標準的な方法の上記の限界を軽減することを目的としたバイラテラル統合およびＮＬ平均統合に基づく結合ノイズ除去方法が開発された。この方法を、結合フィルタ処理アルゴリズムまたは結合フィルタと呼ぶ。結合非局所バイラテラルフィルタ処理は以下の原理に基づく。
【００４６】
結合非局所バイラテラルフィルタ処理の原理。平均化の際、バイラテラル（輝度）および非局所（構造的）の２つのマトリックスに基づいて、現在の画素への近さを用いて、すべての画素の重みを計算する。ＮＬ平均フィルタの場合、現在の画素と比較画素とが非局所平均フィルタの意味で類似していれば、すなわちそれらのブロックが構造的に類似していれば、構造的重みが優勢になる。一方、パッチが互いに類似しているのではなく現在の画素および比較画素の強度に類似していれば、バイラテラル重みの影響が増大する。残留ノイズを低減させるために、最も単純なローパスフィルタによる原画像のフィルタ処理によって得ることができる平滑化された画像にわたって平均重みを計算する。
【００４７】
【数６】
2013175182.tif 000008
【００４８】
ここで与えられる結合フィルタ処理スキーム（５）、（６）では、さまざまな方法によって決定され得る重み関数Ｖの選択が重要な役割を果たす。結合フィルタ処理の速度効果の高いアルゴリズムを開発するために、輝度および構造的重みの比率が所与である形式で重み関数を形成することにより、単純であるが品質の点で成功的な結合フィルタのバージョンを生成する可能性がある。
【００４９】
【数７】
2013175182.tif 000009
【００５０】
結合フィルタのベクトルバージョンでは、
【００５１】
【数８】
2013175182.tif 000010
【００５２】
である。
【００５３】
【数９】
2013175182.tif 000011
【００５４】
分離可能なバイラテラルフィルタを適用した場合の結合フィルタ処理の分割可能な形式（９）により、フィルタ処理速度がさらに向上して、ＮＬ平均フィルタのベクトル化されたバージョンの速度にフィルタ処理速度を本質的に近付けることができる。
【００５５】
重み（８）を備えた結合フィルタ処理方法（６）の効果が図５（Ａ，Ｂ，Ｃ，Ｄ，Ｅ）に示されている。図５ＡはＸ線試験パターンの断片を示し、図５Ｂは標準的なＮＬ平均フィルタのベクトル化されたバージョン（４）によって処理されたＸ線試験パターンの断片を示し、図５Ｃは結合フィルタによって処理されたＸ線試験パターンの断片を示し、図５Ｄは図５Ａと図５Ｂとの差を示し、図５Ｅは図５Ａと図５Ｃとの差を示す。結合フィルタを適用した結果、標準的なＮＬ平均フィルタを使用した時に現れる数字のエッジに沿って、目に見えるアーティファクトのノイズ除去が行なわれる。
【００５６】
段階３。ＬＦノイズ低減。定数分散のノイズ実現と同一強度の平滑な画像の領域内のバイラテラルフィルタおよびＮＬ平均フィルタは、最も単純な線形近隣平均化フィルタの特徴を有する（すなわち、移動平均の役割を果たす）。平均化の有限領域の適用により、ＬＦに対応するスペクトルノイズ成分の一部は、低減されることはなく、平滑な画像の領域において十分に見ることができる。
【００５７】
クレームに記載の発明では、線形回帰を用いたバイラテラルフィルタを使用して、画像のピラミッドにおいて残留ＬＦノイズの低減が行なわれる[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]。画像のピラミッドは以下のように構築される。
【００５８】
【数１０】
2013175182.tif 000012
【００５９】
ここで、Ｉ（ｐ）は原画像であり、Ａｋ（ｐ），Ｈｋ（ｐ）は、画像ピラミッドのｋレベルのローパス（近似）画像および高周波数（詳細）画像であり、Ｒｅｄｕｃｅは、画像サンプリング周波数をフィルタ処理して２度低減させる（スケーリング）演算子であり、すなわちＲｅｄｕｃｅ（Ａ）＝（Ａ＊Ｌｒ）↓２（Ｌｒは線形ローパスフィルタであり、↓２はスケーリングである）であり、Ｅｘｐａｎｄは、クラウディング（crowding）サンプルおよび画像フィルタ処理演算子であり、すなわちＥｘｐａｎｄ（Ａ）＝（Ａ↑２）＊Ｌｅ（↑２はクラウディングサンプルであり、Ｌｅはローパスフィルタである）である。
【００６０】
ノイズ分散分布は、ローパスピラミッド成分（１０）にわたって推定される。そのため、標準偏差σを有するホワイトガウスノイズが結合フィルタの入力に付加される。近隣平均化フィルタとの平滑な画像の領域内のフィルタの上記の類似性のために、ピラミッドレベルの数字に対するローパスノイズ二乗平均平方根誤差依存度が得られる。
【００６１】
【数１１】
2013175182.tif 000013
【００６２】
ここで、σ０は、結合フィルタによって処理された画像ノイズの二乗平均平方根誤差レベルであり、σｋは、ピラミッド（１０）の画像ノイズＡｋ（ｐ）の二乗平均平方根誤差レベルであり、ｑ（ｋ，Ｌｒ）は、画像Ａｋ（ｐ）における二乗平均平方根誤差レベルの減少の係数である。
【００６３】
残留ノイズを低減させるために、ピラミッド（１０）の画像Ａｋ（ｐ）はバイラテラルフィルタによってフィルタ処理される。ぎざぎざのアーティファクトを低減させるために[Buades他、The staircasing effect in neighborhood filters and its solution, IEEE Transactions on Image Processing, Vol. 15 (6), pp. 1499-1505, 2006]、線形回帰を用いたバイラテラルフィルタが使用される。
【００６４】
【数１２】
2013175182.tif 000014
【００６５】
段階４。過剰平滑化補正。クレームに記載の発明では、有用なＨＦ情報の冗長な平滑化に対しては、フィルタ処理された画像とその元のバージョンとを混合することによって補正が行なわれる。フィルタ処理された画像およびその元のバージョンの画素類似性および構造的類似性について混合重みが演算される。構造的類似性を演算するために、画像画素パッチ間のユークリッド距離が演算される。画素類似性（すなわち、強度類似性）を演算するために、適切な画像の画素強度間の距離が計算される。
【００６６】
【数１３】
2013175182.tif 000015
【００６７】
最終的な画素、および画像類似性が計算された画素の構造的尺度を演算するために、構造的類似性が以下の形式で組合わせられる。
【００６８】
【数１４】
2013175182.tif 000016
【００６９】
そして、過剰平滑化補正は、以下の規則に従ってフィルタ処理された画像と元のノイズ画像との一体化を伴う。
【００７０】
【数１５】
2013175182.tif 000017
【００７１】
補正（１７）による（１６）のような類似性尺度の定義により、画像の単一の急な変化（いわゆるバースト）の除去が可能になる。
【００７２】
発明の最良の実施例
【００７３】
【数１６】
2013175182.tif 000018
【００７４】
ノイズ画像Ｎｅ（ｐ）に基づく分散演算からエッジを除去する際、正および負の変化を有する画像が形成される。
【００７５】
【数１７】
2013175182.tif 000019
【００７６】
エッジに対応する大きなエリアを選択するために、所与の二値画像に対して、２×２マスクを使用することによって、その後の拡張を伴う形態的侵食が行なわれる。
【００７７】
【数１８】
2013175182.tif 000020
【００７８】
そして、これらの画像は、原画像エッジのマップを得るために一体化される。
【００７９】
【数１９】
2013175182.tif 000021
【００８０】
区間ノイズ分散推定値を演算するために、画像Ｉｅ（ｐ）の最小強度Ｉｍｉｎおよび最大強度Ｉｍａｘが得られ、ピッチｈＭが選択され、次いで強度範囲がピッチｈＭ（たとえば、ｈＭ＝３２）を有する区間Ｍｉに分割される。各画像画素Ｉｅ（ｐ）について、この画素値を含む区間が求められ、画像Ｎｅ（ｐ）の適切な画素値を用いて所与の区間Ｍｉにおけるノイズ分散推定値σ２（ｉ）が計算され、エッジマップの値がＥ（ｐ）＝１である画素値は除外される。区間ノイズ分散推定値を計算する際、不偏推定値の式が用いられる。
【００８１】
【数２０】
2013175182.tif 000022
【００８２】
得られた区間ノイズ分散推定値σ２（ｉ）は、各区間Ｍｉについて、この区間におけるノイズ分散推定値のその後の再計算によるノイズの３標準偏差と等しいしきい値の大きさを超えるノイズ画素値を反復的に除外するように特定される。
【００８３】
【数２１】
2013175182.tif 000023
【００８４】
強度に対するノイズ分散の依存度を推定するために、得られたノイズ分散の区間推定値に基づいて補間表関数が作成される。この表関数は、ノイズ分散区間推定値のロバストな区分的線形近似を使用することによって作成される。そのため、区間Ｍｉのグリッド上では、（区間Ｍｉにおける多数の点ｎｉに依存可能にされ得る）近似化ステップｈＩおよび半径ｒＩが選択される。表関数の以前のステップ中に得られた値の近似化のために以下の規則が使用される。
【００８５】
【数２２】
2013175182.tif 000024
【００８６】
ここで、ｋは近似クノーテ（knote）の数であり、ｍｋは区間Ｍｉ（ｋ・ｈＩ）の中心であり、パラメータａ，ｂは大きさが最小の偏差に基づいて計算される[Press他、Numerical Recipes in C: The Art of Scientific Computing, Second Edition. New York: Cambridge University Press, 1999]。
【００８７】
【数２３】
2013175182.tif 000025
【００８８】
評価画像および補間表に基づいてノイズマップが生成され、画像の各画素は原画像の適切な画素におけるノイズ分散を評価する。
【００８９】
フィルタ処理品質とフィルタ速度との優れた対応関係をもたらすために、式（９）を使用して結合バイラテラル非局所フィルタ処理が行なわれ、ここで、λ＝０．８５であり、平均化領域の半径はｔ＝５であり、画素ブロックの半径はｆ＝４であるように選択され、１つの画素の重複を有する状態でブロックが選択される。強度および構造について画素類似性を計算するために、連続的な遷移を有するしきい値関数Ｆ（ｒ，σ，ｋ，ｈ）＝ｅｘｐ（−ｍａｘ（ｒ２−ｋ・σ２，０）／ｈ２）（ｋおよびｈはパラメータである）が使用される。
【００９０】
【数２４】
2013175182.tif 000026
【００９１】
構造的類似性を計算するための過剰平滑化補正では、半径５のブロックが使用される。重み関数の形態を規定するパラメータは、結合フィルタ重み関数のパラメータと一致する。強度について画素類似性の重み関数を演算する際、幾何学的距離は考慮に入れられない。
【符号の説明】
【００９２】
１  Ｘ線管、２  Ｘ線ビーム、３  検出器。
[参照文献]
【００９３】
【表１】
2013175182.tif 000027
【００９４】
【表２】
2013175182.tif 000028
【００９５】
【表３】
2013175182.tif 000029
【図１】
2013175182.tif 000030
【図２】
2013175182.tif 000031
【図３】
2013175182.tif 000032
【図４】
2013175182.tif 000033
【図５Ａ】
2013175182.tif 000034
【図５Ｂ】
2013175182.tif 000035
【図５Ｃ】
2013175182.tif 000036
【図５Ｄ】
2013175182.tif 000037
【図５Ｅ】
2013175182.tif 000038
【外国語明細書】
2013175182000001.pdf