画像再構成装置及びプログラム
【課題】種々の投影像の撮像態様に容易に対応でき、応用性を向上した画像再構成装置及びプログラムを提供する。
【解決手段】観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れ、観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、投影像を生成する確率的モデルを表す演算式を生成し、ベイズ推定を用いて、確率的モデルを表す演算式と、画像データとから観測対象の断面像を推定する画像再構成装置である。
【解決手段】観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れ、観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、投影像を生成する確率的モデルを表す演算式を生成し、ベイズ推定を用いて、確率的モデルを表す演算式と、画像データとから観測対象の断面像を推定する画像再構成装置である。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、画像再構成装置及びプログラムに係り、特に画質の向上に関する。
【背景技術】
【０００２】
観測対象のＸ線撮像画像を、観測対象の外部から様々な角度θで撮影し、撮影によって得られた複数の画像を用いて、観測対象の内部の画像を再構成する技術が知られている。この技術は、ＣＴ（Computed Tomography）と呼ばれ、医療や製品製造などの種々の場面で利用されている。
【０００３】
ＣＴ技術は、その画像を再構成する方法によって、解析的方法と、逐次近似法とに分類され、また逐次近似法はさらに、代数的方法と統計的方法とに分類される。
【０００４】
図７は、従来のＣＴ技術における画像の再構成法の例を表す説明図である。観測対象ｘをＣＴにより撮影した投影像ｙは、元の像ｘのラドン変換（Radon変換）Ｗの結果として得られる。すなわちｙ＝Ｗｘ。そこで理論的には、この投影像ｙを逆ラドン変換することにより（ｘ＝Ｗtｙ）、元の像が得られることとなる。しかしながら、現実には撮影角度θが離散的であるため、観測対象の中心部から外側へ向かうに連れて、投影像間のサンプリング間隔が広がり（Ａ）、像がぼけてしまう（Ｂ）。
【０００５】
そこでフィルタｆを用い、まず投影像ｙに適用した上でｆ＊ｙ、逆変換Ｗtを行う方法が提案されている。このようなフィルタｆとしては、例えば、Ramachandran-Lakshminarayananフィルタと呼ばれるものがある。この方法は広く知られているので、詳しい説明を省略するが、これにより得られた逆変換後の像は比較的元の像に近いものとなる。
【０００６】
特許文献１にはＣＴ画像に対して超解像処理を行う技術が開示されている。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２００５−０９５３２９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、上記従来の方法では、検出器に由来するノイズや、撮像や再構成時の回転軸ずれ、回転量のゆらぎ等の影響が考慮されず、再構成した像に様々なアーチファクト（artifact：実際には存在しない像）が発生することが知られている。
【０００９】
また逆ラドン変換には、再構成像全体に亘る積分操作が行われるので、情報の一部が欠損している場合、致命的なアーチファクトを発生させる場合がある。例えば、電子顕微鏡ＣＴ等では撮像対象の構造的な制約から、撮影可能な角度に制限があり、制限に係る角度の範囲では撮像ができないので、情報の欠損が生じる。また観測対象の一部に、他の部分とは放射線吸収率が大きく異なる部分があると、放射状のアーチファクトが生じる。
【００１０】
一方、フィルタｆを用いずに、予め定めた更新ベクトルｑ_kと、重みαkとを用いて、元の像の推定像ｘkを、逐次的にｘk＝ｘk-1＋αkｑkとして近似していく方法もある。この方法によれば、フィルタを用いる方法に比べ、情報の欠損に対してはアーチファクトの発生を抑制できるものの、他の要因に基づくアーチファクトを抑制できない。
【００１１】
さらに、これらフィルタｆを用いる方法や、逐次的に近似する方法では、投影像の撮像態様に応じてフィルタｆを設計しなければならないなど、応用性に乏しい。
【００１２】
本発明は上記実情に鑑みて為されたもので、種々の投影像の撮像態様に容易に対応でき、応用性を向上した画像再構成装置及びプログラムを提供することを、その目的の一つとする。
【００１３】
また本発明の他の目的の一つは、種々の要因で発生するアーチファクトを抑制するなど、画質を向上できる画像再構成装置及びプログラムを提供することである。
【課題を解決するための手段】
【００１４】
上記従来例の問題点を解決するための発明は、画像再構成装置であって、観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、を含むこととしたものである。
【００１５】
また本発明の一態様に係る画像再構成装置は、前記確率的モデルを表す演算式を生成する手段は、前記観測対象に基づいて投影像が得られるまでの過程に係る条件として、撮影時の外乱による未知パラメータを含む条件を表す演算式を生成する手段であり、当該未知パラメータを推定する手段をさらに含むこととしたものである。
さらに本発明の別の態様に係る画像再構成装置は、前記投影像は、観測対象を透過する電磁波または粒子線を検出することによって得たものであり前記確率的モデルを表す演算式を生成する手段は、前記観測対象の透過に伴う前記電磁波または粒子線の、波長に依存したエネルギーの減衰が、前記投影像にもたらす影響を表す変換行列を用いて、演算式である投影変換行列を生成することとしたものである。
【００１６】
さらに本発明の別の態様に係るプログラムは、コンピュータを、観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、として機能させることとしたものである。
【図面の簡単な説明】
【００１７】
【図１】本発明の実施の形態に係る画像再構成装置の例を表す構成ブロック図である。
【図２】本発明の実施の形態に係る画像再構成装置の一例に係る機能ブロック図である。
【図３】断面像から投影像を得る過程の一例を表す説明図である。
【図４】本発明の実施の形態に係る画像再構成装置の処理例を表すフローチャート図である。
【図５】本発明の実施の形態に係る画像再構成装置と、従来例の画像構成装置とによりそれぞれ再構成される画像の例を表す説明図である。
【図６】本発明の実施の形態に係る画像再構成装置と、従来例の画像構成装置とにより、それぞれ再構成された画像の評価の例を表す説明図である。
【図７】従来のＣＴ技術における画像の再構成法の例を表す説明図である。
【発明を実施するための形態】
【００１８】
本発明の実施の形態について図面を参照しながら説明する。本発明の実施の形態に係る画像再構成装置１は、図１に例示するように、制御部１１と、記憶部１２と、入力部１３と、出力部１４とを含んで構成されている。
【００１９】
ここで制御部１１は、ＣＰＵ（Central Processing Unit）等のプログラム制御デバイスであり、記憶部１２に格納されたプログラムに従って動作する。本実施の形態では、この制御部１１は、観測対象に対して互いに異なる方向から得た、当該観測対象を透過した複数の投影像に係る画像データを受け入れ、当該観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、投影像を生成する確率的モデルを表す演算式を生成し、ベイズ推定を用いて、当該確率的モデルを表す演算式と、受け入れた画像データとから観測対象の断面像を推定する処理を実行する。この制御部１１の詳しい処理の内容については後に述べる。
【００２０】
記憶部１２は、メモリデバイスやディスクデバイス等であり、制御部１１によって実行されるプログラムを保持している。このプログラムは、例えばＤＶＤ−ＲＯＭ等のコンピュータ可読な記録媒体に格納されて提供され、この記憶部１２に複写されたものであってもよい。また、この記憶部１２は、制御部１１のワークメモリとしても動作する。
【００２１】
入力部１３は、例えばシリアルまたはパラレル入力インタフェースを含み、ＣＴ装置にて撮像されて得られた画像データを受け入れて制御部１１に出力する。出力部１４は、ディスプレイやプリンタ等であり、制御部１１から入力される指示に従って画像データを表示出力する。
【００２２】
本実施の形態の画像再構成装置１が処理の対象とする複数の投影像は、観測対象に対して互いに異なる方向から、観測対象を透過し、または観測対象により散乱された電磁波や粒子線等を、複数の検出素子（channel）を配列した検出器（detector）で撮像して得られるもので、一例としては、観測対象を輪切りにする仮想的な平面上で、観測対象に対して、複数の角度方向からＸ線を照射する一般的なＣＴ装置による投影像である。つまり観測対象を上記平面で切ったときの断面像を、複数の角度から投影した複数の投影像を処理の対象としている。
【００２３】
本実施の形態の画像再構成装置１の制御部１１は、機能的には図２に例示するように、画像データ受入部２１と、輝度ベクトル抽出部２２と、変換行列生成部２３と、ベイズ推定部２４と、再構成部２５とを含んで構成されている。画像データ受入部２１は、入力部１３から複数（例えばＫ個）の投影像の画像データの入力を受け入れ、それぞれを記憶部１２に蓄積して格納する。ここで投影像の画像データはそれぞれ複数の画素値（検出素子ごとの出力値）を表すデータを含む。この画素値のデータは、例えば輝度の値を含む。
【００２４】
輝度ベクトル抽出部２２は、画像データ受入部２１が記憶部１２に格納したｋ番目（ｋ＝１，２，…，Ｋ）の画像データの各々について、予め定めた順に画素値のデータ（例えば輝度の値）を配列して、輝度値ベクトルｙ(k)（ｋ＝１，２，…，Ｋ）を得る。このＫ個の輝度値ベクトルの集合Ｄ＝｛ｙ(k)｝がサイノグラム（sinogram）に相当する。
【００２５】
変換行列生成部２３は、ベイズ推定に用いる事前分布に係る演算式を生成する。すなわち、元の画像（断面像）を表す輝度ベクトルをｘとし、投影像を得るためのラドン変換を表す投影変換行列をＷ(k)、撮像時の条件により生じる誤りをε(k)として、
【数１】
と表すことができる。すなわち、
【数２】
【００２６】
ここで、元の画像である断面像と投影像との対応関係、つまり元の画像の輝度ベクトルｘからサイノグラムＤが得られる確率を表す尤度関数ｐ（Ｄ｜ｘ）は、投影変換行列Ｗ(k)により一意に定まる部分を除いた部分（ε(k)＝Ｗ(k)Ｄ−ｘ）、誤り（ノイズ）ε(k)の生じる確率ｐ（ε(k)）に等しいから、
【数３】
を得る。ここでβは精度（分散の逆数）であり、ＩはＷ(k)と同じ次元の単位行列、Ｍはｙ(k)の次元、Ｎ（ｙ(k)｜ａ，ｂ）は平均ベクトルａ、分散行列をｂとする分布関数である。この分布関数は具体的にはガウス分布やポアソン分布など広く知られた分布関数のうちから、予め利用者が選択しておく。
【００２７】
また変換行列生成部２３は、投影変換行列Ｗ(k)を、少なくとも一つのアフィン変換を表す行列（平行移動、回転、縮小拡大等を表す行列）またはそれらの積（複数の場合）として、例えば次のように定める。すなわち観測対象から投影像が得られるまでの過程に係る条件として、ｋ番目の投影像を得たときの角度をθk、投影像のサイズ（輝度ベクトルの次元）をＭ、観測対象の元の像のサイズ（輝度ベクトルの次元）をＮ、撮像により点像が、投影像においていかなるサイズの円盤像に投影されるかを表す、円盤像の幅、つまり、点広がり関数（point spread function）の大きさγを用い、当該θkの回転を表す行列Ｒ（θk）、ξだけ平行移動をすることを表す行列をＴ（ξ）、ｒ倍の拡大縮小を表す行列をＭ（ｒ）として、まず行列
【数４】
を求める。この行列Ｆ(k)は、位置
【数５】
にある画素を所定の位置に移動したうえで、θkだけ回転させ、ｒ（ｒ＝Ｍ／Ｎ）倍して、
【数５】
を当該ｒ倍した位置まで再度移動することを意味する。これにより、次のように画素ベクトルの要素ｖjを変換すれば、投影像上の画素値ｕj(k)が得られる。
【数６】
【００２８】
そしてこのｕj(k)を用いた
【数７】
により、投影変換行列Ｗ(k)の要素を次のように表すことができる。
【数８】
【００２９】
また変換行列生成部２３は、観測対象の像に関する制約条件に基づく事前分布（つまり元の画像である断面像自体の評価関数）を表す演算式ｐ（ｘ）を生成する。具体的にこの分布は次のように定められる。すなわち観測対象の像（再構成されるべき像）においては画素間に相関があり、互いに隣接する画素同士では輝度の差が所定の閾値を超えないとする。このようなｐ（ｘ）としては、例えばガウス分布がある（次式）。
【００３０】
【数９】
【００３１】
なお、分散を表す固定行列（fixed matrix）Ｚxは、
【数１０】
であり、ここにｖiは位置ｉにある画素値（ここでは輝度の値）、ｒは画素間の距離、Ａは経験的に定められる係数である。変換行列生成部２３は、これら生成したｐ（Ｄ｜ｘ）、及びｐ（ｘ）に係る演算式（２），（４）式を、確率的モデルを表す演算式として出力する。また、このｐ（ｘ）としては、ポアソン分布等、他の分布を用いてもよい。
【００３２】
ベイズ推定部２４は、ベイズ推定を用いて事後分布ｐ（ｘ｜Ｄ）を演算する。これは、再構成されるべき元の画像である断面像自体の評価関数と、断面像と投影像との関係に基づく評価関数とのそれぞれの値を最大化するような画像を断面像の推定結果として探索することに相当する。具体的にはベイズの定理から、
【数１１】
であり、このうち分母のｐ（Ｄ）は、すべてのｘについて分子のｐ（ｘ）・ｐ（Ｄ｜ｘ）の和をとったものに他ならないから、予め規格化して「１」としておく。すると、サイノグラムＤが与えられたとき、像がｘである確率を与える事後分布ｐ（ｘ｜Ｄ）は結局、
【数１２】
となる。ここにμは平均のベクトルである。ベイズ推定部２４は、この（７）式に含まれる平均ベクトルμを演算する。この平均ベクトルμは、
【数１３】
である。ここに、Σは分散行列であり、
【数１４】
と表される。再構成部２５は、平均のベクトルμを像の輝度のベクトルとして演算し、出力する。
【００３３】
本実施の形態の画像再構成装置１は、以上の構成を備えてなり、次のように動作する。すなわち、既知である６４×６４画素の２次元画像Ｓに対し、その中心まわりに、予め定めた角度から、角度θk（ｋ＝１，２，…，Ｋ）だけ傾いた走査方向に沿って、当該走査方向に対し垂直な線分に射影した、それぞれｎ画素の複数の像（この場合は複数の一次元像）Ｙ1，Ｙ2，…，ＹKを得る（図３）。
【００３４】
本実施の形態の画像再構成装置１は、図４に示す処理を行い、これらの、それぞれｎ画素分の複数の一次元像を受け入れる。そして、各一次元像のそれぞれの画素の輝度値を配列した複数のｎ次元のベクトルｙ1，ｙ2，…，ｙkを得る（Ｓ１）。なお、ここでは輝度値としているが、カラーの場合、予め定めた色空間上の値（例えばＲ、Ｇ、Ｂの各値）のベクトルとしてもよい。
【００３５】
また画像再構成装置１は、観測対象から投影像が得られるまでの過程に係る条件として、ｋ番目の投影像を得たときの角度をθk、投影像のサイズ（輝度ベクトルの次元）をＭ、観測対象の元の像のサイズ（輝度ベクトルの次元）をＮ、撮像により点像が、投影像においていかなるサイズの円盤像に投影されるかを表す、円盤像の幅、つまり、点広がり関数（point spread function）の大きさγを用い、当該θkの回転を表す行列Ｒ（θk）、ξだけ平行移動をすることを表す行列をＴ（ξ）として、（３）式で表される行列Ｆ(k)を求める（Ｓ２）。
【００３６】
そして画像再構成装置１は、このＦ(k)により画素ベクトルの要素ｖjを変換して投影像上の画素値ｕj(k)を得て、このｕj(k)を用いて要素を定めた投影変換行列Ｗ(k)を得る（Ｓ３）。そしてこの行列を用いてベイズ推定に用いる尤度関数、ｐ（Ｄ｜ｘ）を（２）式のように求める（Ｓ４）。
【００３７】
また画像再構成装置１は、観測対象の像に関する制約条件に基づく事前分布を表す演算式ｐ（ｘ）を生成し（（４）式）（Ｓ５）、これら生成したｐ（Ｄ｜ｘ）、及びｐ（ｘ）に係る演算式（２），（４）式を確率的モデルを表す演算式とし、ベイズ推定の式（（６）式）を用いて、事後推定の結果である平均ベクトルμを（８）式により演算する（Ｓ６）。画像再構成装置１は、この平均ベクトルμの表す像を再構成した画像として出力する（Ｓ７）。
【００３８】
ここで平均ベクトルμは、再構成した画像に含まれるべき画素の輝度値（画素値）を予め定めた順に配列してベクトルとしたものであるから、当該順に従って元の画像の各画素の値を、当該平均ベクトルμにおける対応する成分の値から定める（例えばその値そのものとする）ことで、元の画像が得られる。
【００３９】
このような本実施の形態によると、元となる観測対象の断面像から投影像が得られるプロセスをそのまま観測画像の生成の確率的モデルとして記述することで、ベイズ統計学に基づく演算により、投影像から最尤断面像（観測対象の画像）を直接的に推定する。また、これより、投影変換行列を観測対象から投影像が得られるまでの過程に係る条件に基づいて設定することで、例えばヘリカルスキャン、コーンビーム、マルチＣＴ、位相ＣＴ、３次元ＣＴ等の種々の投影過程による投影像を対象として再構成像を得ることができる。
【００４０】
さらに本実施の形態の画像再構成装置１は、ラドン変換を表す投影変換行列Ｗ(k)を演算するにあたり、観測対象から投影像が得られるまでの過程に係る条件を表す演算式として、（３）式に定めた行列Ｆ(k)に代えて、次のような行列Ｆ(k)を用いることとしてもよい。すなわち、ｋ番目の投影像を得たときの角度をθk、投影像のサイズ（輝度ベクトルの次元）をＭ、観測対象の元の像のサイズ（輝度ベクトルの次元）をＮ、撮像により点像が、投影像においていかなるサイズの円盤像に投影されるかを表す、円盤像の幅、つまり、点広がり関数（point spread function）の大きさγを用い、当該θkの回転を表す行列Ｒ（θk）、ξだけ平行移動をすることを表す行列をＴ（ξ）、さらにｋ番目の投影像を得たときの位置ズレを表すシフトベクトルＳkによる平行移動Ｔ（Ｓk）を含め、
【数１５】
としてもよい。
【００４１】
この行列Ｆ(k)は、ｋ番目の投影像を撮像したときにＳkだけの位置ズレが生じたとしたものであるが、実際には、この位置ズレの量Ｓkは、外乱であるために測定が困難である。また、θk、γについても予め定めたものから、外乱のために実際の撮影時にはずれている可能性もある。すなわち、Ｓk，θk，γは、未知のパラメータであるということができる。
【００４２】
そこで本実施の形態の画像再構成装置１は、これら未知のパラメータを、いわゆるＥＭアルゴリズム（期待値最大化法）や、直接推定の方法を用いて推定する。ＥＭアルゴリズムを用いる場合、未知パラメータのセットをα
【数１６】
とおく。これにより、最終的に求めたい確率的モデルを表す演算式の一つの対数値
【数１７】
を考え、この期待値を最大化する。
【００４３】
すなわち画像再構成装置１は、当初は、この未知パラメータのセットの初期値αt＝α0を仮に定める。当初はＳk，γを０、θkはｋ番目の投影像を得たときのＸ線の飛跡方向を表す角度、つまり撮影時の設定角度としておいてもよいし、ランダムに決定してもよい。
【００４４】
そして画像再構成装置１は、未知パラメータのセットがαtであるとき、期待値を評価するステップ（Ｅステップ）では、（１１）式のうち、αtに係る項である
【数１８】
を演算する。以下この（１２）式を最大化することを考える。そしてこの（１２）式の期待値を最大化するステップ（Ｍステップ）として、（１２）式を最大化するαを求め、これを次の未知パラメータのセットαt+1とする。すなわち、
【数１９】
とする。
【００４５】
画像再構成装置１は、このαが収束するまで、つまり、ベクトルαt+1−αtの大きさが予め定めたしきい値を下回るまで、ＥステップとＭステップとを交互に繰り返し実行する。そして収束後のαが得られれば、当該αに含まれる未知パラメータの値を用いて、行列Ｆ(k)を求め、さらには投影変換行列Ｗ(k)を求めて、尤度関数ｐ（Ｄ｜ｘ）を定める。
【００４６】
また、この未知パラメータを推定する方法として、画像再構成装置は、周辺尤度を最大化するαを直接的に求めてもよい。すなわち周辺尤度を
【数２０】
として、これを最大化するαを求める
【数２１】
こととしてもよい。
【００４７】
またＣＴ撮影では、金属が含まれる生体組織の投影像を得た場合など、当該金属における高いＸ線吸収率により、投影像にデータが欠損した領域が生じる場合がある。これは例えば金属製のクラウン（歯科で用いる被せ物）が含まれる口腔を撮像した場合等にあり得るが、従来の方法では当該欠損した領域により生じるアーチファクト（いわゆるメタル・アーチファクト）が避け得なかった。本実施の形態の画像再構成装置１では、観測対象の像に関する制約条件に基づく事前分布ｐ（ｘ）として、当該欠損した領域を含んだ推定像のエッジを表現するイジングモデル等を用いることなどとすることで、このメタル・アーチファクトを軽減できる。
【００４８】
そのほか、この投影変換行列Ｗ(k)を求める際の行列Ｆ(k)に含めることのできる未知パラメータとして、投影像の撮像に用いた撮像素子等の特性（感度など）があり、観測対象の像に関する制約条件に基づく事前分布ｐ（ｘ）を定める際に、推定される断面に含まれるべき物体のＸ線吸収率などの物質の特性を表す情報等、種々のものを利用できる。これにより、予め分かっている観測対象の像（断面像）に関する情報を有効に利用して断面像の再構成の処理を行うことができる。また、投影変換行列Ｗ(k)を求める際の行列Ｆ(k)に得られる画像の解像度に係るパラメータを含むので、投影像よりも解像度の高い画像を再構成する、いわゆる超解像処理を行うことも可能となっている。
【００４９】
さらに、撮像角度が制限されている場合や、撮像対象の内部に、検出素子で検出する電磁波や粒子線源がある場合（蛍光体がある場合）、これら検出素子で検出する電磁波や粒子線を散乱ないし吸収する物質が撮像対象の内部に含まれている場合には、観測対象から投影像が得られるまでの過程に係る条件を表現した行列Ｆ(k)を生成するときに、それぞれの条件を表す行列を乗じたものとして生成すればよい。
【００５０】
また撮像時にゾーンプレートを用いている場合など、非線形な散乱が発生する場合は、当該非線形な散乱を線形に近似した行列Ｌを作成する。そしてこの行列Ｌを含む行列の積により、（３）または（１０）式に示したものと同様に、行列Ｆ(k)を得る。また隣接する検出素子間で出力値に相関がある場合も、その旨を表す行列を含む複数の行列の積により、行列Ｆ(k)を得る。
さらに、本実施の形態では、投影像を得るまでの過程がわかっていれば、いままでに得られているＣＴデータから再構成した画像を得ることができる。
【００５１】
さらに電子線やＸ線を利用したＣＴ等ではビームハードニング（Beam Hardening）アーチファクトの影響がある。これら検出素子で検出する電磁波や粒子線は一般に単波長ではなく、比較的波長の長い成分から短い成分まで種々の波長の電磁波や物質波（ド・ブロイ波）を含むものである。このビームハードニングは、これらＸ線等が撮像対象を透過する間に、比較的波長の長い成分がより吸収されて、透過に従って短い波長の成分が多くなり、結果としてエネルギーが高くなる現象をいう。
医療用にＣＴを用いる場合にこの現象が発生すると、体内に実在しない低吸収領域があるように見えてしまう場合がある。
そこで（３）または（１０）式に示したような、行列Ｆ(k)を行列の積のうちに、ビームハードニング（Ｘ線等のエネルギーの減衰）に伴うＣＴデータの変化を表す行列を含める。
【００５２】
具体的には、行列Ｆ(k)を、ｋ番目の投影像の撮像過程で生じる、投影対象の各部分における、波長に依存した電磁波または粒子線のエネルギーの減衰量を未知パラメータとして含む行列Ａ（ａλ1，ａλ2…）と、その他、撮像過程での変換を表す行列Ｍ（ｒ）、Ｔ（Ｓk）等との積により表すこととすればよい。ここで、ａλnは、予め定めた複数の波長（例えば、吸収率が互いに異なると考えられる各周波数帯から少なくとも一つずつの代表波長を選択することにより定めればよい）λ1，λ2，…λnの成分ごとの吸収係数（波長λi（ｉ＝１，２，…）の成分の減衰割合）を意味する。つまり、この行列を含んだ複数の行列の積で表される行列Ｆ(k)は、検出素子ごとの出力値である画素値を、多次元の吸収係数で表現するものとなる。これにより、波長帯ごとの吸収率を推定し、波長帯ごとの吸収率の相違によって発生するアーチファクトなどを軽減する。
このように、観測対象から投影像が得られるまでの過程に係る種々の条件を表した行列Ｆ(k)を生成することで、種々の撮像態様や、撮像対象の性質を反映した画像の再構成を行うことが可能となる。
【実施例】
【００５３】
本発明の実施例として、画像再構成装置１を、制御部１１であるＣＰＵとしてIntel Xeno 2.66GHz（４コア）を採用したコンピュータを用いて実現した例を次に示す。
【００５４】
図５は、元の６４×６４画素の二次元像（Ａ）について得た３２×８、３２×１６、３２×３２、３２×６４画素のサイノグラム（Ｒ）を得て、それぞれ単純に逆ラドン変換を行って得た再構成画像（Ｂ）、フィルタを用いた再構成画像（Ｃ）、逐次近似法による再構成画像（Ｄ）、本実施例の画像再構成装置１により得た再構成画像（Ｅ）を対比したものである。
図５において、各再構成画像（Ｂ，Ｃ，Ｄ，Ｅ）は上から順に、それぞれ３２×８、３２×１６、３２×３２、３２×６４画素のサイノグラムから得られる再構成画像を表す。
【００５５】
この図５の例において、再構成画像Ｃ，Ｄ，Ｅを表すベクトル（ｘにハット（＾）を付している）と、原画像ｘと、単純に逆ラドン変換を行って得た再構成画像を表すベクトル（ｘにチルド（〜）を付している）とを用いて演算される再現性の評価値ＩＳＮＲ（式（１６））
【数２２】
を得た結果を、図６に示す。図６において、横軸はサイノグラムのサイズ、縦軸はＩＳＮＲ値を示す。
【００５６】
図６に例示するように、本実施例の方法によると、フィルタを用いる例や逐次近似を用いる例に比べ、再現性を大幅に改善できる。
【符号の説明】
【００５７】
１ 画像再構成装置、１１ 制御部、１２ 記憶部、１３ 入力部、１４ 出力部、２１ 画像データ受入部、２２ 輝度ベクトル抽出部、２３ 変換行列生成部、２４ ベイズ推定部、２５ 再構成部。
【技術分野】
【０００１】
本発明は、画像再構成装置及びプログラムに係り、特に画質の向上に関する。
【背景技術】
【０００２】
観測対象のＸ線撮像画像を、観測対象の外部から様々な角度θで撮影し、撮影によって得られた複数の画像を用いて、観測対象の内部の画像を再構成する技術が知られている。この技術は、ＣＴ（Computed Tomography）と呼ばれ、医療や製品製造などの種々の場面で利用されている。
【０００３】
ＣＴ技術は、その画像を再構成する方法によって、解析的方法と、逐次近似法とに分類され、また逐次近似法はさらに、代数的方法と統計的方法とに分類される。
【０００４】
図７は、従来のＣＴ技術における画像の再構成法の例を表す説明図である。観測対象ｘをＣＴにより撮影した投影像ｙは、元の像ｘのラドン変換（Radon変換）Ｗの結果として得られる。すなわちｙ＝Ｗｘ。そこで理論的には、この投影像ｙを逆ラドン変換することにより（ｘ＝Ｗtｙ）、元の像が得られることとなる。しかしながら、現実には撮影角度θが離散的であるため、観測対象の中心部から外側へ向かうに連れて、投影像間のサンプリング間隔が広がり（Ａ）、像がぼけてしまう（Ｂ）。
【０００５】
そこでフィルタｆを用い、まず投影像ｙに適用した上でｆ＊ｙ、逆変換Ｗtを行う方法が提案されている。このようなフィルタｆとしては、例えば、Ramachandran-Lakshminarayananフィルタと呼ばれるものがある。この方法は広く知られているので、詳しい説明を省略するが、これにより得られた逆変換後の像は比較的元の像に近いものとなる。
【０００６】
特許文献１にはＣＴ画像に対して超解像処理を行う技術が開示されている。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２００５−０９５３２９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００８】
しかしながら、上記従来の方法では、検出器に由来するノイズや、撮像や再構成時の回転軸ずれ、回転量のゆらぎ等の影響が考慮されず、再構成した像に様々なアーチファクト（artifact：実際には存在しない像）が発生することが知られている。
【０００９】
また逆ラドン変換には、再構成像全体に亘る積分操作が行われるので、情報の一部が欠損している場合、致命的なアーチファクトを発生させる場合がある。例えば、電子顕微鏡ＣＴ等では撮像対象の構造的な制約から、撮影可能な角度に制限があり、制限に係る角度の範囲では撮像ができないので、情報の欠損が生じる。また観測対象の一部に、他の部分とは放射線吸収率が大きく異なる部分があると、放射状のアーチファクトが生じる。
【００１０】
一方、フィルタｆを用いずに、予め定めた更新ベクトルｑ_kと、重みαkとを用いて、元の像の推定像ｘkを、逐次的にｘk＝ｘk-1＋αkｑkとして近似していく方法もある。この方法によれば、フィルタを用いる方法に比べ、情報の欠損に対してはアーチファクトの発生を抑制できるものの、他の要因に基づくアーチファクトを抑制できない。
【００１１】
さらに、これらフィルタｆを用いる方法や、逐次的に近似する方法では、投影像の撮像態様に応じてフィルタｆを設計しなければならないなど、応用性に乏しい。
【００１２】
本発明は上記実情に鑑みて為されたもので、種々の投影像の撮像態様に容易に対応でき、応用性を向上した画像再構成装置及びプログラムを提供することを、その目的の一つとする。
【００１３】
また本発明の他の目的の一つは、種々の要因で発生するアーチファクトを抑制するなど、画質を向上できる画像再構成装置及びプログラムを提供することである。
【課題を解決するための手段】
【００１４】
上記従来例の問題点を解決するための発明は、画像再構成装置であって、観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、を含むこととしたものである。
【００１５】
また本発明の一態様に係る画像再構成装置は、前記確率的モデルを表す演算式を生成する手段は、前記観測対象に基づいて投影像が得られるまでの過程に係る条件として、撮影時の外乱による未知パラメータを含む条件を表す演算式を生成する手段であり、当該未知パラメータを推定する手段をさらに含むこととしたものである。
さらに本発明の別の態様に係る画像再構成装置は、前記投影像は、観測対象を透過する電磁波または粒子線を検出することによって得たものであり前記確率的モデルを表す演算式を生成する手段は、前記観測対象の透過に伴う前記電磁波または粒子線の、波長に依存したエネルギーの減衰が、前記投影像にもたらす影響を表す変換行列を用いて、演算式である投影変換行列を生成することとしたものである。
【００１６】
さらに本発明の別の態様に係るプログラムは、コンピュータを、観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、として機能させることとしたものである。
【図面の簡単な説明】
【００１７】
【図１】本発明の実施の形態に係る画像再構成装置の例を表す構成ブロック図である。
【図２】本発明の実施の形態に係る画像再構成装置の一例に係る機能ブロック図である。
【図３】断面像から投影像を得る過程の一例を表す説明図である。
【図４】本発明の実施の形態に係る画像再構成装置の処理例を表すフローチャート図である。
【図５】本発明の実施の形態に係る画像再構成装置と、従来例の画像構成装置とによりそれぞれ再構成される画像の例を表す説明図である。
【図６】本発明の実施の形態に係る画像再構成装置と、従来例の画像構成装置とにより、それぞれ再構成された画像の評価の例を表す説明図である。
【図７】従来のＣＴ技術における画像の再構成法の例を表す説明図である。
【発明を実施するための形態】
【００１８】
本発明の実施の形態について図面を参照しながら説明する。本発明の実施の形態に係る画像再構成装置１は、図１に例示するように、制御部１１と、記憶部１２と、入力部１３と、出力部１４とを含んで構成されている。
【００１９】
ここで制御部１１は、ＣＰＵ（Central Processing Unit）等のプログラム制御デバイスであり、記憶部１２に格納されたプログラムに従って動作する。本実施の形態では、この制御部１１は、観測対象に対して互いに異なる方向から得た、当該観測対象を透過した複数の投影像に係る画像データを受け入れ、当該観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、投影像を生成する確率的モデルを表す演算式を生成し、ベイズ推定を用いて、当該確率的モデルを表す演算式と、受け入れた画像データとから観測対象の断面像を推定する処理を実行する。この制御部１１の詳しい処理の内容については後に述べる。
【００２０】
記憶部１２は、メモリデバイスやディスクデバイス等であり、制御部１１によって実行されるプログラムを保持している。このプログラムは、例えばＤＶＤ−ＲＯＭ等のコンピュータ可読な記録媒体に格納されて提供され、この記憶部１２に複写されたものであってもよい。また、この記憶部１２は、制御部１１のワークメモリとしても動作する。
【００２１】
入力部１３は、例えばシリアルまたはパラレル入力インタフェースを含み、ＣＴ装置にて撮像されて得られた画像データを受け入れて制御部１１に出力する。出力部１４は、ディスプレイやプリンタ等であり、制御部１１から入力される指示に従って画像データを表示出力する。
【００２２】
本実施の形態の画像再構成装置１が処理の対象とする複数の投影像は、観測対象に対して互いに異なる方向から、観測対象を透過し、または観測対象により散乱された電磁波や粒子線等を、複数の検出素子（channel）を配列した検出器（detector）で撮像して得られるもので、一例としては、観測対象を輪切りにする仮想的な平面上で、観測対象に対して、複数の角度方向からＸ線を照射する一般的なＣＴ装置による投影像である。つまり観測対象を上記平面で切ったときの断面像を、複数の角度から投影した複数の投影像を処理の対象としている。
【００２３】
本実施の形態の画像再構成装置１の制御部１１は、機能的には図２に例示するように、画像データ受入部２１と、輝度ベクトル抽出部２２と、変換行列生成部２３と、ベイズ推定部２４と、再構成部２５とを含んで構成されている。画像データ受入部２１は、入力部１３から複数（例えばＫ個）の投影像の画像データの入力を受け入れ、それぞれを記憶部１２に蓄積して格納する。ここで投影像の画像データはそれぞれ複数の画素値（検出素子ごとの出力値）を表すデータを含む。この画素値のデータは、例えば輝度の値を含む。
【００２４】
輝度ベクトル抽出部２２は、画像データ受入部２１が記憶部１２に格納したｋ番目（ｋ＝１，２，…，Ｋ）の画像データの各々について、予め定めた順に画素値のデータ（例えば輝度の値）を配列して、輝度値ベクトルｙ(k)（ｋ＝１，２，…，Ｋ）を得る。このＫ個の輝度値ベクトルの集合Ｄ＝｛ｙ(k)｝がサイノグラム（sinogram）に相当する。
【００２５】
変換行列生成部２３は、ベイズ推定に用いる事前分布に係る演算式を生成する。すなわち、元の画像（断面像）を表す輝度ベクトルをｘとし、投影像を得るためのラドン変換を表す投影変換行列をＷ(k)、撮像時の条件により生じる誤りをε(k)として、
【数１】
と表すことができる。すなわち、
【数２】
【００２６】
ここで、元の画像である断面像と投影像との対応関係、つまり元の画像の輝度ベクトルｘからサイノグラムＤが得られる確率を表す尤度関数ｐ（Ｄ｜ｘ）は、投影変換行列Ｗ(k)により一意に定まる部分を除いた部分（ε(k)＝Ｗ(k)Ｄ−ｘ）、誤り（ノイズ）ε(k)の生じる確率ｐ（ε(k)）に等しいから、
【数３】
を得る。ここでβは精度（分散の逆数）であり、ＩはＷ(k)と同じ次元の単位行列、Ｍはｙ(k)の次元、Ｎ（ｙ(k)｜ａ，ｂ）は平均ベクトルａ、分散行列をｂとする分布関数である。この分布関数は具体的にはガウス分布やポアソン分布など広く知られた分布関数のうちから、予め利用者が選択しておく。
【００２７】
また変換行列生成部２３は、投影変換行列Ｗ(k)を、少なくとも一つのアフィン変換を表す行列（平行移動、回転、縮小拡大等を表す行列）またはそれらの積（複数の場合）として、例えば次のように定める。すなわち観測対象から投影像が得られるまでの過程に係る条件として、ｋ番目の投影像を得たときの角度をθk、投影像のサイズ（輝度ベクトルの次元）をＭ、観測対象の元の像のサイズ（輝度ベクトルの次元）をＮ、撮像により点像が、投影像においていかなるサイズの円盤像に投影されるかを表す、円盤像の幅、つまり、点広がり関数（point spread function）の大きさγを用い、当該θkの回転を表す行列Ｒ（θk）、ξだけ平行移動をすることを表す行列をＴ（ξ）、ｒ倍の拡大縮小を表す行列をＭ（ｒ）として、まず行列
【数４】
を求める。この行列Ｆ(k)は、位置
【数５】
にある画素を所定の位置に移動したうえで、θkだけ回転させ、ｒ（ｒ＝Ｍ／Ｎ）倍して、
【数５】
を当該ｒ倍した位置まで再度移動することを意味する。これにより、次のように画素ベクトルの要素ｖjを変換すれば、投影像上の画素値ｕj(k)が得られる。
【数６】
【００２８】
そしてこのｕj(k)を用いた
【数７】
により、投影変換行列Ｗ(k)の要素を次のように表すことができる。
【数８】
【００２９】
また変換行列生成部２３は、観測対象の像に関する制約条件に基づく事前分布（つまり元の画像である断面像自体の評価関数）を表す演算式ｐ（ｘ）を生成する。具体的にこの分布は次のように定められる。すなわち観測対象の像（再構成されるべき像）においては画素間に相関があり、互いに隣接する画素同士では輝度の差が所定の閾値を超えないとする。このようなｐ（ｘ）としては、例えばガウス分布がある（次式）。
【００３０】
【数９】
【００３１】
なお、分散を表す固定行列（fixed matrix）Ｚxは、
【数１０】
であり、ここにｖiは位置ｉにある画素値（ここでは輝度の値）、ｒは画素間の距離、Ａは経験的に定められる係数である。変換行列生成部２３は、これら生成したｐ（Ｄ｜ｘ）、及びｐ（ｘ）に係る演算式（２），（４）式を、確率的モデルを表す演算式として出力する。また、このｐ（ｘ）としては、ポアソン分布等、他の分布を用いてもよい。
【００３２】
ベイズ推定部２４は、ベイズ推定を用いて事後分布ｐ（ｘ｜Ｄ）を演算する。これは、再構成されるべき元の画像である断面像自体の評価関数と、断面像と投影像との関係に基づく評価関数とのそれぞれの値を最大化するような画像を断面像の推定結果として探索することに相当する。具体的にはベイズの定理から、
【数１１】
であり、このうち分母のｐ（Ｄ）は、すべてのｘについて分子のｐ（ｘ）・ｐ（Ｄ｜ｘ）の和をとったものに他ならないから、予め規格化して「１」としておく。すると、サイノグラムＤが与えられたとき、像がｘである確率を与える事後分布ｐ（ｘ｜Ｄ）は結局、
【数１２】
となる。ここにμは平均のベクトルである。ベイズ推定部２４は、この（７）式に含まれる平均ベクトルμを演算する。この平均ベクトルμは、
【数１３】
である。ここに、Σは分散行列であり、
【数１４】
と表される。再構成部２５は、平均のベクトルμを像の輝度のベクトルとして演算し、出力する。
【００３３】
本実施の形態の画像再構成装置１は、以上の構成を備えてなり、次のように動作する。すなわち、既知である６４×６４画素の２次元画像Ｓに対し、その中心まわりに、予め定めた角度から、角度θk（ｋ＝１，２，…，Ｋ）だけ傾いた走査方向に沿って、当該走査方向に対し垂直な線分に射影した、それぞれｎ画素の複数の像（この場合は複数の一次元像）Ｙ1，Ｙ2，…，ＹKを得る（図３）。
【００３４】
本実施の形態の画像再構成装置１は、図４に示す処理を行い、これらの、それぞれｎ画素分の複数の一次元像を受け入れる。そして、各一次元像のそれぞれの画素の輝度値を配列した複数のｎ次元のベクトルｙ1，ｙ2，…，ｙkを得る（Ｓ１）。なお、ここでは輝度値としているが、カラーの場合、予め定めた色空間上の値（例えばＲ、Ｇ、Ｂの各値）のベクトルとしてもよい。
【００３５】
また画像再構成装置１は、観測対象から投影像が得られるまでの過程に係る条件として、ｋ番目の投影像を得たときの角度をθk、投影像のサイズ（輝度ベクトルの次元）をＭ、観測対象の元の像のサイズ（輝度ベクトルの次元）をＮ、撮像により点像が、投影像においていかなるサイズの円盤像に投影されるかを表す、円盤像の幅、つまり、点広がり関数（point spread function）の大きさγを用い、当該θkの回転を表す行列Ｒ（θk）、ξだけ平行移動をすることを表す行列をＴ（ξ）として、（３）式で表される行列Ｆ(k)を求める（Ｓ２）。
【００３６】
そして画像再構成装置１は、このＦ(k)により画素ベクトルの要素ｖjを変換して投影像上の画素値ｕj(k)を得て、このｕj(k)を用いて要素を定めた投影変換行列Ｗ(k)を得る（Ｓ３）。そしてこの行列を用いてベイズ推定に用いる尤度関数、ｐ（Ｄ｜ｘ）を（２）式のように求める（Ｓ４）。
【００３７】
また画像再構成装置１は、観測対象の像に関する制約条件に基づく事前分布を表す演算式ｐ（ｘ）を生成し（（４）式）（Ｓ５）、これら生成したｐ（Ｄ｜ｘ）、及びｐ（ｘ）に係る演算式（２），（４）式を確率的モデルを表す演算式とし、ベイズ推定の式（（６）式）を用いて、事後推定の結果である平均ベクトルμを（８）式により演算する（Ｓ６）。画像再構成装置１は、この平均ベクトルμの表す像を再構成した画像として出力する（Ｓ７）。
【００３８】
ここで平均ベクトルμは、再構成した画像に含まれるべき画素の輝度値（画素値）を予め定めた順に配列してベクトルとしたものであるから、当該順に従って元の画像の各画素の値を、当該平均ベクトルμにおける対応する成分の値から定める（例えばその値そのものとする）ことで、元の画像が得られる。
【００３９】
このような本実施の形態によると、元となる観測対象の断面像から投影像が得られるプロセスをそのまま観測画像の生成の確率的モデルとして記述することで、ベイズ統計学に基づく演算により、投影像から最尤断面像（観測対象の画像）を直接的に推定する。また、これより、投影変換行列を観測対象から投影像が得られるまでの過程に係る条件に基づいて設定することで、例えばヘリカルスキャン、コーンビーム、マルチＣＴ、位相ＣＴ、３次元ＣＴ等の種々の投影過程による投影像を対象として再構成像を得ることができる。
【００４０】
さらに本実施の形態の画像再構成装置１は、ラドン変換を表す投影変換行列Ｗ(k)を演算するにあたり、観測対象から投影像が得られるまでの過程に係る条件を表す演算式として、（３）式に定めた行列Ｆ(k)に代えて、次のような行列Ｆ(k)を用いることとしてもよい。すなわち、ｋ番目の投影像を得たときの角度をθk、投影像のサイズ（輝度ベクトルの次元）をＭ、観測対象の元の像のサイズ（輝度ベクトルの次元）をＮ、撮像により点像が、投影像においていかなるサイズの円盤像に投影されるかを表す、円盤像の幅、つまり、点広がり関数（point spread function）の大きさγを用い、当該θkの回転を表す行列Ｒ（θk）、ξだけ平行移動をすることを表す行列をＴ（ξ）、さらにｋ番目の投影像を得たときの位置ズレを表すシフトベクトルＳkによる平行移動Ｔ（Ｓk）を含め、
【数１５】
としてもよい。
【００４１】
この行列Ｆ(k)は、ｋ番目の投影像を撮像したときにＳkだけの位置ズレが生じたとしたものであるが、実際には、この位置ズレの量Ｓkは、外乱であるために測定が困難である。また、θk、γについても予め定めたものから、外乱のために実際の撮影時にはずれている可能性もある。すなわち、Ｓk，θk，γは、未知のパラメータであるということができる。
【００４２】
そこで本実施の形態の画像再構成装置１は、これら未知のパラメータを、いわゆるＥＭアルゴリズム（期待値最大化法）や、直接推定の方法を用いて推定する。ＥＭアルゴリズムを用いる場合、未知パラメータのセットをα
【数１６】
とおく。これにより、最終的に求めたい確率的モデルを表す演算式の一つの対数値
【数１７】
を考え、この期待値を最大化する。
【００４３】
すなわち画像再構成装置１は、当初は、この未知パラメータのセットの初期値αt＝α0を仮に定める。当初はＳk，γを０、θkはｋ番目の投影像を得たときのＸ線の飛跡方向を表す角度、つまり撮影時の設定角度としておいてもよいし、ランダムに決定してもよい。
【００４４】
そして画像再構成装置１は、未知パラメータのセットがαtであるとき、期待値を評価するステップ（Ｅステップ）では、（１１）式のうち、αtに係る項である
【数１８】
を演算する。以下この（１２）式を最大化することを考える。そしてこの（１２）式の期待値を最大化するステップ（Ｍステップ）として、（１２）式を最大化するαを求め、これを次の未知パラメータのセットαt+1とする。すなわち、
【数１９】
とする。
【００４５】
画像再構成装置１は、このαが収束するまで、つまり、ベクトルαt+1−αtの大きさが予め定めたしきい値を下回るまで、ＥステップとＭステップとを交互に繰り返し実行する。そして収束後のαが得られれば、当該αに含まれる未知パラメータの値を用いて、行列Ｆ(k)を求め、さらには投影変換行列Ｗ(k)を求めて、尤度関数ｐ（Ｄ｜ｘ）を定める。
【００４６】
また、この未知パラメータを推定する方法として、画像再構成装置は、周辺尤度を最大化するαを直接的に求めてもよい。すなわち周辺尤度を
【数２０】
として、これを最大化するαを求める
【数２１】
こととしてもよい。
【００４７】
またＣＴ撮影では、金属が含まれる生体組織の投影像を得た場合など、当該金属における高いＸ線吸収率により、投影像にデータが欠損した領域が生じる場合がある。これは例えば金属製のクラウン（歯科で用いる被せ物）が含まれる口腔を撮像した場合等にあり得るが、従来の方法では当該欠損した領域により生じるアーチファクト（いわゆるメタル・アーチファクト）が避け得なかった。本実施の形態の画像再構成装置１では、観測対象の像に関する制約条件に基づく事前分布ｐ（ｘ）として、当該欠損した領域を含んだ推定像のエッジを表現するイジングモデル等を用いることなどとすることで、このメタル・アーチファクトを軽減できる。
【００４８】
そのほか、この投影変換行列Ｗ(k)を求める際の行列Ｆ(k)に含めることのできる未知パラメータとして、投影像の撮像に用いた撮像素子等の特性（感度など）があり、観測対象の像に関する制約条件に基づく事前分布ｐ（ｘ）を定める際に、推定される断面に含まれるべき物体のＸ線吸収率などの物質の特性を表す情報等、種々のものを利用できる。これにより、予め分かっている観測対象の像（断面像）に関する情報を有効に利用して断面像の再構成の処理を行うことができる。また、投影変換行列Ｗ(k)を求める際の行列Ｆ(k)に得られる画像の解像度に係るパラメータを含むので、投影像よりも解像度の高い画像を再構成する、いわゆる超解像処理を行うことも可能となっている。
【００４９】
さらに、撮像角度が制限されている場合や、撮像対象の内部に、検出素子で検出する電磁波や粒子線源がある場合（蛍光体がある場合）、これら検出素子で検出する電磁波や粒子線を散乱ないし吸収する物質が撮像対象の内部に含まれている場合には、観測対象から投影像が得られるまでの過程に係る条件を表現した行列Ｆ(k)を生成するときに、それぞれの条件を表す行列を乗じたものとして生成すればよい。
【００５０】
また撮像時にゾーンプレートを用いている場合など、非線形な散乱が発生する場合は、当該非線形な散乱を線形に近似した行列Ｌを作成する。そしてこの行列Ｌを含む行列の積により、（３）または（１０）式に示したものと同様に、行列Ｆ(k)を得る。また隣接する検出素子間で出力値に相関がある場合も、その旨を表す行列を含む複数の行列の積により、行列Ｆ(k)を得る。
さらに、本実施の形態では、投影像を得るまでの過程がわかっていれば、いままでに得られているＣＴデータから再構成した画像を得ることができる。
【００５１】
さらに電子線やＸ線を利用したＣＴ等ではビームハードニング（Beam Hardening）アーチファクトの影響がある。これら検出素子で検出する電磁波や粒子線は一般に単波長ではなく、比較的波長の長い成分から短い成分まで種々の波長の電磁波や物質波（ド・ブロイ波）を含むものである。このビームハードニングは、これらＸ線等が撮像対象を透過する間に、比較的波長の長い成分がより吸収されて、透過に従って短い波長の成分が多くなり、結果としてエネルギーが高くなる現象をいう。
医療用にＣＴを用いる場合にこの現象が発生すると、体内に実在しない低吸収領域があるように見えてしまう場合がある。
そこで（３）または（１０）式に示したような、行列Ｆ(k)を行列の積のうちに、ビームハードニング（Ｘ線等のエネルギーの減衰）に伴うＣＴデータの変化を表す行列を含める。
【００５２】
具体的には、行列Ｆ(k)を、ｋ番目の投影像の撮像過程で生じる、投影対象の各部分における、波長に依存した電磁波または粒子線のエネルギーの減衰量を未知パラメータとして含む行列Ａ（ａλ1，ａλ2…）と、その他、撮像過程での変換を表す行列Ｍ（ｒ）、Ｔ（Ｓk）等との積により表すこととすればよい。ここで、ａλnは、予め定めた複数の波長（例えば、吸収率が互いに異なると考えられる各周波数帯から少なくとも一つずつの代表波長を選択することにより定めればよい）λ1，λ2，…λnの成分ごとの吸収係数（波長λi（ｉ＝１，２，…）の成分の減衰割合）を意味する。つまり、この行列を含んだ複数の行列の積で表される行列Ｆ(k)は、検出素子ごとの出力値である画素値を、多次元の吸収係数で表現するものとなる。これにより、波長帯ごとの吸収率を推定し、波長帯ごとの吸収率の相違によって発生するアーチファクトなどを軽減する。
このように、観測対象から投影像が得られるまでの過程に係る種々の条件を表した行列Ｆ(k)を生成することで、種々の撮像態様や、撮像対象の性質を反映した画像の再構成を行うことが可能となる。
【実施例】
【００５３】
本発明の実施例として、画像再構成装置１を、制御部１１であるＣＰＵとしてIntel Xeno 2.66GHz（４コア）を採用したコンピュータを用いて実現した例を次に示す。
【００５４】
図５は、元の６４×６４画素の二次元像（Ａ）について得た３２×８、３２×１６、３２×３２、３２×６４画素のサイノグラム（Ｒ）を得て、それぞれ単純に逆ラドン変換を行って得た再構成画像（Ｂ）、フィルタを用いた再構成画像（Ｃ）、逐次近似法による再構成画像（Ｄ）、本実施例の画像再構成装置１により得た再構成画像（Ｅ）を対比したものである。
図５において、各再構成画像（Ｂ，Ｃ，Ｄ，Ｅ）は上から順に、それぞれ３２×８、３２×１６、３２×３２、３２×６４画素のサイノグラムから得られる再構成画像を表す。
【００５５】
この図５の例において、再構成画像Ｃ，Ｄ，Ｅを表すベクトル（ｘにハット（＾）を付している）と、原画像ｘと、単純に逆ラドン変換を行って得た再構成画像を表すベクトル（ｘにチルド（〜）を付している）とを用いて演算される再現性の評価値ＩＳＮＲ（式（１６））
【数２２】
を得た結果を、図６に示す。図６において、横軸はサイノグラムのサイズ、縦軸はＩＳＮＲ値を示す。
【００５６】
図６に例示するように、本実施例の方法によると、フィルタを用いる例や逐次近似を用いる例に比べ、再現性を大幅に改善できる。
【符号の説明】
【００５７】
１ 画像再構成装置、１１ 制御部、１２ 記憶部、１３ 入力部、１４ 出力部、２１ 画像データ受入部、２２ 輝度ベクトル抽出部、２３ 変換行列生成部、２４ ベイズ推定部、２５ 再構成部。
【特許請求の範囲】
【請求項１】
観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、
前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、
ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、
を含む画像再構成装置。
【請求項２】
請求項１記載の画像再構成装置であって、
前記確率的モデルを表す演算式を生成する手段は、前記観測対象に基づいて投影像が得られるまでの過程に係る条件として、撮影時の外乱による未知パラメータを含む条件を表す演算式を生成する手段であり、
当該未知パラメータを推定する手段をさらに含む画像再構成装置。
【請求項３】
請求項１または２記載の画像再構成装置であって、
前記投影像は、観測対象を透過する電磁波または粒子線を検出することによって得たものであり、
前記確率的モデルを表す演算式を生成する手段は、前記観測対象の透過に伴う前記電磁波または粒子線の、波長に依存したエネルギーの減衰が、前記投影像にもたらす影響を表す変換行列を用いて、演算式である投影変換行列を生成することを特徴とする画像再構成装置。
【請求項４】
コンピュータを、
観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、
前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、
ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、
として機能させるプログラム。
【請求項１】
観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、
前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、
ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、
を含む画像再構成装置。
【請求項２】
請求項１記載の画像再構成装置であって、
前記確率的モデルを表す演算式を生成する手段は、前記観測対象に基づいて投影像が得られるまでの過程に係る条件として、撮影時の外乱による未知パラメータを含む条件を表す演算式を生成する手段であり、
当該未知パラメータを推定する手段をさらに含む画像再構成装置。
【請求項３】
請求項１または２記載の画像再構成装置であって、
前記投影像は、観測対象を透過する電磁波または粒子線を検出することによって得たものであり、
前記確率的モデルを表す演算式を生成する手段は、前記観測対象の透過に伴う前記電磁波または粒子線の、波長に依存したエネルギーの減衰が、前記投影像にもたらす影響を表す変換行列を用いて、演算式である投影変換行列を生成することを特徴とする画像再構成装置。
【請求項４】
コンピュータを、
観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、
前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、
ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、
として機能させるプログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【公開番号】特開２０１３−５８４０（Ｐ２０１３−５８４０Ａ）
【公開日】平成２５年１月１０日（２０１３．１．１０）
【国際特許分類】
【出願番号】特願２０１１−１３８８３３（Ｐ２０１１−１３８８３３）
【出願日】平成２３年６月２２日（２０１１．６．２２）
【出願人】（５０４１３７９１２）国立大学法人 東京大学 (1,942)
【Ｆターム（参考）】
【公開日】平成２５年１月１０日（２０１３．１．１０）
【国際特許分類】
【出願日】平成２３年６月２２日（２０１１．６．２２）
【出願人】（５０４１３７９１２）国立大学法人 東京大学 (1,942)
【Ｆターム（参考）】
[ Back to top ]