JP4364995B2 - 散乱吸収体内部の光路分布計算方法 - Google Patents
散乱吸収体内部の光路分布計算方法 Download PDFInfo
- Publication number
- JP4364995B2 JP4364995B2 JP2000078534A JP2000078534A JP4364995B2 JP 4364995 B2 JP4364995 B2 JP 4364995B2 JP 2000078534 A JP2000078534 A JP 2000078534A JP 2000078534 A JP2000078534 A JP 2000078534A JP 4364995 B2 JP4364995 B2 JP 4364995B2
- Authority
- JP
- Japan
- Prior art keywords
- time
- medium
- optical path
- light
- path length
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Lifetime
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N21/47—Scattering, i.e. diffuse reflection
- G01N21/4795—Scattering, i.e. diffuse reflection spatially resolved investigating of object in scattering medium
Description
【０００１】
【発明の属する技術分野】
本発明は、散乱吸収体内部の光路分布の計算方法に関する。より詳細には、光入射位置及び光検出位置を測定対象物の表面に沿って移動させて測定対象物の内部の情報を得る装置に適用可能な散乱吸収体内部の時間分解型光路分布の計算方法に関する。
【０００２】
この時間分解型光路分布から、適宜の計測方式における光路分布の計算、及び光減衰量や平均光路長に対する散乱吸収体各部の寄与度（寄与関数あるいは重み関数とも呼ばれる）を計算することができる。
【０００３】
【従来の技術】
従来、内部計測を含む広義の光ＣＴ等の分野においては、重み関数あるいは寄与関数を用いて測定対象物の内部情報を得る計測方法が知られており、これらに適用される重み関数などに関して既に多くの報告がなされている。これらは、例えば、以下の文献に記載されている。
【０００４】
(1) S. Arridge: SPIE Institutes for Advanced Optical Technologies, Vol. IS11, Medical Optical Tomography: Functional Imaging and Monitoring, 35-64(1993); (2) R. L. Barbour and H. L. Graber: ibid. 87-120(1993); (3) H. L. Graber, J. Chang, R. Aronson and R. L. Barbour: ibid. 121-143(1993); (4) J. C. Schotland, J. C. Haselgrove and J. S. Leigh: Applied Optics, 32, 448-5453(1993); (5) Chang, R. Aronson, H. L. Graber and R. L. Barbour: Proc. SPIE, 2389, 448-464(1995); (6) B. W. Pogue, M. S. Patterson, H. Jiang and K. D. Paulsen: Phys. Med. Biol. 40, 1709-1729(1995); (7) S. R. Arridge: Applied Optics, 34, 7395-7409(1995); (8) H. L. Graber, J. Chang, and R. L. Barbour: Proc. SPIE, 2570, 219-234(1995); (9) A. Maki and H. Koizumi: OSA TOPS, Vol.2, 299-304(1996); (10) H. Jiang, K. D. Paulsen and Ulf L. Osterberg: J. Opt. Soc. Am. A13, 253-266(1996); (11) S. R. Arridge and J. C. Hebden: Phys. Med. Biol. 42, 841-853(1997); (12) S. B. Colak, D. G. Papaioannou, G. W. 't Hooft, M. B. van der Mark, H. Schomberg, J. C. J. Paasschens, J. B. M. Melissen and N. A. A. J. van Astten: Applied Optics, 36, 180-213(1997)。
【０００５】
しかしながら、上記の文献（１）〜（１２）に開示されている光ＣＴ等の内部情報計測方法及びそれらに適用されている重み関数は以下のような問題を包含しているため、計測精度及び使い勝手に大きな問題が生じていた。また、実際に、上記分野において十分な性能をもつ光ＣＴの実用例は、まだ報告されていない。
【０００６】
すなわち、従来の光ＣＴ等における第１の問題は、媒体内の光子移動解析あるいは光子移動モデルが、輸送方程式に拡散近似を適用した光拡散方程式に基いているということにある。つまり、拡散近似は媒体中のフォトンの平均自由光路長に対して十分な大きさの媒体でのみ成立するから、比較的小さい媒体、内部の複雑な形状の組織、及び複雑な形状の媒体を取り扱うことができないという問題がある。また、拡散近似は等方散乱を前提としているので、非等方な散乱特性をもつ生体組織のような測定対象物の計測に適用すると、非等方散乱を等方散乱で近似したことに起因する無視できない誤差が発生するという問題もある。
【０００７】
更に、拡散方程式などの微分方程式は、解析的あるいは有限要素法など、いづれの数値計算手法を用いても、前もって境界条件（媒体形状や界面での反射特性など）を設定した後、解を求めなけらばならないという厄介な問題がある。つまり、通常、生体組織のような測定対象物では、境界条件が計測すべき場所や計測に用いる光の波長などによって変わるから、これらの影響を補正して精度を向上させようとすると、境界条件が変わるたびに複雑な計算をやり直す必要が生じ、結局、計算時間が極端に長くなるという大きな問題がある。
【０００８】
第２の問題は、測定対象物の内部情報を演算する際に、狭義の重み関数（寄与関数とも呼ばれる）、つまり媒体のインパルス応答を構成する光子集合に対する平均光路長（加重平均）或いはそれに相当する位相遅れ（周波数領域での計測）を適用していることにある。この場合、狭義の重み関数（平均光路長或いはそれに相当する位相遅れ）は、吸収係数や吸収分布に依存して変化するため、その取り扱いが極めて複雑になる。実際にこのような依存性を考慮した計算法を用いると、繰り返し回数の多い演算が不可避になるため、演算時間が極めて長くなり、実用的でないという問題が生じる。そこで、通常は、吸収係数や吸収分布に対する依存性を無視するが、このような近似に起因して誤差が増大するという重大な問題が生じる。
【０００９】
また、先に述べた重み関数の吸収係数や吸収分布に対する依存性に起因する誤差を低減するために、適宜の吸収があるときの重み関数を計算して用いるという方法がある。しかし、吸収があるときの重み関数の計算に要する時間は、吸収が無いときの重み関数の計算に要する時間に比較して、極端に長くなるという問題が発生する。
【００１０】
以上のようなことから、従来の狭義の重み関数（平均光路長或いはそれに相当する位相遅れ）を使用した計測方法は実用的でないと結論されている。
【００１１】
更に、上記の狭義の重み関数を適用する場合とは別に、輸送方程式の近似式や光拡散方程式に摂動理論を適用して、信号光と散乱吸収媒体の光学特性との関係を用いて測定対象物の内部情報を得る計測方法がある。しかしながら、この方法では、非線形効果（２次以上の項）の取り扱いが極めて複雑になる。この際、2次以上の項を含めた計算は、理論的にはコンピュータで計算することが可能であるが、現在の最速のコンピュータを用いても演算時間が膨大なものになり、実用化は不可能である。このため、常套的には２次以上の項を無視している。そのため、この方法では、複数個の比較的強い吸収領域が存在する媒体の光ＣＴ画像を再構成する際などで、吸収領域間の相互作用に起因する大きな誤差が発生するという問題が生じていた。
【００１２】
以上のように、従来の光ＣＴ等の測定対象物の内部情報を得る計測方法は、十分な精度の再構成画像が得られず、空間解像度、画像ひずみ、定量性、測定感度、所要計測時間などに大きな問題があるため、実用化が困難であった。
【００１３】
本発明者は、以上のような状況を打破するため、次のことが重要であると考えて一連の研究を推進してきた。すなわち、測定対象物の内部情報を得る計測方法、特に光ＣＴを実現させるために重要なことは、強い散乱媒体である生体組織の中を移動する光の振る舞いの詳細を明らかにし、検出される信号光と吸収成分を含む散乱媒体（散乱吸収体）の光学特性との関係をより厳密に記述し、信号光、並びに、信号光と散乱媒体の光学特性との関係を利用して光ＣＴ画像を再構成する新しいアルゴリズムを開発するということである。そして、本発明者は、散乱媒体中を移動する光の振る舞いを解析する際に、マイクロ・べア・ランバート則（Microscopic Beer-Lambert Law、以下「ＭＢＬ」という）を適用し、光ＣＴ画像の再構成に対して光路分布、つまり光子がどこを通過したかという情報を利用することを考えた。
【００１４】
そして、本発明者は、ＭＢＬに基づいて、散乱吸収体に対応する光子伝播モデル（有限格子モデル）を提案すると共に散乱吸収体の光学特性と信号光との関係を表す解析式を導出し、散乱吸収体内の光の振る舞いを解析する方法を開発してきた。これらの結果を本発明者らは、例えば以下の文献に掲載している。
【００１５】
(13) Y. Tsuchiya and T. Urakami: "Photon migration model for turbid biological medium having various shapes", Jpn. J. Appl. Phys. 34. Part2, pp. L79-81(1995); (14) Y. Tsuchiya and T. Urakami, "Frequency domain analysis of photon migration based on the microscopic Beer-Lambert law", Jpn. J. Appl. Phys. 35, Part1, pp.4848-4851(1996); (15) Y. Tsuchiya and T. Urakami: "Non-invasive spectroscopy of various]y shaped turbid media like human tissue based on the microscopic Beer-Lambert law", OSA TOPS, Biomedical Optical Spectroscopy and Diagnostics 1996, 3, pp.98- 100(1996); (16) Y. Tsuchiya and T. Urakami: "Quantitation of absorbing substances in turbid media such as human tissues based on the microscopic Beer-Lambert law". Optics Commun. 144, pp.269-280(1997); (17) Y. Tsuchiya and T. Urakami: "Optical quantitation of absorbers in variously shaped turbid media based on the microscopic Beer-Lambert law: A new approach to optical computerized tomography", Advances in Optical Biopsy and Optical Mammography (Annals of the New York Academy of Sciences), 838, pp.75-94(1998); (18) Y. Ueda, K. Ohta, M. Oda, M. Miwa, Y. Yamashita, and Y. Tsuchiya: "Average value method: A new approach to practical optical computed tomography for a turbid medium such as human tissue", Jpn. J. Appl. Phvs. 37, Part1. 5A, pp.2717-2723(1998); (19) 土屋 裕;「マイクロ・ベール・ランバート則と平均値法に基づく光ＣＴ画像の再構成」,O plus E, Vol.21, No.7, 814-821; (20) H. Zhang, M. Miwa, Y. Yamashita, and Y. Tsuchiya: "Quantitation of absorbers in turbid media using time-integrated spectroscopy based on microscopic Beer-Lambert law", Jpn. J. Appl. Phys. 37, Part1, pp.2724-2727(1998); (21) H. Zhang, Y. Tsuchiya, T. Urakami, M. Miwa, and Y. Yamashita: "Time integrated spectroscopy of turbid media based on the microscopic Beer-Lambert law: Consideration of the wavelength dependence of scattering properties", Optics Commun. 153, pp.314-322(1998); (22) Y. Tsuchiya, H. Zhang, T. Urakami, M. Miwa, and Y. Yamashita:"Time integrated spectroscopy of turbid media such as human tissues based on the microscopic Beer-Lambert law", Proc. JICAST'98/CPST'98, Joint international conference on Advanced science and technology, Hamamatsu, Aug.29-30, pp.237 -240(1998); (23) H. Zhang, T. Urakami, Y. Tsuchiya, Z. Liu, and T. Hiruma: "Time integrated spectroscopy of turbid media based on the microscopic Beer-Lambert law: Application to small-size phantoms having different boundary conditions", J. Biomedical Optics. 4, pp.183-190(1999);(24) Y. Tsuchiya, Y. Ueda, H. Zhang, Y. Yamashita, M. Oda, and T. Urakami: "Analytic expressions for determining the concentrations of absorber in turbid media by time-gating measurements", OSA TOPS, Advances in Optical Imaging and Photon Migration, 21, pp.67-72(1998); (25) H. Zhang, M. Miwa, T. Urakami, Y. Yamashita, and Y. Tsuchiya: "Simple subtraction method for determining the mean path length traveled by photons in turbid media", Jpn. J. Appl. Phys. 37, Part1, pp.700-704 (1998)。
【００１６】
このＭＢＬは「任意の散乱吸収分布をもつ媒体中において微視的にみたとき、吸収係数がμaである部分を伝播する光子は、伝播する長さｌのジグザグ光路に沿って吸収により指数関数的に減衰し、その光子の生存率は媒体の散乱特性や境界条件とは無関係な値exp(-μal)になり、吸収に起因する減衰量はμalになる」という法則であり、例えば下記式で表現される。
【数１】
ここで、f(l)は、実測される各光子の光路長ｌの分布関数、f0(l)は吸収の存在しないときの分布関数である。
【００１７】
この（１）式は、任意の散乱吸収分布をもつ媒体において、吸収と散乱が独立の事象であること、及び吸収に関して重ね合わせの原理が成り立つこと、すなわち、散乱吸収分布をもつ媒体が多成分系である場合において全体の吸光度はそれぞれの成分の吸光度の和で与えられることを示している。
【００１８】
本発明者らは上記（１）式で表現されるＭＢＬから、散乱媒体の各種光応答と光学特性との関係を表す式を導出し、散乱媒体の各種の応答が吸収に依存する減衰項と吸収に依存しない項とに分離して記述できることを明らかにした。このことから、多波長の光を用いてこの減衰項を分光計測すれば、吸収係数と、吸収成分に固有の吸光係数との関係から、散乱媒体内の吸収成分の濃度の絶対値を定量することができるようになった。
【００１９】
以上のようなＭＢＬに基づく計測方法は原理的に、媒体形状、境界条件及び散乱の影響を受けないという大きな特長があり、非等方散乱媒体や小さい媒体への適用も可能である。ただし、ここで計測される吸収係数や吸収成分の濃度は、吸収が均一な散乱媒体に対して正しい計測ができるが、吸収分布が均一でない不均一媒体に対する計測では光が通過した部分の平均吸収係数や平均濃度（光路分布に対する加重平均）となる。従って、高精度の光ＣＴを実現するためには、光路分布を知る必要が生じた。
【００２０】
そこで、本発明者らは、ＭＢＬに基づく散乱媒体内の光の振る舞いの解析方法をさらに発展させて、散乱や吸収が不均一に分布する不均一散乱媒体(不均一系)に対しても適用可能にして、媒体の形状や散乱特性などに影響されずに散乱媒体内における吸収成分の濃度分布を定量計測する種々の計測方法を開発してきた。
【００２１】
本発明者らがこれまでにＭＢＬに基づいて導出した具体的な計測方法としては、例えば、インパルス応答を利用する時間分解計測（ＴＲＳ）、インパルス応答の時間積分を利用する時間積分計測（ＴＩＳ）、インパルス応答の時間ゲート積分を利用する時間分解ゲート積分計測（ＴＧＳ）、周波数領域での位相変調計測（ＰＭＳ）、更には、平均値法に基づく光ＣＴ画像の再構成法（ＡＶＭ）などがある。これらの計測方法は、特願平１０−１４４３００、上記文献(7)、及び上記文献(15)〜(25)に記載されている。
【００２２】
上記のような新しい計測方法の出現によって計測精度は飛躍的に向上したが、光ＣＴ画像の再構成法では、依然として、従来型の重み関数、つまり平均光路長あるいはそれの相当する位相遅れが用いられていた。具体的には、有限格子モデルで記述した均一な吸収係数μνを有する仮想媒体のvoxe1ｉの吸収係数を少しだけ変化させ、その光出力をモンテカルロ計算、経路積分（Path Integral）や輸送方程式の数値計算、光拡散方程式の数値計算などを利用して計算し、変化前後の光出力の差、すなわち、均一な吸収係数μνを有する仮想媒体を基準として、光強度、光路長、減衰量などの偏差に対応するvoxe1ｉの重み関数（あるいは寄与関数）Ｗｉを求めている。
【００２３】
【発明が解決しようとする課題】
しかしながら、有限格子モデルを用いて記述した仮想媒体のvoxe1ｉの吸収係数を変化させることにより各種の計測方法に対応した重み関数（あるいは寄与関数）を求める上記従来の計測方法では、voxe1の総数に等しい回数の光拡散方程式等の順問題計算を行う必要があるので、この計算時間が極めて長くなり、結局、計測時間が非常に長くなるという問題があった。この場合、従来の狭義の重み関数（平均光路長或いはそれに相当する位相遅れ）を使用した計測方法では、計測精度が十分でないという問題もあった。更に、この計測精度を向上させようとすると、さらに計算時間が長くなり、結果として、さらに計測時間が長くなるという問題があった。
【００２４】
本発明は、上記従来の課題に鑑みて成されたものであり、生体などの不均一系の散乱吸収体内部の光路分布を精密、且つ迅速に算出することのできる方法を提供することを目的とする。なお、このようにして算出された光路分布から、散乱体内部の吸収成分の濃度分布を定量計測するための各種の具体的な計測方法に対応する重み関数Ｗｉを迅速に計算することができる。
【００２５】
【課題を解決するための手段】
本発明者は、上記目的を達成すべく鋭意、前記ＭＢＬに基づく解析法等の研究を進めた結果、Ｎ（１≦Ｎ）個のvoxelに分割された被計測媒体のvoxelｉ内の時間分解型光路長liが、この被計測媒体と同一の散乱特性と境界条件をもつが吸収がないと仮定した仮想媒体のvoxelｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）内の時間分解型光路長liと一致すること、この仮想媒体のvoxelｉ内の時間分解型光路長liは、仮想媒体のvoxelｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）内の光子存在確率等で記述することできること、更に、この仮想媒体のvoxelｉ内の時間分解型光路長liは、モンテカルロ計算、あるいは、経路積分、輸送方程式、または光拡散方程式の数値計算などによって仮想媒体内の光子移動を計算した結果を利用して、直接的に然も迅速に計算できることを見出し、本発明に到達した。
【００２６】
なお、従来のモンテカルロ計算は、有限格子モデル内のvoxelｉ内における散乱事象の発生回数を計数することを基本原理としており、従来の重み関数に対してもこの概念が踏襲・適用されている。この際、voxelｉ内の散乱特性が均一であれば、散乱回数と平均自由光路長との関係から光路長を推定することができるが、voxelｉ内の散乱特性が不均一であるときは、散乱回数と平均自由光路長との関係が散乱分布に依存して複雑になるから、本発明の目的である時間分解型光路長を精密に求めることができない。従って、本発明ではこのような従来の問題点を回避するため、媒体のvoxelｉ内の光子存在確率等からのvoxelｉ内の時間分解型光路長liを算出するという新しい手法を用いている。
【００２７】
すなわち、本発明の散乱吸収体内部の光路分布計算方法は、
計測波長領域において、散乱吸収体である被計測媒体を、所定の大きさを有すると共に一様な吸収分布を有するとみなせるＮ（１≦Ｎ）個のvoxe1に分割する第１ステップと、
前記Ｎ個のvoxe1に分割された前記被計測媒体に対して、形状、境界条件及び散乱分布が前記被計測媒体と同一であり、屈折率分布が一様で、且つ吸収がないとみなせる仮想媒体を仮定する第２ステップと、
前記仮想媒体の表面に、インパルス光が入射する光入射位置と、前記光入射位置から入射して前記仮想媒体内を伝播した前記インパルス光の時間ｔにおけるインパルス応答s(t)を検出する光検出位置とを設定する第３ステップと、
前記仮想媒体内へ前記光入射位置からインパルス光を入射する場合に、前記光検出位置において検出される前記インパルス応答s(t)を算出する第４ステップと、
前記光入射位置から前記仮想媒体内にインパルス光を入射する場合に、前記光検出位置において時間ｔの時点で検出される前記インパルス応答s(t)を構成する光子集合が、前記光検出位置において検出される以前の時間t'(0≦t'≦t)の時点で前記仮想媒体内の任意のvoxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）内に存在していた光子存在確率密度Ｕi(t',t)を算出する第５ステップと、
前記光子存在確率密度Ｕi(t',t)を使用して前記光検出位置において時間ｔの時点で検出される前記インパルス応答s(t)を構成する光子集合が前記voxe1ｉ内に存在していた光子存在累積確率Ｕi(t)を算出する第６ステップと、
前記光子存在累積確率Ｕi(t)及び前記インパルス応答s(t)を使用して前記voxe1ｉの時間分解型光路長liを算出する第８ステップと、
を備えることを特徴する方法である。
【００２８】
従来は、有限格子モデルで記述した仮想媒体のvoxeliの吸収係数を少し変化させて、吸収係数変化の前後の出力光強度の比、或いは減衰量の差等の計算値から各種の重み関数Ｗiを計算していた。この方法では、voxe1の総数Ｎに等しい回数の光拡散方程式等の順問題計算を行う必要があるので計算時間が長くなり、結局、計測時間が非常に長くなるという大きな問題があったが、上記ではこの問題が解決されている。更にその上、上記では吸収のない媒体について時間分解型光路長liを計算しているから、従来の吸収がある媒体に対する光拡散方程式等の順問題計算に長い時間を要したという問題も、解決されている。
【００２９】
また、本発明の方法では、一対の光入射位置と検出位置とからなる組に対して、全voxelのパラメータが並列に計算されるので、更に一層、計算時間を短縮することができる。つまり、本発明の散乱吸収体内部の光路分布計算方法は、時間分解型光路長liを直接精密、且つ迅速に計算することをはじめて見出したものである。
【００３０】
具体的にvoxe1ｉ内の時間分解型光路長liを求めるには、後述するように、有限格子モデル内に光入射位置からインパルス光を入射する場合に、光検出位置において時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光検出位置において検出される以前の時間t'(0≦t'≦t)の時点にvoxe1ｉに存在する光子存在確率密度Ｕi(t',t)を考える。
【００３１】
この新しい時間変数t'を導入して記述される時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)は、相反原理(reciprocity theorem)を適用することにより下記式のように表現できる。
【数２】
［式中、Ｕpi(t')は、時間t'＝０において光入射位置ｐから入射した光子が時間t'＝t'においてvoxe1ｉ内に存在する光子存在確率密度を示し、
Ｕqi(t-t')は、時間t-t'＝０において光検出位置ｑから入射した光子が時間t-t'においてvoxe1ｉ内に存在する光子存在確率密度を示す。］
【００３２】
上記のＵi(t',t)は、（２．２）式で示されるように二つの順問題を計算することにより算出することができる。この時、Ｕpi(t')は、有限格子モデルの光入射位置ｐからインパルス光を入射して計算することができる。また、Ｕqi(t-t')も有限格子モデルの光検出位置ｑからインパルス光を入射して計算することができる。なお、（２．２）式で示される二つの順問題の具体的な計算方法については後に詳細に説明する。
【００３３】
更に本発明においては、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて時間分解型光路長liの時間に対する加重平均Ｌi を迅速に算出することができる。
【００３４】
従来は、時間分解型光路長liの具体的な計算方法が知られていなかったため、最初に時間分解型光路長liを求めて、その後で、時間分解型光路長liの加重平均Li及び重み関数Ｗiを算出するという発想がなかった。実際、時間分解型光路長liを直接計算して利用する計測方法の例はこれまでに報告されていない。
【００３５】
更に本発明においては、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の各種の具体的な計測方法に対応する重み関数Ｗiを迅速に計算することができる。
【００３６】
つまり、被計測媒体及びこれに対応する仮想媒体をＮ個のvoxelに分割した有限格子モデルを利用する場合において、ＭＢＬに基づくＴＲＳ、ＴＩＳ、ＴＧＳ、ＰＭＳ、ＡＶＭ等の具体的な計測方法において、被計測媒体内のvoxel i内の時間分解型光路長liと各種の計測法で計測される計測値および対応する重み関数との関係は、既に本発明者らのよって明らかにされている。
【００３７】
例えば、減衰量に対する重み関数は、ＴＲＳでは時間分解型光路長liそのものに等しく、ＴＩＳでは時間分解型光路長liの時間に対する加重平均Li（平均光路長Liともよばれる）や時間分解型光路長liの分散等を用いて記述される。より具体的には、減衰量と平均光路長を利用するＴＩＳやＡＶＭでは、被計測媒体内のvoxeli内の時間分解型光路長liと、その加重平均である平均光路長Li、更には時間分解型光路長liの分散等との関係を用いて記述される重み関数Ｗｉを用いるが、これらは後述するような（７）、（９）、（１２）式で表現されることが知られている。
【００３８】
ここで、本発明の散乱吸収体内部の光路分布計算方法における「時間分解型光路分布」とは、図1に示すように、３次元の被計測媒体をＮ（１≦Ｎ）個のvoxe1（体積素）に分割した有限格子モデルとし、この有限格子モデルの光入射位置ｐからインパルス光を入射した場合に、有限格子モデルの光検出位置ｑにおいて計測されるインパルス応答の光路長がｌとなる光子集合の光路分布を示す。具体的には、この時間分解型光路分布はvoxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）の時間分解型光路長liを要素とするＮ次元ベクトル［li］で表される。上記の定義から、この時間分解型光路分布は光路長ｌの関数、つまり時間tの関数である。各種の具体的な計測法における出力光信号の減衰量は、この時間分解型光路長liとvoxe1ｉの吸収係数の積の関数で表される。また、この時間分解型光路長liの時間に対する加重平均Liは、voxe1ｉの平均光路長とよばれ、これを要素とするＮ次元ベクトル［Li］は平均光路長分布とよばれる。
【００３９】
更に、重み関数は注目している物理量に対する吸収情報の寄与度として定義される。例えば、ＴＲＳで得られる時間分解型の減衰量に対する重み関数は、吸収や吸収分布に無関係となり、［時間分解型減衰量］＝［重み関数］・［吸収係数差］となり、この重み関数は上記で求めた時間分解型光路長liに等しい。また、インパルス応答の時間積分（ＴＲＳ）で得られる減衰量に対に重み関数は［ＴＲＳ減衰量］＝［重み関数］・［吸収係数差］で定義され、平均光路長差に対に対する重み関数は［平均光路長差］＝［重み関数］・［吸収係数差］で定義される。ＴＲＳの場合、これらの重み関数は吸収係数の関数になる。
【００４０】
このような重み関数を導入することによって、光ＣＴにおける画像再構成のための行列式が行列の積、つまり［重み関数］・［吸収係数差］として記述され、この式から［吸収係数差］すなわち吸収分布を計算することができる。
【００４１】
なお、上記の被計測媒体と仮想媒体との散乱特性は不均一であってもよく、順問題の具体的な計算には、モンテカルロ計算、経路積分（Path Integral）や輸送方程式の数値計算、光拡散方程式の数値計算などを利用することができる。
【００４２】
ここで、本発明の散乱吸収体内部の光路分布計算方法により求めた時間分解型光路長liが適用される「所定の計測方式」とは、特に限定されるものではなく、例えば、発明者らがＭＢＬに基づいて導出した時間分解計測（ＴＲＳ）、時間積分計測（ＴＩＳ）、時間分解ゲート積分計測（ＴＧＳ）、更には周波数領域での位相変調計測（ＰＭＳ）、平均値法に基づく光ＣＴ画像の再構成法（ＡＶＭ）等の散乱吸収体の内部計測方法でもよく、先に文献(1)〜(12)に示したようなＭＢＬに基づかない計測方法でもよい。
【００４３】
また、本発明の散乱吸収体内部の光路分布計算方法は、
計測波長領域において、散乱吸収体である被計測媒体を、所定の大きさを有すると共に一様な吸収分布を有するとみなせるＮ（１≦Ｎ）個のvoxe1に分割する第１ステップと、
前記Ｎ個のvoxe1に分割された前記被計測媒体に対して、形状、境界条件及び散乱分布が前記被計測媒体と同一であり、屈折率分布が一様で、吸収がないとみなせる仮想媒体を仮定する第２ステップと、
前記仮想媒体の表面に、インパルス光が入射する光入射位置と、前記光入射位置から入射して前記仮想媒体内を伝播した前記インパルス光の時間ｔにおけるインパルス応答ｓ(t)を検出する光検出位置とを設定する第３ステップと、
前記光入射位置から前記仮想媒体内にインパルス光を入射する場合に、前記光検出位置において時間ｔの時点で検出される前記インパルス応答ｓ(t)を構成する光子集合が、前記光検出位置において検出される以前の時間t'(0≦t'≦t)の時点で前記仮想媒体内の任意のvoxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）内に存在していた光子存在確率密度Ｕi(t',t)を算出する第５ステップと、
前記光子存在確率密度Ｕi(t',t)を使用して前記光検出位置において時間ｔの時点で検出される前記インパルス応答ｓ(t)を構成する光子集合が前記voxe1ｉ内に存在していた光子存在累積確率Ｕi(t)を算出する第６ステップと、
前記光子存在累積確率Ｕi(t)の全voxe1に対する総和Ｕ(t)を算出する第７ステップと、
前記光子存在累積確率Ｕi(t)及び前記光子存在累積確率Ｕi(t)の全voxe1に対する総和Ｕ(t)を使用して前記voxe1ｉの時間分解型光路長liを算出する第１１ステップと、
を備えることを特徴とする方法である。
【００４４】
上記の場合にも、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて時間分解型光路長liの時間に対する加重平均Ｌi を迅速に算出することができる。また、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の各種の具体的な計測方法に対応する重み関数Ｗiを迅速に計算することができる。
【００４５】
【発明の実施の形態】
以下、図面を参照しながら本発明による散乱吸収体内部光路分布計算方法の好適な実施形態について詳細に説明する。なお、以下の説明では、同一または相当部分には同一符号を付し、重複する説明は省略する。また、散乱吸収体内を散乱されつつ伝播する光子に関しては３次元座標を用いて考える必要があるが、以下では説明を簡単にするために場合により２次元座標を用いて説明する。
【００４６】
先ず、本発明の原理について説明する。
［本発明の原理］
先に述べたように、本発明者らがＭＢＬに基づいて開発してきた散乱媒体内の光の振る舞いを解析する解析法は、散乱や吸収が不均一に分布する不均一散乱媒体(不均一系)に適用することができる。そして、ＭＢＬに基づく解析法により、不均一散乱媒体の種々の光応答（光出力）が吸収に依存する減衰項と吸収に依存しない項とに分離して記述できること、更に、この光応答の減衰項が被計測媒体内の光路分布と吸収分布で記述できることが演繹される。従って、前もって被計測媒体内の光路分布を知ることができれば、被計測媒体内部の吸収成分の濃度分布を定量する光ＣＴなどが可能になる。つまり、光ＣＴを実現するには、被計測媒体内の光路分布を知ることと、光路分布を迅速に計算する方法を開発することが重要である。
【００４７】
本発明は、被計測媒体及びこれに対応する仮想媒体の有限格子モデルにおいて、以下の（i）〜（v）の事柄が成立することを出発点としている。
【００４８】
すなわち、（i）被計測媒体内の屈折率分布が一様であると仮定すれば、被計測媒体のインパルス応答（時間分解波形）では、光路長をｌとし、屈折率で決まる被計測媒体中の光速度をｃ、時間をｔとすれば、常にｌ＝ｃｔが成立する。ここで、被計測媒体内の屈折率分布が一様であると仮定するのは、被計測媒体として生体などの不均一散乱吸収体を選択した場合には、生体の主成分は水であるためその屈折率分布は一様であると見なしてよいからである。従って、以下の説明においては特にことわりがない限り屈折率分布が一様な被計測媒体を考える。
【００４９】
（ii）散乱と吸収が独立であることから、散乱と吸収が任意に分布する不均一媒体に対して光入射位置、及び光検出位置を定めると、インパルス応答を構成する(インパルス応答として検出される)光子集合の光路分布、つまり時間分解型光路分布が一意的に決まる。
【００５０】
（iii）上記の光路分布は、時間の関数となり、入射光子数（つまり光子密度）、被計測媒体の吸収や吸収分布に無関係であり、被計測媒体に吸収がないと仮定したときの光路分布に等しい。そして、３次元の被計測媒体を所定の大きさを有するＮ（１≦Ｎ）個のvoxe1（体積素）に分割した場合において、光路長がｌとなる光子集合の光路分布は、各voxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）の光路長liを要素とするＮ次元ベクトル［li］で表すことができる。
【００５１】
（iv）この時間分解型光路分布［li］は、被計測媒体と同一の散乱特性（不均一であってよい）と境界条件をもつが、屈折率分布が一様で、且つ吸収がない仮想媒体を仮想して計算することができる。
【００５２】
（v）各種の光ＣＴ画像再構成アルゴリズムに用いる各種の重み関数Ｗiを上記に定義した各voxe1ｉの時間分解型光路長liを用いて記述することができる。
なお、本発明者は、この時間分解型光路分布liを用いて記述される重み関数Ｗiと被計測媒体に対する実測値を用いて不均一媒体内の吸収成分の濃度分布を計測する光ＣＴを提案した（特願平１０−１４４３００、文献(19)）。
【００５３】
図１は被計測媒体である不均一散乱吸収体の内部光子移動を解析するための有限格子モデルを示す。
【００５４】
図１に示す有限格子モデルは、被計測媒体に対して形状、散乱分布、及び境界条件が同一であり、屈折率分布が一様、且つ吸収がないとして仮定される仮想媒体を示したものである。この際、３次元の被計測媒体をＮ（１≦Ｎ）個のvoxe1（体積素）に分割すると共にそれぞれのvoxe1に番号ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）を付けているので、被計測媒体に対応して仮定される仮想媒体も同様にＮ個のvoxe1を有している。
【００５５】
また、図１に示す仮想媒体の有限格子モデルのvoxe1ｉの大きさ及び形状は、対応する被計測媒体のvoxe1ｉ内の吸収分布が一様であると見なせるか、一様であると見なしてさしつかえない範囲及び形状で設定される。これは、実際の計測によって得られる各voxe1ｉ内の吸収分布データは、各voxe1ｉ内の吸収分布の平均値として得られるからである。voxe1ｉの大きさ及び形状が決まることにより全voxe1の総数も決まる。すなわち、仮想媒体のvoxe1ｉの大きさ、形状、及び全voxe1の総数は、先に述べたのようにvoxe1ｉ内の吸収が一様であると見なせる限りにおいて任意である。実際のvoxe1ｉの大きさ及び形状は、使用される計測方式、或いはそれらの計測方式に要求されるＣＴ画像の解像度に応じて決定される。
【００５６】
本発明の散乱吸収体内部の光路分布計算方法においては、図１に示す仮想媒体の有限格子モデルの光入射位置ｐからインパルス光を入射した場合に、有限格子モデルの光検出位置ｑにおいて計測されるインパルス応答の中で、光路長がｌとなる光子集合を考える。ここで、光子入射位置ｐ、光子検出位置ｑは任意である。
【００５７】
図１に示す有限格子モデルに多数の光子からなるインパルス光を入射した場合のインパルス応答には、種々の光路長をもつ光子が含まれるが、観測時間ｔを決めるとｌ＝ｃｔとして光路長ｌが一意的に決まる。そして、光路長がｌとなる光子は多数存在する。そこで、この光路長がｌとなる光子集合の光路分布、すなわちvoxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）における光路長liを、光子が取り得る光路の確率分布で示すことと定義する。これに伴い有限格子モデル内のvoxe1ｉの光路分布は、voxe1ｉ内の時間分解型光路長liを要素とするN次元ベクトル［li］で表される。なお、先に述べたように光路長がｌとなる光子は多数存在するから、voxe1ｉ内の時間分解型光路長liは、これらの個々の光子のvoxe1ｉ内の光路長の平均（集合平均）である。また、このような系では、光入射位置ｐと光検出位置ｑを定めると、インパルス応答を構成する光子集合の時間分解型光路分布［li］が一意的に決まる。
【００５８】
ここで、上記の被計測媒体のインパルス応答を構成する光子集合の時間分解型光路分布［li］は被計測媒体の吸収に無関係であるので、吸収がないと仮定したときの仮想媒体のインパルス応答を構成する光子集合の時間分解型光路分布［li］に等しい。すなわち、仮想媒体のインパルス応答を構成する光子集合の時間分解型光路分布［li］を求めることにより、被計測媒体のインパルス応答を構成する光子集合の時間分解型光路分布［li］を求めることができる。
【００５９】
１. 解析モデルと解析の概要
まず、不均一散乱媒体である被計測媒体における減衰量、光路長などの関係を説明する。吸収と散乱は独立であるから、ここでは、図１に示す有限格子モデルが被計測媒体と同一の吸収分布をもつと仮定すれば良い。ただし、その他の条件は仮想媒体の有限格子モデルと同じである。
【００６０】
すると、まず上記の（i）ないし（v）の知見と、ＭＢＬから、光路長ｌと光路長分布liとの関係は、下記式で表現される。
【数３】
［式中、ｉは被計測媒体中のN個のvoxe1の番号（ｉ＝１，２，…，Ｎ；１≦Ｎ）を示し、liは、voxe1ｉ内の時間分解型光路長を示す。］
【００６１】
また、不均一散乱媒体を被計測媒体としたときの光路長がｌとなる光子集合のインパルス応答は、（１）式に対応する下記式で表現される。
【数４】
【数５】
【数６】
［式中、μaiはvoxe1ｉの吸収係数、s(t)は被計測媒体に吸収がないと仮定したときに検出位置ｑにおいて検出される被計測媒体のインパルス応答、h(t)は検出位置において検出される被計測媒体のインパルス応答を示す。また、（５）式の最後の辺（右辺）は時間分解型光路長li が吸収や吸収分布に依存しないことから導かれる。］
【００６２】
（４）式及び（５）式は、それぞれ（３）式で表現されるインパルス応答を微分形と積分形で表示したものである。（５）式より不均一散乱媒体のインパルス応答h(t)が吸収に依存しない項ln s(t)と下記式で表される吸収に依存する減衰項とに分離して記述できることがわかる。
【数７】
【００６３】
更に、（５）式より、この光応答の減衰量 ln[h(t)/s(t)] が被計測媒体内の時間分解型光路分布 [li] と吸収分布 [μai] とで記述できることがわかる。なお、ここで、被計測媒体内の時間分解型光路分布は、前記仮想媒体内、つまり吸収がないときの時間分解型光路分布に等しいことに注意する必要がある。従って、時間分解型の計測（ＴＲＳ）では、減衰量に対する重み関数は時間分解型光路長分布 [li] に等しく、この関係を利用した光ＣＴ画像再構が可能になる。なお、このような時間分解型光ＣＴは、画像再構成アルゴリズムが最も簡単になる。
【００６４】
より一般的な光ＣＴとしては、上記インパルス応答の時間積分値を利用するＴＩＳ、インパルス応答の所定時間範囲内の積分値を利用するＴＧＳ、正弦波変調光を利用するＰＭＳ、更には、これら各種の方法において、仮想媒体に基準となる均一な吸収を与え、その基準値からの偏差を計測するＡＶＭがる。ところが、これらの全ての方法に対して、減衰量あるいは減衰量偏差が吸収に依存する減衰項と吸収に依存しない項とに分離して記述できること、更に、各種の光応答の減衰量あるいは減衰量偏差が被計測媒体内の時間分解型光路分布と吸収分布で記述できることが導出され、この関係から減衰量あるいは減衰量偏差に対する重み関数を求めることができる。また、後述するように、減衰量あるいは減衰量偏差がキュムラント関数であることを利用すると、ＴＩＳやＰＭＳにおける重み関数が平均光路長や位相遅延およびキュムラントで記述できることが導出される。
【００６５】
例えば、ＴＩＳの場合における光応答は、前述のインパルス応答h(t)の時間積分に相当し、下記式のように表現できる。
【数８】
［式中、Ｉは不均一系のインパルス応答h(t)の時間積分を示す。］
【００６６】
また、（６）式の微分形及び積分形は、（３）式に対応する（４）式及び（５）式と同様に下記式のように表現できる。
【数９】
【数１０】
［（７）及び（８）式中、Liは検出位置において検出された時間積分Ｉを構成する光子集合のvoxelｉにおける平均光路長Li、つまり前記時間分解型光路長liの加重平均〈li〉を示わす。］
【００６７】
この平均光路長（時間分解型光路長liの加重平均）Liは、散乱体の散乱特性、境界条件、及び吸収や吸収分布に依存する。この際、（８）式の右辺第２項の積分は、（５）式のように積で表示することができないことに注意する必要がある。従って、減衰量に対する重み関数は、単なる平均光路長Liではなく、吸収に依存することになり、このことが光ＣＴ再構成アルゴリズムを複雑にしている。
【００６８】
また、仮想媒体に基準となる均一な吸収を与え、その基準値からの偏差として被計測媒体を計測するＡＶＭ計測では、voxe1ｉにおける被計測媒体と仮想媒体との減衰量差（減衰量偏差）(Ｂi は、吸収係数差（吸収係数偏差）(μai 、重み関数 Ｗi を用いて下記式のように表現される。
【数１１】
【数１２】
【数１３】
【数１４】
【数１５】
【数１６】
【数１７】
［（９）〜（１５）式中、μνは、仮想媒体の有する一定の吸収係数を示し、
Δμaiは、μνに対するμaiの吸収係数偏差を示し、
Ｂi(μai)は、被計測媒体のvoxe1ｉのインパルス応答の時間積分Iに関する減衰量を示し、
Ｂi(μν)は、仮想媒体のvoxe1ｉのインパルス応答の時間積分Iに関する減衰量を示し、
ΔＢi(μai)は、Ｂi(μν)に対するＢi(μai)の減衰量偏差を示し、
Ｌi(μν)は、仮想媒体のvoxe1ｉの平均光路長の差（時間積分Iに関する平均光路長の差）を示し、voxe1ｉ内の時間分解型光路長liの加重平均に等しい、
Ｗiは、減衰量偏差ΔＢi(μai)に対する重み関数であり、吸収係数偏差Δμai及びキュムラント(L, L',L",…)の関数をで表される。］
ここで、キュムラントは、仮想媒体に対して求めた時間分解波形から直接算出することができる。
【００６９】
また、上記の計測方法の他に各種の光ＣＴに応じて誘導される他の諸量、例えば、光路長平均の差と光路長分散に対応する重み関数との組み合わせを使用しても光路長平均の差と光路長分散に対応する重み関数との間には上記の（９）式と同様の関係式が成立し、吸収係数分布を定量計測することができる。具体的には波長の異なる光に対してこの計測を行い、得られた吸収係数と物質に固有の吸光係数から、吸収成分の濃度分布(絶対値)を定量計測する。
【００７０】
２．時間分解型光路長の計算方法
次に、有限格子モデルで記述される仮想媒体に対して、時間分解型光路分布［li］の要素である時間分解型光路長liを算出する式を導出する。
更に、この時間分解型光路長liから、ＴＩＳ（時間分解積分方式）、ＴＧＳ（時間分解ゲート積分方式）、及びＰＭＳ（位相変調方式）の光ＣＴに用いる計測方法に対応する各種の重み関数を計算することができるので、これらの具体例についても説明する。
【００７１】
被計測媒体のvoxe1ｉの有する吸収係数はμaiであるが、仮想媒体の有限格子モデルのvoxe1ｉの有する吸収係数は、μai=0である。このとき、図１に示す仮想媒体の有限格子モデルのインパルス応答は、（３）式にμai=0を代入して、下記式のように表現される。
【数１８】
すなわち、光検出位置ｑにおいて検出される仮想媒体のインパルス応答h(t)は、被計測媒体に吸収がないと仮定したときに光検出位置ｑにおいて検出される被計測媒体のインパルス応答s(t)に等しくなる。
【００７２】
このとき、時間分解型光路長liと光路長ｌとの関係は下記式で表現できる。
【数１９】
［式中、ｉは有限格子モデルを構成するvoxe1の番号ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ））を示す。］
【００７３】
上記の時間分解型光路長liと光路長ｌとの関係は、一個の光子がジグザグに伝播して同じvoxe1を2回以上通過した場合でも成立する。なお、上式においては、総光路長lと各voxe1ｉの時間分解型光路長liが時間ｔの関数であることを明示するため、それぞれを「ｌ(t)」、「li(t)」と表記した。
【００７４】
以下、このようなｔをパラメトリック変数として、voxe1ｉの時間分解型光路長liを算出する原理とその過程について詳細に説明する。
【００７５】
図２は、図１に示す有限格子モデルにおいて、光入射位置ｐからインパルス光を入射して伝播光を光検出位置ｑで検出する場合、インパルス応答ｓ(t)を構成し且つ光路長がｌとなる光子集合（すなわち、時間ｔの時点で検出される光子集合）のうち、voxe1ｉを通過した光子集合の経路と確率密度、およびそのような光子集合が光検出位置ｑにおいて時間ｔの時点で検出される値Ｕｐiｑ(t)とを示す模式図である。
【００７６】
ここで、上記有限格子モデルにおいて、光入射位置ｐからインパルス光を入射して伝播光を光検出位置ｑで検出する場合、インパルス応答ｓ(t)を構成し且つ光路長がｌとなる光子集合（すなわち、時間ｔの時点で検出される光子集合）が、光検出位置ｑにおいて検出される以前の時間t'(0≦t'≦t)の時点にvoxe1ｉに存在する光子存在確率密度Ｕi(t',t)を考える。なお、以下の説明において、このような光子存在確率密度Ｕi(t',t)を、単に「時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)」として記述する。
【００７７】
このとき、上記の新しい時間変数t'を導入して記述される時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)は、相反原理(reciprocity theorem)を適用することにより下記式のように表現できる。
【数２０】
［式中、Ｕpi(t')は、時間t'＝０において光入射位置ｐから入射した光子が時間t'＝t'においてvoxe1ｉ内に存在する光子存在確率密度を示し、
Ｕqi(t-t')は、時間t-t'＝０において光検出位置ｑから入射した光子が時間t-t'においてvoxe1ｉ内に存在する光子存在確率密度を示す。］
【００７８】
生体のような散乱媒体では、コリメート光入射やペンシルビーム入射に対して近似的に相反原理が成立するが、より厳密に相反原理を適用する場合には、等方光入射と呼ばれる光入射方法を利用することが好ましい。これは等方光入射法においては、相反原理がより精度よく成立するからである。なお、この等方光入射法は、発明者が既に報告した光入射方法であり、例えば、特願平５−３０１９７９、Y. Tsuchiya, K.Ohta, and T. Urakami: Jpn.J.Appl. Phys.34, Part 1, pp.2495-2501(1995)に開示されている。
【００７９】
また、以下の説明において、上記の（２．２）式で示されるＵpi(t')とＵqi(t-t')とを、それぞれ、単に「時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕpi(t')」、「時間t-t'におけるvoxe1ｉ内の光子存在確率密度Ｕqi(t-t')」として記述する。
【００８０】
図３は、上記の（２．２）式で示される有限格子モデル内の時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)、時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕpi(t')、及び、時間t-t'におけるvoxe1ｉ内の光子存在確率密度Ｕqi(t-t')の一例を示すグラフである。
【００８１】
上記のＵi(t',t)は、（２．２）式で示される二つの順問題を計算することにより算出することができる。この時、Ｕpi(t')は、有限格子モデルの光入射位置ｐからインパルス光を入射して計算することができる。また、Ｕqi(t-t')も有限格子モデルの光検出位置ｑからインパルス光を入射して計算することができる。なお、上記の順問題の具体的な計算には、モンテカルロ計算、経路積分（Path Integral）や輸送方程式の数値計算、光拡散方程式の数値計算などを利用することができる。この際、上記有限格子モデルでは吸収がないから、従来の吸収がある場合の順問題計算より計算アルゴリズムが簡単になり、計算時間も短くなる。
【００８２】
上記の時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)を算出することにより、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光入射位置ｐから光検出位置ｑへと有限格子モデル内を伝播する時間0≦t'≦tの間にvoxe1ｉに存在する確率密度Ｕi(t)を下記式に従って算出することができる。この確率密度Ｕi(t)には、有限格子モデル内（仮想媒体内、被計測媒体内）をジグザグに伝播する光子が異なる時間に同一のvoxelに存在した場合の確率が累積されている。
【数２１】
【００８３】
上記の（２．３） 式で示されるＵi(t)はインパルス応答ｓ(t)を構成し、且つ光路長がｌとなる光子集合のvoxe1ｉ内の光子存在累積確率を表しており、以下の説明において、単に「時間0≦t'≦tにおけるvoxe1ｉ内の光子存在累積確率Ｕi(t)」として記述する。
【００８４】
ここで、有限格子モデル内の全voxe1について（２．３）式で示される光子存在累積確率Ｕi(t)の総和をとり、有限格子モデル内の光子存在累積確率Ｕ(t)を求めると、次式が得られる。
【数２２】
【００８５】
この光子存在累積確率Ｕ(t)は、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光入射位置ｐから光検出位置ｑへと有限格子モデル内を伝播する時間0≦t'≦tの間に有限格子モデル内（全voxe1内）に存在する累積確率を表わす。
【００８６】
以下の説明において、上記の（２．４） 式で示されるＵ(t)を、「時間0≦t'≦tにおける全voxe1内の光子存在累積確率Ｕ(t)」として記述する。
【００８７】
一方、先に述べた時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)に対して、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光検出位置ｑにおいて検出される以前の時間t'(0≦t'≦t)の時点に有限格子モデル内（仮想媒体内）に存在する確率密度Ｕpq (t',t)を考える。すなわち、「時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)」に対して、「時間t'における全voxe1内の光子存在確率密度Ｕpq(t',t)」を考える。この場合、時間t'における光子存在確率密度を考えているから、全voxe1に対する光子存在確率密度の総和はＵpq(t',t)に等しく、このＵpq(t',t)はＵi(t',t)を用いて下記式で表現することができる。
【数２３】
【００８８】
この（２．５）式を時間t'= 0 から時間 t'= t の範囲で定積分したものは、当然、先に求めた時間 0≦t'≦tにおける全voxe1内の光子存在累積確率Ｕ(t)に等しく、下記式で表現することができる。
【数２４】
【００８９】
ここで、本発明の散乱吸収体内部の光路分布計算方法においては、吸収がない仮想媒体を考えているから、時間t'における全voxe1内の光子存在確率密度Ｕpq(t',t)は、時間t'＝0から時間t'＝ｔまで一定である。そこで、累積確率Ｕ(t)は、(２．６)式を用いて下記式のように表現することができる。
【数２５】
【００９０】
上記の事実に加えて、吸収がない仮想媒体において時間t'における全voxe1内の光子存在確率密度Ｕpq(t',t)は、時間t'＝0から時間t'＝ｔまで一定であることから、このＵpq(t',t)は、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答、つまり（１．１）式に求めたｓ(t)に等しいことがわかる。すなわち、Ｕpq(t',t)とｓ(t)との関係が、下記式のように表現される。
【数２６】
【００９１】
以上から、（２．４）式の時間0≦t'≦tにおける全voxe1内の光子存在累積確率Ｕ(t)は（２．８）式の結果を（２．７）式に代入することにより最終的に下記式で表現される。
【数２７】
【００９２】
更に、前出の図２で説明したＵｐiｑ(t)の全voxe1に対する総和は、有限格子モデル内の光子存在累積確率Ｕ(t)に等しく、下記式で示す関係が成り立つことが理解される。
【数２８】
【００９３】
なお、図４に上記の時間0≦t'≦tにおけるvoxe1ｉ内の光子存在累積確率Ｕi(t)、時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)、及び、時間t'における全voxe1内の光子存在確率密度Ｕpq(t',t)の関係を模式的に示す。
また、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光子集合の中で光路長がlとなる光子集合の時間t'におけるvoxeli内の時間分解型光路長をli(t',t)とすれば、このli(t',t)を時間t'= 0 から時間 t'= t まで積分したものが、時間分解型光路長li(t)に等しく、下記式で表現することができる。
【数２９】
【００９４】
ところで、有限格子モデル内の光速度は一定であるから、任意の時間あるいは任意の部分を考えた場合、光路長と格子存在確率は常に比例する。
【００９５】
具体的には、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光子集合の中で光路長がlとなる光子集合を考えた場合、任意の時間t'におけるvoxeli内の時間分解型光路長をli(t',t)は、時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)に比例する。そして、li(t',t)の時間積分値である時間分解型光路長li(t)もＵi(t',t)の時間積分値である光子存在累積確率Ｕi(t)に比例する。更には、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光子集合の中で光路長がlとなる光子集合を考えた場合、光路長lは光子存在累積確率Ｕ(t)に比例する。
【００９６】
この比例係数をＡとおけば、（２．１）式、（２．９）式、および（２．１１）式などの関係から下記式を得ることができる。
【数３０】
［式中、cは媒体中を伝播する光子の光速度である。］
【００９７】
（２．１２）式の関係を同値変形して時間分解型光路長li(t)について解くと下記式が得られる。
【数３１】
この（２．１３） 式により、時間分解型光路分布をＮ次元ベクトル[li] として表現する要素、つまり各voxelごとの時間分解型光路長liは、Ｕi(t',t) を（２．２）式に示す順問題を計算して求め、更に、このＵi(t',t)から計算されるＵi(t)、Ｕ(t)、および時間tを定めたときl = ctとして決まるl、或いはs(t)を用いて算出することができる。
【００９８】
なお、このような時間分解型光路長liは、時間分解応答を利用する時間分解型光ＣＴの画像再構成に直接的に用いられる。つまり、時間分解型光ＣＴにおける重み関数は、時間分解型光路長光路長liに等しいからである。また、後述するように、時間分解型光路長liは、各種の計測法式における重み関数を記述するための基本量となる。
【００９９】
３.時間分解型光路長liの具体的計算法
図５は、本発明の散乱吸収体内部の光路分布計算方法の一実施形態を示すフローチャートである。以下、図５に示すフローチャートについて説明する。
【０１００】
先ず、計測波長領域において、散乱吸収体である被計測媒体を、所定の大きさを有すると共に一様な吸収分布を有するとみなせるＮ（１≦Ｎ）個のvoxe1に分割する（Ｓ１）。
【０１０１】
次に、Ｎ個のvoxe1に分割された被計測媒体に対して、形状、境界条件及び散乱分布が被計測媒体と同一であり、屈折率分布が一様で、且つ吸収がないとみなせる仮想媒体を仮定する（Ｓ２）。
【０１０２】
次に、仮想媒体の表面に、インパルス光が入射する光入射位置ｐと、光入射位置ｐから入射して仮想媒体内を伝播したインパルス光の時間ｔにおけるインパルス応答ｓ(t)を検出する光検出位置ｑとを設定する（Ｓ３）。
【０１０３】
次に、仮想媒体内へ光入射位置ｐからインパルス光を入射させ、光検出位置において検出されるインパルス応答s(t)を算出する（Ｓ４）。
【０１０４】
次に、光入射位置ｐから仮想媒体内にインパルス光を入射し、時間t'＝０において光入射位置ｐから入射した光子が時間t'＝t'においてvoxe1ｉ内に存在する光子存在確率密度Ｕpi(t')を計算する。続いて、光検出位置ｑから仮想媒体内にインパルス光を入射し、時間t-t'＝０において光検出位置ｑから入射した光子が時間 t-t'においてvoxe1ｉ内に存在する光子存在確率密度Ｕqi(t-t')を計算する。そして、（２．２）式に従って、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置において時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光検出位置において検出される以前の時間t'(0≦t'≦t)の時点にvoxe1ｉに存在する光子存在確率密度Ｕi(t',t)を計算する（Ｓ５）。
【０１０５】
以上の場合、一個の光入射位置に対して、有限格子モデルを構成する全voxe1の光子存在確率密度を並列あるいは同時に計算することができる。そのため、各voxelの吸収係数を変化させて、voxe1の総数Ｎに等しい回数だけ順問題を繰り返して解くという時間のかかる従来型の計算が不要になり、計算時間が著しく短縮される。
【０１０６】
次に、（２．３）式に従って、光子存在確率密度Ｕi(t',t)から光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答を構成する光子集合がvoxe1ｉ内に存在していた光子存在累積確率Ｕi(t)を算出する（Ｓ６）。
【０１０７】
最後に（２．１３）式に従って光子存在累積確率Ｕi(t)及びインパルス応答s(t)を使用し、voxe1ｉの時間分解型光路長liを計算する（Ｓ８）。
【０１０８】
以上の本発明の第一実施形態において説明されるように、本発明の散乱吸収体内部の光路分布計算方法は、被計測媒体と同一の散乱特性と境界条件をもつが、吸収がないと仮定した仮想媒体に対して（２．２）式に示す順問題を解くので、従来の吸収があるときの順問題計算比べて、計算時間が短くなる。更に、各voxelの吸収係数を変化させてvoxe1の総数に等しい回数だけ順問題を繰り返して解くという時間のかかる従来の題計算が不要になり、voxe1ｉ内の時間分解型光路長liを直接且つ迅速に計算することができる。
【０１０９】
また、上記の散乱吸収体内部の光路分布計算方法を各種の光応答を用いた光ＣＴに適用する場合には、図５に示すフローチャートの後段において、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liを使用して、被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の計測方式に対応する重み関数Ｗiを算出する。この場合には、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liから所定の計測方式に対応する重み関数Ｗiを直接算出する方法と、はじめに仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liの時間に対する加重平均Liを求め、この加重平均Liを使用して所定の計測方式に対応する重み関数Ｗiを算出する方法とがある。
【０１１０】
先にも述べたように、従来は、時間分解型光路長liの具体的な計算方法が知られていなかったため、最初に時間分解型光路長liを求めて、その後で、時間分解型光路長liの加重平均Li及び重み関数Ｗiを算出するという発想がなかったのに対し、この重み関数の計算方法は、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて時間分解型光路長liの加重平均Li或いは被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の各種の具体的な計測方法に対応する重み関数Ｗiを迅速に計算することができる。なお、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liを使用した各種の光応答を用いた光ＣＴに対応する重み関数の具体的計算方法については後述する。
【０１１１】
図６は、本発明の散乱吸収体内部の光路分布計算方法の第二の実施形態を示すフローチャートである。以下、図６に示すフローチャートについて説明する。
【０１１２】
先ず、計測波長領域において、散乱吸収体である被計測媒体を、所定の大きさを有すると共に一様な吸収分布を有するとみなせるＮ（１≦Ｎ）個のvoxe1に分割する（Ｓ１）。
【０１１３】
次に、Ｎ個のvoxe1に分割された被計測媒体に対して、形状、境界条件及び散乱分布が被計測媒体と同一であり、屈折率分布が一様で、且つ吸収がないとみなせる仮想媒体を仮定する（Ｓ２）。
【０１１４】
次に、仮想媒体の表面に、インパルス光が入射する光入射位置ｐと、光入射位置ｐから入射して仮想媒体内を伝播したインパルス光の時間ｔにおけるインパルス応答ｓ(t)を検出する光検出位置ｑとを設定する（Ｓ３）。
【０１１５】
次に、光入射位置ｐから仮想媒体内にインパルス光を入射し、時間t'＝０において光入射位置ｐから入射した光子が時間t'＝t'においてvoxe1ｉ内に存在する光子存在確率密度Ｕpi(t')を計算する。続いて、光検出位置ｑから仮想媒体内にインパルス光を入射し、時間t-t'＝０において光検出位置ｑから入射した光子が時間t'-tにおいてvoxe1ｉ内に存在する光子存在確率密度Ｕqi(t-t')を計算する。そして、（２．２）式に従って、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置において時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光検出位置において検出される以前の時間t'(0≦t'≦t)の時点にvoxe1ｉに存在する光子存在確率密度Ｕi(t',t)を計算する（Ｓ５）。
【０１１６】
以上の場合にも、一個の光入射位置に対して、有限格子モデルを構成する全voxe1の光子存在確率密度を並列あるいは同時に計算することができる。そのため、各voxelの吸収係数を変化させて、voxe1の総数Ｎに等しい回数だけ順問題を繰り返して解くという時間のかかる従来型の計算が不要になり、計算時間が著しく短縮される。
【０１１７】
次に、（２．３）式に従って、光子存在確率密度Ｕi(t',t)から光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答を構成する光子集合がvoxe1ｉ内に存在していた光子存在累積確率Ｕi(t)を算出する（Ｓ６）。
【０１１８】
続いて、（２．４）式に従って、光子存在累積確率Ｕi(t)の全voxe1に対する総和Ｕ(t)を算出する（Ｓ７）。
【０１１９】
最後に（２．１３）式に従って光子存在累積確率Ｕi(t)及びインパルス応答Ｕ(t)を使用し、voxe1ｉの時間分解型光路長liを計算する（Ｓ１１）。
【０１２０】
以上の本発明の第二実施形態の説明から、先に述べた本発明の第一実施形態の説明と同様に、本発明の散乱吸収体内部の光路分布計算方法は、被計測媒体と同一の散乱特性と境界条件をもつが、吸収がないと仮定した仮想媒体に対して（２．２）式に示す順問題を解くので、従来の吸収があるときの順問題計算比べて、計算時間が短くなることがわかる。更に、各voxelの吸収係数を変化させてvoxe1の総数に等しい回数だけ順問題を繰り返して解くという時間のかかる従来の題計算が不要になり、voxe1ｉ内の時間分解型光路長liを直接且つ迅速に計算することができることがわかる。
【０１２１】
また、先に説明した本発明の第一実施形態と同様、上記の散乱吸収体内部の光路分布計算方法を各種の光応答を用いた光ＣＴに適用する場合には、図６に示すフローチャートの後段において、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liを使用して、被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の計測方式に対応する重み関数Ｗiを算出する。この場合にも、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liから所定の計測方式に対応する重み関数Ｗiを直接算出する方法と、はじめに仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liの時間に対する加重平均Liを求め、この加重平均Liを使用して所定の計測方式に対応する重み関数Ｗiを算出する方法とがある。
【０１２２】
この重み関数の計算方法は、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて時間分解型光路長liの加重平均Li或いは被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の各種の具体的な計測方法に対応する重み関数Ｗiを迅速に計算することができることになる。
【０１２３】
以上説明したように、本発明の散乱吸収体内部の光路分布計算方法は各種の光応答を用いた光ＣＴに適用することができる。この場合には、図５或いは図６に示すフローチャートの後段において、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liを使用して、被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の計測方式に対応する重み関数Ｗiを算出する。
【０１２４】
例えば、具体例としては発明者らがＭＢＬに基づいて導出した時間積分計測（ＴＩＳ）、及び時間分解ゲート積分計測（ＴＧＳ）、周波数領域での位相変調計測（ＰＭＳ）、更には平均値法（ＴＩＳ）等がある。これらの場合、各種の光応答は、先に述べたインパルス応答の数学変換で定義されるが、それぞれに対応する光ＣＴでは、それぞれに異なった重み関数を用いる。
【０１２５】
以下、時間積分計測（ＴＩＳ）、時間分解ゲート積分計測（ＴＧＳ）、及び周波数領域での位相変調計測（ＰＭＳ）等の光応答を用いた光ＣＴに本発明の散乱吸収体内部の光路分布計算方法により求めた時間分解型光路長liを適用する例について説明する。
【０１２６】
４．ＴＩＳにおける散乱吸収体内部の光路分布と重み関数
先に述べたように、インパルス応答h(t)の時間積分Ｉを用いるＴＩＳに基づく光ＣＴでは、インパルス応答h(t)に対するvoxel i の平均光路長Liや光路長分散などを用いて記述される重み関数を用いる。なお、この重み関数Ｗｉは、吸収分布が一様な基準媒体に対する重み関数として定義される。
【０１２７】
まず、吸収が不均一に分布する散乱媒体のインパルス応答h(t)の時間積分Iを考える。この際に考えるモデルは、図1に示したモデルのvoxel i の吸収係数をμaiとしたものである。これらを用いて、インパルス応答の時間積分Ｉは、下記式で与えられる。
【数３２】
［式中、μaiはvoxe1ｉの吸収係数を示し、
s(t)は、被計測媒体に吸収がないと仮定したときに検出位置において検出されるインパルス応答を示し、
Liは、検出位置において検出された時間積分Ｉを構成する光子集合のvoxel i における平均光路長（voxelｉにおける時間分解型光路長liの加重平均）を示す。］
この式では、ｌｎＩが吸収に無関係な項と、各voxelの吸収係数μaiに依存する項(減衰量)とに分離して記述されている。
【０１２８】
ここで、voxel i の平均光路長Liは、散乱体の散乱特性、境界条件、及び吸収や吸収分布に依存する量であり、下記式により与えられる。
【数３３】
【０１２９】
対応する光ＣＴに用いる重み関数Ｗｉは、均一な吸収係数μνをもつ仮想媒体に対して定義され、下記式により与えれる。この場合のモデルは、図1に示したモデルのvoxel i の吸収係数をμνとしたものである。
【数３４】
【数３５】
【数３６】
【数３７】
［（４．３）〜（４．６）式中、μνは、仮想媒体に与えた一定の吸収係数を示し、
Δμaiは、μνに対するμaiの吸収係数偏差を示し、
Ｌi(μν)は、一定の吸収係数μνをもつ仮想媒体のvoxe1ｉにおける平均光路長を示し（voxe1ｉ内の時間分解型光路長liの加重平均）、
Ｌ(μν)は、一定の吸収係数μνを与えた仮想媒体のインパルス応答の平均光路長を示し（仮想媒体内の光路長lの加重平均）、
Ｗiは、先に述べた（９）式中の減衰量偏差ΔＢi(μai)に対する重み関数であり、吸収係数偏差Δμai及びキュムラント(L, L',L",…)の関数で表される。］
ここで、キュムラントは、仮想媒体に対して求めた時間分解波形から直接算出することができる。
【０１３０】
５．ＴＧＳにおける散乱吸収体内部の光路分布と重み関数
まず、吸収が不均一に分布する散乱媒体のインパルス応答の時間分解ゲート積分信号Ｉgは、下記式で与えられる。
【数３８】
［式中、Ｌgiはインパルス応答の時間ゲート積分Ｉgを構成する光子集合のvoxel i における平均光路長（voxe1ｉ内のゲート時間内における時間分解型光路長liの加重平均）を示す。］
ここで、ゲート積分の時間範囲は[t1,t2]であり、これを[0,∞]にすると、前節に求めたインパルス応答h(t)の時間積分信号Iになる。
【０１３１】
対応する光ＣＴに用いる重み関数Ｗgiは、均一な吸収係数μνをもつ仮想媒体に対して定義され、下記式で与えられる。
【数３９】
【数４０】
【数４１】
［（５．２）〜（５．４）式中、μνは、仮想媒体に与えた一定の吸収係数を示し、
Δμaiは、前述の（４．４）式に示したように、μνに対するμaiの吸収係数偏差を示し、
Ｌgi(μν)は、仮想媒体のインパルス応答の時間ゲート積分Ｉgを構成する光子集合のvoxel i における平均光路長（voxe1ｉ内のゲート時間内における時間分解型光路長liの加重平均）を示し、
Ｌg(μν)は、一定の吸収係数μνを与えた仮想媒体のインパルス応答の時間ゲート積分Ｉgを構成する光子集合の平均光路長（加重平均）を示し
Ｗiは、先に述べた（９）式の減衰量偏差ΔＢi(μai)をＴＧＳに適用したときの減衰量偏差ΔＢgi(μai)に対する重み関数であり、吸収係数偏差Δμai及びキュムラント(Lg, Lg',Lg",…)の関数で表される。］
【０１３２】
６．ＰＭＳにおける散乱吸収体内部の光路分布と重み関数
吸収が不均一に分布する散乱媒体に対して、ＰＭＳで得られる応答は、吸収が不均一に分布する散乱媒体に対するインパルス応答h(t)をフーリエ変換したものであり下記式で与えられる。
【数４２】
【数４３】
ここで、Ｒ、Ｘ、Ａ、φはそれぞれ、
【数４４】
【数４５】
【数４６】
【数４７】
である。
［式（６．１）〜（６．６）中、Ｒ及びＸはそれぞれ実数部と虚数部を示し、Ａ及びφはそれぞれ振幅及び位相遅れを示す。］
Ｒ、Ｘ、Ａ及びφは、ロックインアンプなどで容易に計測することができる。
【０１３３】
また、ここで、
【数４８】
【数４９】
【数５０】
【数５１】
であるから、
【数５２】
【数５３】
【数５４】
【数５５】
が得られる。（６．１１）式〜（６．１４）式の各式の右辺第1項は、吸収がないと仮定したときの応答であり、第２項は吸収に依存する項である。
【０１３４】
以上のように、ＰＭＳでは（６．７）式〜（６．１０）式に示す4種類の記述方法があり、それぞれの式に対応する光ＣＴ画像再構成法が可能である。その際、先に述べたのと同様にして、（６．７）式〜（６．１０）式にμai=μνを代入して求めた値と、（６．１１）式〜（６．１４）式の関係を用いて、それぞれの重み関数を算出することができる。
【０１３５】
以上、本発明の好適な実施形態について詳細に説明したが、本発明は上記実施形態に限定されない。
【０１３６】
【発明の効果】
以上説明したように、本発明によれば、被計測媒体と同一の散乱特性と境界条件をもつが、吸収がないと仮定した仮想媒体に対して順問題を解いて、voxe1ｉ内の時間分解型光路長liを直接且つ迅速に計算することができるので、従来の吸収があるときの順問題計算比べて、計算時間が短くなる。更に、各voxelの吸収係数を変化させてvoxe1の総数に等しい回数だけ順問題を繰り返して解くという時間のかかる従来の順問題計算が不要になる。従って、生体などの不均一系の散乱吸収体内部の光路分布を迅速に算出することのできる方法を提供することが可能となる。
【図面の簡単な説明】
【図１】被計測媒体の内部光子移動を解析するための有限格子モデルを示す模式図である。
【図２】図１に示す有限格子モデルにおいて、光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答を構成する光子集合（光路長がlになる）のうち、voxe1ｉを通過した光子集合の経路とその確率を示す模式図である。
【図３】図２に示す有限格子モデル内の時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)、時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕpi(t')、及び、時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕqi(t-t')の時間波形の一例を示すグラフである。
【図４】時間0≦t'≦tにおけるvoxe1ｉ内の光子存在累積確率Ｕi(t)、 時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)、及び、時間t'における全voxe1内の光子存在確率密度Ｕpq(t',t)の関係を模式的に示すグラフである。
【図５】本発明の散乱吸収体内部の光路分布計算方法の一実施形態を示すフローチャートである。
【図６】本発明の散乱吸収体内部の光路分布計算方法の第二の実施形態を示すフローチャートである。
【符号の説明】
Ｓ１…第１ステップ、Ｓ２…第２ステップ、Ｓ３…第３ステップ、Ｓ４…第４ステップ、Ｓ５…第５ステップ、Ｓ６…第６ステップ、Ｓ７…第７ステップ、Ｓ８…第８ステップ、Ｓ１１…第１１ステップ。
【発明の属する技術分野】
本発明は、散乱吸収体内部の光路分布の計算方法に関する。より詳細には、光入射位置及び光検出位置を測定対象物の表面に沿って移動させて測定対象物の内部の情報を得る装置に適用可能な散乱吸収体内部の時間分解型光路分布の計算方法に関する。
【０００２】
この時間分解型光路分布から、適宜の計測方式における光路分布の計算、及び光減衰量や平均光路長に対する散乱吸収体各部の寄与度（寄与関数あるいは重み関数とも呼ばれる）を計算することができる。
【０００３】
【従来の技術】
従来、内部計測を含む広義の光ＣＴ等の分野においては、重み関数あるいは寄与関数を用いて測定対象物の内部情報を得る計測方法が知られており、これらに適用される重み関数などに関して既に多くの報告がなされている。これらは、例えば、以下の文献に記載されている。
【０００４】
(1) S. Arridge: SPIE Institutes for Advanced Optical Technologies, Vol. IS11, Medical Optical Tomography: Functional Imaging and Monitoring, 35-64(1993); (2) R. L. Barbour and H. L. Graber: ibid. 87-120(1993); (3) H. L. Graber, J. Chang, R. Aronson and R. L. Barbour: ibid. 121-143(1993); (4) J. C. Schotland, J. C. Haselgrove and J. S. Leigh: Applied Optics, 32, 448-5453(1993); (5) Chang, R. Aronson, H. L. Graber and R. L. Barbour: Proc. SPIE, 2389, 448-464(1995); (6) B. W. Pogue, M. S. Patterson, H. Jiang and K. D. Paulsen: Phys. Med. Biol. 40, 1709-1729(1995); (7) S. R. Arridge: Applied Optics, 34, 7395-7409(1995); (8) H. L. Graber, J. Chang, and R. L. Barbour: Proc. SPIE, 2570, 219-234(1995); (9) A. Maki and H. Koizumi: OSA TOPS, Vol.2, 299-304(1996); (10) H. Jiang, K. D. Paulsen and Ulf L. Osterberg: J. Opt. Soc. Am. A13, 253-266(1996); (11) S. R. Arridge and J. C. Hebden: Phys. Med. Biol. 42, 841-853(1997); (12) S. B. Colak, D. G. Papaioannou, G. W. 't Hooft, M. B. van der Mark, H. Schomberg, J. C. J. Paasschens, J. B. M. Melissen and N. A. A. J. van Astten: Applied Optics, 36, 180-213(1997)。
【０００５】
しかしながら、上記の文献（１）〜（１２）に開示されている光ＣＴ等の内部情報計測方法及びそれらに適用されている重み関数は以下のような問題を包含しているため、計測精度及び使い勝手に大きな問題が生じていた。また、実際に、上記分野において十分な性能をもつ光ＣＴの実用例は、まだ報告されていない。
【０００６】
すなわち、従来の光ＣＴ等における第１の問題は、媒体内の光子移動解析あるいは光子移動モデルが、輸送方程式に拡散近似を適用した光拡散方程式に基いているということにある。つまり、拡散近似は媒体中のフォトンの平均自由光路長に対して十分な大きさの媒体でのみ成立するから、比較的小さい媒体、内部の複雑な形状の組織、及び複雑な形状の媒体を取り扱うことができないという問題がある。また、拡散近似は等方散乱を前提としているので、非等方な散乱特性をもつ生体組織のような測定対象物の計測に適用すると、非等方散乱を等方散乱で近似したことに起因する無視できない誤差が発生するという問題もある。
【０００７】
更に、拡散方程式などの微分方程式は、解析的あるいは有限要素法など、いづれの数値計算手法を用いても、前もって境界条件（媒体形状や界面での反射特性など）を設定した後、解を求めなけらばならないという厄介な問題がある。つまり、通常、生体組織のような測定対象物では、境界条件が計測すべき場所や計測に用いる光の波長などによって変わるから、これらの影響を補正して精度を向上させようとすると、境界条件が変わるたびに複雑な計算をやり直す必要が生じ、結局、計算時間が極端に長くなるという大きな問題がある。
【０００８】
第２の問題は、測定対象物の内部情報を演算する際に、狭義の重み関数（寄与関数とも呼ばれる）、つまり媒体のインパルス応答を構成する光子集合に対する平均光路長（加重平均）或いはそれに相当する位相遅れ（周波数領域での計測）を適用していることにある。この場合、狭義の重み関数（平均光路長或いはそれに相当する位相遅れ）は、吸収係数や吸収分布に依存して変化するため、その取り扱いが極めて複雑になる。実際にこのような依存性を考慮した計算法を用いると、繰り返し回数の多い演算が不可避になるため、演算時間が極めて長くなり、実用的でないという問題が生じる。そこで、通常は、吸収係数や吸収分布に対する依存性を無視するが、このような近似に起因して誤差が増大するという重大な問題が生じる。
【０００９】
また、先に述べた重み関数の吸収係数や吸収分布に対する依存性に起因する誤差を低減するために、適宜の吸収があるときの重み関数を計算して用いるという方法がある。しかし、吸収があるときの重み関数の計算に要する時間は、吸収が無いときの重み関数の計算に要する時間に比較して、極端に長くなるという問題が発生する。
【００１０】
以上のようなことから、従来の狭義の重み関数（平均光路長或いはそれに相当する位相遅れ）を使用した計測方法は実用的でないと結論されている。
【００１１】
更に、上記の狭義の重み関数を適用する場合とは別に、輸送方程式の近似式や光拡散方程式に摂動理論を適用して、信号光と散乱吸収媒体の光学特性との関係を用いて測定対象物の内部情報を得る計測方法がある。しかしながら、この方法では、非線形効果（２次以上の項）の取り扱いが極めて複雑になる。この際、2次以上の項を含めた計算は、理論的にはコンピュータで計算することが可能であるが、現在の最速のコンピュータを用いても演算時間が膨大なものになり、実用化は不可能である。このため、常套的には２次以上の項を無視している。そのため、この方法では、複数個の比較的強い吸収領域が存在する媒体の光ＣＴ画像を再構成する際などで、吸収領域間の相互作用に起因する大きな誤差が発生するという問題が生じていた。
【００１２】
以上のように、従来の光ＣＴ等の測定対象物の内部情報を得る計測方法は、十分な精度の再構成画像が得られず、空間解像度、画像ひずみ、定量性、測定感度、所要計測時間などに大きな問題があるため、実用化が困難であった。
【００１３】
本発明者は、以上のような状況を打破するため、次のことが重要であると考えて一連の研究を推進してきた。すなわち、測定対象物の内部情報を得る計測方法、特に光ＣＴを実現させるために重要なことは、強い散乱媒体である生体組織の中を移動する光の振る舞いの詳細を明らかにし、検出される信号光と吸収成分を含む散乱媒体（散乱吸収体）の光学特性との関係をより厳密に記述し、信号光、並びに、信号光と散乱媒体の光学特性との関係を利用して光ＣＴ画像を再構成する新しいアルゴリズムを開発するということである。そして、本発明者は、散乱媒体中を移動する光の振る舞いを解析する際に、マイクロ・べア・ランバート則（Microscopic Beer-Lambert Law、以下「ＭＢＬ」という）を適用し、光ＣＴ画像の再構成に対して光路分布、つまり光子がどこを通過したかという情報を利用することを考えた。
【００１４】
そして、本発明者は、ＭＢＬに基づいて、散乱吸収体に対応する光子伝播モデル（有限格子モデル）を提案すると共に散乱吸収体の光学特性と信号光との関係を表す解析式を導出し、散乱吸収体内の光の振る舞いを解析する方法を開発してきた。これらの結果を本発明者らは、例えば以下の文献に掲載している。
【００１５】
(13) Y. Tsuchiya and T. Urakami: "Photon migration model for turbid biological medium having various shapes", Jpn. J. Appl. Phys. 34. Part2, pp. L79-81(1995); (14) Y. Tsuchiya and T. Urakami, "Frequency domain analysis of photon migration based on the microscopic Beer-Lambert law", Jpn. J. Appl. Phys. 35, Part1, pp.4848-4851(1996); (15) Y. Tsuchiya and T. Urakami: "Non-invasive spectroscopy of various]y shaped turbid media like human tissue based on the microscopic Beer-Lambert law", OSA TOPS, Biomedical Optical Spectroscopy and Diagnostics 1996, 3, pp.98- 100(1996); (16) Y. Tsuchiya and T. Urakami: "Quantitation of absorbing substances in turbid media such as human tissues based on the microscopic Beer-Lambert law". Optics Commun. 144, pp.269-280(1997); (17) Y. Tsuchiya and T. Urakami: "Optical quantitation of absorbers in variously shaped turbid media based on the microscopic Beer-Lambert law: A new approach to optical computerized tomography", Advances in Optical Biopsy and Optical Mammography (Annals of the New York Academy of Sciences), 838, pp.75-94(1998); (18) Y. Ueda, K. Ohta, M. Oda, M. Miwa, Y. Yamashita, and Y. Tsuchiya: "Average value method: A new approach to practical optical computed tomography for a turbid medium such as human tissue", Jpn. J. Appl. Phvs. 37, Part1. 5A, pp.2717-2723(1998); (19) 土屋 裕;「マイクロ・ベール・ランバート則と平均値法に基づく光ＣＴ画像の再構成」,O plus E, Vol.21, No.7, 814-821; (20) H. Zhang, M. Miwa, Y. Yamashita, and Y. Tsuchiya: "Quantitation of absorbers in turbid media using time-integrated spectroscopy based on microscopic Beer-Lambert law", Jpn. J. Appl. Phys. 37, Part1, pp.2724-2727(1998); (21) H. Zhang, Y. Tsuchiya, T. Urakami, M. Miwa, and Y. Yamashita: "Time integrated spectroscopy of turbid media based on the microscopic Beer-Lambert law: Consideration of the wavelength dependence of scattering properties", Optics Commun. 153, pp.314-322(1998); (22) Y. Tsuchiya, H. Zhang, T. Urakami, M. Miwa, and Y. Yamashita:"Time integrated spectroscopy of turbid media such as human tissues based on the microscopic Beer-Lambert law", Proc. JICAST'98/CPST'98, Joint international conference on Advanced science and technology, Hamamatsu, Aug.29-30, pp.237 -240(1998); (23) H. Zhang, T. Urakami, Y. Tsuchiya, Z. Liu, and T. Hiruma: "Time integrated spectroscopy of turbid media based on the microscopic Beer-Lambert law: Application to small-size phantoms having different boundary conditions", J. Biomedical Optics. 4, pp.183-190(1999);(24) Y. Tsuchiya, Y. Ueda, H. Zhang, Y. Yamashita, M. Oda, and T. Urakami: "Analytic expressions for determining the concentrations of absorber in turbid media by time-gating measurements", OSA TOPS, Advances in Optical Imaging and Photon Migration, 21, pp.67-72(1998); (25) H. Zhang, M. Miwa, T. Urakami, Y. Yamashita, and Y. Tsuchiya: "Simple subtraction method for determining the mean path length traveled by photons in turbid media", Jpn. J. Appl. Phys. 37, Part1, pp.700-704 (1998)。
【００１６】
このＭＢＬは「任意の散乱吸収分布をもつ媒体中において微視的にみたとき、吸収係数がμaである部分を伝播する光子は、伝播する長さｌのジグザグ光路に沿って吸収により指数関数的に減衰し、その光子の生存率は媒体の散乱特性や境界条件とは無関係な値exp(-μal)になり、吸収に起因する減衰量はμalになる」という法則であり、例えば下記式で表現される。
【数１】
【００１７】
この（１）式は、任意の散乱吸収分布をもつ媒体において、吸収と散乱が独立の事象であること、及び吸収に関して重ね合わせの原理が成り立つこと、すなわち、散乱吸収分布をもつ媒体が多成分系である場合において全体の吸光度はそれぞれの成分の吸光度の和で与えられることを示している。
【００１８】
本発明者らは上記（１）式で表現されるＭＢＬから、散乱媒体の各種光応答と光学特性との関係を表す式を導出し、散乱媒体の各種の応答が吸収に依存する減衰項と吸収に依存しない項とに分離して記述できることを明らかにした。このことから、多波長の光を用いてこの減衰項を分光計測すれば、吸収係数と、吸収成分に固有の吸光係数との関係から、散乱媒体内の吸収成分の濃度の絶対値を定量することができるようになった。
【００１９】
以上のようなＭＢＬに基づく計測方法は原理的に、媒体形状、境界条件及び散乱の影響を受けないという大きな特長があり、非等方散乱媒体や小さい媒体への適用も可能である。ただし、ここで計測される吸収係数や吸収成分の濃度は、吸収が均一な散乱媒体に対して正しい計測ができるが、吸収分布が均一でない不均一媒体に対する計測では光が通過した部分の平均吸収係数や平均濃度（光路分布に対する加重平均）となる。従って、高精度の光ＣＴを実現するためには、光路分布を知る必要が生じた。
【００２０】
そこで、本発明者らは、ＭＢＬに基づく散乱媒体内の光の振る舞いの解析方法をさらに発展させて、散乱や吸収が不均一に分布する不均一散乱媒体(不均一系)に対しても適用可能にして、媒体の形状や散乱特性などに影響されずに散乱媒体内における吸収成分の濃度分布を定量計測する種々の計測方法を開発してきた。
【００２１】
本発明者らがこれまでにＭＢＬに基づいて導出した具体的な計測方法としては、例えば、インパルス応答を利用する時間分解計測（ＴＲＳ）、インパルス応答の時間積分を利用する時間積分計測（ＴＩＳ）、インパルス応答の時間ゲート積分を利用する時間分解ゲート積分計測（ＴＧＳ）、周波数領域での位相変調計測（ＰＭＳ）、更には、平均値法に基づく光ＣＴ画像の再構成法（ＡＶＭ）などがある。これらの計測方法は、特願平１０−１４４３００、上記文献(7)、及び上記文献(15)〜(25)に記載されている。
【００２２】
上記のような新しい計測方法の出現によって計測精度は飛躍的に向上したが、光ＣＴ画像の再構成法では、依然として、従来型の重み関数、つまり平均光路長あるいはそれの相当する位相遅れが用いられていた。具体的には、有限格子モデルで記述した均一な吸収係数μνを有する仮想媒体のvoxe1ｉの吸収係数を少しだけ変化させ、その光出力をモンテカルロ計算、経路積分（Path Integral）や輸送方程式の数値計算、光拡散方程式の数値計算などを利用して計算し、変化前後の光出力の差、すなわち、均一な吸収係数μνを有する仮想媒体を基準として、光強度、光路長、減衰量などの偏差に対応するvoxe1ｉの重み関数（あるいは寄与関数）Ｗｉを求めている。
【００２３】
【発明が解決しようとする課題】
しかしながら、有限格子モデルを用いて記述した仮想媒体のvoxe1ｉの吸収係数を変化させることにより各種の計測方法に対応した重み関数（あるいは寄与関数）を求める上記従来の計測方法では、voxe1の総数に等しい回数の光拡散方程式等の順問題計算を行う必要があるので、この計算時間が極めて長くなり、結局、計測時間が非常に長くなるという問題があった。この場合、従来の狭義の重み関数（平均光路長或いはそれに相当する位相遅れ）を使用した計測方法では、計測精度が十分でないという問題もあった。更に、この計測精度を向上させようとすると、さらに計算時間が長くなり、結果として、さらに計測時間が長くなるという問題があった。
【００２４】
本発明は、上記従来の課題に鑑みて成されたものであり、生体などの不均一系の散乱吸収体内部の光路分布を精密、且つ迅速に算出することのできる方法を提供することを目的とする。なお、このようにして算出された光路分布から、散乱体内部の吸収成分の濃度分布を定量計測するための各種の具体的な計測方法に対応する重み関数Ｗｉを迅速に計算することができる。
【００２５】
【課題を解決するための手段】
本発明者は、上記目的を達成すべく鋭意、前記ＭＢＬに基づく解析法等の研究を進めた結果、Ｎ（１≦Ｎ）個のvoxelに分割された被計測媒体のvoxelｉ内の時間分解型光路長liが、この被計測媒体と同一の散乱特性と境界条件をもつが吸収がないと仮定した仮想媒体のvoxelｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）内の時間分解型光路長liと一致すること、この仮想媒体のvoxelｉ内の時間分解型光路長liは、仮想媒体のvoxelｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）内の光子存在確率等で記述することできること、更に、この仮想媒体のvoxelｉ内の時間分解型光路長liは、モンテカルロ計算、あるいは、経路積分、輸送方程式、または光拡散方程式の数値計算などによって仮想媒体内の光子移動を計算した結果を利用して、直接的に然も迅速に計算できることを見出し、本発明に到達した。
【００２６】
なお、従来のモンテカルロ計算は、有限格子モデル内のvoxelｉ内における散乱事象の発生回数を計数することを基本原理としており、従来の重み関数に対してもこの概念が踏襲・適用されている。この際、voxelｉ内の散乱特性が均一であれば、散乱回数と平均自由光路長との関係から光路長を推定することができるが、voxelｉ内の散乱特性が不均一であるときは、散乱回数と平均自由光路長との関係が散乱分布に依存して複雑になるから、本発明の目的である時間分解型光路長を精密に求めることができない。従って、本発明ではこのような従来の問題点を回避するため、媒体のvoxelｉ内の光子存在確率等からのvoxelｉ内の時間分解型光路長liを算出するという新しい手法を用いている。
【００２７】
すなわち、本発明の散乱吸収体内部の光路分布計算方法は、
計測波長領域において、散乱吸収体である被計測媒体を、所定の大きさを有すると共に一様な吸収分布を有するとみなせるＮ（１≦Ｎ）個のvoxe1に分割する第１ステップと、
前記Ｎ個のvoxe1に分割された前記被計測媒体に対して、形状、境界条件及び散乱分布が前記被計測媒体と同一であり、屈折率分布が一様で、且つ吸収がないとみなせる仮想媒体を仮定する第２ステップと、
前記仮想媒体の表面に、インパルス光が入射する光入射位置と、前記光入射位置から入射して前記仮想媒体内を伝播した前記インパルス光の時間ｔにおけるインパルス応答s(t)を検出する光検出位置とを設定する第３ステップと、
前記仮想媒体内へ前記光入射位置からインパルス光を入射する場合に、前記光検出位置において検出される前記インパルス応答s(t)を算出する第４ステップと、
前記光入射位置から前記仮想媒体内にインパルス光を入射する場合に、前記光検出位置において時間ｔの時点で検出される前記インパルス応答s(t)を構成する光子集合が、前記光検出位置において検出される以前の時間t'(0≦t'≦t)の時点で前記仮想媒体内の任意のvoxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）内に存在していた光子存在確率密度Ｕi(t',t)を算出する第５ステップと、
前記光子存在確率密度Ｕi(t',t)を使用して前記光検出位置において時間ｔの時点で検出される前記インパルス応答s(t)を構成する光子集合が前記voxe1ｉ内に存在していた光子存在累積確率Ｕi(t)を算出する第６ステップと、
前記光子存在累積確率Ｕi(t)及び前記インパルス応答s(t)を使用して前記voxe1ｉの時間分解型光路長liを算出する第８ステップと、
を備えることを特徴する方法である。
【００２８】
従来は、有限格子モデルで記述した仮想媒体のvoxeliの吸収係数を少し変化させて、吸収係数変化の前後の出力光強度の比、或いは減衰量の差等の計算値から各種の重み関数Ｗiを計算していた。この方法では、voxe1の総数Ｎに等しい回数の光拡散方程式等の順問題計算を行う必要があるので計算時間が長くなり、結局、計測時間が非常に長くなるという大きな問題があったが、上記ではこの問題が解決されている。更にその上、上記では吸収のない媒体について時間分解型光路長liを計算しているから、従来の吸収がある媒体に対する光拡散方程式等の順問題計算に長い時間を要したという問題も、解決されている。
【００２９】
また、本発明の方法では、一対の光入射位置と検出位置とからなる組に対して、全voxelのパラメータが並列に計算されるので、更に一層、計算時間を短縮することができる。つまり、本発明の散乱吸収体内部の光路分布計算方法は、時間分解型光路長liを直接精密、且つ迅速に計算することをはじめて見出したものである。
【００３０】
具体的にvoxe1ｉ内の時間分解型光路長liを求めるには、後述するように、有限格子モデル内に光入射位置からインパルス光を入射する場合に、光検出位置において時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光検出位置において検出される以前の時間t'(0≦t'≦t)の時点にvoxe1ｉに存在する光子存在確率密度Ｕi(t',t)を考える。
【００３１】
この新しい時間変数t'を導入して記述される時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)は、相反原理(reciprocity theorem)を適用することにより下記式のように表現できる。
【数２】
Ｕqi(t-t')は、時間t-t'＝０において光検出位置ｑから入射した光子が時間t-t'においてvoxe1ｉ内に存在する光子存在確率密度を示す。］
【００３２】
上記のＵi(t',t)は、（２．２）式で示されるように二つの順問題を計算することにより算出することができる。この時、Ｕpi(t')は、有限格子モデルの光入射位置ｐからインパルス光を入射して計算することができる。また、Ｕqi(t-t')も有限格子モデルの光検出位置ｑからインパルス光を入射して計算することができる。なお、（２．２）式で示される二つの順問題の具体的な計算方法については後に詳細に説明する。
【００３３】
更に本発明においては、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて時間分解型光路長liの時間に対する加重平均Ｌi を迅速に算出することができる。
【００３４】
従来は、時間分解型光路長liの具体的な計算方法が知られていなかったため、最初に時間分解型光路長liを求めて、その後で、時間分解型光路長liの加重平均Li及び重み関数Ｗiを算出するという発想がなかった。実際、時間分解型光路長liを直接計算して利用する計測方法の例はこれまでに報告されていない。
【００３５】
更に本発明においては、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の各種の具体的な計測方法に対応する重み関数Ｗiを迅速に計算することができる。
【００３６】
つまり、被計測媒体及びこれに対応する仮想媒体をＮ個のvoxelに分割した有限格子モデルを利用する場合において、ＭＢＬに基づくＴＲＳ、ＴＩＳ、ＴＧＳ、ＰＭＳ、ＡＶＭ等の具体的な計測方法において、被計測媒体内のvoxel i内の時間分解型光路長liと各種の計測法で計測される計測値および対応する重み関数との関係は、既に本発明者らのよって明らかにされている。
【００３７】
例えば、減衰量に対する重み関数は、ＴＲＳでは時間分解型光路長liそのものに等しく、ＴＩＳでは時間分解型光路長liの時間に対する加重平均Li（平均光路長Liともよばれる）や時間分解型光路長liの分散等を用いて記述される。より具体的には、減衰量と平均光路長を利用するＴＩＳやＡＶＭでは、被計測媒体内のvoxeli内の時間分解型光路長liと、その加重平均である平均光路長Li、更には時間分解型光路長liの分散等との関係を用いて記述される重み関数Ｗｉを用いるが、これらは後述するような（７）、（９）、（１２）式で表現されることが知られている。
【００３８】
ここで、本発明の散乱吸収体内部の光路分布計算方法における「時間分解型光路分布」とは、図1に示すように、３次元の被計測媒体をＮ（１≦Ｎ）個のvoxe1（体積素）に分割した有限格子モデルとし、この有限格子モデルの光入射位置ｐからインパルス光を入射した場合に、有限格子モデルの光検出位置ｑにおいて計測されるインパルス応答の光路長がｌとなる光子集合の光路分布を示す。具体的には、この時間分解型光路分布はvoxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）の時間分解型光路長liを要素とするＮ次元ベクトル［li］で表される。上記の定義から、この時間分解型光路分布は光路長ｌの関数、つまり時間tの関数である。各種の具体的な計測法における出力光信号の減衰量は、この時間分解型光路長liとvoxe1ｉの吸収係数の積の関数で表される。また、この時間分解型光路長liの時間に対する加重平均Liは、voxe1ｉの平均光路長とよばれ、これを要素とするＮ次元ベクトル［Li］は平均光路長分布とよばれる。
【００３９】
更に、重み関数は注目している物理量に対する吸収情報の寄与度として定義される。例えば、ＴＲＳで得られる時間分解型の減衰量に対する重み関数は、吸収や吸収分布に無関係となり、［時間分解型減衰量］＝［重み関数］・［吸収係数差］となり、この重み関数は上記で求めた時間分解型光路長liに等しい。また、インパルス応答の時間積分（ＴＲＳ）で得られる減衰量に対に重み関数は［ＴＲＳ減衰量］＝［重み関数］・［吸収係数差］で定義され、平均光路長差に対に対する重み関数は［平均光路長差］＝［重み関数］・［吸収係数差］で定義される。ＴＲＳの場合、これらの重み関数は吸収係数の関数になる。
【００４０】
このような重み関数を導入することによって、光ＣＴにおける画像再構成のための行列式が行列の積、つまり［重み関数］・［吸収係数差］として記述され、この式から［吸収係数差］すなわち吸収分布を計算することができる。
【００４１】
なお、上記の被計測媒体と仮想媒体との散乱特性は不均一であってもよく、順問題の具体的な計算には、モンテカルロ計算、経路積分（Path Integral）や輸送方程式の数値計算、光拡散方程式の数値計算などを利用することができる。
【００４２】
ここで、本発明の散乱吸収体内部の光路分布計算方法により求めた時間分解型光路長liが適用される「所定の計測方式」とは、特に限定されるものではなく、例えば、発明者らがＭＢＬに基づいて導出した時間分解計測（ＴＲＳ）、時間積分計測（ＴＩＳ）、時間分解ゲート積分計測（ＴＧＳ）、更には周波数領域での位相変調計測（ＰＭＳ）、平均値法に基づく光ＣＴ画像の再構成法（ＡＶＭ）等の散乱吸収体の内部計測方法でもよく、先に文献(1)〜(12)に示したようなＭＢＬに基づかない計測方法でもよい。
【００４３】
また、本発明の散乱吸収体内部の光路分布計算方法は、
計測波長領域において、散乱吸収体である被計測媒体を、所定の大きさを有すると共に一様な吸収分布を有するとみなせるＮ（１≦Ｎ）個のvoxe1に分割する第１ステップと、
前記Ｎ個のvoxe1に分割された前記被計測媒体に対して、形状、境界条件及び散乱分布が前記被計測媒体と同一であり、屈折率分布が一様で、吸収がないとみなせる仮想媒体を仮定する第２ステップと、
前記仮想媒体の表面に、インパルス光が入射する光入射位置と、前記光入射位置から入射して前記仮想媒体内を伝播した前記インパルス光の時間ｔにおけるインパルス応答ｓ(t)を検出する光検出位置とを設定する第３ステップと、
前記光入射位置から前記仮想媒体内にインパルス光を入射する場合に、前記光検出位置において時間ｔの時点で検出される前記インパルス応答ｓ(t)を構成する光子集合が、前記光検出位置において検出される以前の時間t'(0≦t'≦t)の時点で前記仮想媒体内の任意のvoxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）内に存在していた光子存在確率密度Ｕi(t',t)を算出する第５ステップと、
前記光子存在確率密度Ｕi(t',t)を使用して前記光検出位置において時間ｔの時点で検出される前記インパルス応答ｓ(t)を構成する光子集合が前記voxe1ｉ内に存在していた光子存在累積確率Ｕi(t)を算出する第６ステップと、
前記光子存在累積確率Ｕi(t)の全voxe1に対する総和Ｕ(t)を算出する第７ステップと、
前記光子存在累積確率Ｕi(t)及び前記光子存在累積確率Ｕi(t)の全voxe1に対する総和Ｕ(t)を使用して前記voxe1ｉの時間分解型光路長liを算出する第１１ステップと、
を備えることを特徴とする方法である。
【００４４】
上記の場合にも、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて時間分解型光路長liの時間に対する加重平均Ｌi を迅速に算出することができる。また、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の各種の具体的な計測方法に対応する重み関数Ｗiを迅速に計算することができる。
【００４５】
【発明の実施の形態】
以下、図面を参照しながら本発明による散乱吸収体内部光路分布計算方法の好適な実施形態について詳細に説明する。なお、以下の説明では、同一または相当部分には同一符号を付し、重複する説明は省略する。また、散乱吸収体内を散乱されつつ伝播する光子に関しては３次元座標を用いて考える必要があるが、以下では説明を簡単にするために場合により２次元座標を用いて説明する。
【００４６】
先ず、本発明の原理について説明する。
［本発明の原理］
先に述べたように、本発明者らがＭＢＬに基づいて開発してきた散乱媒体内の光の振る舞いを解析する解析法は、散乱や吸収が不均一に分布する不均一散乱媒体(不均一系)に適用することができる。そして、ＭＢＬに基づく解析法により、不均一散乱媒体の種々の光応答（光出力）が吸収に依存する減衰項と吸収に依存しない項とに分離して記述できること、更に、この光応答の減衰項が被計測媒体内の光路分布と吸収分布で記述できることが演繹される。従って、前もって被計測媒体内の光路分布を知ることができれば、被計測媒体内部の吸収成分の濃度分布を定量する光ＣＴなどが可能になる。つまり、光ＣＴを実現するには、被計測媒体内の光路分布を知ることと、光路分布を迅速に計算する方法を開発することが重要である。
【００４７】
本発明は、被計測媒体及びこれに対応する仮想媒体の有限格子モデルにおいて、以下の（i）〜（v）の事柄が成立することを出発点としている。
【００４８】
すなわち、（i）被計測媒体内の屈折率分布が一様であると仮定すれば、被計測媒体のインパルス応答（時間分解波形）では、光路長をｌとし、屈折率で決まる被計測媒体中の光速度をｃ、時間をｔとすれば、常にｌ＝ｃｔが成立する。ここで、被計測媒体内の屈折率分布が一様であると仮定するのは、被計測媒体として生体などの不均一散乱吸収体を選択した場合には、生体の主成分は水であるためその屈折率分布は一様であると見なしてよいからである。従って、以下の説明においては特にことわりがない限り屈折率分布が一様な被計測媒体を考える。
【００４９】
（ii）散乱と吸収が独立であることから、散乱と吸収が任意に分布する不均一媒体に対して光入射位置、及び光検出位置を定めると、インパルス応答を構成する(インパルス応答として検出される)光子集合の光路分布、つまり時間分解型光路分布が一意的に決まる。
【００５０】
（iii）上記の光路分布は、時間の関数となり、入射光子数（つまり光子密度）、被計測媒体の吸収や吸収分布に無関係であり、被計測媒体に吸収がないと仮定したときの光路分布に等しい。そして、３次元の被計測媒体を所定の大きさを有するＮ（１≦Ｎ）個のvoxe1（体積素）に分割した場合において、光路長がｌとなる光子集合の光路分布は、各voxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）の光路長liを要素とするＮ次元ベクトル［li］で表すことができる。
【００５１】
（iv）この時間分解型光路分布［li］は、被計測媒体と同一の散乱特性（不均一であってよい）と境界条件をもつが、屈折率分布が一様で、且つ吸収がない仮想媒体を仮想して計算することができる。
【００５２】
（v）各種の光ＣＴ画像再構成アルゴリズムに用いる各種の重み関数Ｗiを上記に定義した各voxe1ｉの時間分解型光路長liを用いて記述することができる。
なお、本発明者は、この時間分解型光路分布liを用いて記述される重み関数Ｗiと被計測媒体に対する実測値を用いて不均一媒体内の吸収成分の濃度分布を計測する光ＣＴを提案した（特願平１０−１４４３００、文献(19)）。
【００５３】
図１は被計測媒体である不均一散乱吸収体の内部光子移動を解析するための有限格子モデルを示す。
【００５４】
図１に示す有限格子モデルは、被計測媒体に対して形状、散乱分布、及び境界条件が同一であり、屈折率分布が一様、且つ吸収がないとして仮定される仮想媒体を示したものである。この際、３次元の被計測媒体をＮ（１≦Ｎ）個のvoxe1（体積素）に分割すると共にそれぞれのvoxe1に番号ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）を付けているので、被計測媒体に対応して仮定される仮想媒体も同様にＮ個のvoxe1を有している。
【００５５】
また、図１に示す仮想媒体の有限格子モデルのvoxe1ｉの大きさ及び形状は、対応する被計測媒体のvoxe1ｉ内の吸収分布が一様であると見なせるか、一様であると見なしてさしつかえない範囲及び形状で設定される。これは、実際の計測によって得られる各voxe1ｉ内の吸収分布データは、各voxe1ｉ内の吸収分布の平均値として得られるからである。voxe1ｉの大きさ及び形状が決まることにより全voxe1の総数も決まる。すなわち、仮想媒体のvoxe1ｉの大きさ、形状、及び全voxe1の総数は、先に述べたのようにvoxe1ｉ内の吸収が一様であると見なせる限りにおいて任意である。実際のvoxe1ｉの大きさ及び形状は、使用される計測方式、或いはそれらの計測方式に要求されるＣＴ画像の解像度に応じて決定される。
【００５６】
本発明の散乱吸収体内部の光路分布計算方法においては、図１に示す仮想媒体の有限格子モデルの光入射位置ｐからインパルス光を入射した場合に、有限格子モデルの光検出位置ｑにおいて計測されるインパルス応答の中で、光路長がｌとなる光子集合を考える。ここで、光子入射位置ｐ、光子検出位置ｑは任意である。
【００５７】
図１に示す有限格子モデルに多数の光子からなるインパルス光を入射した場合のインパルス応答には、種々の光路長をもつ光子が含まれるが、観測時間ｔを決めるとｌ＝ｃｔとして光路長ｌが一意的に決まる。そして、光路長がｌとなる光子は多数存在する。そこで、この光路長がｌとなる光子集合の光路分布、すなわちvoxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）における光路長liを、光子が取り得る光路の確率分布で示すことと定義する。これに伴い有限格子モデル内のvoxe1ｉの光路分布は、voxe1ｉ内の時間分解型光路長liを要素とするN次元ベクトル［li］で表される。なお、先に述べたように光路長がｌとなる光子は多数存在するから、voxe1ｉ内の時間分解型光路長liは、これらの個々の光子のvoxe1ｉ内の光路長の平均（集合平均）である。また、このような系では、光入射位置ｐと光検出位置ｑを定めると、インパルス応答を構成する光子集合の時間分解型光路分布［li］が一意的に決まる。
【００５８】
ここで、上記の被計測媒体のインパルス応答を構成する光子集合の時間分解型光路分布［li］は被計測媒体の吸収に無関係であるので、吸収がないと仮定したときの仮想媒体のインパルス応答を構成する光子集合の時間分解型光路分布［li］に等しい。すなわち、仮想媒体のインパルス応答を構成する光子集合の時間分解型光路分布［li］を求めることにより、被計測媒体のインパルス応答を構成する光子集合の時間分解型光路分布［li］を求めることができる。
【００５９】
１. 解析モデルと解析の概要
まず、不均一散乱媒体である被計測媒体における減衰量、光路長などの関係を説明する。吸収と散乱は独立であるから、ここでは、図１に示す有限格子モデルが被計測媒体と同一の吸収分布をもつと仮定すれば良い。ただし、その他の条件は仮想媒体の有限格子モデルと同じである。
【００６０】
すると、まず上記の（i）ないし（v）の知見と、ＭＢＬから、光路長ｌと光路長分布liとの関係は、下記式で表現される。
【数３】
【００６１】
また、不均一散乱媒体を被計測媒体としたときの光路長がｌとなる光子集合のインパルス応答は、（１）式に対応する下記式で表現される。
【数４】
【００６２】
（４）式及び（５）式は、それぞれ（３）式で表現されるインパルス応答を微分形と積分形で表示したものである。（５）式より不均一散乱媒体のインパルス応答h(t)が吸収に依存しない項ln s(t)と下記式で表される吸収に依存する減衰項とに分離して記述できることがわかる。
【数７】
更に、（５）式より、この光応答の減衰量 ln[h(t)/s(t)] が被計測媒体内の時間分解型光路分布 [li] と吸収分布 [μai] とで記述できることがわかる。なお、ここで、被計測媒体内の時間分解型光路分布は、前記仮想媒体内、つまり吸収がないときの時間分解型光路分布に等しいことに注意する必要がある。従って、時間分解型の計測（ＴＲＳ）では、減衰量に対する重み関数は時間分解型光路長分布 [li] に等しく、この関係を利用した光ＣＴ画像再構が可能になる。なお、このような時間分解型光ＣＴは、画像再構成アルゴリズムが最も簡単になる。
【００６４】
より一般的な光ＣＴとしては、上記インパルス応答の時間積分値を利用するＴＩＳ、インパルス応答の所定時間範囲内の積分値を利用するＴＧＳ、正弦波変調光を利用するＰＭＳ、更には、これら各種の方法において、仮想媒体に基準となる均一な吸収を与え、その基準値からの偏差を計測するＡＶＭがる。ところが、これらの全ての方法に対して、減衰量あるいは減衰量偏差が吸収に依存する減衰項と吸収に依存しない項とに分離して記述できること、更に、各種の光応答の減衰量あるいは減衰量偏差が被計測媒体内の時間分解型光路分布と吸収分布で記述できることが導出され、この関係から減衰量あるいは減衰量偏差に対する重み関数を求めることができる。また、後述するように、減衰量あるいは減衰量偏差がキュムラント関数であることを利用すると、ＴＩＳやＰＭＳにおける重み関数が平均光路長や位相遅延およびキュムラントで記述できることが導出される。
【００６５】
例えば、ＴＩＳの場合における光応答は、前述のインパルス応答h(t)の時間積分に相当し、下記式のように表現できる。
【数８】
【００６６】
また、（６）式の微分形及び積分形は、（３）式に対応する（４）式及び（５）式と同様に下記式のように表現できる。
【数９】
【００６７】
この平均光路長（時間分解型光路長liの加重平均）Liは、散乱体の散乱特性、境界条件、及び吸収や吸収分布に依存する。この際、（８）式の右辺第２項の積分は、（５）式のように積で表示することができないことに注意する必要がある。従って、減衰量に対する重み関数は、単なる平均光路長Liではなく、吸収に依存することになり、このことが光ＣＴ再構成アルゴリズムを複雑にしている。
【００６８】
また、仮想媒体に基準となる均一な吸収を与え、その基準値からの偏差として被計測媒体を計測するＡＶＭ計測では、voxe1ｉにおける被計測媒体と仮想媒体との減衰量差（減衰量偏差）(Ｂi は、吸収係数差（吸収係数偏差）(μai 、重み関数 Ｗi を用いて下記式のように表現される。
【数１１】
Δμaiは、μνに対するμaiの吸収係数偏差を示し、
Ｂi(μai)は、被計測媒体のvoxe1ｉのインパルス応答の時間積分Iに関する減衰量を示し、
Ｂi(μν)は、仮想媒体のvoxe1ｉのインパルス応答の時間積分Iに関する減衰量を示し、
ΔＢi(μai)は、Ｂi(μν)に対するＢi(μai)の減衰量偏差を示し、
Ｌi(μν)は、仮想媒体のvoxe1ｉの平均光路長の差（時間積分Iに関する平均光路長の差）を示し、voxe1ｉ内の時間分解型光路長liの加重平均に等しい、
Ｗiは、減衰量偏差ΔＢi(μai)に対する重み関数であり、吸収係数偏差Δμai及びキュムラント(L, L',L",…)の関数をで表される。］
ここで、キュムラントは、仮想媒体に対して求めた時間分解波形から直接算出することができる。
【００６９】
また、上記の計測方法の他に各種の光ＣＴに応じて誘導される他の諸量、例えば、光路長平均の差と光路長分散に対応する重み関数との組み合わせを使用しても光路長平均の差と光路長分散に対応する重み関数との間には上記の（９）式と同様の関係式が成立し、吸収係数分布を定量計測することができる。具体的には波長の異なる光に対してこの計測を行い、得られた吸収係数と物質に固有の吸光係数から、吸収成分の濃度分布(絶対値)を定量計測する。
【００７０】
２．時間分解型光路長の計算方法
次に、有限格子モデルで記述される仮想媒体に対して、時間分解型光路分布［li］の要素である時間分解型光路長liを算出する式を導出する。
更に、この時間分解型光路長liから、ＴＩＳ（時間分解積分方式）、ＴＧＳ（時間分解ゲート積分方式）、及びＰＭＳ（位相変調方式）の光ＣＴに用いる計測方法に対応する各種の重み関数を計算することができるので、これらの具体例についても説明する。
【００７１】
被計測媒体のvoxe1ｉの有する吸収係数はμaiであるが、仮想媒体の有限格子モデルのvoxe1ｉの有する吸収係数は、μai=0である。このとき、図１に示す仮想媒体の有限格子モデルのインパルス応答は、（３）式にμai=0を代入して、下記式のように表現される。
【数１８】
【００７２】
このとき、時間分解型光路長liと光路長ｌとの関係は下記式で表現できる。
【数１９】
【００７３】
上記の時間分解型光路長liと光路長ｌとの関係は、一個の光子がジグザグに伝播して同じvoxe1を2回以上通過した場合でも成立する。なお、上式においては、総光路長lと各voxe1ｉの時間分解型光路長liが時間ｔの関数であることを明示するため、それぞれを「ｌ(t)」、「li(t)」と表記した。
【００７４】
以下、このようなｔをパラメトリック変数として、voxe1ｉの時間分解型光路長liを算出する原理とその過程について詳細に説明する。
【００７５】
図２は、図１に示す有限格子モデルにおいて、光入射位置ｐからインパルス光を入射して伝播光を光検出位置ｑで検出する場合、インパルス応答ｓ(t)を構成し且つ光路長がｌとなる光子集合（すなわち、時間ｔの時点で検出される光子集合）のうち、voxe1ｉを通過した光子集合の経路と確率密度、およびそのような光子集合が光検出位置ｑにおいて時間ｔの時点で検出される値Ｕｐiｑ(t)とを示す模式図である。
【００７６】
ここで、上記有限格子モデルにおいて、光入射位置ｐからインパルス光を入射して伝播光を光検出位置ｑで検出する場合、インパルス応答ｓ(t)を構成し且つ光路長がｌとなる光子集合（すなわち、時間ｔの時点で検出される光子集合）が、光検出位置ｑにおいて検出される以前の時間t'(0≦t'≦t)の時点にvoxe1ｉに存在する光子存在確率密度Ｕi(t',t)を考える。なお、以下の説明において、このような光子存在確率密度Ｕi(t',t)を、単に「時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)」として記述する。
【００７７】
このとき、上記の新しい時間変数t'を導入して記述される時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)は、相反原理(reciprocity theorem)を適用することにより下記式のように表現できる。
【数２０】
Ｕqi(t-t')は、時間t-t'＝０において光検出位置ｑから入射した光子が時間t-t'においてvoxe1ｉ内に存在する光子存在確率密度を示す。］
【００７８】
生体のような散乱媒体では、コリメート光入射やペンシルビーム入射に対して近似的に相反原理が成立するが、より厳密に相反原理を適用する場合には、等方光入射と呼ばれる光入射方法を利用することが好ましい。これは等方光入射法においては、相反原理がより精度よく成立するからである。なお、この等方光入射法は、発明者が既に報告した光入射方法であり、例えば、特願平５−３０１９７９、Y. Tsuchiya, K.Ohta, and T. Urakami: Jpn.J.Appl. Phys.34, Part 1, pp.2495-2501(1995)に開示されている。
【００７９】
また、以下の説明において、上記の（２．２）式で示されるＵpi(t')とＵqi(t-t')とを、それぞれ、単に「時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕpi(t')」、「時間t-t'におけるvoxe1ｉ内の光子存在確率密度Ｕqi(t-t')」として記述する。
【００８０】
図３は、上記の（２．２）式で示される有限格子モデル内の時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)、時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕpi(t')、及び、時間t-t'におけるvoxe1ｉ内の光子存在確率密度Ｕqi(t-t')の一例を示すグラフである。
【００８１】
上記のＵi(t',t)は、（２．２）式で示される二つの順問題を計算することにより算出することができる。この時、Ｕpi(t')は、有限格子モデルの光入射位置ｐからインパルス光を入射して計算することができる。また、Ｕqi(t-t')も有限格子モデルの光検出位置ｑからインパルス光を入射して計算することができる。なお、上記の順問題の具体的な計算には、モンテカルロ計算、経路積分（Path Integral）や輸送方程式の数値計算、光拡散方程式の数値計算などを利用することができる。この際、上記有限格子モデルでは吸収がないから、従来の吸収がある場合の順問題計算より計算アルゴリズムが簡単になり、計算時間も短くなる。
【００８２】
上記の時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)を算出することにより、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光入射位置ｐから光検出位置ｑへと有限格子モデル内を伝播する時間0≦t'≦tの間にvoxe1ｉに存在する確率密度Ｕi(t)を下記式に従って算出することができる。この確率密度Ｕi(t)には、有限格子モデル内（仮想媒体内、被計測媒体内）をジグザグに伝播する光子が異なる時間に同一のvoxelに存在した場合の確率が累積されている。
【数２１】
上記の（２．３） 式で示されるＵi(t)はインパルス応答ｓ(t)を構成し、且つ光路長がｌとなる光子集合のvoxe1ｉ内の光子存在累積確率を表しており、以下の説明において、単に「時間0≦t'≦tにおけるvoxe1ｉ内の光子存在累積確率Ｕi(t)」として記述する。
【００８４】
ここで、有限格子モデル内の全voxe1について（２．３）式で示される光子存在累積確率Ｕi(t)の総和をとり、有限格子モデル内の光子存在累積確率Ｕ(t)を求めると、次式が得られる。
【数２２】
この光子存在累積確率Ｕ(t)は、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光入射位置ｐから光検出位置ｑへと有限格子モデル内を伝播する時間0≦t'≦tの間に有限格子モデル内（全voxe1内）に存在する累積確率を表わす。
【００８６】
以下の説明において、上記の（２．４） 式で示されるＵ(t)を、「時間0≦t'≦tにおける全voxe1内の光子存在累積確率Ｕ(t)」として記述する。
【００８７】
一方、先に述べた時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)に対して、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光検出位置ｑにおいて検出される以前の時間t'(0≦t'≦t)の時点に有限格子モデル内（仮想媒体内）に存在する確率密度Ｕpq (t',t)を考える。すなわち、「時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)」に対して、「時間t'における全voxe1内の光子存在確率密度Ｕpq(t',t)」を考える。この場合、時間t'における光子存在確率密度を考えているから、全voxe1に対する光子存在確率密度の総和はＵpq(t',t)に等しく、このＵpq(t',t)はＵi(t',t)を用いて下記式で表現することができる。
【数２３】
この（２．５）式を時間t'= 0 から時間 t'= t の範囲で定積分したものは、当然、先に求めた時間 0≦t'≦tにおける全voxe1内の光子存在累積確率Ｕ(t)に等しく、下記式で表現することができる。
【数２４】
ここで、本発明の散乱吸収体内部の光路分布計算方法においては、吸収がない仮想媒体を考えているから、時間t'における全voxe1内の光子存在確率密度Ｕpq(t',t)は、時間t'＝0から時間t'＝ｔまで一定である。そこで、累積確率Ｕ(t)は、(２．６)式を用いて下記式のように表現することができる。
【数２５】
上記の事実に加えて、吸収がない仮想媒体において時間t'における全voxe1内の光子存在確率密度Ｕpq(t',t)は、時間t'＝0から時間t'＝ｔまで一定であることから、このＵpq(t',t)は、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答、つまり（１．１）式に求めたｓ(t)に等しいことがわかる。すなわち、Ｕpq(t',t)とｓ(t)との関係が、下記式のように表現される。
【数２６】
以上から、（２．４）式の時間0≦t'≦tにおける全voxe1内の光子存在累積確率Ｕ(t)は（２．８）式の結果を（２．７）式に代入することにより最終的に下記式で表現される。
【数２７】
更に、前出の図２で説明したＵｐiｑ(t)の全voxe1に対する総和は、有限格子モデル内の光子存在累積確率Ｕ(t)に等しく、下記式で示す関係が成り立つことが理解される。
【数２８】
なお、図４に上記の時間0≦t'≦tにおけるvoxe1ｉ内の光子存在累積確率Ｕi(t)、時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)、及び、時間t'における全voxe1内の光子存在確率密度Ｕpq(t',t)の関係を模式的に示す。
また、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光子集合の中で光路長がlとなる光子集合の時間t'におけるvoxeli内の時間分解型光路長をli(t',t)とすれば、このli(t',t)を時間t'= 0 から時間 t'= t まで積分したものが、時間分解型光路長li(t)に等しく、下記式で表現することができる。
【数２９】
ところで、有限格子モデル内の光速度は一定であるから、任意の時間あるいは任意の部分を考えた場合、光路長と格子存在確率は常に比例する。
【００９５】
具体的には、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光子集合の中で光路長がlとなる光子集合を考えた場合、任意の時間t'におけるvoxeli内の時間分解型光路長をli(t',t)は、時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)に比例する。そして、li(t',t)の時間積分値である時間分解型光路長li(t)もＵi(t',t)の時間積分値である光子存在累積確率Ｕi(t)に比例する。更には、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光子集合の中で光路長がlとなる光子集合を考えた場合、光路長lは光子存在累積確率Ｕ(t)に比例する。
【００９６】
この比例係数をＡとおけば、（２．１）式、（２．９）式、および（２．１１）式などの関係から下記式を得ることができる。
【数３０】
【００９７】
（２．１２）式の関係を同値変形して時間分解型光路長li(t)について解くと下記式が得られる。
【数３１】
【００９８】
なお、このような時間分解型光路長liは、時間分解応答を利用する時間分解型光ＣＴの画像再構成に直接的に用いられる。つまり、時間分解型光ＣＴにおける重み関数は、時間分解型光路長光路長liに等しいからである。また、後述するように、時間分解型光路長liは、各種の計測法式における重み関数を記述するための基本量となる。
【００９９】
３.時間分解型光路長liの具体的計算法
図５は、本発明の散乱吸収体内部の光路分布計算方法の一実施形態を示すフローチャートである。以下、図５に示すフローチャートについて説明する。
【０１００】
先ず、計測波長領域において、散乱吸収体である被計測媒体を、所定の大きさを有すると共に一様な吸収分布を有するとみなせるＮ（１≦Ｎ）個のvoxe1に分割する（Ｓ１）。
【０１０１】
次に、Ｎ個のvoxe1に分割された被計測媒体に対して、形状、境界条件及び散乱分布が被計測媒体と同一であり、屈折率分布が一様で、且つ吸収がないとみなせる仮想媒体を仮定する（Ｓ２）。
【０１０２】
次に、仮想媒体の表面に、インパルス光が入射する光入射位置ｐと、光入射位置ｐから入射して仮想媒体内を伝播したインパルス光の時間ｔにおけるインパルス応答ｓ(t)を検出する光検出位置ｑとを設定する（Ｓ３）。
【０１０３】
次に、仮想媒体内へ光入射位置ｐからインパルス光を入射させ、光検出位置において検出されるインパルス応答s(t)を算出する（Ｓ４）。
【０１０４】
次に、光入射位置ｐから仮想媒体内にインパルス光を入射し、時間t'＝０において光入射位置ｐから入射した光子が時間t'＝t'においてvoxe1ｉ内に存在する光子存在確率密度Ｕpi(t')を計算する。続いて、光検出位置ｑから仮想媒体内にインパルス光を入射し、時間t-t'＝０において光検出位置ｑから入射した光子が時間 t-t'においてvoxe1ｉ内に存在する光子存在確率密度Ｕqi(t-t')を計算する。そして、（２．２）式に従って、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置において時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光検出位置において検出される以前の時間t'(0≦t'≦t)の時点にvoxe1ｉに存在する光子存在確率密度Ｕi(t',t)を計算する（Ｓ５）。
【０１０５】
以上の場合、一個の光入射位置に対して、有限格子モデルを構成する全voxe1の光子存在確率密度を並列あるいは同時に計算することができる。そのため、各voxelの吸収係数を変化させて、voxe1の総数Ｎに等しい回数だけ順問題を繰り返して解くという時間のかかる従来型の計算が不要になり、計算時間が著しく短縮される。
【０１０６】
次に、（２．３）式に従って、光子存在確率密度Ｕi(t',t)から光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答を構成する光子集合がvoxe1ｉ内に存在していた光子存在累積確率Ｕi(t)を算出する（Ｓ６）。
【０１０７】
最後に（２．１３）式に従って光子存在累積確率Ｕi(t)及びインパルス応答s(t)を使用し、voxe1ｉの時間分解型光路長liを計算する（Ｓ８）。
【０１０８】
以上の本発明の第一実施形態において説明されるように、本発明の散乱吸収体内部の光路分布計算方法は、被計測媒体と同一の散乱特性と境界条件をもつが、吸収がないと仮定した仮想媒体に対して（２．２）式に示す順問題を解くので、従来の吸収があるときの順問題計算比べて、計算時間が短くなる。更に、各voxelの吸収係数を変化させてvoxe1の総数に等しい回数だけ順問題を繰り返して解くという時間のかかる従来の題計算が不要になり、voxe1ｉ内の時間分解型光路長liを直接且つ迅速に計算することができる。
【０１０９】
また、上記の散乱吸収体内部の光路分布計算方法を各種の光応答を用いた光ＣＴに適用する場合には、図５に示すフローチャートの後段において、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liを使用して、被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の計測方式に対応する重み関数Ｗiを算出する。この場合には、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liから所定の計測方式に対応する重み関数Ｗiを直接算出する方法と、はじめに仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liの時間に対する加重平均Liを求め、この加重平均Liを使用して所定の計測方式に対応する重み関数Ｗiを算出する方法とがある。
【０１１０】
先にも述べたように、従来は、時間分解型光路長liの具体的な計算方法が知られていなかったため、最初に時間分解型光路長liを求めて、その後で、時間分解型光路長liの加重平均Li及び重み関数Ｗiを算出するという発想がなかったのに対し、この重み関数の計算方法は、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて時間分解型光路長liの加重平均Li或いは被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の各種の具体的な計測方法に対応する重み関数Ｗiを迅速に計算することができる。なお、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liを使用した各種の光応答を用いた光ＣＴに対応する重み関数の具体的計算方法については後述する。
【０１１１】
図６は、本発明の散乱吸収体内部の光路分布計算方法の第二の実施形態を示すフローチャートである。以下、図６に示すフローチャートについて説明する。
【０１１２】
先ず、計測波長領域において、散乱吸収体である被計測媒体を、所定の大きさを有すると共に一様な吸収分布を有するとみなせるＮ（１≦Ｎ）個のvoxe1に分割する（Ｓ１）。
【０１１３】
次に、Ｎ個のvoxe1に分割された被計測媒体に対して、形状、境界条件及び散乱分布が被計測媒体と同一であり、屈折率分布が一様で、且つ吸収がないとみなせる仮想媒体を仮定する（Ｓ２）。
【０１１４】
次に、仮想媒体の表面に、インパルス光が入射する光入射位置ｐと、光入射位置ｐから入射して仮想媒体内を伝播したインパルス光の時間ｔにおけるインパルス応答ｓ(t)を検出する光検出位置ｑとを設定する（Ｓ３）。
【０１１５】
次に、光入射位置ｐから仮想媒体内にインパルス光を入射し、時間t'＝０において光入射位置ｐから入射した光子が時間t'＝t'においてvoxe1ｉ内に存在する光子存在確率密度Ｕpi(t')を計算する。続いて、光検出位置ｑから仮想媒体内にインパルス光を入射し、時間t-t'＝０において光検出位置ｑから入射した光子が時間t'-tにおいてvoxe1ｉ内に存在する光子存在確率密度Ｕqi(t-t')を計算する。そして、（２．２）式に従って、有限格子モデル内に光入射位置ｐからインパルス光を入射する場合に、光検出位置において時間ｔの時点で検出されるインパルス応答ｓ(t)を構成する光路長がｌとなる光子集合が、光検出位置において検出される以前の時間t'(0≦t'≦t)の時点にvoxe1ｉに存在する光子存在確率密度Ｕi(t',t)を計算する（Ｓ５）。
【０１１６】
以上の場合にも、一個の光入射位置に対して、有限格子モデルを構成する全voxe1の光子存在確率密度を並列あるいは同時に計算することができる。そのため、各voxelの吸収係数を変化させて、voxe1の総数Ｎに等しい回数だけ順問題を繰り返して解くという時間のかかる従来型の計算が不要になり、計算時間が著しく短縮される。
【０１１７】
次に、（２．３）式に従って、光子存在確率密度Ｕi(t',t)から光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答を構成する光子集合がvoxe1ｉ内に存在していた光子存在累積確率Ｕi(t)を算出する（Ｓ６）。
【０１１８】
続いて、（２．４）式に従って、光子存在累積確率Ｕi(t)の全voxe1に対する総和Ｕ(t)を算出する（Ｓ７）。
【０１１９】
最後に（２．１３）式に従って光子存在累積確率Ｕi(t)及びインパルス応答Ｕ(t)を使用し、voxe1ｉの時間分解型光路長liを計算する（Ｓ１１）。
【０１２０】
以上の本発明の第二実施形態の説明から、先に述べた本発明の第一実施形態の説明と同様に、本発明の散乱吸収体内部の光路分布計算方法は、被計測媒体と同一の散乱特性と境界条件をもつが、吸収がないと仮定した仮想媒体に対して（２．２）式に示す順問題を解くので、従来の吸収があるときの順問題計算比べて、計算時間が短くなることがわかる。更に、各voxelの吸収係数を変化させてvoxe1の総数に等しい回数だけ順問題を繰り返して解くという時間のかかる従来の題計算が不要になり、voxe1ｉ内の時間分解型光路長liを直接且つ迅速に計算することができることがわかる。
【０１２１】
また、先に説明した本発明の第一実施形態と同様、上記の散乱吸収体内部の光路分布計算方法を各種の光応答を用いた光ＣＴに適用する場合には、図６に示すフローチャートの後段において、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liを使用して、被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の計測方式に対応する重み関数Ｗiを算出する。この場合にも、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liから所定の計測方式に対応する重み関数Ｗiを直接算出する方法と、はじめに仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liの時間に対する加重平均Liを求め、この加重平均Liを使用して所定の計測方式に対応する重み関数Ｗiを算出する方法とがある。
【０１２２】
この重み関数の計算方法は、voxelｉの時間分解型光路長liを迅速且つ直接計算することができるので、これを用いて時間分解型光路長liの加重平均Li或いは被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の各種の具体的な計測方法に対応する重み関数Ｗiを迅速に計算することができることになる。
【０１２３】
以上説明したように、本発明の散乱吸収体内部の光路分布計算方法は各種の光応答を用いた光ＣＴに適用することができる。この場合には、図５或いは図６に示すフローチャートの後段において、仮想媒体の任意のvoxe1ｉごとに算出される時間分解型光路長liを使用して、被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の計測方式に対応する重み関数Ｗiを算出する。
【０１２４】
例えば、具体例としては発明者らがＭＢＬに基づいて導出した時間積分計測（ＴＩＳ）、及び時間分解ゲート積分計測（ＴＧＳ）、周波数領域での位相変調計測（ＰＭＳ）、更には平均値法（ＴＩＳ）等がある。これらの場合、各種の光応答は、先に述べたインパルス応答の数学変換で定義されるが、それぞれに対応する光ＣＴでは、それぞれに異なった重み関数を用いる。
【０１２５】
以下、時間積分計測（ＴＩＳ）、時間分解ゲート積分計測（ＴＧＳ）、及び周波数領域での位相変調計測（ＰＭＳ）等の光応答を用いた光ＣＴに本発明の散乱吸収体内部の光路分布計算方法により求めた時間分解型光路長liを適用する例について説明する。
【０１２６】
４．ＴＩＳにおける散乱吸収体内部の光路分布と重み関数
先に述べたように、インパルス応答h(t)の時間積分Ｉを用いるＴＩＳに基づく光ＣＴでは、インパルス応答h(t)に対するvoxel i の平均光路長Liや光路長分散などを用いて記述される重み関数を用いる。なお、この重み関数Ｗｉは、吸収分布が一様な基準媒体に対する重み関数として定義される。
【０１２７】
まず、吸収が不均一に分布する散乱媒体のインパルス応答h(t)の時間積分Iを考える。この際に考えるモデルは、図1に示したモデルのvoxel i の吸収係数をμaiとしたものである。これらを用いて、インパルス応答の時間積分Ｉは、下記式で与えられる。
【数３２】
s(t)は、被計測媒体に吸収がないと仮定したときに検出位置において検出されるインパルス応答を示し、
Liは、検出位置において検出された時間積分Ｉを構成する光子集合のvoxel i における平均光路長（voxelｉにおける時間分解型光路長liの加重平均）を示す。］
この式では、ｌｎＩが吸収に無関係な項と、各voxelの吸収係数μaiに依存する項(減衰量)とに分離して記述されている。
【０１２８】
ここで、voxel i の平均光路長Liは、散乱体の散乱特性、境界条件、及び吸収や吸収分布に依存する量であり、下記式により与えられる。
【数３３】
対応する光ＣＴに用いる重み関数Ｗｉは、均一な吸収係数μνをもつ仮想媒体に対して定義され、下記式により与えれる。この場合のモデルは、図1に示したモデルのvoxel i の吸収係数をμνとしたものである。
【数３４】
Δμaiは、μνに対するμaiの吸収係数偏差を示し、
Ｌi(μν)は、一定の吸収係数μνをもつ仮想媒体のvoxe1ｉにおける平均光路長を示し（voxe1ｉ内の時間分解型光路長liの加重平均）、
Ｌ(μν)は、一定の吸収係数μνを与えた仮想媒体のインパルス応答の平均光路長を示し（仮想媒体内の光路長lの加重平均）、
Ｗiは、先に述べた（９）式中の減衰量偏差ΔＢi(μai)に対する重み関数であり、吸収係数偏差Δμai及びキュムラント(L, L',L",…)の関数で表される。］
ここで、キュムラントは、仮想媒体に対して求めた時間分解波形から直接算出することができる。
【０１３０】
５．ＴＧＳにおける散乱吸収体内部の光路分布と重み関数
まず、吸収が不均一に分布する散乱媒体のインパルス応答の時間分解ゲート積分信号Ｉgは、下記式で与えられる。
【数３８】
ここで、ゲート積分の時間範囲は[t1,t2]であり、これを[0,∞]にすると、前節に求めたインパルス応答h(t)の時間積分信号Iになる。
【０１３１】
対応する光ＣＴに用いる重み関数Ｗgiは、均一な吸収係数μνをもつ仮想媒体に対して定義され、下記式で与えられる。
【数３９】
Δμaiは、前述の（４．４）式に示したように、μνに対するμaiの吸収係数偏差を示し、
Ｌgi(μν)は、仮想媒体のインパルス応答の時間ゲート積分Ｉgを構成する光子集合のvoxel i における平均光路長（voxe1ｉ内のゲート時間内における時間分解型光路長liの加重平均）を示し、
Ｌg(μν)は、一定の吸収係数μνを与えた仮想媒体のインパルス応答の時間ゲート積分Ｉgを構成する光子集合の平均光路長（加重平均）を示し
Ｗiは、先に述べた（９）式の減衰量偏差ΔＢi(μai)をＴＧＳに適用したときの減衰量偏差ΔＢgi(μai)に対する重み関数であり、吸収係数偏差Δμai及びキュムラント(Lg, Lg',Lg",…)の関数で表される。］
【０１３２】
６．ＰＭＳにおける散乱吸収体内部の光路分布と重み関数
吸収が不均一に分布する散乱媒体に対して、ＰＭＳで得られる応答は、吸収が不均一に分布する散乱媒体に対するインパルス応答h(t)をフーリエ変換したものであり下記式で与えられる。
【数４２】
【数４４】
［式（６．１）〜（６．６）中、Ｒ及びＸはそれぞれ実数部と虚数部を示し、Ａ及びφはそれぞれ振幅及び位相遅れを示す。］
Ｒ、Ｘ、Ａ及びφは、ロックインアンプなどで容易に計測することができる。
【０１３３】
また、ここで、
【数４８】
【数５２】
【０１３４】
以上のように、ＰＭＳでは（６．７）式〜（６．１０）式に示す4種類の記述方法があり、それぞれの式に対応する光ＣＴ画像再構成法が可能である。その際、先に述べたのと同様にして、（６．７）式〜（６．１０）式にμai=μνを代入して求めた値と、（６．１１）式〜（６．１４）式の関係を用いて、それぞれの重み関数を算出することができる。
【０１３５】
以上、本発明の好適な実施形態について詳細に説明したが、本発明は上記実施形態に限定されない。
【０１３６】
【発明の効果】
以上説明したように、本発明によれば、被計測媒体と同一の散乱特性と境界条件をもつが、吸収がないと仮定した仮想媒体に対して順問題を解いて、voxe1ｉ内の時間分解型光路長liを直接且つ迅速に計算することができるので、従来の吸収があるときの順問題計算比べて、計算時間が短くなる。更に、各voxelの吸収係数を変化させてvoxe1の総数に等しい回数だけ順問題を繰り返して解くという時間のかかる従来の順問題計算が不要になる。従って、生体などの不均一系の散乱吸収体内部の光路分布を迅速に算出することのできる方法を提供することが可能となる。
【図面の簡単な説明】
【図１】被計測媒体の内部光子移動を解析するための有限格子モデルを示す模式図である。
【図２】図１に示す有限格子モデルにおいて、光入射位置ｐからインパルス光を入射する場合に、光検出位置ｑにおいて時間ｔの時点で検出されるインパルス応答を構成する光子集合（光路長がlになる）のうち、voxe1ｉを通過した光子集合の経路とその確率を示す模式図である。
【図３】図２に示す有限格子モデル内の時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)、時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕpi(t')、及び、時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕqi(t-t')の時間波形の一例を示すグラフである。
【図４】時間0≦t'≦tにおけるvoxe1ｉ内の光子存在累積確率Ｕi(t)、 時間t'におけるvoxe1ｉ内の光子存在確率密度Ｕi(t',t)、及び、時間t'における全voxe1内の光子存在確率密度Ｕpq(t',t)の関係を模式的に示すグラフである。
【図５】本発明の散乱吸収体内部の光路分布計算方法の一実施形態を示すフローチャートである。
【図６】本発明の散乱吸収体内部の光路分布計算方法の第二の実施形態を示すフローチャートである。
【符号の説明】
Ｓ１…第１ステップ、Ｓ２…第２ステップ、Ｓ３…第３ステップ、Ｓ４…第４ステップ、Ｓ５…第５ステップ、Ｓ６…第６ステップ、Ｓ７…第７ステップ、Ｓ８…第８ステップ、Ｓ１１…第１１ステップ。
Claims (6)
- 計測波長領域において、散乱吸収体である被計測媒体を、所定の大きさを有すると共に一様な吸収分布を有するとみなせるＮ（１≦Ｎ）個のvoxe1に分割する第１ステップと、
前記Ｎ個のvoxe1に分割された前記被計測媒体に対して、形状、境界条件及び散乱分布が前記被計測媒体と同一であり、屈折率分布が一様で、且つ吸収がないとみなせる仮想媒体を仮定する第２ステップと、
前記仮想媒体の表面に、インパルス光が入射する光入射位置と、前記光入射位置から入射して前記仮想媒体内を伝播した前記インパルス光の時間ｔにおけるインパルス応答s(t) を検出する光検出位置とを設定する第３ステップと、
前記仮想媒体内へ前記光入射位置からインパルス光を入射する場合に、前記光検出位置において検出される前記インパルス応答s(t)を算出する第４ステップと、
前記光入射位置から前記仮想媒体内にインパルス光を入射する場合に、前記光検出位置において時間ｔの時点で検出される前記インパルス応答s(t) を構成する光子集合が、前記光検出位置において検出される以前の時間t'(0≦t'≦t)の時点で前記仮想媒体内の任意のvoxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）内に存在していた光子存在確率密度Ｕi(t',t)を算出する第５ステップと、
前記光子存在確率密度Ｕi(t',t)を使用して前記光検出位置において時間ｔの時点で検出される前記インパルス応答s(t) を構成する光子集合が前記voxe1ｉ内に存在していた光子存在累積確率Ｕi(t)を算出する第６ステップと、
前記光子存在累積確率Ｕi(t)及び前記インパルス応答s(t)を使用して前記voxe1ｉの時間分解型光路長liを算出する第８ステップと、
を備えることを特徴する散乱吸収体内部の光路分布計算方法。 - 前記仮想媒体の任意のvoxe1ｉごとに算出される前記時間分解型光路長liの時間に対する加重平均Ｌi を算出する第９ステップを前記第８ステップの後段に更に備えることを特徴とする請求項１に記載の散乱吸収体内部の光路分布計算方法。
- 前記仮想媒体の任意のvoxe1ｉごとに算出される前記時間分解型光路長liを使用して、前記被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の計測方式に対応する重み関数Ｗiを算出する第１０ステップを前記第８ステップの後段に更に備えることを特徴とする請求項１に記載の散乱吸収体内部の光路分布計算方法。
- 計測波長領域において、散乱吸収体である被計測媒体を、所定の大きさを有すると共に一様な吸収分布を有するとみなせるＮ（１≦Ｎ）個のvoxe1に分割する第１ステップと、
前記Ｎ個のvoxe1に分割された前記被計測媒体に対して、形状、境界条件及び散乱分布が前記被計測媒体と同一であり、屈折率分布が一様で、且つ吸収がないとみなせる仮想媒体を仮定する第２ステップと、
前記仮想媒体の表面に、インパルス光が入射する光入射位置と、前記光入射位置から入射して前記仮想媒体内を伝播した前記インパルス光の時間ｔにおけるインパルス応答s(t) を検出する光検出位置とを設定する第３ステップと、
前記光入射位置から前記仮想媒体内にインパルス光を入射する場合に、前記光検出位置において時間ｔの時点で検出される前記インパルス応答s(t) を構成する光子集合が、前記光検出位置において検出される以前の時間t'(0≦t'≦t)の時点で前記仮想媒体内の任意のvoxe1ｉ（ｉ＝１，２，…，Ｎ；１≦Ｎ）内に存在していた光子存在確率密度Ｕi(t',t)を算出する第５ステップと、
前記光子存在確率密度Ｕi(t',t)を使用して前記光検出位置において時間ｔの時点で検出される前記インパルス応答を構成する光子集合が前記voxe1ｉ内に存在していた光子存在累積確率Ｕi(t)を算出する第６ステップと、
前記光子存在累積確率Ｕi(t)の全voxe1に対する総和Ｕ(t)を算出する第７ステップと、
前記光子存在累積確率Ｕi(t)及び前記光子存在累積確率Ｕi(t)の全voxe1に対する総和Ｕ(t)を使用して前記voxe1ｉの時間分解型光路長liを算出する第１１ステップと、
を備えることを特徴する散乱吸収体内部の光路分布計算方法。 - 前記仮想媒体の任意のvoxe1ｉごとに算出される前記時間分解型光路長liの時間に対する加重平均Ｌiを算出する第９ステップを前記第１１ステップの後段に更に備えることを特徴とする請求項４に記載の散乱吸収体内部の光路分布計算方法。
- 前記仮想媒体の任意のvoxe1ｉごとに算出される前記時間分解型光路長liを使用して、前記被計測媒体内部における吸収成分の濃度分布を定量計測するための所定の計測方式に対応する重み関数Ｗiを算出する第１０ステップを前記第１１ステップの後段に更に備えることを特徴とする請求項４に記載の散乱吸収体内部の光路分布計算方法。
Priority Applications (6)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2000078534A JP4364995B2 (ja) | 2000-03-21 | 2000-03-21 | 散乱吸収体内部の光路分布計算方法 |
| PCT/JP2001/001724 WO2001071319A1 (fr) | 2000-03-21 | 2001-03-06 | Procede de determination de la distribution de parcours optique a l'interieur d'un absorbant de diffusion |
| US10/239,245 US6975401B2 (en) | 2000-03-21 | 2001-03-06 | Method of calculating optical path distribution inside scattering absorber |
| AU2001236105A AU2001236105A1 (en) | 2000-03-21 | 2001-03-06 | Method of calculating optical path distribution inside scattering absorber |
| EP01908345A EP1284416B1 (en) | 2000-03-21 | 2001-03-06 | Method of calculating optical path distribution inside scattering absorber |
| DE60138543T DE60138543D1 (de) | 2000-03-21 | 2001-03-06 | Verfahren zur berechnung der optischen wegverteilung in einem streuenden absorbierer |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2000078534A JP4364995B2 (ja) | 2000-03-21 | 2000-03-21 | 散乱吸収体内部の光路分布計算方法 |
Publications (3)
| Publication Number | Publication Date |
|---|---|
| JP2001264245A JP2001264245A (ja) | 2001-09-26 |
| JP2001264245A5 JP2001264245A5 (ja) | 2006-12-28 |
| JP4364995B2 true JP4364995B2 (ja) | 2009-11-18 |
Family
ID=18595931
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2000078534A Expired - Lifetime JP4364995B2 (ja) | 2000-03-21 | 2000-03-21 | 散乱吸収体内部の光路分布計算方法 |
Country Status (6)
| Country | Link |
|---|---|
| US (1) | US6975401B2 (ja) |
| EP (1) | EP1284416B1 (ja) |
| JP (1) | JP4364995B2 (ja) |
| AU (1) | AU2001236105A1 (ja) |
| DE (1) | DE60138543D1 (ja) |
| WO (1) | WO2001071319A1 (ja) |
Families Citing this family (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20050106744A1 (en) * | 2001-12-26 | 2005-05-19 | Yoshihiko Mizushima | Optical analysis for heterogeneous medium |
| US7286214B2 (en) * | 2003-08-25 | 2007-10-23 | University Of South Florida | Method and program product for determining a radiance field in an optical environment |
| US7057182B2 (en) * | 2004-03-12 | 2006-06-06 | Hewlett-Packard Development Company, L.P. | Method and system for determining distortion in a circuit image |
| CN101313196A (zh) * | 2005-10-17 | 2008-11-26 | 阿而利克斯公司 | 利用空间调制的光学力显微术检测细胞变形性的设备和方法 |
| US8243071B2 (en) * | 2008-02-29 | 2012-08-14 | Microsoft Corporation | Modeling and rendering of heterogeneous translucent materials using the diffusion equation |
| JP5782314B2 (ja) | 2011-07-07 | 2015-09-24 | 浜松ホトニクス株式会社 | 生体計測装置および画像作成方法 |
| JP6154613B2 (ja) * | 2013-01-18 | 2017-06-28 | 浜松ホトニクス株式会社 | 断面画像計測装置及び計測方法 |
| JP6358573B2 (ja) | 2014-05-15 | 2018-07-18 | 国立大学法人浜松医科大学 | 乳房計測装置の作動方法及び乳房計測装置 |
| US11654635B2 (en) | 2019-04-18 | 2023-05-23 | The Research Foundation For Suny | Enhanced non-destructive testing in directed energy material processing |
| CN114204985A (zh) * | 2021-11-12 | 2022-03-18 | 西安理工大学 | 无线紫外光非直视通信中光子检测概率快速估算方法 |
Family Cites Families (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3310782B2 (ja) * | 1994-07-14 | 2002-08-05 | 株式会社日立製作所 | 吸収物質濃度の空間分布画像化装置 |
| JP3771339B2 (ja) * | 1996-01-18 | 2006-04-26 | 浜松ホトニクス株式会社 | 光ｃｔ装置及び光ｃｔによる画像再構成方法 |
| JP3662376B2 (ja) * | 1996-05-10 | 2005-06-22 | 浜松ホトニクス株式会社 | 内部特性分布の計測方法および装置 |
| JP3836941B2 (ja) * | 1997-05-22 | 2006-10-25 | 浜松ホトニクス株式会社 | 光ｃｔ装置及び画像再構成方法 |
| JPH11311569A (ja) * | 1998-04-28 | 1999-11-09 | Hamamatsu Photonics Kk | 内部特性分布の計測方法及び装置 |
| JP3887486B2 (ja) * | 1998-05-26 | 2007-02-28 | 浜松ホトニクス株式会社 | 散乱吸収体の内部特性分布の計測方法及び装置 |
-
2000
- 2000-03-21 JP JP2000078534A patent/JP4364995B2/ja not_active Expired - Lifetime
-
2001
- 2001-03-06 EP EP01908345A patent/EP1284416B1/en not_active Expired - Lifetime
- 2001-03-06 DE DE60138543T patent/DE60138543D1/de not_active Expired - Lifetime
- 2001-03-06 AU AU2001236105A patent/AU2001236105A1/en not_active Abandoned
- 2001-03-06 US US10/239,245 patent/US6975401B2/en not_active Expired - Lifetime
- 2001-03-06 WO PCT/JP2001/001724 patent/WO2001071319A1/ja active Application Filing
Also Published As
| Publication number | Publication date |
|---|---|
| EP1284416A4 (en) | 2004-10-20 |
| DE60138543D1 (de) | 2009-06-10 |
| AU2001236105A1 (en) | 2001-10-03 |
| WO2001071319A1 (fr) | 2001-09-27 |
| EP1284416B1 (en) | 2009-04-29 |
| US20030078503A1 (en) | 2003-04-24 |
| EP1284416A1 (en) | 2003-02-19 |
| US6975401B2 (en) | 2005-12-13 |
| JP2001264245A (ja) | 2001-09-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP3887486B2 (ja) | 散乱吸収体の内部特性分布の計測方法及び装置 | |
| Dehghani et al. | The effects of internal refractive index variation in near-infrared optical tomography: a finite element modelling approach | |
| EP2652708B1 (en) | A method and a system for image integration using constrained optimization for phase contrast imaging with an arrangement of gratings | |
| EP1303832B1 (en) | Imaging of scattering media using relative detector values | |
| US7965389B2 (en) | Method for reconstructing the distribution of fluorophores in a non-homogeneous medium by optical tomography in continuous mode | |
| US10223815B2 (en) | Iterative reconstruction method for spectral, phase-contrast imaging | |
| Florescu et al. | Single-scattering optical tomography | |
| JP4364995B2 (ja) | 散乱吸収体内部の光路分布計算方法 | |
| EP2577272A1 (en) | Low coherence enhanced backscattering tomography and techniques | |
| Hampel et al. | Fast image reconstruction for optical absorption tomography in media with radially symmetric boundaries | |
| CA2175455C (en) | Imaging system and method using direct reconstruction of scattered radiation | |
| Chang et al. | Layer-stripping approach for recovery of scattering media from time-resolved data | |
| US7046832B1 (en) | Imaging of scattering media using relative detector values | |
| Paithankar et al. | Fluorescence lifetime imaging with frequency-domain photon migration measurement | |
| JP2002333400A (ja) | 吸収物質濃度空間分布画像化装置 | |
| Stergar et al. | Assessment of line illumination for hyperspectral imaging by a Monte Carlo simulation | |
| Everitt et al. | Analysis and optimization of a diffuse photon optical tomography of turbid media | |
| JP2000105191A (ja) | 吸収物質濃度空間分布画像装置 | |
| Ueda et al. | Optical imaging reconstruction using the average value as the reference | |
| Jadick et al. | Cramér-Rao lower bound in the context of spectral x-ray imaging with propagation-based phase contrast | |
| Schiffers et al. | Polychromatic Maximum Likelihood Reconstruction for Talbot-Lau X-ray Tomography | |
| Graber et al. | A Layer-Stripping Approach for Recovery of Scattering Media from | |
| Chang et al. | A layer-stripping approach for recovery of scattering media from time-resolved data | |
| Freiberger | Nonlinear reconstruction schemes and experimental design for fluorescence diffuse optical tomography | |
| Tarvainen | Simon R. Arridge⋅ Jari P. Kaipio⋅ Ville Kolehmainen⋅ |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20061107 |
|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20061107 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20090818 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20090820 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120828 Year of fee payment: 3 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 4364995 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120828 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130828 Year of fee payment: 4 |
|
| EXPY | Cancellation because of completion of term |