画像処理装置、画像処理方法、プログラムおよびコンピュータ記録媒体
【課題】 ノイズの影響を考慮することで動画像から動きを精度良く判定し、その判定結果に応じた好適なノイズ低減を行う。
【解決手段】 複数フレームの動画像を時系列に重み付け加算する画像処理装置であって、
前記複数フレームの動画像の中で、現フレームより前のフレーム画像をリカ―シブフィルタ処理するフィルタ手段と、
前記現フレームの画像と前記フィルタ手段で得られた画像との差分画像におけるノイズの統計量に基づいて現フレームの重み付け係数を得る係数取得手段と、
前記係数取得手段で得られた前記複数フレームに対応する重み付け係数に基づき前記複数フレームの動画像を加算する加算手段と、
を有する。
【解決手段】 複数フレームの動画像を時系列に重み付け加算する画像処理装置であって、
前記複数フレームの動画像の中で、現フレームより前のフレーム画像をリカ―シブフィルタ処理するフィルタ手段と、
前記現フレームの画像と前記フィルタ手段で得られた画像との差分画像におけるノイズの統計量に基づいて現フレームの重み付け係数を得る係数取得手段と、
前記係数取得手段で得られた前記複数フレームに対応する重み付け係数に基づき前記複数フレームの動画像を加算する加算手段と、
を有する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像のノイズを低減する技術に関し、特に時間方向に画像を加算することでノイズを低減する技術に関する。
【背景技術】
【０００２】
従来、Ｘ線動画撮影の分野では、入射されたＸ線を蛍光体とイメージインテンシファイアにより可視光に変換し、ＣＣＤを用いたＴＶカメラで可視光像を撮像する２次元のＸ線動画撮影装置が提案されている。また、近年では撮像系をフラットパネルディテクタ（ＦＰＤ）で置き換えたＸ線動画撮影装置も広く普及している。
【０００３】
ところで、上述のＸ線動画撮影装置で得られる動画像（連続的に収集された複数の静止画像のセットのことであり、以下、個々の静止画像をフレームと呼ぶ）には、Ｘ線量子の揺らぎよる量子ノイズや検出器や回路等から発生するシステムノイズなど、様々なノイズが一連の撮影過程において重畳する。特に、医療用のＸ線動画撮影においては被検者への被曝の観点から少ないＸ線量で撮影を行うため、このようなノイズの影響は無視できず、画質を向上させるためにはノイズを低減することが重要となってくる。
【０００４】
そこで、画像処理の分野においては、従来から時間方向もしくは空間方向に画素の値を重み付け加算（所謂フィルタリング）することによって、上述のようなノイズを低減することが行われている。
【０００５】
具体的には、時間方向の重み付け加算は、対象となるフレームと過去のフレームの同一画素の値を平均化することで、時間的に略無相関なノイズの揺らぎを減少させる。また、空間方向の重み付け加算は、対象となる画素とその周辺の画素の値を平均化することで空間的に略無相関なノイズの揺らぎを減少させる。
【０００６】
なお、時間方向の重み付け加算は、フレーム間において被検者の動きによる時間的な信号変化がない場合には理想的な方法である。しかしながら、信号変化がある場合においては、過去のフレームの像が対象のフレーム上に残像として発生するという課題がある。また、空間方向の重み付け加算は、過去のフレームを参照しないため、残像は発生しないが、エッジの境界等の空間的に類似しない構造をノイズと共に不鮮明にするため、効果はあるものの好ましい方法とは言えない。
【０００７】
そこで、時間方向と空間方向を上手く組み合わせた方法も提案されており、例えば特許文献１の方法がある。この方法は、対象となるフレームと過去のフレームの局所平均の差に基づいて動きを判定し、動きが大きい場合は時間方向から空間方向の重み付け加算に切り替えるものである。これにより、動きがない場合は理想的なノイズ低減を行うことができ、動きが有った場合においても、残像を発生させないノイズ低減を行うことができる。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開平６−４７０３６号公報
【発明の概要】
【発明が解決しようとする課題】
【０００９】
ところで、上述のように時間方向と空間方法の重み付け加算を切り替える場合、ノイズが重畳した信号から動きをいかに精度良く判定するかが重要である。しかしながら、特許文献１の方法では、重畳したノイズの影響が十分に考慮されていない。具体的には、特許文献１の方法は、対象となるフレームと過去のフレームの局所平均の差に基づいて動きを判定しているが、その差が、被検者の動きによる時間的な信号変化であるのか、ノイズの時間的な揺らぎによる信号変化であるかを厳密に判定するものではないため、動きを精度良く判定できないという課題がある。
【００１０】
そこで、本発明の例示的な目的は、上記の課題を解決するためになされたもので、ノイズの影響を考慮することで動画像から動きを精度良く判定し、その判定結果に応じた好適なノイズ低減を行うことにある。
【課題を解決するための手段】
【００１１】
上記目的を達成するために、本発明による画像処理装置は複数フレームの動画像を時系列に重み付け加算する画像処理装置であって、
前記複数フレームの動画像の中で、現フレームより前のフレーム画像をリカ―シブフィルタ処理するフィルタ手段と、
前記現フレームの画像と前記フィルタ手段で得られた画像との差分画像におけるノイズの統計量に基づいて現フレームの重み付け係数を得る係数取得手段と、
前記係数取得手段で得られた前記複数フレームに対応する重み付け係数に基づき前記複数フレームの動画像を加算する加算手段と、
を有することを徴とする。
【００１２】
上記目的を達成するために、本発明による画像処理方法は、複数フレームの動画像を時系列に重み付け加算する画像処理方法であって、
前記複数フレームの動画像の中で、現フレームより前のフレーム画像をリカ―シブフィルタ処理するフィルタ工程と、
前記現フレームの画像と前記フィルタ工程で得られた画像との差分画像におけるノイズの統計量に基づいて現フレームの重み付け係数を得る係数取得工程と、
前記係数取得手段で得られた前記複数フレームに対応する重み付け係数に基づき前記複数フレームの動画像を加算する加算工程と、
を有することを徴とする。
【００１３】
また、本発明は、上述した画像処理装置、および画像処理方法をコンピュータに実行させるためのプログラム、並びにプログラムを記録するコンピュータ読み取り可能な記録媒体を含む。
【発明の効果】
【００１４】
ノイズの特性に応じた好適なノイズ低減を行うことができる。
【図面の簡単な説明】
【００１５】
【図１】実施例１および２によるＸ線動画撮影装置全体の構成図である。
【図２】実施例１による画像処理部１１２の処理手順を示すフローチャートである。
【図３】実施例２による画像処理部１１２の処理手順を示すフローチャートである。
【図４】実施例３および４によるＸ線動画撮影装置全体の構成図である。
【図５】実施例３による画像処理部１１２の処理手順を示すフローチャートである。
【図６】実施例４による画像処理部１１２の処理手順を示すフローチャートである。
【図７】ノイズの算出方法を説明する図である。
【図８】ピラミッドアルゴリズムを説明する図である。
【図９】階層画像の周波数特性を説明する図である。
【発明を実施するための形態】
【００１６】
［実施例１］
本発明は、例えば図１に示すようなＸ線動画撮影装置１００に適用される。図１に示すように、Ｘ線動画撮影装置１００は、撮影された動画像をモニタ上に出力する際の効果的な画像処理を行う機能を有するＸ線の動画撮影装置である。Ｘ線動画撮影装置１００は、Ｘ線発生部１０１、Ｘ線検出器１０４、データ収集部１０５、前処理部１０６、ＣＰＵ１０８、メインメモリ１０９、操作部１１０、表示部１１１、画像処理部１１２を備えており、これらはＣＰＵバス１０７を介して互いにデータ授受が可能に接続されている。
【００１７】
上述のような画撮影装置１００において、まず、メインメモリ１０９は、ＣＰＵ１０８での処理に必要な各種のデータなどを記憶すると共に、ＣＰＵ１０８のワーキング・メモリとして機能する。ＣＰＵ１０８は、メインメモリ１０９を用いて、操作部１１０からの操作にしたがった装置全体の動作制御等を行う。これによりＸ線動画撮影装置１００は、以下のように動作する。
【００１８】
まず、操作部１１０を介してユーザから撮影指示が入力されると、この撮影指示はＣＰＵ１０８によりデータ収集部１０５に伝えられる。ＣＰＵ１０８は、撮影指示を受けると、Ｘ線発生部１０１及びＸ線検出器１０４を制御してＸ線動画撮影を実行させる。
【００１９】
Ｘ線動画撮影では、まずＸ線発生部１０１が、被検体１０３に対してＸ線ビーム１０２を照射する。Ｘ線発生部１０１から照射されたＸ線ビーム１０２は、被検体１０３を減衰しながら透過して、Ｘ線検出器１０４に到達する。そしてＸ線検出器１０４によりＸ線画像信号が連続的に出力される。なお、本実施の形態では被検体１０３を人体とする。すなわち、Ｘ線検出器１０４から出力されるＸ線動画像は人体画像となる。
【００２０】
データ収集部１０５は、Ｘ線検出器１０４から連続的に出力されたフレーム単位のＸ線画像信号を所定のデジタル信号に変換してＸ線画像データとして前処理部１０６に供給する。前処理部１０６は、データ収集部１０５からの信号（Ｘ線画像データ）に対して、オフセット補正やゲイン補正等の前処理を行う。この前処理部１０６で前処理が行われたＸ線画像データは、ＣＰＵ１０８の制御により、ＣＰＵバス１０７を介して、メインメモリ１０９、画像処理部１１２に順次転送される。
【００２１】
画像処理部１１２は、複数フレームの動画像を時系列に重み付け加算することで動画像のノイズを低減するものである。複数のフレームで構成される動画像を時間方向にリカ―シブフィルタ処理するフィルタ部１１３と、現フレームの画像データからノイズを取得するノイズ取得部１１４と有する。また、ノイズ取得部１１４で取得したノイズの統計量を取得する統計量取得部１１５と、統計量取得部１１５で取得したノイズの統計量から予め取得されているノイズの分布に従うかを示す評価値を得る検定部１１６と、
この評価値で定められる重み係数を、複数のフレームのそれぞれに対応して取得する図示しない係数取得部１１６１と、を有する。そして、現フレームよりも以前に処理した画像データを図示しないメモリ１１６２から取得し、現フレームの画像データを係数取得部１１６１で取得した重み係数で加算しメモリ１１６２に格納する加算部１１７を備えている。また、これらの各構成部は、ＣＰＵバスに接続されている。
【００２２】
以上のような構成を備えたＸ線動画撮影装置１００において、画像処理部１１２の動作について、図２に示すフローチャートを用いて具体的に説明する。
【００２３】
上述のごとく前処理部１０６によって得られた複数のフレームは、ＣＰＵバス１０７を介して画像処理部１１２に転送されるとともに、処理を開始する現フレームＦｔ（ｔは撮影開始から何枚目のフレームであるかを表す変数であり１以上かつ撮影枚数以下の整数である）としてｔ＝２を設定する。
【００２４】
ノイズ取得部１１４は、フィードバック係数ａを設定する（ｓ２０１）。ここで、フィードバック係数ａは０＜ａ＜１の範囲で任意に設定すれば良く、本実施の形態ではユーザが操作部１１０を介して設定するものとする。
【００２５】
次に、フィルタ部１１３において、設定された現フレームＦｔよりも前のフレームとフィードバック係数ａを用いて、すべての画素に対して時間方向の重み付け加算を行う（ｓ２０２）。具体的には、下記式のように前フレームの処理結果を順次フィードバック加算する、いわゆるリカーシブフィルタ（再帰型）によって時間方向の重み付け加算を行う。
【００２６】
【数１】
ここで、Ｆｒｔ−１が現フレームＦｔよりも前のフレーム画像、すなわち１〜ｔ−１枚目のフレームを時間方向に重み付け加算したデータとなる。なお、本実施の形態では上記式のような再帰型フィルタを用いたがこれに限定されるものではなく、異なる再帰型フィルタを用いてもよい、また非再帰型フィルタを用いても良い。
【００２７】
次に、ノイズ取得部１１４において、現フレームＦｔから時間方向の重み付け加算がなされた前フレームＦｒｔ−１を減算することで差分画像Ｄｔを生成する（ｓ２０３）。具体的には、下記式のように、すべての画素に対してそれぞれ単純な減算を行う。
【００２８】
【数２】
【００２９】
次に、統計量取得部１１５において、差分画像Ｄｔに重畳されるノイズの分散σＤ２を差分画像のノイズの統計量として取得する（ｓ２０４）。
【００３０】
ここで、現フレームＦｔに重畳されるノイズの分散をσＦ２とし、時間方向の重み付け加算がなされた前フレームＦｒｔ−１に重畳されるノイズの分散をσＦｒ２とし、ノイズが時間的に無相関であると仮定すれば分散の加法性により差分画像Ｄｔに重畳されるノイズの分散σＤ２は下記式で表される。
【００３１】
【数３】
【００３２】
また、ノイズの分布（例えば分散の確率分布）が時間的に変化しないと仮定すれば、上記式の右辺第二項は、現フレームＦｔに重畳されるノイズの分散σＦ２および、時間方向の重み付け加算によるノイズ減衰率を考慮して下記式のように変形可能である。
【００３３】
【数４】
ここで、ａはユーザが操作部１１０を介して設定したフィードバック係数である。
【００３４】
なお、現フレームＦｔに重畳されるノイズの分散σＦ２は、検出器への入射Ｘ線量（入射Ｘ線量子数）に比例する量子ノイズを含むため、実際には入射Ｘ線量に依存した値となる。具体的には、検出器への入射Ｘ線量に比例する量子ノイズをσＱ２とし、入射Ｘ線量に依存しないシステムノイズをσＳ２とし、両者が無相関であると仮定すれば、現フレームＦｔに重畳されるノイズの統計量である分散σＦ２は、下記式で表される。
【００３５】
【数５】
ただし、Ｖは現フレームＦｔの画素値であり、ｇは入射Ｘ線量に比例した画素値Ｖを入射Ｘ線量相当の値に変換するためのゲイン係数である。
【００３６】
ここで、上記式を考慮すれば、差分画像Ｄｔに重畳されるノイズの分散σＤ２は下記式となる。
【００３７】
【数６】
従って、実際には現フレームＦｔの画素値Ｖに応じて分散σＤ２を推定する必要がある。そこで、本実施の形態では、対象画素（ｘ，ｙ）の画素値Ｆｔ（ｘ，ｙ）または、対象画素（ｘ，ｙ）を中心とした任意のｎ×ｎ画素の平均画素値をＶ として分散σＤ２（Ｖ）をすべての画素で推定する。
【００３８】
なお、上記式の量子ノイズσＱ２、システムノイズσＳ２、ゲイン係数ｇのパラメータは、撮影過程で一意に決まる値であり、予め取得しておくことが可能である。従って、予め取得したパラメータをメインメモリ１０９に保存しておき、その結果をロードして用いればよい。
【００３９】
ここで、上記パラメータの算出方法については、特に限定するものではないが、例えば、製品検査時に被検体のない状態で検出器への入射Ｘ線量を変えたｎ回の撮影を行い、図７のように各撮影画像の平均画素値Ｖと分散σ２の関係を求めた離散データ｛（Ｖｉ，σｉ２）｜ｉ＝１，２，．．．，ｎ｝を得、この離散データを良く近似するパラメータを算出する。具体的には、下記式で表される残差の２乗和Ｊが、最小となるパラメータを最小二乗法によって求め、このパラメータを予めメインメモリ１０９に保存しておけば良い。
【００４０】
【数７】
【００４１】
次に、検定部１１６において、統計量取得部１１５で求めた差分画像Ｄｔのノイズの統計量が予めえられたノイズの統計分（例えば分散値の確率分布）に従うか否かを表す評価値としての有意確率を求める（ｓ２０５）。ここで、有意確率とは、仮説検定において、ある仮説（帰無仮説）のもとで得られた検定統計量が実現する確率であり、この有意確率が小さい場合は、確率的に仮説が成り立つ可能性が低いことを表すものである。そこで、本実施の形態では「差分画像Ｄｔのノイズの統計量が予め得られたノイズの確率分布Ｎ（μ，σ２）に従う」という帰無仮説に対して、その有意確率Ｐを求める。
【００４２】
なお、「差分画像Ｄｔのノイズの統計量がノイズの確率分布Ｎ（μ，σ２）に従う」ということは、差分画像Ｄｔにはノイズ以外の信号成分が含まれないことを意味し、言い換えれば「差分したフレーム間には被検者の動きによる信号変化がない」ということである。よって、求めた有意確率Ｐは、差分したフレーム間に被検体の動きによる信号変化がないかを確率的に表したものとなる。なお、ノイズの統計量としてノイズの平均値、標準偏差、分散の少なくともいずれかを用いることがえきる。
【００４３】
なお、上述のごとく差分画像Ｄｔに重畳するノイズは画素毎に異なる確率分布（分散）を持つ。そこで、本実施の形態では対象画素（ｘ，ｙ）を中心とした任意のｎ×ｎ画素の範囲を標本集合Ωとし、「標本集合Ωがノイズの確率分布Ｎ（μ，σ２）に従う」という帰無仮説に対して、すべての画素で有意確率Ｐを算出する。ここで、ノイズの母平均μは０であり、母分散σ２は統計量取得部１１５で求めた対象画素（ｘ，ｙ）に対する分散σＤ２（Ｖ）を推定値として用いる。また、ｎは任意のサイズを設定すれば良く、例えば本実施の形態ではｎ＝３とする。
【００４４】
なお、有意確率Ｐの算出方法に関しては、特に限定するものではないが、本実施の形態では、母平均μおよび母分散σ２に対する検定統計量から求めた有意確率を統合することで、上記の帰無仮説に対する有意確率Ｐを算出する。
【００４５】
まず、対象画素（ｘ，ｙ）を中心とした任意のｎ×ｎ画素の範囲の標本集合Ωに対し、母平均μに対する検定統計量Ｔを下記式にて求める。
【００４６】
【数８】
ただし、Ｅは標本平均であり、Ｓ２は不偏分散を表す。また、母平均μは０である。
【００４７】
ここで、上記式にて求めた検定統計量Ｔは、母平均がμの場合に自由度２ｎ−１のＴ分布に従う。よって、「母平均がμである」という帰無仮説おける両側検定（「母平均がμでない」を対立仮説とする検定）に対する有意確率ＰＴは、下記式にて求めることができる。
【００４８】
【数９】
ただし、Γ（・）はガンマ関数を表す。
【００４９】
なお、上記式で求めた有意確率ＰＴは、「標本集合Ωの母平均が０である」という帰無仮説が成り立つ可能性を確率的に求めたものであり、例えばＰＴ＝０．０１であれば、帰無仮説は１％の確率で正しいことになる。言い換えれば、帰無仮説は９９％の確率で誤りであり、「標本集合Ωの母平均が０でない」という対立仮説が成り立つ可能性が非常に高いこと表している。
【００５０】
次に、対象画素（ｘ，ｙ）を中心とした任意のｎ×ｎ画素の範囲の標本集合Ωに対し、母分散σ２に対する検定統計量χ２を下記式にて求める。
【００５１】
【数１０】
ただし、Ｓ２は不偏分散を表す。また、母分散σ２は統計量取得部１１５で求めた対象画素（ｘ，ｙ）に対する分散σＤ２（Ｖ）である。
【００５２】
ここで、上記式にて求めた検定統計量χ２は、母分散がσ２の場合に自由度２ｎ−１のχ２分布に従う。よって、「母分散がσ２である」という帰無仮説における両側検定（「母分散がσ２でない」を対立仮説とする検定）に対する有意確率Ｐχは、下記式にて求めることができる。
【００５３】
【数１１】
ただし、Γ（・）はガンマ関数を表す。
【００５４】
なお、上記式で求めた有意確率Ｐχは、「標本集合Ωの母分散がσＤ２（Ｖ）である」という帰無仮説が成り立つ可能性を確率的に求めたものであり、例えばＰＴ＝０．０１であれば、帰無仮説は１％の確率で正しいことになる。言い換えれば、帰無仮説は９９％の確率で誤りであり、「標本集合Ωの母分散がσＤ２（Ｖ）でない」という対立仮説が成り立つ可能性が非常に高いことを表している。
【００５５】
以上により求めた有意確率ＰＴおよびＰχは、それぞれ標本集合Ωの母平均および母分散に対する有意確率を独立に求めたものであり、本実施の形態ではその両者を満たす、すなわち「標本集合Ωがノイズの確率分布Ｎ（μ，σ２）に従う」という帰無仮説に対して、その有意確率Ｐを求めたい。そこで、求めた２つの有意確率ＰＴおよびＰχを統合し、上記帰無仮説に対する有意確率Ｐを算出する。なお、有意確率の統合方法に関しては特に限定するものではないが、本実施の形態では、Ｆｉｓｈｅｒの方法を用いることとする。具体的には下記式にて有意確率を統合した検定統計量Ｆを求める。
【００５６】
【数１２】
ここで、上記式で統合した検定統計量Ｆは、独立に求めた２つの有意確率ＰＴおよびＰχを統合したものであり、これは自由度４のχ２分布に従う。また、求めた検定統計量Ｆは、２つの有意確率ＰＴおよびＰχが小さくなるほど大きな値をとる。よって、下記式のように、片側検定（上側）に対する有意確率をＰとして算出する。
【００５７】
【数１３】
以上により、対象画素（ｘ，ｙ）に対する有意確率Ｐを算出することができる。ここで、同様の処理をすべての画素に対して実行する。なお、上記式で求めた有意確率Ｐは、「標本集合Ωがノイズの確率分布Ｎ｛０， σＤ２（Ｖ）｝に従う」、すなわち「標本集合Ωには被検者の動きによる信号変化がない」という帰無仮説が成り立つ可能性を確率的に求めたものであり、例えばＰ＝０．０１であれば、帰無仮説は１％の確率で正しいことになる。言い換えれば、帰無仮説は９９％の確率で誤りであり、「標本集合Ωには被検者の動きによる信号変化が有る」という対立仮説が成り立つ可能性が非常に高いことを表している。
【００５８】
次に、加算部１１７において、ｓ２０６〜ｓ２０８のステップを実行することで、有意確率に基づいた時空間方向の重み付け加算を行う。まず、画素毎に求めた有意確率Ｐに従って時間方向の重みづけ係数としてフィードバック係数ａｍを設定する（ｓ２０６）。ここで、上述のごとく有意確率Ｐは「差分したフレーム間には被検者の動きによる信号変化がない」という帰無仮説が成り立つ可能性を確率的に求めたものであり、有意確率Ｐが小さいほど「差分したフレーム間には被検者の動きによる信号変化が有る」という対立仮説が成り立つ可能性が高くなる。よって、このような場合に前フレームの処理結果に対するフィードバック係数を大きく設定すると残像が発生する可能性がある。そこで、下記式のように有意確率Ｐが小さいほどフィードバック係数ａｍが小さくなるように設定する。
【００５９】
【数１４】
ここで、ａは、ユーザが操作部１１０を介して設定したフィードバック係数であり、この値を基準として実際に用いるフィードバック係数ａｍを設定する。
【００６０】
なお、フィードバック係数ａｍの設定方法はこれに限定されるものではなく、例えば下記式のように有意確率Ｐが設定した有意水準αよりも小さい場合、すなわち仮説検定で言う所の対立仮説が採択される場合は０（時間方向の重み付け加算を実行しない）に設定し、有意確率Ｐが設定した有意水準α以上の場合、すなわち仮説検定で言う所の帰無仮説を棄却できない場合はａに設定しても良い。
【００６１】
【数１５】
ここで、有意水準αは仮説検定で一般的に使われる値を用いれば良く、本実施の形態では例えばα＝０．０１とする。
【００６２】
以上、すべての画素に対して上述の処理を実行し、有意確率Ｐに基づいたフィードバック係数ａｍを設定する。なお、このように設定したフィードバック係数ａｍは、各画素の近傍における動きに基づいたものであり、局所的な動きに応じて最適なフィードバック係数が得られるものである。
【００６３】
次に、図示しない空間フィルタ部１８０は、現フレームＦｔのすべての画素に対して空間方向の重み付け加算を空間フィルタ処理として行う（ｓ２０７）。なお、空間方向の重み付け加算は、ＦＩＲ等の単純な空間フィルタで行うこともできるが、エッジの境界等の空間的に類似しない構造をノイズと共に不鮮明にする。そこで、本実施の形態では下記式で表されるεフィルタを空間フィルタとして用いることで、構造のぼけを抑えた重み付け加算を行う。
【００６４】
【数１６】
ただし、ｈはサイズがＮ×Ｎのローパスフィルタであり、本実施の形態では例えば９×９の平均化フィルタを用いる。
【００６５】
なお、上述の処理は、対象画素（ｘ，ｙ）と近傍画素の画素値の差を類似度の指標値とし、この指標値がε値よりも大きいほど近傍画素に対するローパスフィルタの重みを小さくするものである。これにより、類似しない近傍画素を除外した空間方向の重み付け加算がなされ、構造のぼけを抑えることができる。
【００６６】
ここで、ε値は対象画素（ｘ，ｙ）と近傍画素が類似しているか否かを判定する閾値であり、本実施の形態では、下記式にてε値を設定する。
【００６７】
【数１７】
ただし、σ２は現フレームＦｔに重畳するノイズの分散であり、統計量取得部１１５で求めた対象画素（ｘ，ｙ）に対する分散σＦ２（Ｖ）を用いて、各画素でε値を設定する。
【００６８】
なお、上記式の意味するところは、もし対象画素（ｘ，ｙ）と近傍画素がノイズを除けば同じ画素値であると仮定すると、その差はノイズの揺らぎによる信号変化であり、ノイズが正規分布ならば約９９．７％の確率で±３√２σ２の範囲となる。よって、この値をε値として設定することで、もし対象画素（ｘ，ｙ）と近傍画素の画素値の差が想定されるノイズの揺らぎよりも大きい場合、その変化はノイズの揺らぎによる信号変化ではない。すなわち、対象画素（ｘ，ｙ）と近傍画素は類似しない画素であると判定するものである。
【００６９】
以上、すべての画素に対して上述の処理を実行し、ε値に基づいた空間方向の重み付け加算を行う。なお、本実施の形態ではεフィルタを用いたが、これに限定されるものではなくＢｉｌａｔｅｒａｌフィルタ等の類似の非線形フィルタを用いることもできる。
【００７０】
次に、空間方向に重み付け加算した現フレームＦｓｔとｓ２０６のステップで求めたフィードバック係数ａｍを用いて、時空間方向に重み付け加算した現フレームＦｏｔを生成する（ｓ２０８）。具体的には、すべての画素に対して下記式の処理を実行する。
【００７１】
【数１８】
ここで、Ｆｒｔ−１は時間方向の重み付け加算がなされた前フレームであり、加算部１１７で生成したものである。また、ｂは空間方向の重み付け加算の寄与率を表し、フィードバック係数ａｍが小さいほど大きくなるように設定する。すなわち、残像が発生する可能性がある画素に対しては、時間方向の寄与を小さくする代わりに、空間方向の寄与を大きくし空間フィルタ処理の処理効果を上げ、十分なノイズ低減効果を得るものである。
【００７２】
なお、寄与率ｂの設定方法は特に限定するものではないが、本実施の形態では、上記式によるノイズ減衰率Ｒｏ（入力のノイズの分散に対する出力のノイズの分散の比）が、基準値Ｒと略同等となるように寄与率ｂを設定する。ここで、上記式によるノイズ減衰率Ｒｏは、空間方向の重み付け加算（上記式の右辺第二項の波括弧内）におけるノイズ減衰率をＲｓとすれば下記式で表される。
【００７３】
【数１９】
ただし、ａはユーザが操作部１１０を介して設定したフィードバック係数である。
従って、Ｒｏに基準値Ｒを代入し、Ｒｓについて解けば空間方向の重み付け加算が担うノイズ減衰率が求まる。よって、空間方向の重み付け加算によるノイズ減衰率がＲＳとなるような寄与率ｂを設定する。なお、この条件を満たす厳密解を求めるのは困難である。そこで、空間方向に重み付け加算した現フレームＦｓｔに重畳するノイズが無視可能なレベルと仮定し、下記式にて求まる近似値を寄与率ｂとして設定する。
【００７４】
【数２０】
なお、基準値Ｒは任意に設定すればよく、本実施の形態ではユーザが操作部１１０を介して設定したフィードバック係数ａを用いて下記式にて設定する。
【００７５】
【数２１】
ここで、上記式で設定される基準値Ｒは、１〜ｔ枚目のフレームをリカーシブフィルタによって処理した場合のノイズ減衰率そのものである。従って、このように設定した基準値Ｒを用いることで、すべての画素においてフィードバック係数ａを用いたリカーシブフィルタと略同等のノイズ低減効果が得られるものである。
【００７６】
なお、リカーシブフィルタのノイズ減衰率は時間依存性があり、ｔが大きくなるほど小さな値となる。よって、上記式で設定した基準値Ｒを用いるとフレーム間でノイズ低減効果が同等とならない。そこで、下記式のようにｔ→∞の場合のノイズ減衰率を基準値Ｒとして設定しても良い。この場合、すべてのフレームにおいて略同等のノイズ低減効果を得ることができる。
【００７７】
【数２２】
以上、すべての画素に対して上述の処理を実行することで、ノイズ低減がなされた現フレームＦｏｔを生成することができる。
【００７８】
次に、全フレームに対してノイズ低減が実行されたか否かを判定し（ｓ２０９）、全フレームに対してノイズ低減が実行されていない場合は、ｔをインクリメントすることで次のフレームを現フレームとして設定し（ｓ２１０）、ｓ２０２〜ｓ２０８の各ステップを実行する。全フレームに対してノイズ低減が実行された場合は、処理を終了する。
【００７９】
以上、実施例１ではノイズの確率分布に基づいて動きである可能性を確率的に求めることで、精度良く動きを判定することができる。また、動きである可能性が高い画素では時間方向のフィードバック係数を小さくすることで残像を抑え、かつ空間方向の寄与率をノイズ減衰率に基づいて設定することで、すべての画素で略同等のノイズ抑制効果を得ることができる。
【００８０】
［実施例２］
本発明は、Ｘ線動画撮影装置１００において、画像処理部１１２の動作を実施例１とは異なる図３のフローチャートに従った動作とする。なお、図３に示すフローチャートにおいて、図２に示したフローチャートと同様に処理実行するステップは同じ符号を付し、ここでは、上述した実施例１とは異なる構成についてのみ具体的に説明する。
【００８１】
まず、処理を開始する現フレームＦｔとしてｔ＝２を設定するとともに、複数のフィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝を設定する（ｓ３０１）。ここで、フィードバック係数はユーザが操作部１１０を介して設定したフィードバック係数ａを最大値とし、下記式のように０〜ａの範囲をｎ等分したフィードバック係数を設定する。
【００８２】
【数２３】
ここで、ｎは任意の数を設定すれば良く、例えばｎ＝１０とする。
【００８３】
次に、設定した複数のフィードバック係数の１つを選択し（ｓ３０２）、ｓ２０２〜ｓ２０５の各ステップを実施例１と同様に実行することで有意確率を求める。
【００８４】
以上の処理をすべてのフィードバック係数に対して実行するまで処理を繰り返し（ｓ３０３）、フィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝に対する有意確率｛Ｐｉ｜ｉ＝１，２，…ｎ｝を求める。
【００８５】
次に、加算部１１７において、複数のフィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝の中から最適なフィードバック係数ａｍを設定する（ｓ３０４）。ここで、上述した通り有意確率は「差分したフレーム間には被検者の動きによる信号変化がない」という帰無仮説が成り立つ可能性を確率的に求めたものであり、有意確率が小さいほど「差分したフレーム間には被検者の動きによる信号変化が有る」という対立仮説が成り立つ可能性が高くなる。そこで、本実施の形態では有意確率が最も大きくなるフィードバック係数、すなわち対立仮説が成り立つ可能性が最も低いフィードバック係数を最適なフィードバック係数ａｍとして設定する。ただし、最大の有意確率が有意水準αよりも小さい場合、すなわち対立仮説がすべてのフィードバック係数において採択される場合は、すべてのフィードバック係数を棄却し、ａｍ＝０（時間方向の重み付け加算を実行しない）に設定する。
【００８６】
なお、フィードバック係数ａｍの設定方法はこれに限定されるものではなく、例えば、有意確率が有意水準αよりも大きいフィードバック係数の内、最大のものを最適なフィードバック係数ａｍとして設定しても良い。この場合は、設定した有意水準αを満たすものも内、最大のノイズ抑制効果を得られるフィードバック係数が選択されることとなる。
【００８７】
以上、すべての画素に対して上述の処理を実行し、複数の有意確率｛Ｐｉ｜ｉ＝１，２，…ｎ｝に基づいたフィードバック係数ａｍを設定する。
【００８８】
次に、現フレームＦｔのすべての画素に対して空間方向の重み付け加算を行い（ｓ２０７）、空間方向に重み付け加算した現フレームＦｓｔとｓ３０４のステップで求めたフィードバック係数ａｍを用いて、時空間方向に重み付け加算した現フレームＦｏｔを生成する（ｓ３０５）。具体的には、すべての画素に対して下記式を実行する。
【００８９】
【数２４】
ここで、寄与率ｂは実施例１と同様に、上記式によるノイズ減衰率Ｒｏに基づいて設定すれば良い。
【００９０】
なお、上記式ではａｍ＝０の場合は時間方向の重み付け加算が行われず、空間方向の重み付け加算のみが行われる。そのため、検定部１１６において最適なフィードバック係数としてａｍ＝０が設定された場合（すべてのフィードバック係数を棄却しａｍ＝０を設定した場合は除く）、１つ前のフレームＦｔ−１との間に被検者の動きによる信号変化がないにも関わらず、１つ前のフレームＦｔ−１の情報が利用されないことになる。そこで、この場合は、例外的に下記式にて処理を行っても良い。
【００９１】
【数２５】
ここで、ａはユーザが操作部１１０を介して設定したフィードバック係数である。また、Ｆｏｔ−１は時空間方向に重み付け加算された前フレームである。
【００９２】
上記式は、空間方向の重み付け加算のみを実行する代わりに、時空間方向に重み付け加算された前フレームＦｏｔ−１とを重み付け加算するものである。この場合、時間方向の情報も利用した処理となるため、エッジの境界等の構造のぼけを抑えることができる。
【００９３】
次に、全フレームに対してノイズ低減が実行されたか否かを判定し（ｓ２０９）、全フレームに対してノイズ低減が実行されていない場合は、ｔをインクリメントすることで次のフレームを現フレームとして設定し（ｓ２１０）、ｓ３０２〜ｓ３０５の各ステップを実行する。全フレームに対してノイズ低減が実行された場合は、処理を終了する。
【００９４】
以上、実施例２では複数の有意確率に基づき確率的に最もリスクの少ないフィードバック係数、または確率的にリスクが少ないものの内、最もノイズ抑制効果の高いフィードバック係数を選択するものであり、実施例１よりも最適なフィードバック係数を設定できる。
【００９５】
［実施例３］
本発明は、例えば図４に示すようなＸ線動画撮影装置４００に適用される。このＸ線動画撮影装置４００は、Ｘ線動画撮影装置１００に対し、分解部４０１、再構成部４０２を備える構成としている。また、本発明は画像処理部１１２の処理手順を実施例１とは異なる図５に示したフローチャートに従った動作とする。
【００９６】
なお、図４のＸ線動画撮影装置４００において、図１のＸ線動画撮影装置１００と同様に動作する箇所は同じ符号を付し、その詳細は省略する。また、図５に示すフローチャートにおいて、図２に示したフローチャートと同様に処理実行するステップは同じ符号を付し、ここでは、上述した実施例１とは異なる構成についてのみ具体的に説明する。
【００９７】
まず、実施例１と同様に処理を開始する現フレームＦｔとしてｔ＝２を設定するとともに、フィードバック係数ａを設定する（ｓ２０１）。
【００９８】
次に、分解部４０１において、現フレームＦｔを複数の異なる解像度を持つ階層画像に分解する（ｓ５０１）。なお、本実施の形態ではＢｕｒｔとＡｄｅｌｓｏｎによるピラミッドアルゴリズムを用いて階層画像に分解する。
【００９９】
ここで、ピラミッドアルゴリズムの分解について図８の（ａ）及び（ｂ）を参照して説明する。図８（ａ）は、入力画像Ｇｉに対して１レベルの分解を表しており、まず入力画像Ｇｉは図８（ｂ）に示すフィルタでフィルタリングされ（８０）、さらに１／２倍のダウンサンプリング８１によって、解像度が縦横ともに１／２のガウシアン画像Ｇｉ＋１を生成する。また、ガウシアン画像Ｇｉ＋１は、さらに２倍のアップサンプリング８２によって０値が挿入された後、図８（ｂ）に示すフィルタ係数に４倍のゲインを掛けたフィルタによってフィルタリングされ（８３）、解像度が入力画像Ｇｉと同一の画像を生成する。そして、入力画像Ｇｉからこの画像を減算器８４にて減算することでラプラシアン画像Ｌｉを生成する。以上の動作により、入力画像Ｇｉをガウシアン画像Ｇｉ＋１、ラプラシアン画像Ｌｉに分解することができる。さらに、マルチレベルの分解を行う場合は、出力されたガウシアン画像Ｇｉ＋１を入力画像として同一の処理を繰り返せばよい。
【０１００】
なお、本実施の形態では、現フレームＦｔを最初の入力画像とし、上述の処理をＭ回繰り返すことで、解像度の異なる複数のガウシアン画像｛Ｇｔ，ｉ｜ｉ＝０，１，…Ｍ｝およびラプラシアン画像｛Ｌｔ，ｉ｜ｉ＝０，１，…Ｍ−１｝を生成する。ここで、このように生成された画像は図９に示す周波数特性を持つものである。図９のようにガウシアン画像はレベルｉが大きくなるほど、入力画像のより低い周波数のみで構成される画像となる。また、ラプラシアン画像は入力画像の一部の周波数帯域で構成される画像となり、レベルｉが大きくなるほど低い周波数帯域を抽出したものとなる。
【０１０１】
次に、フィルタ部１１３において、設定された現フレームよりも前のガウシアン画像とフィードバック係数ａを用いて、時間方向の重み付け加算をおこなう（ｓ５０２）。具体的には、下記式にてレベルＭを除くすべての画素に対しリカーシブフィルタを用いた重み付け加算を行う。
【０１０２】
【数２６】
【０１０３】
次に、ノイズ取得部１１４において、現フレームのガウシアン画像Ｇｔ，ｉから時間方向の重み付け加算がなされた前フレームのガウシアン画像Ｇｒｔ−１，ｉを減算することで差分画像Ｄｇｔ，ｉを生成する（ｓ５０３）。具体的には、下記式にてレベルＭを除くすべての画素に対してそれぞれ単純な減算を行う。
【０１０４】
【数２７】
【０１０５】
次に、統計量取得部１１５において、差分画像Ｄｇｔ，ｉに重畳されるノイズの分散σｄｉ２を推定する（ｓ５０４）。ここで、現フレームのガウシアン画像Ｇｔ，ｉに重畳されるノイズの分散をσｇｉ２とし、時間方向の重み付け加算によるノイズ減衰率を考慮すれば、差分画像Ｄｇｔ，ｉに重畳されるノイズの分散σｄｉ２は下記式で表される。
【０１０６】
【数２８】
ここで、ａはユーザが操作部１１０を介して設定したフィードバック係数である。
【０１０７】
なお、現フレームのガウシアン画像Ｇｔ，ｉに重畳されるノイズの分散σｇｉ２は、実施例１と同様にして求まる現フレームＦｔに重畳されるノイズの分散σＦ２（Ｖ）および分解過程におけるノイズ減衰率を考慮すれば下記式となる。
【０１０８】
【数２９】
ただし、ｈは図８（ｂ）に示したフィルタ行列であり、（↑ａ） ・はａ倍のアップサンプリングを表す。また、＊は畳み込みを表す。
【０１０９】
ここで、上記式はレベルｉおよび画素値Ｖに依存した値をとる。従って、差分画像Ｄｇｔ，ｉに重畳されるノイズの分散σｄｉ２も同じく、レベルｉおよび画素値Ｖに依存した値をとなる。なお、画素値Ｖは、現フレームＦｔの重み付き平均であるガウシアン画像Ｇｔ，ｉを用いれば良いので、実際の処理としては、レベルｉのガウシアン画像Ｇｔ，ｉの画素値Ｇｔ，ｉ （ｘ，ｙ）に基づいてレベルＭを除くすべての画素に対する推定値を算出すれば良い。
【０１１０】
次に、検定部１１６において、差分画像Ｄｇｔ，ｉが統計量取得部１１５で求めたノイズの統計量に基づいて確率分布に従うか否かを表す有意確率を求める（ｓ５０５）。ここで、上述のごとく差分画像Ｄｇｔ，ｉに重畳するノイズはレベル及び画素毎に異なる確率分布（分散）を持つ。そこで、本実施の形態ではレベルｉの対象画素（ｘ，ｙ）を中心とした任意のｎ×ｎ画素の範囲を標本集合Ωとし、「標本集合Ωがノイズの確率分布Ｎ（μ，σ２）に従う」という帰無仮説に対して、各レベルｉのすべての画素で有意確率Ｐｉを算出する。ここで、ノイズの母平均μは０であり、母分散σ２は統計量取得部１１５で求めたレベルｉの対象画素（ｘ，ｙ）に対する分散σｄｉ２（Ｖ）を推定値として用いる。また、ｎは任意のサイズを設定すれば良く、例えば本実施の形態ではｎ＝３とする。
【０１１１】
なお、標本集合Ωに対する有意確率の算出方法は、母平均μに対する検定統計量Ｔの算出方法を下記式に変更する以外は、実施例１と同様に求めることができる。
【０１１２】
【数３０】
ここで、実施例１と上記式の違いはＣであり、これは空間的に相関のあるノイズを持つガウシアン画像Ｇｔ，ｉから求めた検定統計量Ｔを自由度２ｎ−１のＴ分布に従わせるための補正係数である。具体的には、空間的なノイズの相関を考慮して下記式にて算出する。
【０１１３】
【数３１】
ただし、ｈは図８（ｂ）に示したフィルタ行列であり、（↑ａ） ・はａ倍のアップサンプリングを表す。また、＊は畳み込みを表す。
【０１１４】
なお、ｍは標本サイズと同じサイズの平均化フィルタ行列であり、対象画素（ｘ，ｙ）を中心とした３×３画素の範囲を標本集合Ωとする本実施の形態では、下記式となる。
【０１１５】
【数３２】
以上により求めた母平均μに対する検定統計量Ｔと、実施例１と同様に求めた母分散σ２に対する検定統計量χ２を統合することで、レベルＭを除くすべての画素に対する有意確率Ｐｉを求める。なお、このように算出した有意確率Ｐｉは、画像を複数の周波数帯域に分解し、各々の周波数帯域における動きを確率的に求めたものである。そのため、周波数毎に異なる変化を伴うような非剛体運動においても精度良く動きを判定できるものである。
【０１１６】
次に、加算部１１７において、レベルＭを除くすべての画素に対するフィードバック係数ａｍｉを設定する（ｓ２０６）。ここでは、各レベルｉに対して実施例１と同様にフィードバック係数ａｍｉを設定すれば良い。
【０１１７】
次に、現フレームのラプラシアン画像Ｌｔ，ｉのすべての画素に対して空間方向の重み付け加算を行う（ｓ５０６）。なお、本実施の形態では、ラプラシアン画像Ｌｔ，ｉと同じレベルのガウシアン画像Ｇｔ，ｉからεフィルタと同様の原理にてノイズを抽出し、抽出したノイズを現フレームのラプラシアン画像Ｌｔ，ｉから減算することで、構造のぼけを抑えた重み付け加算を行う。具体的には下記式にて処理を行う。
【０１１８】
【数３３】
ただし、ｈは図８（ｂ）に示したフィルタ行列である。
【０１１９】
上記式において、右辺第二項がε値を基準としてラプラシアン画像Ｌｔ，ｉと同一の周波数帯域からノイズを推定したものとなる。ここで、ε値は、対象画素（ｘ，ｙ）と近傍画素がノイズを除けば同じ画素値であるかを判定する閾値であり、現フレームＦｔに重畳されるノイズの分散σＦ２（Ｖ）および分解過程におけるノイズ減衰率を考慮して下記式にて設定すれば良い。
【０１２０】
【数３４】
ｈは図８（ｂ）に示したフィルタ行列であり、（↑ａ） ・はａ倍のアップサンプリングを表す。 また、＊は畳み込みを表す。
【０１２１】
以上、レベルＭを除くすべての画素に対して上述の処理を実行し、ε値に基づいて空間方向の重み付け加算を行う。
【０１２２】
次に、空間方向の重み付け加算した現フレームのラプラシアン画像Ｌｓｔ，ｉとｓ２０６で求めたフィードバック係数ａｍを用いて、時空間方向に重み付け加算した現フレームのラプラシアン画像Ｌｏｔ，ｉを生成する（ｓ５０７）。具体的には、各レベルｉのすべての画素に対して下記式を実行する。
【０１２３】
【数３５】
ここで、寄与率ｂは実施例１と同様に、上記式によるノイズ減衰率Ｒｏに基づいて設定すれば良い。
【０１２４】
次に、再構成部４０２において、時空間方向に重み付け加算した現フレームのラプラシアン画像Ｌｏｔ，ｉとレベルＭのガウシアン画像Ｇｍ，ｉを用いて時空間方向に重み付け加算した現フレームＦｏｔを生成する（ｓ５０８）。本実施の形態ではＢｕｒｔとＡｄｅｌｓｏｎによるピラミッドアルゴリズムを用いて階層画像に分解したので、その逆変換により再構成を行う。
【０１２５】
ここで、ピラミッドアルゴリズムの再構成について図８の（ｂ）及び（ｃ）を参照して説明する。図８（ｃ）はレベルｉのラプラシアン画像Ｌｉとレベルｉ＋１のガウシアン画像Ｇｉ＋１に対する１レベルの再構成を表しており、まずガウシアン画像Ｇｉ＋１は２倍のアップサンプリング８５によって０値が挿入された後、図８（ｂ）に示すフィルタ係数に４倍のゲインを掛けたフィルタによってフィルタリング（８６）され、解像度が縦横ともに２倍の画像を生成する。そして、ラプラシアン画像Ｌｉとこの画像を加算器８７で加算することによって、再構成されたレベルｉのガウシアン画像Ｇｉを生成する。以上の動作により、レベルｉのラプラシアン画像Ｌｉとレベルｉ＋１のガウシアン画像Ｇｉ＋１からレベルｉのガウシアン画像Ｇｉを再構成することができる。さらにマルチレベルの再構成を行う場合は、出力されたガウシアン画像Ｇｉと１つ下のレベルのラプラシアン画像Ｌｉ−１を入力とし同一の処理を繰り返せば良い。
【０１２６】
なお、本実施の形態では、時空間方向に重み付け加算した現フレームのラプラシアン画像｛Ｌｏｔ，ｉ｜ｉ＝０，１，…Ｍ−１｝と未処理のレベルＭのガウシアン画像Ｇｍ，ｉを入力として上述の処理をＭ回繰り返すことで、ノイズ低減がなされた現フレームＦｏｔを生成する。
【０１２７】
次に、全フレームに対してノイズ低減が実行されたか否かを判定し（ｓ２０９）、全フレームに対してノイズ低減が実行されていない場合は、ｔをインクリメントすることで次のフレームを現フレームとして設定し（ｓ２１０）、ｓ５０１〜ｓ５０８の各ステップを実行する。全フレームに対してノイズ低減が実行された場合は、処理を終了する。
【０１２８】
以上、実施例３では画像を複数の周波数帯域に分解し、各々の周波数帯域に対する有意確率を求めることで、周波数毎に異なる変化を伴うような非剛体運動においても精度良く動きを判定でき、実施例１よりも最適な時空間方向の重み付け加算を行うことができる。
【０１２９】
［実施例４］
本発明は、Ｘ線動画撮影装置４００において、画像処理部１１２の動作を実施例３とは異なる図６のフローチャートに従った動作とする。なお、実施例４は実施例３に対して実施例２と同じ構成を加えたものである。従って、図６に示すフローチャートにおいて、図３または図５と同様に処理実行するステップは同じ符号を付し、ここでは、上述した実施例２または実施例３と異なる構成についてのみ具体的に説明する。
【０１３０】
まず、処理を開始する現フレームＦｔとしてｔ＝２を設定するとともに、複数のフィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝を設定する（ｓ３０１）。
【０１３１】
次に、分解部４０１において、現フレームＦｔを複数の異なる解像度を持つ階層画像に分解する（ｓ５０１）。
【０１３２】
次に、設定した複数のフィードバック係数の１つを選択し（ｓ３０２）、ｓ５０２〜ｓ５０５の各ステップを実施例３と同様に実行することで有意確率を求める。
【０１３３】
以上の処理をすべてのフィードバック係数に対して実行するまで処理を繰り返し（ｓ３０３）、フィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝に対する有意確率｛Ｐｉ｜ｉ＝１，２，…ｎ｝を求める。
【０１３４】
次に、加算部１１７において、複数のフィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝の中から実施例２と同様に最適なフィードバック係数ａｍを設定する（ｓ３０４）。
【０１３５】
以上、すべての画素に対して上述の処理を実行し、複数の有意確率｛Ｐｉ｜ｉ＝１，２，…ｎ｝に基づいたフィードバック係数ａｍを設定する。
【０１３６】
次に、現フレームのラプラシアン画像Ｌｔ，ｉのすべての画素に対して空間方向の重み付け加算を行い（ｓ５０６）、空間方向に重み付け加算したラプラシアン画像Ｌｓｔ，ｉとｓ３０４のステップで求めたフィードバック係数ａｍを用いて、時空間方向に重み付け加算した現フレームのラプラシアン画像Ｌｏｔ，ｉを生成する（ｓ６０１）。具体的には、各レベルｉのすべての画素に対して下記式を実行する。
【０１３７】
【数３６】
ここで、寄与率ｂは実施例１と同様に、上記式によるノイズ減衰率Ｒｏに基づいて設定すれば良い。
【０１３８】
なお、上記式ではａｍ＝０の場合は時間方向の重み付け加算が行われず、空間方向の重み付け加算のみが行われる。そのため、検定部１１６において最適なフィードバック係数としてａｍ＝０が設定された場合（すべてのフィードバック係数を棄却しａｍ＝０を設定した場合は除く）、１つ前のフレームのラプラシアン画像Ｈｔ−１，ｉとの間に被検者の動きによる信号変化がないにも関わらず、１つ前のフレームのラプラシアン画像Ｈｔ−１，ｉの情報が利用されないことになる。そこで、この場合は、例外的に下記式にて処理を行っても良い。
【０１３９】
【数３７】
ここで、ａはユーザが操作部１１０を介して設定したフィードバック係数である。また、Ｈｏｔ−１，ｉは時空間方向に重み付け加算された前フレームのラプラシアン画像である。
【０１４０】
上記式は、空間方向の重み付け加算のみを実行する代わりに、時空間方向に重み付け加算された前フレームのラプラシアン画像Ｈｏｔ−１，ｉとを重み付け加算するものである。この場合、時間方向の情報も利用した処理となるため、エッジの境界等の構造のぼけを抑えることができる。
【０１４１】
次に、再構成部４０２において、時空間方向に重み付け加算した現フレームのラプラシアン画像Ｌｏｔ，ｉとレベルＭのガウシアン画像Ｇｍ，ｉを用いて時空間方向に重み付け加算した現フレームＦｏｔを生成する（ｓ５０８）。
【０１４２】
次に、全フレームに対してノイズ低減が実行されたか否かを判定し（ｓ２０９）、全フレームに対してノイズ低減が実行されていない場合は、ｔをインクリメントすることで次のフレームを現フレームとして設定し（ｓ２１０）、ｓ３０２〜ｓ５０８の各ステップを実行する。全フレームに対してノイズ低減が実行された場合は、処理を終了する。
【０１４３】
以上、実施例４では複数の有意確率に基づき確率的に最もリスクの少ないフィードバック係数、または確率的にリスクが少ないものの内、最もノイズ抑制効果の高いフィードバック係数を選択するものであり、実施例３よりも最適なフィードバック係数を設定できる。
【０１４４】
以上、本発明の好ましい実施例について説明したが、本発明はこれらの実施例に限定されないことはいうまでもなく、その要旨の範囲内で種々の変形及び変更が可能である。
【０１４５】
また、本発明は、前述した実施例の機能を実現するソフトウェアのプログラム（実施例では図に示すフローチャートに対応したプログラム）を、システムあるいは装置に直接あるいは遠隔から供給し、そのシステムあるいは装置のコンピュータが該供給されたプログラムコードを読み出して実行することによっても達成される場合を含む。
【０１４６】
従って、本発明の機能処理をコンピュータで実現するために、該コンピュータにインストールされるプログラムコード自体も本発明を実現するものである。つまり、本発明は、本発明の機能処理を実現するためのコンピュータプログラム自体も含まれる。
【０１４７】
プログラムを供給するための記録媒体としては、例えば、ハードディスク、光ディスク、光磁気ディスク、ＭＯ、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ、磁気テープ、不揮発性のメモリカード、ＲＯＭ、ＤＶＤ（ＤＶＤ−ＲＯＭ，ＤＶＤ−Ｒ）などがある。
【０１４８】
その他、プログラムの供給方法としては、クライアントコンピュータのブラウザを用いてインターネットのホームページに接続し、該ホームページから本発明のコンピュータプログラムそのもの、もしくは圧縮され自動インストール機能を含むファイルをハードディスク等の記録媒体にダウンロードすることによっても供給できる。また、本発明のプログラムを構成するプログラムコードを複数のファイルに分割し、それぞれのファイルを異なるホームページからダウンロードすることによっても実現可能である。つまり、本発明の機能処理をコンピュータで実現するためのプログラムファイルを複数のユーザに対してダウンロードさせるＷＷＷサーバも、本発明に含まれるものである。
【０１４９】
また、本発明のプログラムを暗号化してＣＤ−ＲＯＭ等の記憶媒体に格納してユーザに配布し、所定の条件をクリアしたユーザに対し、インターネットを介してホームページから暗号化を解く鍵情報をダウンロードさせ、その鍵情報を使用することにより暗号化されたプログラムを実行してコンピュータにインストールさせて実現することも可能である。
【０１５０】
また、コンピュータが、読み出したプログラムを実行することによって、前述した実施例の機能が実現される他、そのプログラムの指示に基づき、コンピュータ上で稼動しているＯＳなどが、実際の処理の一部または全部を行ない、その処理によっても前述した実施例の機能が実現され得る。
【符号の説明】
【０１５１】
１１２ 画像処理部
１１３ フィルタ部
１１４ ノイズ取得部
１１５ 統計量取得部
１１６ 検定部
１１７ 加算部
１８０ 空間フィルタ部
４０１ 分解部
４０２ 再構成部
【技術分野】
【０００１】
本発明は、画像のノイズを低減する技術に関し、特に時間方向に画像を加算することでノイズを低減する技術に関する。
【背景技術】
【０００２】
従来、Ｘ線動画撮影の分野では、入射されたＸ線を蛍光体とイメージインテンシファイアにより可視光に変換し、ＣＣＤを用いたＴＶカメラで可視光像を撮像する２次元のＸ線動画撮影装置が提案されている。また、近年では撮像系をフラットパネルディテクタ（ＦＰＤ）で置き換えたＸ線動画撮影装置も広く普及している。
【０００３】
ところで、上述のＸ線動画撮影装置で得られる動画像（連続的に収集された複数の静止画像のセットのことであり、以下、個々の静止画像をフレームと呼ぶ）には、Ｘ線量子の揺らぎよる量子ノイズや検出器や回路等から発生するシステムノイズなど、様々なノイズが一連の撮影過程において重畳する。特に、医療用のＸ線動画撮影においては被検者への被曝の観点から少ないＸ線量で撮影を行うため、このようなノイズの影響は無視できず、画質を向上させるためにはノイズを低減することが重要となってくる。
【０００４】
そこで、画像処理の分野においては、従来から時間方向もしくは空間方向に画素の値を重み付け加算（所謂フィルタリング）することによって、上述のようなノイズを低減することが行われている。
【０００５】
具体的には、時間方向の重み付け加算は、対象となるフレームと過去のフレームの同一画素の値を平均化することで、時間的に略無相関なノイズの揺らぎを減少させる。また、空間方向の重み付け加算は、対象となる画素とその周辺の画素の値を平均化することで空間的に略無相関なノイズの揺らぎを減少させる。
【０００６】
なお、時間方向の重み付け加算は、フレーム間において被検者の動きによる時間的な信号変化がない場合には理想的な方法である。しかしながら、信号変化がある場合においては、過去のフレームの像が対象のフレーム上に残像として発生するという課題がある。また、空間方向の重み付け加算は、過去のフレームを参照しないため、残像は発生しないが、エッジの境界等の空間的に類似しない構造をノイズと共に不鮮明にするため、効果はあるものの好ましい方法とは言えない。
【０００７】
そこで、時間方向と空間方向を上手く組み合わせた方法も提案されており、例えば特許文献１の方法がある。この方法は、対象となるフレームと過去のフレームの局所平均の差に基づいて動きを判定し、動きが大きい場合は時間方向から空間方向の重み付け加算に切り替えるものである。これにより、動きがない場合は理想的なノイズ低減を行うことができ、動きが有った場合においても、残像を発生させないノイズ低減を行うことができる。
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特開平６−４７０３６号公報
【発明の概要】
【発明が解決しようとする課題】
【０００９】
ところで、上述のように時間方向と空間方法の重み付け加算を切り替える場合、ノイズが重畳した信号から動きをいかに精度良く判定するかが重要である。しかしながら、特許文献１の方法では、重畳したノイズの影響が十分に考慮されていない。具体的には、特許文献１の方法は、対象となるフレームと過去のフレームの局所平均の差に基づいて動きを判定しているが、その差が、被検者の動きによる時間的な信号変化であるのか、ノイズの時間的な揺らぎによる信号変化であるかを厳密に判定するものではないため、動きを精度良く判定できないという課題がある。
【００１０】
そこで、本発明の例示的な目的は、上記の課題を解決するためになされたもので、ノイズの影響を考慮することで動画像から動きを精度良く判定し、その判定結果に応じた好適なノイズ低減を行うことにある。
【課題を解決するための手段】
【００１１】
上記目的を達成するために、本発明による画像処理装置は複数フレームの動画像を時系列に重み付け加算する画像処理装置であって、
前記複数フレームの動画像の中で、現フレームより前のフレーム画像をリカ―シブフィルタ処理するフィルタ手段と、
前記現フレームの画像と前記フィルタ手段で得られた画像との差分画像におけるノイズの統計量に基づいて現フレームの重み付け係数を得る係数取得手段と、
前記係数取得手段で得られた前記複数フレームに対応する重み付け係数に基づき前記複数フレームの動画像を加算する加算手段と、
を有することを徴とする。
【００１２】
上記目的を達成するために、本発明による画像処理方法は、複数フレームの動画像を時系列に重み付け加算する画像処理方法であって、
前記複数フレームの動画像の中で、現フレームより前のフレーム画像をリカ―シブフィルタ処理するフィルタ工程と、
前記現フレームの画像と前記フィルタ工程で得られた画像との差分画像におけるノイズの統計量に基づいて現フレームの重み付け係数を得る係数取得工程と、
前記係数取得手段で得られた前記複数フレームに対応する重み付け係数に基づき前記複数フレームの動画像を加算する加算工程と、
を有することを徴とする。
【００１３】
また、本発明は、上述した画像処理装置、および画像処理方法をコンピュータに実行させるためのプログラム、並びにプログラムを記録するコンピュータ読み取り可能な記録媒体を含む。
【発明の効果】
【００１４】
ノイズの特性に応じた好適なノイズ低減を行うことができる。
【図面の簡単な説明】
【００１５】
【図１】実施例１および２によるＸ線動画撮影装置全体の構成図である。
【図２】実施例１による画像処理部１１２の処理手順を示すフローチャートである。
【図３】実施例２による画像処理部１１２の処理手順を示すフローチャートである。
【図４】実施例３および４によるＸ線動画撮影装置全体の構成図である。
【図５】実施例３による画像処理部１１２の処理手順を示すフローチャートである。
【図６】実施例４による画像処理部１１２の処理手順を示すフローチャートである。
【図７】ノイズの算出方法を説明する図である。
【図８】ピラミッドアルゴリズムを説明する図である。
【図９】階層画像の周波数特性を説明する図である。
【発明を実施するための形態】
【００１６】
［実施例１］
本発明は、例えば図１に示すようなＸ線動画撮影装置１００に適用される。図１に示すように、Ｘ線動画撮影装置１００は、撮影された動画像をモニタ上に出力する際の効果的な画像処理を行う機能を有するＸ線の動画撮影装置である。Ｘ線動画撮影装置１００は、Ｘ線発生部１０１、Ｘ線検出器１０４、データ収集部１０５、前処理部１０６、ＣＰＵ１０８、メインメモリ１０９、操作部１１０、表示部１１１、画像処理部１１２を備えており、これらはＣＰＵバス１０７を介して互いにデータ授受が可能に接続されている。
【００１７】
上述のような画撮影装置１００において、まず、メインメモリ１０９は、ＣＰＵ１０８での処理に必要な各種のデータなどを記憶すると共に、ＣＰＵ１０８のワーキング・メモリとして機能する。ＣＰＵ１０８は、メインメモリ１０９を用いて、操作部１１０からの操作にしたがった装置全体の動作制御等を行う。これによりＸ線動画撮影装置１００は、以下のように動作する。
【００１８】
まず、操作部１１０を介してユーザから撮影指示が入力されると、この撮影指示はＣＰＵ１０８によりデータ収集部１０５に伝えられる。ＣＰＵ１０８は、撮影指示を受けると、Ｘ線発生部１０１及びＸ線検出器１０４を制御してＸ線動画撮影を実行させる。
【００１９】
Ｘ線動画撮影では、まずＸ線発生部１０１が、被検体１０３に対してＸ線ビーム１０２を照射する。Ｘ線発生部１０１から照射されたＸ線ビーム１０２は、被検体１０３を減衰しながら透過して、Ｘ線検出器１０４に到達する。そしてＸ線検出器１０４によりＸ線画像信号が連続的に出力される。なお、本実施の形態では被検体１０３を人体とする。すなわち、Ｘ線検出器１０４から出力されるＸ線動画像は人体画像となる。
【００２０】
データ収集部１０５は、Ｘ線検出器１０４から連続的に出力されたフレーム単位のＸ線画像信号を所定のデジタル信号に変換してＸ線画像データとして前処理部１０６に供給する。前処理部１０６は、データ収集部１０５からの信号（Ｘ線画像データ）に対して、オフセット補正やゲイン補正等の前処理を行う。この前処理部１０６で前処理が行われたＸ線画像データは、ＣＰＵ１０８の制御により、ＣＰＵバス１０７を介して、メインメモリ１０９、画像処理部１１２に順次転送される。
【００２１】
画像処理部１１２は、複数フレームの動画像を時系列に重み付け加算することで動画像のノイズを低減するものである。複数のフレームで構成される動画像を時間方向にリカ―シブフィルタ処理するフィルタ部１１３と、現フレームの画像データからノイズを取得するノイズ取得部１１４と有する。また、ノイズ取得部１１４で取得したノイズの統計量を取得する統計量取得部１１５と、統計量取得部１１５で取得したノイズの統計量から予め取得されているノイズの分布に従うかを示す評価値を得る検定部１１６と、
この評価値で定められる重み係数を、複数のフレームのそれぞれに対応して取得する図示しない係数取得部１１６１と、を有する。そして、現フレームよりも以前に処理した画像データを図示しないメモリ１１６２から取得し、現フレームの画像データを係数取得部１１６１で取得した重み係数で加算しメモリ１１６２に格納する加算部１１７を備えている。また、これらの各構成部は、ＣＰＵバスに接続されている。
【００２２】
以上のような構成を備えたＸ線動画撮影装置１００において、画像処理部１１２の動作について、図２に示すフローチャートを用いて具体的に説明する。
【００２３】
上述のごとく前処理部１０６によって得られた複数のフレームは、ＣＰＵバス１０７を介して画像処理部１１２に転送されるとともに、処理を開始する現フレームＦｔ（ｔは撮影開始から何枚目のフレームであるかを表す変数であり１以上かつ撮影枚数以下の整数である）としてｔ＝２を設定する。
【００２４】
ノイズ取得部１１４は、フィードバック係数ａを設定する（ｓ２０１）。ここで、フィードバック係数ａは０＜ａ＜１の範囲で任意に設定すれば良く、本実施の形態ではユーザが操作部１１０を介して設定するものとする。
【００２５】
次に、フィルタ部１１３において、設定された現フレームＦｔよりも前のフレームとフィードバック係数ａを用いて、すべての画素に対して時間方向の重み付け加算を行う（ｓ２０２）。具体的には、下記式のように前フレームの処理結果を順次フィードバック加算する、いわゆるリカーシブフィルタ（再帰型）によって時間方向の重み付け加算を行う。
【００２６】
【数１】
ここで、Ｆｒｔ−１が現フレームＦｔよりも前のフレーム画像、すなわち１〜ｔ−１枚目のフレームを時間方向に重み付け加算したデータとなる。なお、本実施の形態では上記式のような再帰型フィルタを用いたがこれに限定されるものではなく、異なる再帰型フィルタを用いてもよい、また非再帰型フィルタを用いても良い。
【００２７】
次に、ノイズ取得部１１４において、現フレームＦｔから時間方向の重み付け加算がなされた前フレームＦｒｔ−１を減算することで差分画像Ｄｔを生成する（ｓ２０３）。具体的には、下記式のように、すべての画素に対してそれぞれ単純な減算を行う。
【００２８】
【数２】
【００２９】
次に、統計量取得部１１５において、差分画像Ｄｔに重畳されるノイズの分散σＤ２を差分画像のノイズの統計量として取得する（ｓ２０４）。
【００３０】
ここで、現フレームＦｔに重畳されるノイズの分散をσＦ２とし、時間方向の重み付け加算がなされた前フレームＦｒｔ−１に重畳されるノイズの分散をσＦｒ２とし、ノイズが時間的に無相関であると仮定すれば分散の加法性により差分画像Ｄｔに重畳されるノイズの分散σＤ２は下記式で表される。
【００３１】
【数３】
【００３２】
また、ノイズの分布（例えば分散の確率分布）が時間的に変化しないと仮定すれば、上記式の右辺第二項は、現フレームＦｔに重畳されるノイズの分散σＦ２および、時間方向の重み付け加算によるノイズ減衰率を考慮して下記式のように変形可能である。
【００３３】
【数４】
ここで、ａはユーザが操作部１１０を介して設定したフィードバック係数である。
【００３４】
なお、現フレームＦｔに重畳されるノイズの分散σＦ２は、検出器への入射Ｘ線量（入射Ｘ線量子数）に比例する量子ノイズを含むため、実際には入射Ｘ線量に依存した値となる。具体的には、検出器への入射Ｘ線量に比例する量子ノイズをσＱ２とし、入射Ｘ線量に依存しないシステムノイズをσＳ２とし、両者が無相関であると仮定すれば、現フレームＦｔに重畳されるノイズの統計量である分散σＦ２は、下記式で表される。
【００３５】
【数５】
ただし、Ｖは現フレームＦｔの画素値であり、ｇは入射Ｘ線量に比例した画素値Ｖを入射Ｘ線量相当の値に変換するためのゲイン係数である。
【００３６】
ここで、上記式を考慮すれば、差分画像Ｄｔに重畳されるノイズの分散σＤ２は下記式となる。
【００３７】
【数６】
従って、実際には現フレームＦｔの画素値Ｖに応じて分散σＤ２を推定する必要がある。そこで、本実施の形態では、対象画素（ｘ，ｙ）の画素値Ｆｔ（ｘ，ｙ）または、対象画素（ｘ，ｙ）を中心とした任意のｎ×ｎ画素の平均画素値をＶ として分散σＤ２（Ｖ）をすべての画素で推定する。
【００３８】
なお、上記式の量子ノイズσＱ２、システムノイズσＳ２、ゲイン係数ｇのパラメータは、撮影過程で一意に決まる値であり、予め取得しておくことが可能である。従って、予め取得したパラメータをメインメモリ１０９に保存しておき、その結果をロードして用いればよい。
【００３９】
ここで、上記パラメータの算出方法については、特に限定するものではないが、例えば、製品検査時に被検体のない状態で検出器への入射Ｘ線量を変えたｎ回の撮影を行い、図７のように各撮影画像の平均画素値Ｖと分散σ２の関係を求めた離散データ｛（Ｖｉ，σｉ２）｜ｉ＝１，２，．．．，ｎ｝を得、この離散データを良く近似するパラメータを算出する。具体的には、下記式で表される残差の２乗和Ｊが、最小となるパラメータを最小二乗法によって求め、このパラメータを予めメインメモリ１０９に保存しておけば良い。
【００４０】
【数７】
【００４１】
次に、検定部１１６において、統計量取得部１１５で求めた差分画像Ｄｔのノイズの統計量が予めえられたノイズの統計分（例えば分散値の確率分布）に従うか否かを表す評価値としての有意確率を求める（ｓ２０５）。ここで、有意確率とは、仮説検定において、ある仮説（帰無仮説）のもとで得られた検定統計量が実現する確率であり、この有意確率が小さい場合は、確率的に仮説が成り立つ可能性が低いことを表すものである。そこで、本実施の形態では「差分画像Ｄｔのノイズの統計量が予め得られたノイズの確率分布Ｎ（μ，σ２）に従う」という帰無仮説に対して、その有意確率Ｐを求める。
【００４２】
なお、「差分画像Ｄｔのノイズの統計量がノイズの確率分布Ｎ（μ，σ２）に従う」ということは、差分画像Ｄｔにはノイズ以外の信号成分が含まれないことを意味し、言い換えれば「差分したフレーム間には被検者の動きによる信号変化がない」ということである。よって、求めた有意確率Ｐは、差分したフレーム間に被検体の動きによる信号変化がないかを確率的に表したものとなる。なお、ノイズの統計量としてノイズの平均値、標準偏差、分散の少なくともいずれかを用いることがえきる。
【００４３】
なお、上述のごとく差分画像Ｄｔに重畳するノイズは画素毎に異なる確率分布（分散）を持つ。そこで、本実施の形態では対象画素（ｘ，ｙ）を中心とした任意のｎ×ｎ画素の範囲を標本集合Ωとし、「標本集合Ωがノイズの確率分布Ｎ（μ，σ２）に従う」という帰無仮説に対して、すべての画素で有意確率Ｐを算出する。ここで、ノイズの母平均μは０であり、母分散σ２は統計量取得部１１５で求めた対象画素（ｘ，ｙ）に対する分散σＤ２（Ｖ）を推定値として用いる。また、ｎは任意のサイズを設定すれば良く、例えば本実施の形態ではｎ＝３とする。
【００４４】
なお、有意確率Ｐの算出方法に関しては、特に限定するものではないが、本実施の形態では、母平均μおよび母分散σ２に対する検定統計量から求めた有意確率を統合することで、上記の帰無仮説に対する有意確率Ｐを算出する。
【００４５】
まず、対象画素（ｘ，ｙ）を中心とした任意のｎ×ｎ画素の範囲の標本集合Ωに対し、母平均μに対する検定統計量Ｔを下記式にて求める。
【００４６】
【数８】
ただし、Ｅは標本平均であり、Ｓ２は不偏分散を表す。また、母平均μは０である。
【００４７】
ここで、上記式にて求めた検定統計量Ｔは、母平均がμの場合に自由度２ｎ−１のＴ分布に従う。よって、「母平均がμである」という帰無仮説おける両側検定（「母平均がμでない」を対立仮説とする検定）に対する有意確率ＰＴは、下記式にて求めることができる。
【００４８】
【数９】
ただし、Γ（・）はガンマ関数を表す。
【００４９】
なお、上記式で求めた有意確率ＰＴは、「標本集合Ωの母平均が０である」という帰無仮説が成り立つ可能性を確率的に求めたものであり、例えばＰＴ＝０．０１であれば、帰無仮説は１％の確率で正しいことになる。言い換えれば、帰無仮説は９９％の確率で誤りであり、「標本集合Ωの母平均が０でない」という対立仮説が成り立つ可能性が非常に高いこと表している。
【００５０】
次に、対象画素（ｘ，ｙ）を中心とした任意のｎ×ｎ画素の範囲の標本集合Ωに対し、母分散σ２に対する検定統計量χ２を下記式にて求める。
【００５１】
【数１０】
ただし、Ｓ２は不偏分散を表す。また、母分散σ２は統計量取得部１１５で求めた対象画素（ｘ，ｙ）に対する分散σＤ２（Ｖ）である。
【００５２】
ここで、上記式にて求めた検定統計量χ２は、母分散がσ２の場合に自由度２ｎ−１のχ２分布に従う。よって、「母分散がσ２である」という帰無仮説における両側検定（「母分散がσ２でない」を対立仮説とする検定）に対する有意確率Ｐχは、下記式にて求めることができる。
【００５３】
【数１１】
ただし、Γ（・）はガンマ関数を表す。
【００５４】
なお、上記式で求めた有意確率Ｐχは、「標本集合Ωの母分散がσＤ２（Ｖ）である」という帰無仮説が成り立つ可能性を確率的に求めたものであり、例えばＰＴ＝０．０１であれば、帰無仮説は１％の確率で正しいことになる。言い換えれば、帰無仮説は９９％の確率で誤りであり、「標本集合Ωの母分散がσＤ２（Ｖ）でない」という対立仮説が成り立つ可能性が非常に高いことを表している。
【００５５】
以上により求めた有意確率ＰＴおよびＰχは、それぞれ標本集合Ωの母平均および母分散に対する有意確率を独立に求めたものであり、本実施の形態ではその両者を満たす、すなわち「標本集合Ωがノイズの確率分布Ｎ（μ，σ２）に従う」という帰無仮説に対して、その有意確率Ｐを求めたい。そこで、求めた２つの有意確率ＰＴおよびＰχを統合し、上記帰無仮説に対する有意確率Ｐを算出する。なお、有意確率の統合方法に関しては特に限定するものではないが、本実施の形態では、Ｆｉｓｈｅｒの方法を用いることとする。具体的には下記式にて有意確率を統合した検定統計量Ｆを求める。
【００５６】
【数１２】
ここで、上記式で統合した検定統計量Ｆは、独立に求めた２つの有意確率ＰＴおよびＰχを統合したものであり、これは自由度４のχ２分布に従う。また、求めた検定統計量Ｆは、２つの有意確率ＰＴおよびＰχが小さくなるほど大きな値をとる。よって、下記式のように、片側検定（上側）に対する有意確率をＰとして算出する。
【００５７】
【数１３】
以上により、対象画素（ｘ，ｙ）に対する有意確率Ｐを算出することができる。ここで、同様の処理をすべての画素に対して実行する。なお、上記式で求めた有意確率Ｐは、「標本集合Ωがノイズの確率分布Ｎ｛０， σＤ２（Ｖ）｝に従う」、すなわち「標本集合Ωには被検者の動きによる信号変化がない」という帰無仮説が成り立つ可能性を確率的に求めたものであり、例えばＰ＝０．０１であれば、帰無仮説は１％の確率で正しいことになる。言い換えれば、帰無仮説は９９％の確率で誤りであり、「標本集合Ωには被検者の動きによる信号変化が有る」という対立仮説が成り立つ可能性が非常に高いことを表している。
【００５８】
次に、加算部１１７において、ｓ２０６〜ｓ２０８のステップを実行することで、有意確率に基づいた時空間方向の重み付け加算を行う。まず、画素毎に求めた有意確率Ｐに従って時間方向の重みづけ係数としてフィードバック係数ａｍを設定する（ｓ２０６）。ここで、上述のごとく有意確率Ｐは「差分したフレーム間には被検者の動きによる信号変化がない」という帰無仮説が成り立つ可能性を確率的に求めたものであり、有意確率Ｐが小さいほど「差分したフレーム間には被検者の動きによる信号変化が有る」という対立仮説が成り立つ可能性が高くなる。よって、このような場合に前フレームの処理結果に対するフィードバック係数を大きく設定すると残像が発生する可能性がある。そこで、下記式のように有意確率Ｐが小さいほどフィードバック係数ａｍが小さくなるように設定する。
【００５９】
【数１４】
ここで、ａは、ユーザが操作部１１０を介して設定したフィードバック係数であり、この値を基準として実際に用いるフィードバック係数ａｍを設定する。
【００６０】
なお、フィードバック係数ａｍの設定方法はこれに限定されるものではなく、例えば下記式のように有意確率Ｐが設定した有意水準αよりも小さい場合、すなわち仮説検定で言う所の対立仮説が採択される場合は０（時間方向の重み付け加算を実行しない）に設定し、有意確率Ｐが設定した有意水準α以上の場合、すなわち仮説検定で言う所の帰無仮説を棄却できない場合はａに設定しても良い。
【００６１】
【数１５】
ここで、有意水準αは仮説検定で一般的に使われる値を用いれば良く、本実施の形態では例えばα＝０．０１とする。
【００６２】
以上、すべての画素に対して上述の処理を実行し、有意確率Ｐに基づいたフィードバック係数ａｍを設定する。なお、このように設定したフィードバック係数ａｍは、各画素の近傍における動きに基づいたものであり、局所的な動きに応じて最適なフィードバック係数が得られるものである。
【００６３】
次に、図示しない空間フィルタ部１８０は、現フレームＦｔのすべての画素に対して空間方向の重み付け加算を空間フィルタ処理として行う（ｓ２０７）。なお、空間方向の重み付け加算は、ＦＩＲ等の単純な空間フィルタで行うこともできるが、エッジの境界等の空間的に類似しない構造をノイズと共に不鮮明にする。そこで、本実施の形態では下記式で表されるεフィルタを空間フィルタとして用いることで、構造のぼけを抑えた重み付け加算を行う。
【００６４】
【数１６】
ただし、ｈはサイズがＮ×Ｎのローパスフィルタであり、本実施の形態では例えば９×９の平均化フィルタを用いる。
【００６５】
なお、上述の処理は、対象画素（ｘ，ｙ）と近傍画素の画素値の差を類似度の指標値とし、この指標値がε値よりも大きいほど近傍画素に対するローパスフィルタの重みを小さくするものである。これにより、類似しない近傍画素を除外した空間方向の重み付け加算がなされ、構造のぼけを抑えることができる。
【００６６】
ここで、ε値は対象画素（ｘ，ｙ）と近傍画素が類似しているか否かを判定する閾値であり、本実施の形態では、下記式にてε値を設定する。
【００６７】
【数１７】
ただし、σ２は現フレームＦｔに重畳するノイズの分散であり、統計量取得部１１５で求めた対象画素（ｘ，ｙ）に対する分散σＦ２（Ｖ）を用いて、各画素でε値を設定する。
【００６８】
なお、上記式の意味するところは、もし対象画素（ｘ，ｙ）と近傍画素がノイズを除けば同じ画素値であると仮定すると、その差はノイズの揺らぎによる信号変化であり、ノイズが正規分布ならば約９９．７％の確率で±３√２σ２の範囲となる。よって、この値をε値として設定することで、もし対象画素（ｘ，ｙ）と近傍画素の画素値の差が想定されるノイズの揺らぎよりも大きい場合、その変化はノイズの揺らぎによる信号変化ではない。すなわち、対象画素（ｘ，ｙ）と近傍画素は類似しない画素であると判定するものである。
【００６９】
以上、すべての画素に対して上述の処理を実行し、ε値に基づいた空間方向の重み付け加算を行う。なお、本実施の形態ではεフィルタを用いたが、これに限定されるものではなくＢｉｌａｔｅｒａｌフィルタ等の類似の非線形フィルタを用いることもできる。
【００７０】
次に、空間方向に重み付け加算した現フレームＦｓｔとｓ２０６のステップで求めたフィードバック係数ａｍを用いて、時空間方向に重み付け加算した現フレームＦｏｔを生成する（ｓ２０８）。具体的には、すべての画素に対して下記式の処理を実行する。
【００７１】
【数１８】
ここで、Ｆｒｔ−１は時間方向の重み付け加算がなされた前フレームであり、加算部１１７で生成したものである。また、ｂは空間方向の重み付け加算の寄与率を表し、フィードバック係数ａｍが小さいほど大きくなるように設定する。すなわち、残像が発生する可能性がある画素に対しては、時間方向の寄与を小さくする代わりに、空間方向の寄与を大きくし空間フィルタ処理の処理効果を上げ、十分なノイズ低減効果を得るものである。
【００７２】
なお、寄与率ｂの設定方法は特に限定するものではないが、本実施の形態では、上記式によるノイズ減衰率Ｒｏ（入力のノイズの分散に対する出力のノイズの分散の比）が、基準値Ｒと略同等となるように寄与率ｂを設定する。ここで、上記式によるノイズ減衰率Ｒｏは、空間方向の重み付け加算（上記式の右辺第二項の波括弧内）におけるノイズ減衰率をＲｓとすれば下記式で表される。
【００７３】
【数１９】
ただし、ａはユーザが操作部１１０を介して設定したフィードバック係数である。
従って、Ｒｏに基準値Ｒを代入し、Ｒｓについて解けば空間方向の重み付け加算が担うノイズ減衰率が求まる。よって、空間方向の重み付け加算によるノイズ減衰率がＲＳとなるような寄与率ｂを設定する。なお、この条件を満たす厳密解を求めるのは困難である。そこで、空間方向に重み付け加算した現フレームＦｓｔに重畳するノイズが無視可能なレベルと仮定し、下記式にて求まる近似値を寄与率ｂとして設定する。
【００７４】
【数２０】
なお、基準値Ｒは任意に設定すればよく、本実施の形態ではユーザが操作部１１０を介して設定したフィードバック係数ａを用いて下記式にて設定する。
【００７５】
【数２１】
ここで、上記式で設定される基準値Ｒは、１〜ｔ枚目のフレームをリカーシブフィルタによって処理した場合のノイズ減衰率そのものである。従って、このように設定した基準値Ｒを用いることで、すべての画素においてフィードバック係数ａを用いたリカーシブフィルタと略同等のノイズ低減効果が得られるものである。
【００７６】
なお、リカーシブフィルタのノイズ減衰率は時間依存性があり、ｔが大きくなるほど小さな値となる。よって、上記式で設定した基準値Ｒを用いるとフレーム間でノイズ低減効果が同等とならない。そこで、下記式のようにｔ→∞の場合のノイズ減衰率を基準値Ｒとして設定しても良い。この場合、すべてのフレームにおいて略同等のノイズ低減効果を得ることができる。
【００７７】
【数２２】
以上、すべての画素に対して上述の処理を実行することで、ノイズ低減がなされた現フレームＦｏｔを生成することができる。
【００７８】
次に、全フレームに対してノイズ低減が実行されたか否かを判定し（ｓ２０９）、全フレームに対してノイズ低減が実行されていない場合は、ｔをインクリメントすることで次のフレームを現フレームとして設定し（ｓ２１０）、ｓ２０２〜ｓ２０８の各ステップを実行する。全フレームに対してノイズ低減が実行された場合は、処理を終了する。
【００７９】
以上、実施例１ではノイズの確率分布に基づいて動きである可能性を確率的に求めることで、精度良く動きを判定することができる。また、動きである可能性が高い画素では時間方向のフィードバック係数を小さくすることで残像を抑え、かつ空間方向の寄与率をノイズ減衰率に基づいて設定することで、すべての画素で略同等のノイズ抑制効果を得ることができる。
【００８０】
［実施例２］
本発明は、Ｘ線動画撮影装置１００において、画像処理部１１２の動作を実施例１とは異なる図３のフローチャートに従った動作とする。なお、図３に示すフローチャートにおいて、図２に示したフローチャートと同様に処理実行するステップは同じ符号を付し、ここでは、上述した実施例１とは異なる構成についてのみ具体的に説明する。
【００８１】
まず、処理を開始する現フレームＦｔとしてｔ＝２を設定するとともに、複数のフィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝を設定する（ｓ３０１）。ここで、フィードバック係数はユーザが操作部１１０を介して設定したフィードバック係数ａを最大値とし、下記式のように０〜ａの範囲をｎ等分したフィードバック係数を設定する。
【００８２】
【数２３】
ここで、ｎは任意の数を設定すれば良く、例えばｎ＝１０とする。
【００８３】
次に、設定した複数のフィードバック係数の１つを選択し（ｓ３０２）、ｓ２０２〜ｓ２０５の各ステップを実施例１と同様に実行することで有意確率を求める。
【００８４】
以上の処理をすべてのフィードバック係数に対して実行するまで処理を繰り返し（ｓ３０３）、フィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝に対する有意確率｛Ｐｉ｜ｉ＝１，２，…ｎ｝を求める。
【００８５】
次に、加算部１１７において、複数のフィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝の中から最適なフィードバック係数ａｍを設定する（ｓ３０４）。ここで、上述した通り有意確率は「差分したフレーム間には被検者の動きによる信号変化がない」という帰無仮説が成り立つ可能性を確率的に求めたものであり、有意確率が小さいほど「差分したフレーム間には被検者の動きによる信号変化が有る」という対立仮説が成り立つ可能性が高くなる。そこで、本実施の形態では有意確率が最も大きくなるフィードバック係数、すなわち対立仮説が成り立つ可能性が最も低いフィードバック係数を最適なフィードバック係数ａｍとして設定する。ただし、最大の有意確率が有意水準αよりも小さい場合、すなわち対立仮説がすべてのフィードバック係数において採択される場合は、すべてのフィードバック係数を棄却し、ａｍ＝０（時間方向の重み付け加算を実行しない）に設定する。
【００８６】
なお、フィードバック係数ａｍの設定方法はこれに限定されるものではなく、例えば、有意確率が有意水準αよりも大きいフィードバック係数の内、最大のものを最適なフィードバック係数ａｍとして設定しても良い。この場合は、設定した有意水準αを満たすものも内、最大のノイズ抑制効果を得られるフィードバック係数が選択されることとなる。
【００８７】
以上、すべての画素に対して上述の処理を実行し、複数の有意確率｛Ｐｉ｜ｉ＝１，２，…ｎ｝に基づいたフィードバック係数ａｍを設定する。
【００８８】
次に、現フレームＦｔのすべての画素に対して空間方向の重み付け加算を行い（ｓ２０７）、空間方向に重み付け加算した現フレームＦｓｔとｓ３０４のステップで求めたフィードバック係数ａｍを用いて、時空間方向に重み付け加算した現フレームＦｏｔを生成する（ｓ３０５）。具体的には、すべての画素に対して下記式を実行する。
【００８９】
【数２４】
ここで、寄与率ｂは実施例１と同様に、上記式によるノイズ減衰率Ｒｏに基づいて設定すれば良い。
【００９０】
なお、上記式ではａｍ＝０の場合は時間方向の重み付け加算が行われず、空間方向の重み付け加算のみが行われる。そのため、検定部１１６において最適なフィードバック係数としてａｍ＝０が設定された場合（すべてのフィードバック係数を棄却しａｍ＝０を設定した場合は除く）、１つ前のフレームＦｔ−１との間に被検者の動きによる信号変化がないにも関わらず、１つ前のフレームＦｔ−１の情報が利用されないことになる。そこで、この場合は、例外的に下記式にて処理を行っても良い。
【００９１】
【数２５】
ここで、ａはユーザが操作部１１０を介して設定したフィードバック係数である。また、Ｆｏｔ−１は時空間方向に重み付け加算された前フレームである。
【００９２】
上記式は、空間方向の重み付け加算のみを実行する代わりに、時空間方向に重み付け加算された前フレームＦｏｔ−１とを重み付け加算するものである。この場合、時間方向の情報も利用した処理となるため、エッジの境界等の構造のぼけを抑えることができる。
【００９３】
次に、全フレームに対してノイズ低減が実行されたか否かを判定し（ｓ２０９）、全フレームに対してノイズ低減が実行されていない場合は、ｔをインクリメントすることで次のフレームを現フレームとして設定し（ｓ２１０）、ｓ３０２〜ｓ３０５の各ステップを実行する。全フレームに対してノイズ低減が実行された場合は、処理を終了する。
【００９４】
以上、実施例２では複数の有意確率に基づき確率的に最もリスクの少ないフィードバック係数、または確率的にリスクが少ないものの内、最もノイズ抑制効果の高いフィードバック係数を選択するものであり、実施例１よりも最適なフィードバック係数を設定できる。
【００９５】
［実施例３］
本発明は、例えば図４に示すようなＸ線動画撮影装置４００に適用される。このＸ線動画撮影装置４００は、Ｘ線動画撮影装置１００に対し、分解部４０１、再構成部４０２を備える構成としている。また、本発明は画像処理部１１２の処理手順を実施例１とは異なる図５に示したフローチャートに従った動作とする。
【００９６】
なお、図４のＸ線動画撮影装置４００において、図１のＸ線動画撮影装置１００と同様に動作する箇所は同じ符号を付し、その詳細は省略する。また、図５に示すフローチャートにおいて、図２に示したフローチャートと同様に処理実行するステップは同じ符号を付し、ここでは、上述した実施例１とは異なる構成についてのみ具体的に説明する。
【００９７】
まず、実施例１と同様に処理を開始する現フレームＦｔとしてｔ＝２を設定するとともに、フィードバック係数ａを設定する（ｓ２０１）。
【００９８】
次に、分解部４０１において、現フレームＦｔを複数の異なる解像度を持つ階層画像に分解する（ｓ５０１）。なお、本実施の形態ではＢｕｒｔとＡｄｅｌｓｏｎによるピラミッドアルゴリズムを用いて階層画像に分解する。
【００９９】
ここで、ピラミッドアルゴリズムの分解について図８の（ａ）及び（ｂ）を参照して説明する。図８（ａ）は、入力画像Ｇｉに対して１レベルの分解を表しており、まず入力画像Ｇｉは図８（ｂ）に示すフィルタでフィルタリングされ（８０）、さらに１／２倍のダウンサンプリング８１によって、解像度が縦横ともに１／２のガウシアン画像Ｇｉ＋１を生成する。また、ガウシアン画像Ｇｉ＋１は、さらに２倍のアップサンプリング８２によって０値が挿入された後、図８（ｂ）に示すフィルタ係数に４倍のゲインを掛けたフィルタによってフィルタリングされ（８３）、解像度が入力画像Ｇｉと同一の画像を生成する。そして、入力画像Ｇｉからこの画像を減算器８４にて減算することでラプラシアン画像Ｌｉを生成する。以上の動作により、入力画像Ｇｉをガウシアン画像Ｇｉ＋１、ラプラシアン画像Ｌｉに分解することができる。さらに、マルチレベルの分解を行う場合は、出力されたガウシアン画像Ｇｉ＋１を入力画像として同一の処理を繰り返せばよい。
【０１００】
なお、本実施の形態では、現フレームＦｔを最初の入力画像とし、上述の処理をＭ回繰り返すことで、解像度の異なる複数のガウシアン画像｛Ｇｔ，ｉ｜ｉ＝０，１，…Ｍ｝およびラプラシアン画像｛Ｌｔ，ｉ｜ｉ＝０，１，…Ｍ−１｝を生成する。ここで、このように生成された画像は図９に示す周波数特性を持つものである。図９のようにガウシアン画像はレベルｉが大きくなるほど、入力画像のより低い周波数のみで構成される画像となる。また、ラプラシアン画像は入力画像の一部の周波数帯域で構成される画像となり、レベルｉが大きくなるほど低い周波数帯域を抽出したものとなる。
【０１０１】
次に、フィルタ部１１３において、設定された現フレームよりも前のガウシアン画像とフィードバック係数ａを用いて、時間方向の重み付け加算をおこなう（ｓ５０２）。具体的には、下記式にてレベルＭを除くすべての画素に対しリカーシブフィルタを用いた重み付け加算を行う。
【０１０２】
【数２６】
【０１０３】
次に、ノイズ取得部１１４において、現フレームのガウシアン画像Ｇｔ，ｉから時間方向の重み付け加算がなされた前フレームのガウシアン画像Ｇｒｔ−１，ｉを減算することで差分画像Ｄｇｔ，ｉを生成する（ｓ５０３）。具体的には、下記式にてレベルＭを除くすべての画素に対してそれぞれ単純な減算を行う。
【０１０４】
【数２７】
【０１０５】
次に、統計量取得部１１５において、差分画像Ｄｇｔ，ｉに重畳されるノイズの分散σｄｉ２を推定する（ｓ５０４）。ここで、現フレームのガウシアン画像Ｇｔ，ｉに重畳されるノイズの分散をσｇｉ２とし、時間方向の重み付け加算によるノイズ減衰率を考慮すれば、差分画像Ｄｇｔ，ｉに重畳されるノイズの分散σｄｉ２は下記式で表される。
【０１０６】
【数２８】
ここで、ａはユーザが操作部１１０を介して設定したフィードバック係数である。
【０１０７】
なお、現フレームのガウシアン画像Ｇｔ，ｉに重畳されるノイズの分散σｇｉ２は、実施例１と同様にして求まる現フレームＦｔに重畳されるノイズの分散σＦ２（Ｖ）および分解過程におけるノイズ減衰率を考慮すれば下記式となる。
【０１０８】
【数２９】
ただし、ｈは図８（ｂ）に示したフィルタ行列であり、（↑ａ） ・はａ倍のアップサンプリングを表す。また、＊は畳み込みを表す。
【０１０９】
ここで、上記式はレベルｉおよび画素値Ｖに依存した値をとる。従って、差分画像Ｄｇｔ，ｉに重畳されるノイズの分散σｄｉ２も同じく、レベルｉおよび画素値Ｖに依存した値をとなる。なお、画素値Ｖは、現フレームＦｔの重み付き平均であるガウシアン画像Ｇｔ，ｉを用いれば良いので、実際の処理としては、レベルｉのガウシアン画像Ｇｔ，ｉの画素値Ｇｔ，ｉ （ｘ，ｙ）に基づいてレベルＭを除くすべての画素に対する推定値を算出すれば良い。
【０１１０】
次に、検定部１１６において、差分画像Ｄｇｔ，ｉが統計量取得部１１５で求めたノイズの統計量に基づいて確率分布に従うか否かを表す有意確率を求める（ｓ５０５）。ここで、上述のごとく差分画像Ｄｇｔ，ｉに重畳するノイズはレベル及び画素毎に異なる確率分布（分散）を持つ。そこで、本実施の形態ではレベルｉの対象画素（ｘ，ｙ）を中心とした任意のｎ×ｎ画素の範囲を標本集合Ωとし、「標本集合Ωがノイズの確率分布Ｎ（μ，σ２）に従う」という帰無仮説に対して、各レベルｉのすべての画素で有意確率Ｐｉを算出する。ここで、ノイズの母平均μは０であり、母分散σ２は統計量取得部１１５で求めたレベルｉの対象画素（ｘ，ｙ）に対する分散σｄｉ２（Ｖ）を推定値として用いる。また、ｎは任意のサイズを設定すれば良く、例えば本実施の形態ではｎ＝３とする。
【０１１１】
なお、標本集合Ωに対する有意確率の算出方法は、母平均μに対する検定統計量Ｔの算出方法を下記式に変更する以外は、実施例１と同様に求めることができる。
【０１１２】
【数３０】
ここで、実施例１と上記式の違いはＣであり、これは空間的に相関のあるノイズを持つガウシアン画像Ｇｔ，ｉから求めた検定統計量Ｔを自由度２ｎ−１のＴ分布に従わせるための補正係数である。具体的には、空間的なノイズの相関を考慮して下記式にて算出する。
【０１１３】
【数３１】
ただし、ｈは図８（ｂ）に示したフィルタ行列であり、（↑ａ） ・はａ倍のアップサンプリングを表す。また、＊は畳み込みを表す。
【０１１４】
なお、ｍは標本サイズと同じサイズの平均化フィルタ行列であり、対象画素（ｘ，ｙ）を中心とした３×３画素の範囲を標本集合Ωとする本実施の形態では、下記式となる。
【０１１５】
【数３２】
以上により求めた母平均μに対する検定統計量Ｔと、実施例１と同様に求めた母分散σ２に対する検定統計量χ２を統合することで、レベルＭを除くすべての画素に対する有意確率Ｐｉを求める。なお、このように算出した有意確率Ｐｉは、画像を複数の周波数帯域に分解し、各々の周波数帯域における動きを確率的に求めたものである。そのため、周波数毎に異なる変化を伴うような非剛体運動においても精度良く動きを判定できるものである。
【０１１６】
次に、加算部１１７において、レベルＭを除くすべての画素に対するフィードバック係数ａｍｉを設定する（ｓ２０６）。ここでは、各レベルｉに対して実施例１と同様にフィードバック係数ａｍｉを設定すれば良い。
【０１１７】
次に、現フレームのラプラシアン画像Ｌｔ，ｉのすべての画素に対して空間方向の重み付け加算を行う（ｓ５０６）。なお、本実施の形態では、ラプラシアン画像Ｌｔ，ｉと同じレベルのガウシアン画像Ｇｔ，ｉからεフィルタと同様の原理にてノイズを抽出し、抽出したノイズを現フレームのラプラシアン画像Ｌｔ，ｉから減算することで、構造のぼけを抑えた重み付け加算を行う。具体的には下記式にて処理を行う。
【０１１８】
【数３３】
ただし、ｈは図８（ｂ）に示したフィルタ行列である。
【０１１９】
上記式において、右辺第二項がε値を基準としてラプラシアン画像Ｌｔ，ｉと同一の周波数帯域からノイズを推定したものとなる。ここで、ε値は、対象画素（ｘ，ｙ）と近傍画素がノイズを除けば同じ画素値であるかを判定する閾値であり、現フレームＦｔに重畳されるノイズの分散σＦ２（Ｖ）および分解過程におけるノイズ減衰率を考慮して下記式にて設定すれば良い。
【０１２０】
【数３４】
ｈは図８（ｂ）に示したフィルタ行列であり、（↑ａ） ・はａ倍のアップサンプリングを表す。 また、＊は畳み込みを表す。
【０１２１】
以上、レベルＭを除くすべての画素に対して上述の処理を実行し、ε値に基づいて空間方向の重み付け加算を行う。
【０１２２】
次に、空間方向の重み付け加算した現フレームのラプラシアン画像Ｌｓｔ，ｉとｓ２０６で求めたフィードバック係数ａｍを用いて、時空間方向に重み付け加算した現フレームのラプラシアン画像Ｌｏｔ，ｉを生成する（ｓ５０７）。具体的には、各レベルｉのすべての画素に対して下記式を実行する。
【０１２３】
【数３５】
ここで、寄与率ｂは実施例１と同様に、上記式によるノイズ減衰率Ｒｏに基づいて設定すれば良い。
【０１２４】
次に、再構成部４０２において、時空間方向に重み付け加算した現フレームのラプラシアン画像Ｌｏｔ，ｉとレベルＭのガウシアン画像Ｇｍ，ｉを用いて時空間方向に重み付け加算した現フレームＦｏｔを生成する（ｓ５０８）。本実施の形態ではＢｕｒｔとＡｄｅｌｓｏｎによるピラミッドアルゴリズムを用いて階層画像に分解したので、その逆変換により再構成を行う。
【０１２５】
ここで、ピラミッドアルゴリズムの再構成について図８の（ｂ）及び（ｃ）を参照して説明する。図８（ｃ）はレベルｉのラプラシアン画像Ｌｉとレベルｉ＋１のガウシアン画像Ｇｉ＋１に対する１レベルの再構成を表しており、まずガウシアン画像Ｇｉ＋１は２倍のアップサンプリング８５によって０値が挿入された後、図８（ｂ）に示すフィルタ係数に４倍のゲインを掛けたフィルタによってフィルタリング（８６）され、解像度が縦横ともに２倍の画像を生成する。そして、ラプラシアン画像Ｌｉとこの画像を加算器８７で加算することによって、再構成されたレベルｉのガウシアン画像Ｇｉを生成する。以上の動作により、レベルｉのラプラシアン画像Ｌｉとレベルｉ＋１のガウシアン画像Ｇｉ＋１からレベルｉのガウシアン画像Ｇｉを再構成することができる。さらにマルチレベルの再構成を行う場合は、出力されたガウシアン画像Ｇｉと１つ下のレベルのラプラシアン画像Ｌｉ−１を入力とし同一の処理を繰り返せば良い。
【０１２６】
なお、本実施の形態では、時空間方向に重み付け加算した現フレームのラプラシアン画像｛Ｌｏｔ，ｉ｜ｉ＝０，１，…Ｍ−１｝と未処理のレベルＭのガウシアン画像Ｇｍ，ｉを入力として上述の処理をＭ回繰り返すことで、ノイズ低減がなされた現フレームＦｏｔを生成する。
【０１２７】
次に、全フレームに対してノイズ低減が実行されたか否かを判定し（ｓ２０９）、全フレームに対してノイズ低減が実行されていない場合は、ｔをインクリメントすることで次のフレームを現フレームとして設定し（ｓ２１０）、ｓ５０１〜ｓ５０８の各ステップを実行する。全フレームに対してノイズ低減が実行された場合は、処理を終了する。
【０１２８】
以上、実施例３では画像を複数の周波数帯域に分解し、各々の周波数帯域に対する有意確率を求めることで、周波数毎に異なる変化を伴うような非剛体運動においても精度良く動きを判定でき、実施例１よりも最適な時空間方向の重み付け加算を行うことができる。
【０１２９】
［実施例４］
本発明は、Ｘ線動画撮影装置４００において、画像処理部１１２の動作を実施例３とは異なる図６のフローチャートに従った動作とする。なお、実施例４は実施例３に対して実施例２と同じ構成を加えたものである。従って、図６に示すフローチャートにおいて、図３または図５と同様に処理実行するステップは同じ符号を付し、ここでは、上述した実施例２または実施例３と異なる構成についてのみ具体的に説明する。
【０１３０】
まず、処理を開始する現フレームＦｔとしてｔ＝２を設定するとともに、複数のフィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝を設定する（ｓ３０１）。
【０１３１】
次に、分解部４０１において、現フレームＦｔを複数の異なる解像度を持つ階層画像に分解する（ｓ５０１）。
【０１３２】
次に、設定した複数のフィードバック係数の１つを選択し（ｓ３０２）、ｓ５０２〜ｓ５０５の各ステップを実施例３と同様に実行することで有意確率を求める。
【０１３３】
以上の処理をすべてのフィードバック係数に対して実行するまで処理を繰り返し（ｓ３０３）、フィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝に対する有意確率｛Ｐｉ｜ｉ＝１，２，…ｎ｝を求める。
【０１３４】
次に、加算部１１７において、複数のフィードバック係数｛ａｉ｜ｉ＝１，２，．．．，ｎ｝の中から実施例２と同様に最適なフィードバック係数ａｍを設定する（ｓ３０４）。
【０１３５】
以上、すべての画素に対して上述の処理を実行し、複数の有意確率｛Ｐｉ｜ｉ＝１，２，…ｎ｝に基づいたフィードバック係数ａｍを設定する。
【０１３６】
次に、現フレームのラプラシアン画像Ｌｔ，ｉのすべての画素に対して空間方向の重み付け加算を行い（ｓ５０６）、空間方向に重み付け加算したラプラシアン画像Ｌｓｔ，ｉとｓ３０４のステップで求めたフィードバック係数ａｍを用いて、時空間方向に重み付け加算した現フレームのラプラシアン画像Ｌｏｔ，ｉを生成する（ｓ６０１）。具体的には、各レベルｉのすべての画素に対して下記式を実行する。
【０１３７】
【数３６】
ここで、寄与率ｂは実施例１と同様に、上記式によるノイズ減衰率Ｒｏに基づいて設定すれば良い。
【０１３８】
なお、上記式ではａｍ＝０の場合は時間方向の重み付け加算が行われず、空間方向の重み付け加算のみが行われる。そのため、検定部１１６において最適なフィードバック係数としてａｍ＝０が設定された場合（すべてのフィードバック係数を棄却しａｍ＝０を設定した場合は除く）、１つ前のフレームのラプラシアン画像Ｈｔ−１，ｉとの間に被検者の動きによる信号変化がないにも関わらず、１つ前のフレームのラプラシアン画像Ｈｔ−１，ｉの情報が利用されないことになる。そこで、この場合は、例外的に下記式にて処理を行っても良い。
【０１３９】
【数３７】
ここで、ａはユーザが操作部１１０を介して設定したフィードバック係数である。また、Ｈｏｔ−１，ｉは時空間方向に重み付け加算された前フレームのラプラシアン画像である。
【０１４０】
上記式は、空間方向の重み付け加算のみを実行する代わりに、時空間方向に重み付け加算された前フレームのラプラシアン画像Ｈｏｔ−１，ｉとを重み付け加算するものである。この場合、時間方向の情報も利用した処理となるため、エッジの境界等の構造のぼけを抑えることができる。
【０１４１】
次に、再構成部４０２において、時空間方向に重み付け加算した現フレームのラプラシアン画像Ｌｏｔ，ｉとレベルＭのガウシアン画像Ｇｍ，ｉを用いて時空間方向に重み付け加算した現フレームＦｏｔを生成する（ｓ５０８）。
【０１４２】
次に、全フレームに対してノイズ低減が実行されたか否かを判定し（ｓ２０９）、全フレームに対してノイズ低減が実行されていない場合は、ｔをインクリメントすることで次のフレームを現フレームとして設定し（ｓ２１０）、ｓ３０２〜ｓ５０８の各ステップを実行する。全フレームに対してノイズ低減が実行された場合は、処理を終了する。
【０１４３】
以上、実施例４では複数の有意確率に基づき確率的に最もリスクの少ないフィードバック係数、または確率的にリスクが少ないものの内、最もノイズ抑制効果の高いフィードバック係数を選択するものであり、実施例３よりも最適なフィードバック係数を設定できる。
【０１４４】
以上、本発明の好ましい実施例について説明したが、本発明はこれらの実施例に限定されないことはいうまでもなく、その要旨の範囲内で種々の変形及び変更が可能である。
【０１４５】
また、本発明は、前述した実施例の機能を実現するソフトウェアのプログラム（実施例では図に示すフローチャートに対応したプログラム）を、システムあるいは装置に直接あるいは遠隔から供給し、そのシステムあるいは装置のコンピュータが該供給されたプログラムコードを読み出して実行することによっても達成される場合を含む。
【０１４６】
従って、本発明の機能処理をコンピュータで実現するために、該コンピュータにインストールされるプログラムコード自体も本発明を実現するものである。つまり、本発明は、本発明の機能処理を実現するためのコンピュータプログラム自体も含まれる。
【０１４７】
プログラムを供給するための記録媒体としては、例えば、ハードディスク、光ディスク、光磁気ディスク、ＭＯ、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、ＣＤ−ＲＷ、磁気テープ、不揮発性のメモリカード、ＲＯＭ、ＤＶＤ（ＤＶＤ−ＲＯＭ，ＤＶＤ−Ｒ）などがある。
【０１４８】
その他、プログラムの供給方法としては、クライアントコンピュータのブラウザを用いてインターネットのホームページに接続し、該ホームページから本発明のコンピュータプログラムそのもの、もしくは圧縮され自動インストール機能を含むファイルをハードディスク等の記録媒体にダウンロードすることによっても供給できる。また、本発明のプログラムを構成するプログラムコードを複数のファイルに分割し、それぞれのファイルを異なるホームページからダウンロードすることによっても実現可能である。つまり、本発明の機能処理をコンピュータで実現するためのプログラムファイルを複数のユーザに対してダウンロードさせるＷＷＷサーバも、本発明に含まれるものである。
【０１４９】
また、本発明のプログラムを暗号化してＣＤ−ＲＯＭ等の記憶媒体に格納してユーザに配布し、所定の条件をクリアしたユーザに対し、インターネットを介してホームページから暗号化を解く鍵情報をダウンロードさせ、その鍵情報を使用することにより暗号化されたプログラムを実行してコンピュータにインストールさせて実現することも可能である。
【０１５０】
また、コンピュータが、読み出したプログラムを実行することによって、前述した実施例の機能が実現される他、そのプログラムの指示に基づき、コンピュータ上で稼動しているＯＳなどが、実際の処理の一部または全部を行ない、その処理によっても前述した実施例の機能が実現され得る。
【符号の説明】
【０１５１】
１１２ 画像処理部
１１３ フィルタ部
１１４ ノイズ取得部
１１５ 統計量取得部
１１６ 検定部
１１７ 加算部
１８０ 空間フィルタ部
４０１ 分解部
４０２ 再構成部
【特許請求の範囲】
【請求項１】
複数フレームの動画像を時系列に重み付け加算する画像処理装置であって、
前記複数フレームの動画像の中で、現フレームより前のフレーム画像をリカ―シブフィルタ処理するフィルタ手段と、
前記現フレームの画像と前記フィルタ手段で得られた画像との差分画像におけるノイズの統計量に基づいて現フレームの重み付け係数を得る係数取得手段と、
前記係数取得手段で得られた前記複数フレームに対応する重み付け係数に基づき前記複数フレームの動画像を加算する加算手段と、
を有することを徴とする画像処理装置。
【請求項２】
前記係数取得手段は、予めえられたノイズの分布に従うかを示す評価値に基づいて前記重み係数を得ることを特徴とする請求項１に記載の画像処理装置。
【請求項３】
前記ノイズの統計量が予め得られたノイズの確率分布に従うか否かを表す有意確率を求める検定手段を更に備え、
前記加算手段は、前記有意確率に基づいて前記重みづけ係数を得ることを特徴とする請求項１又は２に記載の画像処理装置。
【請求項４】
複数の異なる周波数帯域の画像データに分解する分解手段を更に備え、
前記画像データは前記分解手段で分解された画像であることを特徴とする請求項１乃至３のいずれか一項に記載の画像処理装置。
【請求項５】
複数の異なる周波数帯域の画像データに再構成する再構成手段を更に備え、
前記再構成手段は、前記加算手段で加算した画像データに基づいて再構成することを特徴とする請求項４に記載の画像処理装置。
【請求項６】
前記差分画像におけるノイズの統計量を取得する統計量取得手段を更に備え、
前記検定手段は、前記ノイズの統計量が予め取得された確率分布に従うか否かを示す有意確率を求めることを特徴とする請求項３乃至５のいずれか一項に記載の画像処理装置。
【請求項７】
前記ノイズの統計量は、ノイズの平均値、標準偏差、分散の少なくとも何れか１つであることを特徴とする請求項１乃至６の何れか一項に記載の画像処理装置。
【請求項８】
前記検定手段は、前記現フレームの画像における各画素を中心とした所定の範囲の画素を含む標本集合から有意確率を算出することを特徴とする請求項３乃至７の何れか一項に記載の画像処理装置。
【請求項９】
前記有意確率は、前記標本集合の母平均が前記ノイズの統計量に基づいた確率分布に従うか否かを表すものであることを特徴とする請求項８に記載の画像処理装置。
【請求項１０】
前記有意確率は、前記標本集合の母分散が前記ノイズの統計量に基づいた確率分布に従う否かを表すものであることを特徴とする請求項８に記載の画像処理装置。
【請求項１１】
前記有意確率は、前記標本集合の母平均および母分散が前記ノイズの統計量に基づいた確率分布に従うか否かを表すものであることを特徴とする請求項８に記載の画像処理装置。
【請求項１２】
前記現フレームの画像データを空間方向に空間フィルタ処理する空間フィルタ手段を更に備え、
前記加算手段は前記空間フィルタ処理後の画像データを加算することを特徴とする請求項１乃至１１の何れか一項に記載の画像処理装置。
【請求項１３】
前記空間フィルタ手段の前記空間フィルタ処理の処理効果を前記評価値に応じて変更することを特徴とする請求項１２に記載の画像処理装置。
【請求項１４】
複数フレームの動画像を時系列に重み付け加算する画像処理方法であって、
前記複数フレームの動画像の中で、現フレームより前のフレーム画像をリカ―シブフィルタ処理するフィルタ工程と、
前記現フレームの画像と前記フィルタ工程で得られた画像との差分画像におけるノイズの統計量に基づいて現フレームの重み付け係数を得る係数取得工程と、
前記係数取得手段で得られた前記複数フレームに対応する重み付け係数に基づき前記複数フレームの動画像を加算する加算工程と、
を有することを徴とする画像処理方法。
【請求項１５】
請求項１４に記載の画像処理方法の各工程をコンピュータに実行させるためのプログラム。
【請求項１６】
請求項１４に記載の画像処理方法の各工程をコンピュータに実行させるためのプログラムを記憶した記憶媒体。
【請求項１】
複数フレームの動画像を時系列に重み付け加算する画像処理装置であって、
前記複数フレームの動画像の中で、現フレームより前のフレーム画像をリカ―シブフィルタ処理するフィルタ手段と、
前記現フレームの画像と前記フィルタ手段で得られた画像との差分画像におけるノイズの統計量に基づいて現フレームの重み付け係数を得る係数取得手段と、
前記係数取得手段で得られた前記複数フレームに対応する重み付け係数に基づき前記複数フレームの動画像を加算する加算手段と、
を有することを徴とする画像処理装置。
【請求項２】
前記係数取得手段は、予めえられたノイズの分布に従うかを示す評価値に基づいて前記重み係数を得ることを特徴とする請求項１に記載の画像処理装置。
【請求項３】
前記ノイズの統計量が予め得られたノイズの確率分布に従うか否かを表す有意確率を求める検定手段を更に備え、
前記加算手段は、前記有意確率に基づいて前記重みづけ係数を得ることを特徴とする請求項１又は２に記載の画像処理装置。
【請求項４】
複数の異なる周波数帯域の画像データに分解する分解手段を更に備え、
前記画像データは前記分解手段で分解された画像であることを特徴とする請求項１乃至３のいずれか一項に記載の画像処理装置。
【請求項５】
複数の異なる周波数帯域の画像データに再構成する再構成手段を更に備え、
前記再構成手段は、前記加算手段で加算した画像データに基づいて再構成することを特徴とする請求項４に記載の画像処理装置。
【請求項６】
前記差分画像におけるノイズの統計量を取得する統計量取得手段を更に備え、
前記検定手段は、前記ノイズの統計量が予め取得された確率分布に従うか否かを示す有意確率を求めることを特徴とする請求項３乃至５のいずれか一項に記載の画像処理装置。
【請求項７】
前記ノイズの統計量は、ノイズの平均値、標準偏差、分散の少なくとも何れか１つであることを特徴とする請求項１乃至６の何れか一項に記載の画像処理装置。
【請求項８】
前記検定手段は、前記現フレームの画像における各画素を中心とした所定の範囲の画素を含む標本集合から有意確率を算出することを特徴とする請求項３乃至７の何れか一項に記載の画像処理装置。
【請求項９】
前記有意確率は、前記標本集合の母平均が前記ノイズの統計量に基づいた確率分布に従うか否かを表すものであることを特徴とする請求項８に記載の画像処理装置。
【請求項１０】
前記有意確率は、前記標本集合の母分散が前記ノイズの統計量に基づいた確率分布に従う否かを表すものであることを特徴とする請求項８に記載の画像処理装置。
【請求項１１】
前記有意確率は、前記標本集合の母平均および母分散が前記ノイズの統計量に基づいた確率分布に従うか否かを表すものであることを特徴とする請求項８に記載の画像処理装置。
【請求項１２】
前記現フレームの画像データを空間方向に空間フィルタ処理する空間フィルタ手段を更に備え、
前記加算手段は前記空間フィルタ処理後の画像データを加算することを特徴とする請求項１乃至１１の何れか一項に記載の画像処理装置。
【請求項１３】
前記空間フィルタ手段の前記空間フィルタ処理の処理効果を前記評価値に応じて変更することを特徴とする請求項１２に記載の画像処理装置。
【請求項１４】
複数フレームの動画像を時系列に重み付け加算する画像処理方法であって、
前記複数フレームの動画像の中で、現フレームより前のフレーム画像をリカ―シブフィルタ処理するフィルタ工程と、
前記現フレームの画像と前記フィルタ工程で得られた画像との差分画像におけるノイズの統計量に基づいて現フレームの重み付け係数を得る係数取得工程と、
前記係数取得手段で得られた前記複数フレームに対応する重み付け係数に基づき前記複数フレームの動画像を加算する加算工程と、
を有することを徴とする画像処理方法。
【請求項１５】
請求項１４に記載の画像処理方法の各工程をコンピュータに実行させるためのプログラム。
【請求項１６】
請求項１４に記載の画像処理方法の各工程をコンピュータに実行させるためのプログラムを記憶した記憶媒体。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【公開番号】特開２０１３−４８７８２（Ｐ２０１３−４８７８２Ａ）
【公開日】平成２５年３月１４日（２０１３．３．１４）
【国際特許分類】
【出願番号】特願２０１１−１８９０２４（Ｐ２０１１−１８９０２４）
【出願日】平成２３年８月３１日（２０１１．８．３１）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
【Ｆターム（参考）】
【公開日】平成２５年３月１４日（２０１３．３．１４）
【国際特許分類】
【出願日】平成２３年８月３１日（２０１１．８．３１）
【出願人】（０００００１００７）キヤノン株式会社 (59,756)
【Ｆターム（参考）】
[ Back to top ]