(57)【要約】【課題】デジタル断層撮影において、画質の良い結像画像を高速で生成し表示することができるＸ線診断装置、当該Ｘ線診断装置における結像処理に用いられるフィルタ係数の算出方法を提供すること。【解決手段】複数のＸ線管位置で複数のＸ線画像を撮影し得られる各画像に前処理を適用した後、フィルタリング処理、バックプロジェクションを行うことで結像画像を生成する。このとき、フィルタリングに用いられるフィルタ係数は、スキャン軌道毎、画像フレーム毎、画像内の画素位置毎に応じて好適なものが決定される。【選択図】図３
2008119457.tif 000002
【特許請求の範囲】
【請求項１】
所定のスキャン軌道に沿ってＸ線管を移動させながら被検体に対しＸ線を曝射し取得されたフレーム毎のＸ線画像と、スキャン軌道毎、Ｘ線画像のフレーム毎、検出器の画素位置毎に異なるフィルタ係数とを記憶する記憶ユニットと、
前記所定のスキャン軌道、前記Ｘ線画像のフレーム、前記検出器の画素位置の組み合わせに基づいて各Ｘ線画像の画素毎のフィルタ係数を決定し、決定されたフィルタ係数を用いて前記フレーム毎のＸ線画像又はその処理画像にフィルタ処理を実行するフィルタ処理ユニットと、
フィルタ処理後の前記フレーム毎のＸ線画像を用いて、断層像を生成する断層像生成ユニットと、
を具備することを特徴とする画像処理装置。
【請求項２】
前記所定のスキャン軌道を入力するための入力ユニットをさらに具備することを特徴とする請求項１記載の画像処理装置。
【請求項３】
前記記憶ユニットが記憶するフィルタ係数は、スキャン軌道、Ｘ線画像のフレーム、検出器の画素位置の組み合わせ毎に等価ボケ関数を二次元画像として計算し、
前記等価ボケ関数の逆畳み込み関数を計算し、
前記逆畳み込み関数を用いて、スキャン軌道、Ｘ線画像のフレーム、検出器の画素位置の組み合わせ毎に決定されたものであること、
を特徴とする請求項１又は２記載の画像処理装置。
【請求項４】
前記フィルタユニットは、各画素のフィルタ係数に補間処理を用いるものであることを特徴とする請求項１乃至３のうちいずれか一項記載の画像処理装置。
【請求項５】
前記フィルタ係数のそれぞれは、中央画素において高い値を持ち、中央画素から遠い画素では０に近い値となり、中央画素を基準として決定される対向する二方向に負の値を持つ領域が放射状に延びるものであることを特徴とする請求項１乃至４のうちいずれか一項記載の画像処理装置。
【請求項６】
前記対向する二方向は、投影直線と焦点の移動方向を含む平面と、検出面が交わる線の方向であることを特徴とする請求項５記載の画像処理装置。
【請求項７】
前記負の値を持つ領域においては、中央画素に近い点では絶対値が大きな値を持ち、遠い点では０に近い値となることを特徴とする請求項５記載の画像処理装置。
【請求項８】
所定のスキャン軌道に沿ってＸ線管を移動させながら被検体に対しＸ線を曝射し、検出器の検出面に入射したＸ線を検出する撮影ユニットと、
検出された前記Ｘ線に基づいて、フレーム毎のＸ線画像を生成する画像生成ユニットと、
スキャン軌道毎、Ｘ線画像のフレーム毎、検出器の画素位置毎に異なるフィルタ係数を記憶する記憶ユニットと、
前記所定のスキャン軌道、前記Ｘ線画像のフレーム、前記検出器の画素位置の組み合わせに基づいて各Ｘ線画像の画素毎のフィルタ係数を決定し、決定されたフィルタ係数を用いて前記フレーム毎のＸ線画像又はその処理画像にフィルタ処理を実行するフィルタ処理ユニットと、
フィルタ処理後の前記フレーム毎のＸ線画像を用いて、断層像を生成する断層像生成ユニットと、
を具備することを特徴とするＸ線診断装置。
【請求項９】
前記所定のスキャン軌道を入力するための入力ユニットをさらに具備することを特徴とする請求項８記載のＸ線診断装置。
【請求項１０】
前記記憶ユニットが記憶するフィルタ係数は、スキャン軌道、Ｘ線画像のフレーム、検出器の画素位置の組み合わせ毎に等価ボケ関数を二次元画像として計算し、
前記等価ボケ関数の逆畳み込み関数を計算し、
前記逆畳み込み関数を用いて、スキャン軌道、Ｘ線画像のフレーム、検出器の画素位置の組み合わせ毎に決定されたものであることを特徴とする請求項８又は９記載のＸ線診断装置。
【請求項１１】
前記フィルタユニットは、各画素のフィルタ係数に補間処理を用いるものであることを特徴とする請求項８乃至１０のうちいずれか一項記載のＸ線診断装置。
【請求項１２】
前記フィルタ係数のそれぞれは、中央画素において高い値を持ち、中央画素から遠い画素では０に近い値となり、中央画素を基準として決定される対向する二方向に負の値を持つ領域が放射状に延びるものであることを特徴とする請求項８乃至１１のうちいずれか一項記載の画像処理装置。
【請求項１３】
前記対向する二方向は、投影直線と焦点の移動方向を含む平面と、検出面が交わる線の方向であることを特徴とする請求項１２記載の画像処理装置。
【請求項１４】
前記負の値を持つ領域においては、中央画素に近い点では絶対値が大きな値を持ち、遠い点では０に近い値となることを特徴とする請求項１２記載の画像処理装置。
【請求項１５】
Ｘ線診断装置におけるディジタル断層撮影の再構成処理に用いられるフィルタ係数の算出プログラムであって、
コンピュータに、
スキャン軌道毎、Ｘ線画像のフレーム毎、検出器の画素位置の組み合わせ毎に等価ボケ関数を二次元画像として計算させる機能と、
前記等価ボケ関数の逆畳み込み関数を計算させる機能と、
前記逆畳み込み関数を用いて、スキャン軌道、Ｘ線画像のフレーム、検出器の画素位置の組み合わせ毎のフィルタ係数を計算させる機能と、
を実現させることを特徴とするフィルタ係数の算出プログラム。
【請求項１６】
前記コンピュータに、計算された前記フィルタ係数を、スキャン起動の情報と関連付けて保存手段に保存させる機能をさらに実現させることを特徴とする請求項１５記載のフィルタ係数の算出プログラム。
【請求項１７】
Ｘ線診断装置におけるディジタル断層撮影の再構成処理に用いられるフィルタ係数の算出プログラムであって、
コンピュータに、
逆投影、投影、逆投影、投影の二重積分を実行することにより、等価ボケ関数を計算させる機能と、
前記等価ボケ関数の逆畳み込み関数を計算させる機能と、
前記逆畳み込み関数を用いて、スキャン軌道、Ｘ線画像のフレーム、各Ｘ線画像上の画素位置の組み合わせ毎のフィルタ係数を計算させる機能と、
を実現させることを特徴とするフィルタ係数の算出プログラム。
【請求項１８】
前記コンピュータに、計算された前記フィルタ係数を、スキャン起動の情報と関連付けて保存手段に保存させる機能をさらに実現させることを特徴とする請求項１７記載のフィルタ係数の算出プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、Ｘ線断層画像生成において用いられる画像再構成処理に用いられるフィルタ係数の算出プログラムに関する。
【背景技術】
【０００２】
近年、Ｘ線診断装置等を用いたデジタル断層撮影の結像処理方法として、様々な方法が提案されている（例えば、非特許文献１参照）。その代表的なものに、図１５に示すようなfiltered backprojection(FBP)法と呼ばれるものがある。このFBP法は、複数の方向から撮影したＸ線画像に対数演算などの処理を施した複数の画像にフィルタを適用し、それらをバックプロジェクションすることで多断層の結像画像を得るものである。フィルタとしては、例えばShepp-Logan filter（図１５の式参照）と呼ばれるフィルタなどが知られている。FBP法は、多方向から撮影した全てのＸ線画像に同じフィルタ係数を適用し、さらに1枚の画像の中の位置によっても係数を変えることなく同じフィルタを適用することが特徴である。filtered backprojection法は本来、完全再構成が可能なCTスキャナーのための再構成法であり、完全再構成が不可能なデジタル断層撮影においては適切な結像処理方法ではない。しかしながら、処理時間が短いという特徴があるためfiltered backprojection法をデジタル断層撮影に適用する多くの研究報告がある。
【０００３】
また、デジタル断層撮影への適用において高い画質が得られる結像処理方法がいくつか提案されている。例えば、Maximum Likelihood algorithm(ML法)と呼ばれる反復的結像処理方法とfiltered backprojection法を比較した結果が報告されている（例えば、非特許文献２参照）。この報告では、ML法はfiltered backprojectionよりも画質が良いが、filtered backprojection法の処理時間が15分以下であるのに対し、ML法の処理時間は3時間以上であったことが示されている。さらに、例えばmatrix inversion tomosynthesisと呼ばれる結像処理方法が知られており（例えば、特許文献１参照）、この手法によれば、filtered backprojection法より高い画質が得られることが報告されている（例えば、非特許文献３参照）。
【特許文献１】米国特許第４９０３２０４号公報
【非特許文献１】James T Dobbins III, and Devon J Godfrey: Digital x-ray tomosynthesis: current state of the art and clinical potential, Physics in Medicine and Biology Vol. 48 pp.R65-R106, 2003
【非特許文献２】Tao Wu, Richard H. Moore, Elizabeth A. Rafferty, and Daniel B. Kopans, A comparison of reconstruction algorithms for breast tomosynthesis, Medical Physics Vol. 31, No. 9, September 2004
【非特許文献３】Ying Chen, Joseph Y. Lo, James T. Dobbins III, Noise power spectrum analysis for several digital breast tomosynthesis reconstruction algorithms, Proc. SPIE vol. 6142, pp. 1677-1684, 2006
【発明の開示】
【発明が解決しようとする課題】
【０００４】
臨床の現場において検査のためにＸ線断層画像を生成する場合、異常の有無や状態を読影できるだけの画質が提供されることが必須である。また、処理時間が短かく例えば数分以内で検査画像が提供できることは、実務上は画質が良いことと同程度に重要である。
【０００５】
しかしながら、上述した様に、従来のfiltered backprojection法は処理時間が短いが画質が劣るという問題がある。また、他の画質の良い結像処理法は多くの処理時間がかかり、実用に耐えないという問題がある。
【０００６】
本発明は、上記事情を鑑みてなされたもので、デジタル断層撮影において、画質の良い結像画像を高速で生成し表示することができるＸ線診断装置または画像処理装置、Ｘ線診断装置等における結像処理に用いられるフィルタ係数の算出プログラムを提供することを目的としている。
【課題を解決するための手段】
【０００７】
本発明は、上記目的を達成するため、次のような手段を講じている。
【０００８】
請求項１に記載の発明は、所定のスキャン軌道に沿ってＸ線管を移動させながら被検体に対しＸ線を曝射し取得されたフレーム毎のＸ線画像と、スキャン軌道毎、Ｘ線画像のフレーム毎、検出器の画素位置毎に異なるフィルタ係数とを記憶する記憶ユニットと、前記所定のスキャン軌道、前記Ｘ線画像のフレーム、前記検出器の画素位置の組み合わせに基づいて各Ｘ線画像の画素毎のフィルタ係数を決定し、決定されたフィルタ係数を用いて前記フレーム毎のＸ線画像又はその処理画像にフィルタ処理を実行するフィルタ処理ユニットと、フィルタ処理後の前記フレーム毎のＸ線画像を用いて、断層像を生成する断層像生成ユニットと、を具備する画像処理装置である。
【０００９】
請求項８に記載の発明は、所定のスキャン軌道に沿ってＸ線管を移動させながら被検体に対しＸ線を曝射し、検出器の検出面に入射したＸ線を検出する撮影ユニットと、検出された前記Ｘ線に基づいて、フレーム毎のＸ線画像を生成する画像生成ユニットと、スキャン軌道毎、Ｘ線画像のフレーム毎、検出器の画素位置毎に異なるフィルタ係数を記憶する記憶ユニットと、前記所定のスキャン軌道、前記Ｘ線画像のフレーム、前記検出器の画素位置の組み合わせに基づいて各Ｘ線画像の画素毎のフィルタ係数を決定し、決定されたフィルタ係数を用いて前記フレーム毎のＸ線画像又はその処理画像にフィルタ処理を実行するフィルタ処理ユニットと、フィルタ処理後の前記フレーム毎のＸ線画像を用いて、断層像を生成する断層像生成ユニットと、を具備するＸ線診断装置である。
【００１０】
請求項１５に記載の発明は、Ｘ線診断装置におけるディジタル断層撮影の再構成処理に用いられるフィルタ係数の算出プログラムであって、コンピュータに、スキャン軌道毎、Ｘ線画像のフレーム毎、検出器の画素位置の組み合わせ毎に等価ボケ関数を二次元画像として計算させる機能と、前記等価ボケ関数の逆畳み込み関数を計算させる機能と、前記逆畳み込み関数を用いて、スキャン軌道、Ｘ線画像のフレーム、検出器の画素位置の組み合わせ毎のフィルタ係数を計算させる機能と、を実現させることを特徴とするフィルタ係数の算出プログラムである。
【００１１】
請求項１７に記載の発明は、Ｘ線診断装置におけるディジタル断層撮影の再構成処理に用いられるフィルタ係数の算出プログラムであって、コンピュータに、逆投影、投影、逆投影、投影の二重積分を実行することにより、等価ボケ関数を計算させる機能と、前記等価ボケ関数の逆畳み込み関数を計算させる機能と、前記逆畳み込み関数を用いて、スキャン軌道、Ｘ線画像のフレーム、各Ｘ線画像上の画素位置の組み合わせ毎のフィルタ係数を計算させる機能と、を実現させることを特徴とするフィルタ係数の算出プログラムである。
【発明の効果】
【００１２】
以上本発明によれば、デジタル断層撮影において、画質の良い結像画像を高速で生成し表示することができるＸ線診断装置、画像処理装置及びＸ線診断装置等における結像処理に用いられるフィルタ係数の算出プログラムを実現することができる。
【発明を実施するための最良の形態】
【００１３】
以下、本発明の実施形態を図面に従って説明する。なお、以下の説明において、略同一の機能及び構成を有する構成要素については、同一符号を付し、重複説明は必要な場合にのみ行う。
【００１４】
なお、本実施形態においては、説明を具体的にするため、本発明の技術的思想を乳房撮影用Ｘ線診断装置に適用した場合を例とする。しかしながら、乳房撮影用に限定されず、Ｘ線を用いたデジタル断層撮影など再構成処理を行うもの（例えば循環器撮影用等のＸ線診断装置等）であれば、適用することが可能である。
【００１５】
図１は、本実施形態に係る乳房撮影用Ｘ線診断装置１の外観図を示している。同図に示すように、本乳房撮影用Ｘ線診断装置１は、アーム部１０、支柱部３５から構成されている。アーム部１０には、互いに対向して設置されるＸ線源装置１３及び平面検出器２０、圧迫板１５等が設けられている。また、支柱部３５には、表示パネル２７ａ、フットペダル３１ａ、３１ｂ、タッチパネル３１ｃ等が設けられている。アーム部１０は、支柱部３５に対し上下移動及び軸３６を中心とした回転移動が可能である。このアーム部１０を所定の姿勢にし、平面検出器２０の検出面、圧迫板１５によって乳房を固定することで、例えばトモシンセシス（断層撮影）を行う場合に所望のスキャン軌道によって画像を取得することができる。
【００１６】
図２は、本乳房撮影用Ｘ線診断装置１のブロック構成図を示している。同図に示すように、本乳房撮影用Ｘ線診断装置１は、Ｘ線制御部１１、Ｘ線源装置１３、圧迫板１５、圧迫板制御部１６、圧迫板駆動部１７、平面検出器２０、メモリ２１、前処理部２３、結像処理部２４、画像処理部２５、Ｄ／Ａ変換器２６、表示部２７、中央制御部３０、操作部３１、記憶部３３を具備している。
【００１７】
Ｘ線制御部１１は、中央制御部３０からの指示に従って、所定強度のＸ線が所定レートで曝射されるように、Ｘ線源装置１３を制御する。
【００１８】
Ｘ線源装置１３は、Ｘ線管１３０、Ｘ線照射野限定器１３１を有している。Ｘ線管１３０は、Ｘ線を発生する真空管であり、高電圧発生装置からの高電圧により電子を加速させ、ターゲットに衝突させることでＸ線を発生する。Ｘ線照射野限定器１３１は、Ｘ線管１３０から照射されるＸ線を所定の形状に成形する。
【００１９】
圧迫板１５は、平面検出器２０の検出面とによって乳房を挟み圧迫することで、撮影時において乳房を平たく薄く固定するためのものである。このように乳房を圧迫して固定（ポジショニング）することで、被写体からの散乱Ｘ線を減少させ、乳腺組織の重なりを少なくすることができ、画像コントラストの改善、体動によるノイズ発生等を防止、照射線量の低減等を実現することができる。
【００２０】
圧迫板制御部１６は、中央制御部３０からの指示のもと、圧迫板１５の位置決めを行うために、圧迫板駆動部１７を制御する。また、圧迫板制御部１６は、固定された乳房の圧迫厚（圧迫時の厚さ）を計測し、中央制御部３０に送り出す。
【００２１】
圧迫板駆動部１７は、圧迫板制御部１６の制御のもと、圧迫板１５を駆動する。
【００２２】
平面検出器２０は、シンチレータとフォトダイオードアレイとを有し、被写体を透過したＸ線を光電膜に当てることで電子正孔を生成し、これを半導体スイッチにおいて蓄積し、電気信号として読み出すことでＸ線を電気信号に変換して検出するものである。なお、変換方式は、Ｘ線から電気信号に変換する直接変換であっても良いし、Ｘ線から光を介して電気信号に変換する間接変換であっても良い。
【００２３】
メモリ２１は、平面検出器２０から供給されるデジタル信号を一時的に記憶する。
【００２４】
前処理部２３は、メモリ２１からデジタル信号をフレーム毎に読み出してデジタルＸ線画像を生成し、前処理（例えば、オフセット補正処理、血管画素補正処理、対数演算処理等）を実行する。なお、この前処理においては、必要に応じてピクセル束ねが実施される場合もある。
【００２５】
結像処理部２４は、中央制御部３０による制御のもと、後述するProjection Inversion法による結像処理を実行する。
【００２６】
画像処理部２５は、フレーム毎のＸ線画像を用いて、必要に応じてサブトラクション処理等の所定の画像処理を行う。
【００２７】
Ｄ／Ａ変換器２６は、画像処理部２３から入力した画像データのデジタル信号列を、アナログ信号列に変換する。
【００２８】
表示部２７は、Ｄ／Ａ変換器２６から受け取った信号により、トモシンセシスによって得られるデジタル断層画像等のＸ線画像を表示するＣＲＴ、プラズマディスプレイ、液晶ディスプレイ等の他、装置の動作状態を示す表示パネル２７ａを有する。また、表示部２７は、後述するProjection Inversion法による結像処理において、スキャン軌道に関する条件を入力するための画面を表示する。
【００２９】
中央制御部３０は、画像データの収集に関する制御、及び収集した画像データの画像処理、画像再生処理等に関する制御を行う中央処理装置である。特に、中央制御部３０は、後述するProjection Inversion法による結像処理に関する制御を実行する。
【００３０】
操作部３１は、キーボードや各種スイッチ、マウス、フットペダル３１ａ、３１ｂ、タッチパネル３１ｃ等を備えた入力装置であり、患者情報（患者ＩＤ、検査部位、検査目的等）の入力、圧迫板１５の上下移動、撮影指示、パルスレート選択、画像選択、スキャン軌道に関する条件等を入力する際に使用される。
【００３１】
記憶部３３は、画像処理部２３での画像処理前又は処理後の画像データ等を記憶する。また、記憶装置３３は、後述するProjection Inversion法による結像処理を実現するための専用プログラムを記憶する。さらに、記憶装置３３は、Projection Inversion法による結像処理において用いられるフィルタ係数（すなわち、スキャン軌道毎、画像フレーム毎、画像上の位置毎に決定されるフィルタ係数）を記憶する。
【００３２】
（Projection Inversion法による結像処理）
次に、Projection Inversion法（PI法）による結像処理について説明する。この処理は、大きくフィルタ係数計算処理と断層像生成処理とに分類することができる。以下、それぞれの処理について説明する。
【００３３】
［１．断層像生成処理］
図３は、PI法による結像処理の流れを概念的に示した図である。同図に示すように、
断層像生成処理は、スキャン軌道に沿って複数のＸ線画像を撮影した後、Ｘ線画像に前処理、フィルタリングを施しバックプロジェクションを実施して断層像を得る処理である。フィルタリングには、後述するフィルタ係数計算処理により算出されたフィルタ係数が用いられる。このフィルタ係数は、スキャン軌道毎、画像フレーム毎、画像上の位置毎に異なる。従って、撮影に用いたスキャン軌道によって、適切なフィルタ係数を選択するよう構成される。また、一回の撮影で得られた複数のＸ線画像の各々に対して異なるフィルタ係数が適用されるように、フィルタ係数を選択する構成となっている。さらには、フィルタ係数は、Ｘ線画像の検出器画素位置に応じて異なるため、検出器画素位置に応じて異なるフィルタ係数が適用されるよう構成される。
【００３４】
PI法による結像処理の背景理論は、次のようである。すなわち、結像処理の結果得られるボクセルデータをベクトル形式でｙ、プロジェクション画像をベクトル形式でｇ、プロジェクション演算を行列形式でＷｙと表したとき、プロジェクション過程は、次の式（１）の様に書くことができる。
【数１】
2008119457.tif 000003
【００３５】
εはＸ線画像に加わるノイズである。ここで、プロジェクション画像ｇの各成分は被写体の吸収率分布により減衰したＸ線量の割合（減衰比）の対数をとったものである。
【００３６】
２乗平均推定法(minimum-mean-square-error estimation)によれば、方程式（１）の解は、次の式（２）で得られることが知られている。
【数２】
2008119457.tif 000004
【００３７】
ｇｍｅａｓは、撮影により得られたプロジェクション画像、σ２ｙ、σ２εはそれぞれボリュームの画素値の標準偏差とプロジェクション画像に含まれるノイズの標準偏差である。ここで、次の式（３）によって定義されるボケ補正プロジェクション画像ｘを導入する。
【数３】
2008119457.tif 000005
【００３８】
すると、式（２）は、次の式（４）で表すことができ、ｘがわかればバックプロジェクション演算により解を得ることができる。
【数４】
2008119457.tif 000006
【００３９】
なお、ＷＴｘはｘのバックプロジェクション演算を行うことを表す。式（３）より、ｘは次の式（５）の解である。
【数５】
2008119457.tif 000007
【００４０】
この式（５）によって表現される連立一次方程式を解けばｘは求められるが、この方程式が大次元であるため、直接法により解くことはできない。しかし、次の式（６）の目的関数を反復法を用いて最小化することで、解を得ることができる。
【数６】
2008119457.tif 000008
【００４１】
式（６）を最小化するためには共役勾配法などの最適化手法が用いられる。共役勾配法を用いて式（５）の解を得る場合、Ａｘを計算する処理ルーチンを用意しておけばよく、Ａ、Ｗ、ＷＴなどの大次元の行列を保持しなくても、計算を実行できるという特徴がある。
【００４２】
共役勾配法により式（５）の解が得られれば（これがボケ補正プロジェクション画像ｘである）、式（４）によりバックプロジェクションを行って断層像ｙｅｓｔを得ることができる。
【００４３】
この方法の重要な点の一つは、プロジェクション画像領域でinversionが行われていることである。すなわち、方程式（５）は、右辺、左辺ともプロジェクション領域の画像である。matrix inversion tomosynthesisや他の多くの反復解放は、断層像領域でinversionを行っているのとは対照的である。
【００４４】
式（２）及びこれを変形した式（５）は、デジタル断層撮影法における最適解の方程式であり、これを解く方法が画質の面で最も良い。特に、単純バックプロジェクション法やfiltered backprojection法などと比べ、大きな画質向上がある。しかし、共役勾配法を用いる上記の方法も、数回の反復を実行するため、処理時間が多くかかってしまう。すなわち、上述したすべての方法は、ボケ補正Ｘ線画像ｘを得るために、方程式（５）を直接解いている。そのため、実際にこれを実行するには多くの処理時間がかかるという問題を持っている。
【００４５】
しかしながら、方程式（５）は線形方程式であり右辺、左辺ともプロジェクション領域の画像であることを考えると、ボケ補正Ｘ線画像ｘを得るためのフィルタ演算が存在する可能性がある。フィルタ演算であれば方程式を解くよりはるかに高速に処理を実行することが可能となる。本発明者らは、実際、後述する手法によりそのフィルタを導出し、方程式（５）を解くフィルタ演算が存在することを確かめた。また、フィルタ演算により等価な断層像が得られ、さらに、大幅に処理時間を削減できることを確かめた。
【００４６】
［２．フィルタ係数計算処理］
次に、フィルタ係数計算処理について説明する。この処理では、Ｘ線画像撮影時のスキャン軌道（何枚の画像を撮影するか、及び各々の画像をどの焦点位置、検出器位置、検出器向きで撮影するかを表す）に応じて、Ｘ線画像に適用する適切なフィルタ係数を算出する。この処理は、Ｘ線画像を撮影した後で実施できるが、撮影するスキャン軌道が前もって知られていれば、Ｘ線画像の撮影前に実施しておくこともできる。また、撮影に用いられるスキャン軌道が少数の種類であるならば、その少数のスキャン軌道に対するフィルタ係数をあらかじめ算出しておき、撮影後の断層像生成処理には、あらかじめ算出されたフィルタ係数を適用することができる。
【００４７】
フィルタ係数は、スキャン軌道毎に異なるが、さらに、画像フレーム毎（すなわち、焦点位置、及び検出器位置、検出器向き等の条件毎）や同一画像内の画素位置毎（画像内の領域毎）によっても異なる。従って、スキャン軌道毎、画像フレーム毎、画像内の画素位置（又は検出器画素位置）毎に応じて異なるフィルタ係数が計算される。
【００４８】
図４は、PI法による結像処理に用いられるフィルタ係数の計算処理の流れを概念的に示した図である。同図に示すように、ボケ補正Ｘ線画像を求めるには、一旦、プロジェクション領域の等価ボケ関数ｐｋを求め、その後、ＦＦＴを用いてｐｋを応答関数とするデコンボリューションフィルタを構成する（逆フィルタ算出）という二つステップを採用する。
【００４９】
［２−１．等価ボケ関数の算出処理］
図５乃至図８は、等価ボケ関数を算出する手順を示した図である。プロジェクション画像のあるフレームのある検出画素hが1で他の画素が全て0である状態を想定する（図５中（１））。以降の処理に等価な積分演算が実行できればよいので、画像セットを作る必要はないが、説明のため、このような画像セットｐ１を作るものとする。また、この画像を結像領域にバックプロジェクションし、バックプロジェクション画像ｂ１を得る（図５中（２））。
【００５０】
次に、バックプロジェクション画像ｂ１を全てのフレームにプロジェクションして、プロジェクション画像ｐ２を得る（図６中（３））。その後、プロジェクション画像ｐ２をバックプロジェクションし、ｂ２を得る（図７中（４））。バックプロジェクション画像ｂ２を画素ｈが存在するひとつのフレームにプロジェクションする（図８中（５））。このバックプロジェクションにより得られる１枚のプロジェクション画像から、画素ｈの周辺の画像領域の画像を抜き出したものが等価ボケ関数ｐｋである。
【００５１】
この様なバックプロジェクション−プロジェクション−バックプロジェクション−プロジェクションという演算は、ひとつの２重積分の式で表現することができる。従って、実装上は等価ボケ関数ｐｈはこの２重積分を実行することで算出される。
【００５２】
等価ボケ関数ｐｈの積分式は、上記の逆投影-投影-逆投影-投影の過程をステップごとに追っていくことで得られる。ここではその計算方法のみを示す。
【数７】
2008119457.tif 000009
【００５３】
ｄｈはフィルタ中心の座標値、ｄは検出器上のボケ関数の値を調べたい位置である。上の積分はｄについて陽な式にはなっておらず、ｄはｓ，ｔの変数となっている。しかし、数値上は上式で計算することができ、ｔ，ｓのループ内でｄを決定し、その画素へ積分を累積してゆけばよい。
【００５４】
Ｌは積分領域の大きさを表すパラメータ（積分領域サイズ）で、およそ被写体の厚さの長さを指定する。ｔ１，ｔ２はｔの積分範囲でＬを用いて、次の（ａ３）、（ａ４）で求められる。
【数８】
2008119457.tif 000010
【００５５】
ｒｓｔ，ｒｓｌはそれぞれプロジェクションk, lでの焦点の座標である。ｒｋ（ｓ），ｒｌ（ｓ，ｔ）はプロジェクションk, lの投影直線上の座標であり、ｓ，ｔを用いて次式（ａ５）、（ａ６）で表される。
【数９】
2008119457.tif 000011
【００５６】
ｍｋ（ｓ），ｍｌ（ｓ，ｔ）は幾何学的拡大率で、次の（ａ７）で求められる。
【数１０】
2008119457.tif 000012
【００５７】
θｋ，θｌは投影直線と検出面のなす角で、次の（ａ８）で求められる。
【数１１】
2008119457.tif 000013
【００５８】
ｗｋ，ｗｌはプロジェクションｋ，ｌでの検出器の法線方向を表すベクトルである。ｄはプロジェクションkでの焦点とｒｌ（ｓ，ｔ）を結ぶ投影直線と検出面の交点である。幾何学的拡大率ｍｋ（ｓ，ｔ）を用いて、次の（ａ９）の様に容易に計算できる。
【数１２】
2008119457.tif 000014
【００５９】
［２−２．逆フィルタ算出処理］
まず、最初にベクトルａｋを求める。ａｋは等価ボケ関数ｐｋに画素ｈに対応する画素のみに値γ２＝σ２ε／σ２ｙを加えたベクトルである。ａｋ、ｐｋともに２次元画像を表すことに注意する。
【００６０】
ボケ補正プロジェクション画像ｘのうち、ｋ番目のフレームの画像をｘｋとし、撮影したプロジェクション画像ｇｍｅａｓのｋ番目のフレームの画像をｇｋｍｅａｓとすると、方程式（５）は次の式（８）のように近似的に表される。
【数１３】
2008119457.tif 000015
【００６１】
式（８）において、ａｋとｇｋｍｅａｓが既知でｘｋが未知なのだから、デコンボリューション演算によりｘｋを求めることができる。デコンボリューション演算は、次の式（９）のようにＦＦＴアルゴリズムを用いて実行される。
【数１４】
2008119457.tif 000016
【００６２】
なお、Ｆ［．］は２次元フーリエ変換をあらわす。ここで求めたい逆フィルタは、次の式（１０）によって求められる。
【数１５】
2008119457.tif 000017
【００６３】
［２−３．PI法での断層像生成処理］
フィルタ係数ｆを用いれば断層像の生成処理は、式（１１）によって表されるフィルタリング、式（４）によって表されるバックプロジェクションの二つのステップで実行される。
【数１６】
2008119457.tif 000018
【００６４】
［２−４．算出したフィルタ係数の特徴］
図９は、あるフレームにおける、ある画素位置でのフィルタ係数を算出し、画像化して示した例である。この様なprojection inversion法で算出したフィルタ係数は、画素位置毎の値を持ち、一般に次のような性質を持つ。
【００６５】
・中央画素（画素ｈ）において高い値を持つ。
【００６６】
・中央画素（画素ｈ）から遠い画素では０に近い値となる。
【００６７】
・中央画素（画素ｈ）から対向する２方向に負の値を持つ領域が放射状に延びる。
【００６８】
なお、対向する２方向とは、図１０に示す投影直線（画素ｈと焦点を結ぶ線）と焦点の移動方向を含む平面と、検出面が交わる線の方向である。図１０に示したように、スキャン軌道から等価ボケ関数を計算する方法のためこのようなフィルタ係数が得られる。負の値を持つ領域においては、中央画素に近い点では絶対値が大きな値を持ち、遠い点では０に近い値となる。また、上記の性質は、フィルタ係数の算出に、スキャン軌道を用いて計算した等価ボケ関数を用いることによるものである。逆に上記のような性質を持っていれば、フィルタ係数算出の正確な構成にかかわらず類似の画質を得ることができる。
【００６９】
［２−５．本乳房撮影用Ｘ線診断装置におけるフィルタ係数決定］
以上述べた手法により算出されたフィルタ係数は、断層像生成の度に計算してもよい。しかしながら、断層像生成時におけるTurn around time削減の観点から、本乳房用Ｘ線診断装置１は、事前に（例えば装置の工場出荷時に）コンピュータ等を用いて計算されたスキャン軌道毎、画像フレーム毎、画素位置毎のフィルタ係数を予め記憶しておき、適切なフィルタ係数を選択して断層像生成処理を実行する。
【００７０】
フィルタ係数は、フィルタ係数同定文字列を用いて選択される。フィルタ係数同定文字列は、スキャン軌道同定文字列、スキャン軌道算出パラメータ、画像フレーム番号、画像上の画素位置を特定するためのものである。
【００７１】
スキャン軌道同定文字列は、撮影枚数(フレーム数)、スキャン角度（例えば、円弧スキャンの場合はどの角度範囲を撮影するか）、フレームレート（１秒に何枚のＸ線画像を撮影するか）、焦点・検出器の移動速度、ＳＩＤ（焦点、検出器間距離）の組み合わせによって定義される。例えば、撮影枚数４１、 スキャン角度（円弧スキャンの角度範囲）９０度、 フレームレート０．５［frame/sec］、ＳＩＤ９００ｍｍの場合のスキャン軌道同定文字列は、「SP90_4frms_05dfps_sid900」となる。
【００７２】
なお、本実施形態では、例えば後述するインターフェースを用いて予めプリセットされた値からスキャン軌道に関する条件を選択することで、スキャン軌道同定文字列の入力、設定を行うものとする。しかしながら、これに拘泥されず、直接任意の値や条件を入力することで、スキャン軌道同定文字列を設定するようにしてもよい。
【００７３】
また、スキャン軌道算出パラメータは、式（７）のγなどのパラメータ及びｘ，ｙ方向に何点のフィルタ係数代表点を取るかという数を含むものである。従って、スキャン軌道同定文字列が「SP90_41frms_05dfps_sid900」であり、例えば図１１中の第２フレームの横方向（例えばｘ方向）３番目、縦方向（例えばｙ方向）２番目の検出器画素位置のフィルタ係数同定文字列は、「SP90_41frms_05dfps_sid900/gamma10_n5_3/2_3_2」と定義される。なお、同文字列中の「n5_3」は、横方向に５点,縦方向に３点の代表点がフィルタ計算に使用されることを表す。なお、本実施形態でのフィルタ係数同定文字列は、「スキャン軌道同定文字列+"/"+フィルタ係数算出パラメータ+"/"+検出器画素位置を示す文字列」で表現している。
【００７４】
結像処理部２４は、設定されたフィルタ係数同定文字列に対応するフィルタ係数を記憶部３３から読み出し、それぞれフィルタ係数同定文字列と関連付けてファイル管理して（例えば、フィルタ係数同定文字列をファイル名の一部に用いて）保存する。
【００７５】
（動作）
次に、本乳房撮影用Ｘ線診断装置１のPI法による結像処理における動作について説明する。
【００７６】
図１２は、乳房撮影用Ｘ線診断装置１のPI法による結像処理における動作の流れを示したフローチャートである。同図において、まず、結像処理実行前に、操作部３１からスキャン軌道に関する条件（すなわち、スキャン軌道同定文字列を決定するための条件）、スキャン軌道算出パラメータ、及び他の結像処理条件（断層像の中心位置、断層像の大きさ(field of view)、断層像の画素間隔等）が入力される（ステップＳ１）。
【００７７】
図１３は、スキャン軌道に関する条件を入力するためのインターフェースの一例を示した図である。同図に示すような入力画面から所望の値を選択することにより、スキャン軌道に関する条件が入力され、これに応じたスキャン軌道同定文字列が設定される。なお、図１３の例（すなわち撮影枚数７２、 スキャン角度３０度、 フレームレート７．２［frame/sec］、SID700mmの場合）では、スキャン軌道同定文字列は、「SP90_72frms_7.2dfps_sid700」がスキャン軌道同定文字列となる。
【００７８】
次に、操作部３１から撮影開始が指示されると、中央制御部３０は、指定した枚数のＸ線画像の撮影を、指定したスキャン軌道に沿って実行し、複数のＸ線画像を取得する（ステップＳ２）。
【００７９】
次に、前処理部２３は、撮影した複数のＸ線画像に対して前処理を実行する（ステップＳ３）。なお、この前処理には、オフセット補正、欠陥画素補正処理、対数演算などが含まれ、また、必要に応じ、ピクセル束ねが実施される場合もある。
【００８０】
次に、中央制御部３０は、各画像の画素毎のフィルタ係数同定文字列を生成する（ステップＳ４）。結像処理部２４は、生成されたフィルタ係数同定文字列を用いてフィルタ係数を選択し、フィルタリング処理を実行する（ステップＳ５）。すなわち、結像処理部２４は、フィルタ係数同定文字列からフィルタ係数を計算した代表点の総数を取得し、全てのフィルタ係数をファイルから読み込む。その後、結像処理部２４は、下記の（１）〜（３）までの処理を全てのフレームｋ毎に、全ての検出器画素(i, j)に対し実施する
（１）検出器画素(i, j)がどの代表点(I, J)に近いかを決定する（すなわち、フィルタ係数同定文字列から代表点位置(I, J)を算出し判定する）。
【００８１】
（２）フィルタ係数同定文字列により関連するフィルタ係数を特定する（すなわち、フィルタ係数同定文字列に基づいて、記憶部３３から関連するフィルタ係数を選択する）。なお、今の場合、フィルタ係数同定文字列は「スキャン軌道同定文字列 + "/"+フィルタ係数算出パラメータ+"/"+"k_I_J"」で定義される。
【００８２】
（３）フレームk、検出器画素(i, j)に対して関連付けられたフィルタ係数によるフィルタリングを実施する。
【００８３】
次に、フィルタリング処理が施された画像を用いて、バックプロジェクション処理を実行することで、断層像を生成する（ステップＳ６）。このバックプロジェクション処理は、スキャン軌道同定文字列により特定されたスキャン軌道、及び指定した他の結像処理条件に基づいて実行する。生成された断層像は、表示部２７に表示され、自動的に記憶部３３に記憶される（ステップＳ７）。
【００８４】
以上述べた手法は、検出器の各画素に最も近い代表点でのフィルタ係数を用いる例であるが、補間によりフィルタ係数を決定する方法もある。例えば、検出器の画素(i,j)に近い４つの代表点を選択する。その際、４つの代表点が作る四角形内に(i,j)が含まれるように代表点を選ぶ。４つの代表点のフィルタ係数から、画素(i, j)でのフィルタ係数を双線形補間により求める。なお、双線形補間の場合、選択する代表点は４つであるが、他の補間方法を用いた場合、代表点を選択する数は異なるものになる。
【００８５】
また、フィルタ係数の算出処理において、上述した手法（直接構成法）以外に、フィルタ係数採取法によってもフィルタ係数を求めることができる。この手法によるフィルタ係数算出処理は次のように行う。まず、擬似的なプロジェクション画像を生成する。このプロジェクション画像は、結像領域に点およびstripが存在したと仮定したときのプロジェクション画像である（図１４（ａ））。点の投影像は格子状に配置され、stripの投影像はその隙間に挿入される。stripの向きは全11枚のプロジェクション画像のうち一つの画像に対してstripが点として投影される向きとする。このプロジェクション像を元に前述の式（６）を用いて反復的に解く方法または他の方法により式（５）の解xを求める。その結果、点投影像に対する応答が得られる（図１４（ｂ））。フィルタ係数採取法はこの応答をフィルタ係数として採取するものである。得られたフィルタ係数を用いてプロジェクション画像をフィルタリングし、バックプロジェクションを実施することで、式（６）による反復法と等価な処理をより高速に実行することができる。
【００８６】
以上述べた構成によれば、以下の効果を得ることができる。
【００８７】
本Ｘ線診断装置では、複数のＸ線管位置で複数のＸ線画像を撮影し得られる各画像に前処理を適用した後、フィルタリング処理、バックプロジェクションを行うことで断層像を生成する。このとき、フィルタリングに用いられるフィルタ係数は、スキャン軌道毎、画像フレーム毎、画像内の画素位置毎に応じて好適なものが決定され、また、このフィルタ係数の決定は、デジタル断層撮影法における最適解の方程式（式（２）、式（５）参照）をフィルタ演算によって求めることに対応する。従って、従来に比して画質の高い断層像を生成することができ、診断の質の向上に寄与することができる。
【００８８】
また、本Ｘ線診断装置では、断層像生成において反復演算処理を必要とせず、また、予め算出されたフィルタ係数を記憶しておき、その中からフィルタ係数同定文字列を用いて画素位置毎のフィルタ係数を選択するため、断層像生成の度にフィルタ係数を初めから算出する必要ない。従って、画質の高い断層像を実現しつつ、撮影から断層像表示までのTurn around timeを短縮することができる。
【００８９】
なお、本発明は上記実施形態そのままに限定されるものではなく、実施段階ではその要旨を逸脱しない範囲で構成要素を変形して具体化できる。具体的な変形例としては、例えば次のようなものがある。
【００９０】
本実施形態に係る各機能は、当該処理を実行するプログラムをワークステーション等のコンピュータにインストールし、これらをメモリ上で展開することによっても実現することができる。このとき、コンピュータに当該手法を実行させることのできるプログラムは、磁気ディスク（フロッピー（登録商標）ディスク、ハードディスクなど）、光ディスク（ＣＤ−ＲＯＭ、ＤＶＤなど）、半導体メモリなどの記録媒体に格納して頒布することも可能である。
【００９１】
また、上記実施形態に開示されている複数の構成要素の適宜な組み合わせにより、種々の発明を形成できる。例えば、実施形態に示される全構成要素から幾つかの構成要素を削除してもよい。さらに、異なる実施形態にわたる構成要素を適宜組み合わせてもよい。
【産業上の利用可能性】
【００９２】
以上本発明によれば、デジタル断層撮影において、画質の良い結像画像を高速で生成し表示することができるＸ線診断装置、画像処理装置、Ｘ線診断装置等における結像処理に用いられるフィルタ係数の算出プログラムを実現することができる。
【図面の簡単な説明】
【００９３】
【図１】図１は、本実施形態に係る乳房撮影用Ｘ線診断装置１の外観図を示している。
【図２】図２は、本実施形態に係る乳房撮影用Ｘ線診断装置１のブロック構成図を示している。
【図３】図３は、PI法による結像処理の流れを概念的に示した図である。
【図４】図４は、PI法による結像処理に用いられるフィルタ係数の計算処理の流れを概念的に示した図である。
【図５】図５は、等価ボケ関数を算出する手順を示した図である。
【図６】図６は、等価ボケ関数を算出する手順を示した図である。
【図７】図７は、等価ボケ関数を算出する手順を示した図である。
【図８】図８は、等価ボケ関数を算出する手順を示した図である。
【図９】図９は、あるフレームにおける、ある画素位置でのフィルタ係数を算出し、画像化して示した例である。
【図１０】図１０は、projection inversion法で算出したフィルタ係数により求めたフィルタ係数の特徴を説明するための図である。
【図１１】図１１は、フィルタ係数同定文字列の生成手法を説明するための図である。
【図１２】図１２は、乳房撮影用Ｘ線診断装置１のPI法による結像処理における動作の流れを示したフローチャートである。
【図１３】図１３は、スキャン軌道に関する条件を入力するためのインターフェースの一例を示した図である。
【図１４】図１４（ａ）、（ｂ）は、フィルタ係数採取法を説明するための図である。
【図１５】図１５は、従来技術を説明するための図である。
【符号の説明】
【００９４】
１…乳房撮影用Ｘ線診断装置、１０…アーム部、１１…Ｘ線制御部、１３…Ｘ線源装置、１５…圧迫板、１６…圧迫板制御部、１７…圧迫板駆動部、２０…平面検出器、２１…メモリ、２３…前処理部、２４…結像処理部、２５…画像処理部、２６…Ｄ／Ａ変換器、２７…表示部、３０…中央制御部、３１…操作部、３３…記憶部、３５…支柱部、３６…軸
【図１】
2008119457.tif 000019
【図２】
2008119457.tif 000020
【図３】
2008119457.tif 000021
【図４】
2008119457.tif 000022
【図５】
2008119457.tif 000023
【図６】
2008119457.tif 000024
【図７】
2008119457.tif 000025
【図８】
2008119457.tif 000026
【図９】
2008119457.tif 000027
【図１０】
2008119457.tif 000028
【図１１】
2008119457.tif 000029
【図１２】
2008119457.tif 000030
【図１３】
2008119457.tif 000031
【図１４】
2008119457.tif 000032
【図１５】
2008119457.tif 000033