JP4746761B2 - Radiation image processing apparatus, radiation image processing method, storage medium, and program - Google Patents

Description

【０１２９】
なる式（１）で表される。
【０１３０】
本実施の形態では、上記式（１）で表される空間周波数ｆｍに相当する縞模様がグリッド縞成分として対象画像中に存在することを考慮して、上述した第１処理ステップにおいて、グリッド縞を正確に抽出する。すなわち、グリッド縞の空間周波数ｆｇは、予め判明しているため、サンプリングされたラインデータの中から、対象画像の空間周波数ｆｍ近辺を探索し、その検索結果により、ピーク値を示す空間周波数を以って、対象画像におけるグリッド縞の空間周波数ｆｍと見なす。
【０１３１】
また、上述した第２処理ステップにおいて、空間周波数ｆｍに対し、その空間周波数ｆｍを中心として、できるだけ小さなスパンのＦＩＲフィルタリングを行うことで、有効な画像成分を殆ど除去した状態で、且つ、急峻な変動（例えばエッジ部分）によるアーチファクトの影響が狭い範囲に収まる状態で、グリッド縞成分を粗く抽出する。
【０１３２】
このとき、位相変動をなくすために、例えば、ＦＩＲフィルタの係数系列を偶関数とし、また、上記の狭い範囲を満たすために、３点若しくは５点のＦＩＲフィルタが望ましい。
【０１３３】
具体的には例えば、対称３点のＦＩＲフィルタとし、その係数を（ａ１，ｂ１，ａ１）とすると、当該係数（ａ１，ｂ１，ａ１）を求めるためには、空間周波数ｆｍにおけるレスポンスが“１”であるという条件、及び画像情報の中心値である直流成分を“０”にするという条件の２つの条件を用いることができる。
【０１３４】
すなわち、当該係数演算は、
２^*ａ１＋ｂ１＝０
２^*ａ１^*ｃｏｓ（２πｆｍＴ）＋ｂ１＝１
なる式で表される連立方程式であり、
【０１３５】
【数２】

【０１３６】
なる式（２）で表される解をとる。
【０１３７】
上記のＦＩＲフィルタリングは、空間周波数ｆｍにおいてレスポンスが“１”であるが、空間周波数がそれ以上になると、レスポンスは次第に上昇していく。一般的に、この部分には画像成分が存在しないため、当該ＦＩＲフィルタイリングであってもグリッド縞を充分に抽出できる。
【０１３８】
また、対称５点のＦＩＲフィルタリングとし、その係数を（ａ２，ｂ２，ｃ２，ｂ２，ａ２）とすると、当該係数（ａ２，ｂ２，ｃ２，ｂ２，ａ２）を求めるためには、空間周波数ｆｍにおけるレスポンスが“１”であるという条件、及び画像情報の中心値である直流成分を“０”にするという条件の２つの条件の外に、空間周波数ｆｍにおけるレスポンスの微分値が“０”（ピーク）を示すという条件をも用いることができる。
【０１３９】
すなわち、当該係数演算は、上記式（２）で表される解から、簡単な演算を行うことで、
（−ａ１²，２ａ１（１−ｂ１），１−２ａ１²−（１−ｂ１）²，
２ａ１（１−ｂ１），−ａ１²）
なる解が得られる。
【０１４０】
対称５点のＦＩＲフィルタリングのフィルタの求め方としては、例えば、先ず、上記式（２）で表される係数（ａ１，ｂ１，ａ１）を有するフィルタを、“１”から差し引いた形にすると、空間数端数ｆｍで零点を有するフィルタとなる。このフィルタによるフィルタリングを２回施す処理を考慮すると、やはり空間周波数ｆｍで零点を有するが、位相（符号）の反転がなくなる。このようなフィルタが、対称５点のＦＩＲフィルタリングのフィルタであり、“１”から当該フィルタを差し引くことで、目的とする空間周波数ｆｍにおいてピークを有するフィルタを構成できる。
【０１４１】
図１は、上述したような、対称３点のＦＩＲフィルタリング（以下、「３点ＦＩＲフィルタリング」とも言う）、及び対称５点のＦＩＲフィルタリング（以下、「５点ＦＩＲフィルタリング」とも言う）の形状（空間周波数特性）の例を示したものである。
【０１４２】
上記図１に示されるＦＩＲフィルタによるフィルタリングの結果は、殆どの場合、グリッド縞成分のみを抽出した結果となる。これは、上記図１から明かなように、主に低周波成分からなる有効な画像成分の当該低周波成分の多くが除去されるためである。
【０１４３】
しかしながら、上述したように、FIRフィルタリングにより抽出した成分の中には、かなりの量の有効画像成分が含まれているのも事実である。本来は、空間周波数ｆｍを中心とする急峻な選択特性を有するフィルタによるフィルタリングを行いたいところであるが、これを行ったとしても、対象画像に含まれる急激な変動部分を構成する周波数成分が含まれてしまうことに変りはない。
【０１４４】
そこで、上記の問題を解決するために、本実施の形態では、上述した第３ステップにおいて、グリッド縞成分の局所的な包絡線を求め、その変動から、グリッド縞成分以外の、アーチファクトを発生させる可能性のある成分が含まれる部分を検出することで、グリッド縞成分のみを安定に抽出（作成）する。
【０１４５】
一般の信号の包絡線は、ヒルベルト変換によらなければならないが、単一の正弦波の包絡線は、その周波数におけるレスポンス振幅が“１”であり、９０°（π／２）の位相変動を起こすような空間フィルタを施し、その結果と元信号とのベクトル振幅（２乗和の平方根）をとれば求まる。
【０１４６】
位相が９０°変動するＦＩＲフィルタによるフィルタリングを、離散的なデータに対して施す場合、当該ＦＩＲフィルタの係数を点対称（奇関数）的なものとする。例えば、この係数を（−ａ３，０，ａ３）とすると、空間周波数ｆｍでレスポンスを“１”にするためには、係数（−ａ３，０，ａ３）は、
２^*ａ３^*ｓｉｎ（２πｆｍＴ）＝１
なる式を満たす必要があり、
【０１４７】
【数３】

【０１４８】
なる式（３）で表される解が得られる。
【０１４９】
上記式（３）で表される解の係数（−ａ３，０，ａ３）を有するＦＩＲフィルタにより得られた信号系列と元信号系列との振幅を求める。
【０１５０】
例えば、図２（ａ）は、上記図２５（ａ）で示した元画像信号に対して、上記式（２）で示される解の係数（ａ１，ｂ１，ａ１）を有するフィルタリングを施した結果を示したものである。上記図２（ａ）から明かなように、グリッド縞成分が殆ど抽出されている。
【０１５１】
また、図３（ａ）は、上記図２６（ａ）で示した元画像信号に対して、上記式（２）で示される解の係数（ａ１，ｂ１，ａ１）を有するフィルタリングを施した結果を示したものである。
【０１５２】
上記図２（ｂ）及び図３（ｂ）で示す波形（太線で示す波形）はそれぞれ、上記図２（ａ）及び図３（ａ）で示される元画像信号に対して、上記式（３）で示される解の係数（−ａ３，０，ａ３）を有するフィルタリングを施した結果と、当該元画像信号との２乗和の平方根をとった包絡線を示したものである。
【０１５３】
特に、上記図３（ｂ）に示す包絡線において、同図（ｃ）に示す窪み部分に着目すると、当該窪み部分には明らかに非定常な成分が存在している。これは、抽出されたグリッド縞成分は、単純なフィルタリングでは異常に抽出され（対象画像のエッジ成分等を含み）、対象画像信号から差し引けば、この処理後の対象画像信号にアーチファクトが発生することを意味する。
【０１５４】
そこで、本実施の形態では、上述のようにして求められた包絡線から、異常な数値を示す範囲を特定し、当該範囲のグリッド縞成分をその周辺の数値列からの推定値で補正（置換）する。すなわち、本来有するグリッド縞成分の特徴である、全ての範囲にわたって常に定常な成分を有するという性質を利用して、グリッド縞成分を形成（作成）する。
【０１５５】
補正の際に用いる推定値（予測値）は、異常な数値を示す範囲の周辺のデータの統計的性質から求める。例えば、グリッド縞成分の空間周波数ｆｍが既知であるので、この空間周波数ｆｍを統計的性質として用いることができる。
【０１５６】
例えば、グリッド縞の空間周波数ｆｍ、及び位相φを以って、
【０１５７】
【数４】

【０１５８】
なる式（４）により表される正弦波を用いて、非定常部のグリッド縞成分を形成する。
【０１５９】
例えば、最も簡便な方法としては、フーリエ変換（フーリエ級数展開）を用い、特定の周波数における２つの係数Ａ及びφを周辺画素から求める方法が挙げられる。
【０１６０】
しかしながら、データに欠陥（非定常部分）が存在する等の問題から、通常のフーリエ変換を用いることができない。したがって、ここでは、フーリエ変換を一般化し、最小２乗の意味で、振幅及び位相情報を求める。このため、上記式（４）を、
【０１６１】
【数５】

【０１６２】
なる式（５）に変形する。
【０１６３】
この場合、サンプリング点ｘｉにおけるデータを“ｙｉ”（データ点数ｎ）である時（｛ｘｉ，ｙｉ；ｉ＝０〜ｎ−１｝）の２乗誤差εは、
【０１６４】
【数６】

【０１６５】
なる式（６）で表される。
【０１６６】
このとき、注意すべきは、ここで用いられる成分“ｘｉ，ｙｉ”としては、上述した包絡線の成分の検証から、定常な部分であると判断されたデータのみを選択することである。そして、２乗誤差εを最小化するパラメータＲ，Ｉを、次のようにして求める。
【０１６７】
先ず、
【０１６８】
【数７】

【０１６９】
なる式（７）は、
【０１７０】
【数８】

【０１７１】
なる式（７´）で書き表される。
上記式（７´）の連立方程式を解くことで、パラメータＲ，Ｉが求まり、位相φと振幅Ａを同時に推定できる。
【０１７２】
ここで、データ系列がｋ／（２・ｆｍ）（ｋは正の整数）の区間を等間隔でｍ等分するものであれば、上記式（７´）は、
【０１７３】
【数９】

【０１７４】
なる式（８）で示されるように、特定の周波数の係数を求める離散フーリエ変換（フーリエ級数展開）となる。
【０１７５】
上記式（７´）若しくは上記式（８）により、非定常部分の周辺の定常な適当なデータを用いて、パラメータＲ，Ｉの値を求めることで、不適当であるとして除去された非定常部分の補修（置換）を行う。
【０１７６】
また、その他の補修方法としては、線形予測モデルを考え、グリッド縞の空間周波数を特定せずに、線形予測アルゴリズムによって順次予測して補修を行う方法が挙げられる。
【０１７７】
上述のような補修によって得られた信号波形は、全般的に定常な正弦波であり、グリッド縞成分を非常によく表している成分である。
【０１７８】
しかしながら、上記の信号波形（グリッド縞成分の信号波形）は、もともと上記式（２）で示される係数（ａ１，ｂ１，ａ１）による狭いスパンのＦＩＲフィルタリングの結果であり、上記図１に示したようなフィルタの応答特性を有し、グリッド縞成分以外の画像成分をかなり多く含んでいるものである。
【０１７９】
そこで、本実施の形態では、上記の信号波形に対して、さらに、グリッド縞の空間周波数ｆｍ近辺の成分のみを抽出するフィルタリングを施す。このフィルタリングは、上述したような操作により、既に非定常な成分が補修された成分に対して施されるものであり、当該フィルタリングによるリンギング等のアーチファクトが発生することはない。
【０１８０】
上記のフィルタリング後のグリッド縞成分の信号に対して、上述した包絡線作成を行った際、非常に小さな値（“０”に近い値）が観測された部分が存在した場合、この部分は、元々グリッド縞成分が何らかの理由（例えば、Ｘ線が完全に遮断されている、又はセンサがサチュレーションを起こしている等）で観測されなかった部分であり、グリッド縞成分が元々存在しない部分である。したがって、この部分に関しては、その情報を記録しておき、後段のフィルタリングの後に“０”に置き換える。この処理の結果を、グリッド縞成分として、対象画像信号から差し引く。
【０１８１】
本実施の形態では、対象画像信号から順次対象とするラインデータを１ラインづつ取り出しながら、グリッド縞成分の抽出処理を実行するが、１ラインデータを取り出すときに、その前後の数ラインデータの平均を求め、すなわち画像成分を弱めてから、又はグリッド縞成分を強調してから、グリッド縞成分を抽出することも可能である。
【０１８２】
上記のことは、通常グリッド縞とセンサは、ほぼ平行に配置されており、あるラインのグリッド縞と、その近辺のラインのグリッド縞との各成分は、非常に酷似しているからである。
【０１８３】
したがって、本実施の形態の他の形態として、上記の非常に酷似しているという特徴を利用し、グリッド縞成分の抽出のための計算処理回数を減らすために、処理するラインを間引き、あるラインで抽出されたグリッド縞成分を、その近辺のラインのグリッド縞成分とする。すなわち、グリッド縞成分が抽出されたラインの近辺のラインについてはグリッド縞成分の抽出処理を行なわず、上記のあるラインで抽出されたグリッド縞成分を用いて、その近辺のラインデータからの差し引き処理を行なうこともできる。
【０１８４】
尚、上述の間引きが可能であるか否かを判断するために、上述した第１処理ステップにおいて、グリッド縞の空間周波数を測定する際に、サンプルされた前後のライン若しくはサンプルされたライン同士のグリッド縞の位相差を調べ、グリッド縞がセンサに対して傾いていないことを確認する。
【０１８５】
また、本実施の形態の他の形態としては、上述のようにして、グリッド縞成分を対象画像信号から除去する場合に、信号強度等を検出することで、対象画像信号における照射野を認識し、当該照射野内部の画像データに対してのみ、上述のグリッド縞成分の除去処理を行なうようにしてもよい。
【０１８６】
また、上述したようなグリッド縞成分の除去処理に用いるグリッド縞成分の抽出をＸ線動画を構成する全てのフレーム画像に対して実行すると多くの時間を必要とし、実際には、全てのフレーム画像に当該抽出処理を実行することは実時間では不可能であるため、本実施の形態では、Ｘ線動画を構成する複数のフレーム画像の中の選択フレーム画像（対象フレーム画像）について当該抽出処理を行い、これにより取得したグリッド縞成分を画像として保持し、対象フレーム画像に近隣の時間のフレーム画像に対して、同じグリッド縞成分を用いたグリッド縞成分除去処理を実行する。
【０１８７】
ところで、例えば、固体撮像素子によりＸ線画像を取得する場合、複数の画素を配列してなる当該固体撮像素子特有の問題として、欠陥画素の問題がある。画像情報の冗長性（空間的に低周波成分を主成分とする）から、欠陥画素が少量であれば、その周辺画素値からの平均的な補間により、殆どの場合、欠陥画素についての修復（補正）が可能である。
【０１８８】
しかしながら、一般的には、欠陥画素周辺の統計的な性質により、予測が必要となってくる。例えば、本実施の形態のように、グリッド縞がナイキスト周波数の５０％以上となると、平均補間では予測が逆転してしまう。
【０１８９】
図４は、１次元において、任意の点を欠陥画素とし、その両側の２点の画素の平均で補間する場合の、フィルタリングとしての応答関数を示したものである。
上記図４では、空間周波数を横軸で表している。上記図４に示すように、空間周波数が低く、ナイキスト周波数の５０％以下であれば、応答は“正”、すなわち位相が反転しない。これに対して、ナイキスト周波数の５０％以上となると、位相が反転し、期待される補間結果が得られない。
【０１９０】
図５は、欠陥画素補間の一例を示したものである。
上記図５において、黒丸点は、正常な画素から得られた画素値を表し、矢印で示す点（「欠陥画素位置」）は、データが得られていない欠陥画素を表している。また、上記図５は、グリッド縞が映り込んでいることにより、それぞれの画素データ（黒丸点で示すデータ）が細かく振動している状態を示している。
また、上記図５の「Ａ：平均による補間値」と指している白丸点は、従来の平均補間により得られた画素値であり、同図の「Ｂ：理想的な補間値」と指してある白丸点は、グリッド縞を考慮した補間値である。
【０１９１】
本実施の形態では、上記図５の「Ｂ：理想的な補間値」で示される理想的な補間値を得るため、次の２つの方法を実施する。
（方法１：線形予測の方法）
上記図５において、欠陥画素位置のデータを、その周辺の画素から線形予測して求める。
（方法２）
上述したグリッド縞の除去処理により、元画像信号からグリッド縞成分を除去することで、ナイキスト周波数の５０％以下の主成分を無くした状態で、従来の方法である平均による補間処理を行なう。
【０１９２】
以下に、（方法１）である線形予測の方法について、その概要を説明する。
【０１９３】
まず、処理対象の画像データ（画素データ）として、
データ系列｛Ｘ_n，Ｘ_n-1，Ｘ_n-2，・・・，Ｘ_n-p｝
が与えられ、“ｎ”におけるデータＸ_nが、
【０１９４】
【数１０】

【０１９５】
なる差分方程式（９）で表されるものとする。
上記式（９）において、“ε_n”は白色雑音系列を表し、“ａｉ｛ｉ＝１，・・・，ｐ｝”は線形予測係数を表す。このような系列を「自己回帰過程（ＡＲ過程）Ｘ_n」と呼ぶ。
【０１９６】
上記式（９）を、遅延演算子Ｚ^-1を用いて書き直すと、
【０１９７】
【数１１】

【０１９８】
なる式（１０）となる。
【０１９９】
但し、上記式（１０）は、
【０２００】
【数１２】

【０２０１】
なる式（１０´）で表されるため、ＡＲ過程Ｘ_nは、パルス伝達関数１／Ａ（ｚ_-1）を有する線形フィルタの入力ε_nに対する出力であると定義（スペクトル推定）できる。
【０２０２】
また、上記式（９）は、線形予測係数ａｉ｛ｉ＝１，・・・，ｐ｝が、信頼できるデータ系列から求まれば、（ｎ−１）点目の画素データから、ｎ点目の画素データが予測可能であることを示している。
【０２０３】
線形予測係数ａｉ｛ｉ＝１，・・・，ｐ｝の予測は、装置或いはシステムが定常であると過程して、最尤推定（最小２乗推定）を用いれば行える。すなわち、ε_nのパワー（分散）を最小にするものを求めることができる。ε_nは最小２乗推定によって得られた誤差であるため、その予測次数に対応する以下の相関成分を持たない。 したがって、必要十分な次数ｐで予測されて得られた予測誤差ε_nは式９で定義したように白色雑音となる。
【０２０４】
予測誤差ε_nの分散は、２乗平均（平均値“０”）であるため、当該平均を表す関数Ｅ[＊]は、
【０２０５】
【数１３】

【０２０６】
なる式（１１）で表される。
【０２０７】
上記式（１１）において、
Ｒ（τ）＝Ｅ［Ｘ_mＸ_m+τ］（共分散関数：ｃｏｖａｒｉａｎｃｅ）
であり、最小値を求めるために両辺を係数ａ_kで微分して“０”とおくと、
【０２０８】
【数１４】

【０２０９】
なる連立方程式（１２）が得られる。これは、正規方程式若しくはＹｕｌｅ−Ｗａｌｋｅｒ方程式と呼ばれるものである。
【０２１０】
実際には、自己相関Ｒ（＊）の演算は、全ての画素点で行うことなく、限られた（与えられた）画素数で演算した推定値を用いる。
例えば、高速算法であるＬｅｖｉｎｓｏｎアルゴリズムを用いるが、Ｂｕｒｇアルゴリズムを用いるようにしてもよい。このＢｕｒｇアルゴリズムは、さらに少ない画素数のデータで共分散（自己相関）を直接計算せずに求められる最大エントロピー法によるアルゴリズムである。これらのアルゴリズムでは、予測誤差が正規分布であり画素数が多ければ数学的には一致するが、画素数が少ない場合、Ｂｕｒｇアルゴリズム（最大エントロピー法によるアルゴリズム）が有利である。
【０２１１】
上述のような演算により得られる係数ａ_kを用いて、欠陥画素の前後の画素から予想される画素データを求める。
【０２１２】
[第１の実施の形態]
本発明は、例えば、図６に示すようなＸ線画像取得装置１００に適用される。
【０２１３】
＜Ｘ線画像取得装置１００の全体構成及び動作＞
本実施の形態のＸ線画像取得装置１００は、医療用（画像診断等）のＸ線動画を取得するための装置であり、上記図１に示すように、Ｘ線を被写体１０２（ここでは人体）に対して発生するＸ線発生部１０１と、被写体１０２からの散乱Ｘ線を除去するためのグリッド１０３と、被写体１０２を透過したＸ線量の分布を検出する面状のＸ線センサ１０４と、Ｘ線発生部１０１のコントローラ１０５（ＣＯＮＴ）と、Ｘ線センサ１０４から出力される電気信号をディジタルデータに変換するアナログ／ディジタル（Ａ／Ｄ）変換器１０６と、Ａ／Ｄ変換器１０６から出力されるディジタルデータをフレーム画像データとして一旦蓄積するメモリ１０７と、Ｘ線を放射しない状態での撮影により取得されたディジタルデータを記憶するメモリ１０８と、メモリ１０７内のフレーム画像データに対してメモリ１０８のデータを用いた演算処理を施す演算器１０９と、演算器１０９での演算処理後のフレーム画像データの変換テーブル（参照テーブル：Ｌｏｏｋ ＵｐＴａｂｌｅ、以下、「ＬＵＴ」とも言う）１１０と、Ｘ線センサ１０４を構成する画素毎のゲインのばらつきを補正するためゲインパタンデータを記憶するメモリ１１１と、ＬＵＴ１１０から出力される変換後のフレーム画像データに対してメモリ１１１のゲインパタンデータを用いた演算処理を施す演算器１１２と、演算器１１２での演算結果後のフレーム画像データを一旦記憶するメモリ１１３と、Ｘ線センサ１０４固有の欠陥画素に関しての情報（欠陥画素位置情報等）を記憶するメモリ１１４と、メモリ１１４内の情報を用いてメモリ１１３に記憶されたフレーム画像データに補正処理を施す補正処理部１１５と、メモリ１１３内の補正処理後のフレーム画像データの読み出しを制御するメモリ転送制御部１４２と、メモリ転送制御部１４２の制御によりメモリ１１３から読み出されたフレーム画像データを記憶するメモリ１４３と、メモリ１４３内のフレーム画像データに対してグリッド縞に関する情報を検出するグリッド縞検出部１１６と、グリッド縞検出部１１６で得られた情報に基いてメモリ１４３内のフレーム画像データからグリッド縞成分を抽出するグリッド縞成分抽出部１１７と、グリッド縞成分抽出部１１７で抽出されたグリッド縞成分を一時的に記憶するメモリ１１８と、メモリ１１８内のグリッド縞成分の読み出しを制御するメモリ転送制御部１４０と、メモリ転送制御部１４０の制御によりメモリ１１８から読み出されたグリッド縞成分を記憶するメモリ１４１と、メモリ１１３内のフレーム画像データからメモリ１４１内のグリッド縞成分を差し引く演算器１１９と、演算器１１９での演算結果(グリッド縞成分除去後のフレーム画像データ)を一旦記憶するメモリ１２０と、メモリ１２０内のフレーム画像データに画像処理を施して出力する画像処理部１２１とを備えている。
【０２１４】
上述のようなＸ線画像取得装置１００において、Ｘ線発生部１０１のコントローラ１０５は、不図示の操作部から操作者により発生トリガがかけられると、Ｘ線発生部１０１でのＸ線放射を開始する。
これにより、Ｘ線発生部１０１は、人体である被写体１０２に対して、Ｘ線を放射する。
【０２１５】
Ｘ線発生部１０１から放射されたＸ線は、被写体１０２を透過して、被写体１０２からの散乱Ｘ線を除去するグリッド１０３を介して、Ｘ線センサ１０４へと到達する。
【０２１６】
Ｘ線センサ１０４は、被写体１０２を透過したＸ線量の分布を検出する面（受像面）上に、当該Ｘ線強度を検出する複数の検出器（画素）がマトリックス状に配置された構成とされており、このマトリックス状に配置された複数の検出器（画素）により得られたＸ線強度に対応する電気信号を出力する。
【０２１７】
Ｘ線センサ１０４としては、例えば、次のようなセンサ（１）及び（２）を適用可能である。
センサ（１）：
Ｘ線強度を一旦蛍光に変換し、その蛍光をマトリックス状に配置されている複数の検出器で光電変換して検出するようになされたセンサ。
センサ（２）：
特定の物体に放射されたＸ線が該物体内で光電変換されて遊離した自由電子を、一様な電界によって引き付けて電荷分布を構成し、その電荷分布を、マトリックス状に配置された複数の電荷検出器（キャパシタ）によって電気信号に変換するようになされたセンサ。
【０２１８】
Ａ／Ｄ変換器１０６は、Ｘ線センサ１０４から出力された電気信号をディジタル化して出力する。
具体的には、Ａ／Ｄ変換器１０６は、Ｘ線発生部１０１のＸ線の放射、若しくはＸ線センサ１０４の駆動に同期して、Ｘ線センサ１０４から出力される電気信号を順次ディジタルデータに変換して出力する。
【０２１９】
尚、上記図６では、１つのＡ／Ｄ変換器１０６を設ける構成としているが、例えば、複数のＡ／Ｄ変換器を設け、これらを並行に動作させるように構成してもよい。これにより、ディジタル変換速度を早めることができ、効率よく処理を進めることができる。
【０２２０】
Ａ／Ｄ変換器１０６から出力されたディジタルデータは、フレーム画像データとしてメモリ１０７に一旦記憶される。
したがって、メモリ１０７には、Ｘ線センサ１０４を構成する複数の画素に対応する複数の画素データの集合であるディジタル画像データ（フレーム画像データ）が記憶される。
【０２２１】
メモリ１０８には、Ｘ線を放射しない状態での撮影により取得されたディジタルデータが予め記憶されている。このディジタルデータは、メモリ１０７に記憶されたフレーム画像データから、Ｘ線センサ１０４特有のオフセット的に存在する固定パタンノイズを除去するためのデータである。したがって、予め、Ｘ線画像取得装置１００において、Ｘ線発生部１０１によるＸ線を放射しない状態で撮影を行ない、これにより取得されたディジタルデータを画像データとして、メモリ１０８に記憶させておく。
【０２２２】
演算器１０９は、メモリ１０７に記憶されたフレーム画像データ（被写体１０２を透過したＸ線により得られた画像データ）を構成する複数の画素データのそれぞれに対して、メモリ１０８に記憶された画像データ（Ｘ線なしの撮影により得られた固定パタンノイズの画像データ）を構成する複数の画素データの中の対応する位置の画素データを減算する処理を実行する。
【０２２３】
ＬＵＴ１１０は、演算器１０９での処理後のフレーム画像データを、対数に比例した値に変換して出力する。
【０２２４】
メモリ１１１には、ＬＵＴ１１０による変換後のフレーム画像データに対して、Ｘ線センサ１０４を構成する各画素のゲインのばらつきを補正するためゲインパタンデータが記憶されている。このため、予め、Ｘ線画像取得装置１００において、被写体１０２がない状態でＸ線撮影を行ない、これにより得られた画像データから、メモリ１０８に記憶されたディジタルデータを用いて固定パタンノイズを除去し、さらにＬＵＴ１１０によって対数値に比例した値に変換して得られたデータを、ゲインパタンデータとしてメモリ１１１に記憶させておく。
【０２２５】
演算器１１２は、ＬＵＴ１１０から出力されたフレーム画像データから、メモリ１１１のゲインパタンデータを減算（対数変換されていなければ除算に相当）して出力する。
この演算器１１２での減算処理されたフレーム画像データは、メモリ１１３に一旦記憶される。
【０２２６】
尚、メモリ１１１に記憶させるゲインパタンデータに使用する画像データの取得の際に、グリッド１０３を装着した状態で撮影を行なえば、これにより得られるゲインパタンデータ自体に、グリッド縞が写り込むことになる。予想されるのは、演算器１１２により、対象画像データからゲインパタンデータを減算した際に、被写体１０２に写り込んだグリッド縞自体がゲインの変動に近いものであることにより、グリッド縞成分が除去される可能性があることである。しかしながら、被写体１０２なしの撮影で得られたゲインパタンデータは、画像取得毎（実際の撮影毎）に毎回取得される可能性は少なく、殆どの場合、１日一回、或いはさらに低い頻度で得られるものであり、また、Ｘ線発生部１０１とＸ線センサ１０４の位置関係は撮影毎に変化する可能性があるため、グリッド縞成分は上記の減算によって除去されない。また、上記位置関係が不変であっても、被写体１０２ありの撮影と、被写体１０２なし撮影なしの撮影とでは、一般にＸ線の散乱線量や線質が異なるため、グリッド縞のコントラストが異なり、グリッド縞成分は減算によって除去されない。
何れの場合にもグリッド１０３方向が一致していれば、グリッド縞の空間周波数が変動することはない。好適には、ゲインパタンデータを取得する場合には、グリッド１０３自体を取り外して、グリッド縞がゲインパタンデータに含まれないようにすべきである。
【０２２７】
メモリ１１４には、Ｘ線センサ１０４固有の欠陥画素に関しての情報（欠陥画素位置情報等）が記憶される。
具体的には例えば、一般に平面状のＸ線センサは、半導体製造技術で製造されるが、その歩留まりは１００％ではなく、製造工程での何らかの原因により、複数の検出器（画素）の内のいくらかは、検出器としての意味を無さない、すなわちその出力が意味を持たない欠陥画素である。ここでは、製造工程において、或いは不図示の手段によって、Ｘ線センサ１０４を予め検査し、その結果得られた欠陥画素の位置情報をメモリ１１４に記憶させておく。
【０２２８】
補正処理部１１５は、メモリ１１４に記憶された欠陥画素の位置情報により、メモリ１１３に記憶されたフレーム画像データを構成する複数の画素データの中の欠陥画素データを補正し、当該補正後の画素データを、再びメモリ１１３の該当する位置に記憶させる。
【０２２９】
メモリ転送制御部１４２は、後述するグリッド縞成分抽出部１１７から与えられる演算終了信号に基づいて、メモリ１１３内のフレーム画像データ（補正処理部１１５による補正処理後の画像データ）の読み出しを制御する。
メモリ１４３は、メモリ転送制御部１４２によりメモリ１１３から読み出されたフレーム画像データを記憶する。
【０２３０】
グリッド縞検出部１１６は、メモリ１４３内のフレーム画像データ（補正処理部１１５による補正処理後の画像データ）に対して、グリッド縞の解析を行い、グリッド縞の空間周波数ｆｍ、及びグリッド縞の角度θを検出して出力する。
【０２３１】
グリッド縞成分抽出部１１７は、メモリ１４３内のフレーム画像データ（補正処理部１１５による補正処理後の画像データ）を読み出し、当該読出画像データから、グリッド縞検出部１１６で得られたグリッド縞の空間周波数ｆｍ及びグリッド縞の角度θに基いて、グリッド縞成分を抽出し、その後、当該抽出処理終了を示す演算終了信号を、メモリ転送制御部１４２及びメモリ転送制御部１４０のそれぞれに対して出力する。
グリッド縞成分抽出部１１７で得られたグリッド縞成分は、メモリ１１８に一旦記憶される。
【０２３２】
メモリ転送制御部１４０は、グリッド縞成分抽出部１１７から与えられる演算終了信号に基づいて、メモリ１１８内のグリッド縞成分の読み出しを制御する。
メモリ１４１は、メモリ転送制御部１４０によりメモリ１１８から読み出されたグリッド縞成分を記憶する。
【０２３３】
演算器１１９は、メモリ１１３内のフレーム画像データ（補正処理部１１５による補正処理後の画像データ）から、メモリ１４１に記憶されたグリッド縞成分を差し引く。
演算器１１９によりグリッド縞成分が差し引かれたフレーム画像データは、メモリ１２０に一旦記憶される。
【０２３４】
画像処理部１２１は、メモリ１２０内のフレーム画像データに対して、観察者が観察しやすいように画像処理を施す。
ここでの画像処理としては、例えば、次のような処理が挙げられる。
・フレーム画像からのランダムノイズの除去処理。
・フレーム画像を表示した際に、観察者が見やすい濃度値になるように、階調を変換する、或いは詳細部分を強調する。
・フレーム画像から観察者にとって不要な部分を切り取り、フレーム画像の情報量を減らす、或いはフレーム画像情報を圧縮する。
【０２３５】
画像処理部１２１での処理後のフレーム画像データは、不図示の手段により、外部或いはＸ線画像取得装置１００内において、表示部への表示や、記憶部若しくは記憶媒体への格納、記録媒体への記録、或いは解析処理等が施される。
【０２３６】
＜Ｘ線画像取得装置１００の特徴とする構成及び動作＞
図７（ａ）〜（ｅ）は、Ｘ線画像取得装置１００において、特に、メモリ転送制御部１４０及びメモリ転送制御部１４２によるデータ読出及び当該読出データに対する処理のタイミングに着目して図示したものである。
【０２３７】
上記図７（ａ）は、Ｘ線センサ１０４で順次取得されるフレーム画像データがメモリ１１３へ記憶されるタイミングを示したものである。メモリ１１３へ順次記憶されたフレーム画像データは、演算器１１９により、メモリ１４１内のグリッド縞成分の除去処理が施される。
【０２３８】
上記図７（ｂ）は、グリッド縞成分抽出部１１７が出力する演算終了信号を示したものである。上記図７（ｂ）に示されるように、当該演算終了信号は、ハイレベル信号で、グリッド縞成分の抽出を終了（演算終了）したことを示す。
【０２３９】
上記図７（ｃ）は、グリッド縞成分抽出部１１７が出力する演算終了信号に基づいた、メモリ転送制御部１４２の動作タイミングを示したものである。上記図７（ｃ）に示されるように、メモリ転送制御部１４２は、グリッド縞成分抽出部１１７からの演算終了信号により、グリッド縞成分の抽出を終了したタイミングであれば（当該演算終了信号がハイレベル信号である場合）、メモリ１１３内のフレーム画像データをメモリ１４３に対して読み出す。
【０２４０】
上記図７（ｄ）は、グリッド縞検出部１１６及びグリッド縞成分抽出部１１７の動作タイミングを示したものである。上記図７（ｄ）に示されるように、メモリ１４３へフレーム画像データが転送されると、グリッド縞検出部１１６及びグリッド縞成分抽出部１１７は動作し、この結果得られたグリッド縞成分が、メモリ１１８へ記憶される。この動作が終了したタイミングで、グリッド縞成分抽出部１１７は、上記図７（ｂ）に示したようなハイレベルの演算終了信号を出力する。
【０２４１】
上記図７（ｅ）は、グリッド縞成分抽出部１１７が出力する演算終了信号に基づいた、メモリ転送制御部１４０の動作タイミングを示したものである。上記図７（ｅ）に示されるように、メモリ転送制御部１４０は、グリッド縞成分抽出部１１７からの演算終了信号により、グリッド縞成分の抽出を終了したタイミングであれば（当該演算終了信号がハイレベル信号である場合）、すなわちメモリ１１８へグリッド縞成分が記憶されたタイミングであれば、メモリ１１８内のグリッド縞成分をメモリ１４１に対して読み出す。これにより、メモリ１４１内のグリッド縞成分が更新されることになる。
【０２４２】
すなわち、本実施の形態では、上記図７（ａ）〜（ｅ）に示すように、Ｘ線動画を構成する全てのフレーム画像に対してグリッド縞成分の抽出処理を行うのではなく、１フレーム画像に対するグリッド縞成分の抽出処理が終了したタイミングで、次のグリッド縞成分の抽出処理を実行し、その間のフレーム画像に対して、直前に抽出された同一のグリッド縞成分（メモリ１４１に記憶されたグリッド縞成分）を用いたグリッド縞成分の除去処理を行う。このような動作は、一般的に、連続するフレーム画像では、Ｘ線発生部１０１の位置や、そのＸ線エネルギー、或いは被写体１０２の位置等がほぼ不変であるため、グリッド１０３による縞成分も殆ど変動がないことにより、有効となる動作である。
【０２４３】
尚、本実施の形態では、メモリ転送制御部１４２により、１フレーム画像に対するグリッド縞成分の抽出処理が終了したタイミングで、メモリ１４３へ次の抽出処理を行うフレーム画像を取り込むように構成したが、これに限られることはなく、例えば、他の任意のタイミング制御機構を用いて、定期的にメモリ１４３へフレーム画像を取り込み、当該フレーム画像に対するグリッド縞成分の抽出処理を行うようにしてもよい。
【０２４４】
また、本実施の形態では、Ｘ線動画を対象としているが、例えば、静止画を対象として、以前に抽出したグリッド縞成分を、今回取得した静止画に対するグリッド縞成分の除去処理に用いるようにしてもよい。この場合にも、グリッド縞成分の抽出に要する負担を軽減することができる。
【０２４５】
＜Ｘ線画像取得装置１００の他の具体的な構成及び動作＞
ここでは、上述したＸ線画像取得装置１００において、その他具体的な説明が必要と思われる、次のような構成部分について具体的に説明する。
（１）メモリ１１８に記憶されたグリッド縞成分の画像データ
（２）補正処理部１１５による欠陥画素の補正処理
（３）グリッド縞検出部１１６及びグリッド縞成分抽出部１１７によるグリッド縞成分の検出及び抽出処理
【０２４６】
（１）メモリ１１８に記憶されたグリッド縞成分の画像データ
メモリ１１８に記憶されたグリッド縞成分の画像データは、グリッド縞成分が重畳されたフレーム画像データから減算されるデータであるが、本実施の形態のように、メモリ１１８に記憶されるデータを減算後の対象画像データと対応づけて別途記憶する等の構成にすれば、グリッド縞が除去されたフレーム画像データから、元のグリッド縞が重畳されたフレーム画像データを再現できる。これにより、例えば、グリッド除去処理において、何らかの不具合によりフレーム画像データが損傷を受けた場合であっても、上記の再現処理により、元のフレーム画像データに戻すことが可能となる。
【０２４７】
（２）補正処理部１１５による欠陥画素の補正処理
補正処理部１１５は、以下に説明するような処理を、例えば、マイクロプロセッサを用いたソフトウエアにより実行する。
【０２４８】
図８〜図１１は、Ｘ線センサ１０４における画素欠陥の分布の例を示したものである。
ここでは、画素欠陥は基本的に１画素の幅でしか存在しないものとする。これは、大きなかたまりで隣接する複数の画素欠陥を有するＸ線センサは、欠陥画素の補修が困難であるので一般的に用いないためである。
【０２４９】
上記図８〜図１１において、それぞれの各マス目は画素を表し、黒いマス目は欠陥画素を表している。また、図の下部には、グリッド１０３のグリッド縞の方向（縦方向）を図示している。
【０２５０】
まず、上記図８に示す欠陥画素は基本的な画素欠陥であり、同図に示すように、欠陥画素（黒マス目）の周囲に、８個の隣接する画素成分ａ１〜ａ８が存在している。
【０２５１】
上記図８に示す欠陥画素が存在し、グリッド縞成分が存在する対象画像において、当該グリッド縞成分の空間周波数軸上の信号分布を模式的に示したものが、図１２である。
上記図１２において、横軸は、対象画像の横方向の空間周波数軸ｕを表し、縦軸は、対象画像の縦方向の空間周波数軸ｖを表し、空間周波数軸ｕ及び空間周波数軸ｖの両軸に対して、画素ピッチの逆数である「サンプリング周波数」とその半値である「ナイキスト周波数」を示している。
【０２５２】
グリッド縞は、対象画像の横方向に振動しており、縦方向には一定であるので、グリッド縞の成分は、上記図１２に示すように、空間周波数軸ｕ上に存在することになる（同図白丸参照）。
【０２５３】
通常の画像では、その主成分が、ナイキスト周波数の、さらに半値以下の空間周波数領域に分布しており、グリッド縞成分が存在しなければ、欠陥画素の任意の両側の画素値の平均により補間できる。これは、当該補間の空間スペクトルに与える影響が、例えば、上記図４で示したフィルタリングの応答関数（特性）を示すためである。
【０２５４】
したがって、上記図８に示す欠陥画素の補正の場合、縦方向にはグリッド縞成分がないことにより、縦方向の画素の平均、すなわち画素成分ａ２及画素成分ａ６の平均、或いは何れか一方の画素成分により、ほぼ満足な補正が可能となる。
【０２５５】
上記図９に示す欠陥画素は、横長に連結する画素欠陥である（同図中、黒マス目参照）。このような形状の欠陥画素の場合も、上記図８に示した欠陥画素と同様に、各画素の上下方向がグリッド縞に並行する方向であるため、対象欠陥画素の上下の画素成分により、ほぼ満足な補正が可能となる。
【０２５６】
上記図１０に示す欠陥画素は、縦長に連結する画素欠陥である（同図中、黒マス目参照）。このような形状の欠陥画素の場合、連結欠陥画素の上端の欠陥画素或いは下端の欠陥画素以がいについては、対象画素の上下に信頼できる値を有する画素が存在しない。このような状態の欠陥画素に対して、横方向の単純な平均等で欠陥補正を行えば、上記図５を用いて説明したような、期待しない値の補正結果が得られてしまう。
【０２５７】
そこで、本実施の形態では、対象欠陥画素の右又は左側に連なる正常な画素成分を以って、上記式（１２）で示したような連立方程式を利用し、係数ａ_k（ｋ＝１〜Ｐ）を、対象欠陥画素の左右から求める。
このとき、例えば、使用する画素数を２０程度とし、次数ｋを５程度とする。
【０２５８】
そして、係数ａ_kを以って、上記式（９）により対象欠陥画素値Ｘ_nを予測し、この結果得られた全ての欠陥画素値Ｘ_nの平均を求める。これにより、上記図５に示したような「Ｂ：理想的な補間値」が得られる。
【０２５９】
尚、本実施の形態では、上記式（１２）を用いて、係数ａ_k（ｋ＝１〜Ｐ）を求めるようにしているが、これに限られることはなく、例えば、最大エントロピー法と呼ばれるアルゴリズム等を用いるようにしてもよい。
【０２６０】
上記図１１に示す欠陥画素は、上記図９に示した欠陥画素の状態、及び上記図１０に示した欠陥画素の状態を重ね合わせた状態の欠陥画素である。この状態の欠陥画素の中で問題となる画素は、縦方向の連結欠陥画素と、横方向の連結欠陥画素との交わる部分の画素（十字に重なった部分の画素）、すなわち画素成分ａ４，ａ１２，ａ１８，ａ２４で囲まれた欠陥画素である。
【０２６１】
上記図１１に示すような状態の欠陥画素の補正は、横方向に連なった線状の欠陥画素の補正は、上下画素成分の平均で行い、縦方向に連なった線状の欠陥画素の補正は、上述したような連立方程式を用いて行う。
具体的には、次の３つの方法▲１▼〜▲３▼が挙げられるが、何れの方法を用いても、ほぼ同じ結果が得られる。
【０２６２】
▲１▼画素成分ａ４，ａ１２，ａ１８，ａ２４に囲まれた欠陥画素の補正を、上下の欠陥補正値（補正された欠陥画素の値）の平均値を用いて行なう。
▲２▼画素成分ａ４，ａ１２，ａ１８，ａ２４に囲まれた欠陥画素の補正を、欠陥補正値である左右画素の値を用いて、上記式（１２）の連立方程式を解くことにより行う。
▲３▼画素成分ａ４，ａ１２，ａ１８，ａ２４に囲まれた欠陥画素の補正を、▲１▼の結果と▲２▼の結果の平均値を用いて行なう。
【０２６３】
以上説明したような処理を実行することで、メモリ１１３に記憶されたフレーム画像データを構成する複数の画素データの中の欠陥画素データが補正される。
【０２６４】
（３）グリッド縞検出部１１６及びグリッド縞成分抽出部１１７によるグリッド縞成分の検出及び抽出処理
【０２６５】
グリッド縞検出部１１６は、メモリ１１３内に記憶されたフレーム画像データの一部を読み出し、当該読出データにより、フレーム画像データに含まれるグリッド縞のスペクトルを調べ、当該グリッド縞の空間周波数ｆｍ及び角度θを検出する。このグリッド縞の空間周波数ｆｍ及び角度θ（以下、「角度η」とも言う）の情報は、後段のグリッド縞成分抽出部１１７において、グリッド縞成分の抽出処理に使用される。
【０２６６】
図１３（ａ）〜（ｄ）は、グリッド縞の空間周波数ｆｍ及び角度θの検出処理を説明するための図である。
【０２６７】
上記図１３（ａ）は、対象となるフレーム画像全体のイメージを示したものであり、“Ｌ１”乃至“Ｌ６”は、対象フレーム画像上部からのライン位置を表している。グリッド縞検出部１１６は、ラインＬ１〜Ｌ６をフーリエ変換した結果により、グリッド縞の空間周波数ｆｍを測定する。このとき、グリッド縞検出部１１６は、グリッド縞のスペクトルを検出する際、当該検出能力を上げるために、各ラインＬ１〜Ｌ６の前後数ラインの平均（又はスペクトルの平均）を用いるようにしてもよい。
【０２６８】
上記図１３（ｂ）〜（ｄ）はそれぞれ、ラインＬ１におけるフーリエ変換結果を表したものである。
すなわち、上記図１２（ｂ）は、振幅スペクトル（又はパワースペクトル）を表し、同図（ｃ）は、フーリエ変換結果の余弦波の係数である実数部の値を表し、同図（ｄ）は、正弦波の係数である虚数部の値を表している。
【０２６９】
図１４は、グリッド縞検出部１１６の上記の処理をフローチャートによって示したものである。
【０２７０】
先ず、グリッド縞検出部１１６は、スペクトルの平均を求めるための変数ｃｕｍｕｌａｔｉｏｎをクリアする（ステップＳ２０１）。
また、グリッド縞検出部１１６は、スペクトルの平均を求める際の対象とライン数のカウンタ（変数）ｎをクリアする（ステップＳ２０２）。
また、グリッド縞検出部１１６は、上記図１３（ａ）に示したラインＬ１〜Ｌ６の中から処理対象ライン（選択ライン）を選択する変数ｉを“１”に初期設定する（ステップＳ２０３）。これにより、最初の処理では、ラインＬ１が対象ラインとして選択され処理されることになる。
【０２７１】
そして、グリッド縞検出部１１６は、次のステップＳ２０５〜ステップＳ２１５の処理を、対象フレーム画像のラインＬ１〜Ｌ６の全てのラインについて実行し終えたか否かを判別する（ステップＳ２０４）。
この判別の結果、処理終了した場合のみ、後述するステップＳ２１６へ進み、未だ処理終了していない場合には、次のステップＳ２０５からの処理を実行する。
【０２７２】
ステップＳ２０４の判別の結果、処理未終了の場合、先ず、グリッド縞検出部１１６は、対象フレーム画像のラインＬ１〜Ｌ６の中から、変数ｉで示されるラインＬｉを選択し、そのデータ（ラインデータＬｉ）を取得する（ステップＳ２０５）。
【０２７３】
次に、グリッド縞検出部１１６は、ステップＳ２０５で取得したラインデータＬｉに対して、高速フーリエ変換等のフーリエ変換処理を施す（ステップＳ２０６）。
【０２７４】
次に、グリッド縞検出部１１６は、ステップＳ２０６でのフーリエ変換結果（空間周波数領域のデータ）から、パワースペクトル（又は振幅スペクトル）を取得する（ステップＳ２０７）。
【０２７５】
次に、グリッド縞検出部１１６は、ステップＳ２０７で取得したパワースペクトルにおいて、グリッド縞を示す有意なスペクトル（ピーク値）が存在するか否かを判別する（ステップＳ２０８）。
【０２７６】
具体的には例えば、まず、グリッド縞を発生させる原因となるグリッド鉛の絶対的な空間周波数は、グリッド１０３を設置した段階で既知であることにより、その周波数を“ｆｇ”として用いることで、ステップＳ２０８での判別処理を正確に行える。
【０２７７】
すなわち、Ｘ線センサ１０４のサンプリングピッチを“Ｔｓ”とすると、グリッド縞の発生する大まかな空間周波数ｆｍは、
【０２７８】
【数１５】

【０２７９】
なる条件式（１３）により特定できる。
【０２８０】
このとき、上記式（１３）において、▲３▼に示される条件を満たしている場合には、▲２▼で得られた空間周波数ｆｍ´を用い、▲３▼に示される条件を満たしていない場合には、「Ｊ←Ｊ＋１」として▲１▼を実行する。
【０２８１】
グリッド縞の正確な空間周波数ｆｍは、上記式（１３）で得られる“ｆｍ´”の近辺に存在するはずであり、ピーク値（グリッド縞を示す有意なスペクトル）が存在するか否かを判断する際に、当該近辺のみを検索すれば、画像成分やノイズ成分等の影響で異なるピーク値が存在したとしても、その影響を受けることなく、グリッド縞を示す有意なスペクトルであるピーク値の検出を行なえる。
【０２８２】
また、グリッド１０３のグリッド鉛の周波数ｆｇは、かなり正確に製造されるものであるが、撮影の際に、グリッド１０３とＸ線センサ１０４との間に任意の距離以上あると、Ｘ線発生部１０１からのＸ線ビームがコーンビーム状であることにより、当該Ｘ線ビームが拡大されてＸ線センサ１０４に到達してしまう。このため、正確な空間周波数ｆｍが異なるものになってしまい、簡単に予測することができない。
したがって、上記式（１３）で得られる“ｆｍ´”近辺の周波数のうち、ピーク値を示す周波数を、空間周波数ｆｍとして求める。
【０２８３】
但し、上記の場合、求められたピーク値に相当するものが、実際に安定して存在する有意なピーク値であるか否かを判断する必要がある。この判断は、通常ノイズレベルを基準に行う。このノイズレベルとしては、予め測定されるものでも構わないし、スペクトルの高域のピーク値以外の成分の平均値を代用するようにしてもよい。例えば、ピーク値を示した近隣のスペクトル値のパワースペクトルの総和（又は平均値）と、ノイズレベルとの比が、１０程度以上あれば、通常有意な安定したピーク値であると判断する。
【０２８４】
上述のようなステップＳ２０８の判別の結果、有意なピーク値が存在しない場合、グリッド縞検出部１１６は、次のラインを処理するために、変数ｉをカウントアップして（ステップＳ２１７）、再びステップＳ２０４へと戻り、これ以降の処理ステップを繰り返し実行する。
【０２８５】
一方、ステップＳ２０８の判別の結果、有意なピーク値が存在する場合、グリッド縞検出部１１６は、当該ピーク値を示す空間周波数をＰｉとして（ステップＳ２０９）、これを変数ｃｕｍｕｌａｔｉｏｎに対して加算する（ステップＳｓ２１０）。
【０２８６】
また、グリッド縞検出部１１６は、グリッド縞の位相を求め、変数θｎ及び変数Ｍｎに対して、当該位相及び変数ｉに示されるライン位置を設定し（ステップＳ２１１〜ステップＳ２１４）、変数ｎをインクリメントする（ステップＳ２１５）。
その後、グリッド縞検出部１１６は、次のラインを処理するために、変数ｉをカウントアップして（ステップＳ２１７）、再びステップＳ２０４へと戻り、これ以降の処理ステップを繰り返し実行する。
【０２８７】
上述のようにして、全てのラインＬ１〜Ｌ６に対してステップＳ２０５〜ステップＳ２１５の処理を実行し終えると、グリッド縞検出部１１６は、現在の変数ｃｕｍｕｌａｔｉｏｎの値を、現在の変数ｎの値（有意なピーク値数）で除算することで、平均的なグリッド縞の空間周波数ｆｍを求める（ステップＳ２１６）。
【０２８８】
また、グリッド縞検出部１１６は、上記図１４の処理実行後、ステップＳ２１３で得られた変数θｉ（ｉ＝０〜ｎ−１）を、ステップＳ２１４でのライン位置Ｍｉ（ｉ＝０〜ｎ−１）により、平均的なグリッド縞の角度η（角度θ）を求めることに用いる。
【０２８９】
すなわち、グリッド縞検出部１１６は、ｉ番目のラインＬｉの位相θ_iと、（ｉ＋１）番目のラインＬ（ｉ＋１）の位相θ_i+1との位相差（θ_i−θ_i+1）によるグリッド１０３の位置の差を、グリッド縞の空間周波数ｆｍを以って、
[（θ_i−θ_i+1）/２π]／ｆｍ
なる式により求め、ラインＬｉとラインＬ（ｉ＋１）のライン差を、
（Ｍｉ＋１）−Ｍｉ
なる式で求め、これらの結果を以って、グリッド縞の角度ηを、
【０２９０】
【数１６】

【０２９１】
なる式（１４）により求める。
【０２９２】
そして、グリッド縞検出部１１６は、上記式（１４）により得られた角度ηが異常に傾いていない場合、グリッド縞を抽出する処理を数ライン毎に実行し、これに対して当該角度ηが異常に傾いている場合、グリッド縞を抽出する処理を１ライン毎を実行する。
【０２９３】
尚、上述したグリッド縞検出部１１６が実行する処理では、グリッド縞の方向（縦又は横方向）が既知であることを前提としたが、例えば、グリッド縞の方向も不明である場合、例えば、予め縦方向と横方向の双方に対して同様の処理を実行し、有意なピーク値が検出された方向を、グリッド縞に略直交する方向とする。
【０２９４】
以上説明したような処理をグリッド縞検出部１１６が実行することで、グリッド縞の空間周波数ｆｍ及び角度θ（角度η）が求められる。
グリッド縞成分抽出部１１７は、グリッド縞検出部１１６により得られたグリッド縞の空間周波数ｆｍ及び角度θ（角度η）を用いて、実際にメモリ１４３に記憶されている、グリッド縞成分を含む対象フレーム画像データから、グリッド縞成分を抽出し、そのグリッド縞成分をメモリ１１８ヘ格納する。
【０２９５】
図１４は、グリッド縞成分抽出部１１７でのグリッド縞成分抽出処理をフローチャートによって示したものである。
【０２９６】
先ず、グリッド縞成分抽出部１１７に対しては、処理パラメータとして、グリッド縞検出部１１６にて得られたグリッド縞の角度θ及びグリッド縞の空間周波数ｆｍが与えられる。
グリッド縞成分抽出部１１７は、上記処理パラメータ（グリッド縞の角度θ及びグリッド縞の空間周波数ｆｍ）に基いて、以下に説明するステップＳ３００〜ステップＳ３２１の処理を実行する。
【０２９７】
尚、グリッド縞成分抽出部１１７に対して与えられたグリッド縞の空間周波数ｆｍの値が“０”の場合等は、対象フレーム画像上にグリッド縞が存在しない、すなわちグリッド１０３を使用せずに撮影が行なわれた場合であるので、この場合、グリッド縞成分抽出部１１７は、上記図１５に示される処理を実行しない。また、例えば、メモリ１１８には、この場合のグリッド縞成分として“０”データが格納される。或いは、演算器１１９が機能しないことにより、メモリ１２０に対して、メモリ１１３に格納された対象画像データがそのまま格納される。
【０２９８】
グリッド縞検出部１１６からグリッド縞成分抽出部１１７に対して、グリッド縞の角度θ及びグリッド縞の空間周波数ｆｍが与えられると、先ず、グリッド縞成分抽出部１１７は、上記式（２）を用いて、空間周波数ｆｍから係数ａ１，ａ２（又は対称５点フィルタの係数（ａ２，ｂ２，ｃ２，ｂ２，ａ２））を求める（ステップＳ３００）。ここで得られた係数に対応するＦＩＲフィルタを「ＦＩＲ１」とする。
【０２９９】
次に、グリッド縞成分抽出部１１７は、上記式（３）を用いて、空間周波数ｆｍから係数ａ３を求める（ステップＳ３０１）。ここで得られた係数に対応するＦＩＲフィルタを「ＦＩＲ２」とする。
【０３００】
次に、グリッド縞成分抽出部１１７は、空間周波数ｆｍの領域でのＦＩＲフィルタリングを行うために、空間周波数ｆｍを中心とするウインドウ関数を生成する（ステップＳ３０２）。
ここでのウインドウ関数としては、例えば、空間周波数ｆｍを中心としたガウス分布形状の関数を適用可能である。
【０３０１】
次に、グリッド縞成分抽出部１１７は、グリッド縞の角度θに基いて、対象フレーム画像を構成するラインデータに対してグリッド縞の抽出処理を実行するラインの範囲を決定する（ステップＳ３０３〜ステップＳ３０５）。
【０３０２】
具体的には例えば、グリッド縞成分抽出部１１７は、角度θの基準値を「０．１度」とし、グリッド縞検出部１１６で得られたグリッド縞の角度θが基準値０．１度よりも大きいか或いは小さいかを判別する（ステップＳ３０３）。この判別の結果、角度θが基準値０．１度よりも小さい場合、グリッド縞成分抽出部１１７は、変数ｋｉｐに対して「５」を設定することで、５ライン毎に当該処理を省くことを決定する（ステップＳ３０４）。一方、角度θが基準値０．１度以上の場合、グリッド縞成分抽出部１１７は、変数ｓｋｉｐに対して「０」を設定することで、全てのラインに対して当該処理を実行することを決定する（ステップＳ３０５）。
【０３０３】
尚、ステップＳ３０３〜ステップＳ３０５において、例えば、変数ｓｋｉｐに対する設定だけではなく、角度θに基き、さらに細かな設定を行なうようにしてもよい。
【０３０４】
ステップＳ３０４又はステップＳ３０５の処理後、グリッド縞成分抽出部１１７は、変数ｃｏｕｎｔに対して、ステップＳ３０４又はステップＳ３０５にて設定が行なわれた変数ｓｋｉｐの値を設定することで、変数ｃｏｕｎｔの初期化を行なう（ステップＳ３０６）。
【０３０５】
そして、グリッド縞成分抽出部１１７は、現在の変数ｃｏｕｎｔの値が、変数ｋｉｐの値以上であるか否か、すなわち対象ラインデータに対するグリッド縞の抽出処理を実行すべきであるか否かを判別する（ステップＳ３０７）。
この判別の結果、処理実行する場合には、ステップＳ３１０からの処理に進み、処理実行でない場合には、ステップＳ３０８からの処理に進む。但し、最初に本ステップＳ３０７の実行時には、必ず処理実行と判別されるため、次のステップＳ３１０へ進む。
【０３０６】
ステップＳ３０７の判別の結果、処理実行の場合（ｓｋｉｐ≦ｃｏｕｎｔ）、先ず、グリッド縞成分抽出部１１７は、メモリ１４３に格納されている対象フレーム画像データから処理対象となる１ラインのデータ（対象ラインデータ）を取得する（ステップＳ３１０）。
【０３０７】
尚、ステップＳ３１０において、対象フレーム画像データから対象ラインデータをそのまま取得するようにしてもよいが、例えば、対象ラインデータの前後数ライン分の平均値（移動平均値）を、実際の処理対象のラインデータとして取得するようにしてもよい。
【０３０８】
次に、グリッド縞成分抽出部１１７は、ステップＳ３１０で取得した対象ラインデータに対して、ステップＳ３００で係数を決定したＦＩＲ１を用いたフィルタリングを施し、ラインバッファ１（不図示）へ格納する（ステップＳ３１１）。
ここでの処理により、ラインバッファ１に対しては、グリッド縞成分を含む画像成分のデータが格納される。
【０３０９】
次に、グリッド縞成分抽出部１１７は、ラインバッファ１に格納されたラインデータに対して、ステップＳ３０１で係数を決定したＦＩＲ２を用いたフィルタリングを施し、ラインバッファ２（不図示）へ格納する（ステップＳ３１２）。
ここでの処理により、ラインバッファ２に対しては、グリッド縞成分の包絡線を求めるためのデータが格納される。
【０３１０】
次に、グリッド縞成分抽出部１１７は、グリッド縞成分の包絡線を求める（ステップＳ３１３）。
すなわち、グリッド縞成分抽出部１１７は、ラインバッファ１内のデータと、ラインバッファ２内のデータを成分とするベクトルの振幅（すなわち２乗和の平方根）を求め、その結果を、ラインバッファ３（不図示）へ格納する。ここでの演算としては、平方根の単調増加性により、平方根を取らない演算であっても適用可能であり、同様の効果が得られる。
【０３１１】
次に、グリッド縞成分抽出部１１７は、ラインバッファ３内のデータ、すなわち包絡線データを調査して、その異常データを検出するためのしきい値の上限値ｔｈ１及び下限値ｔｈ２を決定する（ステップＳ３１４）。
上限値ｔｈ１及び下限値ｔｈ２の決定方法としては、様々な方法を適用可能であるが、例えば、平均値と標準偏差値を求め、平均値から標準偏差値のｎ倍（ｎは、例えば、“３”程度の値）以上ずれた値を上限値ｔｈ１及び下限値ｔｈ２とする方法や、ラインバッファ３内の包絡線データのヒストグラムを求め、その最頻値を中心として上限値ｔｈ１及び下限値ｔｈ２を決定する方法等が挙げられる。
【０３１２】
次に、グリッド縞成分抽出部１１７は、ラインバッファ３内の包絡線データにおいて、値ｔｈ１以上若しくは値ｔｈ２以下のデータを異常データ（画像データが急峻に変動していること等によるデータ）と見なし、その異常データに対応するラインバッファ１内のグリッド縞成分のデータを、当該異常データ周辺のデータから推定して書き換える（ステップＳ３１５）。このとき、全体のグリッド縞成分が周期的な変動パタンを示す安定な状態となるように、データ書き換えを行う。
尚、ステップＳ３１５の処理については、上記式（７´）等の説明部分で説明したので、ここではその詳細は省略する。
【０３１３】
次に、グリッド縞成分抽出部１１７は、ステップＳ３１５の処理により、全体的に安定状態となったラインバッファ１内のグリッド縞成分のデータに対して、フーリエ変換処理を施し、グリッド縞成分の空間周波数領域のデータを求める（ステップＳ３１６）。
尚、ステップＳ３１６での変換処理については、フーリエ変換に限らず、例えば、コサイン変換等の他の直交変換をも適用可能である。
【０３１４】
次に、グリッド縞成分抽出部１１７は、ステップＳ３１６で取得した空間周波数領域のデータに対して、ステップＳ３０２で求めた空間周波数ｆｍを中心とするウインドウ関数によるフィルタリングを施す（ステップＳ３１７）。これにより、グリッド縞成分のデータは、より選択的にグリッド縞を表すようになる。
【０３１５】
次に、グリッド縞成分抽出部１１７は、ステップＳ３１７のフィルタリング処理後のグリッド縞成分のデータに対して、ステップＳ３１６の変換の逆変換処理を施し、この結果を、実際のグリッド縞成分のデータとする（ステップＳ３１８）。
【０３１６】
次に、グリッド縞成分抽出部１１７は、ステップＳ３１８で取得したグリッド縞成分のデータを、メモリ１１８の当該位置へ格納する（ステップＳ３１９）。
【０３１７】
そして、グリッド縞成分抽出部１１７は、変数ｃｏｕｎｔに対して“０”を設定し（ステップＳ３２０）、メモリ１４３の対象フレーム画像データを構成する全てのラインデータについて、ステップＳ３０７からの処理を行ったか否かを判別する（ステップＳ３２１）。
この判別の結果、処理が終了した場合、グリッド縞成分抽出部１１７は、上述した演算終了信号を出力し、本処理終了とする。一方、未だ処理が終了していない場合には、再びステップＳ３０７へ戻り、以降の処理ステップを繰り返し実行する。
【０３１８】
上述したステップＳ３０７の判別の結果、処理実行でない場合（ｓｋｉｐ＞ｃｏｕｎｔ）、グリッド縞成分抽出部１１７は、メモリ１１８の該当位置に対して、前段で取得したグリッド縞成分のデータをコピーし（ステップＳ３０８）、コピーするラインを示す変数ｃｏｕｎｔをインクリメントして（ステップＳ３０９）、再びステップＳ３０７へ戻り、以降の処理ステップを繰り返し実行する。
【０３１９】
[第２の実施の形態]
本発明は、例えば、図１６に示すようなＸ線画像取得装置４００に適用される。
本実施の形態のＸ線画像取得装置４００は、上記図６に示したＸ線画像取得装置１００とは、以下の構成が異なる。
【０３２０】
尚、上記図１６のＸ線画像取得装置４００において、上記図６のＸ線画像取得装置１００と同様に動作する構成部には同じ符号を付し、その詳細な説明は省略する。
【０３２１】
上記図６に示したＸ線画像取得装置１００では、補正処理部１１５が、メモリ１１３内のフレーム画像データに対して、メモリ１１４内のデータを用いた欠陥画素の補正処理を施すように構成した。これに対して、本実施の形態のＸ線画像取得装置４００では、上記図１６に示すように、補正処理部１１５が、メモリ１２０内のフレーム画像データ、すなわちグリッド縞成分の除去後のフレーム画像データに対して、メモリ１１４内のデータを用いた欠陥画素の補正処理を施すように構成した。
【０３２２】
したがって、本実施の形態のＸ線画像取得装置４００によれば、グリッド縞を考慮した画素欠陥補正が必要なくなるため、従来から知られるような、周辺の欠陥ではない画素値の平均値を用いて欠陥画素を補正する等のような、単純な欠陥画素補正を適用可能となる。
【０３２３】
[第３の実施の形態]
本発明は、例えば、図１７に示すようなＸ線画像取得装置５００に適用される。
本実施の形態のＸ線画像取得装置５００は、上記図６に示したＸ線画像取得装置１００とは、以下の構成が異なる。
【０３２４】
尚、上記図１７のＸ線画像取得装置５００において、上記図６のＸ線画像取得装置１００と同様に動作する構成部には同じ符号を付し、その詳細な説明は省略する。
【０３２５】
本実施の形態のＸ線画像取得装置５００は、上記図１７に示すように、上記図６のＸ線画像取得装置１００の構成に対して、さらに、グリッド１０３の装着を検知する検知部（スイッチ）１２２を設けた構成としている。
【０３２６】
検知部１２２は、グリッド１０３の装着の検知結果（グリッド装着信号）を、補正処理部１１５及びグリッド縞検出部１１６へそれぞれ供給する。
【０３２７】
補正処理部１１５は、検知部１２２からのグリッド装着信号により、グリッド１０３が装着されている場合、第１の実施の形態で説明したような、グリッド縞を考慮した欠陥画素補正処理を実行する。そうでない場合は、周辺の欠陥でない画素の画素値の平均値等により欠陥画素を補正する。
【０３２８】
グリッド縞検出部１１６も同様に、検知部１２２からのグリッド装着信号により、グリッド１０３が装着されている場合、第１の実施の形態で説明したような、グリッド縞の検出処理（解析処理）を実行する。但し、上記グリッド装着信号により、グリッド１０３が装着されていない場合、グリッド縞検出部１１６は、グリッド縞の検出処理を実行せずに、即座にグリッド縞無しと判断し、これに該当する処理を実行する。
【０３２９】
上述のように、本実施の形態のＸ線画像取得装置５００では、検出部１２２を設け、この検出結果に基づいて、グリッド縞の検出を行うように構成したので、グリッド縞を検出する処理に要する時間を大幅に短縮できる。
【０３３０】
[第４の実施の形態]
本発明は、例えば、図１８に示すようなＸ線画像取得装置６００に適用される。
本実施の形態のＸ線画像取得装置６００は、上記図６に示したＸ線画像取得装置１００とは、以下の構成が異なる。
【０３３１】
尚、上記図１８のＸ線画像取得装置６００において、上記図６のＸ線画像取得装置１００と同様に動作する構成部には同じ符号を付し、その詳細な説明は省略する。
【０３３２】
本実施の形態のＸ線画像取得装置６００は、上記図１８に示すように、上記図６のＸ線画像取得装置１００の構成に対して、さらに、Ｘ線照射領域データを格納するためのメモリ１２３を設けた構成としている。
【０３３３】
メモリ１２３には、メモリ１１３に格納されたフレーム画像データにおいて、Ｘ線が照射された領域部分のみを切り出した画像データ（照射領域データ）が格納され、このメモリ１２３内の照射領域データに対して、グリッド縞成分の検出及び抽出が行なわれる。
【０３３４】
具体的には、まず、Ｘ線撮影では一般に、被写体１０２（ここでは、人体）の目的とする部位以外への被曝を避けるために、Ｘ線発生部１０１のＸ線発生管球の出口に照射野絞りを設けることが行なわれる。これにより、被写体１０２の必要部分のみにＸ線の照射が行える。
【０３３５】
上記の照射野絞り機能を用いた場合、Ｘ線撮影により得られた画像も、Ｘ線センサ１０４から得られる画像信号全てが有効ではなく、照射野絞りによるＸ線照射野に対応する部分画像のみが有効になる。
【０３３６】
そこで、本実施の形態では、不図示の計算機手段（ＣＰＵ等）により、メモリ１１３内のフレーム画像データから、Ｘ線強度分布や形状、或いはその他の情報に基づいて、Ｘ線照射野に対応する有効な部分画像領域（照射領域）を見出し、その照射領域部分のデータ（照射領域データ）のみを、メモリ１２３に格納する。
【０３３７】
上述のように、本実施の形態では、メモリ１２３内の照射領域データ、すなわち対象画像データ全てではなく、情報量を削減した必要部分のデータのみを対象とするので、処理時間の短縮化を実現できる。
【０３３８】
尚、本実施の形態では、欠陥画素補正後のフレーム画像から、照射領域を切り出すように構成したが、例えば、当該照射領域の切り出し後に、当該照射領域に対して欠陥画素補正を行うようにしてもよい。
【０３３９】
[第５の実施の形態]
第１の実施の形態では、上記図６のＸ線画像取得装置１００において、グリッド縞成分抽出部１１７でのグリッド縞成分抽出処理を、上記図１５のフローチャートに従った処理とした。
本実施の形態では、グリッド縞成分抽出部１１７でのグリッド縞成分抽出処理を、例えば、図１９に示すフローチャートに従った処理とする。
【０３４０】
尚、上記図１９のグリッド縞成分抽出処理において、上記図１５のグリッド縞成分抽出処理と同様に処理するステップには同じ符号を付し、その詳細な説明は省略する。
【０３４１】
まず、本実施の形態でのグリッド縞成分抽出処理の説明の前に、図２０（ａ）は、グリッド縞成分を含むフレーム画像の一部分を示したものであり、同図において、“＊”で示す部分は、Ｘ線を遮断する物質が存在する等の原因により、グリッド縞成分が存在しない部分を示している。
【０３４２】
上記図２０（ｂ）は、同図（ａ）のフレーム画像に対して、第１の実施の形態でのフィルタタリングにより、グリッド縞成分の抽出を行った結果を示したものである。上記図２０（ｂ）に示すように、グリッド縞成分において、フレーム画像の“＊”で示す部分に対応する部分にアーチファクトが現れるため、第１の実施の形態で説明したように、当該部分を包絡線から抽出し、その前後のデータから推測して全体を安定したグリッド縞になるようにする。この結果を示したものが、上記図２０（ｃ）の図である。このように安定したグリッド縞成分であれば、フーリエ変換等の変換処理により、新たなアーチファクトは発生しない。
【０３４３】
第１の実施の形態では、上記図２０（ｃ）に示されるようなグリッド縞成分を、同図（ａ）に示されるような元のフレーム画像から差し引いて、グリッド縞成分を除去するように構成したが、この構成の場合、本来グリッド縞成分が存在しない部分（“＊”で示す部分）に新たなグリッド縞成分が現れてしまうことが考えられる。
【０３４４】
そこで、本実施の形態では、上記図２０（ｄ）に示すように、同図（ｃ）に示されるようなグリッド縞成分において、本来グリッド縞成分が存在しない部分（“＊”で示す部分）を“０”に設定する。
【０３４５】
このため、本実施の形態では、グリッド縞成分抽出部１１７は、上記図１９のフローチャートに従ったグリッド縞成分抽出処理を実行する。
すなわち、グリッド縞成分抽出部１１７は、ステップＳ３１８の処理実行後、この処理により取得した、安定したグリッド縞成分のデータに対して、ステップＳ３１５により推測で補った部分を“０”に置換する処理を施し（ステップＳ７００）、その後、次のステップＳ３１９へ進む。これにより、もともとグリッド縞成分が存在しない部分であっても、グリッド縞除去処理により新たなグリッド縞成分が発生してしまうことを確実に防ぐことができる。
【０３４６】
尚、第１〜第５の実施の形態では、ハードウェア的に構成したが、本装置全体をソフトウェアにより制御することで実現することも可能である。
【０３４７】
また、本発明の目的は、第１〜第５の実施の形態のホスト及び端末の機能を実現するソフトウェアのプログラムコードを記憶した記憶媒体を、システム或いは装置に供給し、そのシステム或いは装置のコンピュータ（又はＣＰＵやＭＰＵ）が記憶媒体に格納されたプログラムコードを読みだして実行することによっても、達成されることは言うまでもない。
この場合、記憶媒体から読み出されたプログラムコード自体が第１〜第５の実施の形態の機能を実現することとなり、そのプログラムコード、及びそのプログラムコードを記憶した記憶媒体は本発明を構成することとなる。
プログラムコードを供給するための記憶媒体としては、ＲＯＭ、フレキシブルディスク、ハードディスク、光ディスク、光磁気ディスク、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、磁気テープ、不揮発性のメモリカード等を用いることができる。
また、コンピュータが読みだしたプログラムコードを実行することにより、第１〜第５の実施の形態の機能が実現されるだけでなく、そのプログラムコードの指示に基づき、コンピュータ上で稼動しているＯＳ等が実際の処理の一部又は全部を行い、その処理によって第１〜第５の実施の形態の機能が実現される場合も含まれることは言うまでもない。
さらに、記憶媒体から読み出されたプログラムコードが、コンピュータに挿入された拡張機能ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれた後、そのプログラムコードの指示に基づき、その機能拡張ボードや機能拡張ユニットに備わるＣＰＵなどが実際の処理の一部又は全部を行い、その処理によって第１〜第５の実施の形態の機能が実現される場合も含まれることは言うまでもない。
【０３４８】
図２１は、上記のコンピュータの機能８００の構成例を示したものである。
コンピュータ機能８００は、上記図２１に示すように、ＣＰＵ８０１と、ＲＯＭ８０２と、ＲＡＭ８０３と、キーボード（ＫＢ）８０９のキーボードコントローラ（ＫＢＣ）８０５と、表示部としてのＣＲＴディスプレイ（ＣＲＴ）８１０のＣＲＴコントローラ（ＣＲＴＣ）８０６と、ハードディスク（ＨＤ）８１１及びフレキシブルディスク（ＦＤ）８１２のディスクコントローラ（ＤＫＣ）８０７と、ネットワーク８４０との接続のためのネットワークインターフェースコントローラ（ＮＩＣ）８０８とが、システムバス８０４を介して互いに通信可能に接続された構成としている。
【０３４９】
ＣＰＵ８０１は、ＲＯＭ８０２或いはＨＤ８１１に記憶されたソフトウェア、或いはＦＤ８１２より供給されるソフトウェアを実行することで、システムバス８０４に接続された各構成部を総括的に制御する。
すなわち、ＣＰＵ８０１は、所定の処理シーケンスに従った処理プログラムを、ＲＯＭ８０２、ＨＤ８１１或いはＦＤ８１２から読み出して実行することで、第１〜第５の本実施の形態での動作を実現するための制御を行う。
【０３５０】
ＲＡＭ８０３は、ＣＰＵ８０１の主メモリ或いはワークエリア等として機能する。
ＫＢＣ８０５は、ＫＢ８０９や図示していないポインティングデバイス等からの指示入力を制御する。
ＣＲＴＣ８０６は、ＣＲＴ８１０の表示を制御する。
ＤＫＣ８０７は、ブートプログラム、種々のアプリケーション、編集ファイル、ユーザファイル、ネットワーク管理プログラム、及び第１〜第６の実施の形態における所定の処理プログラム等を記憶するＨＤ８１１及びＦＤ８１２へのアクセス等を制御する。
ＮＩＣ８０８は、ネットワーク８４０上の他の装置或いはシステムとの双方向のデータのやりとりを制御する。
【０３５１】
【発明の効果】
以上説明したように本発明では、グリッドを使用した放射線撮影による放射線画像から、グリッドに起因する画像成分が除去された良好な放射線画像を得ることのできる、放射線画像処理装置、画像処理システム、放射線画像処理方法、記憶媒体及びプログラムを提供することができる。
【図面の簡単な説明】
【図１】第１〜第５の実施の形態において、対象画像からグリッド縞成分を抽出するためのフイルタの空間周波数特性を説明するための図である。
【図２】上記フィルタにより対象画像から抽出されたグリッド縞成分の一例を説明するための図である。
【図３】上記フィルタにより対象画像から抽出されたグリッド縞成分の他の例を説明するための図である。
【図４】対象画像の欠陥画素補正における空間周波数特性を説明するための図である。
【図５】上記グリッド縞成分が存在する対象画像に対する欠陥画素補正における空間周波数特性を説明するための図である。
【図６】第１の実施の形態において、本発明を適用したＸ線画像取得装置の構成を示すブロック図である。
【図７】上記Ｘ線画像取得装置の動作タイミングを説明するための図である。
【図８】上記Ｘ線画像取得装置での対象画像において、欠陥画素の様子の一例（例１）を説明するための図である。
【図９】上記欠陥画素の様子の一例（例２）を説明するための図である。
【図１０】上記欠陥画素の様子の一例（例３）を説明するための図である。
【図１１】上記欠陥画素の様子の一例（例４）を説明するための図である。
【図１２】上記Ｘ線画像取得装置での対象画像において、グリッド縞成分の空間周波数分布を説明するための図である。
【図１３】上記Ｘ線画像取得装置での対象画像に対する、グリッド縞成分の検出（解析）を説明するための図である。
【図１４】上記グリッド縞成分の検出（解析）処理を説明するためのフローチャートである。
【図１５】上記グリッド縞成分の検出（解析）処理の結果に基づいて、対象画像からグリッド縞成分を抽出する処理を説明するためのフローチャートである。
【図１６】第２の実施の形態において、本発明を適用したＸ線画像取得装置の構成を示すブロック図である。
【図１７】第３の実施の形態において、本発明を適用したＸ線画像取得装置の構成を示すブロック図である。
【図１８】第４の実施の形態において、本発明を適用したＸ線画像取得装置の構成を示すブロック図である。
【図１９】第５の実施の形態において、上記グリッド縞成分の検出（解析）処理の結果に基づいて、対象画像からグリッド縞成分を抽出する処理を説明するためのフローチャートである。
【図２０】上記グリッド縞成分を抽出する処理を具体的例を挙げて説明するための図である。
【図２１】第１〜第５の実施の形態の機能をコンピュータに実現させるためのプログラムを記録したコンピュータ読出可能な記憶媒体から当該プログラムを読み出して実行する構成の一例を示すブロック図である。
【図２２】グリッドを用いた放射線撮影を説明するための図である。
【図２３】上記放射線撮影で得られた画像に対してのフィルタリングの効果の一例を説明するための図である。
【図２４】上記放射線撮影で得られた画像に対してのフィルタリングの効果の他の例を説明するための図である。
【図２５】上記放射線撮影で得られた、グリッド縞成分が重畳した画像に対してのフィルタリングの効果の一例を説明するための図である。
【図２６】上記放射線撮影で得られた、グリッド縞成分が重畳した画像に対してのフィルタリングの効果の一例を説明するための図である。
【符号の説明】
１００ Ｘ線画像取得装置
１０１ Ｘ線発生部
１０２ 被写体
１０３ グリッド
１０４ Ｘ線センサ
１０５ コントローラ
１０６ アナログ／デジタル（Ａ／Ｄ）変換器
１０７ メモリ
１０８ メモリ
１０９ 演算器
１１０ 変換テーブル
１１１ メモリ
１１２ 演算器
１１３ メモリ
１１４ メモリ
１１５ 補正処理部
１１６ グリッド縞検出部
１１７ グリッド縞成分抽出部
１１８ メモリ
１１９ 演算器
１２０ メモリ
１２１ 画像処理部
１４０ メモリ転送制御部
１４１ メモリ
１４２ メモリ転送制御部
１４３ メモリ[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a radiographic image processing apparatus used in an apparatus or system for acquiring a radiographic image of a subject by radiography using, for example, a grid for removing scattered radiation. , Release The present invention relates to a ray image processing method, a storage medium, and a program.
[0002]
[Prior art]
Conventionally, for example, a technique that makes it possible to visualize the inside of a subject by irradiating the subject (object) with radiation represented by X-rays and imaging the spatial distribution (radiation distribution) of the radiation component transmitted through the subject. Although used, the radiation generates scattered radiation (scattered rays) inside the subject, so that the scattered radiation (scattered rays) is also imaged together with the directly transmitted radiation (direct rays) transmitted through the subject.
[0003]
The process of generating scattered radiation depends on the type of radiation, the physical properties, structure, etc. of the subject and is usually unpredictable. For this reason, in order to obtain a radiation image from which scattered radiation has been removed, various ideas are required.
[0004]
In a method for easily obtaining a radiation image from which scattered radiation has been removed, a typical method is to provide a radiation shielding material wall such as lead on the radiation receiving surface and regulate the angle of the radiation reaching the receiving surface. Then, there is a method of shielding the scattered radiation component.
[0005]
Specifically, for example, in X-ray photography in the medical field, in order to remove scattered X-rays from a subject such as a human body, an instrument called a “grid” may be installed between the subject and the X-ray image receiving surface. Has been done.
[0006]
The grid is angled so that the X-ray shielding material such as lead and the X-ray transmitting material such as wood, paper, aluminum, or carbon face each other toward the X-ray generation source (focal point). The instrument is configured by alternately arranging with a predetermined width.
[0007]
Since the grid as described above excludes a part of the direct line, a shadow of the grid (a shadow of the grid, hereinafter also referred to as “grid stripe”) is seen on the X-ray receiving surface. By adopting a configuration such as “accurate spatial period in which shielding substances and X-ray transmitting substances are alternately arranged” and “making the period relatively high in spatial frequency”, it is possible for an observer of an X-ray image. On the other hand, the uncomfortable feeling that grid stripes exist on the X-ray image is reduced as much as possible.
[0008]
FIG. 22 schematically shows the configuration of the grid 940.
In FIG. 22, “910” indicates an X-ray generation source, and “920” indicates an X-ray emission direction. “930” indicates a subject such as a human body, and “950” indicates an X-ray image receiving surface.
[0009]
As shown in FIG. 22, the grid 940 generally has a stripe structure only in one direction on the two-dimensional plane (vertical direction indicated by an arrow in the figure) for reasons such as ease of manufacture. It is.
[0010]
Further, as a method of reducing (removing) the contrast of grid stripes, during X-ray exposure, only the grid is moved in a direction orthogonal to the stripes, and the integration effect during X-ray exposure on the X-ray receiving surface is obtained. There is a way to use.
[0011]
Here, the radiation image receiving surface (image receiving portion) mainly means a radiation film for recording a radiation distribution directly on a photosensitive material.
[0012]
By the way, in recent years, a method of handling a radiographic image as digital data is becoming more general than a method of directly imaging a radiographic image for medical use or the like with a radiographic film.
For example, the radiation distribution is once converted into an electrical signal (analog signal), and the analog signal is digitized (A / D conversion) and converted into numerical data (digital data). Thereby, processing such as filing, image processing, and monitor display can be performed at a low cost on a digitized radiation image (digital radiation image).
[0013]
However, in the case of a digital radiographic image, it is necessary to sample an image signal in a two-dimensional space, so naturally the problem of aliasing based on the sampling theorem becomes obvious.
[0014]
Specifically, in the case of a general image, there is no problem in observing the image even if some aliasing occurs by setting the sampling period in the space to an appropriate period (such as a fine period).
[0015]
On the other hand, in the case of digital radiographic images obtained using a grid, the periodic striped pattern due to the grid becomes very low frequency due to aliasing, or aliasing does not occur depending on the sampling period. However, there may be a case where low-frequency amplitude fluctuations occur, which is a problem for the observer of the image.
[0016]
Therefore, the following method has been proposed as a method for removing an inappropriate grid stripe pattern by aliasing or the like.
Here, the method of digitizing a radiographic image is roughly classified into two methods, and grid stripe patterns for these two methods are removed (method 1) and (method 2).
[0017]
(Method 1)
In the digitization method here, first, the radiation distribution (radiation intensity distribution) is once converted into another energy distribution (energy distribution such as fluorescence), and this is scanned to spatially sample the image only in one direction. An analog video signal (time signal) is generated, and the time signal is A / D converted based on a temporal cycle prepared separately.
As a typical example of the digitization method, there is a method in which the energy stored in the photostimulable phosphor is sequentially fluorescentized / condensed by laser scanning to obtain a video signal, and then A / D conversion is sequentially performed. Can be mentioned.
[0018]
In acquiring a digital radiographic image by the above digitization method, as a method of removing an inappropriate grid stripe pattern by aliasing or the like, for example, Japanese Patent No. 2507659, Japanese Patent No. 2754068, and Japanese Patent Laid-Open No. Hei 8- There is a method described in No. 088765.
[0019]
Specifically, first, scanning with respect to another energy distribution corresponding to the radiation distribution is performed in a direction orthogonal to the grid stripe, the grid stripe is converted into a periodic signal on the video signal, and the analog periodic signal is obtained. On the other hand, after performing low-pass filtering, sampling of the time axis is performed. Such a normal anti-aliasing filter configuration can remove inappropriate grid stripes.
[0020]
As a method similar to the above-described method, Japanese Patent No. 2507659 and the like detect the presence and frequency of a grid-striped pattern by performing Fourier transform on a pre-sampled image, and based on the detection result. A method of removing inappropriate grid stripes by selecting a low-pass filter and performing low-pass filtering using the low-pass filter has been proposed.
[0021]
In addition, in Japanese Patent No. 2754068 and Japanese Patent Application Laid-Open No. 8-088765, the analog low-pass filtering processing proposed in Japanese Patent No. 2507659 etc. is performed, and sampling on the time axis is performed shorter than a desired interval. After acquiring the image information including the grid image without eliminating the aliasing of the grid stripe pattern information, digital low-pass filtering is performed to remove the grid image component, and then digital thinning (resampling) is performed. A method for obtaining an image at a desired sampling interval has been proposed.
[0022]
(Method 2)
In the digitization method here, first, the radiation intensity distribution is once converted into another energy distribution (energy distribution such as fluorescence and electric field intensity), and a plurality of two-dimensionally arranged electric signal conversions are performed. Two-dimensional sampling is directly performed by an element (a photodiode, a capacitor, or the like), and signals from each conversion element are sequentially extracted and A / D converted.
As a representative example of the digitization method, an electrical signal is obtained by a conversion element for each of a plurality of pixels in which a fluorescence distribution or a field intensity distribution over a large area due to radiation is two-dimensionally arrayed by a large screen planar sensor which has been developed in recent years. The method using what is called a radiation flat panel sensor is mentioned.
[0023]
By the way, when acquiring a digital radiation image by the digitization method of Method 2, it is very difficult to remove an inappropriate grid stripe pattern by aliasing or the like. The reason for this is that sampling is performed directly in a two-dimensional space by a plurality of electrical signal conversion elements such as a radiation flat panel sensor (hereinafter also simply referred to as “sensor” or “flat panel sensor”). This is because an anti-aliasing filter for a typical electric signal cannot be applied.
[0024]
In order to solve the above-mentioned problem, direct sampling is performed with a sensor so fine as not to cause aliasing in a two-dimensional space, a digital anti-aliasing filter is applied, and the above-described re-sampling is performed. Although such a method is conceivable, it is difficult to perform fine sampling in a two-dimensional space due to the configuration of the electric signal conversion element of the sensor, and the performance of the electric signal conversion element itself is degraded. Invite up.
[0025]
Therefore, a method of moving the grid during X-ray exposure as in the past has been performed.
[0026]
As another method, in Japanese Patent Laid-Open No. 9-75332, etc., when acquiring a digital X-ray image by directly sampling in a two-dimensional space by a sensor, the purpose is to remove the grid stripes. A method has been proposed to prevent the generation of inappropriate grid stripes on the image by matching the pitch (sensor pixel pitch) completely and matching the gap between the pixel and the area that blocks the direct line due to grid stripes. ing.
[0027]
JP-A-9-78970 and U.S. Pat. S. In Patent 5801385 and the like, the grid stripe contrast is reduced by making the grid stripe interval smaller than the sampling pitch so as to be the same as or close to the width of the opening of the light receiving portion of one pixel (one electric signal conversion element). A reduction method has been proposed.
[0028]
In addition, U.S. S. In Patent 5.050.198 and the like, a grid striped pattern image (grid image) is stored in advance under a plurality of shooting conditions, and when a grid image is actually captured, a plurality of grid images stored in advance are stored. Among them, a method of removing a grid image by dividing an image obtained by photographing by a grid image that matches or approximates photographing conditions has been proposed.
[0029]
[Problems to be solved by the invention]
However, in the conventional radiographic image processing as described above, in particular, when radiographic images are acquired using sensors and grids directly sampled in a two-dimensional space represented by (Method 2), There was a problem like this.
[0030]
First, in the configuration proposed in JP-A-9-75332 or the like, it is very difficult to completely match the interval between grid stripes and the sampling pitch. That is, a flat panel sensor generally created in a semiconductor manufacturing process and a grid created by a combination of relatively thick lead plates are each independently created, or the grid itself is Since it must be created so as to be removable depending on the situation at that time, it is very difficult to perfectly match the grid stripe interval and the sensor pixel pitch (sampling pitch) due to these factors.
[0031]
Also, JP-A-9-78970 and U.S. Pat. S. In the configuration proposed in Patent 5801385 and the like, it is effective to make the interval between the grid stripes smaller than the sampling pitch so as to be the same as or close to the width of the opening of the light receiving portion of one pixel, but the sensor (flat panel sensor) However, when the sampling pitch becomes 0.1 mm or less, for example, the interval between grid stripes is required to be very fine, that is, 10 or more per mm. In such a fine-striped grid, the lead plate thickness for blocking scattered radiation is almost fixed, so the area through which the direct line passes must be narrowed, and as a result, the radiation dose utilization efficiency Becomes extremely low, and good radiography cannot be performed.
[0032]
Further, in the conventional configuration as described above, the grid itself is moved during radiation exposure. However, the movement of the grid increases the cost by a drive system for moving the grid, and the apparatus or system. It is necessary to provide a configuration for adjusting and controlling the relationship between the drive timing and the radiation exposure timing and the relationship between the drive speeds. Therefore, the grid movement configuration is effective for removing grid stripes, but has the above-mentioned limitations and cannot always be adopted.
[0033]
Therefore, in order to solve the above problem, since the obtained radiation image is digital data, a method of removing grid stripes by digital filtering can be considered. According to this method, if the spatial frequency component of the grid stripe and the spatial frequency component of the effective image information based on the subject are completely separated, the grid stripe can be removed with a simple filtering configuration.
[0034]
As a specific example of the above method, Japanese Patent Laid-Open No. 3-12785 proposes a method for removing or reducing data corresponding to the spatial frequency of grid stripes using Fourier transform.
[0035]
A method of removing or reducing data corresponding to the spatial frequency of grid stripes using a normal FIR (Finite Impulse Response) filter is also conceivable.
[0036]
By the way, since the image of the grid stripe is a shadow with reduced radiation transmittance due to an X-ray shielding material such as lead, it is superimposed almost multiply in the signal, but is additively superimposed if logarithmic transformation is performed. It will be. Therefore, filtering as described above is possible.
[0037]
In general, the manufacturing process of the grid for removing scattered radiation is managed with very high accuracy, and the grid stripe having a uniform spatial frequency characteristic over the entire image is also widely used. . For this reason, there is a possibility that the filtering as described above may be performed only for a single spatial frequency.
[0038]
Moreover, in reality, the shape of the grid stripe image (shadow) is not an accurate sine wave shape, so there may be a spatial frequency component of double, triple,... In this case, it is considered that only the fundamental wave component is received due to the blur caused by the two-dimensional spatial dependence of the conversion process (energy conversion process) at the sensor.
[0039]
However, the problem with filtering as described above is that it is almost impossible to limit the spatial frequency band of the image component itself.
[0040]
Specifically, as typified by, for example, those proposed in Japanese Patent No. 2754068 and Japanese Patent Laid-Open No. 8-088765, the spatial sampling pitch is made very fine, and the effective Nyquist frequency is increased after the sampling, If the effective bandwidth (bandwidth below the Nyquist frequency) is widened and the image component and the grid stripe component are completely separated, the grid stripe can be removed by normal filtering without any problem. Although it is clear, as a practical matter, increasing the spatial sampling pitch just to remove grid fringes leads to a significant cost increase of the sensor due to semiconductor processes and other factors, and further increases the radiation acquisition efficiency. This is not effective.
[0041]
Therefore, it is effective in terms of cost and performance to configure the sensor itself with a spatial sampling pitch that allows an effective image component in a band below the Nyquist frequency, but if such a configuration is used, the grid Inevitably, there is a certain degree of overlap between the spatial frequency of the fringes and the spatial frequency of the effective image components.
[0042]
Specifically, for example, referring to FIGS. 23A to 23D, first, FIG. 23A shows the image when the image (original image) to be processed is viewed in one dimension. The image signal is composed of 256 numerical values.
FIG. 23 (b) shows the response characteristics of the filter in the spatial frequency domain when the image signal shown in FIG. 23 (a) is filtered. In FIG. 23 (b), in consideration of the discrete Fourier transform, the frequency domain is represented by numerical values from “0” to “128”, and band cut filtering is performed at a spatial frequency value of “100”. .
FIG. 23 (c) shows the result of applying the filtering shown in FIG. 23 (b) to the image signal shown in FIG. 23 (a). As is clear from FIG. 23C, the image signal in FIG. 23C has almost no change from the characteristics of the image signal shown in FIG.
FIG. 23D shows the difference between the image signal shown in FIG. 23A and the image signal shown in FIG. 23C for confirmation. As apparent from FIG. 23 (d), there is almost no signal component removed by filtering.
[0043]
FIG. 24 (a) shows the image signal (original image signal) shown in FIG. 23 (a) added with a portion that rises sharply in the middle (so-called edge portion). .
FIG. 24 (b) shows the response characteristics of the filter in the spatial frequency domain when the image signal shown in FIG. 24 (a) is filtered, as in FIG. 23 (b). .
FIG. 24C shows the result of applying the filtering shown in FIG. 24B to the image signal shown in FIG. In FIG. 24 (c), it is clear from the circled portion that it is out of the original image signal and unstable vibration (artifact).
FIG. 24D shows a difference between the image signal shown in FIG. 24A and the image signal shown in FIG. 24C for confirmation. As is clear from FIG. 24 (d), many vibration components appear in the portion (including both ends of the signal) that fluctuates sharply.
[0044]
As shown in FIGS. 23A to 23D and FIGS. 24A to 24D, in the case of a normal image signal, the Nyquist frequency (spatial frequency is “128” in the same figure) or less The high-frequency component is not a main component of the image signal and is a component that has almost no information. Therefore, even if a sharp filtering process is performed at this position, there are few problems. On the other hand, when the image signal has a steep portion (edge portion), the image signal is expressed using a considerably high frequency component equal to or lower than the above Nyquist frequency. Problems (artifacts) occur in the fluctuating part.
[0045]
25A to 25D show signal states when a sine wave (sin (2π100x / 256)) is applied to the original image signal shown in FIG. Is shown. As apparent from FIGS. 25A to 25D, the grid stripes are substantially removed by the filtering having the response characteristic of the filter shown in FIG. 25B (see FIG. 25C).
[0046]
26A to 26D show signal states when a sine wave (sin (2π100x / 256)) is applied to the original image signal shown in FIG. Is shown. As is clear from FIGS. 26A to 26D, the same artifact as in FIG. 24C is generated by the filtering having the response characteristic of the filter shown in FIG. (See FIG. 26 (c)).
[0047]
That is, in a simple filtering process such as that described in Japanese Patent Laid-Open No. 3-12785 or the like, if the grid stripe component is removed, there is a possibility that the above-described artifacts are strongly generated. Further, if the width of the impulse response of the filter is narrowed in order to reduce artifacts, the response characteristics of filtering will be reduced in a wide range, and the image itself will become stronger.
[0048]
The problem in the radiographic image as described above is also applicable to a moving image obtained by radiography.
[0049]
Specifically, for example, in the case of a configuration in which an X-ray image of a subject is acquired by capturing an X-ray with a solid-state imaging device through the subject, a series of continuous images (X of a moving image) corresponding to the movement of the subject. Line image (hereinafter also referred to as “X-ray moving image”). Even for such an X-ray moving image, the problems related to the grid stripe component as described above occur in the individual frame images constituting the X-ray moving image.
[0050]
In the case of acquiring an X-ray moving image, the X-ray moving image is acquired at a high speed of about 30 to 120 frames per second. Thus, it is very difficult to completely remove the grid stripe component from the X-ray moving image obtained at high speed. Furthermore, since the acquisition of the X-ray moving image is performed at a high speed, if the grid stripe component is removed by simple filtering, the image itself may be damaged.
[0051]
Therefore, the present invention is made to eliminate the above-mentioned drawbacks, and can obtain a good radiographic moving image in which image components due to the grid are removed from a radiographic image obtained by radiography using the grid. Radiation image processing device , Release An object is to provide a ray image processing method, a storage medium, and a program.
[0052]
[Means for Solving the Problems]
Under such an object, the first invention is a radiographic image processing apparatus for processing a radiographic image obtained by radiography using a grid for removing scattered radiation from a subject, wherein the predetermined radiographic image is obtained. A plurality of radiation images obtained after the predetermined radiation image based on the image component obtained from the image component obtained by the creation means and the image component obtained by the creation means; And removing means for removing the components The creating means includes means for creating the image component based on the feature that the image component superimposed on the radiation image is steady throughout the image. It is characterized by that.
[0072]
A second invention is a radiographic image processing method for processing a radiographic image obtained by radiography using a grid for removing scattered radiation from a subject, and an image resulting from the grid from a radiographic image A creation step of creating a component, and a removal step of performing a removal process of the image component on a plurality of radiation images obtained after the radiation image based on the image component obtained by the creation step. Thus, the creating step includes the step of creating the image component based on the feature that the image component superimposed on the radiographic image is steady over the entire image. It is characterized by that.
[0073]
First 3 The invention of the above 2 In the present invention, the removing step comprises Multiple above For the radiographic image, use the same image component obtained in the above creation step. The The method includes a step of performing a removal process.
[0075]
First 4 According to the second aspect of the invention, in the second aspect of the invention, the creating step includes an analysis step of analyzing the radiographic image to obtain a spatial frequency of the image component.
[0076]
According to a fifth invention, in the fourth invention, the creation step is obtained by an extraction step of extracting a predetermined component including an image component from the radiation image based on an analysis result of the analysis step, and the extraction step. A processing step of processing the predetermined component obtained to obtain the image component; And the removing step is the same as the creating step. The image component obtained by the processing step is removed from the radiation image. Ruko And features.
[0077]
First 6 The invention of the above 5 In the invention, the extracting step includes a step of performing filtering for extracting the component having the spatial frequency obtained in the analyzing step from the radiation image.
[0078]
First 7 The invention of the above 5 In the invention, the processing step includes a step of executing processing for estimating and processing an unsteady portion of the predetermined component from a steady portion before and after the predetermined portion.
[0079]
First 8 The invention of the above 7 In the invention described above, the processing step includes a step of detecting the unsteady portion based on envelope information of the predetermined component.
[0080]
First 9 The invention of the above 7 In the invention, the processing step estimates the amplitude and phase of a sine wave corresponding to the image component based on a stationary part before and after the unsteady part of the predetermined component and a spatial frequency, and the unsteady part. And a step of repairing a necessary part.
[0081]
First 10 The invention of the above 9 In the invention, the processing step includes a step of further filtering the predetermined component including the repaired non-stationary portion to obtain the image component.
[0082]
First 11 The invention of the above 7 In the invention, the processing step is characterized in that the image component is obtained by replacing a portion satisfying a predetermined condition in the non-stationary portion with a predetermined value.
[0083]
First 12 According to the second aspect of the invention, in the second aspect of the invention, the creating step includes a step of creating the image component for a predetermined line selected from the radiation image.
[0084]
First 13 The invention of the above 12 In the invention, the creating step includes a step of creating the image component with respect to an average result of a plurality of lines.
[0085]
First 14 The invention of the above 12 In the invention, the creating step includes a step of using, as an image component of a line for which no image component has been created, the image component acquired from a line in which the image component adjacent to the line is created. To do.
[0087]
First 15 According to the second aspect of the invention, in the second aspect of the invention, the method includes a detection step of detecting mounting of the grid, and the creation step includes a step of executing the processing based on a detection result of the detection step. And
[0092]
First 16 The computer-readable storage medium recording the program according to the invention is the second to the second. 15 A program for causing a computer to execute the processing steps of the radiation processing method according to any of the inventions is stored.
[0094]
First 17 The program according to the invention of No. 2 15 A processing step of the radiation processing method according to any of the inventions is caused to be executed by a computer.
[0095]
DETAILED DESCRIPTION OF THE INVENTION
Hereinafter, embodiments of the present invention will be described with reference to the drawings.
[0096]
[Outline of the first to fifth embodiments]
Here, as embodiments of the present invention, the first to fifth embodiments are exemplified. Prior to specific description of the first to fifth embodiments, an outline of these will be described.
[0097]
Here, X-rays are used as an example of radiation, and an X-ray image obtained by X-ray imaging is processed.
The configuration described below is an example configuration to which the present invention is applied, and is not limited thereto.
[0098]
In the first to fifth embodiments, the X-ray moving image in the case where a series of continuous images (X-ray moving images) corresponding to the movement of the subject is acquired by X-ray imaging using a grid is particularly simple. The problem that it is difficult to remove the grid stripe component by simple filtering processing and the problem that it is difficult to remove the grid stripe component at high speed have been solved by the following configuration.
[0099]
In other words, it is usually superimposed on each of a plurality of frame images constituting the X-ray moving image to be processed (target X-ray moving image) or due to the influence of sensor characteristics (characteristics such as saturation characteristics and defective pixels). By estimating and obtaining grid stripe component information (hereinafter also simply referred to as “grid stripe component”) that should exist as a stable stripe pattern over the entire image that cannot be detected by the filtering process of Extract fringe components.
[0100]
In addition, in the plurality of frame images constituting the target X-ray moving image, the grid stripe component does not substantially change in the adjacent frame images, so that the grid stripe component extraction process is not performed on all the plurality of frame images. Don't skip it.
[0101]
For example, when each frame image is a logarithmic image signal, the extracted grid stripe component is subtracted from each frame image. As a result, stable grid stripe information can be removed without affecting the X-ray moving image.
[0102]
Specifically, first, in the first embodiment, one frame image (target frame image) of an X-ray moving image obtained by X-ray imaging is analyzed. Then, for a pixel defect in a direction orthogonal to the grid stripe, a linear prediction type pixel defect correction process is performed. In addition, based on the analysis result, a grid stripe component is extracted from the target frame image and held as image information.
Here, the time required for the above analysis is longer than the time of continuously acquired frame images. Thus, until the analysis of the next frame image is completed, the previously extracted grid stripe component is used. Then, using the target frame image obtained by analyzing and extracting the grid stripe component as a start frame, the grid stripe component removal process is performed on each frame image until the next frame image to be analyzed. At this time, grid stripe components may be periodically extracted based on control from another timing control configuration corresponding to the time required for the analysis.
[0103]
Further, in the first embodiment, when extracting the grid stripe component from the target frame image, it is first extracted by low-order FIR filtering, and then envelope information is calculated by vector amplitude characteristic calculation with another FIR filtering. get. Then, based on the envelope information, the unsteady portion is extracted from the grid stripe component, and the unsteady portion is repaired using the surrounding steady portion, so that the grid stripe component as a whole Convert to a steady signal sequence. Further, in order to extract only the grid stripe component more appropriately, a filtering process for selectively extracting the spatial frequency corresponding to the grid stripe component is executed, and the result is set as the grid stripe component.
The grid stripe component acquired as described above is subtracted from each frame image. Thereby, the X-ray animation after the removal of the grid stripe component is obtained.
[0104]
In the second embodiment, pixel defect correction is performed on the frame image after removing the grid stripe component. As the pixel defect correction here, for example, pixel defect correction by averaging can be applied.
[0105]
In the third embodiment, detection means for detecting grid mounting is provided, and when the grid mounting is confirmed by the detection means, the frame image is analyzed, and based on this result, predictive pixel defect correction, and Remove grid stripe components.
[0106]
In the fourth embodiment, when the irradiation field of X-ray exposure is narrowed down so as to correspond to the imaging region of the subject, only the partial image corresponding to the irradiation field in the frame image is applied. The processing in the first embodiment is executed.
[0107]
In the fifth embodiment, there is a case where grid stripes do not exist for each part of the image depending on the state of the subject. Take.
Specifically, with respect to the extracted grid stripe component, the information of the portion where the grid stripe component does not exist is replaced with zero (“0”) data, so that the grid stripe removal process is not performed on the portion.
[0108]
In the first to fifth embodiments, for example, since the extracted grid stripe component exists as image information, the grid stripe component is removed from the frame image by holding the image. Even in such a case, the original frame image, that is, the frame image in a state where the grid stripes are not removed, can be restored using the stored image information of the grid stripe components.
[0109]
[Grid Stripe Component Removal Processing in First to Fifth Embodiments]
Here, the grid stripe component removal processing method used in the first to fifth embodiments will be described.
In the following description, the first to fifth embodiments are collectively referred to as “this embodiment”. The target frame image in the X-ray moving image is also referred to as “target image”.
[0110]
For example, a component substantially including the grid stripe component is separated from the frequency indicated by the grid stripe component, the separated component is processed into feature information that the grid stripe will indicate, and the processed information is converted into a grid. This is regarded as a fringe component and removed from the target image signal.
[0111]
The grid fringe component has a fairly strong component according to the spatial spectrum expression. If the spatial frequency of the grid stripe is appropriately selected, the Nyquist frequency at the time of sampling (spatial frequency corresponding to 1/2 of the sampling frequency) It can exist in the vicinity. As a result, such a state as shown in FIGS. 25A to 25D that does not overlap with the main component of a normal image signal can be easily obtained.
[0112]
Here, as described with reference to FIGS. 26A to 26D, only when a steep fluctuation component exists in the target image signal, it is difficult to separate the grid stripe component and the target image signal. It becomes.
[0113]
In some cases, the grid stripes themselves may not exist. This is because when a subject that almost completely blocks X-rays is taken, or when a strong X-ray that exceeds the dynamic range of the sensor reaches the sensor, the grid stripe component is generated by saturation. This is the case.
[0114]
Normally, when an X-ray image is acquired, since the X-ray dose that passes through the inside of the subject is emphasized, the X-ray dose is several hundred times as large as the subject transmission amount outside the subject that does not pass through. In general, it is meaningless to extend the dynamic range of a sensor or sensor amplifier to the outside of a subject without information. In most cases, the outside of the subject exhibits nonlinearity due to saturation and also has a grid stripe component. Or the contrast is lowered.
[0115]
Therefore, in the present embodiment, when there is a steep change in the target image that makes it difficult to separate the grid stripe component and the target image signal in the data in one line, and when the grid stripe itself does not exist. It was possible to adaptively cope with both, and to eliminate the grid stripe component without generating artifacts.
[0116]
Specifically, the grid stripe component removal processing in the present embodiment is mainly processing including the following first processing step to third processing step.
[0117]
First processing step:
Line data in a direction orthogonal to the grid stripe is sampled from the target image obtained by including the grid stripe component, and the spatial frequency of the grid stripe is determined.
[0118]
Second processing step:
The grid fringe component is sequentially extracted from the target image, and the result (grid fringe component) is subtracted from the target image. At this time, considering the occurrence of the artifact, even if the artifact occurs, the comparison has a smaller influence range. A component mainly composed of grid stripes is extracted by FIR filtering with a small span.
[0119]
Third processing step:
Based on the spatial frequency of the grid stripes obtained in the first processing step, the envelope of the component mainly composed of the grid stripes obtained in the second processing step is obtained, and the phase of the component is obtained by FIR filtering different from that component. It is obtained by calculating the vector amplitude with the component moved by 90 °.
[0120]
The grid stripe component removal processing including the first processing step to the third processing step will be described more specifically. First, the envelope information obtained in the third processing step always takes a positive value, The features include the following features (1) and (2).
(1) A very large value is shown for a steep fluctuation portion (for example, an edge portion).
(2) For a portion where there is no grid stripe, a small value that is substantially “0” is shown.
[0121]
In the present embodiment, based on the above-described envelope information, a value portion (very large value portion) indicated by the feature (1) and a value portion (very small value) indicated by the feature (2). By repairing the grid stripe component, the creation of a more stable grid stripe component is realized.
[0122]
As a method of repairing the grid stripe component, for example, there is a method of replacing the portion of the feature (1) with a component predicted from a stable grid stripe portion in the vicinity thereof to obtain a stable grid stripe component as a whole. Can be mentioned.
[0123]
A normal filtering process is performed on the entire line with respect to a component mainly composed of grid stripes obtained only as described above and including only a stable component. At this time, filtering is performed to extract only a narrower range of spatial frequencies centered on the spatial frequency of the grid stripes.
[0124]
The filtering result (extracted component) is the grid stripe component in the target image. When the envelope information has the characteristic (2), that is, the grid stripe component is mainly included in the grid stripe component. When there is a non-existing portion, the grid stripe does not exist originally, so the component mainly composed of the grid stripe of the portion is replaced with “0”.
[0125]
In addition, when performing a filtering process on a target image, a fast Fourier transform algorithm may be used in order to perform a steep filtering process stably and at high speed. In this case, the number of data points is “2” to the nth power. (N is a positive integer). For this reason, "0" is padded around the normal data so as to match the points. Data in a range filled with “0” may be considered as data corresponding to a portion where the envelope information indicates the feature (2).
[0126]
The effective spatial frequency of the grid itself is 30% to 40% of the sampling frequency (reciprocal of the spatial sampling pitch) proposed in Japanese Patent Application No. 2000-028161 and the like. It is effective to be selected from such a spatial frequency (60% to 80% of the Nyquist frequency). This is because the main component of the image is generally concentrated at 30% or less of the sampling frequency, and the component of strong grid stripes having a spatial frequency of 40% or more and 60% or less of the sampling frequency is an interpolation similar to linear interpolation after sampling. This is because when the processing is performed, it seems that another periodic amplitude fluctuation occurs, and the grid stripe itself is not stable.
[0127]
When the spatial frequency of the grid itself is “fg [cyc / mm]” and the sampling pitch of the sensor is “T”, the spatial frequency fm of the grid stripe is
[0128]
[Expression 1]

[0129]
It is represented by the following formula (1).
[0130]
In the present embodiment, in consideration of the fact that a striped pattern corresponding to the spatial frequency fm represented by the above formula (1) exists in the target image as a grid striped component, the grid striped in the first processing step described above. Is extracted accurately. That is, since the spatial frequency fg of the grid stripes is known in advance, the vicinity of the spatial frequency fm of the target image is searched from the sampled line data, and the spatial frequency indicating the peak value is obtained from the search result. Thus, the spatial frequency fm of the grid stripes in the target image is considered.
[0131]
Further, in the second processing step described above, the FIR filtering of the smallest possible span with the spatial frequency fm as the center is performed on the spatial frequency fm, so that the effective image components are almost removed and the steepness is sharp. The grid stripe component is roughly extracted in a state where the influence of the artifact due to the fluctuation (for example, the edge portion) is within a narrow range.
[0132]
At this time, in order to eliminate the phase fluctuation, for example, the coefficient series of the FIR filter is an even function, and a three-point or five-point FIR filter is desirable to satisfy the above narrow range.
[0133]
Specifically, for example, when a FIR filter having three symmetrical points is used and the coefficients are (a1, b1, a1), the response at the spatial frequency fm is “1” in order to obtain the coefficients (a1, b1, a1). Two conditions can be used: a condition of “0” and a condition that the DC component, which is the center value of the image information, is set to “0”.
[0134]
That is, the coefficient calculation is
2 ^* a1 + b1 = 0
2 ^* a1 ^* cos (2πfmT) + b1 = 1
It is a simultaneous equation represented by the formula
[0135]
[Expression 2]

[0136]
The solution represented by the following formula (2) is taken.
[0137]
In the FIR filtering described above, the response is “1” at the spatial frequency fm, but the response gradually increases when the spatial frequency becomes higher. In general, there is no image component in this portion, so that the grid stripes can be sufficiently extracted even with the FIR filtering.
[0138]
In addition, when FIR filtering is performed with five symmetrical points and the coefficients are (a2, b2, c2, b2, a2), the coefficients (a2, b2, c2, b2, a2) are obtained at the spatial frequency fm. In addition to the two conditions of the condition that the response is “1” and the condition that the DC component that is the center value of the image information is “0”, the differential value of the response at the spatial frequency fm is “0” (peak ) Can also be used.
[0139]
That is, the coefficient calculation is performed by performing a simple calculation from the solution represented by the above formula (2).
(-A1 ² , 2a1 (1-b1), 1-2a1 ² -(1-b1) ² ,
2a1 (1-b1), -a1 ² )
The following solution is obtained.
[0140]
As a method of obtaining a filter for symmetric five-point FIR filtering, for example, first, when a filter having coefficients (a1, b1, a1) represented by the above formula (2) is subtracted from “1”, The filter has a zero point at a space fraction fm. Considering the process of performing the filtering by this filter twice, it still has a zero at the spatial frequency fm, but the phase (sign) is not inverted. Such a filter is a symmetric 5-point FIR filtering filter, and a filter having a peak at the target spatial frequency fm can be configured by subtracting the filter from “1”.
[0141]
FIG. 1 shows the shapes of symmetric 3-point FIR filtering (hereinafter also referred to as “3-point FIR filtering”) and symmetric 5-point FIR filtering (hereinafter also referred to as “5-point FIR filtering”) as described above. An example of spatial frequency characteristics) is shown.
[0142]
In most cases, the result of filtering by the FIR filter shown in FIG. 1 is a result of extracting only the grid stripe component. This is because, as is clear from FIG. 1, most of the low-frequency components of the effective image components mainly composed of the low-frequency components are removed.
[0143]
However, as described above, it is a fact that a significant amount of effective image components are included in the components extracted by FIR filtering. Originally, we would like to perform filtering with a filter having a steep selection characteristic centered on the spatial frequency fm, but even if this is done, it will contain frequency components that make up the rapidly changing part included in the target image. It will not change.
[0144]
Therefore, in order to solve the above problem, in the present embodiment, in the third step described above, the local envelope of the grid stripe component is obtained, and artifacts other than the grid stripe component are generated from the fluctuation. By detecting a part including a possible component, only the grid stripe component is stably extracted (created).
[0145]
A general signal envelope must be based on the Hilbert transform, but a single sinusoidal envelope has a response amplitude of “1” at that frequency and a phase variation of 90 ° (π / 2). It can be obtained by applying a spatial filter that causes the occurrence and taking the vector amplitude (square root of the sum of squares) of the result and the original signal.
[0146]
When filtering with discrete data using a FIR filter whose phase varies by 90 °, the coefficients of the FIR filter are point-symmetric (odd functions). For example, when this coefficient is (−a3, 0, a3), in order to set the response to “1” at the spatial frequency fm, the coefficient (−a3, 0, a3) is
2 ^* a3 ^* sin (2πfmT) = 1
Must satisfy the formula
[0147]
[Equation 3]

[0148]
The solution represented by the following formula (3) is obtained.
[0149]
The amplitude of the signal sequence obtained by the FIR filter having the solution coefficient (−a3, 0, a3) represented by the above equation (3) and the original signal sequence is obtained.
[0150]
For example, FIG. 2A shows the result of filtering the original image signal shown in FIG. 25A with the solution coefficients (a1, b1, a1) shown in the equation (2). Is shown. As apparent from FIG. 2A, most of the grid stripe components are extracted.
[0151]
FIG. 3A shows the result of filtering the original image signal shown in FIG. 26A with the solution coefficients (a1, b1, a1) shown in the above equation (2). Is shown.
[0152]
The waveforms shown in FIG. 2 (b) and FIG. 3 (b) (the waveforms shown by bold lines) are the above-described equations (3) for the original image signal shown in FIG. 2 (a) and FIG. The envelope obtained by taking the square root of the square sum of the result of filtering having the solution coefficient (−a3, 0, a3) and the original image signal is shown.
[0153]
In particular, in the envelope shown in FIG. 3B, when attention is paid to the depression shown in FIG. 3C, an unsteady component is clearly present in the depression. This is because the extracted grid stripe component is abnormally extracted by simple filtering (including the edge component of the target image), and if subtracted from the target image signal, an artifact occurs in the target image signal after this processing. Means that.
[0154]
Therefore, in the present embodiment, a range indicating an abnormal numerical value is identified from the envelope obtained as described above, and the grid stripe component of the range is corrected (replaced) with an estimated value from the surrounding numerical sequence. ) In other words, the grid stripe component is formed (created) by utilizing the characteristic of having a steady component over the entire range, which is a characteristic of the inherent grid stripe component.
[0155]
The estimated value (predicted value) used in the correction is obtained from the statistical properties of the data around the range showing abnormal numerical values. For example, since the spatial frequency fm of the grid stripe component is known, this spatial frequency fm can be used as a statistical property.
[0156]
For example, with the spatial frequency fm of the grid stripe and the phase φ,
[0157]
[Expression 4]

[0158]
The grid stripe component of the unsteady part is formed using the sine wave expressed by the following equation (4).
[0159]
For example, as the simplest method, there is a method in which two coefficients A and φ at a specific frequency are obtained from peripheral pixels using Fourier transform (Fourier series expansion).
[0160]
However, normal Fourier transform cannot be used due to problems such as the presence of defects (unsteady portions) in the data. Therefore, here, Fourier transform is generalized, and amplitude and phase information is obtained in the sense of least squares. Therefore, the above equation (4) is
[0161]
[Equation 5]

[0162]
The following equation (5) is transformed.
[0163]
In this case, the square error ε when the data at the sampling point xi is “yi” (number of data points n) ({xi, yi; i = 0 to n−1}) is
[0164]
[Formula 6]

[0165]
It is represented by the following formula (6).
[0166]
At this time, it should be noted that as the component “xi, yi” used here, only data determined to be a stationary part from the above-described verification of the envelope component is selected. Then, parameters R and I that minimize the square error ε are obtained as follows.
[0167]
First,
[0168]
[Expression 7]

[0169]
Equation (7)
[0170]
[Equation 8]

[0171]
This is expressed by the following equation (7 ′).
By solving the simultaneous equations of the above equation (7 ′), the parameters R and I can be obtained, and the phase φ and the amplitude A can be estimated simultaneously.
[0172]
Here, if the data series divides the section of k / (2 · fm) (k is a positive integer) equally into m, the above equation (7 ′) is
[0173]
[Equation 9]

[0174]
As shown by the following equation (8), a discrete Fourier transform (Fourier series expansion) for obtaining a coefficient of a specific frequency is obtained.
[0175]
The unsteady state removed as inappropriate by obtaining the values of the parameters R and I by using the appropriate data around the unsteady portion according to the above formula (7 ') or the above formula (8). Repair (replace) the part.
[0176]
As another repair method, there is a method in which a linear prediction model is considered and repair is performed by predicting sequentially using a linear prediction algorithm without specifying the spatial frequency of grid stripes.
[0177]
The signal waveform obtained by the above-described repair is generally a steady sine wave, and is a component that very well represents the grid stripe component.
[0178]
However, the above signal waveform (grid fringe component signal waveform) was originally the result of narrow span FIR filtering with the coefficients (a1, b1, a1) shown in equation (2) above, and is shown in FIG. It has such a response characteristic of the filter as described above and contains a large amount of image components other than the grid stripe components.
[0179]
Therefore, in the present embodiment, filtering for extracting only components near the spatial frequency fm of the grid stripes is performed on the signal waveform. This filtering is performed on a component in which an unsteady component has already been repaired by the above-described operation, and artifacts such as ringing due to the filtering do not occur.
[0180]
When the envelope generation described above is performed on the signal of the grid fringe component after filtering, when there is a portion where a very small value (a value close to “0”) is observed, this portion is: Originally, the grid stripe component is a portion that has not been observed for some reason (for example, X-rays are completely blocked or the sensor is saturated), and the grid stripe component does not originally exist. Therefore, information on this part is recorded and replaced with “0” after subsequent filtering. The result of this processing is subtracted from the target image signal as a grid stripe component.
[0181]
In the present embodiment, grid stripe component extraction processing is executed while sequentially extracting target line data from the target image signal one line at a time. When one line data is extracted, the average of several line data before and after that is extracted. It is also possible to extract the grid stripe component after obtaining the image, that is, after weakening the image component or enhancing the grid stripe component.
[0182]
The above is because the grid stripe and the sensor are usually arranged almost in parallel, and the components of the grid stripe of a certain line and the grid stripe of a line in the vicinity thereof are very similar.
[0183]
Therefore, as another form of the present embodiment, in order to reduce the number of calculation processes for extracting grid stripe components, the processing line is thinned out by using the above-mentioned very similar feature. The grid stripe component extracted in step 1 is set as the grid stripe component of the line in the vicinity thereof. In other words, the grid fringe component extraction process is not performed on the line in the vicinity of the line from which the grid fringe component is extracted, and the subtraction process from the line data in the vicinity is performed using the grid fringe component extracted in the certain line. Can also be performed.
[0184]
In order to determine whether or not the above-described thinning is possible, in the first processing step described above, when measuring the spatial frequency of the grid stripes, the lines before and after being sampled or between the sampled lines are measured. Check the phase difference of the grid stripe and confirm that the grid stripe is not tilted with respect to the sensor.
[0185]
As another form of the present embodiment, as described above, when the grid stripe component is removed from the target image signal, the irradiation field in the target image signal is recognized by detecting the signal intensity or the like. The above-described grid stripe component removal process may be performed only on the image data inside the irradiation field.
[0186]
Moreover, if the extraction of the grid stripe component used for the grid stripe component removal processing as described above is performed on all the frame images constituting the X-ray moving image, a lot of time is required. Since it is impossible to execute the extraction process in real time, in the present embodiment, the extraction process is performed on a selected frame image (target frame image) among a plurality of frame images constituting the X-ray moving image. The grid stripe component acquired by this is held as an image, and a grid stripe component removal process using the same grid stripe component is performed on a frame image at a time close to the target frame image.
[0187]
By the way, for example, when acquiring an X-ray image with a solid-state image sensor, there is a problem of defective pixels as a problem peculiar to the solid-state image sensor formed by arranging a plurality of pixels. Due to the redundancy of image information (spatially having a low frequency component as a main component), if there are a small number of defective pixels, in most cases, the defective pixels are repaired by average interpolation from the surrounding pixel values ( Correction) is possible.
[0188]
However, in general, prediction is necessary due to the statistical properties around defective pixels. For example, when the grid stripe is 50% or more of the Nyquist frequency as in the present embodiment, the prediction is reversed in the average interpolation.
[0189]
FIG. 4 shows a response function as filtering in a case where an arbitrary point is a defective pixel in one dimension and interpolation is performed with the average of two pixels on both sides thereof.
In FIG. 4 above, the spatial frequency is represented on the horizontal axis. As shown in FIG. 4, when the spatial frequency is low and is 50% or less of the Nyquist frequency, the response is “positive”, that is, the phase is not reversed. On the other hand, when it becomes 50% or more of the Nyquist frequency, the phase is inverted, and an expected interpolation result cannot be obtained.
[0190]
FIG. 5 shows an example of defective pixel interpolation.
In FIG. 5, black circle points represent pixel values obtained from normal pixels, and points indicated by arrows (“defective pixel positions”) represent defective pixels for which data is not obtained. Further, FIG. 5 shows a state in which each pixel data (data indicated by a black dot) vibrates finely due to reflection of grid stripes.
Also, the white circle point “A: average interpolation value” in FIG. 5 is a pixel value obtained by the conventional average interpolation, and “B: ideal interpolation value” in FIG. 5 is indicated. A certain white circle point is an interpolation value in consideration of grid stripes.
[0191]
In the present embodiment, the following two methods are performed in order to obtain an ideal interpolation value indicated by “B: ideal interpolation value” in FIG.
(Method 1: Linear prediction method)
In FIG. 5, the data of the defective pixel position is obtained by linear prediction from the surrounding pixels.
(Method 2)
By removing the grid stripe component from the original image signal by the above-described grid stripe removal process, the interpolation process by averaging which is a conventional method is performed in a state where the main component of 50% or less of the Nyquist frequency is eliminated.
[0192]
The outline of the linear prediction method (method 1) will be described below.
[0193]
First, as image data (pixel data) to be processed,
Data series {X _n , X _n-1 , X _n-2 , ..., X _np }
And data X in “n” _n But,
[0194]
[Expression 10]

[0195]
It is assumed that the difference equation (9)
In the above formula (9), “ε _n “Represents a white noise sequence, and“ ai {i = 1,..., P} ”represents a linear prediction coefficient. Such a sequence is expressed as“ autoregressive process (AR process) X ”. _n "
[0196]
The above equation (9) is transformed into the delay operator Z ^-1 If you rewrite using
[0197]
[Expression 11]

[0198]
Equation (10) is obtained.
[0199]
However, the above formula (10) is
[0200]
[Expression 12]

[0201]
Since it is expressed by the following equation (10 ′), the AR process X _n Is the pulse transfer function 1 / A (z _-1 ) Input of a linear filter with _n Can be defined (spectrum estimation).
[0202]
In addition, the above equation (9) is obtained by calculating the linear prediction coefficient ai {i = 1,..., P} from the reliable (n−1) -th pixel data. This indicates that the pixel data can be predicted.
[0203]
The prediction of the linear prediction coefficient ai {i = 1,..., P} can be performed by using maximum likelihood estimation (least square estimation) while assuming that the apparatus or system is stationary. That is, ε _n Can be obtained that minimizes the power (dispersion). ε _n Is an error obtained by least square estimation, and does not have the following correlation component corresponding to the predicted order. Therefore, the prediction error ε obtained by prediction with the necessary and sufficient order p _n Becomes white noise as defined in Equation 9.
[0204]
Prediction error ε _n Is a mean square (mean value “0”), the function E [*] representing the mean is
[0205]
[Formula 13]

[0206]
It is represented by the following formula (11).
[0207]
In the above formula (11),
R (τ) = E [X _m X _{m +} τ] (covariance function: covariance)
And both sides have a coefficient a to obtain the minimum value. _k If you differentiate it and set it to “0”,
[0208]
[Expression 14]

[0209]
The following simultaneous equations (12) are obtained. This is what is called a normal equation or Yule-Walker equation.
[0210]
Actually, the autocorrelation R (*) is not calculated at all pixel points, but an estimated value calculated with a limited (given) number of pixels is used.
For example, the Levinson algorithm which is a high-speed algorithm is used, but the Burg algorithm may be used. The Burg algorithm is an algorithm based on the maximum entropy method that can be obtained without directly calculating covariance (autocorrelation) with data having a smaller number of pixels. These algorithms mathematically match if the prediction error is a normal distribution and the number of pixels is large, but if the number of pixels is small, the Burg algorithm (an algorithm based on the maximum entropy method) is advantageous.
[0211]
The coefficient a obtained by the above calculation _k Is used to obtain expected pixel data from pixels before and after the defective pixel.
[0212]
[First embodiment]
The present invention is applied to, for example, an X-ray image acquisition apparatus 100 as shown in FIG.
[0213]
<Overall Configuration and Operation of X-ray Image Acquisition Apparatus 100>
The X-ray image acquisition apparatus 100 according to the present embodiment is an apparatus for acquiring a medical (image diagnosis or the like) X-ray moving image. As shown in FIG. 1, the X-ray image is acquired from a subject 102 (here, a human body). ), An X-ray generation unit 101 for generating scattered X-rays from the subject 102, a planar X-ray sensor 104 for detecting the distribution of X-ray dose transmitted through the subject 102, A controller 105 (CONT) of the X-ray generator 101, an analog / digital (A / D) converter 106 that converts an electric signal output from the X-ray sensor 104 into digital data, and an output from the A / D converter 106 A memory 107 for temporarily storing the digital data to be processed as frame image data, and a memory 10 for storing the digital data acquired by imaging without emitting X-rays And an arithmetic unit 109 that performs arithmetic processing using the data in the memory 108 on the frame image data in the memory 107, and a conversion table of frame image data after the arithmetic processing in the arithmetic unit 109 (reference table: Look UpTable, (Hereinafter also referred to as “LUT”) 110, a memory 111 for storing gain pattern data for correcting gain variations among pixels constituting the X-ray sensor 104, and converted frame image data output from the LUT 110. On the other hand, an arithmetic unit 112 that performs arithmetic processing using the gain pattern data of the memory 111, a memory 113 that temporarily stores frame image data after the arithmetic result of the arithmetic unit 112, and a defective pixel unique to the X-ray sensor 104 Memory 114 for storing information (defective pixel position information, etc.) and information in the memory 114 A correction processing unit 115 that performs correction processing on the frame image data stored in the memory 113, a memory transfer control unit 142 that controls reading of the frame image data after the correction processing in the memory 113, and a memory transfer control unit A memory 143 that stores frame image data read from the memory 113 under the control of 142, a grid stripe detection unit 116 that detects information about grid stripes from the frame image data in the memory 143, and a grid stripe detection unit 116. The grid stripe component extraction unit 117 that extracts the grid stripe component from the frame image data in the memory 143 based on the information obtained in the above, and the memory that temporarily stores the grid stripe component extracted by the grid stripe component extraction unit 117 118 and a memory transfer that controls reading of grid stripe components in the memory 118. The control unit 140, the memory 141 that stores the grid stripe component read from the memory 118 under the control of the memory transfer control unit 140, and the calculator 119 that subtracts the grid stripe component in the memory 141 from the frame image data in the memory 113. And a memory 120 that temporarily stores a calculation result (frame image data after removal of the grid stripe component) in the calculator 119, and an image processing unit 121 that performs image processing on the frame image data in the memory 120 and outputs the result. ing.
[0214]
In the X-ray image acquisition apparatus 100 as described above, the controller 105 of the X-ray generation unit 101 starts X-ray emission at the X-ray generation unit 101 when a generation trigger is applied by an operator from an operation unit (not shown). To do.
Thereby, the X-ray generation unit 101 emits X-rays to the subject 102 that is a human body.
[0215]
X-rays radiated from the X-ray generation unit 101 pass through the subject 102 and reach the X-ray sensor 104 via the grid 103 that removes scattered X-rays from the subject 102.
[0216]
The X-ray sensor 104 has a configuration in which a plurality of detectors (pixels) for detecting the X-ray intensity are arranged in a matrix on a surface (image receiving surface) for detecting the distribution of the X-ray dose transmitted through the subject 102. An electric signal corresponding to the X-ray intensity obtained by the plurality of detectors (pixels) arranged in a matrix is output.
[0217]
As the X-ray sensor 104, for example, the following sensors (1) and (2) are applicable.
Sensor (1):
A sensor that converts X-ray intensity into fluorescence once, and detects the fluorescence by photoelectric conversion with a plurality of detectors arranged in a matrix.
Sensor (2):
X-rays radiated to a specific object are photoelectrically converted in the object and free electrons are attracted by a uniform electric field to form a charge distribution, and the charge distribution is divided into a plurality of matrix-arranged charges. A sensor adapted to be converted into an electrical signal by a charge detector (capacitor).
[0218]
The A / D converter 106 digitizes and outputs the electrical signal output from the X-ray sensor 104.
Specifically, the A / D converter 106 sequentially converts the electric signal output from the X-ray sensor 104 into digital data in synchronization with the X-ray emission of the X-ray generation unit 101 or the driving of the X-ray sensor 104. Convert to and output.
[0219]
In FIG. 6, one A / D converter 106 is provided. However, for example, a plurality of A / D converters may be provided and operated in parallel. As a result, the digital conversion speed can be increased, and the processing can proceed efficiently.
[0220]
The digital data output from the A / D converter 106 is temporarily stored in the memory 107 as frame image data.
Accordingly, the memory 107 stores digital image data (frame image data) that is a set of a plurality of pixel data corresponding to a plurality of pixels constituting the X-ray sensor 104.
[0221]
The memory 108 stores in advance digital data acquired by imaging without emitting X-rays. This digital data is data for removing fixed pattern noise existing in an offset manner unique to the X-ray sensor 104 from the frame image data stored in the memory 107. Therefore, in the X-ray image acquisition apparatus 100, imaging is performed in a state where the X-ray generation unit 101 does not emit X-rays, and the digital data acquired thereby is stored in the memory 108 as image data.
[0222]
The computing unit 109 stores image data stored in the memory 108 for each of a plurality of pieces of pixel data constituting frame image data (image data obtained by X-rays transmitted through the subject 102) stored in the memory 107. A process of subtracting pixel data at corresponding positions from among a plurality of pixel data constituting (fixed pattern noise image data obtained by imaging without X-rays) is executed.
[0223]
The LUT 110 converts the frame image data processed by the computing unit 109 into a value proportional to the logarithm and outputs it.
[0224]
The memory 111 stores gain pattern data for correcting variation in gain of each pixel constituting the X-ray sensor 104 with respect to the frame image data converted by the LUT 110. For this reason, in the X-ray image acquisition apparatus 100, X-ray imaging is performed in the absence of the subject 102, and fixed pattern noise is removed from the obtained image data using digital data stored in the memory 108. Further, the data obtained by converting into a value proportional to the logarithmic value by the LUT 110 is stored in the memory 111 as gain pattern data.
[0225]
The computing unit 112 subtracts the gain pattern data in the memory 111 from the frame image data output from the LUT 110 (corresponding to division if logarithmic conversion is not performed), and outputs the result.
The frame image data subjected to the subtraction processing by the calculator 112 is temporarily stored in the memory 113.
[0226]
Note that, when image data used for gain pattern data to be stored in the memory 111 is acquired, if image capturing is performed with the grid 103 attached, grid stripes are reflected in the gain pattern data obtained thereby. Become. What is expected is that when the gain pattern data is subtracted from the target image data by the computing unit 112, the grid stripes reflected in the subject 102 are close to gain fluctuations, so that the grid stripe components are removed. It is possible that However, the gain pattern data obtained by shooting without the subject 102 is unlikely to be acquired every time an image is acquired (every actual shooting), and in most cases, is acquired once a day or less frequently. In addition, since the positional relationship between the X-ray generation unit 101 and the X-ray sensor 104 may change for each imaging, the grid stripe component is not removed by the above subtraction. Even if the positional relationship is not changed, since the X-ray scattered dose and the radiation quality are generally different between the imaging with the subject 102 and the imaging without the subject 102, the grid stripe contrast is different. The fringe component is not removed by subtraction.
In any case, if the grid 103 direction matches, the spatial frequency of the grid stripes will not fluctuate. Preferably, when gain pattern data is acquired, the grid 103 itself should be removed so that grid stripes are not included in the gain pattern data.
[0227]
The memory 114 stores information regarding defective pixels unique to the X-ray sensor 104 (defective pixel position information and the like).
Specifically, for example, a flat X-ray sensor is generally manufactured by a semiconductor manufacturing technique, but the yield is not 100%, and a plurality of detectors (pixels) are not included due to some cause in the manufacturing process. Some are defective pixels that do not make sense as detectors, that is, their output has no meaning. Here, the X-ray sensor 104 is inspected in advance in the manufacturing process or by means (not shown), and the position information of the defective pixel obtained as a result is stored in the memory 114.
[0228]
The correction processing unit 115 corrects the defective pixel data in the plurality of pixel data constituting the frame image data stored in the memory 113 based on the position information of the defective pixel stored in the memory 114, and the corrected pixel The data is stored again in the corresponding position in the memory 113.
[0229]
The memory transfer control unit 142 controls reading of frame image data in the memory 113 (image data after correction processing by the correction processing unit 115) in the memory 113 based on a calculation end signal given from a grid stripe component extraction unit 117 described later. .
The memory 143 stores the frame image data read from the memory 113 by the memory transfer control unit 142.
[0230]
The grid stripe detection unit 116 analyzes grid stripes on the frame image data in the memory 143 (image data after the correction processing by the correction processing unit 115), and the grid stripe spatial frequency fm and the grid stripe angle. Detect and output θ.
[0231]
The grid stripe component extraction unit 117 reads frame image data (image data after the correction processing by the correction processing unit 115) in the memory 143, and the grid stripe space obtained by the grid stripe detection unit 116 from the read image data. Based on the frequency fm and the grid stripe angle θ, a grid stripe component is extracted, and then an operation end signal indicating the end of the extraction process is output to each of the memory transfer control unit 142 and the memory transfer control unit 140. .
The grid stripe component obtained by the grid stripe component extraction unit 117 is temporarily stored in the memory 118.
[0232]
The memory transfer control unit 140 controls reading of the grid stripe component in the memory 118 based on the calculation end signal given from the grid stripe component extraction unit 117.
The memory 141 stores the grid stripe component read from the memory 118 by the memory transfer control unit 140.
[0233]
The computing unit 119 subtracts the grid stripe component stored in the memory 141 from the frame image data in the memory 113 (image data after the correction processing by the correction processing unit 115).
The frame image data from which the grid stripe component has been subtracted by the calculator 119 is temporarily stored in the memory 120.
[0234]
The image processing unit 121 performs image processing on the frame image data in the memory 120 so that the observer can easily observe the frame image data.
Examples of the image processing here include the following processing.
-Random noise removal processing from frame images.
When the frame image is displayed, the gradation is converted or the detailed part is emphasized so that the density value is easy for the observer to see.
A part unnecessary for the observer is cut out from the frame image to reduce the information amount of the frame image or to compress the frame image information.
[0235]
The frame image data processed by the image processing unit 121 is displayed on a display unit, stored in a storage unit or a storage medium, and recorded on a recording medium by an unillustrated means, externally or in the X-ray image acquisition apparatus 100 Recording or analysis processing is performed.
[0236]
<Configuration and Operation Characteristic of X-ray Image Acquisition Apparatus 100>
7A to 7E illustrate the X-ray image acquisition apparatus 100 with particular attention paid to the timing of data reading by the memory transfer control unit 140 and the memory transfer control unit 142 and the processing for the read data. It is.
[0237]
FIG. 7A shows the timing at which the frame image data sequentially acquired by the X-ray sensor 104 is stored in the memory 113. The frame image data sequentially stored in the memory 113 is subjected to processing for removing grid stripe components in the memory 141 by the calculator 119.
[0238]
FIG. 7B shows the calculation end signal output from the grid stripe component extraction unit 117. As shown in FIG. 7B, the calculation end signal is a high-level signal and indicates that the extraction of the grid stripe component is ended (computation is completed).
[0239]
FIG. 7C shows the operation timing of the memory transfer control unit 142 based on the calculation end signal output from the grid stripe component extraction unit 117. As shown in FIG. 7C, the memory transfer control unit 142 has a timing at which the extraction of the grid stripe component is completed by the calculation end signal from the grid stripe component extraction unit 117 (the calculation end signal is When it is a high level signal), the frame image data in the memory 113 is read out to the memory 143.
[0240]
FIG. 7D shows the operation timing of the grid stripe detection unit 116 and the grid stripe component extraction unit 117. As shown in FIG. 7D, when the frame image data is transferred to the memory 143, the grid stripe detection unit 116 and the grid stripe component extraction unit 117 operate, and the grid stripe component obtained as a result is Stored in the memory 118. At the timing when this operation ends, the grid stripe component extraction unit 117 outputs a high-level calculation end signal as shown in FIG.
[0241]
FIG. 7E shows the operation timing of the memory transfer control unit 140 based on the calculation end signal output from the grid stripe component extraction unit 117. As shown in FIG. 7 (e), the memory transfer control unit 140 has a timing at which the extraction of the grid fringe component is completed by the calculation end signal from the grid stripe component extraction unit 117 (the calculation end signal is In the case of a high level signal), that is, when the grid stripe component is stored in the memory 118, the grid stripe component in the memory 118 is read out to the memory 141. Thereby, the grid stripe component in the memory 141 is updated.
[0242]
That is, in this embodiment, as shown in FIGS. 7A to 7E, grid fringe component extraction processing is not performed on all frame images constituting an X-ray moving image, but one frame. At the timing when the grid stripe component extraction processing for the image is completed, the next grid stripe component extraction processing is executed, and the same grid stripe component extracted immediately before (stored in the memory 141) is extracted for the frame image in the meantime. The grid stripe component is removed using the grid stripe component). In general, in such a frame image, the position of the X-ray generation unit 101, its X-ray energy, or the position of the subject 102 is almost unchanged in continuous frame images. This operation is effective when there is no fluctuation.
[0243]
In the present embodiment, the memory transfer control unit 142 is configured to capture the frame image to be subjected to the next extraction process into the memory 143 at the timing when the grid stripe component extraction process for one frame image is completed. However, the present invention is not limited to this. For example, a frame image may be periodically taken into the memory 143 using another arbitrary timing control mechanism, and the grid stripe component extraction process may be performed on the frame image.
[0244]
In this embodiment, X-ray moving images are targeted. However, for example, for a still image, a previously extracted grid stripe component is used for the grid stripe component removal processing for the still image acquired this time. May be. In this case as well, the burden required for extracting grid stripe components can be reduced.
[0245]
<Other Specific Configuration and Operation of X-ray Image Acquisition Apparatus 100>
Here, in the X-ray image acquisition apparatus 100 described above, the following components that are considered to require other specific description will be specifically described.
(1) Grid stripe component image data stored in the memory 118
(2) Correction processing of defective pixels by the correction processing unit 115
(3) Grid stripe component detection and extraction processing by the grid stripe detection unit 116 and the grid stripe component extraction unit 117
[0246]
(1) Grid stripe component image data stored in the memory 118
The image data of the grid stripe component stored in the memory 118 is data that is subtracted from the frame image data on which the grid stripe component is superimposed, but the data stored in the memory 118 is subtracted as in the present embodiment. If the configuration is such that it is stored separately in association with the subsequent target image data, the frame image data on which the original grid stripes are superimposed can be reproduced from the frame image data from which the grid stripes have been removed. Thereby, for example, even if the frame image data is damaged due to some trouble in the grid removal process, it is possible to restore the original frame image data by the above reproduction process.
[0247]
(2) Correction processing of defective pixels by the correction processing unit 115
The correction processing unit 115 executes processing as described below, for example, by software using a microprocessor.
[0248]
8 to 11 show examples of the distribution of pixel defects in the X-ray sensor 104. FIG.
Here, it is assumed that the pixel defect basically exists only in the width of one pixel. This is because an X-ray sensor having a plurality of adjacent pixel defects in a large lump is generally not used because it is difficult to repair defective pixels.
[0249]
8 to 11, each square represents a pixel, and the black square represents a defective pixel. Further, in the lower part of the figure, the direction (vertical direction) of the grid stripes of the grid 103 is illustrated.
[0250]
First, the defective pixel shown in FIG. 8 is a basic pixel defect. As shown in FIG. 8, there are eight adjacent pixel components a1 to a8 around the defective pixel (black cells). Yes.
[0251]
FIG. 12 schematically shows the signal distribution on the spatial frequency axis of the grid stripe component in the target image in which the defective pixel shown in FIG. 8 exists and the grid stripe component exists.
In FIG. 12, the horizontal axis represents the horizontal spatial frequency axis u of the target image, the vertical axis represents the vertical spatial frequency axis v of the target image, and both the spatial frequency axis u and the spatial frequency axis v. A “sampling frequency” that is the reciprocal of the pixel pitch and a “Nyquist frequency” that is a half value thereof are shown on the axis.
[0252]
Since the grid stripe vibrates in the horizontal direction of the target image and is constant in the vertical direction, the component of the grid stripe exists on the spatial frequency axis u as shown in FIG. (See white circle in the figure).
[0253]
In a normal image, the principal component is distributed in a spatial frequency region below the half value of the Nyquist frequency, and if there is no grid stripe component, interpolation can be performed by averaging pixel values on both sides of a defective pixel. . This is because the influence of the interpolation on the spatial spectrum indicates, for example, the filtering response function (characteristic) shown in FIG.
[0254]
Therefore, in the case of the correction of the defective pixel shown in FIG. 8, since there is no grid stripe component in the vertical direction, the average of the pixels in the vertical direction, that is, the average of the pixel component a2 and the pixel component a6, or any one of the pixels. Depending on the component, almost satisfactory correction is possible.
[0255]
The defective pixel shown in FIG. 9 is a pixel defect that is connected horizontally (see black squares in the figure). In the case of a defective pixel having such a shape as well, the vertical direction of each pixel is a direction parallel to the grid stripes as in the case of the defective pixel shown in FIG. Satisfactory correction is possible.
[0256]
The defective pixels shown in FIG. 10 are pixel defects that are vertically connected (see the black squares in the figure). In the case of a defective pixel having such a shape, there is no pixel having a reliable value above and below the target pixel with respect to the defective pixel at the upper end or the defective pixel at the lower end of the connected defective pixel. If defect correction is performed on a defective pixel in such a state with a simple horizontal average or the like, a correction result with an unexpected value as described with reference to FIG. 5 is obtained.
[0257]
Therefore, in the present embodiment, using the simultaneous equations as shown in the above equation (12) with normal pixel components connected to the right or left of the target defective pixel, the coefficient a _k (K = 1 to P) is obtained from the left and right of the target defective pixel.
At this time, for example, the number of pixels to be used is about 20, and the order k is about 5.
[0258]
And coefficient a _k Therefore, the target defective pixel value X is expressed by the above equation (9). _n All the defective pixel values X obtained as a result _n Find the average of. As a result, “B: ideal interpolation value” as shown in FIG. 5 is obtained.
[0259]
In the present embodiment, the coefficient a is calculated using the above equation (12). _k (K = 1 to P) is obtained, but the present invention is not limited to this. For example, an algorithm called a maximum entropy method may be used.
[0260]
The defective pixel shown in FIG. 11 is a defective pixel in which the state of the defective pixel shown in FIG. 9 and the state of the defective pixel shown in FIG. Among the defective pixels in this state, the pixel in question is a pixel at the intersection of the vertically connected defective pixel and the horizontally connected defective pixel (pixels overlapping the cross), that is, the pixel components a4 and a12. , A18, a24 are defective pixels.
[0261]
The correction of the defective pixels in the state as shown in FIG. 11 is performed by correcting the linear defective pixels continuous in the horizontal direction by averaging the upper and lower pixel components, and the correction of the linear defective pixels continuous in the vertical direction is performed. , Using the simultaneous equations as described above.
Specifically, the following three methods {circle around (1)} to {circle around (3)} can be mentioned, and almost the same result can be obtained by any method.
[0262]
(1) Correction of a defective pixel surrounded by the pixel components a4, a12, a18, a24 is performed using the average value of the upper and lower defect correction values (corrected defective pixel values).
(2) Correction of the defective pixel surrounded by the pixel components a4, a12, a18, and a24 is performed by solving the simultaneous equations of the above equation (12) using the left and right pixel values as the defect correction values.
(3) Correction of defective pixels surrounded by the pixel components a4, a12, a18, and a24 is performed using the average value of the results of (1) and (2).
[0263]
By executing the processing as described above, defective pixel data in a plurality of pixel data constituting frame image data stored in the memory 113 is corrected.
[0264]
(3) Grid stripe component detection and extraction processing by the grid stripe detection unit 116 and the grid stripe component extraction unit 117
[0265]
The grid stripe detection unit 116 reads a part of the frame image data stored in the memory 113, examines the spectrum of the grid stripe included in the frame image data based on the read data, and determines the spatial frequency fm and angle of the grid stripe. θ is detected. Information on the spatial frequency fm and the angle θ (hereinafter also referred to as “angle η”) of the grid stripes is used in the grid stripe component extraction process in the grid stripe component extraction unit 117 in the subsequent stage.
[0266]
FIGS. 13A to 13D are diagrams for explaining the detection processing of the spatial frequency fm and the angle θ of the grid stripes.
[0267]
FIG. 13A shows an image of the entire target frame image, and “L1” to “L6” represent line positions from the upper part of the target frame image. The grid stripe detection unit 116 measures the spatial frequency fm of the grid stripe based on the result of Fourier transform of the lines L1 to L6. At this time, when detecting the spectrum of the grid stripe, the grid stripe detection unit 116 may use the average of several lines before and after each line L1 to L6 (or the average of the spectrum) in order to increase the detection capability. Good.
[0268]
FIGS. 13B to 13D each show a Fourier transform result in the line L1.
That is, FIG. 12B shows the amplitude spectrum (or power spectrum), FIG. 12C shows the value of the real part that is the coefficient of the cosine wave of the Fourier transform result, and FIG. , The value of the imaginary part, which is the coefficient of the sine wave.
[0269]
FIG. 14 is a flowchart showing the above processing of the grid stripe detection unit 116.
[0270]
First, the grid stripe detection unit 116 clears a variable cumulation for obtaining an average of spectra (step S201).
In addition, the grid fringe detection unit 116 clears the target for obtaining the average of the spectrum and the counter (variable) n of the number of lines (step S202).
Further, the grid stripe detection unit 116 initializes a variable i for selecting a processing target line (selected line) from the lines L1 to L6 shown in FIG. 13A to “1” (step S203). Thereby, in the first process, the line L1 is selected and processed as the target line.
[0271]
Then, the grid stripe detection unit 116 determines whether or not the processing of the next step S205 to step S215 has been executed for all the lines L1 to L6 of the target frame image (step S204).
As a result of this determination, the process proceeds to step S216 described later only when the process is completed. If the process has not been completed yet, the process from the next step S205 is executed.
[0272]
If the processing is not completed as a result of the determination in step S204, first, the grid stripe detection unit 116 selects the line Li indicated by the variable i from the lines L1 to L6 of the target frame image, and the data (line data). Li) is acquired (step S205).
[0273]
Next, the grid stripe detection unit 116 performs a Fourier transform process such as a fast Fourier transform on the line data Li acquired in step S205 (step S206).
[0274]
Next, the grid stripe detection unit 116 obtains a power spectrum (or amplitude spectrum) from the Fourier transform result (spatial frequency domain data) in step S206 (step S207).
[0275]
Next, the grid stripe detection unit 116 determines whether or not there is a significant spectrum (peak value) indicating the grid stripe in the power spectrum acquired in step S207 (step S208).
[0276]
Specifically, for example, first, the absolute spatial frequency of grid lead that causes grid stripes is known at the stage where the grid 103 is installed, so that the frequency is used as “fg”. The discrimination process in step S208 can be performed accurately.
[0277]
That is, when the sampling pitch of the X-ray sensor 104 is “Ts”, the approximate spatial frequency fm where the grid stripes are generated is
[0278]
[Expression 15]

[0279]
It can be specified by the following conditional expression (13).
[0280]
At this time, in the above formula (13), when the condition shown in (3) is satisfied, the spatial frequency fm ′ obtained in (2) is used and the condition shown in (3) is not satisfied. In this case, (1) is executed as “J ← J + 1”.
[0281]
The accurate spatial frequency fm of the grid stripe should be present in the vicinity of “fm ′” obtained by the above equation (13), and it is determined whether or not a peak value (significant spectrum indicating the grid stripe) exists. In this case, if only the vicinity is searched, even if there are different peak values due to the influence of image components, noise components, etc., detection of peak values that are significant spectra showing grid stripes without being affected Can be done.
[0282]
Further, the frequency fg of the grid lead of the grid 103 is manufactured fairly accurately. However, if there is an arbitrary distance or more between the grid 103 and the X-ray sensor 104 during imaging, the X-ray generation unit Since the X-ray beam from 101 has a cone beam shape, the X-ray beam is expanded and reaches the X-ray sensor 104. For this reason, the exact spatial frequency fm becomes different and cannot be easily predicted.
Accordingly, the frequency indicating the peak value among the frequencies in the vicinity of “fm ′” obtained by the above equation (13) is obtained as the spatial frequency fm.
[0283]
However, in the above case, it is necessary to determine whether or not the value corresponding to the obtained peak value is a significant peak value that actually exists stably. This determination is normally made based on the noise level. The noise level may be measured in advance, or an average value of components other than the peak value in the high band of the spectrum may be substituted. For example, if the ratio between the sum (or average value) of the power spectra of neighboring spectrum values indicating peak values and the noise level is about 10 or more, it is determined that the peak value is usually significant and stable.
[0284]
As a result of the determination in step S208 as described above, if there is no significant peak value, the grid fringe detection unit 116 counts up the variable i to process the next line (step S217), and the step again Returning to S204, the subsequent processing steps are repeatedly executed.
[0285]
On the other hand, if there is a significant peak value as a result of the determination in step S208, the grid fringe detection unit 116 sets Pi as the spatial frequency indicating the peak value (step S209), and adds this to the variable accumulation (step S209). Step Ss210).
[0286]
Further, the grid stripe detection unit 116 obtains the phase of the grid stripe, sets the line position indicated by the phase and the variable i for the variable θn and the variable Mn (step S211 to step S214), and increments the variable n. (Step S215).
Thereafter, in order to process the next line, the grid stripe detection unit 116 counts up the variable i (step S217), returns to step S204 again, and repeatedly executes the subsequent processing steps.
[0287]
As described above, when the processing of steps S205 to S215 is completed for all the lines L1 to L6, the grid fringe detection unit 116 changes the value of the current variable cumulation to the value of the current variable n ( By dividing by the number of significant peak values), the average spatial frequency fm of the grid pattern is obtained (step S216).
[0288]
In addition, after executing the processing of FIG. 14, the grid stripe detection unit 116 uses the variable θi (i = 0 to n−1) obtained in step S213 as the line position Mi (i = 0 to n−) in step S214. This is used to obtain the average grid stripe angle η (angle θ) according to 1).
[0289]
That is, the grid stripe detection unit 116 determines the phase θ of the i-th line Li. _i And the phase θ of the (i + 1) th line L (i + 1) _{i + 1} And the phase difference (θ _i −θ _{i + 1} ) By the spatial frequency fm of the grid pattern,
[(Θ _i −θ _{i + 1} ) / 2π] / fm
The line difference between the line Li and the line L (i + 1) is
(Mi + 1) -Mi
Using these results, the angle η of the grid stripe is
[0290]
[Expression 16]

[0291]
It calculates | requires by Formula (14) which becomes.
[0292]
Then, when the angle η obtained by the above equation (14) is not abnormally inclined, the grid stripe detection unit 116 performs processing for extracting grid stripes every several lines, and the angle η When it is inclined abnormally, a process of extracting grid stripes is executed for each line.
[0293]
In the process executed by the grid stripe detection unit 116 described above, it is assumed that the direction (vertical or horizontal direction) of the grid stripe is known. For example, when the direction of the grid stripe is unknown, for example, The same processing is executed in advance in both the vertical direction and the horizontal direction, and the direction in which a significant peak value is detected is set as a direction substantially orthogonal to the grid stripes.
[0294]
When the grid stripe detection unit 116 executes the processing as described above, the spatial frequency fm and the angle θ (angle η) of the grid stripe are obtained.
The grid stripe component extraction unit 117 uses the grid stripe spatial frequency fm and the angle θ (angle η) obtained by the grid stripe detection unit 116, and includes the grid stripe component that is actually stored in the memory 143. A grid stripe component is extracted from the frame image data, and the grid stripe component is stored in the memory 118.
[0295]
FIG. 14 is a flowchart showing grid stripe component extraction processing in the grid stripe component extraction unit 117.
[0296]
First, the grid stripe component extraction unit 117 is given the grid stripe angle θ obtained by the grid stripe detection unit 116 and the spatial frequency fm of the grid stripe as processing parameters.
The grid fringe component extraction unit 117 performs the processing of step S300 to step S321 described below based on the processing parameters (grid fringe angle θ and grid fringe spatial frequency fm).
[0297]
If the value of the spatial frequency fm of the grid stripe given to the grid stripe component extraction unit 117 is “0”, the grid stripe does not exist on the target frame image, that is, the grid 103 is not used. In this case, the grid stripe component extraction unit 117 does not execute the processing shown in FIG. For example, the memory 118 stores “0” data as the grid stripe component in this case. Alternatively, since the computing unit 119 does not function, the target image data stored in the memory 113 is stored in the memory 120 as it is.
[0298]
When the grid stripe angle θ and the spatial frequency fm of the grid stripe are given from the grid stripe detection unit 116 to the grid stripe component extraction unit 117, first, the grid stripe component extraction unit 117 uses the above formula (2). Then, the coefficients a1 and a2 (or the coefficients (a2, b2, c2, b2, a2) of the symmetric 5-point filter) are obtained from the spatial frequency fm (step S300). The FIR filter corresponding to the coefficient obtained here is “FIR1”.
[0299]
Next, the grid fringe component extraction unit 117 calculates the coefficient a3 from the spatial frequency fm using the above equation (3) (step S301). The FIR filter corresponding to the coefficient obtained here is “FIR2”.
[0300]
Next, the grid fringe component extraction unit 117 generates a window function centered on the spatial frequency fm in order to perform FIR filtering in the region of the spatial frequency fm (step S302).
As the window function here, for example, a Gaussian distribution shape function centered on the spatial frequency fm can be applied.
[0301]
Next, the grid fringe component extraction unit 117 determines a range of lines on which the grid fringe extraction processing is performed on the line data constituting the target frame image based on the grid stripe angle θ (step S303 to step S303). S305).
[0302]
Specifically, for example, the grid stripe component extraction unit 117 sets the reference value of the angle θ to “0.1 degree”, and the grid stripe angle θ obtained by the grid stripe detection unit 116 is less than the reference value of 0.1 degree. Is also larger or smaller (step S303). As a result of this determination, when the angle θ is smaller than the reference value 0.1 degree, the grid stripe component extraction unit 117 sets “5” for the variable “kip”, thereby omitting the processing for every five lines. Is determined (step S304). On the other hand, when the angle θ is equal to or greater than the reference value 0.1 degree, the grid stripe component extraction unit 117 sets the variable skip to “0” to execute the process for all the lines. Determine (step S305).
[0303]
In step S303 to step S305, for example, not only the setting for the variable skip but also a finer setting may be made based on the angle θ.
[0304]
After the process of step S304 or step S305, the grid stripe component extraction unit 117 initializes the variable count by setting the value of the variable skip set in step S304 or step S305 for the variable count. Is performed (step S306).
[0305]
Then, the grid stripe component extraction unit 117 determines whether or not the current variable count value is equal to or greater than the variable kip value, that is, whether or not the grid stripe extraction process for the target line data should be executed. (Step S307).
If it is determined that the process is to be executed, the process proceeds to step S310. If not, the process proceeds to step S308. However, when this step S307 is first executed, it is always determined that the process is executed, so the process proceeds to the next step S310.
[0306]
As a result of the determination in step S307, in the case of processing execution (skip ≦ count), first, the grid stripe component extraction unit 117 reads one line data (target line) to be processed from the target frame image data stored in the memory 143. Data) is acquired (step S310).
[0307]
In step S310, the target line data may be acquired as it is from the target frame image data. For example, an average value (moving average value) for several lines before and after the target line data is obtained as an actual processing target. You may make it acquire as line data.
[0308]
Next, the grid fringe component extraction unit 117 performs filtering using the FIR1 whose coefficient is determined in step S300 on the target line data acquired in step S310, and stores it in the line buffer 1 (not shown) (step). S311).
As a result of this processing, the line buffer 1 stores image component data including grid stripe components.
[0309]
Next, the grid stripe component extraction unit 117 performs filtering using the FIR2 whose coefficient has been determined in step S301 on the line data stored in the line buffer 1, and stores the filtered data in the line buffer 2 (not shown) ( Step S312).
As a result of this processing, the line buffer 2 stores data for obtaining the envelope of the grid stripe component.
[0310]
Next, the grid stripe component extraction unit 117 obtains an envelope of the grid stripe component (step S313).
That is, the grid stripe component extraction unit 117 obtains the amplitude (that is, the square root of the sum of squares) of the vectors having the data in the line buffer 1 and the data in the line buffer 2 as components, and the result is obtained as the line buffer 3 ( (Not shown). As the calculation here, even a calculation that does not take the square root can be applied due to the monotonic increase of the square root, and the same effect can be obtained.
[0311]
Next, the grid fringe component extraction unit 117 examines data in the line buffer 3, that is, envelope data, and determines an upper limit value th1 and a lower limit value th2 for detecting the abnormal data ( Step S314).
Various methods can be applied to determine the upper limit value th1 and the lower limit value th2. For example, an average value and a standard deviation value are obtained, and n times the standard deviation value from the average value (n is, for example, “ 3) (a value about 3 ″) or a method in which the value shifted by more than the upper limit value th1 and the lower limit value th2 is obtained, or a histogram of envelope data in the line buffer 3 is obtained, and the upper limit value th1 and the lower limit value th2 are centered on the mode value. The method etc. of determining are mentioned.
[0312]
Next, in the envelope data in the line buffer 3, the grid fringe component extraction unit 117 regards data having a value greater than or equal to th1 or less than or equal to value th2 as abnormal data (data resulting from a sharp change in image data, etc.). Then, the data of the grid stripe component in the line buffer 1 corresponding to the abnormal data is estimated from the data around the abnormal data and rewritten (step S315). At this time, data rewriting is performed so that the entire grid stripe component is in a stable state showing a periodic variation pattern.
Note that the processing in step S315 has been described in the explanation part such as the above formula (7 '), and therefore the details thereof are omitted here.
[0313]
Next, the grid fringe component extraction unit 117 performs a Fourier transform process on the data of the grid fringe component in the line buffer 1 that has become an overall stable state by the process of step S315, and thereby the space of the grid fringe component. Frequency domain data is obtained (step S316).
Note that the transformation processing in step S316 is not limited to Fourier transformation, and other orthogonal transformations such as cosine transformation can be applied.
[0314]
Next, the grid fringe component extraction unit 117 performs filtering by the window function centered on the spatial frequency fm obtained in step S302 on the spatial frequency domain data acquired in step S316 (step S317). As a result, the data of the grid stripe component more selectively represents the grid stripe.
[0315]
Next, the grid stripe component extraction unit 117 performs the inverse conversion process of the conversion in step S316 on the grid stripe component data after the filtering process in step S317, and uses the result as the actual grid stripe component data. (Step S318).
[0316]
Next, the grid stripe component extraction unit 117 stores the data of the grid stripe component acquired in step S318 at the corresponding position in the memory 118 (step S319).
[0317]
Then, the grid stripe component extraction unit 117 sets “0” for the variable count (step S320), and has performed the processing from step S307 on all the line data constituting the target frame image data in the memory 143. It is determined whether or not (step S321).
As a result of the determination, when the process is completed, the grid stripe component extraction unit 117 outputs the above-described calculation end signal and ends the process. On the other hand, if the processing has not yet been completed, the process returns to step S307 again, and the subsequent processing steps are repeatedly executed.
[0318]
As a result of the determination in step S307 described above, when the processing is not executed (skip> count), the grid stripe component extraction unit 117 copies the data of the grid stripe component acquired in the previous stage to the corresponding position in the memory 118 (step In step S308, the variable count indicating the line to be copied is incremented (step S309), the process returns to step S307, and the subsequent processing steps are repeatedly executed.
[0319]
[Second Embodiment]
The present invention is applied to, for example, an X-ray image acquisition apparatus 400 as shown in FIG.
The X-ray image acquisition apparatus 400 of the present embodiment is different from the X-ray image acquisition apparatus 100 shown in FIG.
[0320]
In the X-ray image acquisition apparatus 400 of FIG. 16, the same reference numerals are given to components that operate in the same manner as the X-ray image acquisition apparatus 100 of FIG. 6, and detailed description thereof will be omitted.
[0321]
In the X-ray image acquisition apparatus 100 illustrated in FIG. 6, the correction processing unit 115 is configured to perform correction processing of defective pixels using the data in the memory 114 on the frame image data in the memory 113. . On the other hand, in the X-ray image acquisition apparatus 400 of the present embodiment, as shown in FIG. 16, the correction processing unit 115 performs frame image data in the memory 120, that is, a frame image after removal of grid stripe components. The correction processing of the defective pixel using the data in the memory 114 is performed on the data.
[0322]
Therefore, according to the X-ray image acquisition apparatus 400 of the present embodiment, pixel defect correction in consideration of grid stripes is not necessary, and therefore, using an average value of pixel values that are not peripheral defects as is conventionally known. Simple defective pixel correction such as correcting defective pixels can be applied.
[0323]
[Third embodiment]
The present invention is applied to, for example, an X-ray image acquisition apparatus 500 as shown in FIG.
The X-ray image acquisition apparatus 500 of the present embodiment is different from the X-ray image acquisition apparatus 100 shown in FIG.
[0324]
In the X-ray image acquisition apparatus 500 of FIG. 17, the same reference numerals are given to the components that operate in the same manner as the X-ray image acquisition apparatus 100 of FIG. 6, and detailed description thereof will be omitted.
[0325]
As shown in FIG. 17, the X-ray image acquisition apparatus 500 of the present embodiment further includes a detector (switch) that detects the mounting of the grid 103 in addition to the configuration of the X-ray image acquisition apparatus 100 of FIG. ) 122 is provided.
[0326]
The detection unit 122 supplies the detection result (grid mounting signal) of the mounting of the grid 103 to the correction processing unit 115 and the grid stripe detection unit 116, respectively.
[0327]
When the grid 103 is mounted according to the grid mounting signal from the detection unit 122, the correction processing unit 115 executes a defective pixel correction process in consideration of grid stripes as described in the first embodiment. If not, the defective pixel is corrected by the average value of the pixel values of the peripheral non-defective pixels.
[0328]
Similarly, when the grid 103 is attached by the grid attachment signal from the detection unit 122, the grid stripe detection unit 116 performs the grid stripe detection process (analysis process) as described in the first embodiment. Execute. However, when the grid 103 is not attached according to the grid attachment signal, the grid stripe detection unit 116 immediately determines that there is no grid stripe without executing the grid stripe detection process, and performs processing corresponding to this. Execute.
[0329]
As described above, in the X-ray image acquisition apparatus 500 according to the present embodiment, the detection unit 122 is provided, and the grid stripe is detected based on the detection result. The time required can be greatly reduced.
[0330]
[Fourth embodiment]
The present invention is applied to, for example, an X-ray image acquisition apparatus 600 as shown in FIG.
The X-ray image acquisition apparatus 600 of the present embodiment is different from the X-ray image acquisition apparatus 100 shown in FIG.
[0331]
In the X-ray image acquisition apparatus 600 of FIG. 18, the same reference numerals are given to the components that operate in the same manner as the X-ray image acquisition apparatus 100 of FIG. 6, and detailed description thereof will be omitted.
[0332]
As shown in FIG. 18, the X-ray image acquisition apparatus 600 of the present embodiment further includes a memory for storing X-ray irradiation area data in addition to the configuration of the X-ray image acquisition apparatus 100 of FIG. 123 is provided.
[0333]
The memory 123 stores image data (irradiation area data) obtained by cutting out only the area irradiated with X-rays from the frame image data stored in the memory 113. The irradiation area data in the memory 123 is stored in the memory 123. The grid stripe component is detected and extracted.
[0334]
Specifically, first, in X-ray imaging, in general, the X-ray generation tube 101 of the X-ray generation unit 101 is irradiated to the exit of the X-ray generation unit 101 in order to avoid exposure of the subject 102 (here, the human body) to a part other than the target part. A field stop is provided. Thereby, only the necessary part of the subject 102 can be irradiated with X-rays.
[0335]
When the above-mentioned irradiation field stop function is used, all image signals obtained from the X-ray sensor 104 are not effective for images obtained by X-ray imaging, and only a partial image corresponding to the X-ray irradiation field by the irradiation field stop is used. Becomes effective.
[0336]
Therefore, in the present embodiment, an unillustrated computer means (CPU or the like) responds to the X-ray irradiation field from the frame image data in the memory 113 based on the X-ray intensity distribution, shape, or other information. An effective partial image area (irradiation area) is found, and only the data of the irradiation area part (irradiation area data) is stored in the memory 123.
[0337]
As described above, in the present embodiment, not only the irradiation area data in the memory 123, that is, not all the target image data, but only the necessary part data with a reduced amount of information, so the processing time can be shortened. it can.
[0338]
In this embodiment, the irradiation area is cut out from the frame image after the defective pixel correction. For example, after the irradiation area is cut out, the defective pixel correction is performed on the irradiation area. Also good.
[0339]
[Fifth embodiment]
In the first embodiment, in the X-ray image acquisition apparatus 100 of FIG. 6 described above, the grid stripe component extraction processing in the grid stripe component extraction unit 117 is performed according to the flowchart of FIG.
In the present embodiment, the grid stripe component extraction process in the grid stripe component extraction unit 117 is, for example, a process according to the flowchart shown in FIG.
[0340]
Note that, in the grid stripe component extraction process of FIG. 19, steps that are the same as the grid stripe component extraction process of FIG. 15 are denoted by the same reference numerals, and detailed description thereof is omitted.
[0341]
First, before the description of the grid stripe component extraction processing in the present embodiment, FIG. 20A shows a part of a frame image including a grid stripe component. In FIG. The portion shown indicates a portion where the grid stripe component does not exist due to the presence of a substance that blocks X-rays.
[0342]
FIG. 20B shows the result of extracting the grid stripe component from the frame image of FIG. 20A by the filtering in the first embodiment. As shown in FIG. 20B, in the grid stripe component, an artifact appears in a portion corresponding to the portion indicated by “*” in the frame image. Therefore, as described in the first embodiment, the portion is It is extracted from the envelope and inferred from the data before and after it so that the whole becomes a stable grid stripe. This result is shown in FIG. 20 (c). With such a stable grid stripe component, no new artifact is generated by a transformation process such as Fourier transform.
[0343]
In the first embodiment, the grid stripe component as shown in FIG. 20C is subtracted from the original frame image as shown in FIG. In the case of this configuration, it is conceivable that a new grid stripe component appears in a portion where no grid stripe component originally exists (portion indicated by “*”).
[0344]
Therefore, in the present embodiment, as shown in FIG. 20 (d), in the grid stripe component as shown in FIG. 20 (c), the portion where the grid stripe component originally does not exist (the portion indicated by “*”). Is set to “0”.
[0345]
For this reason, in this Embodiment, the grid stripe component extraction part 117 performs the grid stripe component extraction process according to the flowchart of the said FIG.
That is, the grid fringe component extraction unit 117 replaces the portion of the stable grid fringe component data obtained by this process after estimation by step S315 with “0” after performing the process of step S318. (Step S700), and then proceeds to the next step S319. Thereby, it is possible to reliably prevent a new grid stripe component from being generated by the grid stripe removal process even in a portion where the grid stripe component does not exist originally.
[0346]
In the first to fifth embodiments, hardware is used. However, the entire apparatus can be realized by controlling the software.
[0347]
Another object of the present invention is to supply a storage medium storing software program codes for realizing the functions of the host and terminal according to the first to fifth embodiments to a system or apparatus, and a computer of the system or apparatus. Needless to say, this can also be achieved by (or CPU or MPU) reading and executing the program code stored in the storage medium.
In this case, the program code itself read from the storage medium realizes the functions of the first to fifth embodiments, and the program code and the storage medium storing the program code constitute the present invention. It will be.
As a storage medium for supplying the program code, ROM, flexible disk, hard disk, optical disk, magneto-optical disk, CD-ROM, CD-R, magnetic tape, nonvolatile memory card, and the like can be used.
Further, by executing the program code read by the computer, not only the functions of the first to fifth embodiments are realized, but also an OS running on the computer based on the instruction of the program code. Needless to say, the present invention includes a case where the functions of the first to fifth embodiments are realized by performing part or all of the actual processing.
Further, after the program code read from the storage medium is written to the memory provided in the extension function board inserted in the computer or the function extension unit connected to the computer, the function extension is performed based on the instruction of the program code. It goes without saying that the CPU or the like provided in the board or the function expansion unit performs part or all of the actual processing, and the functions of the first to fifth embodiments are realized by the processing.
[0348]
FIG. 21 shows an example of the configuration of the computer function 800 described above.
As shown in FIG. 21, the computer function 800 includes a CPU 801, a ROM 802, a RAM 803, a keyboard controller (KBC) 805 of a keyboard (KB) 809, and a CRT controller (CRT) 810 as a display unit (CRT controller 810). CRTC) 806, a disk controller (DKC) 807 of a hard disk (HD) 811 and a flexible disk (FD) 812, and a network interface controller (NIC) 808 for connection to the network 840 via a system bus 804 It is the structure connected so that communication was possible mutually.
[0349]
The CPU 801 generally controls each component connected to the system bus 804 by executing software stored in the ROM 802 or the HD 811 or software supplied from the FD 812.
That is, the CPU 801 reads out and executes a processing program according to a predetermined processing sequence from the ROM 802, the HD 811 or the FD 812, thereby performing control for realizing the operations in the first to fifth embodiments. .
[0350]
The RAM 803 functions as a main memory or work area for the CPU 801.
The KBC 805 controls an instruction input from the KB 809 or a pointing device (not shown).
The CRTC 806 controls the display of the CRT 810.
The DKC 807 controls access to the HD 811 and the FD 812 that store a boot program, various applications, an edit file, a user file, a network management program, a predetermined processing program in the first to sixth embodiments, and the like.
The NIC 808 controls bidirectional data exchange with other devices or systems on the network 840.
[0351]
【The invention's effect】
As described above, in the present invention, a radiographic image processing apparatus, an image processing system, and a radiation that can obtain a good radiographic image from which an image component due to the grid is removed from a radiographic image obtained by radiography using the grid. An image processing method, a storage medium, and a program can be provided.
[Brief description of the drawings]
FIG. 1 is a diagram for describing spatial frequency characteristics of a filter for extracting a grid stripe component from a target image in the first to fifth embodiments.
FIG. 2 is a diagram for explaining an example of a grid stripe component extracted from a target image by the filter.
FIG. 3 is a diagram for explaining another example of grid stripe components extracted from a target image by the filter.
FIG. 4 is a diagram for explaining a spatial frequency characteristic in defective pixel correction of a target image.
FIG. 5 is a diagram for describing a spatial frequency characteristic in defective pixel correction for a target image in which the grid stripe component exists.
FIG. 6 is a block diagram showing a configuration of an X-ray image acquisition apparatus to which the present invention is applied in the first embodiment.
FIG. 7 is a diagram for explaining the operation timing of the X-ray image acquisition apparatus.
FIG. 8 is a diagram for explaining an example (example 1) of a defective pixel in a target image in the X-ray image acquisition apparatus.
FIG. 9 is a diagram for explaining an example (Example 2) of a state of the defective pixel.
FIG. 10 is a diagram for explaining an example (Example 3) of the defective pixel.
FIG. 11 is a diagram for explaining an example (Example 4) of the state of the defective pixel;
FIG. 12 is a diagram for explaining a spatial frequency distribution of grid stripe components in a target image in the X-ray image acquisition apparatus.
FIG. 13 is a diagram for explaining detection (analysis) of a grid stripe component for a target image in the X-ray image acquisition apparatus.
FIG. 14 is a flowchart for explaining detection (analysis) processing of the grid stripe component.
FIG. 15 is a flowchart for explaining a process of extracting a grid stripe component from a target image based on the result of the grid stripe component detection (analysis) process.
FIG. 16 is a block diagram showing a configuration of an X-ray image acquisition apparatus to which the present invention is applied in the second embodiment.
FIG. 17 is a block diagram showing a configuration of an X-ray image acquisition apparatus to which the present invention is applied in the third embodiment.
FIG. 18 is a block diagram showing a configuration of an X-ray image acquisition apparatus to which the present invention is applied in the fourth embodiment.
FIG. 19 is a flowchart for explaining a process of extracting a grid stripe component from a target image based on a result of the grid stripe component detection (analysis) process in the fifth embodiment.
FIG. 20 is a diagram for explaining the process of extracting the grid stripe component with a specific example;
FIG. 21 is a block diagram illustrating an example of a configuration for reading and executing a program from a computer-readable storage medium storing a program for causing a computer to implement the functions of the first to fifth embodiments.
FIG. 22 is a diagram for explaining radiation imaging using a grid.
FIG. 23 is a diagram for explaining an example of the effect of filtering on an image obtained by the radiography.
FIG. 24 is a diagram for explaining another example of the effect of filtering on the image obtained by the radiography.
FIG. 25 is a diagram for describing an example of the effect of filtering on an image obtained by superimposing grid stripe components obtained by radiography.
FIG. 26 is a diagram for explaining an example of the effect of filtering on an image obtained by superimposing grid stripe components obtained by the radiography.
[Explanation of symbols]
100 X-ray image acquisition apparatus
101 X-ray generator
102 subjects
103 grid
104 X-ray sensor
105 controller
106 Analog / digital (A / D) converter
107 memory
108 memory
109 Calculator
110 conversion table
111 memory
112 Calculator
113 memory
114 memory
115 Correction processing unit
116 Grid stripe detector
117 Grid stripe component extraction unit
118 memory
119 arithmetic unit
120 memory
121 Image processing unit
140 Memory transfer control unit
141 memory
142 Memory transfer control unit
143 memory

Claims (17)

A radiographic image processing apparatus that processes a radiographic image obtained by radiography using a grid for removing scattered radiation from a subject,
Creating means for creating an image component resulting from the grid from the predetermined radiation image;
Removing means for performing a removal process of the image component on a plurality of the radiation images obtained after the predetermined radiation image based on the image component obtained by the creating means;
The radiographic image processing apparatus, wherein the creating means includes means for creating the image component based on a feature that an image component superimposed on the radiographic image is steady over the entire image. A radiographic image processing method for processing a radiographic image obtained by radiography using a grid for removing scattered radiation from a subject,
A creation step for creating an image component resulting from the grid from a radiation image;
A removal step of performing a removal process of the image component on a plurality of radiation images obtained after the radiation image based on the image component obtained by the creation step;
The radiographic image processing method, wherein the creating step includes a step of creating the image component based on a feature that an image component superimposed on the radiographic image is steady over the entire image.   The radiographic image processing method according to claim 2, wherein the removing step includes a step of performing a removing process on the plurality of radiographic images using the same image component obtained by the creating step.   The radiographic image processing method according to claim 2, wherein the creating step includes an analysis step of analyzing the radiographic image to obtain a spatial frequency of the image component. The creation step above is
An extraction step for extracting a predetermined component including an image component from the radiation image based on the analysis result of the analysis step;
A processing step of processing the predetermined component obtained by the extraction step to obtain the image component;
Have
The removal step, a radiation image processing method according to claim 4, wherein the image components obtained by the processing steps of the creating step, wherein the benzalkonium be removed from the radiation image.   6. The radiographic image processing method according to claim 5, wherein the extraction step includes a step of performing filtering to extract a component having the spatial frequency obtained in the analysis step from the radiographic image.   6. The radiographic image processing method according to claim 5, wherein the processing step includes a step of executing processing for estimating and processing an unsteady portion of the predetermined component from a steady portion before and after the predetermined portion.   The radiographic image processing method according to claim 7, wherein the processing step includes a step of detecting the unsteady portion based on envelope information of the predetermined component.   The processing step estimates the amplitude and phase of a sine wave corresponding to the image component based on the stationary portion before and after the unsteady portion of the predetermined component and the spatial frequency, and repairs the unsteady portion. The radiation image processing method according to claim 7, further comprising the step of:   The radiographic image processing method according to claim 9, wherein the processing step includes a step of further filtering the predetermined component including the repaired unsteady portion to obtain the image component.   8. The radiographic image processing method according to claim 7, wherein the processing step substitutes a portion satisfying a predetermined condition in the non-stationary portion with a predetermined value to obtain the image component.   The radiographic image processing method according to claim 2, wherein the creating step includes a step of creating the image component for a predetermined line selected from the radiographic image.   The radiographic image processing method according to claim 12, wherein the creating step includes a step of creating the image component for an average result of a plurality of lines.   13. The creation step includes a step of using, as an image component of a line for which no image component has been created, the image component acquired from a line in which the image component in the vicinity of the line is created. The radiation image processing method as described. A detection step for detecting the mounting of the grid,
The radiographic image processing method according to claim 2, wherein the creating step includes a step of executing the processing based on a detection result of the detection step.   A computer-readable storage medium storing a program for causing a computer to execute the processing steps of the radiographic image processing method according to any one of claims 2 to 15.   The program for making a computer perform the process step of the radiographic image processing method of any one of Claims 2-15.

Images