(57)【要約】【課題】簡単なアルゴリズムで復元処理を行うことができる画像再構成方法を提供することを目的とする。【解決手段】ステップＳ１で不完全投影データを取得し、ステップＴ１で被検体の存在範囲を決定し、ステップＴ２で、被検体の存在範囲に発生する基底関数を定義することにより、当該存在範囲に制限して基底関数を決める。ステップＴ３，Ｔ４で不完全投影データと基底関数とをそれぞれフーリエ変換し、不完全投影データと同一になるように基底関数の自由度の数からなる連立方程式をガウスザイデル法で解き、個々の基底関数の係数を決定する。ステップＳ４で、決定された係数を有する基底関数に基づいて復元画像を求める。このように、先見情報として用いる条件も、被検体の存在範囲のみであり、当該存在範囲に制限して基底関数を決めるので、簡単なアルゴリズムで復元処理を行うことができる。【選択図】図２
2013183977.tif 000002
【特許請求の範囲】
【請求項１】
取得された不完全投影データから被検体の形状を復元した画像を再構成する画像再構成方法であって、
前記不完全投影データを取得する不完全投影データ取得工程と、
前記被検体が存在する範囲である被検体の存在範囲を決定する存在範囲決定工程と、
前記被検体の存在範囲に発生する信号のパターンを定義するパターン定義工程と、
前記不完全投影データと前記信号のパターンとをそれぞれ空間周波数上でフーリエ変換し、その空間周波数上の不完全投影データが存在する範囲で、不完全投影データと同一になるように信号のパターンの自由度の数からなる連立方程式を解き、それぞれの信号のパターンの係数を決定する係数決定工程と、
決定された前記係数を有する信号のパターンに基づいて復元画像を求める復元画像算出工程と
を備えることを特徴とする画像再構成方法。
【請求項２】
請求項１に記載の画像再構成方法において、
前記信号のパターンは、基底関数であることを特徴とする画像再構成方法。
【請求項３】
請求項１または請求項２に記載の画像再構成方法において、
前記不完全投影データは、投影角度が３６０°未満の範囲で取得されたデータであることを特徴とする画像再構成方法。
【請求項４】
請求項１から請求項３のいずれかに記載の画像再構成方法において、
存在範囲決定工程では、前記被検体の存在範囲を演算により決定することを特徴とする画像再構成方法。
【請求項５】
請求項１から請求項４のいずれかに記載の画像再構成方法において、
前記連立方程式を解く際に、ガウスザイデル法を用いて繰り返し近似計算を行うことを特徴とする画像再構成方法。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、取得された不完全投影データから被検体の形状を復元した画像を再構成する画像再構成方法に関する。
【背景技術】
【０００２】
不完全投影データが得られる場合としては、以下のケースが考えられる。例えば、Ｘ線ＣＴ装置のように投影角度が３６０°以上の範囲で取得されたデータであっても、被検体の体動や被検体の内部に不透明な物質が存在する場合には、投影データに欠落が生じて不完全投影データが取得されるケースがある。また、低被曝化を図るために、投影角度を３６０°未満の範囲で制限した結果、不完全投影データが取得されるケースがある。
【０００３】
投影角度を３６０°未満の範囲で制限する場合には、以下のような装置を用いて制限する。例えば、Ｘ線ＣＴ装置のガントリ内のＸ線管およびＸ線検出器を３６０°未満の範囲で被検体の回りを回転させることで投影角度を単に制限する。また、Ｘ線管とＸ線検出器とを互いに逆方向に平行移動させて撮影する、Ｘ線管およびＸ線検出器を同方向に平行移動させて撮影する、あるいはＸ線管またはＸ線検出器のいずれか一方のみを移動させて他方を固定させて撮影する、いわゆる「トモシンセシス」によって投影角度を制限する。あるいは、湾曲状に形成されたＣアームの一端がＸ線管を保持し、他端がＸ線検出器を保持し、Ｃアームに沿って回転させて撮影することによって投影角度を制限する。
【０００４】
このような不完全投影データから被検体の形状を復元した画像を再構成するのに、「凸射影法」と呼ばれる手法により復元化する（例えば、非特許文献１参照）。図５のデータのフローに示すように、先見情報として５つの条件と、撮影された不完全投影データとを用いて、凸射影法によりデータの欠落部分を推定することで、元の形状の画像に復元する。
【０００５】
ところで、通信の分野においては、演算の高速化を図るために、連立方程式を解く際に「ガウスザイデル法」と呼ばれる緩和法を用いる手法がある（例えば、特許文献１参照）。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】特開平６−１１９３６８号公報
【非特許文献】
【０００７】
【非特許文献１】工藤博幸、斎藤恒雄，「不完全投影データからのCT画像再構成実験」，Med. Image. Tech. Vol.9 No.1 March 1991 p.43-50
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、従来技術では、このような凸射影法という複雑なアルゴリズムを用いて、先見情報も５つの条件が必要で、先見情報として用いる情報が多いという問題がある。
【０００９】
この発明は、このような事情に鑑みてなされたものであって、簡単なアルゴリズムで復元処理を行うことができる画像再構成方法を提供することを目的とする。
【課題を解決するための手段】
【００１０】
この発明は、このような目的を達成するために、次のような構成をとる。
すなわち、この発明に係る画像再構成方法は、取得された不完全投影データから被検体の形状を復元した画像を再構成する画像再構成方法であって、前記不完全投影データを取得する不完全投影データ取得工程と、前記被検体が存在する範囲である被検体の存在範囲を決定する存在範囲決定工程と、前記被検体の存在範囲に発生する信号のパターンを定義するパターン定義工程と、前記不完全投影データと前記信号のパターンとをそれぞれ空間周波数上でフーリエ変換し、その空間周波数上の不完全投影データが存在する範囲で、不完全投影データと同一になるように信号のパターンの自由度の数からなる連立方程式を解き、それぞれの信号のパターンの係数を決定する係数決定工程と、決定された前記係数を有する信号のパターンに基づいて復元画像を求める復元画像算出工程とを備えることを特徴とするものである。
【００１１】
［作用・効果］この発明に係る画像再構成方法によれば、不完全投影データ取得工程で、不完全投影データを取得し、存在範囲決定工程で、被検体が存在する範囲である被検体の存在範囲を決定する。パターン定義工程で、被検体の存在範囲に発生する信号のパターンを定義することにより、当該存在範囲に制限して信号のパターンを決める。そして、係数決定工程で、不完全投影データと信号のパターンとをそれぞれ空間周波数上でフーリエ変換し、その空間周波数上の不完全投影データが存在する範囲で、不完全投影データと同一になるように信号のパターンの自由度の数からなる連立方程式を解き、それぞれの信号のパターンの係数を決定する。さらに、復元画像算出工程で、決定された係数を有する信号のパターンに基づいて復元画像を求める。このように、先見情報として用いる条件も、被検体の存在範囲のみである。また、当該存在範囲に制限して信号のパターンを決めるので、簡単なアルゴリズムで復元処理を行うことができる。
【００１２】
上述した信号のパターンの一例は、基底関数である。被検体の存在範囲に制限して、制限された存在範囲で個々の基底関数に展開する。個々の基底関数は互いに独立の関係（１次独立）であるので、個々の基底関数が決定し易く、最小限の数で個々の基底関数を定義することができる。もちろん、被検体の存在範囲に制限しても、基底関数以外として、個々の信号のパターンが１次従属であってもよい。ここで、基底関数の係数をＡｋｌとし、基底関数をＢｋｌとしたときに、ΣＡｋｌＢｋｌ＝０（ｋ，ｌは自由度、Σはｋ，ｌを動かしたときの総和）ならばＡ１１＝Ａ１２＝…＝Ａ１ｌ＝Ａ２１＝Ａ２２＝…＝Ａ２ｌ＝…＝Ａｋ１＝Ａｋ２＝…＝Ａｋｌ＝０を満たす（すなわちｋ，ｌを動かしたときに全てのＡｋｌ＝０を満たす）ときに個々の基底関数Ｂｋｌは１次独立であり、そうでなければ１次従属であると定義する。
【００１３】
上述した不完全投影データは、投影角度が３６０°未満の範囲で取得されたデータであってもよい。投影角度を３６０°未満の範囲で制限することで、低被曝化を図ることができ、そのようなデータであっても、欠落部分を推定して復元画像を求めることができる。もちろん、Ｘ線ＣＴ装置で得られたデータのように投影角度が３６０°以上の範囲に取得されたデータを不完全投影データとして用いてもよい。
【００１４】
上述した存在範囲決定工程では、被検体の存在範囲を演算により決定するのが好ましい。例えば他のモダリティ（modality）装置で得られた形態画像（アトラス画像）に基づいて被検体の存在範囲を演算により決定することで、当該存在範囲を自動で決定することができる。もちろん、経験則などのように当該存在範囲を手動で決定してもよい。
【００１５】
なお、当該存在範囲に制限して基底関数等に代表される信号のパターンを決めても、連立方程式は膨大な数となる。そこで、連立方程式を解く際に、ガウスザイデル法を用いて繰り返し近似計算を行うことで、演算の高速化を図ることができる。
【発明の効果】
【００１６】
この発明に係る画像再構成方法によれば、存在範囲決定工程で、被検体が存在する範囲である被検体の存在範囲を決定し、パターン定義工程で、被検体の存在範囲に発生する信号のパターンを定義することにより、当該存在範囲に制限して信号のパターンを決める。このように、先見情報として用いる条件も、被検体の存在範囲のみである。また、当該存在範囲に制限して信号のパターンを決めるので、簡単なアルゴリズムで復元処理を行うことができる。
【図面の簡単な説明】
【００１７】
【図１】実施例に係る断層撮影装置のブロック図である。
【図２】実施例に係る一連の画像再構成方法の流れを示したフローチャートである。
【図３】（ａ）は矩形基底関数の形状、（ｂ）はドーナツ形基底関数の形状である。
【図４】３つの基底関数による復元結果である。
【図５】従来における凸射影法のデータのフローである。
【実施例】
【００１８】
以下、図面を参照してこの発明の実施例を説明する。
図１は、実施例に係る断層撮影装置のブロック図である。本実施例では放射線としてＸ線を例に採って説明するとともに、被検体として人体を例に採って説明する。
【００１９】
断層撮影装置は、図１に示すように、被検体Ｍを載置する天板１と、その被検体Ｍに向けてＸ線を照射するＸ線管２と、被検体Ｍを透過したＸ線を検出するフラットパネル型Ｘ線検出器（以下、「ＦＰＤ」と略記する）３とを備えている。
【００２０】
断層撮影装置は、他に、天板１の昇降および水平移動を制御する天板制御部４や、ＦＰＤ３の走査を制御するＦＰＤ制御部５や、Ｘ線管２の管電圧や管電流を発生させる高電圧発生部６を有するＸ線管制御部７や、ＦＰＤ３から電荷信号であるＸ線検出信号をディジタル化して取り出すＡ／Ｄ変換器８や、Ａ／Ｄ変換器８から出力されたＸ線検出信号に基づいて種々の処理を行う画像処理部９や、これらの各構成部を統括するコントローラ１０や、処理された画像などを記憶するメモリ部１１や、オペレータが入力設定を行う入力部１２や、処理された画像などを表示するモニタ１３などを備えている。
【００２１】
天板制御部４は、天板１を水平移動させて被検体Ｍを撮像位置にまで収容したり、昇降、回転および水平移動させて被検体Ｍを所望の位置に設定したり、水平移動させながら撮像を行ったり、撮像終了後に水平移動させて撮像位置から退避させる制御などを行う。これらの制御は、モータやエンコーダ（図示省略）などからなる天板駆動機構（図示省略）を制御することで行う。
【００２２】
ＦＰＤ制御部５は、ＦＰＤ３を被検体Ｍの長手方向である体軸ｚ方向に沿って平行移動させる制御を行う。この制御は、ラックやピニオンやモータやエンコーダ（図示省略）などからなるＦＰＤ駆動機構（図示省略）を制御することで行う。
【００２３】
高電圧発生部６は、Ｘ線を照射させるための管電圧や管電流を発生してＸ線管２に与える。Ｘ線管制御部７は、Ｘ線管２をＦＰＤ３の平行移動とは逆方向に平行移動させる制御を行う。この制御は、支柱やネジ棒やモータやエンコーダ（図示省略）などからなるＸ線管駆動部（図示省略）を制御することで行う。
【００２４】
また、Ｘ線管制御部７は、Ｘ線管２側のコリメータ（図示省略）の照視野の設定の制御を行う。本実施例では、体軸ｚ方向に広がりを有するファンビーム状のＸ線を照射するようにコリメータを制御して照視野を設定する。
【００２５】
画像処理部９やコントローラ１０は、中央演算処理装置（ＣＰＵ）などで構成されており、メモリ部１１は、ＲＯＭ（Read-only Memory）やＲＡＭ（Random-Access Memory）などに代表される記憶媒体などで構成されている。また、入力部１２は、マウスやキーボードやジョイスティックやトラックボールやタッチパネルなどに代表されるポインティングデバイスで構成されている。
【００２６】
なお、各種の画像処理を行うためのプログラム等をＲＯＭなどに代表される記憶媒体に書き込んで記憶し、その記憶媒体からプログラム等を読み出して画像処理部９のＣＰＵが実行することでそのプログラムに応じた画像処理を行う。特に、不完全投影データ取得や存在範囲決定工程やパターン定義工程や係数決定や復元画像算出に関するプログラムを実行することで、そのプログラムに応じた不完全投影データ取得や存在範囲決定工程やパターン定義工程や係数決定や復元画像算出をそれぞれ行う。
【００２７】
メモリ部１１は、画像処理部９で処理された各々の画像を書き込んで記憶するように構成されている。特に、本実施例では、他のモダリティ装置で得られた形態画像（アトラス画像）をメモリ部１１に書き込んで記憶し、後述する被検体Ｍの存在範囲を演算により決定する際に形態画像を供する。ＦＰＤ制御部５やＸ線管制御部７も、画像処理部９やコントローラ１０と同様にＣＰＵなどで構成されている。
【００２８】
画像処理部９による画像再構成の流れについて、図２、図３を参照して説明する。図２は、実施例に係る一連の画像再構成方法の流れを示したフローチャートであり、図２（ａ）は先見情報のデータを含むフローであり、図２（ｂ）は図２（ａ）のステップＳ２における周波数空間の近似の具体的なフローチャートであり、図３（ａ）は、矩形基底関数の形状である、図３（ｂ）はドーナツ形基底関数の形状である。図２の画像再構成方法については、画像処理部９のＣＰＵがプログラムを実行することにより行う。
【００２９】
（ステップＳ１）不完全投影データ取得
図１に示すように、Ｘ線管２およびＦＰＤ３を互いに逆方向に平行移動させつつ、Ｘ線管２から体軸ｚ方向に広がりを有するファンビーム状のＸ線を照射することで、被検体Ｍに対して互いに異なる方向（投影角度）からそれぞれ照射されたＸ線ビームをＦＰＤ３が検出する。ＦＰＤ３が検出することによって、被検体Ｍに対して互いに異なる方向からそれぞれ照射されたＸ線ビームによる複数のＸ線画像を取得する。Ｘ線画像は、ＦＰＤ３の検出面に投影された投影データである。
【００３０】
図１に示す撮影は、いわゆる上述した「トモシンセシス」であり、これによって得られたデータは、投影角度が制限された不完全投影データである。図１では、Ｘ線管２とＦＰＤ３とを互いに逆方向に平行移動させて撮影する場合を図示したが、Ｘ線管２およびＦＰＤ３を同方向に平行移動させて撮影する、あるいはＸ線管２またはＦＰＤ３のいずれか一方のみを移動させて他方を固定させて撮影する場合に得られたデータに適用してもよい。なお、本実施例ではフラットパネル型Ｘ線検出器（ＦＰＤ）３をＸ線検出器として用いたがイメージインテンシファイアなどのＸ線検出器であってもよい。ＦＰＤ３で検出されたデータを画像処理部９に送り込み、画像処理部９はそのデータを不完全投影データとして取得する。このステップＳ１は、この発明における不完全投影データ取得工程に相当する。
【００３１】
（ステップＳ２）周波数空間での近似
先見情報として被検体Ｍが存在する範囲である被検体Ｍの存在範囲を用いて、周波数空間での近似を行う。図２（ａ）のステップＳ２における周波数空間での近似の具体的なフローは、図２（ｂ）のステップＴ２〜Ｔ４である。
【００３２】
（ステップＴ１）存在範囲決定
先ず、先見情報として被検体Ｍの存在範囲を決定するのに、本実施例では演算により行う。演算については、画像処理部９またはコントローラ１０のＣＰＵにより実行する。他のモダリティ装置で得られた形態画像（アトラス画像）に基づいて被検体Ｍの存在範囲を演算により決定する。例えばある骨の復元を行う場合には、骨のアトラス画像などから決める。なお、モダリティ装置については、形態画像を得る装置であれば、特に限定されず、Ｘ線に関する形態画像を取得する、例えばＸ線ＣＴ装置などに代表されるＸ線撮影装置であってもよいし、Ｘ線以外の形態画像を取得する、例えば磁気共鳴診断（MRI: magnetic resonance image）装置などに代表される装置であってもよい。
【００３３】
被検体Ｍの存在範囲を演算により決定するには、例えば形態画像の画素値を二値化して、二値化された画素値のうち“１”または“０”のいずれかを選択することにより行う。もちろん、二値化に限定されず多値化して、被検体Ｍが存在し得る画素値の閾値で区切ることにより当該存在範囲を決定してもよい。その他にも、存在範囲を簡略化すべく実際に二値化あるいは多値化された画像の不要領域を除去するためにモルフォロジー（morphology）演算のオープニング処理やエロージョン（erosion）処理を行ってもよい。オープニング処理やエロージョン処理については形態学の公知の演算手法であるので、その説明を省略する。
【００３４】
もちろん、経験則などのように当該存在範囲を手動で決定するために、例えば入力部１２に当該存在範囲を入力設定してもよい。また、入力による手動と、演算による自動とを両方組み合わせもよい。例えば、存在範囲の慨要の形を演算により決定することで、当該存在範囲の慨要の形を自動で決定して、当該存在範囲を簡略化するために不要領域を手動で除去してもよい。このステップＴ１は、この発明における存在範囲決定工程に相当する。
【００３５】
（ステップＴ２）基底関数の定義
被検体Ｍの存在範囲に発生する信号のパターンとして、本実施例では基底関数を定義する。本明細書では、被検体Ｍの存在範囲に信号が存在する基底関数として、図３（ａ）に示す矩形基底関数、円形基底関数および図３（ｂ）に示すドーナツ形基底関数をそれぞれ例に採って説明する。
【００３６】
［矩形基底関数］
図３（ａ）に示すようなＭ×Ｎの範囲に信号が存在する矩形の基底関数として、下記（１）式のように定義する。δはデルタ関数であり、デルタ関数により信号の範囲を制御している。また、ｋ，ｌは、基底関数ｅｋｌの自由度であり、ｉは虚数である。
【００３７】
【数１】
2013183977.tif 000003
【００３８】
［円形基底関数］
円形の基底関数として、Zernike多項式を用いて定義する。これはレンズ系の波面収差を表現するときに用いられる関数で、下記（２）式のように表される。なお、極座標（ρ，θ）に変換して使用する。また、ｍ，ｎは、基底関数Ｕｍｎの自由度である。
【００３９】
【数２】
2013183977.tif 000004
【００４０】
ここで、Ｒｍｎ（ρ）は下記（３）式のように表される。
【００４１】
【数３】
2013183977.tif 000005
【００４２】
上記（３）式中の！は階乗である（例えば（ｎ−ｓ）！＝（ｎ−ｓ）×（ｎ−ｓ−１）×…×２×１）。ｎ，ｍは、−ｎ≦ｍ≦ｎかつｎ−｜ｍ｜が偶数である。また、Ｒｍｎ（１）＝１であるので、実空間上でＲ０の範囲に制限するならば、その中で半径ｒの基底を表す場合、ｒ／Ｒ０に正規化して使用する。さらに、本明細書では、ｍ＜０になるとき、下記（４）式のようにｍが正の値の共役をとって基底関数を計算する。これは、実空間での画像の虚部が“０”になるという考えによる。
【００４３】
【数４】
2013183977.tif 000006
【００４４】
［ドーナツ形基底関数］
ドーナツ形基底関数は、図３（ａ）の矩形基底関数と同様に、図３（ｂ）に示すように外径Ｒ０と内径ｒ０とで決まるデルタ関数により、信号の範囲を制御して、極座標表示ｒ，Ｒで円形の基底関数を、下記（５）式のように定義する。ここでは極座標変換を離散化しておく必要がある。
【００４５】
【数５】
2013183977.tif 000007
【００４６】
ここで、上記（５）式中のｅ（ｉ２πｋｒ／Ｒ０）の項は、ｋ≧０になるときはｃｏｓ（２πｋｒ／Ｒ０）とし、ｋ＜０になるときｓｉｎ（２πｋｒ／Ｒ０）として計算する。これは、円形基底関数と同様に、実空間での画像の虚部を“０”にするためである。
【００４７】
［まとめ］
基底関数については、上述した矩形基底関数、円形基底関数およびドーナツ形基底関数には、もちろん限定されない。被検体Ｍの形状に合わせて適切な基底関数を選択することで精度の向上が行える。このステップＴ２は、この発明におけるパターン定義工程に相当する。
【００４８】
（ステップＴ３）フーリエ変換
ステップＳ１で取得された不完全投影データとステップＴ２で定義された基底関数とをそれぞれ空間周波数上でフーリエ変換する。具体的なフーリエ変換については、離散フーリエ変換や高速フーリエ変換など特に限定されない。その空間周波数上の不完全投影データが存在する範囲で、不完全投影データと同一になるように基底関数の自由度の数からなる連立方程式を解く。
【００４９】
先ず、基底関数の係数をＡｋｌとし、基底関数をＢｋｌとし、復元する画像（復元画像）をＦとする。そして、制限角投影でのデータが取得可能な範囲（空間周波数上の不完全投影データが存在する範囲）におけるＦとＢｋｌとの画像の差をＳとすると、下記（６）式のように表される。
【００５０】
【数６】
2013183977.tif 000008
【００５１】
このＳが最小となるＡｋｌを求めるために、下記（７）式のように展開して、最小値を求める。
【００５２】
【数７】
2013183977.tif 000009
【００５３】
つまり、下記（８）式を満たすＡｋｌを求める。
【００５４】
【数８】
2013183977.tif 000010
【００５５】
（ステップＴ４）ガウスザイデル法
上記（８）式は、ｋ×ｌ個の連立方程式を解くことであるが、このｋ，ｌは基底関数の自由度であるので、基底関数の形状に依存するものの制限する存在範囲の画素数に等しいことが多い。そのため、膨大な数の連立方程式を解くことになるので、ガウスザイデル法を用いて計算を行うのが好ましい。
【００５６】
ガウスザイデル法は、連立一次方程式を反復法で解く手法であり、下記（９）式を用いて繰り返し近似計算を行い、解を求める。本明細書での計算では、アルゴリズムを簡略化したヤコビ法による計算（下記（９）式の最右辺）で行う。式からわかるように対角行列が“０”になると収束しないので、基底関数の形状に注意する。
【００５７】
【数９】
2013183977.tif 000011
【００５８】
これらのステップＴ３，Ｔ４は、この発明における係数決定工程に相当する。これらのステップＴ１〜Ｔ３により、図２（ａ）のステップＳ２における周波数空間での近似を行い、図２（ａ）のフローに戻る。
【００５９】
（ステップＳ３）復元画像算出
ステップＴ３で決定された係数Ａを有する基底関数Ｂに基づいて復元画像Ｆを求める。すなわち、復元画像Ｆは、決定された係数Ａを用いて個々の基底関数Ｂの１次結合（すなわちＦ＝ΣＡ×Ｂ）で表されるので、復元画像を容易に求めることができる。このステップＳ３は、この発明における復元画像算出工程に相当する。
【００６０】
［復元結果］
基底関数を用いた復元結果について、図４を参照して説明する。図４は、３つの基底関数による復元結果である。
【００６１】
原画像として、円筒図形を制限角投影で撮影した場合に得られるであろうモデル画像（制限角投影モデル）を不完全投影データとして用いる。また、アトラス画像として、大腿骨のサジタル（Sagittal）画像に前処理として骨成分を抽出した画像を用いる。
【００６２】
先ず、円筒図形の制限角投影モデルを下記のように取得する。円筒図形をフーリエ変換し、周波数空間上で、ある角度の範囲で信号を残し、残りの部分を“０”にする。これを逆フーリエ変換することで、制限角投影モデルを得る。一方、大腿骨のサジタル画像では前処理として白黒変換をした後に適当な閾値で、骨成分のみの画像を得る。また、骨の中心が画像の中心になるように合わせておく。
【００６３】
入力画像として、制限角投影モデルを用いて、矩形基底関数、円形基底関数およびドーナツ形基底関数による復元処理をそれぞれ行い、図４に示すような復元結果を得た。
【００６４】
図４の実空間で示すように、矩形基底関数では円の形状に近づいたが、やや縦長となっており、制限範囲の窓の上下端に信号が張り付くような結果となっている。図４の実空間でも示すように、円形基底関数では矩形基底関数のときよりも円に近い成分ができたが、その外側で高輝度のパターンが発生した結果となっている。ドーナツ形基底関数では、図４の実空間でも示すように、円の上下部がややボケたものの、他の基底関数に比べ、アーティファクトもなく円の形状を最も復元することができたことが確認されている。
【００６５】
本実施例に係る画像再構成方法によれば、不完全投影データ取得工程（ステップＳ１）で、不完全投影データを取得し、存在範囲決定工程（ステップＴ１）で、被検体Ｍが存在する範囲である被検体Ｍの存在範囲を決定する。パターン定義工程（ステップＴ２）で、被検体Ｍの存在範囲に発生する信号のパターン（本実施例では基底関数）を定義することにより、当該存在範囲に制限して信号のパターン（基底関数）を決める。そして、係数決定工程（ステップＴ３，Ｔ４）で、不完全投影データと信号のパターン（基底関数）とをそれぞれ空間周波数上でフーリエ変換し、その空間周波数上の不完全投影データが存在する範囲で、不完全投影データと同一になるように信号のパターン（基底関数）の自由度の数からなる連立方程式を解き（本実施例ではガウスザイデル法を用いて連立方程式を解き）、それぞれの信号のパターン（基底関数）の係数を決定する。さらに、復元画像算出工程（ステップＳ４）で、決定された係数を有する信号のパターン（基底関数）に基づいて復元画像を求める。このように、先見情報として用いる条件も、被検体Ｍの存在範囲のみである。また、当該存在範囲に制限して信号のパターン（基底関数）を決めるので、簡単なアルゴリズムで復元処理を行うことができる。
【００６６】
本実施例では、信号のパターンとして基底関数を用いている。つまり、被検体Ｍの存在範囲に制限して、制限された存在範囲で個々の基底関数に展開する。個々の基底関数は互いに独立の関係（１次独立）であるので、個々の基底関数が決定し易く、最小限の数で個々の基底関数を定義することができる。
【００６７】
本実施例では、不完全投影データは、投影角度が３６０°未満（図１のトモシンセシスでは、例えば垂直な軸から±２０°の投影角度であるので、投影角度は４０°）の範囲で取得されたデータである。投影角度を３６０°未満の範囲で制限することで、低被曝化を図ることができ、そのようなデータであっても、欠落部分を推定して復元画像を求めることができる。
【００６８】
本実施例では、存在範囲決定工程（ステップＴ１）では、被検体Ｍの存在範囲を演算により決定するのが好ましい。例えば他のモダリティ装置で得られた形態画像（アトラス画像）に基づいて被検体Ｍの存在範囲を演算により決定することで、当該存在範囲を自動で決定することができる。
【００６９】
本実施例では、当該存在範囲に制限して基底関数等に代表される信号のパターンを決めても、連立方程式は膨大な数となる。そこで、連立方程式を解く際に、ガウスザイデル法を用いて繰り返し近似計算を行うことで、演算の高速化を図ることができる。
【００７０】
この発明は、上記実施形態に限られることはなく、下記のように変形実施することができる。
【００７１】
（１）上述した実施例では、放射線としてＸ線を例に採って説明したが、Ｘ線以外の放射線（例えばα線、β線、γ線など）や電磁波や光などのように、取得された不完全投影データから被検体の形状を復元した画像を再構成する画像再構成方法の手法であれば、特に限定されない。
【００７２】
（２）上述した実施例では、被検体として人体を用いたが、非破壊検査などに用いられう基板などであってもよい。
【００７３】
（３）上述した実施例では、信号のパターンとして１次独立の基底関数を例に採って説明したが、必ずしも基底関数に限定されない。被検体の存在範囲に制限しても、基底関数以外として、個々の信号のパターンが１次従属であってもよい。
【００７４】
（４）上述した実施例では、投影角度が３６０°未満の範囲で取得されたデータとして、図１のトモシンセシスで得られたデータを例に採って説明したが、Ｘ線ＣＴ装置のガントリ内のＸ線管およびＸ線検出器を３６０°未満の範囲で被検体の回りを回転させることで投影角度を単に制限して得られたデータや、湾曲状に形成されたＣアームの一端がＸ線管を保持し、他端がＸ線検出器を保持し、Ｃアームに沿って回転させて撮影することによって投影角度を制限して得られたデータであってもよい。また、上述したようにデータはＸ線に関する画像に限定されず、Ｘ線以外の放射線や電磁波や光などで得られた画像であってもよい。
【００７５】
（５）上述した実施例では、不完全投影データは、投影角度が３６０°未満の範囲で取得されたデータであったが、Ｘ線ＣＴ装置で得られたデータのように投影角度が３６０°以上の範囲に取得されたデータを不完全投影データとして用いてもよい。同じく、データはＸ線に関する画像に限定されず、Ｘ線以外の放射線や電磁波や光などで得られた画像であってもよい。
【００７６】
（６）上述した実施例では、存在範囲決定工程（ステップＴ１）では、被検体の存在範囲を演算により決定したが、上述したように経験則などのように当該存在範囲を手動で決定してもよい。また、上述したように、手動と自動とを両方組み合わせてもよい。
【００７７】
（７）上述した実施例では、ガウスザイデル法を用いて連立方程式を解いたが、ガウスザイデル法に限定されない。上述の特許文献１（特開平６−１１９３６８号公報）でも示唆されているＳＯＲ法に例示されているように、通常の連立方程式を解くのに用いられる手法であれば、特に限定されない。
【００７８】
（８）個々に基底関数を用いなくとも、複数の基底関数（例えば矩形基底関数、円形基底関数およびドーナツ形基底関数）を組み合わせて、最も最小値をとり得る係数を有した基底関数を選択してもよい。
【符号の説明】
【００７９】
９  …  画像処理部
Ａｋｌ  …  係数
Ｂｋｌ，ｅｋｌ，Ｕｍｎ  …  基底関数
ｋ，ｌ，ｍ，ｎ  …  自由度
Ｍ  …  被検体
【図１】
2013183977.tif 000012
【図２】
2013183977.tif 000013
【図３】
2013183977.tif 000014
【図５】
2013183977.tif 000015
【図４】
2013183977.tif 000016