JP3903027B2 - 放射線画像処理方法及び装置並びにグリッドの選別方法及び装置 - Google Patents

Description

本発明は散乱線除去用グリッドを有する放射線画像撮影装置によって得られた放射線画像の処理に関する。特に、医療用Ｘ線画像の処理に好適な画像処理技術に関する。

放射線撮影画像とは、物体もしくは人体を被写体として、被写体を透過した放射線の空間強度分布すなわち放射線に対する透過率分布を画像化したものである。このような放射線撮影画像を観察する場合に問題となる項目に、被写体などから発せられる散乱線の存在があげられる。一般に、放射線撮影画像は、放射線源と受像面の間に被写体を介在させた状態で得られる受像面の任意位置の放射線強度を、その位置と放射線源を結ぶ直線上の放射線の減衰率とみなして用いて作成される。すなわち、その直線を包含する被写体部分の放射線透過率が得られるため、被写体の局所的な内部情報が平面的に得られる。これを直接線成分と称する。

しかしながら、受像面に到達する放射線には、直接線成分以外にも被写体そのものから２次的に発せられる散乱線成分が存在し、放射線強度はこれら直接線成分と散乱線成分が加算されたものとなる。一般に放射線の強度は、ランダムに到達するエネルギー粒子の数で表される。従って、受像面で得られる放射線強度の分布は必然的にポアソン分布になり、そのばらつき（分散値）はエネルギー粒子の数と等しくなることが知られている。このようなばらつきは受像面で得られるランダムノイズとなるが、これは、散乱線の加算により被写体情報に無関係なノイズが受像面で得られる画像に加算されていることと等価になる。散乱線の加算がランダムノイズの加算と等価であるので、基本的にノイズに埋もれた微小な被写体のコントラスト情報は完全に無くなる場合もある。ここで、散乱線の加算された画像から散乱線成分のみを計算（画像処理）で除去することは非常に困難であり、現状では不可能に近いと考える。

Ｘ線発見当初から散乱線を有効に除去する手段のひとつとして、散乱線除去グリッドの使用がある。図１０は散乱線除去グリッドの使用方法を説明する図であり、５０１は受像面、５０２は散乱線除去グリッド（以下、グリッドという）、５０３は被写体である。グリッド５０２は、Ｘ線源５０４（Ｘ線焦点）へ向けて並べられた鉛製の複数の板で構成される。鉛製の板はＸ線の遮断率が高いため、Ｘ線焦点から到来する直接線成分は破線群５０５で示すように、板と板の間を通り抜けて受像面に到達するが、それ以外の方向性をもつ散乱線成分は、矢印５０６で示すように板で遮断されて受像面に到達しないという仕組みである。この方法により、かなりの量の散乱線成分を除去でき、一部の直接線成分は犠牲になるが、散乱線の影響の少ないコントラストの高い画像を得ることが可能である。

このようなグリッド５０２の弊害は、受像面５０１上に板の陰影である縞状の模様（以下グリッド縞）が現れることである。画像上のグリッド縞は、画像の観察において邪魔となる場合もある。グリッド縞を見えなくする有効な手段としては、Ｘ線曝射中にグリッドをグリッド縞を横切る方向に移動もしくは振動させることである。この方法によれば、グリッド縞のみがぼかされ、被写体の陰影のみが受像面に現れる。しかし、この方法にはグリッドそのものをＸ線曝射中に移動させる機構と、その移動速度を適切におこなうための速度制御、振動対策、移動機構が発生するノイズ対策など、コスト高になるという弊害がある。さらに、撮影するタイミングもグリッド機構の移動に合わせなければならないため自由度が少なくなるという問題点もある。

一方、グリッド縞は画像としては、被写体陰影に乗算されているため、対数変換により加算情報に変換できる。したがって、グリッド縞を含む被写体画像をデジタル画像として取得できれば、画像処理でグリッド縞情報を除去することも可能である。

一般にグリッドの構造は、一方向に鉛板を非常に正確な間隔で並べられたものになっており、画像上のグリッド縞の陰影は、正確な空間周波数ピークを持つ空間的に振動する縞模様になる。このような縞情報を除去する画像処理方法としてはたとえば特許文献１にあるように、ウエーブレット変換を用いて特定の空間周波数領域の応答を低下させる方法や特許文献２にあるようなグリッド縞の特徴に対応したグリッド縞成分のみを予測して作成し、放射線画像より減算させる方法があった。
特開２００１−１８９８６６ 特開２００２−３３０３４４ 特開２００３−３８４８１

近年、たとえば入射するＸ線粒子が半導体内に発生する自由電子を直接捕らえるいわゆる直接型と呼ばれるＸ線センサも実用化されつつあり、Ｘ線像の受像系（以下Ｘ線センサ）の空間周波数応答（以下ＭＴＦ；Modulation Transfer Function）は向上している。また、Ｘ線強度を一旦エネルギの低い蛍光に変換し、その蛍光を捕らえる間接型においても、その空間周波数応答は向上しつつある。サンプリング定理によれば、空間サンプリングピッチの半分であるナイキスト周波数以上の成分はすべてナイキスト周波数以下に折り返して観察されるため、必要以上のＭＴＦの向上は無意味にも思える。しかしながら、ナイキスト周波数以下のＭＴＦはできるだけ向上させたほうが良いため、被写体にナイキスト周波数以上の成分がほとんど無いものであるとの仮定のもとで、Ｘ線センサのＭＴＦを向上させているのである。

図１１（ａ）は、グリッド縞の陰影をあらわした図で、図示するように透過部分と遮断部分が正確に間隔Tgでもって周期的な陰影を作っている。図１１（ｂ）は、（ａ）の空間スペクトルを模式的に示したものである。グリッド縞の陰影は周期関数であるので、基本波である1／Tgの線スペクトルピークに加え、そのｎ倍（ｎ＝２，３，…）の線スペクトルピーク（ｎ倍高調波）が現れている。上述のようにＸ線センサのＭＴＦが向上して来ると、このｎ倍の線スペクトルピークも減少することなく伝わってしまうため、それらが画像上に現れることになる。従来のグリッド縞を低減する画像処理の開示では、このｎ倍の線スペクトルに相当する縞成分を除去することはほとんど考慮されていない。唯一、特許文献３において、グリッド縞の基本波と２倍高調波成分がサンプリングされた画像信号上でほぼ同一周波数になるようにグリッド縞の周波数を設定する方法が開示されているといった程度である。

本発明は上記の課題に鑑みてなされたものであり、グリッド縞の基本波成分に加えて高調波成分を放射線画像より除去可能とし、ＭＴＦの高いＸ線センサにおいてもグリッド縞の影響の少ない放射線画像を取得することを目的とする。

上記の目的を達成するための本発明による放射線画像処理装置の放射線画像処理方法は、
放射線画像からグリッド縞を構成する基本波周波数および少なくとも一つの高調波周波数を取得する取得工程と、
前記取得工程において取得された基本波周波数および少なくとも一つの高調波周波数に基づいて、前記グリッド縞を横切る第１の方向の前記放射線画像の画像信号成分のパワーと、該グリッド縞に沿った第２の方向の前記放射線画像の画像信号成分のパワーとが等しくなるようにグリッド縞予測データを生成する生成工程と、
前記生成工程において生成されたグリッド縞予測データに基づいて、前記放射線画像よりグリッド縞を除去する除去工程とを備える。

また、本発明によれば、
グリッドのみを撮影して得られた放射線画像中のグリッド縞の基本波周波数と少なくとも一つの高調波周波数を測定する測定工程と、
前記放射線画像における、前記グリッド縞を横切る第１の方向と該グリッド縞に沿った第２の方向のそれぞれについてパワースペクトルを求め、前記グリッド縞の基本波周波数部分と高調波周波数部分のパワースペクトルを除いてそれぞれの和を求め、これらパワースペクトルの和の比を取得する取得工程と、
前記取得工程で取得された比に基づいて、前記グリッドの良否を判定する判定工程とを備えることを特徴とするグリッドの選別装置の選別方法が提供される。

以上の本発明の構成によれば、グリッド縞の基本波成分に加えて高調波成分を放射線画像から除去することが可能となる。このため、ＭＴＦの高いＸ線センサにおいてもグリッド縞の影響の少ない放射線画像を取得することができる。

以下、添付の図面を参照して本発明の実施形態を説明する。

＜実施形態の概要＞
まず、各実施形態の詳細な説明を行なう前に、各実施形態の概要を説明する。上述した２倍高調波以上のグリッド縞の高調波成分にも対応するために、実施形態では、２倍高調波以上の成分をも考慮したモデルを採用してグリッド縞をモデル化し、該モデルのパラメータ調整によってグリッド縞を予測する。そして、予測されたグリッド縞を用いてグリッド縞成分を除去する。

また、上記モデルを使って適切なグリッド縞モデルのパラメータを決定する方法としては、グリッド縞を適切に除去された後の画像成分のグリッド縞を横切る方向の分散値とグリッド縞に沿った方向の分散値が概略等しくなることを基準として行う。これは、通常の被写体画像の画像信号成分の局所的なパワーには方向性がなく等方性であるという基本概念を利用している。なお、分散（２乗和の平均）とパワーとは同じ意味であり、数学的にもパワースペクトルの総和は信号成分の分散値に一致することは証明されている。

その基本的な過程は以下の（１）〜（７）の操作による。
（１）画像上のグリッド縞基本波成分である第１の線スペクトルピークFh1を正確に測定する。
（２）第１の線スペクトルピークの周波数とＸ線センサのサンプリングピッチから、２倍高調波成分周波数である第２の線スペクトルピーク周波数Fh2を計算する。
（３）上記第１の線スペクトルピークの周波数と第２の線スペクトルピークの周波数を抽出するフィルタを形成する。
（４）画像情報からグリッド縞を横切る方向のライン情報を取り出し、上記フィルタを通して、第１の線スペクトルピークと第２の線スペクトルピーク両者を含む信号を粗抽出する。
（５）粗抽出された信号成分を複数の小ブロックに分割し、それぞれのブロックのデータを第１のグリッド縞成分と第２のグリッド縞成分を兼ね備えるグリッド縞信号モデルをフィッティングする。このとき、フィッティング誤差が、当該ブロックに対応するグリッド縞に沿った方向の画像情報の分散値とほぼ等価になるようにパラメータ調整を行う。
（６）上記フィッティングで得られた信号情報がグリッド縞情報であり、元の画像情報から減算することによりグリッド縞が適切に除去された画像情報が得られる。
（７）（４）〜（６）を画像の各ラインごとに繰り返す。

ここで特に、２倍高調波の周波数Fh2の発生過程について図１２を参照して説明する。図１２（ａ）は、横軸を空間周波数とし、各種スペクトルピークの位置を表している。すなわち、サンプリングキャリアの空間周波数Fs（-Fs）および、グリッドの空間周波数Fg（-Fg）および２倍高調波の空間周波数周波数2Fg（-2Fg）である。このグリッド信号がサンプリングキャリアによって変調を受ける（乗算）のがサンプリングであると考えられる。よって、図１２（ｂ）のようにグリッド縞の周波数は、両者のコンボリューションにより変化し、サンプリング周波数の半分であるナイキスト周波数Fn（＝Fs／2）以下になり観察される（サンプリング定理）。

ここで、Fh1、Fh2の具体的な計算方法を計算の流れ図である図１３を用いて説明する。図１３では、入力としてグリッド縞が存在する画像情報が与えられる。この画像に対してフーリエ変換を用いてスペクトルピークを探索することにより、画像上の基本グリッド周波数Fh1が求まる（図１３の（ａ））。次に、このFh1を生じたグリッドの基本周波数Fgを計算しなければならない（図１３の（ｂ））。この場合グリッドの空間周波数Fgはあらかじめ大体どの位置に存在するのかを知っておく必要があるがこれは容易である。即ち、通常使用するグリッドには定格があり、その定格を知っておけば十分である。

なお、使用するグリッドの定格（製造された鉛板の間隔）がわかっていたとしても、通常グリッドとＸ線センサの間には空隙があり、Ｘ線撮影によってグリッド縞は拡大されるので、Ｘ線センサ上のグリッド縞の陰影は正確には定格どおりにはならない。また製造誤差も存在する。したがって、図１３の（ａ）の処理ではグリッドの縞周波数を測定した。しかし、おおよそのグリッド縞の周波数範囲は規定できる。本例の場合は、Fn＜Fg＜Fsの位置に存在することはあらかじめわかっているとする。 上記（１）で画像情報を解析することで測定されるのはFh1であり、これはもっとも強いグリッドの基本空間周波数によるものである。Fh1はグリッドの基本周波数信号FgをサンプリングキャリアFsで変調したものであるので下の数式で計算できる。
|Fh1|＝|K・Fs−Fg|；K=1,2,3…、|Fh1|≦Fs/2となるKを選択 …（式１）。

実施形態においては、先にFh1が画像上から測定されているので、Fn＜Fg＜Fsの関係を用いてFgを逆算しなければならない。Fgの計算式は以下になる。
Fg＝|K・Fs−Fh1|；K=1,2,3…、Fn＜Fg＜FsとなるKを選択 …（式２）。

このようにして求められたグリッド縞周波数Fgにおいて、２倍高調波の周波数は正確に「2Fg」になることは明白である（図１３（ｃ））。この2Fgの空間周波数ピークがサンプリングによってどの空間周波数として観察されるかの計算は下式で計算できる（図１３（ｄ））。
|Fh2|＝|K・Fs−2Fg|；K=1,2,3…、|Fh2|≦Fs/2となるKを選択 …（式３）。

以下の第１実施形態では、こうして、算出されたFh1、Fh2を用いてフィルタを形成し、形成されたフィルタを通過させて得られた画像情報を、グリッド縞の基本空間周波数Fh1とその２倍高調波周波数Fh2を考慮したグリッド縞信号モデルでフィッティングし、フィッティングで得られた信号情報をグリッド縞情報として元の画像情報から減算する。こうしてグリッド縞が適切に除去された画像情報が得られる。

上述した基本的な過程における（５）のフィッティング処理は、グリッド縞のパワー（分散値）が等方性を有すると仮定の下で行われるものである。第２実施形態では、グリッド自体の特性（グリッドに沿った方向と横切る方向のパワー（分散）の非等方性）について考慮することにより、更に精度の高いグリッド縞除去を行なう。すなわち、上述の基本的過程において、被写体画像情報のパワーは基本的に方向性がなく等方性であるという概念を用いているが、グリッド自体の陰影には必ずしも等方性が無い場合があるので、この非等方性に着目したグリッド縞除去を行なう。

グリッドは、前述のように鉛板を等間隔にならべたものであり、それを適切に固定するために中間物質としてアルミニウム、樹脂、木材、紙などのＸ線透過性の高いものを間に入れた構造を有している。この場合、鉛板の加工精度、中間物質の素材・加工精度、組み立て精度により、得られるグリッド陰影にはグリッド縞以外の信号成分が存在することになる。しかも、グリッドが一方向の鉛板で構成されているため、その陰影の画像信号成分には等方性がなくなってしまうおそれがある。

図１４は、グリッドをＸ線で撮影した場合の画像信号のパワースペクトルを模式的に示した図である。図１４（ａ）で示すライン５１０はグリッド縞を横切る方向の画像信号のパワースペクトルであり、グリッド基本周波数および２倍高調波縞のピークを示す２本の線スペクトルを包含している。ライン５１１はグリッド縞に沿った方向の画像信号のパワースペクトルであり、ここにはグリッド縞は存在しない。図１４の（ａ）はグリッド縞を横切る方向および沿った方向のパワースペクトルがグリッド縞に関係するところ以外はほぼ同等になっている例を示している。

一方、図１４の（ｂ）のライン５１２は別のグリッドをＸ線で撮影した場合の画像信号における、グリッド縞を横切る方向のパワースペクトルである。図１４（ｂ）の場合、グリッドの製造精度および中間物質の材質などの影響で、本来のグリッド縞には無いはずのスペクトルピーク５１４が発生している。また、５１３はグリッド縞に沿った方向の画像信号のパワースペクトルを示す。図１４（ｂ）においては、グリッドの縞に沿った方向とグリッド縞を横切る方向のパワースペクトルに大きな差が現れている。

第２実施形態では、グリッド陰影のノイズ特性の等方性を表すファクターとして、「グリッド品質値；Grid Quality Value;ＧＱ」をあらかじめ求めておき、この値を用いてグリッド縞を除去するときのパラメータとする。

図１５はＧＱを求める方法を図示したものである。図１５において（ａ）のグラフはグリッド縞に沿った方向のパワースペクトルであり、（ｂ）はグリッド縞を横切る方向のパワースペクトルである。グリッドを横切る方向のパワースペクトルには前述のグリッド基本周波数ピークFh1と２倍高調波ピークFh2が存在する。ここで、ＧＱ値をグリッド縞を横切る方向の分散値とグリッド縞に沿った方向の分散値の比によって定義する。この場合グリッド縞に起因する成分による差およびＸ線撮影時のＸ線放射分布形状の差を無くすために、低空間周波数領域（0.5cyc/mm以下）および、Fh1の周辺およびFh2の周辺の成分を無視して、それ以外の部分のパワースペクトルの加算によってそれぞれの方向の分散値とする。図１５ではそれぞれVarVとVarHが加算値として得られ、
ＧＱ＝VarH／VarV …（式４）
としてＧＱ値が求まる。

したがって、以下に説明する第２実施形態では、上で説明した（１）〜（７）の手続きにおいて、特に（５）の手続きを、
（５）粗抽出された信号成分を複数の小ブロックに分割し、それぞれのブロックのデータを第１のグリッド縞成分と第２のグリッド縞成分を兼ね備え持つグリッド縞信号モデルでフィッティングする。このとき、フィッティング誤差が、当該ブロックに対応するグリッド縞に沿った方向の画像情報の分散値のＧＱ倍になるようにパラメータ調整を行う、ものとする。

また、グリッドモデルのフィッティング方法に関しては、データ点数を調整する方法（第１〜第２実施形態）と複数のモデルを用意しておいて適切なものを選ぶ方法が挙げられる（第３実施形態）。さらに、第４実施形態では、ＧＱ値を積極的に用いて、グリッド自体の品質を定義し、品質の悪いグリッド（ＧＱ値が１から遠ざかるもの）をあらかじめ使わない用にグリッドを選別することを提案する。

以下、各実施形態について具体的に説明する。なお、以下の各実施形態では、２倍高調波のみを考慮する場合について述べているが、それ以上の第３、第４…第ｎ倍の高調波についても一般的に行うことが可能であることは当業者には明らかであろう。

＜第１実施形態＞
図１は第１実施形態によるＸ線撮影装置のブロック図である。図１において、１は散乱線除去グリッド（以下、グリッド）であり、２は複数の画素がマトリックス状に配置されたＸ線画像センサである。なお、Ｘ線画像センサ２の画素マトリックスの画素間隔をTとし、サンプリング空間周波数Fs＝1／Tであるとする。３はＸ線発生装置であり、４は被写体であり、本例では人体である。Ｘ線発生装置３は不図示の高電圧発生装置を含むＸ線発生制御装置によって制御され、操作者の操作によりＸ線を被写体４、グリッド１およびＸ線画像センサ２へ向かってＸ線を発する。Ｘ線画像センサ２はＸ線発生に同期して駆動され、Ｘ線放射されている間、被写体を透過したＸ線量に応じた電荷を各画素に蓄積する。各画素の蓄積電荷値はデジタル値に変換され、画像データとして出力される。

５は画像取得補正部であり、Ｘ線画像センサ２から出力された画像データを取得し、各種補正を行う。ここでいう各種補正とは、Ｘ線画像センサから出力される画像データの暗出力値の補正および各画素のゲインばらつきの補正および欠陥画素の補正を意味する。暗出力値の補正は、Ｘ線を放射しない状態の画像データをあらかじめもしくはＸ線画像取得直後に取得し、被写体を撮影して得られたＸ線画像から減算するものである。また、ゲインばらつきの補正は、被写体が無い状態でＸ線放射を行った場合の画像データを用いて被写体を撮影したＸ線画像を除算するものである。更に、欠陥画素の補正は、製造段階であらかじめ欠陥画素位置を把握しておき、当該画素値を周辺の画素値から平均計算などで予測して補正するものである。

画像取得補正部５で補正された画像データはグリッド縞除去部６に提供され、グリッド縞除去の処理がなされる。続いて、ＱＡ処理部７にて、診断しやすいようにＱＡ処理（Quality Assurance 処理）がなされた後、外部記憶装置もしくは外部の画像印刷装置などへ出力される。ここでいうＱＡ処理とは、階調変換処理、周波数強調処理、切り出し、縮小・拡大処理などの一般的な画像処理を表す。

以上説明した画像取得補正部５、グリッド縞除去部６、ＱＡ処理部７は、例えばメモリに格納されたプログラムを実行して処理の処理を実行するコンピュータ装置によって実現される。以下、本実施形態の特徴的な構成であるところのグリッド縞除去部６について詳細に説明する。

図２はグリッド縞除去部６の内部構成をさらに詳細なブロックに分割して示した図である。グリッド縞除去部６において、グリッド縞推定部８は入力画像からグリッド縞を推定し、これをグリッド画像メモリ９へ記憶する。続いて、減算器１０により入力画像からグリッド縞メモリ９の内容が減算されることにより、グリッド縞を除去した画像が出力される。

続いて図３に示すフローチャートを用いて図２の処理を詳細に説明する。まず、ステップＳ１０１において入力画像に対してグリッドの存在および方向の判定を行う。グリッドの存在や方向などは、機械的な機構をもちいて判定することも可能であるが、本例では画像データを解析することによりグリッドの存在および方向の判定を行う。具体的な方法としては、画像から横方向のラインおよび縦方向のラインを任意に選択し、高速フーリエ変換（ＦＦＴ）を用いて選択したラインの周波数特性を得る。得られた周波数特性において、強いピークの存在するラインの方向にグリッド縞が存在すると判断する。縦横両者にそのような強いピークが存在しない場合にはグリッドが存在しないと判定し、本処理を終了させる（不図示）。

ステップＳ１０２では、ステップＳ１０１で得られた結果より、グリッド縞の存在するライン方向についてさらに詳細に（点数を増やして）フーリエ変換を行い、画像上のグリッド縞周波数Fh1を測定する。次にステップＳ１０３において、２倍高調波成分の周波数Fh2を計算する。２倍高調波成分の周波数Fh2の計算方法は、前述のように、式２によってグリッド縞の周波数Fgを求めた後、それを２倍して、式３により求めることができる。次いで、ステップＳ１０４において、グリッド縞信号を粗抽出するためのＦＩＲ（Finite Impulse Response ）フィルタを形成する。このときＦＩＲフィルタによって信号成分の位相が変化しないように、ＦＩＲフィルタは左右対称な係数を持つものを選択する。図４は本実施形態で用いる５点ＦＩＲフィルタの係数a₀〜a₂の配置を示したものである。このフィルタの周波数特性G(f)は、Ｘ線センサのサンプリングピッチをTとすると以下の式５のようになる。

（式５）中のパラメータa₀〜a₂は以下のように決定される。すなわち、
G(0)=0 ; 画像の直流成分および画像の低周波成分を排除、
G(Fh1) = 1; 画像上のグリッド縞成分を正確に抽出、
G(Fh2)=2; 画像上のグリッド縞2倍高調波成分を正確に抽出、
の３条件を（式５）に代入して連立方程式を解くことにより下記の（式６）が得られる。

（式６）により（式５）に示したＦＩＲフィルタの係数を求め、ＦＩＲフィルタを形成する。

以上のようにしてＦＩＲフィルタを形成したならば、ステップＳ１０５〜Ｓ１１３により、グリッド縞除去の処理を行なう。まず、ステップＳ１０５では、すべての画像情報を処理したかどうかを判断する。未処理のラインがあればステップＳ１０６へ進み、グリッド縞を横切る方向の注目ラインのライン情報を抽出する。本実施形態では、このライン情報はグリッド縞に沿った方向に並ぶｍラインの画素列の平均値で構成される。また同時にｍラインの各画素列の分散値を計算する。すなわち、上下ｍラインの平均画素値とその平均に供した画素のグリッド縞に沿った方向の分散値がライン中の画素数分だけ得られ、平均画素値列によりライン情報が形成される。

図１６はこの様子を説明する図である。グリッド縞を横切る方向（Ｘ方向）に沿ったｍ本のライン（注目ラインを真中に挟む）の、グリッド縞に沿った方向（Ｙ方向）に並ぶｍ個の画素からなる画素列について画素値の平均をとり（式（Ａ））、これを注目ラインの画素値としてライン情報を生成する。また、この画素列に関して式（Ｂ）により分散を求めておく。

次に、ステップＳ１０７において、ステップＳ１０６で得られたライン情報を上述の（式６）で求めたＦＩＲフィルタに通すことにより、グリッド縞情報の粗抽出を行う。そして、ステップＳ１０８において、粗抽出されたグリッド縞情報の中で、画像エッジなどの影響で、後段のフィッティング演算に提供することができない部分（以下不良部分）を判定する。具体的には、粗抽出された信号の振幅値が異常に大きいか、位相変動が異常に大きいかを適当な閾値を用いて判断し、大きい場合にはその位置の情報を以下の演算で用いないようにする。なお、振幅値、位相変動の大きさは、例えば位相が９０度移動するフィルタ（微分フィルタ）を用い、元の信号Ａと微分フィルタの出力信号Ｂと２つの信号により振幅値と位相値の変動を見る。より具体的には、元の信号Ａと微分フィルタの出力信号Ｂの2乗和のルートをとることにより振幅値の変動がわかる。またTan^-1(B/A)の変動を見ることで位相値が得られるため、その変動をチェックすることで位相変動を見ることができる（さらに正確にはTan^-1演算の出力は演算の性格上一般に−π〜＋πの出力値域を持つので、それらをベクトルの回転を考慮した通常の値（範囲を絞らない）に変換、具体的にはつながった値にする必要がある（phase unwrappingと呼ばれる処理である））。

次に、ステップＳ１０９では、上記の不良部分を除く部分をさらに小さなブロック（通常１０画素〜３０画素程度）に分割する。本実施形態では、図１６のＢに示されるように、１つのブロックにｎ個の画素が含まれるものとする。グリッド縞自体は全体で安定したものであるが、被写体の散乱線成分などの影響を受けて、部分的には定常な信号であるが、大きく見ると非定常な信号になっている場合があり、非定常な信号をフィッティングするのを防ぐために、定常と見なせる小部分に分割する。なお、本実施形態では不良部分を除去してから小ブロックに分割する（無駄を省くため）が、小ブロックに分割してから、不良部分を含むブロックを除去するようにしてもよい。次いで、ステップＳ１１０では、得られた各ブロックについてグリッド情報を取り出すためのフィッティングを行う。

本実施形態で用いる、グリッドの数式モデルG1(x)を下式で示す。このモデルによれば、グリッド縞の測定された基本周波数と高調波周波数が考慮される。

このモデルは、n点の入力データ系列｛xi , yi；i=1,2…,n｝が与えられた場合に以下の（式８）の連立方程式を解くことで各パラメータc₁，c₂，d₁，d₂を求めることができる。

ここでω₁=2πFh1、ω₂=2πFh2、Σはｎ個の加算を意味する。

ステップＳ１１０では、上記（式８）を用いて各ブロックのフィッティング演算を行うが、各ブロックのフィッティング演算の詳細について図５のフローチャート及び図１６を用いて説明する。
まず、ステップＳ２０１において、当該ブロック内における、グリッド縞に沿った方向に並ぶ画素列の分散の平均値（以下、平均分散値Vhという）を計算する。即ち、図１６において、v₁〜v_nの平均を求め（式（Ｃ））、これを平均分散値Vhとする。本実施形態のフィッティング演算ではこの平均分散値Vhが基本となり、フィッティング誤差がこの値に近づくような演算を行う。ステップＳ２０２において、ブロック内のデータを用いて（式８）を用いてフィッティング演算を行う。次いで、ステップＳ２０３において、ステップＳ２０３でフィッティングされた関数と元データとの間でフィッティング誤差（２乗誤差の平均）Vsを求める（図１６の式（Ｄ））。式（Ｄ）から明らかなように、ここで求まるVsは、当該注目ラインのデータからグリッド縞を除去した後の分散に等しくなる。よって、このVsがVhに近づくようにフィッティングを調整することにより、グリッド縞を横切る方向とこれに沿った方向の分散が等しくなるようにグリッド縞が除去されることがわかる。

ステップＳ２０４では、フィッティング誤差VsがステップＳ２０１で算出された平均分散値Vhを超えたかどうかを判断する。フィッティング誤差Vsが平均分散値Vhを超えていない場合、処理はステップＳ２０５へ進み、フィッティングの状態を変化させる。本実施形態では、元のデータ（当該ブロック内のライン情報）とフィッティングされた関数との間で誤差の最も大きい画素値をフィッティング演算から除外する。これは、（式７）のグリッド縞モデルから外れた画素値は、グリッド縞成分以外に被写体画像やノイズの影響を多く受けていると思われるためである。すなわち、そのようなグリッド縞成分以外の成分を多く含むと想定される画素をグリッド縞のフィッティングから外してフィッティング演算を行なうことにより、より純粋なグリッド縞情報を得ようとするものである。

上記の画素値除外により、ブロック内のフィッティング誤差Vsはより大きくなる。ステップＳ２０６では、新たなフィッティング演算を行い、ステップＳ２０７でフィッティング誤差Vsを計算しなおす。そしてステップＳ２０４へ戻り、再び平均分散値Vhとフィッティング誤差Vsを比較する。以上の処理を繰り返し、ステップＳ２０４においてVsがVhを超えた時点で、ステップＳ２０８へ進み、それまでの演算のなかでもっともVhに近いフィッティング誤差Vsを持つフィッティングパラメータを選択し、当該ブロックのグリッド縞予測データとする。

以上のようにして、図３のステップＳ１１０では各ブロックについて図５の演算を行い、各ブロックのグリッド縞予測データを作成する。続いて、ステップＳ１１１では、ステップＳ１０８で排除した、フィッティングに提供されなかった不良画素部分のグリッド縞予測データを隣接するフィッティングに提供されたブロックのフィッティングパラメータによって外挿計算する。こうしてライン全体のグリッド縞予測データを完成させる。ステップＳ１１２では完成されたグリッド縞予測データを図２のグリッド画像メモリ９の該当するラインに記憶する。こうして入力画像の全てのラインについてグリッド縞予測データが得られたならば、ステップＳ１０５からステップＳ１１３へ進み、減算器１０により入力画像からグリッド画像を減算し、グリッド縞除去が終了する。

なお、（式５）〜（式８）においては、２倍高調波のみに着目しているが、３倍、４倍などさらに高次の高調波についても同様の演算が可能であり、必要であれば同様の対応が可能であるのは、上記開示から当業者には明らかであろう。

＜第２実施形態＞
第２実施形態では、上述したＧＱ値を用いて、グリッドの品質にあわせたグリッド縞除去を行なう。まず、ＧＱ値を算出するための構成及び処理について説明する。

図６は、ＧＱ値を算出するための構成を示すブロック図であり、図１と比較して被写体４が存在せず、新たにＧＱ値計算部１１が追加されている。この装置では、あらかじめグリッドを被写体なしで撮影し、ＧＱ値計算部１１によりその陰影を解析することによりＧＱ値を求めておく。その後、図１で示した構成を用いてグリッド縞除去を行う。第２実施形態ではグリッド縞除去部６においてＧＱ値を用いる点が第１実施形態と異なる点である。

図６において図１に示した構成と同じ構成には同一の符号を付してある。被写体なしで、Ｘ線源３とグリッド１及びＸ線画像センサ２と適切な距離、位置を保った状態で、グリッド１の陰影を撮影する。そして、ＧＱ値計算部１１によって当該グリッド１のＧＱ値が算出される。以下、図７のフローチャートに従ってＧＱ値計算部１１の動作を詳細に説明する。なお、ＧＱ値の原理的な説明は「実施形態の概要」で上述したとおりである。

図７において、ステップＳ１０１〜Ｓ１０４の処理は第１実施形態（図３）で説明したとおりである。ステップＳ３０２〜Ｓ３０７においてグリッド縞を横切る方向のパワースペクトルが計算される。即ち、あらかじめ定められた部分のｎラインのデータに対してステップＳ１０４で形成したＦＩＲフィルタを用いてフィルタリングを施し（ステップＳ３０２、Ｓ３０３）、フィルタリング後のデータのパワースペクトルを計算し（ステップＳ３０４）、その平均を求める（ステップＳ３０１、Ｓ３０５）。そして、図１５を参照して上述したように、低周波部分とFh1付近、Fh2付近を除去し（ステップＳ３０６）、残ったパワースペクトルの総和（VarH）を求める（ステップＳ３０７）。なお、上記のように低周波部分とFh1付近、Fh2付近を適切に除去するため、パワースペクトルを計算するにあたってフーリエ変換の点数は１０２４点以上であることがのぞましい。

続いて、ステップＳ３０８〜Ｓ３１４により、グリッド縞に沿った方向のパワースペクトルを計算する。まず、あらかじめ定められた部分のグリッド縞に沿った方向のｍラインについて、ＦＩＲフィルタリングを施し（ステップＳ３０９、Ｓ３１０）、パワースペクトルを計算し（ステップＳ３１１）、その平均を求める（ステップＳ３０８、Ｓ３１２）。そして、図１５を参照して上述したように、低周波部分とFh1付近、Fh2付近を除去し（ステップＳ３１３）、残ったパワースペクトルの総和（VarV）を求める（ステップＳ３１４）。なお、上記のように低周波部分とFh1付近、Fh2付近を適切に除去するため、パワースペクトルを計算するにあたってフーリエ変換の点数は１０２４点以上であることがのぞましい。そして、ステップＳ３１５でVarH／VarVを計算し、これをＧＱ値とする。

なお、上記の計算方法において、ＦＩＲフィルタを用いた理由は、低周波成分をより効率よく除去するためであり、ＦＩＲフィルタは必須ではない。また、上記処理ではラインごとにパワースペクトル計算を行って平均処理を行っているが、平均処理後にパワースペクトル演算を行っても同様の結果が得られる。以上のようにして得られたＧＱ値は、グリッド縞を除去された後の画像の品質を表し、フィッティングの目標となる。なお、“画像の品質”と“画像信号成分のパワーの比”とは密接に関連する。通常の人体である被写体を撮影した場合の画像信号の分散（パワー）の縦横の比は局所的にみればほぼ同じになるはずであるが、それが大きく異なるとすると、何らかの強い方向性をもつ画像が被写体に重なっていると考えられる。本実施形態の場合、その重なっているものがグリッドであり、基本の縞成分を除去しても残留しているとすれば、なんらかの不要成分がグリッドに存在していることになる。この不要成分は特性が不明であるため除去できず残留しつづけることになるため、結果的に画像の品質が悪くなるのである。

次に、第２実施形態によるグリッド縞除去処理について説明する。第２実施形態によるグリッド縞除去処理の全体の流れは第１実施形態（図３）と同様である。第１実施形態と異なるのは、ステップＳ１１０におけるフィッティング処理であり、以下この処理について図８を用いて説明する。

図８は、図５と同様ではあるが、ＧＱ値を用いたグリッド縞のフィッティング方法を示している。なお、図５に示した処理と同じ処理には同一のステップ番号を付してある。ステップＳ２０１で、ブロック内のグリッド縞に沿った方向の平均分散値Vhが計算されるが、分散値はグリッドの品質によりＧＱ倍大きくなるはずである。従って第２実施形態では、グリッド縞を横切る方向のフィッティング誤差の目標として平均分散値Vhを用いる代わりに、このVhにＧＱ値を乗じた値Vx（＝Vh×ＧＱ）を用いる（ステップＳ４０１、Ｓ４０２）。そして、上述のステップＳ２０５〜Ｓ２０７を繰り返すことにより、フィッティング誤差VsがVxを越えた時点でステップＳ４０２からステップＳ４０３へ進み、フィッティング誤差VsがVxに最も近いものを採用してグリッド縞の予測を行う。

以上のように第２実施形態によれば、グリッドの特性を考慮したグリッド縞の予測が行なわれるので、より高精度にグリッド縞を除去することが可能となる。

なお、上記第２実施形態では、画像全体に対して一つのＧＱ値を決定し、目標値Vxを設定したが、画像を複数の部分に分け、各部分についてＧＱ値を求めることにより画像上のＧＱ値の分布を作成し、ブロックが属する部分に応じてＧＱ値を切り換えてグリッド縞を予測するようにしてもよい。

＜第３実施形態＞
上記第１及び第２実施形態では、適切なフィッティングを行うために、使用するデータ点数を変化させたが（ステップＳ２０５）、フィッティングの状態を変更する方法はこれに限られるものではない。例えば、モデルそのものを変化させる方法を用いることもできる。例えば、上述の（式７）で示したモデル以外に、下記の（式９）〜（式１１）モデルも考慮する。

ここでω₁=2πFh1、ω₂=2πFh2である。

（式９）は２倍高調波を考慮しないモデルであり、（式１０）はブロック内での多少の非定常性があることを加味し、直線的な振幅変動成分を考慮したモデルであり、（式１１）は２倍高調波成分においても直線的な振幅変動成分を考慮したモデルである。これらモデルの各々においても、（数式８）で示したような連立方程式によって与えられたデータ系列によるフィッティングが可能である。なお、これらのモデルの他にも、モデルの類型としてさまざまなモデル構成が可能であることは言うまでもない。第３実施形態では、ｎ個のモデルについてフィッティングを行い、残留誤差Vsが目標値に最も近いモデルを、フィッティングに最も適したモデルとして選択する。なお、目標値としてはグリッドに沿った方向の平均分散値Vh（第１実施形態）あるいは平均分散値VhのＧＱ倍（第２実施形態）のいずれかを採用することができる。

図９は第３実施形態による、ブロック内におけるグリッド縞の予測をフィッティングにより行うフィッティング処理を説明するフローチャートであり、第２実施形態の図８に置き換わる処理を示した。図８で説明したように、ブロック内のグリッド縞に沿った方向の平均分散値Vhを求め（ステップＳ２０１）、これをＧＱ倍して目標値Vxを得る（ステップＳ４０１）。そして、ステップＳ５０１〜Ｓ５０３において、ｎ個のモデルG₁〜G_nを用いて、同じデータ（同じライン）に関してフィッティングを行い、それぞれの残留誤差Vs₁〜Vs_nを計算する。そして、ステップＳ５０４において、最もVxに近い値の残留誤差を有するモデルをグリッド縞として採択することにより、適切なグリッド縞予測を実現する。

なお、第３実施形態においては、第２実施形態の処理を適用したが、第１実施形態のよう仁平均分散値Vhをそのまま目標値としてもよい。この場合、ステップＳ４０１は省略され、ステップＳ５０４では、最もVhに近い値の残留誤差を有するモデルを採択すればよい。また、また、それぞれのモデルのフィッティングにおいて、図５、図８と同様にデータ点数を変化させつつ、もっともVxに近いフィッティングを行った後に、それぞれのモデルの結果比較を行うという複合された手法を用いることも可能である。

＜第４実施形態＞
第２実施形態で説明したＧＱ値は、グリッド縞除去だけではなく、グリッドそのものの品質管理に用いることが可能である。すなわち、ＧＱ値は１であることが理想であるが、グリッドの製造過程でグリッドの基本周波数以外の固定ノイズが混在する場合には、品質の悪いグリッドであるといえる。たとえば、グリッドを使用する前段階において、図６の装置および図７の計算方法でＧＱ値を求めておき、ＧＱ値が３以上もしくは０．８以下のグリッドは不良なグリッドであるとして排除するなどの処置をとることが可能である。すなわち、不良グリッドの選別にＧＱ値を適用することができる。

以上説明したように上記各実施形態によれば、グリッド縞の基本波成分に加えて高調波成分も除去することにより、ＭＴＦの高いＸ線センサにおいても適切なグリッド縞除去が可能になる。また、その際、画像の横方向と縦方向の分散値を等価もしくはグリッド品質値（ＧＱ値）倍に設定することで、より適切なグリッド縞の予測が行える。

なお、第１乃至第３実施形態において説明したグリッド縞除去、及び第４実施形態で説明したグリッド選別は、コンピュータによってＸ線画像情報を処理することにより実現可能である。従って、本発明の目的は、前述した実施形態の機能を実現するソフトウェアのプログラムコードを記録した記憶媒体を、システムあるいは装置に供給し、そのシステムあるいは装置のコンピュータ（またはＣＰＵやＭＰＵ）が記憶媒体に格納されたプログラムコードを読出し実行することによっても、達成されるものである。

この場合、記憶媒体から読出されたプログラムコード自体が前述した実施形態の機能を実現することになり、そのプログラムコードを記憶した記憶媒体は本発明を構成することになる。

プログラムコードを供給するための記憶媒体としては、例えば、フレキシブルディスク，ハードディスク，光ディスク，光磁気ディスク，ＣＤ−ＲＯＭ，ＣＤ−Ｒ，磁気テープ，不揮発性のメモリカード，ＲＯＭなどを用いることができる。

また、コンピュータが読出したプログラムコードを実行することにより、前述した実施形態の機能が実現されるだけでなく、そのプログラムコードの指示に基づき、コンピュータ上で稼働しているＯＳ（オペレーティングシステム）などが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。

さらに、記憶媒体から読出されたプログラムコードが、コンピュータに挿入された機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに書込まれた後、そのプログラムコードの指示に基づき、その機能拡張ボードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。

第１実施形態によるＸ線撮像装置の構成を示すブロック図である。 グリッド除去部の詳細な構成を示すブロック図である。 第１実施形態によるグリッド除去部の動作を示すフローチャートである。 第１実施形態で採用されるＦＩＲフィルタの係数を説明する図である。 第１実施形態によるグリッド縞の予測処理を示すフローチャートである。 第２実施形態によるＧＱ値測定のための構成を示すブロック図である。 ＧＱ値の計算処理を示すフローチャートである。 第２実施形態によるグリッド縞の予測処理を示すフローチャートである。 第３実施形態によるグリッド縞の予測処理を示すフローチャートである。 グリッドによる散乱線除去を説明する図である。 画像上におけるグリッド縞とその高調波を説明する図である。 グリッド縞とその２倍高調波がサンプリングされたときの周波数変化を説明する図である。 グリッド縞とその２倍高調波がサンプリングされたときの周波数計算の流れを説明する図である。 グリッド画像の縦・横双方のスペクトルの違いを説明する図である。 ＧＱ値の計算方法を説明する図である。 フィッティング処理時の画素平均値及び分散の計算について説明する図である。

Claims (16)

1. 放射線画像からグリッド縞を構成する基本波周波数および少なくとも一つの高調波周波数を取得する取得工程と、
   前記取得工程において取得された基本波周波数および少なくとも一つの高調波周波数に基づいて、前記グリッド縞を横切る第１の方向の前記放射線画像の画像信号成分のパワーと、該グリッド縞に沿った第２の方向の前記放射線画像の画像信号成分のパワーとが等しくなるようにグリッド縞予測データを生成する生成工程と、
   前記生成工程において生成されたグリッド縞予測データに基づいて、前記放射線画像よりグリッド縞を除去する除去工程とを備えることを特徴とする放射線画像処理装置の放射線画像処理方法。
2. 前記生成工程は、
   前記取得された基本波周波数および少なくとも一つの高調波周波数に基づいてフィルタを形成する形成工程と、
   前記形成工程において形成されたフィルタを用いて前記放射線画像からグリッド縞情報を抽出する抽出工程とを有し、
   前記抽出されたグリッド縞情報に基づいて前記グリッド縞予測データを生成することを特徴とする請求項１に記載の放射線画像処理装置の放射線画像処理方法。
3. 前記取得工程は、前記放射線画像をフーリエ変換して得られた周波数特性に基づいて前記基本波周波数と前記少なくとも一つの高調波周波数を得ることを特徴とする請求項１に記載の放射線画像処理装置の放射線画像処理方法。
4. 前記生成工程は、
   前記放射線画像の前記第１の方向に沿った注目ラインを含む複数のラインについて前記第２の方向に並ぶ画素群毎に平均値と分散を算出する算出工程と、
   前記平均値を用いてグリッド縞モデルをフィッティングした場合のフィッティング誤差を前記第１の方向の分散とし、前記算出工程で算出された分散の平均を前記第２の方向の分散とし、該第１及び第２の方向の分散が等しくなるように前記フィッティングを調整して前記グリッド縞予測データを生成する調整工程とを備えることを特徴とする請求項１に記載の放射線画像処理装置の放射線画像処理方法。
5. 前記注目ラインを更に所定画素数毎に分割し、各分割単位で前記算出工程と前記調整工程を実行することを特徴とする請求項４に記載の放射線画像処理装置の放射線画像処理方法。
6. 前記前記調整工程は、フィッティング処理に用いる平均値を間引くことによりフィッティングを調整することを特徴とする請求項４又は５に記載の放射線画像処理装置の放射線画像処理方法。
7. 前記調整工程は、予め用意された複数種類のグリッド縞モデルについてフィッティングを行ない、該第１及び第２の方向の分散が等しくなるようなグリッド縞モデルを採用することを特徴とする請求項４又は５に記載の放射線画像処理装置の放射線画像処理方法。
8. 前記生成工程は、
   前記抽出工程で得られたグリッド縞情報に基づいて、前記測定工程で測定された基本周波数と高調波周波数を考慮したグリッド縞モデルをフィッティングすることにより前記グリッド縞予測データを生成し、
   前記抽出工程で抽出されたグリッド縞情報から信号の振幅値が所定値より大きい及び／又は位相変動が所定値より大きい部分を除去し、
   前記フィッティングの結果に基づいて当該除去された部分を補間することにより、前記グリッド縞予測データを生成することを特徴とする請求項２に記載の放射線画像処理装置の放射線画像処理方法。
9. グリッドのみを撮影して得られたグリッド縞画像に基づいて前記第１及び第２の方向の画像信号成分のパワー比を取得する取得工程を更に備え、
   前記生成工程は、前記第１の方向の画像信号成分のパワーと、前記第２の方向の画像信号成分のパワーに前記パワー比を作用させた値が最も近くなるようにグリッド縞モデルのフィッティングを行ない、当該フィッティングされたグリッド縞モデルを用いてグリッド縞データを生成することを特徴とする請求項１に記載の放射線画像処理装置の放射線画像処理方法。
10. 前記取得工程は、前記グリッド縞画像における前記第１及び第２方向のそれぞれについてパワースペクトルを求め、前記グリッド縞の基本波周波数部分と高調波周波数部分のパワースペクトルを除いてそれぞれの和を求め、該第１及び第２方向のそれぞれについて求めたパワースペクトルの和の比を前記パワー比とすることを特徴とする請求項９に記載の放射線画像処理装置の放射線画像処理方法。
11. グリッドのみを撮影して得られた放射線画像中のグリッド縞の基本波周波数と少なくとも一つの高調波周波数を測定する測定工程と、
    前記放射線画像における、前記グリッド縞を横切る第１の方向と該グリッド縞に沿った第２の方向のそれぞれについてパワースペクトルを求め、前記グリッド縞の基本波周波数部分と高調波周波数部分のパワースペクトルを除いてそれぞれの和を求め、これらパワースペクトルの和の比を取得する取得工程と、
    前記取得工程で取得された比に基づいて、前記グリッドの良否を判定する判定工程とを備えることを特徴とするグリッドの選別装置の選別方法。
12. 放射線画像からグリッド縞を構成する基本波周波数および少なくとも一つの高調波周波数を取得する取得手段と、
    前記取得手段によって取得された基本波周波数および少なくとも一つの高調波周波数に基づいて、前記グリッド縞を横切る第１の方向の前記放射線画像の画像信号成分のパワーと、該グリッド縞に沿った第２の方向の前記放射線画像の画像信号成分のパワーとが等しくなるようにグリッド縞予測データを生成する生成手段と、
    前記生成手段によって生成されたグリッド縞予測データに基づいて、前記放射線画像よりグリッド縞を除去する除去手段とを備えることを特徴とする放射線画像処理装置。
13. グリッドのみを撮影して得られた放射線画像中のグリッド縞の基本波周波数と少なくとも一つの高調波周波数を測定する測定手段と、
    前記放射線画像における、前記グリッド縞を横切る第１の方向と該グリッド縞に沿った第２の方向のそれぞれについてパワースペクトルを求め、前記グリッド縞の基本波周波数部分と高調波周波数部分のパワースペクトルを除いてそれぞれの和を求め、これらパワースペクトルの和の比を取得する取得手段と、
    前記取得手段で取得された比に基づいて、前記グリッドの良否を判定する判定手段とを備えることを特徴とするグリッドの選別装置。
14. 請求項１乃至１０のいずれかに記載の放射線画像処理方法をコンピュータに実行させるための制御プログラム。
15. 請求項１１に記載のグリッド選別方法をコンピュータに実行させるための制御プログラム。
16. 請求項１４又は１５に記載の制御プログラムを格納した記憶媒体。

Images