"観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、を含む画像再構成装置。請求項1記載の画像再構成装置であって、前記確率的モデルを表す演算式を生成する手段は、前記観測対象に基づいて投影像が得られるまでの過程に係る条件として、撮影時の外乱による未知パラメータを含む条件を表す演算式を生成する手段であり、当該未知パラメータを推定する手段をさらに含む画像再構成装置。請求項1または2記載の画像再構成装置であって、前記投影像は、観測対象を透過する電磁波または粒子線を検出することによって得たものであり、前記確率的モデルを表す演算式を生成する手段は、前記観測対象の透過に伴う前記電磁波または粒子線の、波長に依存したエネルギーの減衰が、前記投影像にもたらす影響を表す変換行列を用いて、演算式である投影変換行列を生成することを特徴とする画像再構成装置。コンピュータを、観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、として機能させるプログラム。しかしながら、上記従来の方法では、検出器に由来するノイズや、撮像や再構成時の回転軸ずれ、回転量のゆらぎ等の影響が考慮されず、再構成した像に様々なアーチファクトartifact:実際には存在しない像が発生することが知られている。また逆ラドン変換には、再構成像全体に亘る積分操作が行われるので、情報の一部が欠損している場合、致命的なアーチファクトを発生させる場合がある。例えば、電子顕微鏡CT等では撮像対象の構造的な制約から、撮影可能な角度に制限があり、制限に係る角度の範囲では撮像ができないので、情報の欠損が生じる。また観測対象の一部に、他の部分とは放射線吸収率が大きく異なる部分があると、放射状のアーチファクトが生じる。一方、フィルタfを用いずに、予め定めた更新ベクトルqkと、重みαkとを用いて、元の像の推定像xkを、逐次的にxk=xk1αkqkとして近似していく方法もある。この方法によれば、フィルタを用いる方法に比べ、情報の欠損に対してはアーチファクトの発生を抑制できるものの、他の要因に基づくアーチファクトを抑制できない。さらに、これらフィルタfを用いる方法や、逐次的に近似する方法では、投影像の撮像態様に応じてフィルタfを設計しなければならないなど、応用性に乏しい。本発明は上記実情に鑑みて為されたもので、種々の投影像の撮像態様に容易に対応でき、応用性を向上した画像再構成装置及びプログラムを提供することを、その目的の一つとする。また本発明の他の目的の一つは、種々の要因で発生するアーチファクトを抑制するなど、画質を向上できる画像再構成装置及びプログラムを提供することである。上記従来例の問題点を解決するための発明は、画像再構成装置であって、観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、を含むこととしたものである。また本発明の一態様に係る画像再構成装置は、前記確率的モデルを表す演算式を生成する手段は、前記観測対象に基づいて投影像が得られるまでの過程に係る条件として、撮影時の外乱による未知パラメータを含む条件を表す演算式を生成する手段であり、当該未知パラメータを推定する手段をさらに含むこととしたものである。さらに本発明の別の態様に係る画像再構成装置は、前記投影像は、観測対象を透過する電磁波または粒子線を検出することによって得たものであり前記確率的モデルを表す演算式を生成する手段は、前記観測対象の透過に伴う前記電磁波または粒子線の、波長に依存したエネルギーの減衰が、前記投影像にもたらす影響を表す変換行列を用いて、演算式である投影変換行列を生成することとしたものである。さらに本発明の別の態様に係るプログラムは、コンピュータを、観測対象に対して互いに異なる方向から得た、観測対象を透過した複数の投影像に係る画像データを受け入れる手段と、前記観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、前記投影像を生成する確率的モデルを表す演算式を生成する手段と、ベイズ推定を用いて、前記確率的モデルを表す演算式と、前記画像データとから観測対象の断面像を推定する手段と、として機能させることとしたものである。本発明の実施の形態について図面を参照しながら説明する。本発明の実施の形態に係る画像再構成装置1は、図1に例示するように、制御部11と、記憶部12と、入力部13と、出力部14とを含んで構成されている。ここで制御部11は、CPUCentralProcessingUnit等のプログラム制御デバイスであり、記憶部12に格納されたプログラムに従って動作する。本実施の形態では、この制御部11は、観測対象に対して互いに異なる方向から得た、当該観測対象を透過した複数の投影像に係る画像データを受け入れ、当該観測対象の像の画素間に、予め定めた相関があることを表す分布関数と、当該観測対象に基づいて投影像が得られるまでの過程に係る条件とに基づいて、投影像を生成する確率的モデルを表す演算式を生成し、ベイズ推定を用いて、当該確率的モデルを表す演算式と、受け入れた画像データとから観測対象の断面像を推定する処理を実行する。この制御部11の詳しい処理の内容については後に述べる。記憶部12は、メモリデバイスやディスクデバイス等であり、制御部11によって実行されるプログラムを保持している。このプログラムは、例えばDVD-ROM等のコンピュータ可読な記録媒体に格納されて提供され、この記憶部12に複写されたものであってもよい。また、この記憶部12は、制御部11のワークメモリとしても動作する。入力部13は、例えばシリアルまたはパラレル入力インタフェースを含み、CT装置にて撮像されて得られた画像データを受け入れて制御部11に出力する。出力部14は、ディスプレイやプリンタ等であり、制御部11から入力される指示に従って画像データを表示出力する。本実施の形態の画像再構成装置1が処理の対象とする複数の投影像は、観測対象に対して互いに異なる方向から、観測対象を透過し、または観測対象により散乱された電磁波や粒子線等を、複数の検出素子channelを配列した検出器detectorで撮像して得られるもので、一例としては、観測対象を輪切りにする仮想的な平面上で、観測対象に対して、複数の角度方向からX線を照射する一般的なCT装置による投影像である。つまり観測対象を上記平面で切ったときの断面像を、複数の角度から投影した複数の投影像を処理の対象としている。本実施の形態の画像再構成装置1の制御部11は、機能的には図2に例示するように、画像データ受入部21と、輝度ベクトル抽出部22と、変換行列生成部23と、ベイズ推定部24と、再構成部25とを含んで構成されている。画像データ受入部21は、入力部13から複数例えばK個の投影像の画像データの入力を受け入れ、それぞれを記憶部12に蓄積して格納する。ここで投影像の画像データはそれぞれ複数の画素値検出素子ごとの出力値を表すデータを含む。この画素値のデータは、例えば輝度の値を含む。輝度ベクトル抽出部22は、画像データ受入部21が記憶部12に格納したk番目k=1,2,...,Kの画像データの各々について、予め定めた順に画素値のデータ例えば輝度の値を配列して、輝度値ベクトルyk=1,2,...,Kを得る。このK個の輝度値ベクトルの集合D={y}がサイノグラムsinogramに相当する。変換行列生成部23は、ベイズ推定に用いる事前分布に係る演算式を生成する。すなわち、元の画像断面像を表す輝度ベクトルをxとし、投影像を得るためのラドン変換を表す投影変換行列をW、撮像時の条件により生じる誤りをεとして、と表すことができる。すなわち、ここで、元の画像である断面像と投影像との対応関係、つまり元の画像の輝度ベクトルxからサイノグラムDが得られる確率を表す尤度関数pD|xは、投影変換行列Wにより一意に定まる部分を除いた部分ε=WD-x、誤りノイズεの生じる確率pεに等しいから、を得る。ここでβは精度分散の逆数であり、IはWと同じ次元の単位行列、Mはyの次元、Ny|a,bは平均ベクトルa、分散行列をbとする分布関数である。この分布関数は具体的にはガウス分布やポアソン分布など広く知られた分布関数のうちから、予め利用者が選択しておく。また変換行列生成部23は、投影変換行列Wを、少なくとも一つのアフィン変換を表す行列平行移動、回転、縮小拡大等を表す行列またはそれらの積複数の場合として、例えば次のように定める。すなわち観測対象から投影像が得られるまでの過程に係る条件として、k番目の投影像を得たときの角度をθk、投影像のサイズ輝度ベクトルの次元をM、観測対象の元の像のサイズ輝度ベクトルの次元をN、撮像により点像が、投影像においていかなるサイズの円盤像に投影されるかを表す、円盤像の幅、つまり、点広がり関数pointspreadfunctionの大きさγを用い、当該θkの回転を表す行列Rθk、ξだけ平行移動をすることを表す行列をTξ、r倍の拡大縮小を表す行列をMrとして、まず行列を求める。この行列Fは、位置にある画素を所定の位置に移動したうえで、θkだけ回転させ、rr=M/N倍して、を当該r倍した位置まで再度移動することを意味する。これにより、次のように画素ベクトルの要素vjを変換すれば、投影像上の画素値ujが得られる。そしてこのujを用いたにより、投影変換行列Wの要素を次のように表すことができる。また変換行列生成部23は、観測対象の像に関する制約条件に基づく事前分布つまり元の画像である断面像自体の評価関数を表す演算式pxを生成する。具体的にこの分布は次のように定められる。すなわち観測対象の像再構成されるべき像においては画素間に相関があり、互いに隣接する画素同士では輝度の差が所定の閾値を超えないとする。このようなpxとしては、例えばガウス分布がある次式。なお、分散を表す固定行列fixedmatrixZxは、であり、ここにviは位置iにある画素値ここでは輝度の値、rは画素間の距離、Aは経験的に定められる係数である。変換行列生成部23は、これら生成したpD|x、及びpxに係る演算式2,4式を、確率的モデルを表す演算式として出力する。また、このpxとしては、ポアソン分布等、他の分布を用いてもよい。ベイズ推定部24は、ベイズ推定を用いて事後分布px|Dを演算する。これは、再構成されるべき元の画像である断面像自体の評価関数と、断面像と投影像との関係に基づく評価関数とのそれぞれの値を最大化するような画像を断面像の推定結果として探索することに相当する。具体的にはベイズの定理から、であり、このうち分母のpDは、すべてのxについて分子のpxpD|xの和をとったものに他ならないから、予め規格化して1としておく。すると、サイノグラムDが与えられたとき、像がxである確率を与える事後分布px|Dは結局、となる。ここにμは平均のベクトルである。ベイズ推定部24は、この7式に含まれる平均ベクトルμを演算する。この平均ベクトルμは、である。ここに、Σは分散行列であり、と表される。再構成部25は、平均のベクトルμを像の輝度のベクトルとして演算し、出力する。本実施の形態の画像再構成装置1は、以上の構成を備えてなり、次のように動作する。すなわち、既知である64×64画素の2次元画像Sに対し、その中心まわりに、予め定めた角度から、角度θkk=1,2,...,Kだけ傾いた走査方向に沿って、当該走査方向に対し垂直な線分に射影した、それぞれn画素の複数の像この場合は複数の一次元像Y1,Y2,...,YKを得る図3。本実施の形態の画像再構成装置1は、図4に示す処理を行い、これらの、それぞれn画素分の複数の一次元像を受け入れる。そして、各一次元像のそれぞれの画素の輝度値を配列した複数のn次元のベクトルy1,y2,...,ykを得るS1。なお、ここでは輝度値としているが、カラーの場合、予め定めた色空間上の値例えばR、G、Bの各値のベクトルとしてもよい。また画像再構成装置1は、観測対象から投影像が得られるまでの過程に係る条件として、k番目の投影像を得たときの角度をθk、投影像のサイズ輝度ベクトルの次元をM、観測対象の元の像のサイズ輝度ベクトルの次元をN、撮像により点像が、投影像においていかなるサイズの円盤像に投影されるかを表す、円盤像の幅、つまり、点広がり関数pointspreadfunctionの大きさγを用い、当該θkの回転を表す行列Rθk、ξだけ平行移動をすることを表す行列をTξとして、3式で表される行列Fを求めるS2。そして画像再構成装置1は、このFにより画素ベクトルの要素vjを変換して投影像上の画素値ujを得て、このujを用いて要素を定めた投影変換行列Wを得るS3。そしてこの行列を用いてベイズ推定に用いる尤度関数、pD|xを2式のように求めるS4。また画像再構成装置1は、観測対象の像に関する制約条件に基づく事前分布を表す演算式pxを生成し4式S5、これら生成したpD|x、及びpxに係る演算式2,4式を確率的モデルを表す演算式とし、ベイズ推定の式6式を用いて、事後推定の結果である平均ベクトルμを8式により演算するS6。画像再構成装置1は、この平均ベクトルμの表す像を再構成した画像として出力するS7。ここで平均ベクトルμは、再構成した画像に含まれるべき画素の輝度値画素値を予め定めた順に配列してベクトルとしたものであるから、当該順に従って元の画像の各画素の値を、当該平均ベクトルμにおける対応する成分の値から定める例えばその値そのものとすることで、元の画像が得られる。このような本実施の形態によると、元となる観測対象の断面像から投影像が得られるプロセスをそのまま観測画像の生成の確率的モデルとして記述することで、ベイズ統計学に基づく演算により、投影像から最尤断面像観測対象の画像を直接的に推定する。また、これより、投影変換行列を観測対象から投影像が得られるまでの過程に係る条件に基づいて設定することで、例えばヘリカルスキャン、コーンビーム、マルチCT、位相CT、3次元CT等の種々の投影過程による投影像を対象として再構成像を得ることができる。さらに本実施の形態の画像再構成装置1は、ラドン変換を表す投影変換行列Wを演算するにあたり、観測対象から投影像が得られるまでの過程に係る条件を表す演算式として、3式に定めた行列Fに代えて、次のような行列Fを用いることとしてもよい。すなわち、k番目の投影像を得たときの角度をθk、投影像のサイズ輝度ベクトルの次元をM、観測対象の元の像のサイズ輝度ベクトルの次元をN、撮像により点像が、投影像においていかなるサイズの円盤像に投影されるかを表す、円盤像の幅、つまり、点広がり関数pointspreadfunctionの大きさγを用い、当該θkの回転を表す行列Rθk、ξだけ平行移動をすることを表す行列をTξ、さらにk番目の投影像を得たときの位置ズレを表すシフトベクトルSkによる平行移動TSkを含め、としてもよい。この行列Fは、k番目の投影像を撮像したときにSkだけの位置ズレが生じたとしたものであるが、実際には、この位置ズレの量Skは、外乱であるために測定が困難である。また、θk、γについても予め定めたものから、外乱のために実際の撮影時にはずれている可能性もある。すなわち、Sk,θk,γは、未知のパラメータであるということができる。そこで本実施の形態の画像再構成装置1は、これら未知のパラメータを、いわゆるEMアルゴリズム期待値最大化法や、直接推定の方法を用いて推定する。EMアルゴリズムを用いる場合、未知パラメータのセットをαとおく。これにより、最終的に求めたい確率的モデルを表す演算式の一つの対数値を考え、この期待値を最大化する。すなわち画像再構成装置1は、当初は、この未知パラメータのセットの初期値αt=α0を仮に定める。当初はSk,γを0、θkはk番目の投影像を得たときのX線の飛跡方向を表す角度、つまり撮影時の設定角度としておいてもよいし、ランダムに決定してもよい。そして画像再構成装置1は、未知パラメータのセットがαtであるとき、期待値を評価するステップEステップでは、11式のうち、αtに係る項であるを演算する。以下この12式を最大化することを考える。そしてこの12式の期待値を最大化するステップMステップとして、12式を最大化するαを求め、これを次の未知パラメータのセットαt1とする。すなわち、とする。画像再構成装置1は、このαが収束するまで、つまり、ベクトルαt1-αtの大きさが予め定めたしきい値を下回るまで、EステップとMステップとを交互に繰り返し実行する。そして収束後のαが得られれば、当該αに含まれる未知パラメータの値を用いて、行列Fを求め、さらには投影変換行列Wを求めて、尤度関数pD|xを定める。また、この未知パラメータを推定する方法として、画像再構成装置は、周辺尤度を最大化するαを直接的に求めてもよい。すなわち周辺尤度をとして、これを最大化するαを求めることとしてもよい。またCT撮影では、金属が含まれる生体組織の投影像を得た場合など、当該金属における高いX線吸収率により、投影像にデータが欠損した領域が生じる場合がある。これは例えば金属製のクラウン歯科で用いる被せ物が含まれる口腔を撮像した場合等にあり得るが、従来の方法では当該欠損した領域により生じるアーチファクトいわゆるメタルアーチファクトが避け得なかった。本実施の形態の画像再構成装置1では、観測対象の像に関する制約条件に基づく事前分布pxとして、当該欠損した領域を含んだ推定像のエッジを表現するイジングモデル等を用いることなどとすることで、このメタルアーチファクトを軽減できる。そのほか、この投影変換行列Wを求める際の行列Fに含めることのできる未知パラメータとして、投影像の撮像に用いた撮像素子等の特性感度などがあり、観測対象の像に関する制約条件に基づく事前分布pxを定める際に、推定される断面に含まれるべき物体のX線吸収率などの物質の特性を表す情報等、種々のものを利用できる。これにより、予め分かっている観測対象の像断面像に関する情報を有効に利用して断面像の再構成の処理を行うことができる。また、投影変換行列Wを求める際の行列Fに得られる画像の解像度に係るパラメータを含むので、投影像よりも解像度の高い画像を再構成する、いわゆる超解像処理を行うことも可能となっている。さらに、撮像角度が制限されている場合や、撮像対象の内部に、検出素子で検出する電磁波や粒子線源がある場合蛍光体がある場合、これら検出素子で検出する電磁波や粒子線を散乱ないし吸収する物質が撮像対象の内部に含まれている場合には、観測対象から投影像が得られるまでの過程に係る条件を表現した行列Fを生成するときに、それぞれの条件を表す行列を乗じたものとして生成すればよい。また撮像時にゾーンプレートを用いている場合など、非線形な散乱が発生する場合は、当該非線形な散乱を線形に近似した行列Lを作成する。そしてこの行列Lを含む行列の積により、3または10式に示したものと同様に、行列Fを得る。また隣接する検出素子間で出力値に相関がある場合も、その旨を表す行列を含む複数の行列の積により、行列Fを得る。さらに、本実施の形態では、投影像を得るまでの過程がわかっていれば、いままでに得られているCTデータから再構成した画像を得ることができる。さらに電子線やX線を利用したCT等ではビームハードニングBeamHardeningアーチファクトの影響がある。これら検出素子で検出する電磁波や粒子線は一般に単波長ではなく、比較的波長の長い成分から短い成分まで種々の波長の電磁波や物質波ドブロイ波を含むものである。このビームハードニングは、これらX線等が撮像対象を透過する間に、比較的波長の長い成分がより吸収されて、透過に従って短い波長の成分が多くなり、結果としてエネルギーが高くなる現象をいう。医療用にCTを用いる場合にこの現象が発生すると、体内に実在しない低吸収領域があるように見えてしまう場合がある。そこで3または10式に示したような、行列Fを行列の積のうちに、ビームハードニングX線等のエネルギーの減衰に伴うCTデータの変化を表す行列を含める。具体的には、行列Fを、k番目の投影像の撮像過程で生じる、投影対象の各部分における、波長に依存した電磁波または粒子線のエネルギーの減衰量を未知パラメータとして含む行列Aaλ1,aλ2...と、その他、撮像過程での変換を表す行列Mr、TSk等との積により表すこととすればよい。ここで、aλnは、予め定めた複数の波長例えば、吸収率が互いに異なると考えられる各周波数帯から少なくとも一つずつの代表波長を選択することにより定めればよいλ1,λ2,...λnの成分ごとの吸収係数波長λii=1,2,...の成分の減衰割合を意味する。つまり、この行列を含んだ複数の行列の積で表される行列Fは、検出素子ごとの出力値である画素値を、多次元の吸収係数で表現するものとなる。これにより、波長帯ごとの吸収率を推定し、波長帯ごとの吸収率の相違によって発生するアーチファクトなどを軽減する。このように、観測対象から投影像が得られるまでの過程に係る種々の条件を表した行列Fを生成することで、種々の撮像態様や、撮像対象の性質を反映した画像の再構成を行うことが可能となる。本発明の実施例として、画像再構成装置1を、制御部11であるCPUとしてIntelXeno266GHz4コアを採用したコンピュータを用いて実現した例を次に示す。図5は、元の64×64画素の二次元像Aについて得た32×8、32×16、32×32、32×64画素のサイノグラムRを得て、それぞれ単純に逆ラドン変換を行って得た再構成画像B、フィルタを用いた再構成画像C、逐次近似法による再構成画像D、本実施例の画像再構成装置1により得た再構成画像Eを対比したものである。図5において、各再構成画像B,C,D,Eは上から順に、それぞれ32×8、32×16、32×32、32×64画素のサイノグラムから得られる再構成画像を表す。この図5の例において、再構成画像C,D,Eを表すベクトルxにハット^を付していると、原画像xと、単純に逆ラドン変換を行って得た再構成画像を表すベクトルxにチルド~を付しているとを用いて演算される再現性の評価値ISNR式16を得た結果を、図6に示す。図6において、横軸はサイノグラムのサイズ、縦軸はISNR値を示す。図6に例示するように、本実施例の方法によると、フィルタを用いる例や逐次近似を用いる例に比べ、再現性を大幅に改善できる。"