JP6524497B2 - 画像処理装置、画像処理方法、及び、画像処理プログラム - Google Patents

Description

本発明は、画像処理装置、画像処理方法、及び、画像処理プログラムに関する。

複数の異なる時点にてそれぞれ撮影された複数の非透視画像を処理する画像処理装置が知られている（例えば、非特許文献１及び非特許文献２を参照）。非透視画像に含まれる画素は、前景又は背景を表す。例えば、前景は、車両等の移動する物体であり、背景は、建築物及び道路等の静止した物体である。

この画像処理装置は、混合正規分布に基づいて、画素が、前景及び背景のいずれを表すかを判定する。
具体的には、画像処理装置は、複数の非透視画像における、画素の輝度が複数の値のうちの各値を有する非透視画像の出現頻度に基づいて、当該画素の輝度の確率分布を取得する。換言すると、画素の輝度の確率分布は、複数の非透視画像のうちの、当該画素の輝度が複数の値のうちの各値を有する非透視画像の数の、非透視画像の総数に対する割合を表すとともに、当該画素の輝度を確率変数として用いる。
更に、画像処理装置は、取得した確率分布を表す混合正規分布を推定し、推定した混合正規分布に基づいて、画素が、前景及び背景のいずれを表すかを判定する。

例えば、図１の（Ａ）に表されるように、画像に含まれる画素の輝度は、画像が撮影された時点の変化に伴って変化する。この場合、画素の輝度の確率分布は、図１の（Ｂ）のように表される。

ところで、非透視画像を撮影する方向に沿って延びる直線上に複数の物体が存在する場合、非透視画像に含まれる画素は、最も手前側に位置する物体を表す。従って、前景が移動する物体である場合、画素が背景を表す時間は、当該画素が前景を表す時間よりも長くなりやすい。

そこで、画像処理装置は、画素の輝度が、推定した混合正規分布を構成する複数の正規分布Ｇ１，Ｇ２のうちの、確率が最大である正規分布Ｇ１に対応する場合、当該画素が背景を表すと判定する。

Ｊ．Ｃｈｏｉ、外２名、「Ａｄａｐｔｉｖｅ Ｓｈａｄｏｗ Ｅｓｔｉｍａｔｏｒ ｆｏｒ Ｒｅｍｏｖｉｎｇ Ｓｈａｄｏｗ ｏｆ Ｍｏｖｉｎｇ Ｏｂｊｅｃｔ」、Ｃｏｍｐｕｔｅｒ Ｖｉｓｉｏｎ ａｎｄ Ｉｍａｇｅ Ｕｎｄｅｒｓｔａｎｄｉｎｇ、Ｅｌｓｅｖｉｅｒ Ｉｎｃ．、２０１０年、第１１４号、ｐ．１０１７−１０２９ Ｓ．Ｃａｌｄｅｒａｒａ、外２名、「Ｖｉｓｉｏｎ Ｂａｓｅｄ Ｓｍｏｋｅ Ｄｅｔｅｃｔｉｏｎ Ｓｙｓｔｅｍ Ｕｓｉｎｇ Ｉｍａｇｅ Ｅｎｅｒｇｙ ａｎｄ Ｃｏｌｏｒ Ｉｎｆｏｒｍａｔｉｏｎ」、Ｍａｃｈｉｎｅ Ｖｉｓｉｏｎ ａｎｄ Ａｐｐｌｉｃａｔｉｏｎｓ、Ｓｐｒｉｎｇｅｒ、２０１０年、第２２号、ｐ．７０５−７１９

この画像処理装置を用いて透視画像を処理することにより、透視画像において対象物体（例えば、腫瘍、軟部組織、患部、又は、病変組織等）を識別することが考えられる。しかしながら、透視画像を撮影する方向に沿って延びる直線上に複数の物体が存在する場合、透視画像に含まれる画素の輝度は、当該複数の物体のそれぞれに起因する輝度の成分の重畳である。例えば、複数の物体のそれぞれに起因する輝度の成分の重畳は、複数の物体のそれぞれに起因する輝度の成分の和である。

このため、上述した画像処理装置を用いて透視画像を処理しても、透視画像において対象物体を高い精度にて識別することが困難であった。

本発明の目的の一つは、透視画像において対象物体を高い精度にて識別することが可能な画像処理装置を提供することにある。

一つの側面では、画像処理装置は、
複数の物体のそれぞれに起因する輝度の成分の重畳である画素の輝度を含む、複数の透視画像における、上記画素の輝度が複数の値のうちの各値を有する透視画像の出現頻度に基づいて、当該画素の輝度の確率分布を取得する取得部と、
上記取得された確率分布を表す混合分布を推定する分布推定部と、
上記推定された混合分布を構成する複数の部分分布のうちの期待値が最小である期待値最小分布に対応し且つ上記画素の輝度の少なくとも一部を構成する成分を推定する成分推定部と、
上記透視画像に含まれる上記画素の輝度から、上記推定された成分を除去する除去部と、を備える。

他の一つの側面では、画像処理方法は、
複数の物体のそれぞれに起因する輝度の成分の重畳である画素の輝度を含む、複数の透視画像における、前記画素の輝度が複数の値のうちの各値を有する透視画像の出現頻度に基づいて、当該画素の輝度の確率分布を取得し、
前記取得された確率分布を表す混合分布を推定し、
前記推定された混合分布を構成する複数の部分分布のうちの期待値が最小である期待値最小分布に対応し且つ前記画素の輝度の少なくとも一部を構成する成分を推定し、
前記透視画像に含まれる前記画素の輝度から、前記推定された成分を除去する。

他の一つの側面では、画像処理プログラムは、
複数の物体のそれぞれに起因する輝度の成分の重畳である画素の輝度を含む、複数の透視画像における、前記画素の輝度が複数の値のうちの各値を有する透視画像の出現頻度に基づいて、当該画素の輝度の確率分布を取得し、
前記取得された確率分布を表す混合分布を推定し、
前記推定された混合分布を構成する複数の部分分布のうちの期待値が最小である期待値最小分布に対応し且つ前記画素の輝度の少なくとも一部を構成する成分を推定し、
前記透視画像に含まれる前記画素の輝度から、前記推定された成分を除去する、
処理をコンピュータに実行させる。

透視画像において対象物体を高い精度にて識別できる。

輝度の時間に対する変化の一例と、輝度に対する確率分布の一例と、を表すグラフである。 第１実施形態の放射線照射システムの構成の一例を表すブロック図である。 図２の制御装置の構成の一例を表すブロック図である。 図２の制御装置の機能の一例を表すブロック図である。 図２の制御装置に入力される透視画像の一例を表す図である。 輝度の時間に対する変化の一例を表すグラフである。 輝度に対する確率分布の一例を表すグラフである。 推定される除去対象成分の時間に対する変化の一例を表すグラフである。 入力透視画像からの除去対象成分の除去の一例を表す説明図である。 図２の制御装置が実行する処理の一例を表すフローチャートである。 模擬透視画像の一例を表す説明図である。 模擬透視画像からの除去対象成分の除去結果の一例を表す図である。 模擬透視画像の各画素の輝度に対する二乗平均平方根誤差の時間に対する変化の一例を表すグラフである。 臨床透視画像に対する対象物体の位置の推定結果の一例を表す図である。 臨床透視画像の対象物体の位置に対する二乗平均平方根誤差の時間に対する変化の一例を表すグラフである。

本願の発明者は、透視画像に含まれる画素の輝度の確率分布を表す混合分布を構成する複数の分布のうちの期待値が最小である期待値最小分布が、対象物体と異なる物体に対応する確率が高い、という知見を得た。

そこで、後述する各実施形態において、画像処理装置は、期待値最小分布に対応し且つ画素の輝度の少なくとも一部を構成する成分を推定し、透視画像に含まれる画素の輝度から、推定した成分を除去する。
これによれば、透視画像に含まれる画素の輝度の中から、対象物体に起因する成分を高い精度にて抽出できる。この結果、透視画像において対象物体を高い精度にて識別できる。

以下、本発明の、画像処理装置、画像処理方法、及び、画像処理プログラム、に関する各実施形態について図２乃至図１５を参照しながら説明する。

＜第１実施形態＞
（構成）
図２に表されるように、第１実施形態の放射線照射システム１は、撮影装置１０と、制御装置２０と、放射線照射装置３０と、を備える。
撮影装置１０は、所定の撮影周期が経過する毎に生体ＬＢの透視画像を撮影する。本例では、撮影装置１０は、Ｘ線を用いて透視画像を撮影する。本例では、透視画像において、物体に起因する輝度は、当該物体に対するＸ線の透過率が低くなるほど高くなる。

撮影装置１０は、透視画像を撮影する毎に、撮影された透視画像を制御装置２０へ出力する。本例では、制御装置２０へ出力される透視画像は、Ｘ軸方向においてＮ_ｘｍａｘ個の画素が並ぶとともに、Ｘ軸方向に直交するＹ軸方向においてＮ_ｙｍａｘ個の画素が並ぶように、格子状に配列された、Ｎ_ｘｍａｘＮ_ｙｍａｘ個の画素のそれぞれの輝度を含む。Ｎ_ｘｍａｘ及びＮ_ｙｍａｘは、自然数を表す。

制御装置２０には、撮影装置１０から、上記撮影周期の時間間隔を有する複数の異なる時点にてそれぞれ撮影された複数の透視画像が入力される。制御装置２０に入力された透視画像は、入力透視画像とも表される。制御装置２０は、複数の入力透視画像を処理し、処理後の透視画像に基づいて生体ＬＢにおける対象物体の位置を推定する。例えば、対象物体は、腫瘍、軟部組織、患部、又は、病変組織等である。対象物体は、特定の物体と表されてもよい。
なお、複数の入力透視画像の処理については後述する。

制御装置２０は、推定した対象物体の位置に基づいて、放射線照射装置３０による生体ＬＢへの放射線の照射を制御する。
本例では、制御装置２０は、推定した対象物体の位置に基づいて、放射線を生体ＬＢに照射する期間を制御する。なお、制御装置２０は、推定した対象物体の位置に基づいて、生体ＬＢに照射する放射線の強度を制御してもよい。また、制御装置２０は、推定した対象物体の位置に基づいて、放射線を照射する位置を制御してもよい。

放射線照射装置３０は、制御装置２０による制御に従って、放射線を生体ＬＢに照射する。

ここで、制御装置２０について説明を加える。
制御装置２０は、図３に表されるように、処理装置２１と、記憶装置２２と、を備える。
処理装置２１は、記憶装置２２に記憶されたプログラムを実行することにより、後述する機能を実現する。本例では、処理装置２１は、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）である。なお、処理装置２１は、ＤＳＰ（Ｄｉｇｉｔａｌ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｏｒ）、又は、プログラム可能な論理回路装置（ＰＬＤ；Ｐｒｏｇｒａｍｍａｂｌｅ Ｌｏｇｉｃ Ｄｅｖｉｃｅ）により構成されていてもよい。

記憶装置２２は、情報を読み書き可能に記憶する。例えば、記憶装置２２は、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）、ＨＤＤ（Ｈａｒｄ Ｄｉｓｋ Ｄｒｉｖｅ）、ＳＳＤ（Ｓｏｌｉｄ Ｓｔａｔｅ Ｄｉｓｋ）、半導体メモリ、及び、有機メモリの少なくとも１つを備える。なお、記憶装置２２は、フレキシブルディスク、光ディスク、光磁気ディスク、及び、半導体メモリ等の記録媒体と、記録媒体から情報を読み取り可能な読取装置と、を備えていてもよい。
なお、制御装置２０は、集積回路（例えば、ＬＳＩ（Ｌａｒｇｅ Ｓｃａｌｅ Ｉｎｔｅｇｒａｔｉｏｎ）等）により実現されてよい。
制御装置２０は、画像処理装置の一例である。

図４に表されるように、制御装置２０の機能は、補正部２０１と、取得部２０２と、分布推定部２０３と、成分推定部２０４と、除去部２０５と、位置推定部２０６と、照射制御部２０７と、を備える。

本例では、制御装置２０は、透視画像に含まれるＮ_ｘｍａｘＮ_ｙｍａｘ個の画素のうちの、所定の対象領域に含まれる複数の画素のそれぞれの輝度から除去対象成分を除去する除去処理を、所定の反復回数Ｍだけ繰り返し実行する。Ｍは、自然数を表す。例えば、除去対象成分は、輝度のうちの、対象物体と異なる物体（例えば、骨又は筋肉等）に起因する成分である。
なお、対象物体は、骨又は筋肉であってもよい。この場合、対象物体と異なる物体は、軟部組織であってよい。

対象領域は、Ｘ軸方向における範囲が、ｘ_０＋１番目の画素から、ｘ_０＋Ｎ_ｘ番目の画素までであり、且つ、Ｙ軸方向における範囲が、ｙ_０＋１番目の画素から、ｙ_０＋Ｎ_ｙ番目の画素までである。ｘ_０及びｙ_０は、１以上の整数を表す。Ｎ_ｘ及びＮ_ｙは、２以上の整数を表す。従って、対象領域に含まれる画素の数は、Ｎ_ｘＮ_ｙ個である。本例では、対象領域は、放射線照射システム１のユーザによって設定される。なお、対象領域は、除去処理が実行された回数に応じて変更されてよい。例えば、対象領域は、除去処理が実行される毎に狭くなるように変更されてよい。

１回目の除去処理において、制御装置２０は、数式１に表されるように、除去処理の対象となる透視画像として、撮影装置１０からの入力透視画像を用いる。

Ｉ_ｍ（ｘ，ｙ，ｔ）は、ｍ回目の除去処理の対象となるＴ個の透視画像のうちのｔ番目の透視画像に含まれる複数の画素のうちの、Ｘ軸方向においてｘ番目であり且つＹ軸方向においてｙ番目である画素の輝度を表す。ｍは、１からＭまでの各整数を表す。

数式１において、ｘは、１からＮ_ｘｍａｘまでの各整数を表す。ｙは、１からＮ_ｙｍａｘまでの各整数を表す。Ｔは、２以上の整数を表す。ｔは、１からＴまでの各整数を表す。本例では、ｔ＋１番目の透視画像は、ｔ番目の透視画像が撮影された時点から、撮影周期だけ後の時点にて撮影された透視画像である。

Ｉ_ｐ（ｘ，ｙ，ｔ）は、Ｔ個の入力透視画像のうちのｔ番目の入力透視画像に含まれる複数の画素のうちの、Ｘ軸方向においてｘ番目であり且つＹ軸方向においてｙ番目である画素の輝度を表す。

本例では、２回目以降の除去処理において、制御装置２０は、除去処理の対象となる透視画像として、直前の除去処理において除去対象成分が除去された後の透視画像を用いる。なお、制御装置２０は、２回目以降の少なくとも１つの除去処理において、除去処理の対象となる透視画像として、直前の除去処理における除去対象成分からなる透視画像を用いてもよい。

補正部２０１は、ｍ回目の除去処理の対象となるＴ個の透視画像の中から基準透視画像を選択する。本例では、補正部２０１は、１番目の透視画像を基準透視画像として選択する。そして、補正部２０１は、選択した基準透視画像に対する、Ｔ個の透視画像のそれぞれにおける被写体の移動距離及び移動方向を推定する。本例では、被写体は、１番目の透視画像の、対象領域に含まれる複数の画素により表される。被写体は、テンプレート画像と表されてもよい。

本例では、補正部２０１は、テンプレートマッチング法を用いることにより、被写体の移動距離及び移動方向を推定する。補正部２０１は、テンプレートマッチング法において、相関係数を評価関数として用いてよい。本例では、被写体の移動距離及び移動方向は、Ｘ軸方向における移動量ｕ_ｍ（ｔ）と、Ｙ軸方向における移動量ｖ_ｍ（ｔ）と、により表される。ｕ_ｍ（ｔ）及びｖ_ｍ（ｔ）は、整数を表す。

補正部２０１は、Ｔ個の透視画像のそれぞれに対して、推定した移動方向と逆方向へ、推定した移動距離だけ被写体を移動させるように当該透視画像を被写体の位置に関して補正する。本例では、補正部２０１は、数式２に表されるように、Ｔ個の透視画像のそれぞれに含まれる複数の画素のそれぞれの輝度を補正する。Ｉ_ｍ’（ｘ，ｙ，ｔ）は、補正後の輝度を表す。なお、画素の輝度の補正は、当該画素を含む画像の補正の一例である。

例えば、被写体の移動距離及び移動方向は、対象領域に含まれる複数の画素のそれぞれに対する移動量ベクトルに基づく主要移動量ベクトルにより表されてよい。例えば、画素に対する移動量ベクトルは、画素の移動距離及び移動方向を表す。例えば、主要移動量ベクトルは、平均移動量ベクトル、又は、最頻移動量ベクトルである。平均移動量ベクトルは、対象領域に含まれる複数の画素に対して、各画素の移動量ベクトルを平均することにより得られる。最頻移動量ベクトルは、対象領域に含まれる複数の画素において、最も出現頻度が高い移動量ベクトルである。

取得部２０２は、補正部２０１による補正後の輝度Ｉ_ｍ’（ｘ，ｙ，ｔ）に基づいて、対象領域に含まれる複数の画素のそれぞれに対して、Ｔ個の透視画像のうちの、当該画素の輝度が複数の値のうちの各値を有する透視画像の数を計数することにより、輝度に対する度数分布を取得する。

取得部２０２は、対象領域に含まれる複数の画素のそれぞれに対して、取得した度数分布に基づいて、当該画素の輝度を確率変数として用いる確率分布を取得する。本例では、確率分布は、Ｔ個の透視画像のうちの、画素の輝度が複数の値のうちの各値を有する透視画像の数の、透視画像の総数Ｔに対する割合を表す。

図５は、透視画像の一例を表す。透視画像の中央部に位置し、円形状を有し、且つ、輝度が他の部分よりも高い部分ＣＣの位置が、透視画像間で変化する場合を想定する。更に、ある画素の輝度が、図６に表されるように、時間の経過に伴って変化する場合を想定する。この場合、取得部２０２は、図７に表されるように、確率分布ＰＤを取得する。

分布推定部２０３は、対象領域に含まれる複数の画素のそれぞれに対して、取得部２０２により取得された確率分布を表す混合分布ｐ（ｉ_ｍ（ｘ，ｙ））を推定する。混合分布は、混合確率分布と表されてもよい。本例では、分布推定部２０３は、ＥＭ（Ｅｘｐｅｃｔａｔｉｏｎ Ｍａｘｉｍｉｚａｔｉｏｎ）アルゴリズムに従って数式３に表されるような混合正規分布を推定する。混合正規分布は、混合分布の一例である。

ｉ_ｍ（ｘ，ｙ）は、Ｘ軸方向においてｘ番目であり且つＹ軸方向においてｙ番目である画素の輝度を表す確率変数である。数式３において、ｘは、ｘ_０＋１からｘ_０＋Ｎ_ｘまでの各整数を表す。ｙは、ｙ_０＋１からｙ_０＋Ｎ_ｙまでの各整数を表す。

Ｎ（ｉ｜μ，σ^２）は、期待値がμであり、分散がσ^２であり、且つ、確率変数がｉである正規分布を表す。期待値は、平均値と表されてもよい。正規分布は、部分分布の一例である。なお、部分分布として正規分布と異なる分布が用いられてもよい。Ｊ（ｘ，ｙ）は、Ｘ軸方向においてｘ番目であり且つＹ軸方向においてｙ番目である画素の輝度に対する混合分布を構成する部分分布の数を表す。ｊは、１からＪ（ｘ，ｙ）までの各整数を表す。π_ｊ（ｘ，ｙ）は、上記混合分布を構成するＪ（ｘ，ｙ）個の部分分布のうちのｊ番目の部分分布に対する混合係数を表す。μ_ｊ（ｘ，ｙ）は、上記混合分布を構成するＪ（ｘ，ｙ）個の部分分布のうちのｊ番目の部分分布の期待値を表す。σ_ｊ ^２（ｘ，ｙ）は、上記混合分布を構成するＪ（ｘ，ｙ）個の部分分布のうちのｊ番目の部分分布の分散を表す。

図７に表されるように、取得部２０２が確率分布ＰＤを取得した場合、分布推定部２０３は、２個の部分分布Ｇ１，Ｇ２により構成される混合分布ＭＧを推定する。

成分推定部２０４は、数式４に表されるように、対象領域に含まれる複数の画素のそれぞれに対して、分布推定部２０３により推定された混合分布を構成するＪ（ｘ，ｙ）個の部分分布のうちの期待値が最小である期待値最小分布を決定する。

θ（ｘ，ｙ）は、Ｘ軸方向においてｘ番目であり且つＹ軸方向においてｙ番目である画素の輝度に対する混合分布を構成するＪ（ｘ，ｙ）個の部分分布のうちのθ（ｘ，ｙ）番目の部分分布が期待値最小分布であることを表す。θ（ｘ，ｙ）は、期待値最小分布のインデックスと表されてもよい。数式４において、ｘは、ｘ_０＋１からｘ_０＋Ｎ_ｘまでの各整数を表す。ｙは、ｙ_０＋１からｙ_０＋Ｎ_ｙまでの各整数を表す。

更に、成分推定部２０４は、Ｔ個の透視画像のそれぞれの対象領域に含まれる複数の画素のそれぞれに対して、決定した期待値最小分布に対応し、且つ、当該画素の輝度の少なくとも一部を構成する成分を推定する。画素の輝度に対する期待値最小分布に対応し、且つ、当該画素の輝度の少なくとも一部を構成する成分は、除去対象成分と表されてもよい。

本例では、成分推定部２０４は、除去対象成分を下記のように推定する。
先ず、成分推定部２０４は、補正部２０１による補正後の輝度Ｉ_ｍ’（ｘ，ｙ，ｔ）に基づいて、Ｔ個の透視画像のそれぞれの対象領域に含まれる複数の画素のそれぞれに対して、当該画素の輝度に対する混合分布を構成するＪ（ｘ，ｙ）個の部分分布のそれぞれに、当該画素の輝度が起因する確率を算出する。

本例では、成分推定部２０４は、上記確率として、数式５に表される負担率γ_ｋ（ｘ，ｙ，ｔ）を用いる。数式５において、ｘは、ｘ_０＋１からｘ_０＋Ｎ_ｘまでの各整数を表す。ｙは、ｙ_０＋１からｙ_０＋Ｎ_ｙまでの各整数を表す。

負担率γ_ｋ（ｘ，ｙ，ｔ）は、ｔ番目の透視画像の、Ｘ軸方向においてｘ番目であり且つＹ軸方向においてｙ番目である画素の輝度がＩ_ｍ’（ｘ，ｙ，ｔ）であるという条件の下で、当該輝度がｋ番目の部分分布に起因する確率である。ｋは、１からＪ（ｘ，ｙ）までの各整数を表す。

ｚ_ｋは、０又は１を有する２値確率変数である。ｚ_ｋが１を有することは、輝度がｋ番目の部分分布に起因することを表す。｛ｚ_１，・・・，ｚ_{Ｊ（ｘ，ｙ）}｝は、１−ｏｆ−Ｊ（ｘ，ｙ）表現に従ったＪ（ｘ，ｙ）次元のベクトルである。換言すると、ｚ_１，・・・，ｚ_{Ｊ（ｘ，ｙ）}は、いずれか１つが１を有するとともに、他のすべてが０を有する。

成分推定部２０４は、数式６に表されるように、Ｔ個の透視画像のそれぞれの対象領域に含まれる複数の画素のそれぞれに対して、当該画素の輝度に対する混合分布を構成するＪ（ｘ，ｙ）個の部分分布のうちの、期待値最小分布と異なる部分分布に対する負担率γ_ｊ（ｘ，ｙ，ｔ）のそれぞれが、当該期待値最小分布に対する負担率γ_{θ（ｘ，ｙ）}（ｘ，ｙ，ｔ）よりも小さい場合、除去対象成分ρ（ｘ，ｙ，ｔ）として、当該画素の輝度Ｉ_ｍ’（ｘ，ｙ，ｔ）を推定する。ここで、ｊは、１からＪ（ｘ，ｙ）までの整数のうちのθ（ｘ，ｙ）と異なる各整数を表す。

一方、成分推定部２０４は、数式７に表されるように、Ｔ個の透視画像のそれぞれの対象領域に含まれる複数の画素のそれぞれに対して、当該画素の輝度に対する混合分布を構成するＪ（ｘ，ｙ）個の部分分布のうちの、期待値最小分布と異なる部分分布に対する負担率γ_ｊ（ｘ，ｙ，ｔ）の少なくとも１つが、当該期待値最小分布に対する負担率γ_{θ（ｘ，ｙ）}（ｘ，ｙ，ｔ）以上である場合、除去対象成分ρ（ｘ，ｙ，ｔ）として、当該期待値最小分布の期待値μ_{θ（ｘ，ｙ）}（ｘ，ｙ）を推定する。ここで、ｊは、１からＪ（ｘ，ｙ）までの整数のうちのθ（ｘ，ｙ）と異なる各整数を表す。

期待値最小分布の期待値μ_{θ（ｘ，ｙ）}（ｘ，ｙ）は、期待値最小分布に基づく値の一例である。
なお、この場合、成分推定部２０４は、期待値最小分布の期待値μ_{θ（ｘ，ｙ）}（ｘ，ｙ）に代えて、期待値最小分布の期待値μ_{θ（ｘ，ｙ）}（ｘ，ｙ）及び分散σ_{θ（ｘ，ｙ）} ^２（ｘ，ｙ）に基づいて決定された値を除去対象成分ρ（ｘ，ｙ，ｔ）として推定してもよい。
このようにして、成分推定部２０４は、除去対象成分を推定する。

図７に表されるように、分布推定部２０３が、２つの部分分布Ｇ１，Ｇ２により構成される混合分布ＭＧを推定した場合、成分推定部２０４は、部分分布Ｇ１を期待値最小分布として決定する。更に、この場合、成分推定部２０４は、図８に表されるように、期間Ｔ１，Ｔ３において、除去対象成分ρ（ｘ，ｙ，ｔ）として、期待値最小分布の期待値μ_{θ（ｘ，ｙ）}（ｘ，ｙ）を推定する。また、この場合、成分推定部２０４は、期間Ｔ２，Ｔ４において、除去対象成分ρ（ｘ，ｙ，ｔ）として輝度Ｉ_ｍ’（ｘ，ｙ，ｔ）を推定する。

除去部２０５は、Ｔ個の透視画像のそれぞれの対象領域に含まれる複数の画素のそれぞれに対して、当該画素の輝度Ｉ_ｍ’（ｘ，ｙ，ｔ）から、成分推定部２０４により推定された除去対象成分ρ（ｘ，ｙ，ｔ）を除去する。

本例では、除去部２０５は、数式８に表されるように、Ｔ個の透視画像のそれぞれの対象領域に含まれる複数の画素のそれぞれに対して、当該画素の輝度Ｉ_ｍ’（ｘ，ｙ，ｔ）から、成分推定部２０４により推定された除去対象成分ρ（ｘ，ｙ，ｔ）を減算する。Ｉ_ｍ’’（ｘ，ｙ，ｔ）は、減算後の輝度を表す。

更に、本例では、除去部２０５は、減算後の輝度Ｉ_ｍ’’（ｘ，ｙ，ｔ）を、ｍ＋１回目の除去処理の対象となる透視画像の輝度Ｉ_ｍ＋１（ｘ，ｙ，ｔ）として設定する。なお、除去部２０５は、ｍが２以上である場合において、成分推定部２０４により推定された除去対象成分ρ（ｘ，ｙ，ｔ）を、ｍ＋１回目の除去処理の対象となる透視画像の輝度Ｉ_ｍ＋１（ｘ，ｙ，ｔ）として設定してもよい。

本例では、除去部２０５は、Ｍ回の除去処理によって除去対象成分が除去された後のＴ個の透視画像を位置推定部２０６へ出力する。なお、除去部２０５は、Ｍが２以上である場合、Ｍ回目の除去処理における除去対象成分ρ（ｘ，ｙ，ｔ）からなる透視画像を位置推定部２０６へ出力してもよい。

このように、除去処理は、補正部２０１、取得部２０２、分布推定部２０３、成分推定部２０４、及び、除去部２０５により実行される。

図９の（Ａ）は、除去処理の対象となる透視画像の一例を表す。図９の（Ｂ）は、対象物体と異なる物体に起因する輝度の成分としての除去対象成分の一例を表す。図９の（Ｃ）は、除去処理によって除去対象成分が除去された後の透視画像の一例を表す。この場合、図９の（Ｃ）に表されるように、除去対象成分が除去された後の透視画像によれば、対象物体を高い精度にて識別できる。また、除去対象成分が対象物体に起因する輝度の成分である場合、除去対象成分が除去された後の透視画像に代えて、除去対象成分からなる透視画像を用いることにより、対象物体を高い精度にて識別できる。

位置推定部２０６には、Ｍ回の除去処理によって除去対象成分が除去された後のＴ個の透視画像が除去部２０５から入力される。なお、位置推定部２０６には、Ｍが２以上である場合、Ｍ回目の除去処理における除去対象成分ρ（ｘ，ｙ，ｔ）からなる透視画像が除去部２０５から入力されてもよい。

位置推定部２０６は、除去部２０５から入力されたＴ個の透視画像に基づいて、生体ＬＢにおける対象物体の位置を推定する。例えば、位置推定部２０６は、テンプレートマッチング法を用いることにより、対象物体の位置を推定する。

照射制御部２０７は、位置推定部２０６により推定された対象物体の位置に基づいて、放射線照射装置３０による生体ＬＢへの放射線の照射を制御する。

（動作）
次に、放射線照射システム１の動作について説明する。
撮影装置１０は、撮影周期が経過する毎に生体ＬＢの透視画像を撮影し、透視画像を撮影する毎に、撮影された透視画像を制御装置２０へ出力する。

制御装置２０は、撮影装置１０から入力された複数の入力透視画像を処理する。本例では、制御装置２０は、図１０にフローチャートにより表される処理を、所定の実行周期が経過する毎に実行する。

具体的には、制御装置２０は、処理回数ｍを初期値１に設定する（図１０のステップＳ１０１）。次いで、制御装置２０は、ｍ回目の除去処理の対象となる透視画像に含まれる画素の輝度Ｉ_ｍ（ｘ，ｙ，ｔ）を、入力透視画像に含まれる画素の輝度Ｉ_ｐ（ｘ，ｙ，ｔ）に設定する（図１０のステップＳ１０２）。

次いで、制御装置２０は、Ｔ個の透視画像の中から選択された基準透視画像に対する、Ｔ個の透視画像のそれぞれにおける被写体の移動距離及び移動方向を推定する（図１０のステップＳ１０３）。そして、制御装置２０は、Ｔ個の透視画像のそれぞれに対して、推定した移動方向と逆方向へ、推定した移動距離だけ被写体を移動させるように当該透視画像を補正する（図１０のステップＳ１０４）。

次いで、制御装置２０は、ステップＳ１０４にて補正された輝度Ｉ_ｍ’（ｘ，ｙ，ｔ）に基づいて、対象領域に含まれる複数の画素のそれぞれに対して、当該画素の輝度を確率変数として用いる確率分布を取得する（図１０のステップＳ１０５）。

そして、制御装置２０は、対象領域に含まれる複数の画素のそれぞれに対して、ステップＳ１０５にて取得された確率分布を表す混合分布ｐ（ｉ_ｍ（ｘ，ｙ））を推定する（図１０のステップＳ１０６）。

次いで、制御装置２０は、対象領域に含まれる複数の画素のそれぞれに対して、ステップＳ１０６にて推定された混合分布を構成するＪ（ｘ，ｙ）個の部分分布のうちの期待値が最小である期待値最小分布を決定する（図１０のステップＳ１０７）。

そして、制御装置２０は、Ｔ個の透視画像のそれぞれの対象領域に含まれる複数の画素のそれぞれに対して、当該画素の輝度に対する混合分布を構成するＪ（ｘ，ｙ）個の部分分布のそれぞれに、当該画素の輝度Ｉ_ｍ’（ｘ，ｙ，ｔ）が起因する確率を表す負担率γ_ｋ（ｘ，ｙ，ｔ）を算出する（図１０のステップＳ１０８）。

次いで、制御装置２０は、Ｔ個の透視画像のそれぞれの対象領域に含まれる複数の画素のそれぞれに対して、除去対象成分ρ（ｘ，ｙ，ｔ）を推定する（図１０のステップＳ１０９）。

本例では、制御装置２０は、期待値最小分布に対する負担率γ_{θ（ｘ，ｙ）}（ｘ，ｙ，ｔ）が、期待値最小分布と異なる部分分布に対する負担率γ_ｊ（ｘ，ｙ，ｔ）よりも高い場合、当該画素の輝度Ｉ_ｍ’（ｘ，ｙ，ｔ）を除去対象成分ρ（ｘ，ｙ，ｔ）として推定する。
更に、本例では、制御装置２０は、期待値最小分布に対する負担率γ_{θ（ｘ，ｙ）}（ｘ，ｙ，ｔ）が、期待値最小分布と異なる部分分布に対する負担率γ_ｊ（ｘ，ｙ，ｔ）以下である場合、当該期待値最小分布の期待値μ_{θ（ｘ，ｙ）}（ｘ，ｙ）を除去対象成分ρ（ｘ，ｙ，ｔ）として推定する。

次いで、制御装置２０は、Ｔ個の透視画像のそれぞれの対象領域に含まれる複数の画素のそれぞれに対して、当該画素の輝度Ｉ_ｍ’（ｘ，ｙ，ｔ）から、ステップＳ１０９にて推定された除去対象成分ρ（ｘ，ｙ，ｔ）を除去することにより、ｍ＋１回目の除去処理の対象となる透視画像に含まれる画素の輝度Ｉ_ｍ＋１（ｘ，ｙ，ｔ）を設定する（図１０のステップＳ１１０）。透視画像に含まれる画素の輝度の設定は、透視画像の生成と表されてもよい。

そして、制御装置２０は、処理回数ｍに１を加算する（図１０のステップＳ１１１）。次いで、制御装置２０は、処理回数ｍが反復回数Ｍよりも大きいか否かを判定する（図１０のステップＳ１１２）。

処理回数ｍが反復回数Ｍ以下である場合、制御装置２０は、「Ｎｏ」と判定し、ステップＳ１０３へ戻り、処理回数ｍが反復回数Ｍよりも大きくなるまで、ステップＳ１０３からステップＳ１１２までの処理を繰り返し実行する。

処理回数ｍが反復回数Ｍよりも大きい場合、制御装置２０は、ステップＳ１１２にて「Ｙｅｓ」と判定し、ステップＳ１１３へ進む。そして、制御装置２０は、直前の除去処理において推定された除去対象成分ρ（ｘ，ｙ，ｔ）が除去された後の透視画像に基づいて、対象物体の位置を推定する（図１０のステップＳ１１３）。なお、制御装置２０は、直前の除去処理において推定された除去対象成分ρ（ｘ，ｙ，ｔ）からなる透視画像に基づいて、対象物体の位置を推定してもよい。

次いで、制御装置２０は、ステップＳ１１３にて推定された対象物体の位置に基づいて、放射線照射装置３０による放射線の照射を制御する（図１０のステップＳ１１４）。
その後、制御装置２０は、図１０の処理を終了する。

以上、説明したように、第１実施形態の制御装置２０は、画素毎の輝度の確率分布を表す混合分布を推定する。更に、制御装置２０は、推定した混合分布を構成する複数の部分分布のうちの期待値が最小である期待値最小分布に対応し且つ当該画素の輝度の少なくとも一部を構成する成分を推定する。加えて、制御装置２０は、透視画像に含まれる画素の輝度から、推定した成分を除去する。

これによれば、透視画像に含まれる画素の輝度の中から、対象物体に起因する成分を高い精度にて抽出できる。この結果、透視画像において対象物体を高い精度にて識別できる。

更に、第１実施形態の制御装置２０は、画素の輝度が期待値最小分布に起因する確率が、推定された混合分布を構成する複数の部分分布のうちの、期待値最小分布と異なる部分分布に当該画素の輝度が起因する確率よりも高い場合、当該画素の輝度を上記成分として推定する。

画素の輝度が期待値最小分布に起因する確率が、混合分布を構成する複数の部分分布のうちの、期待値最小分布と異なる部分分布に画素の輝度が起因する確率よりも高い場合、当該画素の輝度が、期待値最小分布に対応する物体に起因する成分のみを有する確率が高い。

制御装置２０は、この場合、画素の輝度を、期待値最小分布に対応する成分として推定する。これにより、透視画像に含まれる画素の輝度から、期待値最小分布に対応する物体に起因する成分を高い精度にて除去できる。この結果、透視画像において対象物体を高い精度にて識別できる。

更に、第１実施形態の制御装置２０は、画素の輝度が期待値最小分布に起因する確率が、推定された混合分布を構成する複数の部分分布のうちの、期待値最小分布と異なる部分分布に当該画素の輝度が起因する確率よりも低い場合、期待値最小分布に基づく値を上記成分として推定する。

画素の輝度が期待値最小分布に起因する確率が、混合分布を構成する複数の部分分布のうちの、期待値最小分布と異なる部分分布に画素の輝度が起因する確率よりも低い場合、当該画素の輝度が、期待値最小分布に対応する物体に起因する成分に加えて、他の物体に起因する成分も含む確率が高い。

制御装置２０は、この場合、期待値最小分布に基づく値（例えば、期待値最小分布の期待値）を、期待値最小分布に対応する成分として推定する。これにより、透視画像に含まれる画素の輝度から、期待値最小分布に対応する物体に起因する成分を適切に除去できる。この結果、透視画像において対象物体を高い精度にて識別できる。

更に、第１実施形態の制御装置２０は、複数の透視画像のそれぞれに対して、上記成分の除去を実行する。加えて、制御装置２０は、上記成分が除去された複数の透視画像に対して、確率分布の取得、混合分布の推定、成分の推定、及び、成分の除去を、繰り返し実行する。

これによれば、画素の輝度が、対象物体と異なる、複数の物体のそれぞれに起因する成分を含む場合に、透視画像に含まれる画素の輝度の中から、対象物体に起因する成分を高い精度にて抽出できる。この結果、透視画像において対象物体を高い精度にて識別できる。

更に、第１実施形態の制御装置２０は、複数の透視画像の中から選択された基準透視画像に対する、複数の透視画像のそれぞれにおける被写体の移動距離及び移動方向を推定する。加えて、制御装置２０は、複数の透視画像のそれぞれに対して、推定した移動方向と逆方向へ、推定した移動距離だけ被写体を移動させるように当該透視画像を補正する。

これによれば、透視画像に含まれる画素の輝度に対して被写体の移動が及ぼす影響を抑制できる。この結果、透視画像に含まれる画素の輝度の中から、対象物体に起因する成分を高い精度にて抽出できる。この結果、透視画像において対象物体を高い精度にて識別できる。

次に、臨床透視画像を模擬した模擬透視画像に対して、第１実施形態の制御装置２０が処理した結果について説明する。
模擬透視画像は、図１１に表されるように、第１乃至第３の基本画像を重ね合わせた画像である。図１１の（Ａ）は、模擬透視画像を表す。図１１の（Ｂ）乃至（Ｄ）は、第１乃至第３の基本画像をそれぞれ表す。本例では、複数の基本画像を重ね合わせた画像の各画素の輝度は、数式９に表されるように、当該複数の基本画像の当該画素の輝度の和である。

Ｉ_ｐ（ｘ，ｙ，ｔ）は、模擬透視画像に含まれる画素の輝度を表す。Ｌ_ｊ（ｘ，ｙ，ｔ）は、第ｊの基本画像に含まれる画素の輝度を表す。Ｊは、３を表す。ｊは、１からＪまでの各整数を表す。
第１の基本画像は、臨床透視画像に、期待値が０であり且つ分散が４であるガウスノイズを付加した画像である。
第２の基本画像は、長方形に、期待値が１５であり且つ分散が４であるガウスノイズを付加した画像である。
第３の基本画像は、腫瘍を模擬した円形に、期待値が１０であり且つ分散が１であるガウスノイズを付加した画像である。

第２の基本画像における長方形は、時間の経過に伴って、長方形のＸ軸方向における位置が正弦波に従って変化するように、Ｘ軸方向にて移動する。
第３の基本画像における円形は、時間の経過に伴って、円形のＹ軸方向における位置が正弦波に従って変化するように、Ｙ軸方向にて移動する。

図１２の（Ａ）は、各基本画像の真値を表す。図１２の（Ｂ）は、第１実施形態の制御装置２０による処理の結果を表す。図１２の（Ｃ）は、比較例の画像処理装置による処理の結果を表す。比較例の画像処理装置は、非特許文献（Ｙ．Ｗｅｉｓｓ、「Ｄｅｒｉｖｉｎｇ Ｉｎｔｒｉｎｓｉｃ Ｉｍａｇｅｓ ｆｒｏｍ Ｉｍａｇｅ Ｓｅｑｕｅｎｃｅｓ」、Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｃｏｎｆｅｒｅｎｃｅ ｏｎ Ｃｏｍｐｕｔｅｒ Ｖｉｓｉｏｎ、２００１年、第２巻、ｐ．６８−７５）に記載の固有画像抽出法（ＤＩＩ；Ｄｅｒｉｖｉｎｇ Ｉｎｔｒｉｎｓｉｃ Ｉｍａｇｅ）を用いる。

図１２の（Ｂ）の図面における最上段は、１回目の除去処理において推定された除去対象成分を表す。図１２の（Ｂ）の図面における中段は、２回目の除去処理において推定された除去対象成分を表す。図１２の（Ｂ）の図面における最下段は、２回目の除去処理において除去対象成分が除去された後の透視画像を表す。

第１実施形態の制御装置２０による２回目の除去処理において除去対象成分が除去された後の透視画像に含まれる画素の輝度の、第３の基本画像に含まれる画素の輝度に対する二乗平均平方根誤差Ｌ_ｒｍｓｅは、数式１０により表される。図１３は、二乗平均平方根誤差Ｌ_ｒｍｓｅの時間に対する変化を表す。第１実施形態の制御装置２０に対する二乗平均平方根誤差Ｌ_ｒｍｓｅは、比較例の画像処理装置よりも小さい。

このように、第１実施形態の制御装置２０によれば、透視画像に含まれる画素の輝度の中から、対象物体に起因する成分を高い精度にて抽出できる。この結果、透視画像において対象物体を高い精度にて識別できる。

次に、臨床透視画像に対して、第１実施形態の制御装置２０が処理した結果について説明する。
図１４の（Ａ）は、入力透視画像を表す。図１４の（Ｂ）は、比較例の画像処理装置による処理の結果を表す。図１４の（Ｃ）は、第１実施形態の制御装置２０による処理の結果を表す。図１４において、点線の長方形は、対象物体の真の位置を表す。図１４の（Ｂ）において、実線の長方形は、比較例の画像処理装置により推定された対象物体の位置を表す。図１４の（Ｃ）において、実線の長方形は、第１実施形態の制御装置２０により推定された対象物体の位置を表す。

第１実施形態の制御装置２０により推定された対象物体の位置の、対象物体の真の位置に対する二乗平均平方根誤差Ｐ_ｒｍｓｅは、数式１１により表される。ω_ｅ（ｔ）及びλ_ｅ（ｔ）は、対象物体のＸ軸方向及びＹ軸方向における、推定された位置をそれぞれ表す。ω_０（ｔ）及びλ_０（ｔ）は、対象物体のＸ軸方向及びＹ軸方向における真の位置をそれぞれ表す。図１５は、二乗平均平方根誤差Ｐ_ｒｍｓｅの時間に対する変化を表す。第１実施形態の制御装置２０に対する二乗平均平方根誤差Ｐ_ｒｍｓｅは、比較例の画像処理装置よりも小さい。

このように、第１実施形態の制御装置２０によれば、透視画像において対象物体の位置を高い精度にて推定できる。

第１実施形態において、制御装置２０は、放射線照射装置３０を制御する。ところで、制御装置２０は、放射線照射装置３０を制御しない画像処理装置であってもよい。この場合、制御装置２０は、除去対象成分が除去された後の透視画像を出力又は送信してもよい。

なお、本発明は、上述した実施形態に限定されない。例えば、上述した実施形態に、本発明の趣旨を逸脱しない範囲内において当業者が理解し得る様々な変更が加えられてよい。例えば、本発明の趣旨を逸脱しない範囲内において、上述した実施形態の他の変形例として、上述した実施形態及び変形例の任意の組み合わせが採用されてもよい。

１ 放射線照射システム
１０ 撮影装置
２０ 制御装置
２１ 処理装置
２２ 記憶装置
２０１ 補正部
２０２ 取得部
２０３ 分布推定部
２０４ 成分推定部
２０５ 除去部
２０６ 位置推定部
２０７ 照射制御部
３０ 放射線照射装置
ＬＢ 生体

Claims (7)

1. 複数の物体のそれぞれに起因する輝度の成分の重畳である画素の輝度を含む、複数の透視画像における、前記画素の輝度が複数の値のうちの各値を有する透視画像の出現頻度に基づいて、当該画素の輝度の確率分布を取得する取得部と、
   前記取得された確率分布を表す混合分布を推定する分布推定部と、
   前記推定された混合分布を構成する複数の部分分布のうちの期待値が最小である期待値最小分布に対応し且つ前記画素の輝度の少なくとも一部を構成する成分を推定する成分推定部と、
   前記透視画像に含まれる前記画素の輝度から、前記推定された成分を除去する除去部と、
   を備える、画像処理装置。
2. 請求項１に記載の画像処理装置であって、
   前記成分推定部は、
   前記画素の輝度が前記期待値最小分布に起因する確率が、前記推定された混合分布を構成する複数の部分分布のうちの、前記期待値最小分布と異なる部分分布に前記画素の輝度が起因する確率よりも高い場合、当該画素の輝度を前記成分として推定する、画像処理装置。
3. 請求項１又は請求項２に記載の画像処理装置であって、
   前記成分推定部は、
   前記画素の輝度が前記期待値最小分布に起因する確率が、前記推定された混合分布を構成する複数の部分分布のうちの、前記期待値最小分布と異なる部分分布に前記画素の輝度が起因する確率よりも低い場合、当該期待値最小分布に基づく値を前記成分として推定する、画像処理装置。
4. 請求項１乃至請求項３のいずれか一項に記載の画像処理装置であって、
   前記除去部は、前記複数の透視画像のそれぞれに対して、前記成分の除去を実行し、
   前記画像処理装置は、
   前記成分が除去された前記複数の透視画像に対して、前記確率分布の取得、前記混合分布の推定、前記成分の推定、及び、前記成分の除去を、繰り返し実行する、画像処理装置。
5. 請求項１乃至請求項４のいずれか一項に記載の画像処理装置であって、
   前記複数の透視画像の中から選択された基準透視画像に対する、前記複数の透視画像のそれぞれにおける被写体の移動距離及び移動方向を推定し、前記複数の透視画像のそれぞれに対して、前記推定した移動方向と逆方向へ、前記推定した移動距離だけ前記被写体を移動させるように当該透視画像を補正する補正部を備える、画像処理装置。
6. 複数の物体のそれぞれに起因する輝度の成分の重畳である画素の輝度を含む、複数の透視画像における、前記画素の輝度が複数の値のうちの各値を有する透視画像の出現頻度に基づいて、当該画素の輝度の確率分布を取得し、
   前記取得された確率分布を表す混合分布を推定し、
   前記推定された混合分布を構成する複数の部分分布のうちの期待値が最小である期待値最小分布に対応し且つ前記画素の輝度の少なくとも一部を構成する成分を推定し、
   前記透視画像に含まれる前記画素の輝度から、前記推定された成分を除去する、画像処理方法。
7. 複数の物体のそれぞれに起因する輝度の成分の重畳である画素の輝度を含む、複数の透視画像における、前記画素の輝度が複数の値のうちの各値を有する透視画像の出現頻度に基づいて、当該画素の輝度の確率分布を取得し、
   前記取得された確率分布を表す混合分布を推定し、
   前記推定された混合分布を構成する複数の部分分布のうちの期待値が最小である期待値最小分布に対応し且つ前記画素の輝度の少なくとも一部を構成する成分を推定し、
   前記透視画像に含まれる前記画素の輝度から、前記推定された成分を除去する、
   処理をコンピュータに実行させる画像処理プログラム。