(57)【要約】【課題】  デジタルＸ線画像ノイズ推定の速度および品質の観点から、より効率的な方法を策定する。【解決手段】  デジタルＸ線フィルムノイズアセスメントの方法は、原画像の取得、推定画像を取得するための原画像の低周波フィルタリング、原画像および推定画像の間の差分としてのノイズ画像作成、形態的フィルタリングによるノイズ画像画素の除去、推定画像の強度範囲の区間への分割を含む。さらに、方法は、各区間について、推定画像の画素に対応するノイズ画像の画素を蓄積すること、ノイズ分散の区間推定値を算出すること、σ３基準に従った除去ノイズ画素によって区間推定値を改善すること、信号強度に対するノイズの依存度を表す表形式関数をもたらすノイズ分散の区間推定値のロバストな局所線形近似、推定画像および表形式関数に基づいて、デジタル原画像の画素ごとのノイズ分散推定値としてノイズマップを算出することを含む。【選択図】図１
2012150802.tif 000002
【特許請求の範囲】
【請求項１】
デジタルＸ線フィルムノイズアセスメントの方法であって、
原画像を取得するステップと、
原画像の低周波フィルタリングによって、推定された画像を取得するステップと、
原画像および推定された画像の間の差分としてノイズ画像を作成するステップと、
原画像における急激な変化に対応するノイズ画像画素を除去するステップと、
推定された画像の強度範囲を区間に分割するステップとを含み、これにより、推定された画像の各画素は適切な区間と関連し、さらに、
各区間について、推定された画像の画素に対応するいくつかのノイズ画像の画素を蓄積するステップと、
そのような区間で蓄積されたノイズ画像の画素を使用して、ノイズ分散の区間推定値を算出するステップと、
σ３基準に従った除去ノイズ画素の使用によって区間推定値を改善するステップとを含み、
原画像中のエッジに対応するいくつかのノイズ画像画素の形態的抽出の使用によって、原画像における急激な変化に対応するいくつかのノイズ画像画素の除去が行なわれ、
ノイズ分散の区間推定値のロバストな局所線形近似が行なわれ、それにより、信号強度に対するノイズの依存度を表す表形式関数が得られ、
推定された画像および、取得された、信号強度に対するノイズの依存度を表す表形式関数に基づいて、デジタル原画像の画素ごとのノイズ分散推定値としてノイズマップが算出される、方法。
【発明の詳細な説明】
【技術分野】
【０００１】
発明の分野
本発明はデジタル画像処理に関し、Ｘ線を含む高エネルギ放射の使用によって取得されるデジタル画像の処理に関する課題を解決するために使用することができる。本発明は、特にデジタルＸ線フィルムのノイズアセスメントに関する。
【背景技術】
【０００２】
背景
現在、デジタルＸ線撮影においてさまざまな画像処理アルゴリズムが使用されている。クリスプニングまたは解剖的組織分割などの方法は、画像ノイズレベルデータを使用することができる。また、実際には、あらゆる定性ノイズ抑制方法はノイズ分散をパラメータとして使用する。したがって、画像のノイズレベルを判定するためにのみデジタル原画像を有するという問題を考慮することが重要である。この問題は、デジタル画像のノイズ分散が有用な信号強度に大きく依存することによって複雑化する。
【０００３】
【数１】
2012150802.tif 000003
【０００４】
吸収された光子のランダムな変動は、量子ノイズまたは光子ノイズと呼ばれる。現代の検出器では、光子ノイズが画像ノイズの主な源である。追加的なノイズ源は、検出器システムのノイズ：読出しノイズ、熱ノイズ、増幅器のノイズ、量子化ノイズ等に関連する［５］。所与のノイズ源の全体的な影響は、ガウス分布ランダム量によってモデル化することができる［１〜３］。線形の電子回路における従来のモデルによれば、デジタル画像の全体的なノイズ（光子ノイズおよび追加的なノイズ源）の分散は、有用な信号の線形関数である［１１］。
【０００５】
【数２】
2012150802.tif 000004
【０００６】
ここで、Ｉ（ｐ）は画素ｐの信号強度レベルである。
【発明の概要】
【発明が解決しようとする課題】
【０００７】
ノイズに関するデジタル画像アセスメントの問題を扱っている多くの出版物がある［１，３，４，６〜８，１０］。論文［６］は、非パラメトリックノイズアセスメント方法を提示し、主な強調点は、リアルタイムアルゴリズムを策定することに集中している。ノイズアセスメントへの２パラメータ手法が論文［３］で考慮されている。（画像検出器から直接得られ、かつガンマ補正等といった非線形変換を受けていない）生データデジタル画像のノイズは、信号ランダム値に関連した付加物としてモデル化され、その分散は法則（１）に従った信号に依存する。モデルノイズ分散曲線をグラフ化するこの方法は、露出不足および露出過度、すなわちダイナミックレンジの境界での法則（１）の妨害を招く検出器動作における非線形性を考慮に入れる。
【課題を解決するための手段】
【０００８】
発明の概要
本請求項に最も密接に関連する方法は、論文［２，６］に提示されている方法であり、それによれば原画像に対する信号依存ノイズの推定は、以下の段階を実行するのに十分である。
・生データデジタル画像の低周波フィルタリングによって有用な信号を推定し、原画像と、それにより得られた推定ノイズ画像との間の差分を算出すること。
・いずれかの方法の使用によって、生データデジタル画像における急激な変化（境界、単一の「ホット」画素）と相関するノイズ画像を有する画素値を除去すること。
・推定された画像の強度範囲を区間に分割し、次いで各区間について、推定された画像の画素に対応するノイズ画像の画素を記憶すること。
・そのような区間で記憶されたノイズ画像の画素値に基づいて、すべての区間におけるノイズ分散を算出すること。
【０００９】
本発明は、上記の論文に述べられているノイズ推定原理を使用する。請求項に記載の発明の課題は、デジタルＸ線画像ノイズ推定の速度および品質の観点から、より効率的な方法を策定することである。
【００１０】
原画像を取得すること、原画像の低周波フィルタリングによって、推定された画像を取得すること、生データデジタル画像およびその推定値の間の差分として、信号依存ノイズを作成すること、生データデジタル画像における急激な変化に対応するノイズ画像を有する画素値を除去すること、および推定された画像の強度範囲を区間に分割することを含み、推定された画像の各画素は適切な区間に属し、さらに、各区間について、推定された画像の画素に対応するノイズ画像の画素を記憶すること、そのような区間で記憶されたノイズ画像の画素値に基づいて、すべての区間においてノイズ分散を算出すること、およびσ３手法を用いてノイズ画素を除去することによって、すべての区間における分散を形成することを含むデジタルＸ線フィルムのための請求項に記載のノイズアセスメント方法における技術的成果は、以下の手段によって実現される。
・生データデジタル画像における急激な変化と相関する、ノイズ画像を有する画素値は、原画像のエッジに対応する、ノイズ画像を有する画素値の形態的抽出によって除去される。
・ノイズ分散の区間推定値のロバストな局所線形近似が行なわれ、信号強度に対するノイズの依存度を示す表形式関数が得られる。
・推定された画像と計算された表形式関数とに基づいて、原画像ノイズ分散の画素ごとの推定値である画像のように見えるノイズマップを算出する。
【００１１】
アルゴリズムの実行は、いくつかの段階を含む。
１．  生データデジタル画像の低周波フィルタリングの使用による有用な信号の推定、および原画像とその推定値との間の差分を算出することによるノイズ画像の取得。
２．  生データデジタル画像のエッジに対応する、ノイズ画像を有する画素値の形態的抽出。
３．  推定された画像の強度範囲を区間に分割すること、およびそのような各区間について、この区間において記憶されたノイズ画像の画素に基づいてノイズ分散を算出すること。
４．  σ３手法に基づいてノイズ画素を除去することによって、各区間において分散を形成すること。
５．  信号強度に対するノイズの依存度を示す表形式関数を計算するためのノイズ分散の区間推定値のロバストな局所線形近似。
６．  推定された画像と、ノイズおよび強度の間の決定された表形式の依存度とに基づいてノイズマップを算出すること。
【００１２】
添付の図面は発明を例示し、書面での説明と合わせて、発明の原理について説明するように機能する。実施例の同じまたは同様の要素を参照するために、可能な限り、図面全体にわたって同じ参照番号が使用される。
【図面の簡単な説明】
【００１３】
【図１】Ｘ線フィルムを得るための装置を示す図である。
【図２】実際のＸ線フィルムの断片を示す図である。
【図３】境界（非ノイズ画素）のマップを示す図である。
【図４】表形式関数プロットを示す図である。
【図５】ノイズマップを示す図である。
【発明を実施するための形態】
【００１４】
発明の詳細な説明
Ｘ線フィルムを得るための可能な変形例の１つは、図１に示される装置の使用によって実行される。当該装置は、Ｘ線ビーム３を発生するＸ線管１を含む。Ｘ線ビーム３はＸ線検査下にある物体２を透過する。Ｘ線は検出器４によって知覚される。検出器４は、Ｘ線発光スクリーンおよび感光性素子のマトリックスを含む。上記スクリーンは、感光性面の能動部品と光学的に接続される。入射光３は、まず発光層において可視光に変換され、その強度が感光性素子のマトリックスによって記録される。感光性素子のマトリックスからデジタルデータが読み出され、画像としてモニタに表示される。
【００１５】
推定された画像作成およびノイズ画像取得の段階において、原画像の有用な信号推定が、たとえば原画像の低周波フィルタリングの使用（線形低周波数畳み込み、メジアンフィルタリング等）によって実行される［２，６］。リアルタイム用途における実行速度に厳しい要件を与えるために、最も単純な線形フィルタリング（たとえば二項フィルタ）を使用することが合理的である。得られた平滑化画像を使用して、ノイズ画像、すなわち原画像とフィルタリングされた画像との差分が発生する。
【００１６】
最も単純なフィルタによる推定は理想的ではなく、エッジが過剰に平滑化される。その結果、原画像と平滑化画像との差分は、平滑領域のノイズ画素に加えて、急激な変化（解剖組織輪郭、以下非ノイズ画素）に対応するある数の画素を含むノイズ画像をもたらす。これらの画素は、ノイズ分散値を大幅に歪ませる可能性があるため、ノイズ統計値を算出する際には除去されるべきである。この目的のために、概して原画像からしきい値設定平滑化導関数に至る異なる方法［３，８］が使用され、信号対ノイズ局所推定値によってしきい値が決定される。このような推定値は、多数の細部を含む画像領域では不十分であることを証明する。このため、本発明は、エッジの検出および除去の段階において、標準偏差の導関数および局所推定値を算出することを必要としない、エッジを選択するためのより容易な手法を提示する。非ノイズ画素を除去するためのこの形態的手法の本質は以下のとおりである。
１．  ノイズ画像は２つの成分、正の変化および負の変化の二値画像に分割される。
２．  エッジに対応する領域を選択するために、取得された画像に対して侵食（ｅｒｏｓｉｏｎ）および拡張（ｄｉｌａｔｉｏｎ）などの形態的操作を行なう。
３．  処理された画像を組み合わせて１つの二値画像、原画像の外形マップを取得する。
より完全な記憶薄肉構造を目的として、侵食および拡張の形態的操作は、小型のマスク（ウィンドウ２×２）の使用によって行なわれる。図３は、図２に表示されている実際の放射線画像のエッジ評価へのここに述べる手法によって得られた二値画像を示す。
【００１７】
ノイズ分散の区間推定の段階において、推定された画像の最小強度および最大強度（強度範囲限界）を算出し、たとえば３２階調に等しい分割ピッチを選択する。さらに、推定された画像の各画素について、この画素値が属する区間と、所与の区間においてノイズ分散評価を算出するのに使用するノイズ画像の対応画素値とを求める（同時に、境界に対応する画素を除去する）。ノイズ分散の区間推定値を算出する際、異なる公式、たとえば標準的な不偏推定値、または絶対偏差のメジアンの式を使用するロバストな推定値を使用することができる［２，３，６］。このプロシージャは、信号強度に対するノイズ分散の依存度を表わす表形式関数をもたらす。
【００１８】
エッジのマップ生成において出現する不正確さは、統計の推定値から非ノイズ画素を正確に除去するための障害であり、分散の区間推定値を算出する際に総誤差に繋がる可能性がある。したがって、分散の区間推定値を改善する段階において、推定データは、反復異常値除去技術［６］によって改善される。各区間についてノイズ画素が反復して除去され、その絶対値は３標準ノイズ偏差と等しいしきい値を上回り、その後この区間においてノイズ分散推定値が再算出される。
【００１９】
ノイズ分散の区間推定値が算出された後、信号強度に対するノイズ分散の依存度の推定の段階となる。パラメータ評価の際、いずれかの方法（たとえば最小二乗法、尤度関数の最小化、ディレクティブ最適化）によってノイズモデルパラメータ（１）を推定することができる。センサの非線形性を考慮に入れることもできる［２］。しかし文献［６］で言及されているように、ノイズ分散と信号強度との関係を十分に説明するパラメトリックモデルの生成は、多くの要因から重大な困難を招き得る。これらの要因は、センサの非線形性、原画像データの非線形的な前処理（たとえば対数を取ること）を含み得る。したがって、実施の容易さおよびアプリケーションの汎用性の点から、区間推定値に基づいて信号強度に対するノイズ分散の依存度の非パラメトリック推定を行なう手法がより有益と思われる。
【００２０】
本発明では、所要の依存度を生じさせるために、取得されたノイズ分散の区間推定値に基づいて補間表形式関数が生成される非パラメトリック手法が提供される。この表形式関数は、ノイズ分散の区間推定値のロバストな局所線形近似に基づいて作成される。ロバストな技術の利用は、異常値作用（ノイズ分散の区間推定値における総誤差）の低減を追加的に可能とし、一方で近似の局所性は、強度に対するノイズの依存度を示す曲線の複雑な傾向の反復をもたらす。したがって、原画像の各強度について得られた表形式関数は、ノイズ分散の推定値と一致する。表形式入力点は、たとえば、推定された画像の強度を使用し得る。図４は、図２に示される画像の標準ノイズ偏差（円）の区間推定値に基づく、強度に対する標準ノイズ偏差の依存度を表わす表形式関数（実線）のプロットを示す。Ｘ軸は図２に示された画像の信号強度値を表し、Ｙ軸は所与の強度値に対応する標準ノイズ偏差値を表す。
【００２１】
実際には、パラメトリックノイズ推定の場合、原画像のノイズ分散を安定化させる変換が使用される手法が使用され得る［２，３，９，１１］。これにより、信号に依存するノイズフィルタリングの問題は、付加的な、ノイズから独立した分散定数（所与）を抑制するためのタスクに至る。本請求項は非パラメトリックノイズ推定を提供し、したがってノイズマップ作成の段階において以下の手法を使用することが提示される。すなわち、推定された画像および補間表を使用してノイズマップが作成され、その画像の各画素は、原画像の適切な画素におけるノイズ偏差の平均二乗値を推定する。ノイズマップは、実施に十分な精度で画素ごとのノイズ推定をもたらす。画像ノイズレベルの可能な限り正確なデータを使用することで、大幅に向上したノイズ減少アルゴリズムが得られる。図５はノイズマップを示し、その画像の各画素はノイズ標準の推定値をもたらす。
【００２２】
発明の最良の実施例
推定された画像の作成およびノイズ画像取得の段階において、サイズ３×３の低周波線形二項フィルタによって原画像Ｉ（ｘ，ｙ）がフィルタリングされる。
【００２３】
【数３】
2012150802.tif 000005
【００２４】
これにより、平滑化画像Ｉｅ（ｘ，ｙ）＝Ｉ＊Ｈが得られる。さらに、ノイズ画像Ｎｅ（ｘ，ｙ）＝Ｉ（ｘ，ｙ）−Ｉｅ（ｘ，ｙ）が算出される。
【００２５】
ノイズ画像Ｎｅ（ｘ，ｙ）を使用したエッジの除去の段階において、正の変化および負の変化の２つの画像が形成される。
【００２６】
【数４】
2012150802.tif 000006
【００２７】
エッジに対応する大領域を選択するために、サイズ２×２のマスクを使用した侵食および拡張の形態的操作を連続的に二値画像データに対して行なう。
【００２８】
【数５】
2012150802.tif 000007
【００２９】
次いで、これらの画像を組み合わせて、原画像のエッジのマップを取得する。
【００３０】
【数６】
2012150802.tif 000008
【００３１】
ノイズ分散の区間推定の段階は、画像Ｉｅ（ｘ，ｙ）の最小強度Ｉｍｉｎおよび最大強度Ｉｍａｘの算出と、ステップｈＭの選択と、ステップｈＭ（ｈＭ＝３２）による区間Ｍｉへの強度範囲の分割とを含む。画像Ｉｅ（ｘ，ｙ）の各画素について、これらの画素値を含む区間が判定され、画像Ｎｅ（ｘ，ｙ）の適切な画素値を使用して所与の区間Ｍｉにおけるノイズ分散推定値σ２（ｉ）を算出する。
【００３２】
【数７】
2012150802.tif 000009
【００３３】
各区間Ｍｉについてσ２（ｉ）の区間分散推定値を改善する段階においてノイズ画素が反復的に除去され、その絶対値は３標準ノイズ偏差と等しいしきい値を上回り、その後この区間においてノイズ分散推定値が再算出される。
【００３４】
【数８】
2012150802.tif 000010
【００３５】
取得された区間ノイズ分散推定値に基づく、信号強度に対するノイズ分散の依存度の推定の段階により、補間表形式関数が得られた。この表形式関数は、区間分散推定値のロバストな局所線形近似に基づく。このため、区間グリッドＭｉでは、（区間Ｍｉにおける点の数ｎｉに依存し得る）ステップｈＩおよび近似半径ｒＩが選択される。先の段階で取得された表形式関数の値は、以下の式にしたがって近似される。
【００３６】
【数９】
2012150802.tif 000011
【００３７】
推定された画像と補間表とに基づくノイズマップ作成のこの段階によってノイズマップが得られ、その画像の各画素は、原画像の適切な画素におけるノイズ分散を推定する。
【００３８】
参考文献
１．Bosco A., Battiato S., Bruna A.R., Rizzo R. Noise reduction for cfa image sensors exploiting hvs behavior. Sensors 9(3), 1692-1713 (2009)
２．Foi A., Trimeche M., Katkovnik V., Egiazarian K. Practical poissonian-gaussian noise modeling and fitting for single-image raw-data. Image Processing, IEEE Transactions on 17, 1737-1754 (October 2008)
３．Foi А. Practical denoising of clipped or overexposed noisy images. Proc. 16th European Signal Process. Conf., EUSIPCO 2008, Lausanne, Switzerland, August 2008
４．Forstner W. Image preprocessing for feature extraction in digital intensity, color and range images. In Springer Lecture Notes on Earth Sciences, 1998
５．Gino M. Noise, Noise, Noise
http://www.astrophys-assist.com/educate/noise/noise.htm
６．Hensel M., Lundt B., Pralow T., Grigat R.-R. Robust and Fast Estimation of Signal-Dependent Noise in Medical X-Ray Image Sequences. In Handels, H., et al., ed.: Bildverarbeitung fur die Medizin 2006: Algorithmen, Systeme, Anwendun-gen. Springer (2006) 46-50
７．Liu C., Szeliski R., Kang S.B., Zitnick C.L., Freeman W.T. Automatic estimation and removal of noise from a single image. IEEE Transactions on Pattern Analysis and Machine Intelligence 30(2), 299{314 (2008)
８．Salmeri M., Mencattini A., Rabottino G., Lojacono R. Signal-dependent Noise Characterization for Mammographic Images Denoising. IMEKO TC4 Symposium (IMEKOTC4 '08), Firenze, Italy, September 2008
９．Starck J., Murtagh F., Bijaoui A. Image Processing and Data Analysis: The Multiscale Approach. Cambridge University Press, 1998
１０．Waegli B. Investigations into the Noise Characteristics of Digitized Aerial Images. In: Int. Arch. For Photogr. And Remote Sensing, Vol. 32-2, pages 341-348, 1998
１１．Яне Б. Цифровая обработка изображений. М., Техносфера, 2007, 583 c
１２．W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery. Numerical Recipes in C: The Art of Scientific Computing, Second Edition. New York: Cambridge University Press, 1992
【符号の説明】
【００３９】
１  Ｘ線管、２  物体、３  Ｘ線ビーム、４  検出器。
【図１】
2012150802.tif 000012
【図２】
2012150802.tif 000013
【図３】
2012150802.tif 000014
【図４】
2012150802.tif 000015
【図５】
2012150802.tif 000016