JP2013127773A

JP2013127773A - デジタルｘ線フレームシリーズにおけるノイズ低減方法

Info

Publication number: JP2013127773A
Application number: JP2012216792A
Authority: JP
Inventors: Vasil Evich Merkur Ev Sergey; セルゲイ・バシリエビッチ・メルクリエフ
Original assignee: ZAKRYTOE AKCIONERNOE OBSHCHESTVO IMPUL'S
Current assignee: ZAKRYTOE AKCIONERNOE OBSHCHESTVO IMPUL'S
Priority date: 2011-10-07
Filing date: 2012-09-28
Publication date: 2013-06-27
Also published as: CN103177423A; KR20130038794A; EP2579207A2; EA201101502A1; US20130089247A1; EA017302B1; EP2579207A3

Abstract

【課題】デジタルＸ線透視画像のノイズ低減方法を提供する。
【解決手段】動き推定ブロックは重複によって考慮され、動き補償の間、ブロックの重複は平坦な重み付けウィンドウで平均化される。フリッカ効果補償は、基準フレームおよび現在のフレームを、所与のフレームのＬＦフィルタ処理によって得られる平均値の適切なフレームで除算することによって実施される。フリッカ効果を補償するために、基準フレームの局所的な平均強度が現在のフレームの局所的な平均強度となる。基準フレームはその平均値の動き補償されたフレームで除算され、現在のフレームの平均値のフレームで乗算される。動き補償アーティファクト補正をもたらす帰納的平均化は、所与のフレームの画素および構造的類似性の計算に基づいて現在のフィルタ処理されたフレームをその当初のバージョンと混合することを伴う。
【選択図】図１

Description

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

従来の技術水準
現在、臨床実務は、シャープネスの強調、解剖学的構造の強調、減法等の、デジタルＸ線フレームシリーズ処理のさまざまな方法を含む。デジタルＸ線フレームシリーズ品質改良の最も重要な方法は、画像ノイズ低減方法と考えられる。血管造影法プロシージャ中の医療機器の制御はリアルタイムに実行され、カテーテルを通して造影剤を注入することによって脈管系が検査される。患者／医療スタッフの放射線吸収線量を減少させるために、低用量が与えられると、低い信号対ノイズ（Ｓ／Ｎ比）比を有する医療画像の取得がもたらされる［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］。そのため、デジタルＸ線撮影では、ノイズが診断情報を守ることができる状態でのフィルタ処理方法を開発し、高い解像度値、フレーム速度および相当なレベルの信号強度依存ノイズにてリアルタイムに動作するための課題について検討することが急務である。

医療用デジタルフレームシリーズのフィルタ処理のためのアルゴリズムを開発する際に研究者が遭遇する問題がいくつかある。第１の問題は、そのレベルが信号強度に本質的に依存する画像ノイズの推定に関連する。第２の問題は、自動強度調節のシステムが使用されるときに起こる、またはＸ線管の不安定さによる、信号強度の考えられる急激な変化を補償する必要性である。第３の問題は、異なる種類の動作、つまり表またはＸ線管／検出器の動き、患者の移動または患者の体内における生理学的プロセス、血管造影剤の流れによるフレームシリーズ可変性の検討に関連する。最後に、第４の問題は、リアルタイムにフィルタ処理を実施するために、フィルタ品質と要件との均衡を観察することが必要であるという点である。この問題は、フレームの積層内の時空間平均、変換空間内の時間的フィルタ処理（ピラミッド形変換、ウェーブレット変換等）、組合わせ平均法等のフィルタ処理方法の選択も含む。これらの問題の各々について詳細に説明する。

信号依存ノイズ推定の問題について検討する。実際には、あらゆる技術水準のノイズフィルタ処理方法は、単一のデジタルフレームについては、したがってそれらのシリーズについては、ノイズスペクトル特性についての情報を使用する。比較的多数のノイズ低減方法は、ノイズ分散値をパラメータとして利用する。多くの著者の出版物は、ノイズ分散値がより正確に分かればノイズ低減アルゴリズムの品質がより高まることを実証している［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］。

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

吸収された光子のランダムな変動は光子ノイズと呼ばれる。現代の検出器では、光子ノイズが主なノイズ源である。

付加的なノイズ源は、読出しノイズ、熱的ノイズ、アンプノイズ、量子化ひずみ等の検出器ノイズを含む［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］。

σ^２（Ｉ（ｐ））＝ａＩ（ｐ）＋ｂ（１）
ここでＩ（ｐ）は、画像の画素ｐ＝（ｘ，ｙ）における信号強度レベルである。

したがって、デジタル画像ノイズは信号強度に直線的に依存する。多くの出版物が信号依存ノイズ推定の問題を扱っている［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］の著作には、デジタル医療フレームシリーズのフレームごとのノイズ推定の非パラメータ方法が記載されており、リアルタイムに動作するアルゴリズムに特に重点が置かれている。デジタルフレームノイズ推定への２パラメータ手法が［Foi他、Practical poissonian-gaussian noise modeling and fitting for single-image raw-data, Image Processing, IEEE Transactions on 17, 1737-1754, 2008］の論文に記載されている。検出器から得られ、ガンマ補正などといった非線形変換を経ていない当初のデジタル画像のノイズは、信号ランダム変数への加算性によってモデル化される：
Ｉ_ｎ（ｐ）＝Ｉ（ｐ）＋σ^２（Ｉ（ｐ））ｎ（ｐ）（２）
ここで、Ｉ_ｎ（ｐ）は観察中のノイズ画像の画素ｐ中の信号強度レベルであり、σ^２（Ｉ（ｐ））はタイプ（１）の信号強度へのノイズ分散の依存度であり、ｎ（ｐ）∈Ｎ（０，１）は正常なランダム変数である。露出不足および露出過度、つまりダイナミックレンジの境界上の線形関数（１）の妨害を引起すセンサの非線形性を考慮に入れたノイズ分散のモデル曲線を作成する方法が示唆されている。

ダイナミックフィルタを開発する際の第２の問題は、フィルタ処理中のフレームシリーズの強度の変化、いわゆるフリッカを考慮に入れ、補償することが必要な点である。この問題は、ビデオフレーム処理の一連の異なる方法の第一段階であるフリッカ効果補償タスクに類似している。［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］。自動強度調節システムまたはＸ線管の不安定さにより、デジタルＸ線フレームシリーズにフリッカが生じることがある。時間的画像シリーズ平均化においてフリッカ補償が欠如すると、フィルタ処理品質が低下する。したがって、フリッカ補償のない動きフレームにおける静止の単純な時間的平均化により、現在のフレームの推定値の相当な歪みがもたらされる。動きとブロック整合技術に基づくその補償とがある場合、フレームシーケンス強度における相当な変化により、動き推定における概算誤差と、したがってその不十分な補償とが生じる。

フリッカ推定および補償には、２つの主な局面がある。一方で、フリッカ効果は１シリーズ中の１フレームの視野に対する局所的性質を有していない。他方、それは動的過程のみにおいて、つまりフレームのシーケンスを考慮するときに生じる。これらの特徴は、フレームコントラスト／強度追従モーションにおけるいずれかの局所的な変化を識別することを伴う、フリッカ効果補償における主な問題を含むか、またはＸ線管強度変化によって引起こされる。［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］に広範に記載されているフリッカ効果モデリングの自然な方法は、以下の式に基づく。

Ｉ（ｐ）＝Ａ（ｐ）Ｉ_０（ｐ）＋Ｂ（ｐ）（３）
ここで、Ｉ_０（ｐ）は画素ｐの信号強度レベルであり、Ａ（ｐ），Ｂ（ｐ）はフレームの視野において特定される滑らかで緩やかに変動する関数であり、Ｉ（ｐ）は観察中のフレームの強度である。このモデルは、当初画像強度が観察中の画像強度に線形従属することを示唆するとともに、上記のフリッカ効果空間不均一性を伴う。相当数のフリッカ推定および補償方法はグリッドへのフレームセグメンテーションに基づき、グリッドは、重複するブロックと、考えられる線外値フィルタ処理によるＡ（ｐ），Ｂ（ｐ）関数のパラメータのブロックごとの推定とを含み、より正確な推定のために、モデルパラメータが動き推定と同時に算出される。

正確なノイズフィルタ処理方法を開発する際の第３の問題は、動きによって引起こされるフレームシーケンスの変化を考慮に入れる必要性を伴う。動きは、診断上重要ないずれかのデジタル医療用フレームシリーズの必須の属性である。非固定フレームを平均化すると、得られる画像において動きブレアーティファクトが生じる。実際には、フレームシリーズ時空間平均における動きの問題に対処するには２つの解決策がある。第１の種類の手法は、動きに関連した画像領域の検出と、これらの領域における時間的ノイズ低減の力の低下とに基づく。［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］。第２の種類の手法は、連続するフレーム間の画素の動き、したがって動き推定および補償を推定する［Brailean他、Noise Reduction Filters for Dynamic Image Sequences: A review. Proc. IEEE, vol. 83, pp. 1272-1292, 1995］。動きフレームにおける空電は、動き補償プロシージャの使用によって動き推定後に明示的に生成することができる。

第１の種類の最も単純なアルゴリズムは、リアルタイムでの実行のために実現可能な合理的な設備による。そのような方法では、動き検出は、たとえば、現在の平均化されたフレームと基準フレームとの差分をしきい値化することであり得る。ある現在のフレームの画素について、上記の差分が有限のしきい値より高い場合、時間的平均化は、たとえばその現在のフレームの空間において平均化することで置換される。しきい値は、当然ながらノイズ推定段階で得られたノイズレベルに相当することになる。そのような動き推定方法は、いくつかの欠点、インパルスノイズの視認性、ぎざぎざのあるエッジを有している。他方で、（たとえばブロック整合に基づいて）動き推定および補償に使用される第２の種類のアルゴリズムは、より高品質であると考えられる［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］、また、高い計算要求により、これらの方法はリアルタイムで実現するのが困難である。

第４の問題は、フィルタ品質とリアルタイム実行要求との間に均衡をもたらす必要性を意味する。高周波および大きなフレームサイズ（たとえば３０フレーム／秒、解像度５１２×５１２画素、１５フレーム／秒、解像度１０２４×１０２４画素／フレーム）でリアルタイムに医療用フレームシリーズ処理を実施する必要性は、現代のＸ線診断システムにおいて高品質フィルタ処理結果をもたらす要求と矛盾する。ノイズ推定または動き推定であるそのようなプロシージャそのものが、現代のコンピュータ容量に左右される大きな計算要求を有する。

発明の開示
クレームに記載の発明は、医療用途のデジタルフレームシリーズの品質を高めることを目的とする。

クレームに記載の発明の技術的成果は、ノイズ低減しやすい画質向上方法を呈示することである。

デジタルＸ線フレームの取得と、フレームごとの信号依存ノイズ分散推定と、動き推定および補償と、フリッカ効果推定および補償と、帰納的フレーム平均化とを含む、クレームに記載のデジタルフレームシリーズにおけるノイズ低減方法の技術的成果は、基準フレームにおける急激な変化と相関するノイズ画素値の形態的消去によるノイズ分散推定中に実現される。次いで、ノイズ分散の区間推定値のロバストな区分的線形近似を使用することによって、信号強度に対するノイズの依存度を表形式にし、基準デジタルフレームの画素ごとのノイズ分散推定の形態のノイズマップを決定し、動き推定および補償の際、ピラミッド方式のブロック整合が使用され、動き補償ブロック重複が平均化される際に推定ブロックが重複によって考慮され、フリッカ効果は、基準フレームおよび現在のフレームを、所与のフレームの線形ＬＦフィルタ処理によって得られる平均値の対応フレームで除算することによって得られ、フリッカ効果を補償する際、局所的な平均基準フレーム強度は局所的な平均現在フレーム強度となり、したがって基準フレームはその平均値の動き補償されたフレームで除算され、現在のフレームの平均値のフレームで乗算され、帰納的平均化の際、動き補償アーティファクト補正は、それらの画素および構造的類似性の算出に基づいて、現在のフィルタ処理されたフレームをその当初のバージョンと混合することによって実施される。

クレームに記載の方法の考えられる実施例バージョンでは、フレームのエッジに対応する画素の算出された推定値に対する影響を減少させる（エッジの除去）ためのノイズ分散推定において、現在のフレームと動き補償された基準フレームとの差分が使用される。

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

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

フリッカ効果補償を考慮に入れるクレームに記載の方法に最も近いものは、［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］。動きベクトルの推定品質に対して大幅に影響するフリッカ効果の考察は、動き推定段階において現在のフレームおよび基準フレームがそれらの平均値フレームで除算される最も単純な手法、つまり一種の輝度正規化によって実施される。フリッカフレーム強度の考察に加え、動きの推定は、ピラミッド方式を用いることによってなされる。その理由は、医療用フレームシリーズのリアルタイムのフィルタ処理中は、速い動きの推定および補償（たとえば、検査対象物に密接に接近する際の心臓検査−冠動脈造影法中）が特に困難であるという事実にある。その結果、特に低いフレーム速度での動き補償に成功するためには、適切なブロック整合領域が大きくなりすぎる。したがって、速い動き補償中の計算の複雑度が高まるという問題を解決するために、本発明では動き推定にピラミッド形アルゴリズムが使用され、
・現在のフレームおよび動き補償されるべき基準フレームは、いくつかのレベルのピラミッド状のフレーム（検査種類に応じてせいぜい３レベル）に分解され、
・ブロック動きベクトル推定は、ピラミッドレベルの数を考慮してその半径が選択される小領域上を走査することによってＬＦピラミッド構成要素において行なわれ、
・ピラミッドレベルを最低から最高に変換する際、動き補償ブロックのサイズが変倍され、小半径の領域におけるサーチによってその動きベクトルが特定される。

ピラミッド方式の利用は、基準フレームにおける好適なブロック整合領域を大幅に拡張する。矩形ブロックに基づく動き補償に共通の過度なブロッキングを減少させるために、クレームに記載の発明は、重複ブロック補償を提案する［Orchard他、Overlapped Block Compensation: An Estimation-Theoretic Approach. IEEE Transactions On Image Processing, Vol. 3, No. 5, September 1994, pp. 693-699］。典型的な重複サイズは、ブロックの半径の４分の１に等しい。動き補償されたフレームを作成する際、重複したブロック領域に属する画素が平滑重み付け関数によって平均化される。

クレームに記載の発明の示差的特徴は、動き補償アルゴリズムのアーティファクトの改善の特別な方法の利用である。サーチ（整合）領域の半径が小さすぎるために概算動き補償誤差が発生する可能性があり、重複する対象物またはフレームコンテンツの急激な変化は、フィルタ処理されたフレームと、その空間で平均化されたバージョンとの混合によって減少する。混合重みは、画素および、現在のフレームの、動き補償されたフレームとの構造的類似性に基づいて計算される。画素類似性の尺度として、相互の距離が用いられる［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］。画素および構造的類似性は、線形的に混合される。

効果的なノイズ低減をリアルタイムに提供することを目的としたクレームに記載の発明では、信号依存ノイズの自動推定、フレームシリーズの考えられる急激な強度変化の推定、動き推定および補償、動き補償アーティファクトの改善、ならびに時間的過剰平滑化によるカルマンフィルタが使用される。与えられたタスクを、共有ＰＣのマルチコアプロセッサとｎｖｉｄｉａ（登録商標）のｃｕｄａ（登録商標）対応のグラフィックスカードとに振り分けると、リアルタイムで高品質なフィルタ処理要求に適した結果を得ることに成功した。

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

Ｘ線フレーム取得のための装置を示す図である。図４に表される血管の研究のフレームシリーズについて、信号強度に対する正常なノイズ偏差の依存度を表す表関数を示す図である。図４に表されるフレームのための正常なノイズ偏差マップを示す図である。フレームシリーズの当初のフレームを示す図である。ノイズ除去レベルが８０％（ｒ＝０．８、式（１２））である、所与のフレームシリーズフィルタ処理の結果を示す図である。図４および図５の理解を向上させるためにシャープネス強調フィルタを使用した。図４および図５に表されるフレームの差分画像を示す図である。

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

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

段階１。信号強度に依存するノイズ分散のフレームごとの推定。この段階では、現在のフレームのＬＦ線形フィルタ処理によってフレーム推定が実施される［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］。したがって、信号依存ノイズフィルタ処理の問題は、付加的な信号非依存ノイズ分散を減少させるという課題にある。クレームに記載の発明では、非パラメータ的ノイズ推定を実施しノイズマップを構築しつつ、以下の手法が提案される。推定された画像および補間表に基づいて、画像であるノイズマップが構築され、その各画素は、適切な当初の画像画素における平均平方ノイズ偏差を推定する。ノイズマップは、実用化に十分な精度の画素ごとのノイズ推定値をもたらす。ノイズ分散推定中に、フレームエッジに対応するフレーム画素の推定値の算出に対する影響を低下させる（エッジの除去）ために、現在の動き補償されたフレームと基準の動き補償されたフレームとの差分が使用される、クレームに記載の発明の実施例バージョンが可能である。

段階２。動きおよびフリッカ効果の推定および補償。クレームに記載の発明では、フレームシリーズ強度の考えられる可変性を考慮したブロック対応のサーチに基づいて動き推定が実施される。この目的のため、動きベクトルサーチプロシージャを使用する前に、現在のフレームおよび基準フレームがそれらの平均値のフレームで除算される。対応するシリーズフレームの平均値のフレームは、動き推定ブロックのサイズ、たとえば１６×１６画素に匹敵するフィルタ孔を備えた低域通過フィルタ処理によって得られる。重複がブロックサイズの１／４に等しいとき、動き推定ブロックが重複して与えられる。ブロック整合の際、動き推定ピラミッド形アルゴリズムが使用され、その間、
・現在のフレームおよび現在のフレームに対して動き補償されるべき基準フレームがいくつかのピラミッドレベルに分割され、
・ブロック動きベクトル推定がＬＦフレームピラミッドで実施され、
・下方から上方へのピラミッドレベルの移行中、ブロックサイズおよびそれらの重複が変倍され、小半径領域に対するサーチによって、所与のブロックの動きベクトルが特定され、
・動き補償されたフレームを構築する際、重複領域が平均化される。

基準フレームから現在のフレームへの推定は、ＬＦピラミッドのフレームにおいて実施される。せいぜい３つの分解レベルが使用される。原則として２つである。なぜなら、より多くの分解レベルを微細な物体（薄い容器）の動き推定に使用すると、動きベクトル推定値において概算誤差が存在し得るからである。下方から上方へのピラミッドレベルの移行中、ブロックサイズおよびそれらの重複が倍増され、小領域、たとえば３×３画素に対するサーチによって動きベクトルが特定される。ブロック類似性の尺度であるユークリッド平方メートル法を計算しつつ、早く算出を停止させる技術が使用される［Wang他、Fast Algorithms for the Estimation of Motion Vectors. IEEE Trans. Image Processing, vol. 8, no. 3, pp. 435-438, March 1999］。

動き推定後、補償された基準フレームが形成され、過度のブロッキングを減少させるために、ブロックの重複が平均化される。ブロックの重複を平均化するために、平滑重み付けウィンドウが使用される。さらに、フリッカ効果を補償するために、基準フレームの局所的手段が現在のフレームの局所的手段となる。

段階３。動き補償アーティファクト補正をもたらす帰納的な平均化。整合領域の半径が小さすぎることにより生じる動き補償概算誤差、重複する物体、またはフレームコンテンツの急激な変化は、クレームに記載の発明では、フィルタ処理されたフレームとその空間平均化されたバージョンとの混合によって減少される。混合重みは、現在のフレームおよび動き補償されたフレームの画素および構造的類似性に基づいて算出される。画素類似性の尺度として、相互距離が使用される［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］。

Ｉ_ｒ（ｐ），Ｉ_ｃ（ｐ）を、画素ｐ中の基準フレームおよび現在のフレームの強度値とする。

画素と構造的類似性とを組み合わせるフレームの、結果として得られる類似性の尺度の計算については、以下の式：
Ｗ（ｐ）＝λ・Ｗ_ｓ（ｐ）＋（１-λ）・Ｗ_ｐ（ｐ）（５）
が使用される。

λは、適切な画素強度に基づいた重みに対して優勢な構造的重みを特定する。フレーム混合中、λ≒１（たとえば０．９）での混合された重み計算は、画素強度の単一インパルス値の影響を減少させることと、同時に（当初の現在の補償されていないフレームのコンテンツと同様の）正確な構造的フレームコンテンツを再構築することとをもたらす。

時間的平均化は、以下の式にしたがって、動き補償された基準フレーム（フレーム予測子）および現在のフレーム（観察されたフレーム）について一次元カルマンフィルタ処理の使用によって反復的に実施される。

発明の最良の実施例

区間ノイズ分散推定値を計算するために、フレーム強度Ｉ_ｅ（ｘ，ｙ）のＩ_ｍｉｎおよびＩ_ｍａｘが求められ、強度範囲をピッチｈ_Ｍ（ｈ_Ｍ＝３２）で区間Ｍ_ｉに分割する。フレームＩ_ｅ（ｘ，ｙ）の各画素について、この画素を含む区間が求められ、Ｎ_ｅ（ｘ，ｙ）の適切な画素値を用いて所与の区間Ｍ_ｉにおけるノイズ分散推定値σ^２（ｉ）が算出され、エッジマップ値Ｅ（ｘ，ｙ）＝１を有する画素値は除外される。

区間ノイズ分散推定値を計算する際、不偏分散推定値の式が適用される。

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

取得されたノイズ分散の区間推定値を使用した、強度に対するノイズ分散の依存度の推定について、補間表関数が構築される（図２の円）。この表関数は、区間分散推定値のロバストな線形近似に基づく。

各画素が当初のフレーム（図３）の適切な画素におけるノイズ分散を推定するフレームであるノイズマップは、推定されたフレームおよび補間表に基づいて構築される。

動き推定および補償は、ブロック整合およびフレームシリーズ強度の考えられる可変性を伴う。そのため、動き推定の前に、現在のフレームおよび基準フレームがそれらの平均値のフレームで除算される。平均値のフレームは、動き推定ブロックの寸法、たとえば１６×１６画素に匹敵するフィルタ孔を備えた低域通過フィルタ処理により取得される。動き推定ブロックは、重複がたとえば４に等しい場合には、重複によって与えられる。ブロック整合の際、動き推定ピラミッド形アルゴリズムが使用され、その間、
・現在のフレームおよび現在のフレームに対して動き補償されるべき基準フレームがいくつかのピラミッドレベルに分割され、
・ブロック動きベクトル推定がＬＦ画像ピラミッドで実施され、
・下方から上方へのピラミッドレベルの移行中、ブロックサイズおよびそれらの重複が変倍され、小半径領域に対するサーチによって所与のブロックの動きベクトルが特定され、
・動き補償されたフレームを構築する際、重複領域が平均化される。

基準フレームから現在のフレームへの推定は、ＬＦピラミッドのフレームにおいて実施される。せいぜい３つの分解レベルが使用される。原則として２つである。なぜなら、微細な物体（薄い容器）の動き推定のためにより多くの分解レベルを使用すると、動きベクトル推定において概算誤差が存在し得るからである。下方から上方へのピラミッドレベルの移行中、ブロックサイズおよびそれらの重複が倍増され、小領域たとえば３×３画素に対するサーチによって動きベクトルが特定される。ブロック類似性の尺度であるユークリッド平方メートル法を計算しつつ、早く算出を停止させる技術が使用される［Wang他、Fast Algorithms for the Estimation of Motion Vectors. IEEE Trans. Image Processing, vol. 8, no. 3, pp. 435-438, March 1999］。

動き推定後、補償された基準フレームが形成され、過度のブロッキングを減少させるために、ブロックの重複が平均化される。基準フレームノイズマップが同様の方法で補償される。ブロックの重複を平均化するために、平滑重み付けウィンドウ、たとえばコサインが使用される。

さらに、フリッカ効果を補償するために、基準フレームの局所的手段が現在のフレームの局所的手段となる。そのため、補償された基準フレームは、その動き補償された平均フレームで除算され、式４に従って現在のフレームの平均値のフレームで乗算される。基準フレームのノイズマップに対して、同様の平均化を行う。

現在のフレームの帰納的平均化。時間的平均は、式（６）〜（１２）に従って動きおよび強度が補償されたフレームについて一次元（画素ごとの）カルマンフィルタ処理の使用によって帰納的に実施される。

動き補償誤差によって引起こされる時間的平均化アーティファクトを修正するための、フィルタ処理された正規化されたフレームと、その当初のバージョンとの混合（式（５）、（９）および（１０））。現在のフレームと動きおよび強度が補償されたフレームとに基づいて混合重み（５）が算出され、構造的類似性の尺度として、現在のフレームおよび補償されたフレームの画素ブロック間のユークリッド距離が使用される。

構造的類似性の算出には、その寸法が選択されるたとえば１１×１１画素のスライディングブロックが与えられる。画素類似性の算出は、現在のフレームおよび基準の動き補償されたフレーム間の強度空間の距離として計算される。実際のシリーズにおける式（１１）、（１２）中のノイズ低減度ｒは、たとえばｒ＝０．８が選択される。これは５分の１のノイズ低減に相当する。図４〜図６は、実際のシリーズ（血管造影法検査）へのフィルタの適用の結果を示す。図４は、シリーズの当初のフレームを示し、図５はフィルタ処理されたフレームを示し、図６は、図４と図５とのフレームの差分を示す。よりよく理解するため、当初のフレームおよびフィルタ処理されたフレームに対してシャープネス改善プロシージャを行う。

１Ｘ線管、２Ｘ線ビーム、３検出器。

Claims

デジタルＸ線フレームシリーズのノイズ低減方法であって、デジタルＸ線フレームシリーズの取得と、信号依存ノイズ分散のフレームごとの推定と、動き推定および補償と、フリッカ効果推定および補償と、帰納的フレーム平均化とを含み、
信号強度依存ノイズ分散のフレームごとの推定は、当初のフレームのエッジに対応するノイズフレームの画素値の形態的消去を伴い、
信号強度に対するノイズの表形式依存度は、ノイズ分散の区間推定値のロバストな区分的線形近似の使用によって得られ、
ノイズマップは、当初のデジタルフレームのノイズ分散の画素ごとの推定として計算され、
動きおよびフリッカ効果の推定および補償の間、ブロック整合のためのピラミッド方式が利用され、
動き推定ブロックは重複によって考慮され、動き補償中に、ブロック重複が平滑重み付けウィンドウを適用することによって平均化され、
フリッカ効果補償は、基準フレームおよび現在のフレームを、所与のフレームのＬＦフィルタ処理によって得られる平均値の適切なフレームで除算することによって実施され、
フリッカ効果を補償するために、基準フレームの局所的な平均強度は現在のフレームの局所的な平均強度となることを特徴とし、
基準フレームは、その平均値の動き補償されたフレームで除算され、現在のフレームの平均値のフレームで乗算され、
動き補償アーティファクト補正による帰納的平均化は、所与のフレームの画素および構造的類似性の計算に基づいて、現在のフィルタ処理されたフレームをその当初のバージョンと混合することを伴う、方法。
フレームのエッジに対応する計算された推定画像画素への影響を減少させることを目的とするノイズ分散推定は、現在のフレームと基準の動き補償されたフレームとの差分を要する、請求項１に記載の方法。