JP3730927B2 - 吸収物質濃度空間分布画像化装置 - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明は、散乱物質中に吸収物質が存在してなる被検体において、吸収物質濃度の空間分布を求める方法に関する。特に、光を用いて生体内部あるいは大気中の吸収物質濃度の空間分布を計測する方法、装置に関する。
【０００２】
【従来の技術】
現在、医療用の画像診断装置として、ＭＲＩ（核磁気共鳴イメージング装置）Ｘ線ＣＴ、超音波診断装置等がある。それぞれ、臨床医療の現場で多用され有用性が認められている。現在、臨床医療に使用されている画像診断装置は、主に生体内の組織、構造を断層像として得るものが多く、生体内に分布する化学物質の代謝活動に起因する生体機能を計測するための簡便で実用的な画像診断装置は少ない。生体機能の計測により得られる情報は、脳疾患やガンの予防診断、脳機能の解明に有用と期待されており、簡便、実用的な生体機能画像診断装置の開発が医療の分野では待たれている。散乱物質中に、光の波長により吸収特性が異なる吸収物質が存在してなる被検体中の吸収物質の濃度の空間分布を画像化する光計測装置として、生体内の特定の色素濃度の空間分布を画像化する生体光計測装置がある。生体中には、生体機能に関与する吸収物質が多く存在し、吸収物質の濃度の空間分布を計測できる生体光計測装置は、生体機能を画像化して診断に有効となる装置の一つとして開発が急がれている。
【０００３】
生体機能の中でも、生体内酸素飽和度は、生体中に多く含まれる特定色素の濃度に対応しており、この色素濃度は可視から近赤外領域における光吸収の量から計測することが可能である。特定色素として、酸化型あるいは還元型の別により、異なる光吸収の分光特性を持つことが知られている、血液中に存在するヘモグロビン（以下では、Ｈｂと略記する）、細胞中に存在するチトクロームａａ３、筋肉中に存在するミオグロビン等がある。さらに、光による計測は、光ファイバーや小型光源を用いることにより、簡便、小型な装置として開発される可能性を持っている。
【０００４】
可視から近赤外の光を用いた生体機能を計測する装置が、例えば、特開昭５７−１１５２３２号公報、特開昭６３−２７５３２３号公報に記載されている。生体光計測装置で重要な課題の１つとして、生体内に存在する吸収物質の空間分布を画像化する方法の確立が挙げられる。この課題は吸収物質の濃度と位置を求める問題である。従来技術での画像化方法として、（１）Ｘ線ＣＴのアルゴリズムを用いる方法、（２）輸送方程式を用いる方法、（３）モンテカルロシミュレーションの結果を用いる方法、に大別される。
【０００５】
（１）の方法の具体的例として、時間ゲート法を用いて得た計測データから、Ｘ線ＣＴにおいて使用されてきた画像再構成のアルゴリズムにより画像を得る方法（特開平０１−２０９３４２号公報に記載）が挙げられる。
（２）の方法の具体的例として、アリッジ( Aridge )等の方法（" Reconstruction methods for infrared absorption imaging ", Proc. SPIE 1431, 204-215(1991)に記載）や、シンガー( Singer )等の方法（" Image reconstruction of the interior of bodies that diffuse radiation ", Science, 248, 990-993(1990)に記載）が挙げられる。両者の方法共に、散乱と吸収現象を表わす数学モデルに基づき、照射光強度と散乱を表わす変数（以下では、散乱変数と略記する）と吸収物質の吸収を表わす変数（以下では、吸収変数と略記する）とから、輸送方程式を立て検出光強度を計算する。ここで、散乱変数と吸収変数は散乱体内部の位置の関数であり、また、吸収変数は吸収物質の濃度を反映する値である。上記の輸送方程式の散乱変数と吸収変数の値を散乱体内部の任意の位置で独立に変化させ、計算された出光強度と実際に計測された検出光強度との差を評価し、その差が最小となるまで計算を繰り返す。この繰り返しの結果得られた、吸収変数の空間分布が吸収物質濃度の空間分布を表わしており、その分布を画像化して表示する。アリッジ( Aridge )等の方法では、散乱変数として散乱係数を、吸収変数として吸収係数を用いており、シンガー( Singer )等の方法では、散乱変数として散乱体内に設定された面積要素もしくは体積要素から出て行く光子の確率を、吸収変数として要素内で吸収される光子の確率を用いている。
【０００６】
（３）の方法の具体的例として、バーバー ランダル エル等によるモンテカルロ法を利用した画像化方法がある（特表平３−５０５９２２号公報に記載）。この方法の概要は以下の通りである。媒体中の全ての点（位置及び方向）に対し、フラックスによって与えられる重みを割当て、その点で検出器の吸収体の存在しない媒体に対する応答に対する光子の期待される寄与を割り当てることによって、複雑な媒体を画像化する技術に関する。照射される輻射線の測定された強度が、予定される強度と比較されることにより、対応する散乱体積の透過係数が決定される。予定される強度は、予め規定されたモデル媒体に対するモンテカルロシュミレーションを用いた計算により得られ、透過強度係数が重み関数に適用される。媒体中における微小な吸収の検出器応答特性に対する影響は、影響を受ける体積要素に対する吸収散乱断面積と対応する重み関数との積となるものと仮定される。重み関数は、モンテカルロ法を用いて計算され、領域内の全ての点の、与えられた光源−検出器の配置に対する検出器応答性への相対寄与の３次元マップを与える。画像化は、透過係数と、光源−検出器の各配置ごとに、各体積要素に対応する重み関数との積を計算することにより達成される。これらの積の総和は、光源−検出器の配置のすべてに対するすべての体積要素の重ね合わせに対応する。この計算は、全ての光源の位置の関数として、全ての領域にわたる相対吸収断面積の３次元マップを与える。
【０００７】
【発明が解決しようとする課題】
生体光計測装置は、Ｈｂ（ヘモグロビン）のように、生体内で機能する吸収物質濃度の空間分布を、画像として表示して、例えば癌の様な疾病の予防診断への応用が期待されている。しかし、生体は光に対して強い散乱特性を有しているため（散乱係数は１．０〔ｍｍ^{- １}〕程度である）、生体内部の吸収物質濃度の空間分布の精度良い画像化は非常に困難である。生体光計測装置を臨床医療で実用化するためには、生体の様な散乱体の内部に存在する吸収物質濃度の空間分布を精度良く画像化する方法が重要である。現在、このような生体内部の吸収物質濃度の画像化方法として、実用化レベルに達していると考えられる方法はまだ存在しない。従来技術における方法では、様々な問題があり、上記（１）の方法では、３次元的な散乱現象を２次元空間上で扱うため再構成された画像の空間分解能が低い。上記（２）の方法では、繰り返し計算をするため長時間の演算が必要で、現状の計算機速度では高分解能の画像を得ることは実質的に困難であり、また局所最小値（ local minimun ）に落ちこみ真の最適解が得られにくい場合がある。上記（３）の方法では、不透明な媒体中にある標的物体を３次元的に画像化することの記載はあるが、画像化する具体的な手順は数式的に明確には示されておらず、さらに媒体中の吸収物質の濃度を計算することの記載はない。
【０００８】
本願発明の目的は、３次元的な散乱現象を考慮して、特に生体内で機能する吸収物質濃度を精度良く、短時間に空間分解能で求め、吸収物質の濃度の３次元空間分布を求める方法、及び装置を提供することにある。
【０００９】
【課題を解決するための手段】
本願発明は、被検体を通過して検出された光強度が、被検体中に存在する吸収物質の濃度と、散乱に基づく光子の平均光路長との線形結合により表わされると仮定し、平均光路長の空間分布を使用して、被検体中に存在する吸収物質の濃度の３次元分布を求めることに特徴を有する。即ち、光吸収物質を含んでなる光散乱体からなる被検体に、所定の光照射位置から所定の波長を有するパルス光又は連続光を照射する照射ステップと、被検体を通過する光の強度を所定の光検出位置で検出する検出ステップと、を有し被検体内に分布する吸収物質濃度の空間分布画像化方法において、光照射位置から光検出位置までの複数の光路を求めるステップと、複数の光路から平均光路長を求めるステップと、平均光路長と、被検体に照射される光の強度と、光検出位置で検出される検出光強度と、被検体に照射される光に対する吸収物質の光学定数とから、各体積要素毎の吸収物質の濃度を求めることを特徴とする。平均光路長の空間分布は、モンテカルロ法により模擬して計算により求める。この計算では、被検体と実質的に同一の散乱係数を有し、仮想被検体を包含する大きさの仮想散乱体内の一点から発する複数光子のそれぞれのたどる複数光路を求める。仮想被検体は被検体と同一又は類似した形状を有し吸収物質を含まないものとする。
【００１０】
より詳細に説明すると、本願発明の吸収物質濃度の空間分布画像化方法では、被検体と同一又は類似した形状を有し、吸収物質を含まない仮想被検体を想定して、仮想被検体において設定される光照射位置と光検出位置の複数の対のそれぞれで、光子がたどる光照射位置から光検出位置に至るまでの複数の光路を求めるステップと、複数の光路から光照射位置から光検出位置までの仮想被検体内における平均光路長を求めるステップと、平均光路長と、被検体に照射される照射光強度と、光検出位置で検出される検出光強度と、被検体に照射される光に対する吸収物質の光学定数（吸光係数）と、から吸収物質の濃度の空間分布を求めるステップと、吸収物質の濃度の空間分布を表示するステップと、を有することに特徴を有する。被検体の形状を表す空間は複数の体積要素に分割され、各体積要素ごとにに平均光路長が求められ、各体積要素内の吸収物質濃度が一様であると仮定し、各体積要素毎の吸収物質濃度が求められる。
【００１１】
より具体的な吸収物質濃度の空間分布画像化のステップは、上記対における照射光強度と検出光強度との比に基づいて計測値行列を生成するステップと、上記対における平均光路長に基づいて平均光路長の空間分布行列を生成するステップと、計測値行列と平均光路長の空間分布行列から吸収物質濃度を求めるステップ、とからなる。予め与えられている所定の形状と所定の散乱係数を有する被検体に本発明の手順を適用する場合は、平均光路長の空間分布行列を生成するステップで求められ記憶手段に予め記憶されている平均光路長の空間分布行列を使用する。
【００１２】
光路を求めるステップは、被検体と実質的に同一の散乱係数を有し、仮想被検体を包含する大きさの仮想散乱体内の一点から発する光子のたどる光路を模擬するモンテカルロ法により模擬して求められ、モンテカルロ法では、光子が順次散乱される複数の散乱点の内の連続する２つの散乱点のそれぞれの空間座標位置を関連づける散乱角度特性を表す関数、及び乱数を使用する。モンテカルロ法により求められた複数の光路のそれぞれを表す複数の散乱点の空間座標位置が記憶され、複数の光路データが記憶手段に記憶される。光路データは、仮想散乱体内の一点から発する単一光子のたどる有限光路長を有する光路データを求めるか、あるいは仮想被検体の２点を結ぶ最大距離以上の半径を有する球を仮想散乱体として、球の中心から複数光子を発して、複数光子のそれぞれの光子のたどる光路のうち、球の内部に存在する光路データを求める。
【００１３】
仮想被検体での光照射位置から光検出位置までの光路は、仮想被検体での光照射位置と仮想散乱体内の１点とを一致させ、光検出位置とモンテカルロ法により求められた光路とが交差する点を、光路を表す散乱点の空間座標、又は仮想被検体の形状を表す空間座標のいずれか一方の空間座標を、他方の空間座標に対して座標変換（回転移動及び／又は並進移動）して求める。求められた仮想被検体での光照射位置から仮想被検体での光検出位置までの複数光路のそれぞれのを、被検体での対応する光照射位置から光検出位置までの光路とする。なお、仮想被検体の形状を表す空間は、被検体の形状を表す空間を複数の体積要素に分割したと同様に体積要素に分割される。また、光照射位置と光検出位置の対の数は、被検体、仮想被検体を分割した体積要素の総数以上の値に設定される。上記の複数光路データから仮想被検体の各体積要素を通過する平均光路長を体積要素ごとに求め、被検体の各体積要素を通過する平均光路長とする。光照射位置（ｉ）と光検出位置（ｊ）の対（ｉ−ｊ）に対して得られた光路に関するデータに基づいて、光照射位置と光検出位置の位置関係が入れ替わった光照射位置（ｊ）と光検出位置（ｉ）の対（ｊ−ｉ）に対する光路に関するデータを、対（ｉ−ｊ）に対して得られた光路に関するデータから、座標変換により求めることができる。
【００１４】
照射ステップにおいて被検体に照射される光は、４００ｎｍから２０００ｎｍの範囲にある複数の波長を有する光であり、複数の波長の差が１００ｎｍ以下とする。被検体を通過する光の強度が時間分解計測される場合には、検出光強度は時間分解計測された光強度の所定の時間内での光強度の積分値から求められ、所定の時間内に光照射位置から光検出位置に到達する光子に関する光路から求められる平均光路長に基づいて平均光路長の空間分布行列が生成される。
【００１５】
本発明では、計測された検出光強度から計測値行列Ｙを生成し、散乱体内部の吸収物質の吸光係数と、吸収物質の存在しない散乱体を模擬して得る散乱体内における光が通過する光路データをモンテカルロ法により求め、光路データから、平均光路長の空間分布行列Ａを生成し、行列Ｙと行列Ａとから、未知数行列Ｘを導出して、求められた未知数行列Ｘから被検体内部の吸収物質濃度を求めこの空間分布を表示する。さらに、異なる複数位置に配置される光照射位置と光検出位置の対の数は、被検体、仮想被検体を分割した体積要素の総数以上の値に設定するので、ある体積要素に注目するとき、複数方向からの散乱の影響が考慮されるので、従来方法に比較して、約３倍以上の空間分解能を有する３次元的な吸収物質の空間分布を得ることができ、１０〜１００倍程度計算時間が短縮され、さらに定量性の改善が可能である。本方法では、局所の最適解に落ち込むことなく解を得ることができる。さらに、モンテカルロ法による光路の計算が、予め計算された光路データが記憶保存されている場合には、記憶保存されている光路データを使用するので、従来のモンテカルロ法を適用する方法と比較して、平均光路長の空間分布を約１００倍以上の速度で求めることができる。
【００１６】
【発明の実施の形態】
従来技術として、光を用いて生体内部に分布する特定色素濃度の空間分布を計測する生体光計測装置があり、また大気中に存在する吸収物質濃度の空間分布の計測として、レーザーレーダがある。この様に、散乱物質中に存在する吸収物質の濃度の空間分布を求め画像化する方法を必要とする分野は広い。特に以下では、生体光計測装置を代表的な例として説明をする。
【００１７】
本願発明は、光計測装置で計測された検出光強度から、被検体内部に存在する吸収物質の濃度を画像化する方法に係り、本願発明における手順は要約すると、計測された検出光強度から計測値行列Ｙを生成するｓｔｅｐ１、散乱体内部の吸収物質の吸光係数と、吸収物質の存在しない散乱体を模擬して得る散乱体内における光が通過する光路データをモンテカルロ法により求め、光路データから、平均光路長の空間分布行列Ａを生成するｓｔｅｐ２、行列Ｙと行列Ａとから未知数行列Ｘを導出するｓｔｅｐ３、求められた未知数行列Ｘから被検体内部の吸収物質濃度を求めこの空間分布を表示するｓｔｅｐ４からなる。モンテカルロ法による光路の計算が、ｓｔｅｐ２の処理開始以前に既になされ、光路データが記憶保存されている場合には、記憶保存されている光路データを使用する。
【００１８】
まず最初に、被検体によって照射光が散乱されない場合（即ち、被検体が非散乱体である場合）に関して説明し、次いで、被検体によって照射光が散乱される場合（即ち、被検体が散乱体である場合）に関して説明する。図１に示す様な、光源１から被検体２へ所定の波長の光を照射し、光路３を通過してきた光の強度を検出器４で検出する光計測系を考える。図１は、光路３を含む被検体２の断面を示す。被検体２中に１成分の吸収物質が一様濃度ｃで存在する場合に、光源１からの照射光強度をＩ₀、検出器４で検出される検出光強度をＩ_d、吸収物質の吸光係数をε、光照射位置と光検出位置間の幾何学的距離で決まる光路長をＬ’とすると、（数１）で定義されるベール・ランバート則が成り立つ。
Ｉ_d／Ｉ₀＝ｅｘｐ（−εｃＬ’） …（数１）
図１と同様の光計測系を考え、図２に示すように、被検体２が、吸収物質濃度の異なる２個の体積要素５−１、５−２に分けられる場合に、体積要素５−１中での吸収物質濃度をｃ₁、光路長をＬ₁’、体積要素５−２中での吸収物質濃度をｃ₂、光路長をＬ₂’とし、光源１からの照射光強度をＩ₀、体積要素５−１から５−２に入射する光の強度をＩ_C、検出器４で検出される検出光強度Ｉ_dとすると、（数１）より、体積要素５−１、体積要素５−２に関し、それぞれ（数２）、（数３）が成り立つ。
Ｉ_C／Ｉ₀＝ｅｘｐ（−εｃ₁Ｌ₁’） …（数２）
Ｉ_d／Ｉ_C＝ｅｘｐ（−εｃ₂Ｌ₂’） …（数３）
体積要素５−１から５−２に入射する光の強度Ｉ_Cは計測されないので、（数２）、（数３）を用いてＩ_Cを消去すると、（数４）が成り立つ。
Ｉ_d／Ｉ₀＝ｅｘｐ（−ε（ｃ₁Ｌ₁’＋ｃ₂Ｌ₂’）） …（数４）
次に、被検体が散乱体である場合について記述する（以下の説明では、特に明記しない場合は、被検体は散乱体とする）。
【００１９】
図３に示すような光計測系で、被検体２中に１成分の吸収物質が一様な濃度ｃで存在する場合に、光源１からの照射光強度をＩ₀、検出器４で検出される検出光強度をＩ_d、吸収物質の吸光係数をε、被検体内部に吸収物質が存在しないときの光路長（以下光路長と略す）をＬ、散乱によって被検体２の外に出て行く成分を減衰定数Ｄ_Sとすると、（数５）が成り立つ（ Victor Twersky,"Absorption and Multiple Scattering by Biological Suspensions ", J. Opt. Soc. Am., 60, 1089(1970)に記載）。
Ｉ_d／Ｉ₀＝Ｄ_S・ｅｘｐ（−εｃＬ） …（数５）
ここで、光を光子として扱うと、光源１から照射された光子は被検体２中で散乱を繰り返して検出器４に到達する。光子は散乱の度に進行方向が変わるので、その光路は、図３中に示す光路３の様にジグザグを描く。なお、図３では、図１、図２と同様に被検体２は立体であり、簡単のために、２本の光路３が１つの断面内にある場合を例示しており、一般に、光路３は３次元的にジグザグを描く。散乱による進行方向の変化はランダムであり光路３は各光子によって異なり、（数５）中の光路長Ｌは複数個存在し一義的に決定できない。そこで、複数個存在する光路長の平均値（以下では、平均光路長と略記する）を近似値として用いる。検出器に到達する光子数をｍ本とすると、光路の数はｍ本存在し、ｍ本ある光路のうちｋ番目の光路長をＬ_Kとして、平均光路長Ｌ_Aを（数６）で定義し、さらに、平均光路長Ｌ_Aを用いて（数５）を書き換えると（数７）になる。
Ｌ_A＝（ΣＬ_K）／ｍ …（数６）
（数６）において、加算Σはｋ＝１〜ｍについて行なう。
Ｉ_d／Ｉ₀＝Ｄ_S・ｅｘｐ（−εｃＬ_A） …（数７）
図３に示す光計測系を、図４に示すように、被検体２を、ｎ個の体積要素に分割して考え、光源から発した光子は、全ての体積要素を通過して検出器４に至るものと仮定する。さらに各体積要素における吸収物質濃度は一般に異なると考える。図４は、簡単のために、図３で例示したのと同様に、光源１から照射された光子が被検体２中で散乱を繰り返して検出器４に到達する光路が全て、厚さが１体積要素からなる１断面内にある場合を示し、各体積要素の中心部分に示される太い線分は、この１断面内に存在する複数個の光路を各体積要素に分割して考え、各体積要素において複数個の光路のそれぞれ進行方向での距離を加算し一定値で割算して得た光路長を示し（得られる光路長の最大長を体積要素の辺長さ以下とするために一定値で割算した）、光源１と検出器４を結ぶ方向に平行に示したものである。言いかえると、図４は、被検体２として２次元物体を想定し、被検体２を空間（面積）要素に分割して考え、光源１からの複数光子のそれぞれが、２次元平面内で散乱の度に進行方向を変えて光路を形成し進む場合に、各面積要素において複数光子の描くそれぞれの光路を、それぞれの光子の進行方向での距離を加算し一定値で割算して得た光路長として太い線分により、光源１と検出器４を結ぶ方向に平行に示したものである。吸収物質が存在せず、かつ一様な散乱体では、図４に示すように、この太い線分の長さの分布は光源１と検出器４を結ぶ直線に対して対称となり、この直線から距離が離れる程太い線分の長さは短くなる。しかし、一般に被検体２は一様な散乱体ではなくしかも、被検体２の各部位は異なる濃度の吸収物質が存在するため、このような場合には図４に示す太い線分の長さは、各要素で異なることになる。
【００２０】
本願発明では、光源１からの複数光子のそれぞれは、被検体２を構成する全ての体積要素を一般に通過していくと想定し、各体積要素が図３に示す被検体であると見做して（数７）を適用し、さらに複数の異なる濃度の吸収物質が存在する場合に適用される（数４）を適用する。即ち、光源１からの複数光子のそれぞれは、散乱のためジグザグを描いて進行するが、このとき各光子は、全ての体積要素を通過して検出器４に至るものと仮定すると、（数４）、（数７）を考慮して（数８）が成立するものとする。但し、（数８）において、各体積要素中での、吸収物質濃度をｃ_j、平均光路長をＬ_Ajとし、添字ｊは、ｊ番目の体積要素内における吸収物質濃度、平均光路長の値であることを示し、加算Σはｊ＝１〜ｎについて行なう。ｎは、体積要素の総数であり、図４に示す例では、５×５＝２５である。

なお、光源からの光子は、被検体を構成する全ての体積要素を一般に通過していくと想定するが、光子が通過しない体積要素では、Ｌ_Aj＝０であり、体積要素毎の平均光路長Ｌ_Aj（ｊ＝１〜ｎ）を求める方法に関しては後述するので、Ｌ_Aj（ｊ＝１〜ｎ）は既知として扱う。（数８）において、減衰定数Ｄ_Sを直接求めることは困難であるため、照射光として波長の異なる２種類の光を用い、各波長毎に計測を行ないＤ_Sを消去する。波長λ₁の照射光、波長λ₂の照射光を用いた計測に関して、それぞれ（数９）、（数１０）が成り立つ。ここで、平均光路長Ｌ_Ajは波長に依存しないものとする。Ｄ_S、εは波長に依存する。

（数９）、（数１０）において、加算Σはｊ＝１〜ｎについて行なう。
【００２１】
（数９）及び（数１０）を（数１１）及び（数１２）を用いて変数変換をすると、（数１３）及び（数１４）が得られる。

例えば、被検体が生体組織の場合には、波長λ₁、λ₂の波長差の絶対値（│λ₁−λ₂│）が、少なくとも１００ｎｍ、特に３０ｎｍ以内であれば、近似的に次の（数１５）が成り立つ。

と置いて、（数１６）の両辺の自然対数（Ｌｏｇ_e）とり（数１８）を得て、

さらに、

と置いて、
ｙ＝Σ（χ（ｊ，λ₁，λ₂）Ｌ_Aj） …（数２０）
を得る。（数２０）において、加算Σはｊ＝１〜ｎについて行なう。
【００２２】
（数２０）は、線形方程式であり、少なくともｎ組の同様の方程式を用意すればｎ個ある未知数χ（ｊ，λ₁，λ₂）（ｊ＝１〜ｎ）を求めることができる。１つの方程式は被検体に関する１つの計測位置（光源１と検出器４の位置関係）毎に成立するので、図５に示す様に、１つの光源１に対して検出器４を被検体２に複数配置して複数位置で検出するか、あるいは複数の照射位置から光を照射し、ｎ個の検出位置で光検出を行うと、ｎ組の（数２０）の方程式が得られ、未知数χ（ｊ，λ₁，λ₂）（ｊ＝１〜ｎ）と同数の連立方程式（数２１）を得ることができる。但し、ｉはｉ番目の光検出位置を示す添え字である。なお図５は、被検体２が２次元体である例を示す。

（数２２）、（数２３）、（数２４）により、行列Ｙ、Ａ、Ｘを定義すると、（数２１）で表される連立方程式は（数２５）となる。以下の説明においては、Ｙを計測値行列、Ａを光路長の空間分布行列、Ｘを未知数行列と呼ぶ。
Ｙ＝│ｙ₁…ｙ_i…ｙ_n│_T …（数２２）
（数２２）において、Ｔは転置を表し、Ｙはｎ行１列の行列である。
Ａ＝│Ｌ_Aij│ …（数２３）
（数２３）において、Ａはｉ行ｊ列の要素の値をＬ_Aijとするｎ行ｎ列の行列である。

（数２４）において、Ｔは転置を表し、Ｘはｎ行１列の行列である。
Ｙ＝ＡＸ …（数２５）
（数２５）から未知数行列Ｘを求めるためには、（数２６）で示すように、Ａの逆行列を両辺にかければよい。
Ａ^-1Ｙ＝Ａ^-1ＡＸ＝Ｘ …（数２６）
ここで、求められた未知数行列Ｘの各成分は、（数２７）で表される。

（数２７）に照射光として用いた各波長λ₁、λ₂に対する吸光係数を代入して、各体積要素中の吸収物質濃度ｃ_jを得る。
【００２３】
以上の説明では、被検体２をｎ個の体積要素に分割して、ｎ個の検出位置（光源１（光照射位置）と検出器４の対）で光検出を行ない、ｎ組の（数２０）の方程式が得たが、ｎをｎ’≧ｎなるｎ’に置き換え、即ちｎ’個の検出位置（光源１（光照射位置）と検出器４の対）で光検出を行ない、ｎ’組の（数２０）の方程式が得ても、以下のようにして未知数行列Ｘを得ることができる。以上の説明において（数２２）、（数２３）を以下のように置き換えればよい。即ち、（数２２）は、

とする。（数２８）においてＴは転置を表し、Ｙはｎ’行１列の行列である。
【００２４】
（数２３）は、

とする。（数２９）において、Ａはｉ行ｊ列の要素の値をＬ_Aijとするｎ’行ｎ列の行列である。このような置き換えを行なった（数２５）の解は一般に、
Ｘ＝（（Ａ_T）ＷＡ）^-1（（Ａ_T）Ｗ）Ｙ …（数３０）
として得られる。なお（数３０）において、Ｔは反転を表し、Ａ_TはＡの反転行列であり、Ｗは重み行列を表すｎ’行ｎ’列の行列である。重み行列Ｗをｎ’行ｎ’列の単位行列とするとき（数３０）は、
Ｘ＝（（Ａ_T）Ａ）^-1（Ａ_T）Ｙ …（数３１）
となる。
【００２５】
このように、被検体をｎ個の体積要素に分割して、ｎ’≧ｎなるｎ’個の各種方向での検出位置（光源１（光照射位置）と検出器４の対）で光検出を行なうことにより、精度よく未知行列Ｘ、即ち被検体の各体積要素中の吸収物質濃度ｃ_jを得ることができる。
【００２６】
さらに以上の説明において、計測対象である被検体が多成分の吸収物質を含む場合には、照射光に用いる波長の種類を増やして、多成分の吸収物質の濃度を求めることができる。例えば、被検体がαとβの２種類の吸収物質を含む場合、（数５）を拡張して（数３２）が成り立つ。但し、吸収物質α、βの分子吸光係数と濃度をそれぞれ、ε_α、ε_β、ｃ_α、ｃ_βで表す。
Ｉ_d／Ｉ₀＝Ｄ_S・ｅｘｐ（−（ε_αｃ_α＋ε_βｃ_β）Ｌ_A） …（数３２）
図４に示すように、被検体２をｎ個の体積要素に分割する場合には、吸収物質が１種類ある場合と同様にして、（数３２）から（数３３）が得られる。

（数３３）において、加算Σはｊ＝１…ｎについて行なう。
【００２７】
減衰定数Ｄ_Sを２種類の波長λ₁、λ₂の照射光を用いて消去すると（数３４）が得られ、（数３５）の置換を行なうと、（数３４）は（数３６）で表される。

同様にして２種類の波長λ₁、λ₃の照射光を用いて、（数３６）と同様の式を得る。
【００２８】
（数３４）、（数３６）において、加算Σはｊ＝１…ｎについて行なう。（数３６）の両辺の自然対数（Ｌｏｇ_e）をとると線形な式となるので、以下、先に説明した１種類の吸収物質がある場合と同様に、ｎ’≧ｎなるｎ’個の検出位置（光源１（光照射位置）と検出器４の対）で光検出を行ない、（数３６）を解くことができ、各体積要素内における変数χ（ｊ，λ₁，λ₂）が求まる。ここで、波長λ₁と波長λ₃から得られる各体積要素内における変数χ（ｊ，λ₁，λ₃）と前記の変数χ（ｊ，λ₁，λ₂）により連立方程式（数３７）が成り立つ。

ここで、吸収物質α、βの吸光係数が既知のとき、連立方程式（数３７）に吸光係数の値を代入して、各体積要素内の２種類の吸収物質の濃度ｃ_α、ｃ_βを求めることができる。従って、被検体内部にγ種類の吸収物質が存在する場合、（γ＋１）種類の波長の照射光を用いて、各吸収物質の濃度分布を求めることができる。
【００２９】
以上の説明した被検体を複数の体積要素に分割した場合において、特に被検体を分割しない場合（即ち、分割の数をｎ＝１とする）にも、以上の方法は適用できることは言うまでもない。この場合、（数２９）はｎ’行１列の行列となり、異なる各種の方向位置に配置される光照射位置と光検出位置の対により計測されるｎ’個の計測値からなる計測値行列と、上記の対におけるそれぞれの平均光路長とから吸収物質の濃度が求められる。もちろん、被検体が複数種類の吸収物質を含む場合にも、被検体を分割せず分割の数をｎ＝１として、（数３７）を適用できる。このようにして求められる吸収物質の濃度は、被検体の２次元的空間（被検体の特定の断面の周囲に２次元的に上記の対が配置されるとき）、あるいは３次元的空間（被検体の外周に３次元的に上記の対が配置されるとき）での平均的な値が求められることになる。
【００３０】
以上の説明において、行列Ｘを求めるためには、行列Ａが既知でなければならない。行列Ａを構成する成分Ｌ_Aijは、ｉ番目の計測位置（光源と検出器の位置関係）におけるｊ番目の体積要素内の平均光路長を表している。即ち、Ｌ_Aijは、各計測位置毎に検出器に到達する光子の光路が既知であれば求めることができる。図６は、光源１から被検体２を通過し、検出器４に到達する光子の光路３を、被検体２が２次元体である例をとって模式的に表した図である。
【００３１】
本発明では、光源から検出器に至る複数の光路を求める方法として、被検体内の光散乱を模擬するモンテカルロ法シミュレーションを用いる。被検体内の光散乱を模擬するモンテカルロ法シミュレーションの例として、Ｂ．Ｃ．Ｗｉｌｓｏｎ等の方法（" A Monte Carlo model for the absorption and flux distributions of light in tissue ", Med. Phys., 10(6), 824ー830(1983)に記載）が知られている。本願発明で使用するモンテカルロ法は、よく知られた方法であるが、この方法の概要を以下に説明する。光源から発した１光子のたどる経路は各散乱点を表す、極座標（ｒ，θ，φ）で示される。各散乱点では、１光子が次に達する散乱点を与える、ｒ、θ、φの各座標値が確率的に定められ、散乱点間の距離は散乱距離Ｌ_Sから、光子の進行方向を定めるための角度方向は、θ方向の散乱角度θ_Sとφ方向の散乱角度φ_Sとから決定される。なお、光子がたどる物体中には吸収物質は含まれないものとする。各散乱毎の散乱距離Ｌ_Sは散乱係数μ_Sを用いて（数３８）で表わされ、
Ｌ_S＝−（ｌｏｇ_e（Ｒ））／μ_S …（数３８）
である。Ｒは一様乱数（０＜Ｒ＜１）である。θ方向の散乱角度θ_Sは、例えばＭｉｅ散乱理論から近似して得られる位相関数の近似式である、Ｈｅｎｅｙ−Ｇｒｅｅｎｓｔｅｉｎの式を用いて、一様乱数Ｒ（０＜Ｒ＜１）と位相関数( phase function )（即ち１回散乱における散乱強度の角度分布を示す関数）のコサイン平均ｇとを使用して（数３９）で表わされ、

である。なお、ｇは、位相関数をＰ_fとし、積分∫を立体角ω＝０〜４πの範囲について行なうものとして、ｇ＝（∫Ｐ_f・ｃｏｓθｄω）／（∫Ｐ_fｄω）により与えれる。一般の生体関連物質では、ｇ≒０．９であることが知られている。特に、ｇ＝０の場合には等方的散乱を表し、（数３９）は、
θ_S＝ｃｏｓ^-1（２Ｒ−１） …（数４０）
となる。φ方向の散乱角度φ_Sは一様な確率分布とみなすことができ、一様乱数Ｒを用いた（数４１）で表わされる。
φ_S＝２πＲ …（数４１）
光子のたどる各散乱点は、このように、一様乱数（０＜Ｒ＜１）Ｒにより、確率的に与えられる（Ｌ_S，θ_S，φ_S）により、順次与えられる。以上が、本願発明で使用するモンテカルロ法の概要である。
【００３２】
このモンテカルロ法による光子のたどる経路シミュレーションを用いて、ｉ番目の計測位置（光照射位置（光源）と光検出位置の組み合わせ）における１光路を求めた時に、ｊ番目の体積要素内での光路長をｐ（ｉ，ｊ）とする。実際の被検体内での光の伝播を実質的に再現するに十分な光路数ｍ本（被検体形状に依存するが、少なくとも１０００本程度とする）を、上記で説明したモンテカルロ法シミュレーションを用いて計算し記憶しておき、各光路毎にｐ（ｉ，ｊ）を求める。ｋ（ｋ＝１〜ｍ）番目の光路におけるｐ（ｉ，ｊ）をｐ（ｉ．ｊ．ｋ）とすると、Ｌ_Aijは（数４２）で計算される。
Ｌ_Aij＝（Σｐ（ｉ，ｊ，ｋ））／ｍ …（数４２）
（数４２）において、加算Σはｋ＝１〜ｍについて行なう。全ての計測位置（ｉ）と、全ての体積要素（ｊ）に関してＬ_Aijを計算すれば、平均光路長の空間分布行列Ａが求められる。
【００３３】
以上説明した本発明の手順は、照射光として連続光を用い被検体を通過する光強度を計測信号とする場合、あるいは、照射光としてパルス光を用い、被検体を通過するパルス光に対応する全光強度を計測信号とする場合であるが、被検体を通過する光強度を時間分解計測し、計測された光強度の１部を計測信号とする場合にも、同様の手順で吸収物質の濃度分布を求めることができる。
【００３４】
図７に示す様な時間スペクトルを有するδ関数状のパルス光を、被検体、例えば生体の所定部分に照射し、照射位置と対向する側の位置で通過光強度を、ストリークカメラを用いて時間分解計測すると、図８に示す様に時間的に拡がった通過光強度の時間スペクトルＩ_d（t）が得られる。ただし、横軸はパルス光が照射されてから通過光が検出されるまでの時間（検出時間）を表す。検出光強度の時間スペクトルＩ_d（t）が、時間的に広がる理由は、被検体が散乱体である場合には、各光子の光路長に応じて検出される時間が変化するためである。図６の様な系で、ｉ番目の各計測位置毎に検出光強度の時間スペクトルＩ_di（t）を計測し、検出時間ｔ₁からｔ₂までＩ_di（t）を積分した値を検出光強度Ｉ_diとして用い計測値行列Ｙを求め、検出時間ｔ₁からｔ₂までの平均光路長の空間分布行列Ａを求めれば、同様の手順で未知数行列Ｘを求めることができる。検出時間ｔ₁からｔ₂までの平均光路長の空間分布行列Ａは、検出時間ｔ₁からｔ₂までの間に検出器に到達する光子の光路のみについて（数４２）を用いて計算すれば求めることができる。このように検出時間ｔ₁からｔ₂まで、時間ゲートを設定してＩ_di（t）を積分した値を求める際に、時間ゲート幅、（ｔ₂−ｔ₁）の大きさを、小さくすればする程、最終的に再構成される被検体中の吸収体の３次元空間分布の空間分解能、及び定量性を向上させることが可能となる。
【００３５】
以下、以上説明した内容を実際に適用した実施例について詳細に説明する。まず簡単のため、本発明を、生体の散乱特性を模擬した立方体の被検体に適用した例を示す。なお、本願発明は被検体が立方体の形状に限定されるものではなく、任意の形状に適用できることは言うまでもない。任意形状の被検体では、被検体の形状を小さい立方体の集合として近似して表してもよいし、被検体の内部を小さい立方体の集合とし、被検体の外面を含む部分を、一部を欠損する小さい立方体で近似して表してもよい。さらに、検査対象物体の特定部分のみの空間、例えば、特定組織部分を含む空間部分、生体表面もしくは表面近傍にある腫瘍を含む空間部分等を、被検体部位として上記と同様の方法で表すことも可能である。
【００３６】
図９は、本実施例に用いた被検体の形状を示す図である。被検体２は１辺２５ｍｍの立方体で、被検体内部の位置を表す座標系（Ｘ０，Ｙ０，Ｚ０）は図９に示す通りである。被検体２の各面を図９に示す様に、Ｐ１〜Ｐ６とする。ここで、Ｐ１はＺ０＝０（ｍｍ）の面、Ｐ２はＸ０＝２５（ｍｍ）の面、Ｐ３はＺ０＝２５（ｍｍ）の面、Ｐ４はＸ０＝０（ｍｍ）の面、Ｐ５はＹ０＝２５（ｍｍ）の面、Ｐ６はＹ０＝０（ｍｍ）の面である。被検体２の散乱特性として、（数３８）で使用する散乱係数を生体に近い値であるμ_S＝１．０〔ｍｍ^-1〕とし、（数３９）で使用する先に説明した位相関数のコサイン平均は、ｇ＝０とした。即ち、１回散乱における散乱光強度の角度分布は一様であり等方的であるとした。
【００３７】
また、図１０に示す様に、被検体２内部を５×５×５の仮想の体積要素に分割し、それぞれの体積要素の位置を（Ｖ_x，Ｖ_y，Ｖ_z）で表し、各体積要素番号ｊは（数４３）で定義する。例えば、図１０の体積要素５の位置は（５，５，１）であり、体積要素番号ｊ＝２５となる。
ｊ＝Ｖ_x＋５（Ｖ_y−１）＋２５（Ｖ_z−１） …（数４３）
図１１は、照射光７を図９のＰ１の面上にある光照射位置６−１から照射していることを示している。本実施例では、全ての光照射位置、及び光検出位置は、被検体表面に接する任意の体積要素の面の重心位置に設定されている。図９の様に、光を照射するｉ番目の面Ｐｉに固有の座標（Ｘ，Ｙ）_Piを定義すると、光照射位置６−１は（２，２）_P1、光照射位置６−２は（２，２）_P2、光照射位置６−３は（３，３）_P5で表される。また、本実施例では、光検出位置を光照射位置が対向する側の被検体外面にある全ての体積要素面の重心位置に設けた。例えば、図１１において、Ｐ１面上の光照射位置６−１から照射光７を照射した時には、対向側のＰ３面上に２５点の光検出位置を設けており、それぞれの光検出位置で検出光強度が計測される（面Ｐｉの定義については図９を参照）。本実施例では、被検体２を合計１２５個の体積要素で分割しているので、それぞれの体積要素内で一様な吸収物質濃度を未知数とすると、求めるべき１２５個の未知数があることになる。未知数を決定するためには、光計測位置（光源と光検出器の位置関係）も１２５組だけ必要となる。そこで、本実施例では、照射点を（２，２）_P1、（２，２）_P2、（４，４）_P3、（４，２）_P4、（３，３）_P5の５点とし、光検出位置は各光照射位置の対向側の面上に２５点設定し、合計１２５（５×２５）組の光計測位置を用意した（面Ｐｉの定義については図９を参照）。
【００３８】
本実施例では、被検体として、吸収物質の位置の異なる２通りを用意した。これら被検体内部の吸収物質の位置を図１２と図１３に示す。図１２と図１３における吸収物質８の位置を（Ｖ_x，Ｖ_y，Ｖ_z）で表すと、図１２の被検体では、吸収物質が（３，３，３）の体積要素内に均一濃度で存在し、図１３の被検体では、吸収物質が（３，３，３）と（１，３，１）の体積要素内に均一濃度でそれぞれ存在している。照射光として２波長λ₁、λ₂のパルス光を用い、各検出位置において通過光強度を時間分解計測した。波長λ₁、λ₂のそれぞれに対する吸収物質の吸光係数は、ε（λ₁）＝１．０〔(ｍｍ・ｍＭ)^-1〕、ε（λ₂）＝０．０〔 (ｍｍ・ｍＭ)^-1)〕とし、吸収物質濃度は全て０．０１〔ｍＭ〕とした。図１４に、先に説明したモンテカルロシュミレーションにより得られた、透過光強度（検出光強度）の時間スペクトルＩ_di（λ₁，ｔ）、Ｉ_di（λ₂，ｔ）の１例を示す。横軸は検出時間を表しており、光子が散乱されずに直接検出位置に到達する時刻を原点としている。なお、図７、図８と同様に、図１４の横軸は、それぞれの横軸の値に、光子の進行する媒体での光速度を乗算して、光子の進行距離として示すこともできる。本実施例では、ｉ番目の計測位置における検出光強度Ｉ_di（λ₁）、Ｉ_di（λ₂）として、各計測位置（ｉ）毎に得られる検出光強度の時間スペクトルを検出時間ｔ₁＝０〔ｐｓ〕から検出時間ｔ₂＝１５００〔ｐｓ〕まで積分した値を用いた。なお、ｔ₂を１５００〔ｐｓ〕より、小さい値に設定することにより、最終的に得られる吸収体の３次元分布の空間分解能、及び定量精度は、更に向上する。
【００３９】
計測された検出光強度から、被検体内部に存在する吸収物質の濃度を画像化する手順の概要を図１５に示す。この画像化の手順は大きく分けると、計測値行列Ｙを生成するｓｔｅｐ１、平均光路長の空間分布行列Ａを生成するｓｔｅｐ２、行列Ｙと行列Ａ逆行列から未知数行列Ｘを導出するｓｔｅｐ３、求められた未知数行列Ｘから被検体内部の吸収物質濃度を求めこの空間分布を表示するｓｔｅｐ４からなる。なお、モンテカルロ法による光路の計算が、ｓｔｅｐ２の処理開始以前に既になされ、光路データが記憶保存されている場合には、記憶保存されている光路データを使用することは言うまでもない。以下、各ｓｔｅｐにおける詳細な手順を、図１６と図１７にフローチャートにより説明する。
【００４０】
〔ｓｔｅｐ１〕： 計測値行列Ｙの生成。
【００４１】
ｓｔｅｐ１−１： 各計測位置（ｉ）毎に通過光強度を計測して、それぞれ検出光強度Ｉ_di（λ₁）、Ｉ_di（λ₂）を求める。
【００４２】
ｓｔｅｐ１−２： 各計測位置（ｉ）毎に照射光に用いた２波長の照射光強度Ｉ_d0（λ₁）、Ｉ_d0（λ₂）を計測する。
【００４３】
ｓｔｅｐ１−３： 全計測位置（ｉ）に関して計測終了後、ｓｔｅｐ１−１、及びｓｔｅｐ１−２で求められた検出光強度Ｉ_di（λ₁）、Ｉ_di（λ₂）、及びＩ_d0（λ₁）、Ｉ_d0（λ₂）を（数１９）に代入してｙ_iを得る。
【００４４】
ｓｔｅｐ１−４： ｓｔｅｐ１−３で得たｙ_i（i＝１〜１２５）から計測値行列Ｙ（数２２）を生成する。
【００４５】
〔ｓｔｅｐ２〕： 平均光路長の空間分布行列Ａの生成。
【００４６】
ｓｔｅｐ２−１： 計測対象である被検体の３次元的な外形形状を計測して、外形形状を表す座標データを演算処理装置に入力し、仮想被検体データを生成する。外形形状を計測する方法としては、光切断法、ＭＲＩ、Ｘ線ＣＴ等その他方法を用いる。本実施例では被検体は１辺２５ｍｍの立方体であり、この外形形状データを入力する。
【００４７】
ｓｔｅｐ２−２： 計測時に設定される計測位置（光照射位置と光検出位置）の座標を入力する。本実施例では先に説明した１２５組の計測位置を演算処理装置に入力する。
【００４８】
ｓｔｅｐ２−３： 計測対象となる被検体の散乱特性として、被検体の散乱係数μ_S、及び位相関数のコサイン平均ｇを演算処理装置に入力する。本実施例では、散乱係数μ_S＝１．０〔ｍｍ^-1〕、指向性を表わす変数としてｇ＝０を用いる。
【００４９】
ｓｔｅｐ２−４： ｓｔｅｐ２−１〜ｓｔｅｐ２−３で入力したパラメータを基に、モンテカルロ法を用いて吸収物質の存在しない被検体内の光散乱を模擬する計算を行なう。設定された各計測位置（ｉ）毎に複数の光路を計算し、各計測位置（ｉ）毎の各体積要素（ｊ）毎における平均光路長Ｌ_Aij（平均光路長の空間分布）をメモリに記憶しておく。この時、通過光強度を時間分解計測する場合には、計測信号として用いる検出時間ｔ₁からｔ₂の間に到達する光子の光路のみから、平均光路長の空間分布Ｌ_Aijを求めて、メモリに記憶しておく。本実施例では、各計測位置毎に３００００本の光路を計算し、検出時間は０〔ｐｓｅｃ〕から１５００〔ｐｓｅｃ〕とした。本ステップで用いるモンテカルロ法に関しては先に概要を説明した。なお、モンテカルロ法による光路の計算が、ｓｔｅｐ２−４の処理開始以前に既になされ、光路データが記憶保存されている場合には、記憶保存されている光路データを使用することは言うまでもない。光路データが記憶保存されている場合の、各計測位置（ｉ）毎の各体積要素（ｊ）毎における平均光路長Ｌ_Aijを求めるための処理方法に関しては、後で詳細に説明する。
【００５０】
ｓｔｅｐ２−５： ｓｔｅｐ２−４で全計測位置（ｉ）と全体積要素内（ｊ）に関して平均光路長Ｌ_Aij（ｉ＝１〜１２５、ｊ＝１〜１２５）を求めた後、平均光路長の空間分布行列Ａ（数２３）を生成する。ただし、計測対象である被検体の形状が予想される場合、あるいは被検体の形状と類似した形状を用いる場合には、あらかじめ平均光路長の空間分布行列Ａを計算しておく。例えば、生体頭部や乳房が被検体である場合には、複数種類の頭部形状や乳房形状について、平均光路長の空間分布行列Ａを計算しておき、実際に計測対象となる生体頭部や乳房と最も近い形状に関して計算した平均光路長の空間分布行列Ａを選択する。また、生体頭部等は、球体あるいは楕円体に近似して計算する。
【００５１】
〔ｓｔｅｐ３〕： 未知数行列Ｘの導出。
【００５２】
ｓｔｅｐ３−１： ｓｔｅｐ１−４で生成された計測値行列Ｙとｓｔｅｐ２−５で生成された平均光路長の空間分布行列Ａの逆行列Ａ^-1を求め（数２６）に代入して、未知数行列Ｘを求める。もちろん（数３０）により未知数行列Ｘを求めてもよい。
【００５３】
ｓｔｅｐ３−２： 未知数行列Ｘの各成分（数２７）より、照射光の各波長に対する吸収物質の吸光係数を用いて、各体積要素毎の吸収物質濃度を求める。
【００５４】
〔ｓｔｅｐ４〕： 画像表示。
【００５５】
ｓｔｅｐ３で得られた吸収物質濃度の空間分布から、任意断層面の吸収物質濃度分布を画像表示する。
【００５６】
図１８から図２３に、上記で説明したフローチャートに従って求められた吸収物質濃度の空間分布の結果を表示する。図１８から図２０は、図１２に示した被検体における吸収物質濃度の空間分布を求めた結果である。図１８はＶ_y＝２の面、図１９はＶ_y＝３の面、図２０はＶ_y＝４の面を表示している。図２１から図２３は、図１３に示した被検体における吸収物質濃度の空間分布を求めた結果である。図２１はＶ_y＝２の面、図２２はＶ_y＝３の面、図２３はＶ_y＝４の面を表示している。図１８から図２３において、ｘ＝、ｙ＝、ｚ＝のそれぞれに続いて記される数は対応する座標の値を示す。以上の計算では、被検体の散乱係数として、生体に近い散乱係数、μ＝１．０〔ｍｍ^-1〕を設定し、（数３９）中の位相関数のコサイン平均をｇ＝０（等方散乱）と仮定している。一般の生体関連物質では、ｇ≒０．９であることが知られているが、上記の計算ではｇ＝０の仮定にもかかわらず、得られた結果は、明確に予め設定された被検体における吸収物質濃度の空間分布を再現しており、吸収物質濃度が異なる５×５×５のサイズの空間が明確に識別されており、が空間分解能も良好である。
【００５７】
図２４に、Ｘ線ＣＴのアルゴリズム（ Filtered Back Projection 法）を用いた結果例を示す。図２４は、外径２５ｍｍ、高さ１００ｍｍを有する円柱（散乱係数１．０ｍｍ^-1である）の中心軸位置に、外径５．ｍｍ、高さ１００ｍｍを有する円柱（空間分解能の評価において広く仮定される完全吸収体とし、吸収系数は無限大とする）の吸収体を配置した被検体に関する再構成画像である。画像の再構成に使用した計測データは、被検体に照射する光をδ関数と見做して、先に概要を説明したモンテカルロ法を使用したシミュレーションにより、図１４に示すような検出光強度の時間スペクトルを生成した。画像の空間分解能を上げるために、時間ゲート幅を小さく設定して、検出光強度の時間スペクトルの０〔ｐｓｅｃ〕から３００〔ｐｓｅｃ〕までの積分信号を計測信号とした。断層像は、計算により円柱の高さ方向の中心位置で Filtered Back Projection 法で求めた。
【００５８】
図２４では、再構成画像の結果を立体図で示しており、高さ方向の値の最大値を１．０、高さ方向の表示の最小間隔を０．０５として示している。図１８〜図２３と図２４との比較から、本発明による画像化方法は従来の方法と比較して明らかに空間分解能が向上していることがわかる。 Filtered Back Projection 法では、散乱の影響を十分考慮することができず、検出光強度の時間スペクトルの積分時間間隔を仮に短くした計測信号を使用しても、Ｓ／Ｎ比が劣化してしまい高い分解能を得ることは困難である。図２４の結果では、外径５ｍｍの吸収体は、空間分解能は半値幅で評価すると１５ｍｍであり、５ｍｍ×５ｍｍ×５ｍｍの大きさの吸収体を再構成した、図１８〜図２３に示す結果では、吸収体は明確に鋭く（シャープ）に検出されている。このように、本発明での空間分解能は、 Filtered Back Projection 法に比較して約３倍以上改善されている。
【００５９】
図２５は、本発明の方法が適用される装置構成を説明する図である。光源９は少なくとも２波長の光を発する光源で、例えばチタンサファイアレーザ、もしくは複数の半導体レーザで構成される。波長の切り替えは、計算機２４で制御する。まず、光源９から１波長（λ₁）の光が選択され、光ファイバー、レンズ等の光学系で構成される光源用導波路１０(以下、特に断りの無い場合、導波路は光ファイバーあるいは光学系で構成される)を通り、光は分波器１１で参照光用導波路１２と光源用光ファイバー１３の２系統に分岐される。ストリークカメラ等の光検出器２０により、参照光用導波路１２に分岐した光は照射光強度Ｉ₀をもつ参照光として計測される。光源用光ファイバー１３に分岐した光は光スイッチ１４に入射する。光スイッチ１４は、被検体２の表面上に配置された複数の光照射用光ファイバー１５の１本を選択し、光源用光ファイバー１３と接続され、光源用光ファイバー１３からの光は被検体２の任意位置から被検体２内部に照射される。
【００６０】
ここで、光照射用光ファイバー１５は、光ファイバー固定ホルダー１６で位置を固定維持されている。被検体２を通過した光は、被検体２の表面上に配置された複数の光検出用光ファイバー１７に入射する。光スイッチ１８は、複数の光検出用光ファイバー１７の中から必要な光ファイバーのみを選択し、検出器用光ファイバー１９と接続する。選択された光検出用光ファイバー１７を通過する検出光は光検出器２０に導かれ、各検出器用光ファイバー１９毎に検出光強度の時間分解信号が計測される。光源の強度信号と、検出光強度の時間分解信号は信号伝送ケーブル２１でＡ／Ｄ変換器２２に転送されデジタル信号に変換される。デジタル信号に変換された光源の強度信号と検出光強度の時間分解信号は、信号伝送ケーブル２３で計算機２４に転送され、計算機内の記憶装置に保存される。
【００６１】
引続き、上記で用いた波長（λ₁）と異なる波長（λ₂）の光を選択し、上記と同様の計測を行なう。１照射位置に関して計測終了後、前回の照射位置あるいは検出位置と異なる位置における光計測をするため、光スイッチ１４及び１８を光スイッチ制御ケーブル２６及び２７を介して計算機２４により制御して、光源用光ファイバー１３と光照射用光ファイバー１５の接続位置及び検出器用光ファイバー１９と光検出用光ファイバー１７の接続位置を切り替える。必要な計測が終了した後、先に説明した図１５、図１６、図１７のフローチャートの手順に従い、吸収物質濃度の空間分布を求め表示装置２５に表示する。
【００６２】
図２６に、光照射用光ファイバー１５と光検出用光ファイバー１７の配置例を示す。被検体の面の番号Ｐｉを図９と同様に定義し、照射位置と検出位置を図９と同様にする場合には、照射用光ファイバー１３を（２，２）_P1、（２，２）_P2、（４，４）_P3、（４，２）_P4、（３，３）_P5の５点に配置し、検出用光ファイバー１７はＰ１〜Ｐ５の各面の各体積要素面の重心位置に配置する。ただし、前記の５点の光照射位置には、光照射用光ファイバー１５と光検出用光ファイバー１７の２本の光ファイバーが配置されるので、それぞれの光ファイバーの先端位置を多少ずらして配置するものとする。光照射用光ファイバー１５、及び光検出用光ファイバー１７には、例えば径１ｍｍの石英ファイバーを用いる。
【００６３】
図２７に、被検体２を生体頭部とする場合の、光照射用光ファイバー１５と光検出用光ファイバー１７の配置例を示す。被検体が生体頭部の場合には、ヘルメット状の光ファイバー固定具２８を用いて光照射用光ファイバー１５、あるいは光検出用光ファイバー１７を固定維持する。光ファイバー固定具２８には、光照射用光ファイバー１５、あるいは光検出用光ファイバー１７を挿入する穴が設けられ、各光ファイバーを任意の穴位置に挿入する。被検体が生体の様な場合にも、仮想の体積要素５を立方体とするが、図２７に示す様に、表面に接する体積要素５は立方体とはならない。
【００６４】
ここで、被検体を１０００個の体積要素に分割する場合には、計測位置（照射位置と検出位置の関係）を少なくとも１０００組用意する必要があり、計測位置を１０００組設定するために、図２５に示す光照射用ファイバー１５と光検出用光ファイバー１７の配置方法としては、様々の方法がある。例えば、光照射用光ファイバー１５を少なくとも２０本、光検出用光ファイバー１７を少なくとも１０００本、それぞれ配置し、少なくとも２０本の光照射用光ファイバー１５から光を照射した時に、少なくとも１０００本の光検出用光ファイバー１７の中から少なくとも５０本を選択して、検出光強度を計測すればよい。この計測は、図２５と同様な装置構成で実現できる。
【００６５】
生体に照射光として近赤外光を用いると、吸収物質としては酸化ヘモグロビン（以下ＨｂＯ₂と略す）と脱酸素化ヘモグロビン（以下Ｈｂと略す）が支配的であり、これら２成分の吸収物質濃度の空間分布を求めることができる。この場合は、照射光として３種類の波長（λ₁，λ₂，λ₃）の光を用いて、各波長毎に光計測を行い、（λ₁，λ₂）の波長を照射光として用いた計測から図１５、図１６、図１７に記載の、ｓｔｅｐ１からｓｔｅｐ３までと同様のフローチャートの手順に従い、未知行列Ｘ、即ちχ（ｊ，λ₁，λ₂）を求め、次に（λ₁，λ₃）の波長を照射光として用いた計測から同様にしてχ（ｊ，λ₁，λ₃）を求め、得られたχ（ｊ，λ₁，λ₂）、χ（ｊ，λ₁，λ₃）から連立方程式（数３７）を立て、連立方程式（数３７）に各波長の光に対する吸収物質の吸光係数εを代入して解き、各体積要素毎のＨｂの濃度ｃ（Ｈｂ，ｊ）と、ＨｂＯ₂の濃度ｃ（ＨｂＯ₂，ｊ）を求める。
【００６６】
図２８に、連立方程式を解く際に用いるＨｂとＨｂＯ₂の吸光係数の波長スペクトルを示す（ S.Wray 等, " Characterization of the near infrared absorption spectra of cytochrome aa3 and hemoglobin for the non-invasive monitoring of cerebral oxygenation "; Biochim. Biophys. Acta 933, 184-192(1988) より数値を抜粋してグラフ化した）。図２８中で、実線がＨｂの吸光係数の波長スペクトルを、点線がＨｂＯ₂の吸光係数の波長スペクトルを、それぞれ示す。
【００６７】
図２８に示すような吸収物質に関する吸光係数の波長スペクトルをあらかじめ計算機内の記憶装置に記憶しておき、照射光として用いる波長に対する吸光係数を求める。頭部以外の生体組織、例えば乳房等の様な部位を被検体とする場合も、専用の光ファイバー固定具を使用し、照射光として近赤外光を用いて、被検体内部のＨｂとＨｂＯ₂の空間分布を求めることができる。
【００６８】
また、被検体中にγ種類の吸収物質が存在する場合には、照射光として（γ＋１）種類の波長の光を用いて、被検体中のγ種類の吸収物質濃度の空間分布を得ることができる。例えば、近赤外光を用いて筋肉中のＨｂとＨｂＯ₂とミオグロビンの空間分布を求める場合には、４種類の波長の照射光を用いればよい。
【００６９】
被検体の一部の領域についての吸収物質濃度の空間分布を求める方法を、図２９を用いて説明する。図２９では、生体頭部を被検体２とし、被検体の一部である注目領域３１内の吸収物質濃度の空間分布を求めることを示している。図２９において、注目領域３１の形状は２５ｍｍ角の立方体で、５×５×５の合計１２５個の体積要素に分割し、各体積要素内の吸収物質濃度を求める。注目領域３１の近傍に配置した光照射用光ファイバー１５から照射光２９を被検体２に照射し、被検体２内で散乱された光が、注目領域３１の近傍に配置された光検出用光ファイバー１７に到達し検出光３０として検出される（この時検出される光は、光照射位置と同一方向に配置された検出器により検出される、反射光である）。但し、光照射用光ファイバー１５から光検出用光ファイバー１７までの光路３は、注目領域３１外の領域を通過する光路もあるので、光が到達する領域３２（以下では、光到達領域と略記する）の内部全体についての吸収物質濃度の空間分布を求めなければならない。光到達領域３２の形状は、あらかじめモンテカルロ法によるシミュレーションや、拡散方程式を用いたシミュレーション等で決定する。
【００７０】
例えば、光照射用光ファイバー１５から光検出用光ファイバー１７に到達する光の９５％が通過する領域を包絡するように、光到達領域３２を定める。
【００７１】
光到達領域３２内の吸収物質濃度の空間分布を求めるには、光到達領域３２を注目領域３１と同様に体積要素に分割して求める方法と、注目領域３１外の光到達領域３２に含まれる領域を１体積要素と扱って求める方法がある。両方法とも体積要素数が異なるだけであり、本質的には同じであるので、本実施例では注目領域３１外の領域を１体積要素として説明する。従って、光吸収物質濃度の空間分布を求めるべき光到達領域３２内には、体積要素数が１２６個存在することになり、異なる計測位置（照射位置と検出位置の位置関係）を少なくとも１２６組用意しなければならない。計測位置の配置方法は無数にあるが、図２９において、例えば光照射用光ファイバー１５を注目領域３１付近に少なくとも９本配置し、光検出用光ファイバー１７を注目領域付近に少なくとも１４本配置すれば、少なくとも１２６（９×１４）組の計測位置を用意できる。図２５と同様な装置構成で、各計測位置毎に被検体２を通過してくる検出光３０の強度を計測し、図１５、図１６、図１７と同様のフローチャートの手順で、光到達領域３２内の吸収物質濃度の空間分布を求めることができる。ただし、照射光として近赤外光を用い、各体積要素中のＨｂＯ₂とＨｂの濃度を求める場合には、照射光として３種類の波長を用いる。また、平均光路長の空間分布行列Ａを計算する際には、光到達領域３２の外部形状が入力される。光到達領域３２の外部形状は、あらかじめ求めることができるので、平均光路長の空間分布行列Ａも、あらかじめ計算してメモリに記憶しておいても構わない。
【００７２】
以下、被検体の任意の光照射位置から任意の光検出位置までの、複数の光路を計算して一旦メモリに光路データとして記憶しておき、記憶された光路データに基づいて、平均光路長の空間分布行列Ａをもとめるための手順（以下では、簡単のため回転移動型モンテカルロ法と略記する）について図を用いて、詳細に説明する。任意の光照射位置から任意の光検出位置までの複数光路は、被検体内の光散乱を被検体の散乱特性を用いて模擬するモンテカルロ法で計算して求めた。従来の光散乱を模擬し光路を計算する従来のモンテカルロ法は、ウィルソン等の行った方法が一般的であり、先に説明した。被検体の任意の光照射位置から任意の光検出位置までの光路を、従来のモンテカルロ法を用いて、必要に応じてその都度計算してもよいが、その計算のために膨大な計算時間を要するため実用的ではない。本発明では、光路データを記憶しておき、これを使用することより、その都度光路を計算する従来方法より、約１００倍以上高速化できる。
【００７３】
以下では被検体形状として、図１２あるいは図１３と同様な立方体を用いて説明する。図３０は、回転移動型モンテカルロ法を説明するために想定した被検体形状と、照射光２９を照射する光照射位置３３と検出光３０を検出する光検出位置３４の相対位置関係を示す図である。図３０において光照射位置３３から光検出位置３４に至るまでの複数の光路を回転移動型モンテカルロ法で求める手順を、代表例として以下説明する。図３０において、立方体の被検体２は、仮想的に５×５×５の体積要素に分割され、各体積要素５は立方体である。回転移動型モンテカルロ法の概念を図３１に示す。図３１は、簡単のため、光照射位置３３と光検出位置３４の位置を含む断面を示しており、この断面内に存在する光路３７が示される。最初に、半径が被検体２の任意の２点をとりその間の距離の最大値（図３０の被検体形状では、頂点３５−１と頂点３５−２の距離）と等しい、球形の仮想散乱体３６を計算機内で発生する。球形の仮想散乱体３６（図３０では図示せず）の中心Ｏから発射される１光子は、被検体２と等価な散乱特性で多重散乱を繰り返し、仮想散乱体３６の境界面まで到達する。ここで与えている散乱特性とは、先に説明した被検体の散乱係数と位相関数で与えられる。
【００７４】
仮想散乱体３６の境界面に到達するまでの光子の散乱点を、光路データとして記憶保存しておくことにより、１光子の仮想散乱体３６中における光路３７を再現することができる。発射された１光子が境界面に到達した後、同様の過程を繰り返して複数の光子に関する光路を得る。統計的な精度を実質的に補償するに十分な光路数が得られた後、図３１に示す様に、被検体２の外形形状を有する仮想被検体３８を発生する。被検体２上の光照射位置３３と光検出位置３４の位置に対応して、仮想被検体３８上に光照射位置３３と光検出位置３４を設定する。続いて、仮想被検体３８上の光照射位置３３を、球形の仮想散乱体３６の中心位置に配置し、球形の仮想散乱体３６の中心Ｏを回転中心として仮想被検体３８を回転させて、仮想被検体３８と光路の相対位置を変化させ、仮想被検体３８上の光検出位置３４と、上記の過程で記憶保存されている光路とが交差する点を求め、中心Ｏからこの交差点に至るまでの光路を再現する散乱点の座標を仮想被検体３８に設定されている座標系に変換する。即ち、球形の仮想散乱体３６の中心Ｏを回転中心として、仮想被検体３８上の光照射位置３３と光検出位置３４を結ぶ線分を回転させて、この線分の先端点である光検出位置３４と記憶保存されている光路とが交差する点を求める。上記の保存されている光路データのそれぞれに関して、同様の操作を繰り返し、仮想被検体３８上の光照射位置３３から光検出位置３４に至るまでの複数の光路が求められる。例えば、図３１中に示す様に、１光路に関して処理が終了した後、仮想散乱体３８を中心Ｏの周りに矢印の方向に回転し、次の光路データに関する処理を行う。求められた仮想被検体３６上の光照射位置３３から光検出位置３４に至るまでの光路データ（各光路に属する散乱点のデータ）から、各体積要素毎の光路長を計算することができ、各体積要素毎の光路長を記憶しておき、被検体２上の光照射位置３３と光検出位置３４の計測位置における、平均光路長の空間分布を（数４２）から計算することができる。以上が、回転移動型モンテカルロ法の２次元での説明であるが、３次元においても同様の処理ができることは言うまでもない。被検体の形状が一般の場合にも、被検体を複数の空間要素に分割して、被検体の外形形状を有する仮想被検体を処理装置において生成して、以上の説明と同様の処理を行なって、任意形状の被検体に関しての、光照射位置と光検出位置とが設定されたとき、この光照射位置と光検出位置とを通過する平均光路長の空間分布を（数４２）から計算することができる。
【００７５】
以上説明した処理の手順をフローチャートとして図３２、及び図３３に示し、以下にその詳細を各ステップ毎に説明する。ただし、ここでは図１７に示されている様に、既に被検体の外形形状を表すデータが入力されており、さらに、被検体の散乱特性も入力されているものとする。
【００７６】
ｓｔｅｐ１： 被検体の表面上に任意の２点を取り、その距離の最大値が半径となるような球の仮想散乱体３６を発生する。仮想散乱体３６の中心から光子を出射し、仮想散乱体境界面から光子が出ていくまでの光路データを保存する。統計精度を実質的に補償するために、少なくとも１０００本の光路データを保存する。仮想散乱体３６の散乱特性は、被検体の散乱特性を反映する。
【００７７】
ｓｔｅｐ２： 体積要素に分割された仮想被検体を発生し、被検体上の光照射位置に対応して、光照射位置が仮想被検体に設定される。仮想被検体上の光照射位置を球形の仮想散乱体中心に配置する。
【００７８】
ｓｔｅｐ３： ｓｔｅｐ１で保存されている任意の光路データから、光路の１本を選択し、被検体上の光照射位置と同位置に設定される仮想被検体上の光検出位置と、光路が交差するように被検体を回転する。この時、回転中心は仮想散乱体の中心（即ち光照射位置）とする。本ステップで任意光路を選ぶ際に、既に選択された光路は重複して選択しない。
【００７９】
ｓｔｅｐ４： ｓｔｅｐ３で光照射位置と光検出位置を通過した光路が、光路中で仮想被検体の境界面と交差するか否かを判断する。交差する場合はｓｔｅｐ３に戻り、次の光路に関して処理を行ない、交差しない場合はｓｔｅｐ５に進む。
【００８０】
ｓｔｅｐ５： 光照射位置から光検出位置への到達時間ｔを総光路長から計算し、検出時間条件ｔ₁、ｔ₂に対して、ｔ₁≦ｔ≦ｔ₂の条件を満たしているか否かを判断する。条件を満たしていない場合は、ｓｔｅｐ３に戻り次の光路データに関して処理を行ない、条件を満たしている場合は、ｓｔｅｐ６に進む。
ｓｔｅｐ６： 光照射位置から光検出位置に至るまでの光路の散乱点の座標を、仮想散乱体の座標系から仮想被検体の座標系に変換した後、各体積要素毎の光路長ｐ（ｉ，ｊ，ｋ）を計算し、記憶保存する。
ｓｔｅｐ７： 保存されている全光路データに関する以上の各ｓｔｅｐの処理が終了したか否かを判断する。計算終了でない場合はｓｔｅｐ３に戻り次の光路に関して処理を行ない、処理終了の場合はｓｔｅｐ９へ進む。
【００８１】
ｓｔｅｐ８： 本ステップまでに処理された光路の累積数ｍは、（数３８）の右辺の分母に相当するｍであり、検出器への到達時間〔ｔ₁〜ｔ₂〕の時間区間に到達する光子数ｍと等しく、ｓｔｅｐ６で保存されている各体積要素毎の光路長ｐ（ｉ，ｊ，ｋ）から、各体積要素毎の累積光路長（（数３８）の右辺の分子）を計算し、光路の累積数ｍで除算する。本ｓｔｅｐでｉ番目の光計測位置（光照射位置と光検出位置の組）における、ｊ番目の体積要素毎の平均光路長Ｌ_Aij（平均光路長の空間分布）を求めたことになる。異なる計測位置に関して平均光路長の空間分布を計算する場合には、ｓｔｅｐ９に進む。
【００８２】
ｓｔｅｐ９： 全計測位置（光照射位置と光検出位置の組）に関して、処理を終了したか否かを判断する。処理終了していない場合は、ｓｔｅｐ２に戻り次の計測位置に関して処理を行なう。処理終了の場合には、ｓｔｅｐ１０に進む。
【００８３】
ｓｔｅｐ１０： 全計測位置に関する平均光路長の空間分布を記憶保存する。
【００８４】
以上のフローチャートに示す手順では、全計測位置に関して必ずしも処理する必要はない。幾何学的に等価な関係にある異なる計測位置（光照射位置と光検出位置の組）が複数ある場合に関しては、その中から１つの計測位置を選択して、光路を求め、得られた光路が、他の幾何学的に等価な計測位置を通過するように座標変換して、処理時間を短縮できる。この処理は、前記のｓｔｅｐ６で行えばよい。例えば、図３０において、光照射位置３３が（３，３）_P1であり、光検出位置３４が（４，３）_P2である計測位置は、光照射位置が（３，３）_P1であり、光検出位置が面Ｐ２と平行な面Ｐ４上の点である（４，３）_P4、等の計測位置と幾何学的に等価と考える（面Ｐｉの定義は図９参照）。
【００８５】
以上が、回転移動型モンテカルロ法の説明であるが、同様の考え方により、図３４に示すように、仮想散乱体３６（十分大きいサイズを有する空間であれば任意の形状でよい）の１点（Ｏ₀）から、発する１光子の複数の散乱点を、光路データとして記憶保存しておき、光路３７を再現することができる（この方法を、移動回転型モンテカルロ法と呼ぶ）。なお、図３４では、図３１と同様に、簡単のため、図３０に示される光照射位置３３と光検出位置３４の位置を含む断面を示しており、この断面内に存在する光路３７が示される。仮想散乱体３６中における仮想被検体３８上の光照射位置３３を、仮想散乱体３６の任意の位置Ｏに配置し、仮想散乱体３６の中心Ｏを回転中心として仮想被検体３８を回転させて、仮想被検体３８と光路３７の相対位置を変化させ、仮想被検体３８上の光検出位置３４と、保存された光路３７とが交差する点を求め、点Ｏからこの交差点に至るまでの光路を再現する散乱点の座標を仮想被検体３８に設定されている座標系に変換し、各体積要素毎の光路長ｐ（ｉ，ｊ，ｋ）を計算し、記憶保存する。次いで、光照射位置３３の位置を点Ｏから、光路の沿って点Ｏ’に移動させて同様の処理を行なう。点Ｏと点Ｏ’との間の距離は、仮想被検体の光照射位置３３と光検出位置３４との間の距離より大きくすることが望ましい。その他の処理に関しては、先に説明した回転移動型モンテカルロ法の場合と同様である。以上が、移動回転型モンテカルロ法の２次元での説明であるが、３次元においても同様の処理ができることは言うまでもない。同様にして、被検体の形状が一般の場合にも、被検体に光照射位置と光検出位置とが設定されたとき、この光照射位置と光検出位置とを通過する平均光路長の空間分布を（数３８）から計算することができる。
【００８６】
【発明の効果】
従来の画像化方法では、空間分解能、定量性、計算時間の点で問題点があったが、本発明の方法では、従来方法に比較して、空間分解能で約３倍以上の改善がなされた３次元的な吸収物質の空間分布を得ることができ、１０〜１００倍程度計算時間が短縮され、局所の最適解に落ち込むことなく解を得ることができ、さらに定量性の改善が可能である。また、回転移動型、もしくは移動回転型モンテカルロ法を用いることにより、従来のモンテカルロ法を適用する方法と比較して、平均光路長の空間分布を約１００倍以上の速度で求めることができる。
【図面の簡単な説明】
【図１】吸収物質を一様に含む被検体中での光路を含む断面図。
【図２】異なる層状の吸収物質からなる被検体中での光路を含む断面図。
【図３】散乱体に吸収物質を含む被検体中での光路を含む断面図。
【図４】散乱体に吸収物質を含む被検体中を複数体積要素に分割したときの各体積要素での光路を模式的に示す図。
【図５】被検体に光源と検出器を配置する例を示す図。
【図６】検出器に到達する光子の光路例を示す図。
【図７】照射光の時間スペクトルの例を示す図。
【図８】被検体を通過し検出される光の時間スペクトルの例を示す図。
【図９】実施例に用いた被検体の形状を示す図。
【図１０】実施例に用いた被検体を体積要素に分割する例を示す図。
【図１１】実施例に用いた光照射位置の配置例を示す図。
【図１２】実施例に用いた被検体内での吸収体の位置を示す図。
【図１３】実施例に用いた被検体内での吸収体の位置を示す図。
【図１４】実施例に用いた計算により求めた検出光の時間スペクトル例を示す図。
【図１５】被検体内の吸収物質濃度の空間分布を求める処理を示すフローチャート。
【図１６】被検体内の吸収物質濃度の空間分布を求める処理を示すフローチャート。
【図１７】被検体内の吸収物質濃度の空間分布を求める処理を示すフローチャート。
【図１８】本発明により得られた吸収物質濃度の空間分布の例を示す図。
【図１９】本発明により得られた吸収物質濃度の空間分布の例を示す図。
【図２０】本発明により得らえた吸収物質濃度の空間分布の例を示す図。
【図２１】本発明により得らえた吸収物質濃度の空間分布の例を示す図。
【図２２】本発明により得られた吸収物質濃度の空間分布の例を示す図。
【図２３】本発明により得られた吸収物質濃度の空間分布の例を示す図。
【図２４】従来方法により得られた吸収物質濃度の空間分布の例を示す図。
【図２５】本発明が適用される装置の構成例を示す図。
【図２６】被検体に光照射用、及び光検出用光ファイバーを配置する例を示す図。
【図２７】生体頭部に光照射用光、及び光検出用ファイバーを配置する例を示す図。
【図２８】酸素化ヘモグロビンと脱酸素化ヘモグロビンの吸光係数の波長スペクトル。
【図２９】本発明を生体の１部の領域に適用する例を示す図。
【図３０】本発明による回転移動型モンテカルロ法を適用する被検体の例を示す図。
【図３１】本発明による回転移動型モンテカルロ法の概念を説明するための図。
【図３２】本発明による回転移動型モンテカルロ法により平均光路長の空間分布を求める手順を示すフローチャート。
【図３３】本発明による回転移動型モンテカルロ法により平均光路長の空間分布を求める手順を示すフローチャート。
【図３４】本発明による移動回転型モンテカルロ法の概念を説明するための図。
【符号の説明】
Ｐ１…被検体の第１面、Ｐ２…被検体の第２面、Ｐ３…被検体の第３面、Ｐ４…被検体の第４面、Ｐ５…被検体の第５面、Ｐ６…被検体の第６面、１…光源、２…被検体、３…光路、４…検出器、５…体積要素、６…光照射位置、７…照射光、８…吸収物質、９…光源、１０…光源用導波路、１１…分波器、１２…参照光用導波路、１３…光源用光ファイバー、１４…光スイッチ、１５…光照射用光ファイバー、１６…光ファイバー固定具、１７…光検出用光ファイバー、１８…光スイッチ、１９…検出器用光ファイバー、２０…光検出器、２１…信号伝送ケーブル、２２…Ａ／Ｄ変換器、２３…信号伝送ケーブル、２４…計算機、２５…表示装置、２６…光スイッチ制御ケーブル、２７…光スイッチ制御ケーブル、２８…光ファイバー固定具、２９…照射光、３０…検出光、３１…注目領域、３２…注目領域外の領域、３３…光照射位置、３４…光検出位置、３５…被検体の頂点、３６…仮想散乱体、３７…光路、３８…仮想被検体。

Claims (4)

1. 被検体に光を照射し、前記被検体内の光吸収物質の濃度の空間分布を画像化する方法であって、
   前記被検体上の光照射位置と光検出位置の複数の対に対してそれぞれ、前記光照射位置から前記光検出位置に至るまでの散乱光の光路をシミュレーションにより求めるステップと、
   前記光路から前記複数の対ごとに平均光路長を算出するステップと、
   複数の前記平均光路長の空間分布を示す空間分布行列を作るステップと、
   前記空間分布行列と、前記光吸収物質の光学定数と、前記光検出位置で検出された光の強度とから、前記光吸収物質の濃度の空間分布を算出するステップと、
   前記光吸収物質の濃度の空間分布を表示するステップとを有することを特徴とする吸収物質濃度空間分布画像化方法。
2. 前記光路は、モンテカルロ法により求められることを特徴とする請求項１記載の吸収物質濃度空間分布画像化方法。
3. 複数の前記光検出位置で検出された光の強度から計測値行列を作るステップを有し、
   前記空間分布行列と、前記光吸収物質の光学定数と、前記計測値行列とから前記光吸収物質の濃度の空間分布を算出することを特徴とする請求項１記載の吸収物質濃度空間分布画像化方法。
4. 前記空間分布行列を記憶手段に保存するステップを有することを特徴とする請求項１記載の吸収物質濃度空間分布画像化方法。

Images