JP3887486B2 - 散乱吸収体の内部特性分布の計測方法及び装置 - Google Patents
散乱吸収体の内部特性分布の計測方法及び装置 Download PDFInfo
- Publication number
- JP3887486B2 JP3887486B2 JP14430098A JP14430098A JP3887486B2 JP 3887486 B2 JP3887486 B2 JP 3887486B2 JP 14430098 A JP14430098 A JP 14430098A JP 14430098 A JP14430098 A JP 14430098A JP 3887486 B2 JP3887486 B2 JP 3887486B2
- Authority
- JP
- Japan
- Prior art keywords
- absorption coefficient
- light
- medium
- reference value
- value
- 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 - Fee Related
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
Landscapes
- Physics & Mathematics (AREA)
- Optics & Photonics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- General Health & Medical Sciences (AREA)
- General Physics & Mathematics (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
- Photometry And Measurement Of Optical Pulse Characteristics (AREA)
- Length Measuring Devices By Optical Means (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、散乱吸収体の内部特性分布を計測するための方法及びそのための装置に関する。より詳細には、本発明は、光入射位置及び光検出位置を測定対象物の表面に沿って移動させて内部の情報を得る装置に適用可能な内部特性分布の計測方法並びにそのための装置に関する。
【０００２】
【従来の技術】
光ＣＴ(computer tomography)は、生体内部の光学特性あるいは吸収成分の濃度分布を計測する手法又は装置を意味し、光を生体組織に入射したときにその中を移動して検出される光（信号光）を利用する。
【０００３】
現在までに直進光、準直進光、及び散乱光を利用する３種の光ＣＴ技術が報告されている。このうち、直進光を利用する方法は光の利用率が極めて悪く、極小さい媒体にしか適用できない。例えば、近赤外線を用いる場合、標準的な生体組織では１cmの厚さに対して透過光は約1/105倍に減衰する。これに比べて、散乱光を利用する方法は、媒体から出射する全ての光を利用するため、信号対雑音比が向上するので、より大きい媒体への適用が期待される。準直進光を利用する方法は、前記二つの方法の中間にある。そして、実用的な光ＣＴは、生体に対する最大許容光入射パワー、計測感度、所要計測時間などの制限から散乱光を利用する必要があるが、以下の技術上の問題に起因して上記光ＣＴは未だ実用化に至っていない。
【０００４】
第１の問題は、散乱吸収体内部の光あるいはフォトンの振る舞いを十分な精度で記述する方法が開発されていないことである。散乱吸収体内部を移動するフォトンの振る舞いの解析には、従来から、輸送方程式を近似したもの、あるいは輸送理論に拡散近似を適用した光拡散方程式が用いられてきた。しかし、拡散近似は媒体中のフォトンの平均自由光路長に対して十分に大きな媒体においてのみ成立するため、比較的小さい媒体、内部の複雑な形状の組織、及び複雑な形状の媒体を取り扱うことができない。また、拡散近似は等方散乱を前提としているため、非等方な散乱特性をもつ実際の生体組織の計測に適用すると、非等方散乱に起因する無視できない誤差が発生する。さらに、拡散方程式は、解析的あるいは数値的（有限要素法など）手法のいずれを用いても、前もって境界条件を設定しないとその解が求められない。つまり、計測に先立って、光入射・検出する各部の境界条件、つまり媒体形状や界面での反射特性を設定する必要がある。また、これらが個体差などによって変化すれば、それに応じて境界条件を変えて、計算をやりなおす必要がある。従って、輸送方程式の近似式や光拡散方程式から導出される信号光と散乱吸収体の光学特性との関係を利用する光ＣＴは、精度及び使い勝手に大きな問題がある。
【０００５】
また、上記とは別に、輸送方程式の近似式や光拡散方程式に摂動理論を適用して信号光と散乱吸収体の光学特性との関係を導出し、この関係を用いて光ＣＴ画像を再構成する方法がある。しかし、この方法は、非線形効果（２次以上の項）の取り扱いが極めて複雑になる。この際、２次以上の項を含めた計算は、理論的にはコンピュータで計算することが可能であるが、現在の最速のコンピュータを用いても演算時間が膨大なものになり、実用化は不可能である。このため、常套的には２次以上の項を無視している。そのため、この方法においては、複数個の比較的強い吸収領域が存在する媒体の光ＣＴ画像を再構成する際などに吸収領域の相互作用を無視することができなくなり、これに起因する大きな誤差が発生する。
【０００６】
第２の問題は、従来の光ＣＴが、狭義の重み関数、すなわち平均光路長あるいはそれに相当する位相遅れを利用していることである。このため、吸収係数に依存して変化する検出光の平均光路長の取り扱いが極めて複雑になる。そこで通常は近似を用いるが、近似を用いると誤差が増大するという重大な問題がある。このような狭義の重み関数を利用した方法に関しては、例えば以下の文献に記載されている。(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)。
【０００７】
なお、光ＣＴの画像再構成に際しては、信号光が散乱媒体のどの部分を通過してきたか、つまり媒体内の信号光の光路分布を知ることが最も重要である。この観点から、信号光の光路分布をランダムウォーク理論などで導出する方法もあるが、上記の問題は解決されていない。
【０００８】
以上のように、従来の光ＣＴにおいては十分な精度の再構成画像が得られず、空間解像度、画像ひずみ、定量性、測定感度、所要計測時間などに大きな問題がある。
【０００９】
本発明者は、以上のような状況を打破するため、次のことが重要であると考えて一連の研究を推進してきた。すなわち、光ＣＴを実現させるために重要なことは、強い散乱媒体である生体組織の中を移動する光の振る舞いを明らかにし、検出される信号光と吸収成分を含む散乱媒体（散乱吸収体）の光学特性との関係を明らかにし、前記信号光と前記関係とを利用して光ＣＴ画像を再構成するアルゴリズムを開発するということである。
【００１０】
そして本発明者は、▲１▼マイクロ・べア・ランバート則(Microscopic Beer-Lambert Law、以下「ＭＢＬ」という)に基づくモデルを提案すると共に、▲２▼散乱吸収体の光学特性と信号光との関係を表す解析式を導出し、例えば以下の文献に掲載している。(13) Y. Tsuchiya and T. Urakami: Jpn. J. Appl. Phys. 34, L79-81 (1995); (14) Y. Tsuchiya and T. Urakami: Jpn. J. Appl. Phys. 35, 4848-4851 (1996); (15) Y. Tsuchiya and T. Urakami: Optics Communications, 144, 269-280 (1997); (16) H. Zhang, M. Miwa, Y. Yamashita and Y. Tsuchiya: Ext. Abstr. Optics Japan '97, 30aA08 (1997)。このＭＢＬに基づく方法は原理的に、媒体形状、境界条件及び散乱の影響を受けないという大きな特長があり、非等方散乱媒体や小さい媒体への適用も可能である。この結果、散乱吸収体内の光子移動が明らかになり、この理論に基づいて散乱媒体の吸収係数あるいは吸収成分濃度を計測することが可能になった。
【００１１】
【発明が解決しようとする課題】
しかしながら、本発明者は、上記従来の方法はいずれも上記従来の狭義の重み関数を利用していたため、不均一系への適用が困難であり、生体などの散乱吸収体内の吸収成分についての計測精度が未だ充分なものではないことを見出した。
【００１２】
本発明は、上記従来の課題に鑑みてなされたものであり、生体などの不均一系の散乱吸収体について適用した場合であっても吸収成分についてより高い精度の計測が可能な方法及び装置を提供することを目的とする。
【００１３】
【課題を解決するための手段】
本発明者は、上記目的を達成すべく鋭意研究した結果、上記本発明者による知見をさらに発展させて、▲３▼上記マイクロ・べア・ランバート則を不均一系に適用し、▲４▼不均一散乱媒体の光学特性と信号光との関係を表す解析式を導出し、▲５▼前記信号光と不均一散乱吸収体の光学特性との関係から従来とは異なる定義の重み関数を導出し、▲６▼この重み関数を用いて光ＣＴ画像を再構成するためのアルゴリズム及び装置を提供することによって上記目的が達成されることを見出し、本発明に到達した。
【００１４】
すなわち、本発明の散乱吸収体の内部特性分布の計測方法は、
ヒトを除く散乱吸収体である被計測媒体中に一以上の光入射位置から順次光を入射する光入射ステップと、
前記被計測媒体内部を透過した光を複数の光検出位置で検出する光検出ステップと、
検出された各光に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記光の所定パラメータの測定値を取得する測定値取得ステップと、
前記被計測媒体の吸収係数の基準値を設定する基準値設定ステップと、
前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記パラメータの推定値を求める推定値算出ステップと、
前記吸収係数の基準値を用いて、複数のボクセルに分割された前記被計測媒体の各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて求める重み関数演算ステップと、
前記パラメータの測定値、前記パラメータの推定値及び前記重み関数に基づいて、前記各ボクセルにおける前記吸収係数の基準値に対する吸収係数の偏差を算出する吸収係数偏差算出ステップと、
前記吸収係数の基準値及び前記吸収係数の偏差に基づいて前記各ボクセルにおける吸収係数の絶対値を算出して前記被計測媒体における吸収係数の絶対値分布を求める吸収係数絶対値算出ステップと、
を含むことを特徴とする方法である。
【００１５】
また、本発明の散乱吸収体の内部特性分布の計測装置は、
散乱吸収体である被計測媒体中に一以上の光入射位置から順次光を入射する光入射手段と、
前記被計測媒体内部を透過した光を複数の光検出位置で検出する光検出手段と、
検出された各光に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記光の所定パラメータの測定値を取得する測定値取得手段と、
前記被計測媒体の吸収係数の基準値を設定する基準値設定手段と、
前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記パラメータの推定値を求める推定値算出手段と、
前記吸収係数の基準値を用いて、複数のボクセルに分割された前記被計測媒体の各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて求める重み関数演算手段と、
前記パラメータの測定値、前記パラメータの推定値及び前記重み関数に基づいて、前記各ボクセルにおける前記吸収係数の基準値に対する吸収係数の偏差を算出する吸収係数偏差算出手段と、
前記吸収係数の基準値及び前記吸収係数の偏差に基づいて前記各ボクセルにおける吸収係数の絶対値を算出して前記被計測媒体における吸収係数の絶対値分布を求める吸収係数絶対値算出手段と、
を備えることを特徴とする装置である。
【００１６】
本発明の方法及び装置においては、個々の計測に際して本発明によって初めて開示される各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて直接的に求め、その重み関数、所定パラメータの測定値及び前記パラメータの推定値に基づいて各ボクセルにおける吸収係数の偏差が算出される。この際、重み関数は計測状況とは無関係な一種の式で表わされる。このように本発明においては、個々の計測状況が変化しても適切な重み関数に基づいて吸収係数の偏差が算出されるため、近似法の採用に伴う誤差の発生が防止される。そのため、本発明の方法及び装置においては、生体などの不均一系の散乱吸収体について適用した場合であっても、各ボクセルにおける吸収係数の偏差が正確に求まることとなり、このような吸収係数偏差に基づいて得られる内部特性分布（光ＣＴ画像）の測定精度が向上する。
【００１７】
上記本発明の方法及び装置において採用される重み関数として、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける平均光路長と光路長分布の分散との関数が挙げられる。
【００１８】
この場合、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける平均光路長を求める平均光路長取得ステップ（平均光路長取得手段）を更に含んでいることが好ましく、前記重み関数演算ステップ（重み関数演算手段）において、下記の式：
【数９】
［式中、Ｗiは重み関数、Ｚi（μav）は平均光路長、Δμaiは吸収係数の偏差、σiv 2は光路長分布の分散をそれぞれ示す］
に基づいて前記重み関数を求めることが好ましい。
【００１９】
また、上記本発明の方法及び装置において採用される他の重み関数として、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける所定時間領域内の平均光路長と光路長分布の分散との関数が挙げられる。
【００２０】
この場合、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける所定時間領域内の平均光路長を求める平均光路長取得ステップ（平均光路長取得手段）を更に含んでいることが好ましく、前記重み関数演算ステップ（重み関数演算手段）において、下記の式：
【数１０】
［式中、Ｗiは重み関数、Ｚi（μav）は平均光路長、Δμaiは吸収係数の偏差、σiv 2は光路長分布の分散をそれぞれ示す］
に基づいて前記重み関数を求めることが好ましい。
【００２１】
更に、上記本発明の方法及び装置において採用される更に他の重み関数として、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける群遅延と分布の分散との関数が挙げられる。
【００２２】
この場合、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける群遅延を求める群遅延取得ステップ（群遅延取得手段）を更に含んでおり、前記重み関数演算ステップ（重み関数演算手段）において、下記の式：
【数１１】
［式中、Ｗiは重み関数、
【数１２】
はｃ（媒体内の光速）倍の群遅延、Δμaiは吸収係数の偏差、σf 2は分布の分散をそれぞれ示す］
に基づいて前記重み関数を求めることが好ましい。
【００２３】
上記本発明の方法及び装置に好適な重み関数を採用した場合、吸収係数に依存して変化する平均光路長又は群遅延のみならず分布の分散をも考慮されて重み関数が演算されるため、このような重み関数を用いて得られる内部特性分布の測定精度がより向上する傾向にある。
【００２４】
本発明の方法及び装置は、更に、前記吸収係数の絶対値を用いて前記各ボクセルにおける吸収成分の濃度を算出して前記被計測媒体における吸収成分の濃度分布を求める濃度算出ステップ（濃度算出手段）を含んでもよい。このような方法及び装置によれば、前述のように正確に求められた吸収係数偏差に基づいて吸収成分の濃度が求められるため、高精度の濃度分布が得られる。
【００２５】
上記本発明の方法及び装置を少なくとも２つの吸収成分を含有している被計測媒体に適用する場合には、前記光入射ステップ（光入射手段）において前記被計測媒体中に入射される光が、前記吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有していることが好ましい。この場合、前記光検出ステップ（光検出手段）において前記少なくとも２つの波長を有する光をそれぞれ検出し、前記測定値取得ステップ（測定値取得手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記測定値を取得し、前記基準値設定ステップ（基準値設定手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記基準値を設定し、前記推定値算出ステップ（推定値算出手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記推定値を求め、前記重み関数演算ステップ（重み関数演算手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記重み関数を求め、前記吸収係数偏差算出ステップ（吸収係数偏差算出手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収係数の偏差を算出し、前記吸収係数絶対値算出ステップ（吸収係数絶対値算出手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収係数の絶対値を算出し、前記濃度算出ステップ（濃度算出手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収成分の濃度を算出することによって、前記被計測媒体における前記各吸収成分の濃度分布が高精度で求められる。
【００２６】
上述の本発明の方法及び装置は、前記の求められた分布に基づいて、前記被計測媒体における該分布を示す光ＣＴ画像を表示する画像表示ステップ（画像表示手段）を更に含んでもよい。このような本発明の方法及び装置によれば、前述のように精度良く求められた内部特性分布を画像化して高精度な光ＣＴ画像として表示することが可能である。
【００２７】
【発明の実施の形態】
以下、図面を参照しつつ本発明の原理及び好適な実施形態について詳細に説明する。尚、図面中、同一又は相当部分には同一符号を付することとする。なお、散乱されつつ進行する光に関しては３次元座標を用いて考える必要があるが、以下では説明を簡単にするために場合により２次元座標を用いて説明する。先ず、本発明の原理について説明する。
【００２８】
［本発明の原理］
本発明は、マイクロ・ベア・ランバート則に基づいて導出される新規な重み関数を用いる。すなわち、本発明にかかわる重み関数は、従来の重み関数、つまり各ボクセル（voxel、体積素）内の光路長や平均光路長、光子滞在時間、あるいはそれらに対応する位相遅延量などとは異なる。本発明で採用する重み関数は、媒体形状、境界条件、光入射・検出位置及びそれらの間の距離、媒体の散乱係数などを繰り込んだ形で表される。また、この重み関数は散乱媒体に対して一意的に定まる。本発明では上記のような重み関数を用いて光ＣＴの画像再構成を行うため、後述するように、従来の問題点である誤差が大きく低減され、高精度の光ＣＴ画像を得ることができる。
【００２９】
１．光子移動モデル
先ず、本発明の基礎である散乱媒体内の光子移動のモデルを説明する。
【００３０】
Ａ）散乱媒体は再入射不可能な任意の３次元形状をもつ。従って、媒体から出た光子が再度その媒体に入射することはない。この制限内で、散乱媒体と外部媒体との境界条件、例えば形状や屈折率のマッチング条件などは任意である。
【００３１】
Ｂ）散乱媒体内で、光子と媒体との非線形作用や、光子と光子との相互作用はない。従って、媒体の光学定数は光子密度に依存しない。
【００３２】
Ｃ）マクロに見たとき、屈折特性は一様（均一）である。しかし、吸収係数分布は一様でない（不均一である）。以下、このような媒体を便宜的に不均一媒体あるいは不均一系とよぶ。従って、以下で議論する散乱媒体内の光速度は一定値ｃになる。
【００３３】
Ｄ）上記不均一散乱媒体に対して、光入射・検出位置ｐ（ｕ，ｖ）（表面位置ｕに光子を入射し、他の表面位置ｖで散乱光子を検出する）を定めたとき、光路分布（μa＝０）が一意的に決まる。この光路分布（μa＝０）は、散乱媒体に吸収がないと仮定したときに光子が取り得る光路の分布あるいは光路長分布である。なお、このように吸収がないと仮定した散乱媒体を「仮想媒体」という。
【００３４】
Ｅ）上記の光路分布（μa＝０）は、散乱媒体に吸収がないと仮定した仮想媒体に対して定義したから、散乱媒体の吸収あるいは吸収分布には影響されない。また、この光路分布（μa＝０）は、散乱媒体の散乱特性、媒体形状、境界条件、及び光入射・検出位置ｐ（ｕ，ｖ）に依存する。そして、これらの影響は光路分布（μa＝０）の中に繰り込まれている。つまり、散乱等の影響は、吸収がない仮想媒体に対する光路分布（μa＝０）に反映されるため、以下では光路分布（μa＝０）やμa＝０である仮想媒体に対する応答を、単に散乱等の影響とよぶこともある。
【００３５】
Ｆ）また、前記Ｂ）項と関連するが、上記の光路分布（μa＝０）は、複数あるいは多数光子（光子集合）入射や多数光子の移動に対して不変である。
【００３６】
また、媒体内の光子移動と減衰に関しては、以下の関係：
Ｇ）吸収係数がμaである均一媒体内を光路長ｌ(エル)＝ｃｔだけ移動する単一光子に対する減衰量(attenuation)はμaｌとなるという関係（以下、この関係を「ＭＢＬ第１基本則」とよぶ）；及び
Ｈ）散乱媒体内の多数の光子（光子集合）に対する減衰量は、それぞれの光子に対する減衰量の重ね合せ(superposition)で表されるという関係；
が成立する。
【００３７】
２．不均一系内の光子移動に関する基礎事項
最初に不均一系内の光子移動に関する基礎事項を説明する。
【００３８】
2.1 不均一媒体内の光子移動
吸収（又は吸収成分）が不均一に分布する３次元の不均一媒体全体をＮ個のボクセルに分割して番号ｉを付け、それぞれのボクセルｉについて非等方散乱係数μs、散乱角の余弦の平均値ｇ、吸収係数μaiを定義する（図１参照）。光子入射・検出位置はｐ（ｕ，ｖ）である。ここで、μsとμaiは不均一媒体固有の値で、かつμsは定数、またμaiは位置の関数である。この際、ボクセルｉは任意の大きさと形状にすることができ、大きい体積をもつボクセルではそのボクセル内の平均吸収係数や平均光路長を考えればよい。また、空間的に離れた２個以上のボクセルを合わせて、１個のボクセルと考えてもよい。
【００３９】
いま、単一光子を繰り返して入射するとして、ｍ番目の入射光子（m-th光子）に対して時間ｔに検出される出力ｈm(t)を考える。そして、このm-th光子のボクセルｉにおける滞在（飛行）時間をｔim、また対応する光路長をｌim（＝ｃｔim）とする。この場合、検出される光子の光路長ｌとｌimとの関係は、ｌ＝ｃｔ＝Σｌimである。ここで、添字ｍはm-th光子に関する量を表す。上記の系でm-th光子の光路、つまり光路分布（μa＝０）が決まると、ｌimが一意的に決まる。また、ボクセル番号ｉを決めると、そのボクセルｉの吸収係数μaiが一意的に決まる。この場合、吸収係数は任意に選べるから、ボクセルの光路長と吸収係数は独立である。
【００４０】
時間ｔに検出されるm-th光子の減衰量（減光量ともよばれる）Ｂmは、
【数１３】
と表される。但し、光子が通過しなかったボクセルｉに対して、ｌim＝０すなわちｔim＝０を定義した（光路長分布がゼロ）。従って、(2.1.1)式の中のＢmには光子が通過しなかったボクセルの影響は含まれていない。また、μaiとｌimはともに非負の値をもつから、Ｂm＝Σ（μaiｌim）＝０となるのはｌim≠０であるボクセルｉの吸収係数がμai＝０のときに限られる。この際、ｌim＝０となるボクセルの影響は(2.1.1)式の中に含まれていないから、Σ（μaiｌim）＝０なる条件はΣμai＝０と置き換えてよい。
【００４１】
この(2.1.1)式は、不均一系におけるＭＢＬ第１基本則、つまり「吸収係数がμaiである散乱体の中をジグザグ光路に沿って長さｌ＝Σｌimだけ移動した光子の減衰量は、媒体の散乱特性とは無関係な値Σ（μaiｌim）になる」ことを表す。この関係は、光子移動の際に吸収係数が変化する位置でＭＢＬを接続したものである。また、この関係は、ジグザグ光路が交差（光子が２回以上同一場所を通過）する場合にも成立する。
【００４２】
ここで、Σμai＝０のときのm-th光子に対する応答をｓm（μs，ｔ）とおけば、m-th光子に対する応答ｈm(t)は、
【数１４】
となる。各ボクセルの光路長と吸収係数は独立であるから、
【数１５】
が得られる。但し、上式中のＢimは、ボクセルｉにおける減衰量である。この(2.1.3)式は、不均一系内のm-th光子に関するＭＢＬを表わす。
【００４３】
すると、(2.1.3)式から、不均一系内のm-th光子に対して、
【数１６】
が得られる。
【００４４】
不均一媒体に多数の光子（光子集合）を入射したときの応答を構成する光子で、光路長がｌとなる光子を考える。この場合、光路長がｌとなる光路は媒体内に無限といえるほど多数存在する。そして、多数の光子（光子集合）に対する減衰量は、単一光子に対する減衰量の重ね合せ(superposition)で表される。また、通常の応答、例えばインパルス応答には、光路長の異なる光子が含まれている。この場合にも、減衰量は単一光子に対する減衰量の重ね合せで表される。また、先に定義した光路分布（μa＝０）は、多数の光子入射に対して不変である。従って、多数光子入射に対する不均一系の各種応答の減衰量は、単一光子に対する減衰量の重ね合せ、つまり線形問題として解析することができる。
【００４５】
2.2 不均一媒体の光路分布（μa＝０）
再入射不可能な任意の不均一媒体に対して光入射・検出位置ｐ（ｕ，ｖ）を定めたとき、光子が取り得る光路分布は一意的に決まる。この光路分布（μa＝０）は、その不均一媒体に吸収がない（μa＝０）と仮定したとき（仮想媒体、μa＝０）の光路分布で定義されるから、不均一媒体内の吸収やその分布には影響されない。しかし、光路分布（μa＝０）は、不均一媒体の光学特性、媒体形状や境界条件、光入射・検出位置ｐ（ｕ，ｖ）、及び光路長に依存する。
【００４６】
不均一系に対するインパルス応答、インパルス応答の時間分解ゲート積分信号、及び周波数領域の応答を構成する光子集合に対して、それぞれの光路分布（μa＝０）を定義することができる。また、光路分布（μa＝０）が決まれば、各ボクセルの平均光路長（μa＝０）が一意的に決まる。
【００４７】
以上のような光路分布（μa＝０）は、モンテカルロシミュレーションなどによって計算することができる。さらに、モンテカルロシミュレーションでは、光路分布（μa＝０）を基準にして、任意の一様な吸収（例えばμa＝μav）がある場合の検出光の光路分布や平均光路長を計算することもできる。
【００４８】
2.3 検出光の光路分布
インパルス応答には、種々の光路長をもつ光子が含まれる。媒体に吸収があるとき、個々の光子の減衰量はＢm＝Σ（μaiｌim）になるが、光路長が異なる光子集合の平均光路長は、吸収係数の関数になる。従って、応答の時間波形は、吸収に依存して変化することになる。つまり、吸収があると光子移動の途中で光子が消滅し、その消滅の程度が光路長あるいは光路分布（μa＝０）によって異なるため、その分だけ検出光に対する平均光路長や光路分布が変化する。この結果、検出光の光路分布は、通常、吸収に依存することになり、前出の光路分布（μa＝０）に等しくならない。
【００４９】
2.4 重み関数
不均一媒体内をｌだけ移動した単一光子の減衰量は、前述したようにΣ（μaiｌim）で表される。そして、不均一媒体における多数光子の減衰量は、単一光子の減衰量の重ね合せで表される。ここで、多光子で構成される不均一媒体の各種応答の減衰量をボクセルｉにおける寄与係数Ｗiとボクセルｉの吸収係数μaiとの積で表すことを考える。そして、前記の寄与係数Ｗiを「重み関数(Weight Function)」とよぶ。
【００５０】
前出の(2.1.4)式と(2.1.5)式の場合には、光路長ｌ＝ｃｔ＝Σｌimが吸収係数μaiの関数にならないから、光路分布（μa＝０）と重み関数は等しい。しかし、インパルス応答の時間積分値のように平均光路長が吸収係数の関数になる場合、後述するように、重み関数は光路分布（μa＝０）や検出光の光路分布に等しくない。また、一般に重み関数Ｗiは吸収に依存する。
【００５１】
３．不均一系における各種の応答関数
前述のことを踏まえて、不均一系における各種の応答を求めて、不均一系に対するＭＢＬを導出する。この際、不均一媒体に対する光入射・検出位置をｐ（ｕ，ｖ）とする。すると、時間領域や周波数領域の各種計測で得られる各種応答に対して、それぞれの光路分布（μa＝０）が一意的に決まる。そして、各種応答の減衰量は、各種応答を構成する単一光子に対する減衰量の重ね合せになる。なお、光路分布（μa＝０）は吸収分布に対して不変であるから、計測系の設定、つまり光入射・検出位置ｐ（ｕ，ｖ）や計測方式が同一であれば、被計測媒体が不均一媒体であっても均一媒体であっても、同一の光路分布（μa＝０）になる。ＭＢＬに基づく記述・解析方法は、不均一系の応答の場合でも、μa＝０としたときの応答（つまり散乱や境界条件の影響）と吸収による減衰量とを分離した形で記述することができるという点で優れている。
【００５２】
3.1 不均一系のインパルス応答
不均一系のインパルス応答ｈ(t)は、種々の光路長をもつ多数の検出光子で構成される。この場合、観測時間ｔを決めるとｌ＝ｃｔとして光路長ｌが決まる。いま、光路長がｌである検出光子の数Ｍが十分に大きいとする。この場合、全ての光子の光路長はｌm＝ｌであるから、Ｍ個の光子に対する平均光路長もｌ、つまり、
【数１７】
である。このような光路長がｌとなる光子集合の光路分布（μa＝０）は、媒体の形状や散乱特性及び光入射・検出位置ｐ（ｕ，ｖ）に対して一意的に決まり、不均一媒体の吸収係数や吸収分布には関係しない。そこで、μa＝０のときのインパルス応答に対して光路長がｌとなる十分に大きい数Ｍ個の光子（光子集合）を考え、この光子集合の光路分布（μa＝０）を各ボクセルｉ（ｉ＝１〜Ｎ）の平均光路長ｌi（μa＝０）で表す。ここで(3.1.1)式を参照すると、この平均光路長ｌiに関して、
【数１８】
が成立する。つまり、多数光子に対する各ボクセルの平均光路長ｌiの全ボクセルに対する総和は多数光子に対する平均光路長ｌに等しい。
【００５３】
また、インパルス応答を構成する光子（複数）の中で光路長がｌ（エル）となる光子集合に対する減衰量Ｂhとボクセルｉにおける減衰量Ｂihは、
【数１９】
となる。従って、Σμai＝０のときのインパルス応答をｓ（μs，ｔ）とおけば、不均一系のインパルス応答ｈ(t)は、
【数２０】
と表される。すると、不均一系のインパルス応答に対するＭＢＬ、つまり
【数２１】
が得られる。
【００５４】
以上から、不均一系のインパルス応答ｈ(t)は、
【数２２】
となる。これらの式では、先に述べたように、μa＝０としたときのインパルス応答、つまり媒体形状や散乱の影響を表す項ｓ（μs，ｔ）と、吸収に依存する減衰量を表す項Ｂhとが分離されている。
【００５５】
なお、上記(3.1.7)式から、不均一系のインパルス応答ｈ(t)に対する重み関数は、先に定義した光路分布（μa＝０）で決まるボクセルｉ（ｉ＝１〜Ｎ）の光路長ｌi（μa＝０）に等しいことがわかる。従って、時間分解計測による光ＣＴ画像再構成アルゴリズムは、比較的簡単になる。
【００５６】
以下では、上記で求めた不均一系のインパルス応答ｈ(t)を構成する光子集合の光路分布（μa＝０）、つまり前記で定義した各ボクセルの光路長ｌi（μa＝０）、あるいはそれらの時間分布（異なる値をもつそれぞれのｌ（μa＝０）に対するｌi（μa＝０）であり、時間分解型の光路分布（μa＝０）とよぶ）を用いて各種の応答とそれらに対する重み関数を導出する。
【００５７】
3.2 不均一系のインパルス応答の時間分解ゲート積分信号
ここでは、先ず不均一系のインパルス応答の時間分解ゲート積分信号に対する解析式を導出し、減衰量や平均光路長の吸収依存性を考慮した重み関数を導出する。なお、時間域の積分範囲［ｔ1，ｔ2］を［０，∞］にすると、定常（ＣＷ）光に対する応答を表すことになる。
【００５８】
3.2.1 時間分解ゲート積分信号の解析式
インパルス応答の時間分解ゲート積分信号ＩTは、(3.1.8)式を積分して
【数２３】
と表される。
【００５９】
この式から、
【数２４】
が得られる。ここで〈ｌiT〉＝Ｌi（μa）はインパルス応答の時間分解ゲート積分信号に対するボクセルｉの平均光路長である。このＬi（μa）は媒体の散乱特性や境界条件、及びｔ1，ｔ2に依存する。また、後出の(3.2.9)式に示すように、Ｌi（μa）はμaiに対する単調非増加関数である。さらに、上記のＬi（μa）はボクセルｉの吸収係数μaiの関数であるが、ボクセルｉを含めた全ボクセルｉ（ｉ＝１〜Ｎ）の吸収係数μaiの関数でもある。このようなＬi（μa）の全ボクセルの吸収に対する依存性については後述する。
【００６０】
すると、(3.2.2)式から
【数２５】
が得られる。つまり、インパルス応答の時間分解ゲート積分信号ＩTに対する各ボクセルの平均光路長〈ｌiT〉＝Ｌi（μa）の全ボクセルに対する総和は〈ｌT〉に等しい。
【００６１】
次に、(3.2.2)式をμaiで積分して、時間分解ゲート積分信号ＩTに対する基本式を求めると、
【数２６】
が得られる。右辺第１項は積分定数であり、(3.2.1)式にΣμai＝０を代入して得られる。この際、右辺第１項は吸収とは無関係、つまり散乱や境界の影響を表し、右辺第２項は減衰量ＢTを表す。ここで、上記までの結果をまとめると、
【数２７】
が得られる。
【００６２】
ここで、先に述べたＬi（μa）の全ボクセルの吸収係数に対する依存性が問題になる。つまり、この依存性があると、(3.2.5)式をこれ以上簡単にすることができない。この事実は重要であり、光ＣＴ画像再構成アルゴリズムを複雑にしている。
【００６３】
ところで一般には、各ボクセル間の吸収係数μaiの差、つまり平均吸収係数からの偏差はそれほど大きくない場合が多い。この場合、注目しているボクセルｉの平均光路長Ｌi（μa）に対する他部分（又は他のボクセル）の吸収係数偏差の影響は小さい。換言すれば、各ボクセル間の吸収係数の差が小さい場合、注目しているボクセルｉの平均光路長Ｌi（μa）は、全ボクセルに対する吸収係数の平均値で大まかに決まり、他部分の吸収係数の平均吸収係数に対する偏差に依存するＬi（μa）の変化は無視できる程度に小さく、Ｌi（μa）はボクセルｉの吸収係数μaiのみの関数であると見なしてよい。この場合、ボクセルｉの吸収係数の平均吸収係数との差（吸収係数偏差）に対応するＬi（μa）の差分は、線形近似して記述することができる。従って、以下では、Ｌi（μa）はボクセルｉの吸収係数μaiのみの関数であると仮定する。つまり、ボクセルｉの平均光路長Ｌi（μa）を考える際に、ボクセルｉの吸収係数μai以外のパラメータは一定（不変）であると仮定する。なお、２つ以上のボクセルの吸収係数が同時に変化する場合、これらの２つのボクセルを合わせたボクセルを１つの新しいボクセルと考える。このことは、「１つの測定値から複数ボクセルの吸収係数の変化を定量することはできない」ことと等価であり、理論上の矛盾はない。
【００６４】
以上の線形近似を適用すると、ボクセルｉにおける減衰量ＢiTを定義すると、
【数２８】
が得られる。この式は、不均一系の時間分解ゲート積分信号ＩTを構成する光子集合に対するＭＢＬを表す。但し、先に述べた線形近似を用いた。また、前述したように、積分範囲［ｔ1，ｔ2］を［０，∞］にすると、通常の時間分解積分信号、つまり定常（ＣＷ）光に対するＭＢＬを表す。
【００６５】
次に、時間分解ゲート積分信号ＩTに対する重み関数Ｗiについて考える。先ず、平均値定理を用いて(3.2.4)式と(3.2.6)式の減衰量ＢTとＢiTを変形すると、
【数２９】
となる。但し、μx0は０≦μx0≦μaiなる適宜の値である。従って、2.4節に記載した重み関数の定義から、時間分解ゲート積分信号ＩTに対するボクセルｉの重み関数はＷi0＝Ｌi（μx0）であることがわかる。この重み関数Ｗi0は、時間分解ゲート積分信号に対するボクセルｉの平均光路長Ｌi（μai）＝〈ｌiT〉とは異なる。
【００６６】
なお、前記Ｌi（μa）の全ボクセルの吸収に対する依存性が問題になる場合には、(3.2.6)式の代わりに(3.2.5)式を使用する必要がある。
【００６７】
以上のことも、積分範囲［ｔ1，ｔ2］を［０，∞］にすることによって、通常の時間分解積分信号、つまり定常（ＣＷ）光に対する応答に適用することができる。
【００６９】
3.2.2 平均光路長の吸収依存性
各ボクセルの平均光路長は、そのボクセルの吸収に依存する。また、重み関数もそのボクセルの吸収に依存する。本節ではボクセルｉの吸収係数μaiが（μai＋ｈ）に変化（Δμai＝ｈ）したときの減衰量や平均光路長の変化を解析して、それらの変化量を表す解析式を導出する。この結果、光ＣＴの画像再構成に際して問題となる平均光路長と重み関数の差、及びそれらの吸収依存性が明らかになる。
【００７０】
先ず、(3.2.7.2)式に示したボクセルｉの減衰量ＢiTをμaiの周りにテイラー展開すると、
【数３４】
となる。但し、ＢiT (n)（μa）はＢiT（μa）のμaに対するｎ階導関数である。また、(3.2.6)式の関係から、
【数３５】
が得られる。従って、(3.2.2)式を用いてＢiT (2)（μai）を求めると、
【数３６】
となる。つまり、−ＢiT (2)（μai）＝−Ｌi (1)（μai）はｌiTの平均値まわりの分散を示す。従って、減衰量ＢiT（μai）は、
【数３７】
となる。
【００７１】
ここで、吸収係数がμaiから（μai＋ｈ）に変化したときの減衰量の差ΔＢiTを(3.2.7.2)式から求めると、
【数３８】
となる。このμxは前出のμx0と異なり、min（μai，μai＋ｈ）≦μx≦max（μai，μai＋ｈ）なる適宜の値である。また、このＬi（μx）は減衰量の差ΔＢiTに対する重み関数Ｗiである。ここで、(3.2.11)式と(3.2.12)式を比べると、
【数３９】
となる。
【００７２】
他方、平均光路長Ｌi（μai）をμaiの周りにテイラー展開すると、
【数４０】
となる。但し、Ｌi (n)（μa）はＬi（μa）のμaに対するｎ階導関数である。従って、(3.2.13)式及び(3.2.14)式から、
【数４１】
が得られる。なお、(3.2.15)式の最右辺（第４式）に示す近似式は、本発明者によって既に均一散乱媒体内の吸収成分の濃度定量に用いられ、その有効性が確認されている。
【００７３】
以上を踏まえてボクセルｉの吸収係数μaiが（μai＋ｈ）に変化したとき、(3.2.4)式を用いて変化前後の時間分解ゲート積分信号lnＩTの差ΔＩを求めると、
【数４２】
となる。つまり、インパルス応答の時間分解ゲート積分信号lnＩTの差ΔＩと、吸収係数がμaiであるときの重み関数Ｗiから、ボクセルｉの吸収係数変化分ｈ＝Δμaiを求めることができる。この場合、重み関数Ｗiは平均光路長Ｌi（μai）＝〈ｌiT〉及び光路長分布の分散σi 2の関数である。
【００７４】
また、複数ボクセルｉ（１〜Ｎの中の複数の値）の吸収係数μaiが（μai＋ｈi）に変化したとき、(3.2.4)式を用いて変化前後の時間分解ゲート積分信号lnＩTの差ΔＩを求めると、
【数４３】
となる。つまり、インパルス応答の時間分解ゲート積分信号lnＩTの差ΔＩと、吸収係数がμaiであるときの重み関数Ｗiから、複数のボクセルｉ（ｉ＝１〜Ｎ）の吸収係数変化分ｈi＝Δμaiを求めることができる。但し、この場合にはＮ個の連立方程式を解く必要があり、Ｎ個の独立な計測値ΔＩが必要である。また、前記の重み関数Ｗiは平均光路長Ｌi（μai）＝〈ｌiT〉及び光路長分布の分散σi 2の関数である。
【００７５】
以上の場合、吸収係数μaiを被計測媒体の平均吸収係数あるいはそれに近い値に選んでおくと、精度が向上する。つまり、さきに述べたＬi（μa）の全ボクセルの吸収に対する依存性の問題が大きく緩和される。
【００７６】
なお、各ボクセルの吸収係数μaiの平均吸収係数に対する偏差が極端に大きい場合で、他のボクセルの吸収係数偏差に対するＬi（μa）の依存性が問題になる場合には、(3.2.5)式に戻って再計算をするような繰り返し計算が必要になる。
【００７７】
以上のことは、積分範囲［ｔ1，ｔ2］を［０，∞］にすると、通常の時間分解積分信号、つまり定常（ＣＷ）光に対する応答に適用することができる。
【００７８】
3.2.3 アベレージ・バリュー・メソッドへの応用
光ＣＴは、不均一媒体中の吸収成分の濃度分布を計測しようとするものである。この場合、先に述べた２種類の平均光路長の吸収依存性、つまりＬi（μa）の全ボクセルの吸収に対する依存性、及びＬi（μa）のボクセルｉ自身の吸収に対する依存性が問題を複雑にしている。本発明者は、これらの吸収依存性に基づく誤差を低減する一つの方法として、アベレージ・バリュー・メソッド（Average Value Method、以下「ＡＶＭ」という）を開発した。このＡＶＭは、計測すべき不均一媒体に対して近似的な平均吸収係数を推定あるいは計測して、この値を基準にして、媒体各部の吸収係数の偏差を定量する。この場合、上記で推定あるいは計測する平均吸収係数は、真値に等しい必要はない。この方法は、平均光路長の吸収係数依存性の問題を大きく緩和するもので、前記で推定あるいは計測した平均吸収係数が真値に近いほど、再構成光ＣＴ画像の精度が向上する。
【００７９】
いま、被計測媒体の平均吸収係数に近い値の吸収係数μavをもつ仮想媒体（μa＝μav）を仮定する。すると計測対象である不均一媒体の吸収係数はμai＝μav＋Δμai（ｉ＝１〜Ｎ）となる。時間分解ゲート積分計測を考えると、仮想媒体に対する計測値Ｉ（μav）と実際の媒体に対する計測値 Ｉ（μai＝μav＋Δμai）（ｉ＝１〜Ｎ）との関係は(3.2.4)式及び(3.2.15)式、あるいは(3.2.17)式から、
【数４４】
となる。つまり、仮想媒体と実際の媒体に対するインパルス応答の時間分解ゲート積分信号lnＩTの差ΔＩと、吸収係数がμavであるときの重み関数Ｗiから、複数のボクセルｉ（ｉ＝１〜Ｎ）の吸収係数偏差Δμaiを求めることができる。但し、この場合にはＮ個の連立方程式を解く必要があり、Ｎ個の独立な計測値ΔＩが必要である。また、前記重み関数Ｗiは、吸収係数がμavであるときの平均光路長Ｌi（μav）＝〈ｌiT〉とその分散σiv 2との関数である。なお、上記のような一様な吸収がある媒体（仮想媒体）のボクセルｉの平均光路長及び分散は、モンテカルロシミュレーションなどで計算することができる。
【００８０】
前記(3.2.18)式を行列式で表せば、
【数４５】
となる。但し、重み関数Ｗiは、
【数４６】
である。
【００８１】
3.3 不均一系の周波数領域の応答
周波数領域の応答は、前出の(3.1.8)式に示すインパルス応答ｈ(t)をフーリエ変換して、
【数４７】
となる。ここで、Ｒ，Ｘ，Ａ，φはそれぞれ、
【数４８】
である。すなわち、Ｒ及びＸはそれぞれ実部及び虚部であり、Ａ及びφはそれぞれ振幅及び位相遅れであり、これらはロックインアンプなどで容易に計測することができる。
【００８２】
すると、後で詳述するように、
【数４９】
なる関係が得られる。ここで、各式の右辺の添字ｉが付いたパラメータは、ボクセルｉにおける量であることを示し、前節に述べた平均光路長に相当する。また、ここでは、後で詳述するように、
【数５０】
を定義した。以上では、
【数５１】
の関係が成り立つ。
【００８３】
従って、(3.3.3.1)式ないし(3.3.3.4)式の関係式から、前出の(3.2.4)式の導出と同様にして、
【数５２】
が得られる。右辺第１項は積分定数であり、それぞれ、左辺に示すパラメータにΣμai＝０を代入して得られる。
【００８４】
以上の(3.3.6.1)式ないし(3.3.6.4)式は、前節に記載した時間分解ゲート積分信号に対する式と相似であるから、前節と同様にして、これらの式を変形することができる。そこで以下では、(3.3.6.3)式を例として取り上げる。先ず、前節と同様に、周波数領域の応答に対する重み関数を考える、先ず、平均値定理を用いると(3.3.6.3)式の減衰項（右辺第２項）は、
【数５３】
となる。ここで、
【数５４】
は検出光に対するボクセルｉの重み関数（c×群遅延）、つまり周波数領域の応答の重み関数を表す。但し、μxは０≦μx≦μaiなる適宜の値である。この重み関数(3.3.7.3)は、前節に述べた重み関数Ｗi（μx）＝Ｌi（μx）に相当し、ボクセルｉの平均群遅延
【数５５】
のｃ倍の値
【数５６】
とは異なる値をもつ。
【００８５】
さらに、ボクセルｉの減衰量Ｂif及び平均群遅延(3.3.7.4)の吸収係数依存性は、前節に述べた減衰量ＢiT及び平均光路長Ｌi（μx）と同様になる。つまり、群遅延(3.3.7.4)の吸収依存性を考えると、(3.3.3.4)式及び(3.3.5.2)式より、
【数５７】
となる。但し、異なるボクセル間において（ｉ≠ｊのとき）、
【数５８】
を仮定した。つまり、平均群遅延(3.3.7.4)はボクセルｉ以外のボクセルの減群遅延
【数５９】
に依存しないとした。これらは、前節における仮定と同様である。なお、上式中の
【数６０】
は、３種の変調周波数を利用して計測することができる。
【００８６】
以上から、前述の説明（3.2.3節）と同様にして、
【数６１】
が得られる。この際、一様な吸収がある媒体（仮想媒体）（μai＝μav）のボクセルｉの平均群遅延及び分散は、モンテカルロシミュレーションなどで計算することができる。従って、μaiが既知であるとき、吸収係数がμxであるときの平均群遅延
【数６２】
を計算で求めることができる。
【００８７】
［(3.3.3.1)式〜(3.3.3.4)式の導出］
先ず、ｔ＝Σｔiであることを考慮して、(3.3.2.1)式及び(3.3.2.2)式から、
【数６３】
が得られる。そこで、
【数６４】
を定義する。この場合、
【数６５】
の関係が成り立つ。また、(3.3.2.1)式及び(3.3.2.2)式から、
【数６６】
が得られる。以上から、
【数６７】
が得られる。
【００８８】
さらに、(3.3.2.3)式及び上式から、
【数６８】
が得られる。但し、ここでは、
【数６９】
を定義した。この場合、上式に(B.3.1)式及び(B.3.2)式を代入すると、
【数７０】
の関係が得られ、
【数７１】
となる。以上により、(3.3.3.3)式が得られる。
【００８９】
また、同様にして、(3.3.2.4)式、(3.3.3.1)式及び(3.3.3.2)式から
【数７２】
が得られる。但し、ここでは、
【数７３】
を定義した。従って、
【数７４】
が得られる。
【００９０】
なお、上記のとき、
【数７５】
が成立する。従って、
【数７６】
の関係が成り立つ。
【００９１】
４．不均一系に関する一般式
前章までの説明から明らかなように、不均一系におけるＭＢＬは、微分形で表すと有効である。つまり、不均一系に対するインパルス応答ｈ(t)の減衰量Ｂhを表す(3.1.6)式、インパルス応答の時間分解ゲート積分信号ＩT（μs、μa）の減衰量ＢTを表す(3.2.6)式、及び周波数領域の計測で検出される振幅Ａの減衰量Ｂfを表す(3.3.3.3)式等は、微分形で表され、かつ、それらの式は相似である。
【００９２】
これらの応答（検出光の所定パラメータ）をＹ、減衰量をＢ、各ボクセルの平均光路長をＺiとおけば、上記の各種応答は一般式、
【数７７】
で統一して表すことができる。
【００９３】
上記で、対象とする応答が光路長の異なる複数光子からなる場合、平均光路長Ｚはμaiの関数となる。また、対象とする応答が等しい光路長ｌをもつ光子からなる場合、平均光路長は光路分布（μa＝０）で決まるＺi＝ｌiとなり、このｌiはμaiに依存しない。以上から明らかなように、(4.1)式は不均一系に対するＭＢＬの一般式を表す。すなわち、この(4.1)式は、不均一媒体における減衰量が、各ボクセルの吸収係数と平均光路長のみによって決まることを表わしている。従って、不均一系に対するＭＢＬはこの事実を記述したものと考えてよい。
【００９４】
なお、μa＝０のときの応答をＳ（μs）とおけば、上記の(4.1)式から、
【数７８】
が得られる。但し、Ｗiは各ボクセルの重み関数である。
【００９５】
ここで、不均一媒体の平均吸収係数に近い値を基準にして、媒体各部の吸収係数分布を定量するＡＶＭによる光ＣＴ画像再構成との関係を説明する。この場合、不均一媒体の平均吸収係数に近い値μavが既知であるとし、ボクセルｉの吸収係数をμav＋Δμaiとする。また、3.2.3節と同様にして、吸収係数の偏差Δμaiの範囲内で、Ｚi（μai）がμaiのみの関数であると仮定すれば、(3.2.18)式と同様にして、
【数７９】
が得られる。但し、σiv 2は分散を表す。また、(4.4)式を行列で表示すれば、
【数８０】
となる。但し、重み関数Ｗiは、
【数８１】
と表される。
【００９６】
なお、一様な吸収がある媒体（仮想媒体）のボクセルｉの平均光路長及び分散は、モンテカルロシミュレーションなどで計算することができる。また、先に述べたようなＺi（μai）の全ボクセルの吸収に対する依存性が問題になる場合には、次に示す(4.7)式に戻って再計算するような繰り返し計算が必要である。
【００９７】
【数８２】
（第１実施形態）
図２〜図６を参照して本発明の好適な一実施形態である第１実施形態について説明する。図２には、散乱吸収体である被計測媒体１０の内部に分布する吸収成分の濃度分布を計測する内部特性分布計測装置（光ＣＴ画像計測装置）１が示してある。先ず、内部特性分布計測装置１の構成について説明する。
【００９８】
図２に示す装置は、光入射ファイバー２１を備えており、光入射ファイバー２１の先端が被計測媒体１０の表面の異なる位置ｕj（ｊ＝１〜ｐ）に移動可能となっている。そして、光入射ファイバー２１には波長選択器２２を介して光源２３が光学的に接続されており、光源２３から発せられた光は波長選択器２２で波長選択され、光入射ファイバー２１を通して異なる位置ｕj（ｊ＝１〜ｐ）から順次被計測媒体１０中に入射される。光源２３には、発光ダイオード、レーザーダイオード、Ｈｅ−Ｎｅレーザー等種々のものが使用できる。また、光源２３は、パルス光や方形波光、又はそれらの変調光を発生するものでもよい。本実施形態において使用する光源２３は単波長の光を発生するものであってもよいが、２波長以上の光を発生可能なものであってもよい。
【００９９】
また、図２に示す装置は、複数の光検出ファイバー２５を備えており、光検出ファイバー２５の先端が被計測媒体１０の表面の複数の位置ｖk（ｋ＝１〜ｑ）にそれぞれ設置されている。そして、光検出ファイバー２５には光検出器２６がそれぞれ光学的に接続されており、被計測媒体１０中で散乱されつつ透過した光は光検出ファイバー２５を介して光検出器２６に導かれ、光検出器２６で受光信号を検出信号（電気信号）に変換して増幅し、それぞれに対応する検出信号が出力される。光検出器２６は光電子増倍管のほか、光電管、フォトダイオード、アバランシェフォトダイオード、ＰＩＮフォトダイオード等、あらゆる種類の光検出器を使用することができる。光検出器２６の選択に際しては、使用される測定光の波長の光が検出できる分光感度特性をもっていれば良い。また、光信号が微弱であるときは高感度あるいは高利得の光検出器を使用することが好ましい。
【０１００】
なお、被計測媒体１０の表面における光入射ファイバー２１による光入射面及び光検出ファイバー２５による光検出面以外の場所は、光を吸収あるいは遮光する構造にしておくことが望ましい。また、被計測媒体１０の内部を拡散伝搬した光が複数の波長の光を含む場合には、光検出器２６と光検出ファイバー２５との間に波長選択フィルタ（図示せず）を適宜配置してもよい。
【０１０１】
光源２３及び光検出器２６にはＣＰＵ（制御及び演算処理部）３０が電気的に接続されており、光入射に同期した光検出のタイミング等がＣＰＵ３０によって制御されると共に、光検出器２６から出力された検出信号はＣＰＵ３０に導かれる。また、複数の波長を有する測定光を使用する場合は、入射光の波長もＣＰＵ３０によって制御される。具体的な手法としては、異なる波長の光を時分割で入射させて使用する手法と、異なる波長の光を同時に含む光を使用する手法とがある。具体的な波長選択手段としては、ミラーを用いた光ビーム切り換え器、フィルターを用いた波長切り換え器、光スイッチを用いた光切り換え器等がある。
【０１０２】
上記の光入射ファイバー２１、波長選択器２２、光源２３及びＣＰＵ３０が本発明にかかる光入射手段を構成し、上記の光検出ファイバー２５、光検出器２６及びＣＰＵ３０が本発明にかかる光検出手段を構成する。
【０１０３】
図２に示す内部特性分布計測装置１は、更に、オペレーティングシステム（ＯＳ）４１及び後で詳述する内部特性分布計測プログラム４２が記憶された第１記憶装置４０と、各種ファイルが記憶された第２記憶装置５０と、得られた内部特性分布を示す光ＣＴ画像データを記憶する画像メモリ６１と、作業用データを一時的に記憶する作業用メモリ６２と、データの入力を受け付けるキーボード７１及びマウス７２を備える入力装置７０と、得られたデータを出力するディスプレイ８１及びプリンタ８２を備える出力装置８０とを備えており、これらも電気的に接続されているＣＰＵ３０によって制御される。なお、上記の記憶装置及びメモリは、コンピュータの内部メモリ（ハードディスク）であっても、フレキシブルディスクであってもよい。
【０１０４】
第２記憶装置５０は、基本光路長分布データファイル５１と、測定データファイル５２と、入力データファイル５３と、吸収係数データファイル５４と、吸収成分濃度データファイル５５とを備えている。
【０１０５】
この内、基本光路長分布データファイル５１には、前もって用意された基本光路長分布（種々の計測に使う種々の重み関数演算の元になる共通の光路長分布）が格納されている。このような基本光路長分布は、既知の手法、例えばモンテカルロシミュレーションや光拡散方程式を利用して計算することができ、例えば（17）J. Haselgrove, J. Leigh, C. Yee, N-G, Wang, M. Maris and B. Chance: Proc. SPIE, Vol. 1431, 30-41（1991）；（18）J. C. Schotland, J. C. Haselgrove and J. S. Leigh: Appl. Opt. 32, 448-453（1993）；（19）Y. Tsuchiya, K. Ohta and T. Urakami: Jpn. J. Appl. Opt. 34, 2495-2501（1995）；（20）S. R. Arrige: Appl. Opt. 34, 7395-7409（1995）に記載されている。
【０１０６】
また、入力装置７０を用いて予め設定された計測条件や既知値が入力され、このような入力データが入力データファイル５３に格納されている。このような入力データとしては、予め任意に設定されたボクセル（体積素）の数、形状及び大きさ、被計測媒体の形状、光入射位置、光検出位置、散乱係数、平均吸収係数、計測に用いる光の波長、計測の種類（時間分解積分計測、時間分解ゲート積分計測、位相変調計測など）、計測対象となる吸収成分の所定の波長における消光比などがある。
【０１０７】
更に、測定データファイル５２には、内部特性分布計測プログラム４２の実行過程において光検出器２６から発せられた検出信号に基づいて求められる所定パラメータの測定値が光入射位置と光検出位置との組み合わせ毎に記憶される。また、吸収係数データファイル５４及び吸収成分濃度データファイル５５には、内部特性分布計測プログラム４２の実行によって得られた吸収係数データ及び吸収成分濃度データが記憶される。
【０１０８】
上記のＣＰＵ３０、第１記憶装置４０及び第２記憶装置５０が、本発明にかかる測定値取得手段と、基準値設定手段と、推定値算出手段と、重み関数演算手段と、吸収係数偏差算出手段と、吸収係数絶対値算出手段と、濃度算出手段とを構成し、上記の出力装置８０が画像表示手段を構成する。これらの本発明にかかる諸手段については、図３及び図４に示す本発明の方法の一実施形態のフローチャート（図２に示す内部特性分布計測プログラム４２の処理を示すフローチャート）に基づいて以下に詳細に説明する。
【０１０９】
１）図３に示すフローチャートにおいては、先ず、光源２３で生成した光（光ビーム）を光入射ファイバー２１を介して被計測媒体１０の異なる光入射位置ｕj（ｊ＝１〜ｐ）に順次入射し（Ｓ１０１）、被計測媒体内部で散乱されつつ透過した各光（各光入射位置ｕjから入射された光の一部）をそれぞれ複数の光検出位置ｖk（ｋ＝１〜ｑ）に設置した光検出ファイバー２５を介して光検出器２６で検出する（Ｓ１０２）。
【０１１０】
そして、各光入射位置ｕjから入射された光に対して、複数の光検出位置ｖk（ｋ＝１〜ｑ）で検出された光に対応する光検出信号が光検出器２６から発せられ、ＣＰＵ３０において検出光の所定パラメータの測定値Ｙnに変換される。このようにして、光入射位置ｕj（ｊ＝１〜ｐ）と光検出位置ｖk（ｋ＝１〜ｑ）との組み合わせＰn毎に検出光の所定パラメータの測定値Ｙnが得られ、測定データファイル５２にｎセット（ｎ＝ｐ×ｑ）の測定データセットとして記憶される（Ｓ１０３）。この際、ｎセットのデータに対する光入射・検出位置の組み合わせＰnはそれぞれ異なる。
【０１１１】
本発明にかかる所定パラメータの測定値Ｙnとしては、測定光の測定対象物内部での散乱及び吸収に関係する所定パラメータの測定値が好ましく、例えば検出光の光量、位相差（又は位相遅れ）、振幅、時間分解波形等のパラメータの測定値が好適に用いられる。また、ＣＰＵ３０において、光源２３からの光の発生に同期した信号を利用して光検出信号に対する所定時間域（ｔ1→ｔ2）での積分演算を行えば時間分解ゲート積分計測が可能となり、他方、積分範囲を０→∞とすれば通常の時間分解積分計測が可能となる。なお、パルス光等を利用する場合には、この同期信号を省略することができる。また、ＣＰＵ３０において、測定値を平均化フィルタリングや最小二乗フィッティング等を利用して修正してもよい。
【０１１２】
２）次に、被計測媒体１０の吸収係数の基準値μav及び輸送散乱係数の基準値μ'sを設定する（Ｓ１０４）。なお、吸収係数の基準値μavは、仮想媒体に対する真の平均吸収係数に等しくする必要はなく、μav＝０と近似してもよいが、一般的にはμavとして真の平均吸収係数に近い値を用いる方が高い計測精度が得られる傾向にある。従って、吸収係数の基準値μavとして、被計測媒体１０を巨視的に見たとき、すなわち被計測媒体１０が全体的に一様な吸収係数を有するとみなした場合の平均的な吸収係数（近似値）を採用することが好ましい。また同様に、輸送散乱係数の基準値μ'sとしても、被計測媒体１０が全体的に一様な散乱係数を有するとみなした場合の平均的な輸送散乱係数（近似値）を採用することが好ましい。
【０１１３】
このような吸収係数μavと輸送散乱係数μ'sは、前記の複数の光検出位置ｖkで検出された光検出信号（又は所定パラメータの測定値Ｙn）から求めることができる。例えば、前記光検出信号がインパルスレスポンス（光検出信号の時間波形に対して十分に短いと見なせるパルス光入射に対する光検出信号）である場合には、その時間波形を光拡散方程式にフィッティングすることで、所定の光入射・検出位置の組み合わせＰnに対応する媒体内部の吸収係数μavと輸送散乱係数μ'sを求めることができる。従って、各々の光入射・検出位置の組み合わせＰnに対して求めた前記の吸収係数μavの平均値及び輸送散乱係数μ'sの平均値を採用すればよい。さらに、これらの値は、時間分解ゲート積分計測や位相変調計測に利用することができる。また、このような吸収係数μav及び輸送散乱係数μ'sの取得方法に関しては、例えば（21）R. Berg, S. Andersson-Engels, O. Jarlman and S. Svanberg: Proc. SPIE, Vol. 1431, 110-119（1991）；（22）M. Miwa, Y. Ueda and B. Chance: Proc. SPIE, Vol. 2389, 142-149（1995）；（23）R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli and G. Valentini, Phy. Med. Biol. 42, 1971-1979（1997）に記載されている。
【０１１４】
また、上記基準値を設定する際に、被計測媒体１０を巨視的に見たときの平均的な他の光学定数を更に求めてもよく、このような光学定数としては、被計測媒体の屈折率ｎe、散乱係数μs、散乱角の余弦の平均値ｇが挙げられる。なお、被計測媒体の屈折率ｎeは、通常は水の屈折率で近似してよい。
【０１１５】
３）次に、被計測媒体１０が全体的に一様な吸収係数の基準値μavを有するとみなして、吸収係数の基準値μavに基づいて、光入射位置と光検出位置との組み合わせＰn毎に前記の所定パラメータの推定値Ｙavnを求める（Ｓ１０５）。
【０１１６】
すなわち、例えば、上記の平均的な光学定数（吸収係数の基準値μav、輸送散乱係数の基準値μ's等）を与えて光拡散方程式を解けば、各光入射・検出位置の組み合わせＰnにおける前記の所定パラメータ（光検出信号）の推定値Ｙavnが求められる。また、任意の光学定数（例えば、吸収係数）を与えたモンテカルロ計算の結果から、光学定数を変えたときの結果を推定する方法が開発されている（例えば、A. Kienle and M. S. Patterson: Phys. Med. Biol. 41, 2221-2227 (1996)に記載されている）。従って、このようなモンテカルロ計算の結果から、前記吸収係数の基準値μavを用いて各光入射・検出位置の組み合わせＰnにおける前記の所定パラメータ（光検出信号）の推定値Ｙavnを計算することができる。
【０１１７】
４）他方、被計測媒体１０が全体的に一様な吸収係数の基準値μavを有するとみなして、吸収係数の基準値μavに基づいて、複数のボクセルに分割された被計測媒体１０の各ボクセルｉにおける平均光路長Ｚi（μav）を求める（Ｓ１０６）。
【０１１８】
なお、本発明では、個々の計測に際してその状況に合致した重み関数を後述の方法で求める必要がある。この場合、計測のリアルタイム性を考慮して、重み関数をその場で（計測時に）迅速に求める必要があるため、種々の計測に使う種々の重み関数演算の元になる共通の基本光路長分布を前もって用意しておき、計測時に状況に応じてその基本光路長分布を変換することが好ましく、それによって必要な重み関数を迅速に求めることが可能となる傾向にある。
【０１１９】
従って、基本光路長分布データファイル５１からモンテカルロ計算等によって予め求められている前記基本光路長分布を読み出し、前記の吸収係数の基準値μav及び輸送散乱係数の基準値μ'sに基づいて、これらの基準値を有する仮想媒体に対する検出光の光路長分布を求めることができる。これが基本光路長分布の変換であり、例えば（24）A. Kienle and M. Patterson: Phys. Med. Biol. 41, 2221-2227（1996）に示されている方法を利用することができる。
【０１２０】
なお、ここで求められた仮想媒体に対する検出光の光路長分布Ｚi（μav）に基づいて、前記の各光入射・検出位置の組み合わせＰnにおける所定パラメータ（光検出信号）の推定値Ｙavnを計算することもできる。このような場合における処理の流れを図３中に鎖線Ａで示す。
【０１２１】
５）次に、上記のように吸収係数の基準値μavを用いて得られた各ボクセルｉにおける平均光路長Ｚi（μav）に基づいて、各ボクセルｉにおける重み関数Ｗiをマイクロ・ベア・ランバート則に則して(4.6)式にしたがって求める（Ｓ１０７）と共に、前記の各光入射・検出位置の組み合わせＰnにおける所定パラメータの測定値Ｙn及び推定値Ｙavnと前記重み関数Ｗiとに基づいて、各ボクセルｉにおける吸収係数の基準値μavに対する吸収係数の偏差Δμaiを(4.5)式にしたがって算出する（Ｓ１０８）。
【０１２２】
なお、求めるべき重み関数Ｗiは(4.6)式に示したもので、
【数８３】
である。ここで、重み関数Ｗiはそのボクセルのもつ吸収係数偏差Δμai＝μai−μav（ｉ＝１〜Ｎ）の関数になっていることに注意する必要がある。従って、このような重み関数Ｗiは、以下のような繰り返し演算によって吸収係数偏差Δμaiと共に求めることが好ましい。すなわち、Ｓ１０７及びＳ１０８の演算処理（図３中に鎖線Ｂで示す）は、通常、繰り返し演算されることとなる。図４は、このような繰り返し演算によって重み関数Ｗi及び吸収係数偏差Δμaiを求める処理の好適な一実施形態を示すフローチャートである。以下、図４に示すフローチャートについて詳述する。
【０１２３】
先ず、次数Ｋの初期値をＫ＝１、吸収係数偏差Δμaiの初期値をΔμaiK＝０と設定する（Ｓ２０１）。次いで、各ボクセルｉにおける平均光路長Ｚi（μav）及び吸収係数偏差ΔμaiKに基づいて(4.6)式にしたがって重み関数ＷiKを求める（Ｓ２０２）。なお、(4.6)式中のσiv 2は分散を表し、例えば（3.2.10）式や（3.3.8）式を計算して求められる。
【０１２４】
次に、前記の所定パラメータの測定値Ｙn及び推定値Ｙavnと前記重み関数ＷiKとに基づいて、各ボクセルｉにおける吸収係数の偏差Δμai(K+1)を(4.5)式にしたがって算出する（Ｓ２０３）。なお、前記ＹnとＹavnとの比の自然対数値［ln（Ｙavn／Ｙn）］が(4.5)式中のΔＹnに相当する。
【０１２５】
そして、Ｓ２０３によって求められた吸収係数の偏差Δμai(K+1)とＳ２０２で用いた吸収係数偏差ΔμaiKとの残差Ｅを算出し（Ｓ２０４）、その値が所定値（err）以下であるか否かを判断する（Ｓ２０５）。
【０１２６】
その結果、上記残差Ｅが所定値以下である場合はその時の重み関数ＷiKを重み関数Ｗiとして、吸収係数偏差Δμai(K+1)を吸収係数偏差Δμaiとしてそれぞれ出力する（Ｓ２０６）。他方、上記残差Ｅが所定値以下でない場合は、Ｓ２０２に用いる吸収係数偏差ΔμaiKを前記Ｓ２０３によって求められた吸収係数偏差Δμai(K+1)に置き換えて（次数Ｋに１を付加して：Ｓ２０７）、重み関数ＷiKの演算（Ｓ２０２）、吸収係数偏差Δμai(K+1)の演算（Ｓ２０３）及び前記残差Ｅの算出（Ｓ２０４）を再度行なう。そして、得られた残差Ｅが所定値（err）以下であるか否かを判断し（Ｓ２０５）、残差Ｅが所定値以下となるまで上記のＳ２０７→Ｓ２０２→Ｓ２０３→Ｓ２０４→Ｓ２０５を繰り返す。なお、上記のＳ２０１を初期値設定手段、Ｓ２０２を重み関数演算手段、Ｓ２０３を吸収係数偏差演算手段、Ｓ２０４、Ｓ２０５及びＳ２０６を評価手段、Ｓ２０７を吸収係数偏差修正手段でそれぞれ行なってもよい。
【０１２７】
６）以上によって、吸収係数偏差の分布（各ボクセルｉにおける吸収係数の偏差）Δμai＝μai−μav（ｉ＝１〜Ｎ）が得られる。そして、Δμaiは先に述べた吸収係数の基準値（好ましくは平均吸収係数）μavからの偏差を示す。従って、各ボクセルｉにおける吸収係数の偏差Δμai及び吸収係数の基準値μavを用いて、次式：
μai＝Δμai＋μav（ｉ＝１〜Ｎ）
にしたがって各ボクセルｉにおける吸収係数の絶対値μaiを算出することができ（Ｓ１０９）、吸収係数データファイル５４に記憶される。図５に、吸収係数の基準値μavと偏差Δμaiと絶対値μaiとの関係を示す。そして、このようにして求められた各ボクセルｉにおける吸収係数の絶対値μaiに基づいて被計測媒体内部における吸収係数絶対値に関する分布が求められ、その分布を示す光ＣＴ画像が出力装置８０に表示される（Ｓ１１０）。
【０１２８】
７）次に、吸収成分の濃度演算について説明する。ここでは、被計測媒体が一種類の吸収成分を含む場合を考える。すると、その吸収成分のボクセルｉにおける濃度偏差ΔＣiは、ベア・ランバート則から導かれる次式：
εΔＣi＝μai−μav (C.2)
で与えられる。但し、εは吸収成分の単位濃度当たりの吸収係数（又は消光係数）であり、分光光度計で測定することができる。さらに、前記吸収成分のボクセルｉにおける濃度Ｃiは、ベア・ランバート則から導かれる次式：
εＣi＝μai (C.3)
で与えられる。従って、吸収成分の既知の吸収係数を用いて上記の各ボクセルｉにおける吸収係数の絶対値μaiから各ボクセルｉにおける吸収成分の濃度Ｃiを算出することができ（Ｓ１１１）、吸収成分濃度データファイル５５に記憶される。そして、このようにして求められた各ボクセルｉにおける吸収成分の濃度Ｃiに基づいて被計測媒体内部における吸収成分濃度に関する分布が求められ、その分布を示す光ＣＴ画像が出力装置８０に表示される（Ｓ１１２）。
【０１２９】
（第２実施形態）
上記の第１実施形態の方法において２種類以上の波長の光を用いれば、２種類以上の吸収成分の濃度偏差や濃度絶対値が測定できる。本実施形態では、波長λ1とλ2の２つの光で計測する２波長分光計測法について説明する。すなわち、波長λ1及びλ2の光に対応する各吸収成分の単位濃度当たりの吸収係数（又は消光係数）ε1及びε2が既知であれば、前記(C.2)式及び(C.3)式に対して２本の式が成立する。従って、これらの式を連立させることによって、２種類の吸収成分の濃度偏差や濃度絶対値を測定することができる。
【０１３０】
より具体的には、散乱吸収体が少なくとも２つの吸収成分（例えば酸素化型及び脱酸素化型ヘモグロビン）を含有する場合は、それらの吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有する測定光を使用し、各波長を有する測定光のそれぞれに関して測定値Ｙn、推定値Ｙavn、基準値μav及び重み関数Ｗiを求め、それらに基づいて各波長を有する測定光のそれぞれに関して吸収係数偏差Δμai、吸収係数絶対値μaiを求めることによって、各吸収成分の濃度Ｃiの分布が求められる。
【０１３１】
以下に、上記２波長分光法を利用したヘモグロビンの濃度の計測について説明する。哺乳類の脳における吸収成分の主なものは、水、チトクローム（cytochrom)、酸素化型及び脱酸素化型ヘモグロビンである。近赤外線領域での水及びチトクロームの吸収は、酸素化型及び脱酸素化型ヘモグロビンに対して、ほぼ無視することができる程度に少ない。また、酸素化型及び脱酸素化型ヘモグロビンは、図６に示すように、吸収スペクトルが異なる。さらに、頭蓋骨は、近赤外線に対して、散乱体と考えてよい。
【０１３２】
いま前節までに述べた方法で波長λ1 とλ2 の２種の波長の光に対して吸収係数μa1とμa2が求められたとすれば、ランバート・ベール（Lambert-Beer) 則によって、次式が成立する。
【０１３３】
μa1＝εHb,1〔Ｈｂ〕＋εHbO,1 〔ＨｂＯ〕
μa2＝εHb,2〔Ｈｂ〕＋εHbO,2 〔ＨｂＯ〕
但し、
εHb,1；脱酸素化型ヘモグロビンの波長λ1 に対するモル吸収係数〔ｍｍ-1・Ｍ-1〕
εHbO,1 ；酸素化型ヘモグロビンの波長λ1 に対するモル吸収係数〔ｍｍ-1・Ｍ-1〕
εHb,2；脱酸素化型ヘモグロビンの波長λ2 に対するモル吸収係数〔ｍｍ-1・Ｍ-1〕
εHbO,2 ；酸素化型ヘモグロビンの波長λ2 に対するモル吸収係数〔ｍｍ-1・Ｍ-1〕
〔Ｈｂ〕：脱酸素化型ヘモグロビンのモル濃度〔Ｍ〕
〔ＨｂＯ〕：酸素化型ヘモグロビンのモル濃度〔Ｍ〕
である。
【０１３４】
従って、既知のパラメータεHb,1、εHbO,1 、εHb,2、εHbO,2 及び計測値から演算されたμa1とμa2から、脱酸素化型ヘモグロビンのモル濃度〔Ｈｂ〕、及び酸素化型ヘモグロビンのモル濃度〔ＨｂＯ〕を求めることができる。
【０１３５】
また、上記に対してチトクロームを考慮する場合のように、吸収スペクトルが既知である３成分のそれぞれの濃度の定量は、３波長以上の光を使用すればよい。一般的には、吸収スペクトルが既知であるｎ個の成分の濃度の定量計測は、ｎ個又は（ｎ＋１）個の波長に対する吸収係数の計測値から、上記と同様にして求めることができる。
【０１３６】
さらに、飽和度Ｙは、
Ｙ＝〔ＨｂＯ〕／（〔Ｈｂ〕＋〔ＨｂＯ〕）
であるから、
μa1／μa2＝〔εHb,1＋Ｙ（εHbO,1 −εHb,1）〕÷〔εHb,2＋Ｙ（εHbO,2 −εHb,2）〕
を用いて、既知のパラメータεHb,1、εHbO,1 、εHb,2、εHbO,2 と計測値から演算されたμa1及びμa2とから、飽和度Ｙが容易に算出される。
【０１３７】
以上の方法では、本発明によって各波長の光に対する吸収係数μa1とμa2が高精度に求められるため、各濃度も高精度に求められる。なお、酸素化型及び脱酸素化型ヘモグロビンに対して吸収が同一値になる波長（８００ｎｍ、isosbestic wavelength)を使用すれば上記の式はさらに簡単になる。
【０１３８】
（第３実施形態）
本実施形態は本発明を位相変調計測に応用する例を示す。この場合、前記の第１実施形態における入射光を振幅変調光、所定パラメータを検出光の振幅Ａと位相遅れφ、検出光の光路長分布Ｚi（μav）をｃ（媒体内の光速）倍の群遅延
【数８４】
の分布、σiv 2を分散にそれぞれ変更する以外は第１実施形態と同様にして測定値Ｙn、推定値Ｙavn、基準値μav及び重み関数Ｗiを求め、それらに基づいて吸収係数偏差Δμai、吸収係数絶対値μaiを求めることによって、吸収成分濃度Ｃiの分布が求められる。
【０１３９】
なお、この場合は、求めるべき重み関数Ｗiは次式：
【数８５】
に示したものであり、σf 2は分布の分散を表す。また、第１実施形態における基本光路長分布に代えて、ｃ倍の群遅延分布が用いられる。
【０１４０】
以上、本発明の好適な実施形態について説明したが、本発明は勿論上記実施形態に限定されるものではない。
【０１４１】
すなわち、上記実施形態においては吸収係数の基準値及び輸送散乱係数の基準値を光ＣＴ装置自身で得られたデータから求めているが、これらの基準値を別の装置で求めるようにしてもよい。この場合の利点は、例えば光ＣＴ装置で得られるデータはＣＷ（連続光）で測定し、上記基準値を求める装置でのみパルス光や変調光を用いればよくなるため、光ＣＴ装置のシステムの構成が簡単になる。なお、別の装置で上記基準値を求める手法は位相変調法や時間分解分光法であってもよい。
【０１４２】
また、上記実施形態においては光入射位置を走査させていたが、光入射位置と同期して光検出位置を走査させてもよい。また、被計測媒体の周囲に複数の光入射ファイバー及び複数の光検出ファイバーを配置しておき、光入射に使用するファイバーをそれぞれ順次変えていくことによって光入射位置を移動させてもよい。
【０１４３】
更に、以上説明した本発明の装置及び方法は、一般的な断層像を得る光ＣＴだけでなく、３次元媒体内部の吸収成分を３次元計測する装置及び方法に利用できることは明らかである。また、本発明の装置及び方法は、頭部の皮膚、頭蓋骨、灰白質、白質のように、多数の層状の構造をもつ媒体の計測に利用できることも明白である。なお、この場合は、それぞれの層をそれぞれのボクセルに対応させればよい。
【０１４４】
【発明の効果】
以上説明したように、本発明によれば、新規な重み関数を使用して個々の計測に際して各ボクセルにおける重み関数がマイクロ・ベア・ランバート則に基づいて直接的に求められるため、計測状況に応じた適切な重み関数に基づいて吸収係数の偏差が算出され、従って近似法の採用に伴う誤差の発生が防止される。そのため、本発明によれば、生体などの不均一系の散乱吸収体について適用した場合であっても各ボクセルにおける吸収係数の偏差を正確に求めることが可能となり、このような吸収係数偏差に基づいて高精度の内部特性分布を得て光ＣＴ画像化することが可能となる。
【０１４５】
このように、本発明によって再入射不可能な様々な外形をもつ種々の散乱吸収体の内部吸収成分の分布を高精度かつ迅速に計測することが可能となり、例えば、光マンモグラフィー、頭部光ＣＴ、全身用光ＣＴ、臨床モニター、診断や解析、手術や治療への応用も可能となる。
【図面の簡単な説明】
【図１】不均一媒体内の光子移動に関するモデルを示す模式図である。
【図２】本発明の散乱吸収体の内部特性分布の計測装置の一実施形態を示す模式図である。
【図３】本発明の散乱吸収体の内部特性分布の計測方法の一実施形態を示すフローチャートである。
【図４】本発明にかかる重み関数の演算方法の一実施形態を示すフローチャートである。
【図５】各ボクセルにおける吸収係数分布（吸収係数の基準値と吸収係数の偏差との関係）を示すグラフである。
【図６】ヘモグロビン及びミオグロビンの吸収スペクトルを示すグラフである。
【符号の説明】
１…内部特性分布計測装置、１０…被計測媒体、２１…光入射ファイバー、２２…波長選択器、２３…光源、２５…光検出ファイバー、２６…光検出器、３０…ＣＰＵ、４０…第１記憶装置、４１…ＯＳ、４２…内部特性分布計測プログラム、５０…第２記憶装置、５１…基本光路長分布データファイル、５２…測定データファイル、５３…入力データファイル、５４…吸収係数データファイル、５５…吸収成分濃度データファイル、６１…画像メモリ、６２…作業用メモリ、７０…入力装置、７１…キーボード、７２…マウス、８０…出力装置、８１…ディスプレイ、８２…プリンタ。
【発明の属する技術分野】
本発明は、散乱吸収体の内部特性分布を計測するための方法及びそのための装置に関する。より詳細には、本発明は、光入射位置及び光検出位置を測定対象物の表面に沿って移動させて内部の情報を得る装置に適用可能な内部特性分布の計測方法並びにそのための装置に関する。
【０００２】
【従来の技術】
光ＣＴ(computer tomography)は、生体内部の光学特性あるいは吸収成分の濃度分布を計測する手法又は装置を意味し、光を生体組織に入射したときにその中を移動して検出される光（信号光）を利用する。
【０００３】
現在までに直進光、準直進光、及び散乱光を利用する３種の光ＣＴ技術が報告されている。このうち、直進光を利用する方法は光の利用率が極めて悪く、極小さい媒体にしか適用できない。例えば、近赤外線を用いる場合、標準的な生体組織では１cmの厚さに対して透過光は約1/105倍に減衰する。これに比べて、散乱光を利用する方法は、媒体から出射する全ての光を利用するため、信号対雑音比が向上するので、より大きい媒体への適用が期待される。準直進光を利用する方法は、前記二つの方法の中間にある。そして、実用的な光ＣＴは、生体に対する最大許容光入射パワー、計測感度、所要計測時間などの制限から散乱光を利用する必要があるが、以下の技術上の問題に起因して上記光ＣＴは未だ実用化に至っていない。
【０００４】
第１の問題は、散乱吸収体内部の光あるいはフォトンの振る舞いを十分な精度で記述する方法が開発されていないことである。散乱吸収体内部を移動するフォトンの振る舞いの解析には、従来から、輸送方程式を近似したもの、あるいは輸送理論に拡散近似を適用した光拡散方程式が用いられてきた。しかし、拡散近似は媒体中のフォトンの平均自由光路長に対して十分に大きな媒体においてのみ成立するため、比較的小さい媒体、内部の複雑な形状の組織、及び複雑な形状の媒体を取り扱うことができない。また、拡散近似は等方散乱を前提としているため、非等方な散乱特性をもつ実際の生体組織の計測に適用すると、非等方散乱に起因する無視できない誤差が発生する。さらに、拡散方程式は、解析的あるいは数値的（有限要素法など）手法のいずれを用いても、前もって境界条件を設定しないとその解が求められない。つまり、計測に先立って、光入射・検出する各部の境界条件、つまり媒体形状や界面での反射特性を設定する必要がある。また、これらが個体差などによって変化すれば、それに応じて境界条件を変えて、計算をやりなおす必要がある。従って、輸送方程式の近似式や光拡散方程式から導出される信号光と散乱吸収体の光学特性との関係を利用する光ＣＴは、精度及び使い勝手に大きな問題がある。
【０００５】
また、上記とは別に、輸送方程式の近似式や光拡散方程式に摂動理論を適用して信号光と散乱吸収体の光学特性との関係を導出し、この関係を用いて光ＣＴ画像を再構成する方法がある。しかし、この方法は、非線形効果（２次以上の項）の取り扱いが極めて複雑になる。この際、２次以上の項を含めた計算は、理論的にはコンピュータで計算することが可能であるが、現在の最速のコンピュータを用いても演算時間が膨大なものになり、実用化は不可能である。このため、常套的には２次以上の項を無視している。そのため、この方法においては、複数個の比較的強い吸収領域が存在する媒体の光ＣＴ画像を再構成する際などに吸収領域の相互作用を無視することができなくなり、これに起因する大きな誤差が発生する。
【０００６】
第２の問題は、従来の光ＣＴが、狭義の重み関数、すなわち平均光路長あるいはそれに相当する位相遅れを利用していることである。このため、吸収係数に依存して変化する検出光の平均光路長の取り扱いが極めて複雑になる。そこで通常は近似を用いるが、近似を用いると誤差が増大するという重大な問題がある。このような狭義の重み関数を利用した方法に関しては、例えば以下の文献に記載されている。(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)。
【０００７】
なお、光ＣＴの画像再構成に際しては、信号光が散乱媒体のどの部分を通過してきたか、つまり媒体内の信号光の光路分布を知ることが最も重要である。この観点から、信号光の光路分布をランダムウォーク理論などで導出する方法もあるが、上記の問題は解決されていない。
【０００８】
以上のように、従来の光ＣＴにおいては十分な精度の再構成画像が得られず、空間解像度、画像ひずみ、定量性、測定感度、所要計測時間などに大きな問題がある。
【０００９】
本発明者は、以上のような状況を打破するため、次のことが重要であると考えて一連の研究を推進してきた。すなわち、光ＣＴを実現させるために重要なことは、強い散乱媒体である生体組織の中を移動する光の振る舞いを明らかにし、検出される信号光と吸収成分を含む散乱媒体（散乱吸収体）の光学特性との関係を明らかにし、前記信号光と前記関係とを利用して光ＣＴ画像を再構成するアルゴリズムを開発するということである。
【００１０】
そして本発明者は、▲１▼マイクロ・べア・ランバート則(Microscopic Beer-Lambert Law、以下「ＭＢＬ」という)に基づくモデルを提案すると共に、▲２▼散乱吸収体の光学特性と信号光との関係を表す解析式を導出し、例えば以下の文献に掲載している。(13) Y. Tsuchiya and T. Urakami: Jpn. J. Appl. Phys. 34, L79-81 (1995); (14) Y. Tsuchiya and T. Urakami: Jpn. J. Appl. Phys. 35, 4848-4851 (1996); (15) Y. Tsuchiya and T. Urakami: Optics Communications, 144, 269-280 (1997); (16) H. Zhang, M. Miwa, Y. Yamashita and Y. Tsuchiya: Ext. Abstr. Optics Japan '97, 30aA08 (1997)。このＭＢＬに基づく方法は原理的に、媒体形状、境界条件及び散乱の影響を受けないという大きな特長があり、非等方散乱媒体や小さい媒体への適用も可能である。この結果、散乱吸収体内の光子移動が明らかになり、この理論に基づいて散乱媒体の吸収係数あるいは吸収成分濃度を計測することが可能になった。
【００１１】
【発明が解決しようとする課題】
しかしながら、本発明者は、上記従来の方法はいずれも上記従来の狭義の重み関数を利用していたため、不均一系への適用が困難であり、生体などの散乱吸収体内の吸収成分についての計測精度が未だ充分なものではないことを見出した。
【００１２】
本発明は、上記従来の課題に鑑みてなされたものであり、生体などの不均一系の散乱吸収体について適用した場合であっても吸収成分についてより高い精度の計測が可能な方法及び装置を提供することを目的とする。
【００１３】
【課題を解決するための手段】
本発明者は、上記目的を達成すべく鋭意研究した結果、上記本発明者による知見をさらに発展させて、▲３▼上記マイクロ・べア・ランバート則を不均一系に適用し、▲４▼不均一散乱媒体の光学特性と信号光との関係を表す解析式を導出し、▲５▼前記信号光と不均一散乱吸収体の光学特性との関係から従来とは異なる定義の重み関数を導出し、▲６▼この重み関数を用いて光ＣＴ画像を再構成するためのアルゴリズム及び装置を提供することによって上記目的が達成されることを見出し、本発明に到達した。
【００１４】
すなわち、本発明の散乱吸収体の内部特性分布の計測方法は、
ヒトを除く散乱吸収体である被計測媒体中に一以上の光入射位置から順次光を入射する光入射ステップと、
前記被計測媒体内部を透過した光を複数の光検出位置で検出する光検出ステップと、
検出された各光に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記光の所定パラメータの測定値を取得する測定値取得ステップと、
前記被計測媒体の吸収係数の基準値を設定する基準値設定ステップと、
前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記パラメータの推定値を求める推定値算出ステップと、
前記吸収係数の基準値を用いて、複数のボクセルに分割された前記被計測媒体の各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて求める重み関数演算ステップと、
前記パラメータの測定値、前記パラメータの推定値及び前記重み関数に基づいて、前記各ボクセルにおける前記吸収係数の基準値に対する吸収係数の偏差を算出する吸収係数偏差算出ステップと、
前記吸収係数の基準値及び前記吸収係数の偏差に基づいて前記各ボクセルにおける吸収係数の絶対値を算出して前記被計測媒体における吸収係数の絶対値分布を求める吸収係数絶対値算出ステップと、
を含むことを特徴とする方法である。
【００１５】
また、本発明の散乱吸収体の内部特性分布の計測装置は、
散乱吸収体である被計測媒体中に一以上の光入射位置から順次光を入射する光入射手段と、
前記被計測媒体内部を透過した光を複数の光検出位置で検出する光検出手段と、
検出された各光に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記光の所定パラメータの測定値を取得する測定値取得手段と、
前記被計測媒体の吸収係数の基準値を設定する基準値設定手段と、
前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記パラメータの推定値を求める推定値算出手段と、
前記吸収係数の基準値を用いて、複数のボクセルに分割された前記被計測媒体の各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて求める重み関数演算手段と、
前記パラメータの測定値、前記パラメータの推定値及び前記重み関数に基づいて、前記各ボクセルにおける前記吸収係数の基準値に対する吸収係数の偏差を算出する吸収係数偏差算出手段と、
前記吸収係数の基準値及び前記吸収係数の偏差に基づいて前記各ボクセルにおける吸収係数の絶対値を算出して前記被計測媒体における吸収係数の絶対値分布を求める吸収係数絶対値算出手段と、
を備えることを特徴とする装置である。
【００１６】
本発明の方法及び装置においては、個々の計測に際して本発明によって初めて開示される各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて直接的に求め、その重み関数、所定パラメータの測定値及び前記パラメータの推定値に基づいて各ボクセルにおける吸収係数の偏差が算出される。この際、重み関数は計測状況とは無関係な一種の式で表わされる。このように本発明においては、個々の計測状況が変化しても適切な重み関数に基づいて吸収係数の偏差が算出されるため、近似法の採用に伴う誤差の発生が防止される。そのため、本発明の方法及び装置においては、生体などの不均一系の散乱吸収体について適用した場合であっても、各ボクセルにおける吸収係数の偏差が正確に求まることとなり、このような吸収係数偏差に基づいて得られる内部特性分布（光ＣＴ画像）の測定精度が向上する。
【００１７】
上記本発明の方法及び装置において採用される重み関数として、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける平均光路長と光路長分布の分散との関数が挙げられる。
【００１８】
この場合、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける平均光路長を求める平均光路長取得ステップ（平均光路長取得手段）を更に含んでいることが好ましく、前記重み関数演算ステップ（重み関数演算手段）において、下記の式：
【数９】
に基づいて前記重み関数を求めることが好ましい。
【００１９】
また、上記本発明の方法及び装置において採用される他の重み関数として、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける所定時間領域内の平均光路長と光路長分布の分散との関数が挙げられる。
【００２０】
この場合、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける所定時間領域内の平均光路長を求める平均光路長取得ステップ（平均光路長取得手段）を更に含んでいることが好ましく、前記重み関数演算ステップ（重み関数演算手段）において、下記の式：
【数１０】
に基づいて前記重み関数を求めることが好ましい。
【００２１】
更に、上記本発明の方法及び装置において採用される更に他の重み関数として、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける群遅延と分布の分散との関数が挙げられる。
【００２２】
この場合、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける群遅延を求める群遅延取得ステップ（群遅延取得手段）を更に含んでおり、前記重み関数演算ステップ（重み関数演算手段）において、下記の式：
【数１１】
【数１２】
に基づいて前記重み関数を求めることが好ましい。
【００２３】
上記本発明の方法及び装置に好適な重み関数を採用した場合、吸収係数に依存して変化する平均光路長又は群遅延のみならず分布の分散をも考慮されて重み関数が演算されるため、このような重み関数を用いて得られる内部特性分布の測定精度がより向上する傾向にある。
【００２４】
本発明の方法及び装置は、更に、前記吸収係数の絶対値を用いて前記各ボクセルにおける吸収成分の濃度を算出して前記被計測媒体における吸収成分の濃度分布を求める濃度算出ステップ（濃度算出手段）を含んでもよい。このような方法及び装置によれば、前述のように正確に求められた吸収係数偏差に基づいて吸収成分の濃度が求められるため、高精度の濃度分布が得られる。
【００２５】
上記本発明の方法及び装置を少なくとも２つの吸収成分を含有している被計測媒体に適用する場合には、前記光入射ステップ（光入射手段）において前記被計測媒体中に入射される光が、前記吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有していることが好ましい。この場合、前記光検出ステップ（光検出手段）において前記少なくとも２つの波長を有する光をそれぞれ検出し、前記測定値取得ステップ（測定値取得手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記測定値を取得し、前記基準値設定ステップ（基準値設定手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記基準値を設定し、前記推定値算出ステップ（推定値算出手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記推定値を求め、前記重み関数演算ステップ（重み関数演算手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記重み関数を求め、前記吸収係数偏差算出ステップ（吸収係数偏差算出手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収係数の偏差を算出し、前記吸収係数絶対値算出ステップ（吸収係数絶対値算出手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収係数の絶対値を算出し、前記濃度算出ステップ（濃度算出手段）において前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収成分の濃度を算出することによって、前記被計測媒体における前記各吸収成分の濃度分布が高精度で求められる。
【００２６】
上述の本発明の方法及び装置は、前記の求められた分布に基づいて、前記被計測媒体における該分布を示す光ＣＴ画像を表示する画像表示ステップ（画像表示手段）を更に含んでもよい。このような本発明の方法及び装置によれば、前述のように精度良く求められた内部特性分布を画像化して高精度な光ＣＴ画像として表示することが可能である。
【００２７】
【発明の実施の形態】
以下、図面を参照しつつ本発明の原理及び好適な実施形態について詳細に説明する。尚、図面中、同一又は相当部分には同一符号を付することとする。なお、散乱されつつ進行する光に関しては３次元座標を用いて考える必要があるが、以下では説明を簡単にするために場合により２次元座標を用いて説明する。先ず、本発明の原理について説明する。
【００２８】
［本発明の原理］
本発明は、マイクロ・ベア・ランバート則に基づいて導出される新規な重み関数を用いる。すなわち、本発明にかかわる重み関数は、従来の重み関数、つまり各ボクセル（voxel、体積素）内の光路長や平均光路長、光子滞在時間、あるいはそれらに対応する位相遅延量などとは異なる。本発明で採用する重み関数は、媒体形状、境界条件、光入射・検出位置及びそれらの間の距離、媒体の散乱係数などを繰り込んだ形で表される。また、この重み関数は散乱媒体に対して一意的に定まる。本発明では上記のような重み関数を用いて光ＣＴの画像再構成を行うため、後述するように、従来の問題点である誤差が大きく低減され、高精度の光ＣＴ画像を得ることができる。
【００２９】
１．光子移動モデル
先ず、本発明の基礎である散乱媒体内の光子移動のモデルを説明する。
【００３０】
Ａ）散乱媒体は再入射不可能な任意の３次元形状をもつ。従って、媒体から出た光子が再度その媒体に入射することはない。この制限内で、散乱媒体と外部媒体との境界条件、例えば形状や屈折率のマッチング条件などは任意である。
【００３１】
Ｂ）散乱媒体内で、光子と媒体との非線形作用や、光子と光子との相互作用はない。従って、媒体の光学定数は光子密度に依存しない。
【００３２】
Ｃ）マクロに見たとき、屈折特性は一様（均一）である。しかし、吸収係数分布は一様でない（不均一である）。以下、このような媒体を便宜的に不均一媒体あるいは不均一系とよぶ。従って、以下で議論する散乱媒体内の光速度は一定値ｃになる。
【００３３】
Ｄ）上記不均一散乱媒体に対して、光入射・検出位置ｐ（ｕ，ｖ）（表面位置ｕに光子を入射し、他の表面位置ｖで散乱光子を検出する）を定めたとき、光路分布（μa＝０）が一意的に決まる。この光路分布（μa＝０）は、散乱媒体に吸収がないと仮定したときに光子が取り得る光路の分布あるいは光路長分布である。なお、このように吸収がないと仮定した散乱媒体を「仮想媒体」という。
【００３４】
Ｅ）上記の光路分布（μa＝０）は、散乱媒体に吸収がないと仮定した仮想媒体に対して定義したから、散乱媒体の吸収あるいは吸収分布には影響されない。また、この光路分布（μa＝０）は、散乱媒体の散乱特性、媒体形状、境界条件、及び光入射・検出位置ｐ（ｕ，ｖ）に依存する。そして、これらの影響は光路分布（μa＝０）の中に繰り込まれている。つまり、散乱等の影響は、吸収がない仮想媒体に対する光路分布（μa＝０）に反映されるため、以下では光路分布（μa＝０）やμa＝０である仮想媒体に対する応答を、単に散乱等の影響とよぶこともある。
【００３５】
Ｆ）また、前記Ｂ）項と関連するが、上記の光路分布（μa＝０）は、複数あるいは多数光子（光子集合）入射や多数光子の移動に対して不変である。
【００３６】
また、媒体内の光子移動と減衰に関しては、以下の関係：
Ｇ）吸収係数がμaである均一媒体内を光路長ｌ(エル)＝ｃｔだけ移動する単一光子に対する減衰量(attenuation)はμaｌとなるという関係（以下、この関係を「ＭＢＬ第１基本則」とよぶ）；及び
Ｈ）散乱媒体内の多数の光子（光子集合）に対する減衰量は、それぞれの光子に対する減衰量の重ね合せ(superposition)で表されるという関係；
が成立する。
【００３７】
２．不均一系内の光子移動に関する基礎事項
最初に不均一系内の光子移動に関する基礎事項を説明する。
【００３８】
2.1 不均一媒体内の光子移動
吸収（又は吸収成分）が不均一に分布する３次元の不均一媒体全体をＮ個のボクセルに分割して番号ｉを付け、それぞれのボクセルｉについて非等方散乱係数μs、散乱角の余弦の平均値ｇ、吸収係数μaiを定義する（図１参照）。光子入射・検出位置はｐ（ｕ，ｖ）である。ここで、μsとμaiは不均一媒体固有の値で、かつμsは定数、またμaiは位置の関数である。この際、ボクセルｉは任意の大きさと形状にすることができ、大きい体積をもつボクセルではそのボクセル内の平均吸収係数や平均光路長を考えればよい。また、空間的に離れた２個以上のボクセルを合わせて、１個のボクセルと考えてもよい。
【００３９】
いま、単一光子を繰り返して入射するとして、ｍ番目の入射光子（m-th光子）に対して時間ｔに検出される出力ｈm(t)を考える。そして、このm-th光子のボクセルｉにおける滞在（飛行）時間をｔim、また対応する光路長をｌim（＝ｃｔim）とする。この場合、検出される光子の光路長ｌとｌimとの関係は、ｌ＝ｃｔ＝Σｌimである。ここで、添字ｍはm-th光子に関する量を表す。上記の系でm-th光子の光路、つまり光路分布（μa＝０）が決まると、ｌimが一意的に決まる。また、ボクセル番号ｉを決めると、そのボクセルｉの吸収係数μaiが一意的に決まる。この場合、吸収係数は任意に選べるから、ボクセルの光路長と吸収係数は独立である。
【００４０】
時間ｔに検出されるm-th光子の減衰量（減光量ともよばれる）Ｂmは、
【数１３】
【００４１】
この(2.1.1)式は、不均一系におけるＭＢＬ第１基本則、つまり「吸収係数がμaiである散乱体の中をジグザグ光路に沿って長さｌ＝Σｌimだけ移動した光子の減衰量は、媒体の散乱特性とは無関係な値Σ（μaiｌim）になる」ことを表す。この関係は、光子移動の際に吸収係数が変化する位置でＭＢＬを接続したものである。また、この関係は、ジグザグ光路が交差（光子が２回以上同一場所を通過）する場合にも成立する。
【００４２】
ここで、Σμai＝０のときのm-th光子に対する応答をｓm（μs，ｔ）とおけば、m-th光子に対する応答ｈm(t)は、
【数１４】
【数１５】
【００４３】
すると、(2.1.3)式から、不均一系内のm-th光子に対して、
【数１６】
【００４４】
不均一媒体に多数の光子（光子集合）を入射したときの応答を構成する光子で、光路長がｌとなる光子を考える。この場合、光路長がｌとなる光路は媒体内に無限といえるほど多数存在する。そして、多数の光子（光子集合）に対する減衰量は、単一光子に対する減衰量の重ね合せ(superposition)で表される。また、通常の応答、例えばインパルス応答には、光路長の異なる光子が含まれている。この場合にも、減衰量は単一光子に対する減衰量の重ね合せで表される。また、先に定義した光路分布（μa＝０）は、多数の光子入射に対して不変である。従って、多数光子入射に対する不均一系の各種応答の減衰量は、単一光子に対する減衰量の重ね合せ、つまり線形問題として解析することができる。
【００４５】
2.2 不均一媒体の光路分布（μa＝０）
再入射不可能な任意の不均一媒体に対して光入射・検出位置ｐ（ｕ，ｖ）を定めたとき、光子が取り得る光路分布は一意的に決まる。この光路分布（μa＝０）は、その不均一媒体に吸収がない（μa＝０）と仮定したとき（仮想媒体、μa＝０）の光路分布で定義されるから、不均一媒体内の吸収やその分布には影響されない。しかし、光路分布（μa＝０）は、不均一媒体の光学特性、媒体形状や境界条件、光入射・検出位置ｐ（ｕ，ｖ）、及び光路長に依存する。
【００４６】
不均一系に対するインパルス応答、インパルス応答の時間分解ゲート積分信号、及び周波数領域の応答を構成する光子集合に対して、それぞれの光路分布（μa＝０）を定義することができる。また、光路分布（μa＝０）が決まれば、各ボクセルの平均光路長（μa＝０）が一意的に決まる。
【００４７】
以上のような光路分布（μa＝０）は、モンテカルロシミュレーションなどによって計算することができる。さらに、モンテカルロシミュレーションでは、光路分布（μa＝０）を基準にして、任意の一様な吸収（例えばμa＝μav）がある場合の検出光の光路分布や平均光路長を計算することもできる。
【００４８】
2.3 検出光の光路分布
インパルス応答には、種々の光路長をもつ光子が含まれる。媒体に吸収があるとき、個々の光子の減衰量はＢm＝Σ（μaiｌim）になるが、光路長が異なる光子集合の平均光路長は、吸収係数の関数になる。従って、応答の時間波形は、吸収に依存して変化することになる。つまり、吸収があると光子移動の途中で光子が消滅し、その消滅の程度が光路長あるいは光路分布（μa＝０）によって異なるため、その分だけ検出光に対する平均光路長や光路分布が変化する。この結果、検出光の光路分布は、通常、吸収に依存することになり、前出の光路分布（μa＝０）に等しくならない。
【００４９】
2.4 重み関数
不均一媒体内をｌだけ移動した単一光子の減衰量は、前述したようにΣ（μaiｌim）で表される。そして、不均一媒体における多数光子の減衰量は、単一光子の減衰量の重ね合せで表される。ここで、多光子で構成される不均一媒体の各種応答の減衰量をボクセルｉにおける寄与係数Ｗiとボクセルｉの吸収係数μaiとの積で表すことを考える。そして、前記の寄与係数Ｗiを「重み関数(Weight Function)」とよぶ。
【００５０】
前出の(2.1.4)式と(2.1.5)式の場合には、光路長ｌ＝ｃｔ＝Σｌimが吸収係数μaiの関数にならないから、光路分布（μa＝０）と重み関数は等しい。しかし、インパルス応答の時間積分値のように平均光路長が吸収係数の関数になる場合、後述するように、重み関数は光路分布（μa＝０）や検出光の光路分布に等しくない。また、一般に重み関数Ｗiは吸収に依存する。
【００５１】
３．不均一系における各種の応答関数
前述のことを踏まえて、不均一系における各種の応答を求めて、不均一系に対するＭＢＬを導出する。この際、不均一媒体に対する光入射・検出位置をｐ（ｕ，ｖ）とする。すると、時間領域や周波数領域の各種計測で得られる各種応答に対して、それぞれの光路分布（μa＝０）が一意的に決まる。そして、各種応答の減衰量は、各種応答を構成する単一光子に対する減衰量の重ね合せになる。なお、光路分布（μa＝０）は吸収分布に対して不変であるから、計測系の設定、つまり光入射・検出位置ｐ（ｕ，ｖ）や計測方式が同一であれば、被計測媒体が不均一媒体であっても均一媒体であっても、同一の光路分布（μa＝０）になる。ＭＢＬに基づく記述・解析方法は、不均一系の応答の場合でも、μa＝０としたときの応答（つまり散乱や境界条件の影響）と吸収による減衰量とを分離した形で記述することができるという点で優れている。
【００５２】
3.1 不均一系のインパルス応答
不均一系のインパルス応答ｈ(t)は、種々の光路長をもつ多数の検出光子で構成される。この場合、観測時間ｔを決めるとｌ＝ｃｔとして光路長ｌが決まる。いま、光路長がｌである検出光子の数Ｍが十分に大きいとする。この場合、全ての光子の光路長はｌm＝ｌであるから、Ｍ個の光子に対する平均光路長もｌ、つまり、
【数１７】
【数１８】
【００５３】
また、インパルス応答を構成する光子（複数）の中で光路長がｌ（エル）となる光子集合に対する減衰量Ｂhとボクセルｉにおける減衰量Ｂihは、
【数１９】
【数２０】
【数２１】
【００５４】
以上から、不均一系のインパルス応答ｈ(t)は、
【数２２】
【００５５】
なお、上記(3.1.7)式から、不均一系のインパルス応答ｈ(t)に対する重み関数は、先に定義した光路分布（μa＝０）で決まるボクセルｉ（ｉ＝１〜Ｎ）の光路長ｌi（μa＝０）に等しいことがわかる。従って、時間分解計測による光ＣＴ画像再構成アルゴリズムは、比較的簡単になる。
【００５６】
以下では、上記で求めた不均一系のインパルス応答ｈ(t)を構成する光子集合の光路分布（μa＝０）、つまり前記で定義した各ボクセルの光路長ｌi（μa＝０）、あるいはそれらの時間分布（異なる値をもつそれぞれのｌ（μa＝０）に対するｌi（μa＝０）であり、時間分解型の光路分布（μa＝０）とよぶ）を用いて各種の応答とそれらに対する重み関数を導出する。
【００５７】
3.2 不均一系のインパルス応答の時間分解ゲート積分信号
ここでは、先ず不均一系のインパルス応答の時間分解ゲート積分信号に対する解析式を導出し、減衰量や平均光路長の吸収依存性を考慮した重み関数を導出する。なお、時間域の積分範囲［ｔ1，ｔ2］を［０，∞］にすると、定常（ＣＷ）光に対する応答を表すことになる。
【００５８】
3.2.1 時間分解ゲート積分信号の解析式
インパルス応答の時間分解ゲート積分信号ＩTは、(3.1.8)式を積分して
【数２３】
【００５９】
この式から、
【数２４】
【００６０】
すると、(3.2.2)式から
【数２５】
【００６１】
次に、(3.2.2)式をμaiで積分して、時間分解ゲート積分信号ＩTに対する基本式を求めると、
【数２６】
【数２７】
【００６２】
ここで、先に述べたＬi（μa）の全ボクセルの吸収係数に対する依存性が問題になる。つまり、この依存性があると、(3.2.5)式をこれ以上簡単にすることができない。この事実は重要であり、光ＣＴ画像再構成アルゴリズムを複雑にしている。
【００６３】
ところで一般には、各ボクセル間の吸収係数μaiの差、つまり平均吸収係数からの偏差はそれほど大きくない場合が多い。この場合、注目しているボクセルｉの平均光路長Ｌi（μa）に対する他部分（又は他のボクセル）の吸収係数偏差の影響は小さい。換言すれば、各ボクセル間の吸収係数の差が小さい場合、注目しているボクセルｉの平均光路長Ｌi（μa）は、全ボクセルに対する吸収係数の平均値で大まかに決まり、他部分の吸収係数の平均吸収係数に対する偏差に依存するＬi（μa）の変化は無視できる程度に小さく、Ｌi（μa）はボクセルｉの吸収係数μaiのみの関数であると見なしてよい。この場合、ボクセルｉの吸収係数の平均吸収係数との差（吸収係数偏差）に対応するＬi（μa）の差分は、線形近似して記述することができる。従って、以下では、Ｌi（μa）はボクセルｉの吸収係数μaiのみの関数であると仮定する。つまり、ボクセルｉの平均光路長Ｌi（μa）を考える際に、ボクセルｉの吸収係数μai以外のパラメータは一定（不変）であると仮定する。なお、２つ以上のボクセルの吸収係数が同時に変化する場合、これらの２つのボクセルを合わせたボクセルを１つの新しいボクセルと考える。このことは、「１つの測定値から複数ボクセルの吸収係数の変化を定量することはできない」ことと等価であり、理論上の矛盾はない。
【００６４】
以上の線形近似を適用すると、ボクセルｉにおける減衰量ＢiTを定義すると、
【数２８】
【００６５】
次に、時間分解ゲート積分信号ＩTに対する重み関数Ｗiについて考える。先ず、平均値定理を用いて(3.2.4)式と(3.2.6)式の減衰量ＢTとＢiTを変形すると、
【数２９】
【００６６】
なお、前記Ｌi（μa）の全ボクセルの吸収に対する依存性が問題になる場合には、(3.2.6)式の代わりに(3.2.5)式を使用する必要がある。
【００６７】
以上のことも、積分範囲［ｔ1，ｔ2］を［０，∞］にすることによって、通常の時間分解積分信号、つまり定常（ＣＷ）光に対する応答に適用することができる。
【００６９】
3.2.2 平均光路長の吸収依存性
各ボクセルの平均光路長は、そのボクセルの吸収に依存する。また、重み関数もそのボクセルの吸収に依存する。本節ではボクセルｉの吸収係数μaiが（μai＋ｈ）に変化（Δμai＝ｈ）したときの減衰量や平均光路長の変化を解析して、それらの変化量を表す解析式を導出する。この結果、光ＣＴの画像再構成に際して問題となる平均光路長と重み関数の差、及びそれらの吸収依存性が明らかになる。
【００７０】
先ず、(3.2.7.2)式に示したボクセルｉの減衰量ＢiTをμaiの周りにテイラー展開すると、
【数３４】
【数３５】
【数３６】
【数３７】
【００７１】
ここで、吸収係数がμaiから（μai＋ｈ）に変化したときの減衰量の差ΔＢiTを(3.2.7.2)式から求めると、
【数３８】
【数３９】
【００７２】
他方、平均光路長Ｌi（μai）をμaiの周りにテイラー展開すると、
【数４０】
【数４１】
【００７３】
以上を踏まえてボクセルｉの吸収係数μaiが（μai＋ｈ）に変化したとき、(3.2.4)式を用いて変化前後の時間分解ゲート積分信号lnＩTの差ΔＩを求めると、
【数４２】
【００７４】
また、複数ボクセルｉ（１〜Ｎの中の複数の値）の吸収係数μaiが（μai＋ｈi）に変化したとき、(3.2.4)式を用いて変化前後の時間分解ゲート積分信号lnＩTの差ΔＩを求めると、
【数４３】
【００７５】
以上の場合、吸収係数μaiを被計測媒体の平均吸収係数あるいはそれに近い値に選んでおくと、精度が向上する。つまり、さきに述べたＬi（μa）の全ボクセルの吸収に対する依存性の問題が大きく緩和される。
【００７６】
なお、各ボクセルの吸収係数μaiの平均吸収係数に対する偏差が極端に大きい場合で、他のボクセルの吸収係数偏差に対するＬi（μa）の依存性が問題になる場合には、(3.2.5)式に戻って再計算をするような繰り返し計算が必要になる。
【００７７】
以上のことは、積分範囲［ｔ1，ｔ2］を［０，∞］にすると、通常の時間分解積分信号、つまり定常（ＣＷ）光に対する応答に適用することができる。
【００７８】
3.2.3 アベレージ・バリュー・メソッドへの応用
光ＣＴは、不均一媒体中の吸収成分の濃度分布を計測しようとするものである。この場合、先に述べた２種類の平均光路長の吸収依存性、つまりＬi（μa）の全ボクセルの吸収に対する依存性、及びＬi（μa）のボクセルｉ自身の吸収に対する依存性が問題を複雑にしている。本発明者は、これらの吸収依存性に基づく誤差を低減する一つの方法として、アベレージ・バリュー・メソッド（Average Value Method、以下「ＡＶＭ」という）を開発した。このＡＶＭは、計測すべき不均一媒体に対して近似的な平均吸収係数を推定あるいは計測して、この値を基準にして、媒体各部の吸収係数の偏差を定量する。この場合、上記で推定あるいは計測する平均吸収係数は、真値に等しい必要はない。この方法は、平均光路長の吸収係数依存性の問題を大きく緩和するもので、前記で推定あるいは計測した平均吸収係数が真値に近いほど、再構成光ＣＴ画像の精度が向上する。
【００７９】
いま、被計測媒体の平均吸収係数に近い値の吸収係数μavをもつ仮想媒体（μa＝μav）を仮定する。すると計測対象である不均一媒体の吸収係数はμai＝μav＋Δμai（ｉ＝１〜Ｎ）となる。時間分解ゲート積分計測を考えると、仮想媒体に対する計測値Ｉ（μav）と実際の媒体に対する計測値 Ｉ（μai＝μav＋Δμai）（ｉ＝１〜Ｎ）との関係は(3.2.4)式及び(3.2.15)式、あるいは(3.2.17)式から、
【数４４】
【００８０】
前記(3.2.18)式を行列式で表せば、
【数４５】
【数４６】
【００８１】
3.3 不均一系の周波数領域の応答
周波数領域の応答は、前出の(3.1.8)式に示すインパルス応答ｈ(t)をフーリエ変換して、
【数４７】
【数４８】
【００８２】
すると、後で詳述するように、
【数４９】
【数５０】
【数５１】
【００８３】
従って、(3.3.3.1)式ないし(3.3.3.4)式の関係式から、前出の(3.2.4)式の導出と同様にして、
【数５２】
【００８４】
以上の(3.3.6.1)式ないし(3.3.6.4)式は、前節に記載した時間分解ゲート積分信号に対する式と相似であるから、前節と同様にして、これらの式を変形することができる。そこで以下では、(3.3.6.3)式を例として取り上げる。先ず、前節と同様に、周波数領域の応答に対する重み関数を考える、先ず、平均値定理を用いると(3.3.6.3)式の減衰項（右辺第２項）は、
【数５３】
【数５４】
【数５５】
【数５６】
【００８５】
さらに、ボクセルｉの減衰量Ｂif及び平均群遅延(3.3.7.4)の吸収係数依存性は、前節に述べた減衰量ＢiT及び平均光路長Ｌi（μx）と同様になる。つまり、群遅延(3.3.7.4)の吸収依存性を考えると、(3.3.3.4)式及び(3.3.5.2)式より、
【数５７】
【数５８】
【数５９】
【数６０】
【００８６】
以上から、前述の説明（3.2.3節）と同様にして、
【数６１】
【数６２】
【００８７】
［(3.3.3.1)式〜(3.3.3.4)式の導出］
先ず、ｔ＝Σｔiであることを考慮して、(3.3.2.1)式及び(3.3.2.2)式から、
【数６３】
【数６４】
【数６５】
【数６６】
【数６７】
【００８８】
さらに、(3.3.2.3)式及び上式から、
【数６８】
【数６９】
【数７０】
【数７１】
【００８９】
また、同様にして、(3.3.2.4)式、(3.3.3.1)式及び(3.3.3.2)式から
【数７２】
【数７３】
【数７４】
【００９０】
なお、上記のとき、
【数７５】
【数７６】
【００９１】
４．不均一系に関する一般式
前章までの説明から明らかなように、不均一系におけるＭＢＬは、微分形で表すと有効である。つまり、不均一系に対するインパルス応答ｈ(t)の減衰量Ｂhを表す(3.1.6)式、インパルス応答の時間分解ゲート積分信号ＩT（μs、μa）の減衰量ＢTを表す(3.2.6)式、及び周波数領域の計測で検出される振幅Ａの減衰量Ｂfを表す(3.3.3.3)式等は、微分形で表され、かつ、それらの式は相似である。
【００９２】
これらの応答（検出光の所定パラメータ）をＹ、減衰量をＢ、各ボクセルの平均光路長をＺiとおけば、上記の各種応答は一般式、
【数７７】
【００９３】
上記で、対象とする応答が光路長の異なる複数光子からなる場合、平均光路長Ｚはμaiの関数となる。また、対象とする応答が等しい光路長ｌをもつ光子からなる場合、平均光路長は光路分布（μa＝０）で決まるＺi＝ｌiとなり、このｌiはμaiに依存しない。以上から明らかなように、(4.1)式は不均一系に対するＭＢＬの一般式を表す。すなわち、この(4.1)式は、不均一媒体における減衰量が、各ボクセルの吸収係数と平均光路長のみによって決まることを表わしている。従って、不均一系に対するＭＢＬはこの事実を記述したものと考えてよい。
【００９４】
なお、μa＝０のときの応答をＳ（μs）とおけば、上記の(4.1)式から、
【数７８】
【００９５】
ここで、不均一媒体の平均吸収係数に近い値を基準にして、媒体各部の吸収係数分布を定量するＡＶＭによる光ＣＴ画像再構成との関係を説明する。この場合、不均一媒体の平均吸収係数に近い値μavが既知であるとし、ボクセルｉの吸収係数をμav＋Δμaiとする。また、3.2.3節と同様にして、吸収係数の偏差Δμaiの範囲内で、Ｚi（μai）がμaiのみの関数であると仮定すれば、(3.2.18)式と同様にして、
【数７９】
【数８０】
【数８１】
【００９６】
なお、一様な吸収がある媒体（仮想媒体）のボクセルｉの平均光路長及び分散は、モンテカルロシミュレーションなどで計算することができる。また、先に述べたようなＺi（μai）の全ボクセルの吸収に対する依存性が問題になる場合には、次に示す(4.7)式に戻って再計算するような繰り返し計算が必要である。
【００９７】
【数８２】
図２〜図６を参照して本発明の好適な一実施形態である第１実施形態について説明する。図２には、散乱吸収体である被計測媒体１０の内部に分布する吸収成分の濃度分布を計測する内部特性分布計測装置（光ＣＴ画像計測装置）１が示してある。先ず、内部特性分布計測装置１の構成について説明する。
【００９８】
図２に示す装置は、光入射ファイバー２１を備えており、光入射ファイバー２１の先端が被計測媒体１０の表面の異なる位置ｕj（ｊ＝１〜ｐ）に移動可能となっている。そして、光入射ファイバー２１には波長選択器２２を介して光源２３が光学的に接続されており、光源２３から発せられた光は波長選択器２２で波長選択され、光入射ファイバー２１を通して異なる位置ｕj（ｊ＝１〜ｐ）から順次被計測媒体１０中に入射される。光源２３には、発光ダイオード、レーザーダイオード、Ｈｅ−Ｎｅレーザー等種々のものが使用できる。また、光源２３は、パルス光や方形波光、又はそれらの変調光を発生するものでもよい。本実施形態において使用する光源２３は単波長の光を発生するものであってもよいが、２波長以上の光を発生可能なものであってもよい。
【００９９】
また、図２に示す装置は、複数の光検出ファイバー２５を備えており、光検出ファイバー２５の先端が被計測媒体１０の表面の複数の位置ｖk（ｋ＝１〜ｑ）にそれぞれ設置されている。そして、光検出ファイバー２５には光検出器２６がそれぞれ光学的に接続されており、被計測媒体１０中で散乱されつつ透過した光は光検出ファイバー２５を介して光検出器２６に導かれ、光検出器２６で受光信号を検出信号（電気信号）に変換して増幅し、それぞれに対応する検出信号が出力される。光検出器２６は光電子増倍管のほか、光電管、フォトダイオード、アバランシェフォトダイオード、ＰＩＮフォトダイオード等、あらゆる種類の光検出器を使用することができる。光検出器２６の選択に際しては、使用される測定光の波長の光が検出できる分光感度特性をもっていれば良い。また、光信号が微弱であるときは高感度あるいは高利得の光検出器を使用することが好ましい。
【０１００】
なお、被計測媒体１０の表面における光入射ファイバー２１による光入射面及び光検出ファイバー２５による光検出面以外の場所は、光を吸収あるいは遮光する構造にしておくことが望ましい。また、被計測媒体１０の内部を拡散伝搬した光が複数の波長の光を含む場合には、光検出器２６と光検出ファイバー２５との間に波長選択フィルタ（図示せず）を適宜配置してもよい。
【０１０１】
光源２３及び光検出器２６にはＣＰＵ（制御及び演算処理部）３０が電気的に接続されており、光入射に同期した光検出のタイミング等がＣＰＵ３０によって制御されると共に、光検出器２６から出力された検出信号はＣＰＵ３０に導かれる。また、複数の波長を有する測定光を使用する場合は、入射光の波長もＣＰＵ３０によって制御される。具体的な手法としては、異なる波長の光を時分割で入射させて使用する手法と、異なる波長の光を同時に含む光を使用する手法とがある。具体的な波長選択手段としては、ミラーを用いた光ビーム切り換え器、フィルターを用いた波長切り換え器、光スイッチを用いた光切り換え器等がある。
【０１０２】
上記の光入射ファイバー２１、波長選択器２２、光源２３及びＣＰＵ３０が本発明にかかる光入射手段を構成し、上記の光検出ファイバー２５、光検出器２６及びＣＰＵ３０が本発明にかかる光検出手段を構成する。
【０１０３】
図２に示す内部特性分布計測装置１は、更に、オペレーティングシステム（ＯＳ）４１及び後で詳述する内部特性分布計測プログラム４２が記憶された第１記憶装置４０と、各種ファイルが記憶された第２記憶装置５０と、得られた内部特性分布を示す光ＣＴ画像データを記憶する画像メモリ６１と、作業用データを一時的に記憶する作業用メモリ６２と、データの入力を受け付けるキーボード７１及びマウス７２を備える入力装置７０と、得られたデータを出力するディスプレイ８１及びプリンタ８２を備える出力装置８０とを備えており、これらも電気的に接続されているＣＰＵ３０によって制御される。なお、上記の記憶装置及びメモリは、コンピュータの内部メモリ（ハードディスク）であっても、フレキシブルディスクであってもよい。
【０１０４】
第２記憶装置５０は、基本光路長分布データファイル５１と、測定データファイル５２と、入力データファイル５３と、吸収係数データファイル５４と、吸収成分濃度データファイル５５とを備えている。
【０１０５】
この内、基本光路長分布データファイル５１には、前もって用意された基本光路長分布（種々の計測に使う種々の重み関数演算の元になる共通の光路長分布）が格納されている。このような基本光路長分布は、既知の手法、例えばモンテカルロシミュレーションや光拡散方程式を利用して計算することができ、例えば（17）J. Haselgrove, J. Leigh, C. Yee, N-G, Wang, M. Maris and B. Chance: Proc. SPIE, Vol. 1431, 30-41（1991）；（18）J. C. Schotland, J. C. Haselgrove and J. S. Leigh: Appl. Opt. 32, 448-453（1993）；（19）Y. Tsuchiya, K. Ohta and T. Urakami: Jpn. J. Appl. Opt. 34, 2495-2501（1995）；（20）S. R. Arrige: Appl. Opt. 34, 7395-7409（1995）に記載されている。
【０１０６】
また、入力装置７０を用いて予め設定された計測条件や既知値が入力され、このような入力データが入力データファイル５３に格納されている。このような入力データとしては、予め任意に設定されたボクセル（体積素）の数、形状及び大きさ、被計測媒体の形状、光入射位置、光検出位置、散乱係数、平均吸収係数、計測に用いる光の波長、計測の種類（時間分解積分計測、時間分解ゲート積分計測、位相変調計測など）、計測対象となる吸収成分の所定の波長における消光比などがある。
【０１０７】
更に、測定データファイル５２には、内部特性分布計測プログラム４２の実行過程において光検出器２６から発せられた検出信号に基づいて求められる所定パラメータの測定値が光入射位置と光検出位置との組み合わせ毎に記憶される。また、吸収係数データファイル５４及び吸収成分濃度データファイル５５には、内部特性分布計測プログラム４２の実行によって得られた吸収係数データ及び吸収成分濃度データが記憶される。
【０１０８】
上記のＣＰＵ３０、第１記憶装置４０及び第２記憶装置５０が、本発明にかかる測定値取得手段と、基準値設定手段と、推定値算出手段と、重み関数演算手段と、吸収係数偏差算出手段と、吸収係数絶対値算出手段と、濃度算出手段とを構成し、上記の出力装置８０が画像表示手段を構成する。これらの本発明にかかる諸手段については、図３及び図４に示す本発明の方法の一実施形態のフローチャート（図２に示す内部特性分布計測プログラム４２の処理を示すフローチャート）に基づいて以下に詳細に説明する。
【０１０９】
１）図３に示すフローチャートにおいては、先ず、光源２３で生成した光（光ビーム）を光入射ファイバー２１を介して被計測媒体１０の異なる光入射位置ｕj（ｊ＝１〜ｐ）に順次入射し（Ｓ１０１）、被計測媒体内部で散乱されつつ透過した各光（各光入射位置ｕjから入射された光の一部）をそれぞれ複数の光検出位置ｖk（ｋ＝１〜ｑ）に設置した光検出ファイバー２５を介して光検出器２６で検出する（Ｓ１０２）。
【０１１０】
そして、各光入射位置ｕjから入射された光に対して、複数の光検出位置ｖk（ｋ＝１〜ｑ）で検出された光に対応する光検出信号が光検出器２６から発せられ、ＣＰＵ３０において検出光の所定パラメータの測定値Ｙnに変換される。このようにして、光入射位置ｕj（ｊ＝１〜ｐ）と光検出位置ｖk（ｋ＝１〜ｑ）との組み合わせＰn毎に検出光の所定パラメータの測定値Ｙnが得られ、測定データファイル５２にｎセット（ｎ＝ｐ×ｑ）の測定データセットとして記憶される（Ｓ１０３）。この際、ｎセットのデータに対する光入射・検出位置の組み合わせＰnはそれぞれ異なる。
【０１１１】
本発明にかかる所定パラメータの測定値Ｙnとしては、測定光の測定対象物内部での散乱及び吸収に関係する所定パラメータの測定値が好ましく、例えば検出光の光量、位相差（又は位相遅れ）、振幅、時間分解波形等のパラメータの測定値が好適に用いられる。また、ＣＰＵ３０において、光源２３からの光の発生に同期した信号を利用して光検出信号に対する所定時間域（ｔ1→ｔ2）での積分演算を行えば時間分解ゲート積分計測が可能となり、他方、積分範囲を０→∞とすれば通常の時間分解積分計測が可能となる。なお、パルス光等を利用する場合には、この同期信号を省略することができる。また、ＣＰＵ３０において、測定値を平均化フィルタリングや最小二乗フィッティング等を利用して修正してもよい。
【０１１２】
２）次に、被計測媒体１０の吸収係数の基準値μav及び輸送散乱係数の基準値μ'sを設定する（Ｓ１０４）。なお、吸収係数の基準値μavは、仮想媒体に対する真の平均吸収係数に等しくする必要はなく、μav＝０と近似してもよいが、一般的にはμavとして真の平均吸収係数に近い値を用いる方が高い計測精度が得られる傾向にある。従って、吸収係数の基準値μavとして、被計測媒体１０を巨視的に見たとき、すなわち被計測媒体１０が全体的に一様な吸収係数を有するとみなした場合の平均的な吸収係数（近似値）を採用することが好ましい。また同様に、輸送散乱係数の基準値μ'sとしても、被計測媒体１０が全体的に一様な散乱係数を有するとみなした場合の平均的な輸送散乱係数（近似値）を採用することが好ましい。
【０１１３】
このような吸収係数μavと輸送散乱係数μ'sは、前記の複数の光検出位置ｖkで検出された光検出信号（又は所定パラメータの測定値Ｙn）から求めることができる。例えば、前記光検出信号がインパルスレスポンス（光検出信号の時間波形に対して十分に短いと見なせるパルス光入射に対する光検出信号）である場合には、その時間波形を光拡散方程式にフィッティングすることで、所定の光入射・検出位置の組み合わせＰnに対応する媒体内部の吸収係数μavと輸送散乱係数μ'sを求めることができる。従って、各々の光入射・検出位置の組み合わせＰnに対して求めた前記の吸収係数μavの平均値及び輸送散乱係数μ'sの平均値を採用すればよい。さらに、これらの値は、時間分解ゲート積分計測や位相変調計測に利用することができる。また、このような吸収係数μav及び輸送散乱係数μ'sの取得方法に関しては、例えば（21）R. Berg, S. Andersson-Engels, O. Jarlman and S. Svanberg: Proc. SPIE, Vol. 1431, 110-119（1991）；（22）M. Miwa, Y. Ueda and B. Chance: Proc. SPIE, Vol. 2389, 142-149（1995）；（23）R. Cubeddu, A. Pifferi, P. Taroni, A. Torricelli and G. Valentini, Phy. Med. Biol. 42, 1971-1979（1997）に記載されている。
【０１１４】
また、上記基準値を設定する際に、被計測媒体１０を巨視的に見たときの平均的な他の光学定数を更に求めてもよく、このような光学定数としては、被計測媒体の屈折率ｎe、散乱係数μs、散乱角の余弦の平均値ｇが挙げられる。なお、被計測媒体の屈折率ｎeは、通常は水の屈折率で近似してよい。
【０１１５】
３）次に、被計測媒体１０が全体的に一様な吸収係数の基準値μavを有するとみなして、吸収係数の基準値μavに基づいて、光入射位置と光検出位置との組み合わせＰn毎に前記の所定パラメータの推定値Ｙavnを求める（Ｓ１０５）。
【０１１６】
すなわち、例えば、上記の平均的な光学定数（吸収係数の基準値μav、輸送散乱係数の基準値μ's等）を与えて光拡散方程式を解けば、各光入射・検出位置の組み合わせＰnにおける前記の所定パラメータ（光検出信号）の推定値Ｙavnが求められる。また、任意の光学定数（例えば、吸収係数）を与えたモンテカルロ計算の結果から、光学定数を変えたときの結果を推定する方法が開発されている（例えば、A. Kienle and M. S. Patterson: Phys. Med. Biol. 41, 2221-2227 (1996)に記載されている）。従って、このようなモンテカルロ計算の結果から、前記吸収係数の基準値μavを用いて各光入射・検出位置の組み合わせＰnにおける前記の所定パラメータ（光検出信号）の推定値Ｙavnを計算することができる。
【０１１７】
４）他方、被計測媒体１０が全体的に一様な吸収係数の基準値μavを有するとみなして、吸収係数の基準値μavに基づいて、複数のボクセルに分割された被計測媒体１０の各ボクセルｉにおける平均光路長Ｚi（μav）を求める（Ｓ１０６）。
【０１１８】
なお、本発明では、個々の計測に際してその状況に合致した重み関数を後述の方法で求める必要がある。この場合、計測のリアルタイム性を考慮して、重み関数をその場で（計測時に）迅速に求める必要があるため、種々の計測に使う種々の重み関数演算の元になる共通の基本光路長分布を前もって用意しておき、計測時に状況に応じてその基本光路長分布を変換することが好ましく、それによって必要な重み関数を迅速に求めることが可能となる傾向にある。
【０１１９】
従って、基本光路長分布データファイル５１からモンテカルロ計算等によって予め求められている前記基本光路長分布を読み出し、前記の吸収係数の基準値μav及び輸送散乱係数の基準値μ'sに基づいて、これらの基準値を有する仮想媒体に対する検出光の光路長分布を求めることができる。これが基本光路長分布の変換であり、例えば（24）A. Kienle and M. Patterson: Phys. Med. Biol. 41, 2221-2227（1996）に示されている方法を利用することができる。
【０１２０】
なお、ここで求められた仮想媒体に対する検出光の光路長分布Ｚi（μav）に基づいて、前記の各光入射・検出位置の組み合わせＰnにおける所定パラメータ（光検出信号）の推定値Ｙavnを計算することもできる。このような場合における処理の流れを図３中に鎖線Ａで示す。
【０１２１】
５）次に、上記のように吸収係数の基準値μavを用いて得られた各ボクセルｉにおける平均光路長Ｚi（μav）に基づいて、各ボクセルｉにおける重み関数Ｗiをマイクロ・ベア・ランバート則に則して(4.6)式にしたがって求める（Ｓ１０７）と共に、前記の各光入射・検出位置の組み合わせＰnにおける所定パラメータの測定値Ｙn及び推定値Ｙavnと前記重み関数Ｗiとに基づいて、各ボクセルｉにおける吸収係数の基準値μavに対する吸収係数の偏差Δμaiを(4.5)式にしたがって算出する（Ｓ１０８）。
【０１２２】
なお、求めるべき重み関数Ｗiは(4.6)式に示したもので、
【数８３】
【０１２３】
先ず、次数Ｋの初期値をＫ＝１、吸収係数偏差Δμaiの初期値をΔμaiK＝０と設定する（Ｓ２０１）。次いで、各ボクセルｉにおける平均光路長Ｚi（μav）及び吸収係数偏差ΔμaiKに基づいて(4.6)式にしたがって重み関数ＷiKを求める（Ｓ２０２）。なお、(4.6)式中のσiv 2は分散を表し、例えば（3.2.10）式や（3.3.8）式を計算して求められる。
【０１２４】
次に、前記の所定パラメータの測定値Ｙn及び推定値Ｙavnと前記重み関数ＷiKとに基づいて、各ボクセルｉにおける吸収係数の偏差Δμai(K+1)を(4.5)式にしたがって算出する（Ｓ２０３）。なお、前記ＹnとＹavnとの比の自然対数値［ln（Ｙavn／Ｙn）］が(4.5)式中のΔＹnに相当する。
【０１２５】
そして、Ｓ２０３によって求められた吸収係数の偏差Δμai(K+1)とＳ２０２で用いた吸収係数偏差ΔμaiKとの残差Ｅを算出し（Ｓ２０４）、その値が所定値（err）以下であるか否かを判断する（Ｓ２０５）。
【０１２６】
その結果、上記残差Ｅが所定値以下である場合はその時の重み関数ＷiKを重み関数Ｗiとして、吸収係数偏差Δμai(K+1)を吸収係数偏差Δμaiとしてそれぞれ出力する（Ｓ２０６）。他方、上記残差Ｅが所定値以下でない場合は、Ｓ２０２に用いる吸収係数偏差ΔμaiKを前記Ｓ２０３によって求められた吸収係数偏差Δμai(K+1)に置き換えて（次数Ｋに１を付加して：Ｓ２０７）、重み関数ＷiKの演算（Ｓ２０２）、吸収係数偏差Δμai(K+1)の演算（Ｓ２０３）及び前記残差Ｅの算出（Ｓ２０４）を再度行なう。そして、得られた残差Ｅが所定値（err）以下であるか否かを判断し（Ｓ２０５）、残差Ｅが所定値以下となるまで上記のＳ２０７→Ｓ２０２→Ｓ２０３→Ｓ２０４→Ｓ２０５を繰り返す。なお、上記のＳ２０１を初期値設定手段、Ｓ２０２を重み関数演算手段、Ｓ２０３を吸収係数偏差演算手段、Ｓ２０４、Ｓ２０５及びＳ２０６を評価手段、Ｓ２０７を吸収係数偏差修正手段でそれぞれ行なってもよい。
【０１２７】
６）以上によって、吸収係数偏差の分布（各ボクセルｉにおける吸収係数の偏差）Δμai＝μai−μav（ｉ＝１〜Ｎ）が得られる。そして、Δμaiは先に述べた吸収係数の基準値（好ましくは平均吸収係数）μavからの偏差を示す。従って、各ボクセルｉにおける吸収係数の偏差Δμai及び吸収係数の基準値μavを用いて、次式：
μai＝Δμai＋μav（ｉ＝１〜Ｎ）
にしたがって各ボクセルｉにおける吸収係数の絶対値μaiを算出することができ（Ｓ１０９）、吸収係数データファイル５４に記憶される。図５に、吸収係数の基準値μavと偏差Δμaiと絶対値μaiとの関係を示す。そして、このようにして求められた各ボクセルｉにおける吸収係数の絶対値μaiに基づいて被計測媒体内部における吸収係数絶対値に関する分布が求められ、その分布を示す光ＣＴ画像が出力装置８０に表示される（Ｓ１１０）。
【０１２８】
７）次に、吸収成分の濃度演算について説明する。ここでは、被計測媒体が一種類の吸収成分を含む場合を考える。すると、その吸収成分のボクセルｉにおける濃度偏差ΔＣiは、ベア・ランバート則から導かれる次式：
εΔＣi＝μai−μav (C.2)
で与えられる。但し、εは吸収成分の単位濃度当たりの吸収係数（又は消光係数）であり、分光光度計で測定することができる。さらに、前記吸収成分のボクセルｉにおける濃度Ｃiは、ベア・ランバート則から導かれる次式：
εＣi＝μai (C.3)
で与えられる。従って、吸収成分の既知の吸収係数を用いて上記の各ボクセルｉにおける吸収係数の絶対値μaiから各ボクセルｉにおける吸収成分の濃度Ｃiを算出することができ（Ｓ１１１）、吸収成分濃度データファイル５５に記憶される。そして、このようにして求められた各ボクセルｉにおける吸収成分の濃度Ｃiに基づいて被計測媒体内部における吸収成分濃度に関する分布が求められ、その分布を示す光ＣＴ画像が出力装置８０に表示される（Ｓ１１２）。
【０１２９】
（第２実施形態）
上記の第１実施形態の方法において２種類以上の波長の光を用いれば、２種類以上の吸収成分の濃度偏差や濃度絶対値が測定できる。本実施形態では、波長λ1とλ2の２つの光で計測する２波長分光計測法について説明する。すなわち、波長λ1及びλ2の光に対応する各吸収成分の単位濃度当たりの吸収係数（又は消光係数）ε1及びε2が既知であれば、前記(C.2)式及び(C.3)式に対して２本の式が成立する。従って、これらの式を連立させることによって、２種類の吸収成分の濃度偏差や濃度絶対値を測定することができる。
【０１３０】
より具体的には、散乱吸収体が少なくとも２つの吸収成分（例えば酸素化型及び脱酸素化型ヘモグロビン）を含有する場合は、それらの吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有する測定光を使用し、各波長を有する測定光のそれぞれに関して測定値Ｙn、推定値Ｙavn、基準値μav及び重み関数Ｗiを求め、それらに基づいて各波長を有する測定光のそれぞれに関して吸収係数偏差Δμai、吸収係数絶対値μaiを求めることによって、各吸収成分の濃度Ｃiの分布が求められる。
【０１３１】
以下に、上記２波長分光法を利用したヘモグロビンの濃度の計測について説明する。哺乳類の脳における吸収成分の主なものは、水、チトクローム（cytochrom)、酸素化型及び脱酸素化型ヘモグロビンである。近赤外線領域での水及びチトクロームの吸収は、酸素化型及び脱酸素化型ヘモグロビンに対して、ほぼ無視することができる程度に少ない。また、酸素化型及び脱酸素化型ヘモグロビンは、図６に示すように、吸収スペクトルが異なる。さらに、頭蓋骨は、近赤外線に対して、散乱体と考えてよい。
【０１３２】
いま前節までに述べた方法で波長λ1 とλ2 の２種の波長の光に対して吸収係数μa1とμa2が求められたとすれば、ランバート・ベール（Lambert-Beer) 則によって、次式が成立する。
【０１３３】
μa1＝εHb,1〔Ｈｂ〕＋εHbO,1 〔ＨｂＯ〕
μa2＝εHb,2〔Ｈｂ〕＋εHbO,2 〔ＨｂＯ〕
但し、
εHb,1；脱酸素化型ヘモグロビンの波長λ1 に対するモル吸収係数〔ｍｍ-1・Ｍ-1〕
εHbO,1 ；酸素化型ヘモグロビンの波長λ1 に対するモル吸収係数〔ｍｍ-1・Ｍ-1〕
εHb,2；脱酸素化型ヘモグロビンの波長λ2 に対するモル吸収係数〔ｍｍ-1・Ｍ-1〕
εHbO,2 ；酸素化型ヘモグロビンの波長λ2 に対するモル吸収係数〔ｍｍ-1・Ｍ-1〕
〔Ｈｂ〕：脱酸素化型ヘモグロビンのモル濃度〔Ｍ〕
〔ＨｂＯ〕：酸素化型ヘモグロビンのモル濃度〔Ｍ〕
である。
【０１３４】
従って、既知のパラメータεHb,1、εHbO,1 、εHb,2、εHbO,2 及び計測値から演算されたμa1とμa2から、脱酸素化型ヘモグロビンのモル濃度〔Ｈｂ〕、及び酸素化型ヘモグロビンのモル濃度〔ＨｂＯ〕を求めることができる。
【０１３５】
また、上記に対してチトクロームを考慮する場合のように、吸収スペクトルが既知である３成分のそれぞれの濃度の定量は、３波長以上の光を使用すればよい。一般的には、吸収スペクトルが既知であるｎ個の成分の濃度の定量計測は、ｎ個又は（ｎ＋１）個の波長に対する吸収係数の計測値から、上記と同様にして求めることができる。
【０１３６】
さらに、飽和度Ｙは、
Ｙ＝〔ＨｂＯ〕／（〔Ｈｂ〕＋〔ＨｂＯ〕）
であるから、
μa1／μa2＝〔εHb,1＋Ｙ（εHbO,1 −εHb,1）〕÷〔εHb,2＋Ｙ（εHbO,2 −εHb,2）〕
を用いて、既知のパラメータεHb,1、εHbO,1 、εHb,2、εHbO,2 と計測値から演算されたμa1及びμa2とから、飽和度Ｙが容易に算出される。
【０１３７】
以上の方法では、本発明によって各波長の光に対する吸収係数μa1とμa2が高精度に求められるため、各濃度も高精度に求められる。なお、酸素化型及び脱酸素化型ヘモグロビンに対して吸収が同一値になる波長（８００ｎｍ、isosbestic wavelength)を使用すれば上記の式はさらに簡単になる。
【０１３８】
（第３実施形態）
本実施形態は本発明を位相変調計測に応用する例を示す。この場合、前記の第１実施形態における入射光を振幅変調光、所定パラメータを検出光の振幅Ａと位相遅れφ、検出光の光路長分布Ｚi（μav）をｃ（媒体内の光速）倍の群遅延
【数８４】
【０１３９】
なお、この場合は、求めるべき重み関数Ｗiは次式：
【数８５】
【０１４０】
以上、本発明の好適な実施形態について説明したが、本発明は勿論上記実施形態に限定されるものではない。
【０１４１】
すなわち、上記実施形態においては吸収係数の基準値及び輸送散乱係数の基準値を光ＣＴ装置自身で得られたデータから求めているが、これらの基準値を別の装置で求めるようにしてもよい。この場合の利点は、例えば光ＣＴ装置で得られるデータはＣＷ（連続光）で測定し、上記基準値を求める装置でのみパルス光や変調光を用いればよくなるため、光ＣＴ装置のシステムの構成が簡単になる。なお、別の装置で上記基準値を求める手法は位相変調法や時間分解分光法であってもよい。
【０１４２】
また、上記実施形態においては光入射位置を走査させていたが、光入射位置と同期して光検出位置を走査させてもよい。また、被計測媒体の周囲に複数の光入射ファイバー及び複数の光検出ファイバーを配置しておき、光入射に使用するファイバーをそれぞれ順次変えていくことによって光入射位置を移動させてもよい。
【０１４３】
更に、以上説明した本発明の装置及び方法は、一般的な断層像を得る光ＣＴだけでなく、３次元媒体内部の吸収成分を３次元計測する装置及び方法に利用できることは明らかである。また、本発明の装置及び方法は、頭部の皮膚、頭蓋骨、灰白質、白質のように、多数の層状の構造をもつ媒体の計測に利用できることも明白である。なお、この場合は、それぞれの層をそれぞれのボクセルに対応させればよい。
【０１４４】
【発明の効果】
以上説明したように、本発明によれば、新規な重み関数を使用して個々の計測に際して各ボクセルにおける重み関数がマイクロ・ベア・ランバート則に基づいて直接的に求められるため、計測状況に応じた適切な重み関数に基づいて吸収係数の偏差が算出され、従って近似法の採用に伴う誤差の発生が防止される。そのため、本発明によれば、生体などの不均一系の散乱吸収体について適用した場合であっても各ボクセルにおける吸収係数の偏差を正確に求めることが可能となり、このような吸収係数偏差に基づいて高精度の内部特性分布を得て光ＣＴ画像化することが可能となる。
【０１４５】
このように、本発明によって再入射不可能な様々な外形をもつ種々の散乱吸収体の内部吸収成分の分布を高精度かつ迅速に計測することが可能となり、例えば、光マンモグラフィー、頭部光ＣＴ、全身用光ＣＴ、臨床モニター、診断や解析、手術や治療への応用も可能となる。
【図面の簡単な説明】
【図１】不均一媒体内の光子移動に関するモデルを示す模式図である。
【図２】本発明の散乱吸収体の内部特性分布の計測装置の一実施形態を示す模式図である。
【図３】本発明の散乱吸収体の内部特性分布の計測方法の一実施形態を示すフローチャートである。
【図４】本発明にかかる重み関数の演算方法の一実施形態を示すフローチャートである。
【図５】各ボクセルにおける吸収係数分布（吸収係数の基準値と吸収係数の偏差との関係）を示すグラフである。
【図６】ヘモグロビン及びミオグロビンの吸収スペクトルを示すグラフである。
【符号の説明】
１…内部特性分布計測装置、１０…被計測媒体、２１…光入射ファイバー、２２…波長選択器、２３…光源、２５…光検出ファイバー、２６…光検出器、３０…ＣＰＵ、４０…第１記憶装置、４１…ＯＳ、４２…内部特性分布計測プログラム、５０…第２記憶装置、５１…基本光路長分布データファイル、５２…測定データファイル、５３…入力データファイル、５４…吸収係数データファイル、５５…吸収成分濃度データファイル、６１…画像メモリ、６２…作業用メモリ、７０…入力装置、７１…キーボード、７２…マウス、８０…出力装置、８１…ディスプレイ、８２…プリンタ。
Claims (18)
- ヒトを除く散乱吸収体である被計測媒体中に一以上の光入射位置から順次光を入射する光入射ステップと、
前記被計測媒体内部を透過した光を複数の光検出位置で検出する光検出ステップと、
検出された各光に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記光の所定パラメータの測定値を取得する測定値取得ステップと、
前記被計測媒体の吸収係数の基準値を設定する基準値設定ステップと、
前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記パラメータの推定値を求める推定値算出ステップと、
前記吸収係数の基準値を用いて、複数のボクセルに分割された前記被計測媒体の各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて求める重み関数演算ステップと、
前記パラメータの測定値、前記パラメータの推定値及び前記重み関数に基づいて、前記各ボクセルにおける前記吸収係数の基準値に対する吸収係数の偏差を算出する吸収係数偏差算出ステップと、
前記吸収係数の基準値及び前記吸収係数の偏差に基づいて前記各ボクセルにおける吸収係数の絶対値を算出して前記被計測媒体における吸収係数の絶対値分布を求める吸収係数絶対値算出ステップと、を含み、
前記重み関数は、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける平均光路長と光路長分布の分散との関数であることを特徴とする、散乱吸収体の内部特性分布の計測方法。 - 前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける平均光路長を求める平均光路長取得ステップを更に含んでおり、
前記重み関数演算ステップにおいて、下記の式：
に基づいて前記重み関数を求めることを特徴とする、請求項１記載の方法。 - ヒトを除く散乱吸収体である被計測媒体中に一以上の光入射位置から順次光を入射する光入射ステップと、
前記被計測媒体内部を透過した光を複数の光検出位置で検出する光検出ステップと、
検出された各光に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記光の所定パラメータの測定値を取得する測定値取得ステップと、
前記被計測媒体の吸収係数の基準値を設定する基準値設定ステップと、
前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記パラメータの推定値を求める推定値算出ステップと、
前記吸収係数の基準値を用いて、複数のボクセルに分割された前記被計測媒体の各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて求める重み関数演算ステップと、
前記パラメータの測定値、前記パラメータの推定値及び前記重み関数に基づいて、前記各ボクセルにおける前記吸収係数の基準値に対する吸収係数の偏差を算出する吸収係数偏差算出ステップと、
前記吸収係数の基準値及び前記吸収係数の偏差に基づいて前記各ボクセルにおける吸収係数の絶対値を算出して前記被計測媒体における吸収係数の絶対値分布を求める吸収係数絶対値算出ステップと、を含み、
前記重み関数は、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける所定時間領域内の平均光路長と光路長分布の分散との関数であることを特徴とする、散乱吸収体の内部特性分布の計測方法。 - 前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける所定時間領域内の平均光路長を求める平均光路長取得ステップを更に含んでおり、
前記重み関数演算ステップにおいて、下記の式：
に基づいて前記重み関数を求めることを特徴とする、請求項３記載の方法。 - ヒトを除く散乱吸収体である被計測媒体中に一以上の光入射位置から順次光を入射する光入射ステップと、
前記被計測媒体内部を透過した光を複数の光検出位置で検出する光検出ステップと、
検出された各光に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記光の所定パラメータの測定値を取得する測定値取得ステップと、
前記被計測媒体の吸収係数の基準値を設定する基準値設定ステップと、
前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記パラメータの推定値を求める推定値算出ステップと、
前記吸収係数の基準値を用いて、複数のボクセルに分割された前記被計測媒体の各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて求める重み関数演算ステップと、
前記パラメータの測定値、前記パラメータの推定値及び前記重み関数に基づいて、前記各ボクセルにおける前記吸収係数の基準値に対する吸収係数の偏差を算出する吸収係数偏差算出ステップと、
前記吸収係数の基準値及び前記吸収係数の偏差に基づいて前記各ボクセルにおける吸収係数の絶対値を算出して前記被計測媒体における吸収係数の絶対値分布を求める吸収係数絶対値算出ステップと、を含み、
前記重み関数は、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける群遅延と分布の分散との関数であることを特徴とする、散乱吸収体の内部特性分布の計測方法。 - 前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける群遅延を求める群遅延取得ステップを更に含んでおり、
前記重み関数演算ステップにおいて、下記の式：
に基づいて前記重み関数を求めることを特徴とする、請求項５記載の方法。 - 前記吸収係数の絶対値を用いて、前記各ボクセルにおける吸収成分の濃度を算出して前記被計測媒体における吸収成分の濃度分布を求める濃度算出ステップを更に含むことを特徴とする、請求項１〜６のうちのいずれか一項記載の方法。
- 前記被計測媒体が少なくとも２つの吸収成分を含有しており、
前記光入射ステップにおいて前記被計測媒体中に入射される光が、前記吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有しており、
前記光検出ステップにおいて前記少なくとも２つの波長を有する光をそれぞれ検出し、
前記測定値取得ステップにおいて前記少なくとも２つの波長を有する光に関してそれぞれ前記測定値を取得し、
前記基準値設定ステップにおいて前記少なくとも２つの波長を有する光に関してそれぞれ前記基準値を設定し、
前記推定値算出ステップにおいて前記少なくとも２つの波長を有する光に関してそれぞれ前記推定値を求め、
前記重み関数演算ステップにおいて前記少なくとも２つの波長を有する光に関してそれぞれ前記重み関数を求め、
前記吸収係数偏差算出ステップにおいて前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収係数の偏差を算出し、
前記吸収係数絶対値算出ステップにおいて前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収係数の絶対値を算出し、
前記濃度算出ステップにおいて前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収成分の濃度を算出し、前記被計測媒体における前記各吸収成分の濃度分布を求めることを特徴とする、請求項７記載の方法。 - 前記の求められた分布に基づいて、前記被計測媒体における該分布を示す光ＣＴ画像を表示する画像表示ステップを更に含むことを特徴とする、請求項１〜８のうちのいずれか一項記載の方法。
- 散乱吸収体である被計測媒体中に一以上の光入射位置から順次光を入射する光入射手段と、
前記被計測媒体内部を透過した光を複数の光検出位置で検出する光検出手段と、
検出された各光に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記光の所定パラメータの測定値を取得する測定値取得手段と、
前記被計測媒体の吸収係数の基準値を設定する基準値設定手段と、
前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記パラメータの推定値を求める推定値算出手段と、
前記吸収係数の基準値を用いて、複数のボクセルに分割された前記被計測媒体の各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて求める重み関数演算手段と、
前記パラメータの測定値、前記パラメータの推定値及び前記重み関数に基づいて、前記各ボクセルにおける前記吸収係数の基準値に対する吸収係数の偏差を算出する吸収係数偏差算出手段と、
前記吸収係数の基準値及び前記吸収係数の偏差に基づいて前記各ボクセルにおける吸収係数の絶対値を算出して前記被計測媒体における吸収係数の絶対値分布を求める吸収係数絶対値算出手段と、を備え、
前記重み関数は、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける平均光路長と光路長分布の分散との関数であることを特徴とする、散乱吸収体の内部特性分布の計測装置。 - 前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける平均光路長を求める平均光路長取得手段を更に含んでおり、
前記重み関数演算手段において、下記の式：
に基づいて前記重み関数を求めることを特徴とする、請求項１０記載の装置。 - 散乱吸収体である被計測媒体中に一以上の光入射位置から順次光を入射する光入射手段と、
前記被計測媒体内部を透過した光を複数の光検出位置で検出する光検出手段と、
検出された各光に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記光の所定パラメータの測定値を取得する測定値取得手段と、
前記被計測媒体の吸収係数の基準値を設定する基準値設定手段と、
前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記パラメータの推定値を求める推定値算出手段と、
前記吸収係数の基準値を用いて、複数のボクセルに分割された前記被計測媒体の各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて求める重み関数演算手段と、
前記パラメータの測定値、前記パラメータの推定値及び前記重み関数に基づいて、前記各ボクセルにおける前記吸収係数の基準値に対する吸収係数の偏差を算出する吸収係数偏差算出手段と、
前記吸収係数の基準値及び前記吸収係数の偏差に基づいて前記各ボクセルにおける吸収係数の絶対値を算出して前記被計測媒体における吸収係数の絶対値分布を求める吸収係数絶対値算出手段と、を備え、
前記重み関数は、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける所定時間領域内の平均光路長と光路長分布の分散との関数であることを特徴とする、散乱吸収体の内部特性分布の計測装置。 - 前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける所定時間領域内の平均光路長を求める平均光路長取得手段を更に含んでおり、
前記重み関数演算手段において、下記の式：
に基づいて前記重み関数を求めることを特徴とする、請求項１２記載の装置。 - 散乱吸収体である被計測媒体中に一以上の光入射位置から順次光を入射する光入射手段と、
前記被計測媒体内部を透過した光を複数の光検出位置で検出する光検出手段と、
検出された各光に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記光の所定パラメータの測定値を取得する測定値取得手段と、
前記被計測媒体の吸収係数の基準値を設定する基準値設定手段と、
前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記光入射位置と前記光検出位置との組み合わせ毎に前記パラメータの推定値を求める推定値算出手段と、
前記吸収係数の基準値を用いて、複数のボクセルに分割された前記被計測媒体の各ボクセルにおける重み関数をマイクロ・ベア・ランバート則に基づいて求める重み関数演算手段と、
前記パラメータの測定値、前記パラメータの推定値及び前記重み関数に基づいて、前記各ボクセルにおける前記吸収係数の基準値に対する吸収係数の偏差を算出する吸収係数偏差算出手段と、
前記吸収係数の基準値及び前記吸収係数の偏差に基づいて前記各ボクセルにおける吸収係数の絶対値を算出して前記被計測媒体における吸収係数の絶対値分布を求める吸収係数絶対値算出手段と、を備え、
前記重み関数は、前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなしたときの、前記各ボクセルにおける群遅延と分布の分散との関数であることを特徴とする、散乱吸収体の内部特性分布の計測装置。 - 前記被計測媒体が全体的に一様な前記吸収係数の基準値を有するとみなして、前記吸収係数の基準値に基づいて、前記各ボクセルにおける群遅延を求める群遅延取得手段を更に含んでおり、
前記重み関数演算手段において、下記の式：
に基づいて前記重み関数を求めることを特徴とする、請求項１４記載の装置。 - 前記吸収係数の絶対値を用いて、前記各ボクセルにおける吸収成分の濃度を算出して前記被計測媒体における吸収成分の濃度分布を求める濃度算出手段を更に備えることを特徴とする、請求項１０〜１５のうちのいずれか一項記載の装置。
- 前記被計測媒体が少なくとも２つの吸収成分を含有しており、
前記光入射手段において前記被計測媒体中に入射される光が、前記吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有しており、
前記光検出手段において前記少なくとも２つの波長を有する光をそれぞれ検出し、
前記測定値取得手段において前記少なくとも２つの波長を有する光に関してそれぞれ前記測定値を取得し、
前記基準値設定手段において前記少なくとも２つの波長を有する光に関してそれぞれ前記基準値を設定し、
前記推定値算出手段において前記少なくとも２つの波長を有する光に関してそれぞれ前記推定値を求め、
前記重み関数演算手段において前記少なくとも２つの波長を有する光に関してそれぞれ前記重み関数を求め、
前記吸収係数偏差算出手段において前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収係数の偏差を算出し、
前記吸収係数絶対値算出手段において前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収係数の絶対値を算出し、
前記濃度算出手段において前記少なくとも２つの波長を有する光に関してそれぞれ前記吸収成分の濃度を算出し、前記被計測媒体における前記各吸収成分の濃度分布を求めることを特徴とする、請求項１６記載の装置。 - 前記の求められた分布に基づいて、前記被計測媒体における該分布を示す光ＣＴ画像を表示する画像表示手段を更に備えることを特徴とする、請求項１０〜１７のうちのいずれか一項記載の装置。
Priority Applications (6)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP14430098A JP3887486B2 (ja) | 1998-05-26 | 1998-05-26 | 散乱吸収体の内部特性分布の計測方法及び装置 |
| PCT/JP1999/002696 WO1999061893A1 (fr) | 1998-05-26 | 1999-05-24 | Procede et appareil permettant de mesurer la repartition de caracteristique d'un corps d'absorption/diffusion |
| AU38509/99A AU3850999A (en) | 1998-05-26 | 1999-05-24 | Method and device for measuring internal characteristic distribution of scattering/absorbing body |
| EP99921231A EP1081486B1 (en) | 1998-05-26 | 1999-05-24 | Method and device for measuring internal characteristic distribution of scattering/absorbing body |
| DE69928392T DE69928392T2 (de) | 1998-05-26 | 1999-05-24 | Verfahren und vorrichtung zur messung der charakteristischen inneren verteilung von einem streuenden/absorbierenden körper |
| US09/716,264 US6335792B1 (en) | 1998-05-26 | 2000-11-21 | Method and apparatus for measuring internal property distribution in scattering medium |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP14430098A JP3887486B2 (ja) | 1998-05-26 | 1998-05-26 | 散乱吸収体の内部特性分布の計測方法及び装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH11337476A JPH11337476A (ja) | 1999-12-10 |
| JP3887486B2 true JP3887486B2 (ja) | 2007-02-28 |
Family
ID=15358876
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP14430098A Expired - Fee Related JP3887486B2 (ja) | 1998-05-26 | 1998-05-26 | 散乱吸収体の内部特性分布の計測方法及び装置 |
Country Status (6)
| Country | Link |
|---|---|
| US (1) | US6335792B1 (ja) |
| EP (1) | EP1081486B1 (ja) |
| JP (1) | JP3887486B2 (ja) |
| AU (1) | AU3850999A (ja) |
| DE (1) | DE69928392T2 (ja) |
| WO (1) | WO1999061893A1 (ja) |
Families Citing this family (46)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3950243B2 (ja) * | 1998-11-05 | 2007-07-25 | 浜松ホトニクス株式会社 | 散乱吸収体の内部情報の計測方法及び装置 |
| WO2000062462A2 (en) * | 1999-04-09 | 2000-10-19 | Sony Electronics, Inc. | Method for switching signal input based on device capability |
| WO2000065330A1 (fr) * | 1999-04-21 | 2000-11-02 | Hamamatsu Photonics K.K. | Procede d'analyse optique de milieu non homogene |
| US6937884B1 (en) * | 1999-09-14 | 2005-08-30 | The Research Foundation Of State University Of New York | Method and system for imaging the dynamics of scattering medium |
| JP4364995B2 (ja) * | 2000-03-21 | 2009-11-18 | 浜松ホトニクス株式会社 | 散乱吸収体内部の光路分布計算方法 |
| GB2371358A (en) * | 2001-01-22 | 2002-07-24 | Optokem Ltd | Light scattering particle characterisation apparatus and detection means |
| JP3623743B2 (ja) * | 2001-02-26 | 2005-02-23 | 株式会社スペクトラテック | 生体情報測定装置 |
| JP2004537731A (ja) * | 2001-08-02 | 2004-12-16 | ザ リサーチ ファウンデーション オブ ステイト ユニヴァーシティ オブ ニューヨーク | 線形方程式系の解を向上させる方法およびシステム |
| WO2003034032A2 (en) | 2001-10-18 | 2003-04-24 | Cargill Robert L | Method and apparatus for quantifying an 'integrated index' of a material medium |
| US20050106744A1 (en) * | 2001-12-26 | 2005-05-19 | Yoshihiko Mizushima | Optical analysis for heterogeneous medium |
| JP2003202287A (ja) * | 2002-01-08 | 2003-07-18 | Hamamatsu Photonics Kk | 散乱吸収体測定方法及び装置 |
| JP4592245B2 (ja) * | 2002-03-25 | 2010-12-01 | 株式会社日立製作所 | 輻射計算方法及びそれを用いた熱吸収分布図作成方法 |
| US6890155B2 (en) | 2002-04-29 | 2005-05-10 | Thomas Cartwright | Fan blade |
| US6850314B2 (en) * | 2002-08-08 | 2005-02-01 | Board Of Reagents University Of Houston | Method for optical sensing |
| US6992762B2 (en) * | 2002-11-11 | 2006-01-31 | Art Advanced Research Technologies Inc. | Method and apparatus for time resolved optical imaging of biological tissues as part of animals |
| US7720525B2 (en) * | 2003-03-12 | 2010-05-18 | New Art Advanced Research Technologies Inc. | Method and apparatus for combining continuous wave and time domain optical imaging |
| DE102004047417B4 (de) * | 2003-09-29 | 2007-01-04 | Gebauer, Gerd, Dr. | Makromolekül- und Aerosoldiagnostik in gasförmiger und flüssiger Umgebung |
| US7920908B2 (en) | 2003-10-16 | 2011-04-05 | David Hattery | Multispectral imaging for quantitative contrast of functional and structural features of layers inside optically dense media such as tissue |
| US7197404B2 (en) * | 2004-03-01 | 2007-03-27 | Richard Andrew Holland | Computation of radiating particle and wave distributions using a generalized discrete field constructed from representative ray sets |
| GB2429058B (en) * | 2004-03-06 | 2008-12-03 | Michael Trainer | Method and apparatus for determining the size and shape of particles |
| JP4517111B2 (ja) * | 2004-06-07 | 2010-08-04 | 独立行政法人情報通信研究機構 | 脳機能測定装置、脳機能測定方法及び脳機能測定プログラム |
| JP2010512904A (ja) * | 2006-12-19 | 2010-04-30 | コーニンクレッカ フィリップス エレクトロニクス エヌ ヴィ | 不透明媒体のイメージング |
| US20080218727A1 (en) * | 2006-12-22 | 2008-09-11 | Art, Advanced Research Technologies Inc. | Method and apparatus for optical image reconstruction using contour determination |
| US20080177163A1 (en) * | 2007-01-19 | 2008-07-24 | O2 Medtech, Inc. | Volumetric image formation from optical scans of biological tissue with multiple applications including deep brain oxygenation level monitoring |
| US8292809B2 (en) | 2008-03-31 | 2012-10-23 | Nellcor Puritan Bennett Llc | Detecting chemical components from spectroscopic observations |
| US9044601B2 (en) * | 2008-10-20 | 2015-06-02 | Dr. Fred J. Currell | System and methods for accelerating simulation of radiation treatment |
| JP5283525B2 (ja) | 2009-01-30 | 2013-09-04 | 富士フイルム株式会社 | 光断層情報の生成方法、光断層情報生成装置及び光断層情報の生成プログラム |
| US8401608B2 (en) * | 2009-09-30 | 2013-03-19 | Covidien Lp | Method of analyzing photon density waves in a medical monitor |
| WO2011063306A1 (en) | 2009-11-19 | 2011-05-26 | Modulated Imaging Inc. | Method and apparatus for analysis of turbid media via single-element detection using structured illumination |
| US8391943B2 (en) | 2010-03-31 | 2013-03-05 | Covidien Lp | Multi-wavelength photon density wave system using an optical switch |
| WO2012050847A2 (en) | 2010-09-28 | 2012-04-19 | Masimo Corporation | Depth of consciousness monitor including oximeter |
| US9775545B2 (en) | 2010-09-28 | 2017-10-03 | Masimo Corporation | Magnetic electrical connector for patient monitors |
| JP2012237595A (ja) * | 2011-05-10 | 2012-12-06 | Sumitomo Electric Ind Ltd | 光トモグラフィ装置 |
| EP2844983A2 (en) * | 2012-04-30 | 2015-03-11 | Mayo Foundation For Medical Education And Research | Spectrometric systems and methods for improved focus localization of time-and space-varying measurements |
| MX2018016101A (es) | 2012-11-07 | 2022-05-19 | Modulated Imaging Inc | Eficiente formacion modulada de imagenes. |
| JP6154613B2 (ja) * | 2013-01-18 | 2017-06-28 | 浜松ホトニクス株式会社 | 断面画像計測装置及び計測方法 |
| RU2533538C1 (ru) * | 2013-08-19 | 2014-11-20 | Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования "Тверской государственный университет" | Способ раздельного определения вероятностей поглощения и рассеяния фотонов на единицу пути в твердых оптических материалах |
| WO2015040110A1 (en) * | 2013-09-19 | 2015-03-26 | L'oréal Sa | Systems and methods for measuring and categorizing colors and spectra of surfaces |
| US9750413B2 (en) * | 2013-10-03 | 2017-09-05 | National Technology & Engineering Solutions Of Sandia, Llc | Massively parallel diffuse optical tomography |
| USD763938S1 (en) | 2014-04-02 | 2016-08-16 | Cephalogics, LLC | Optical sensor array |
| USD763939S1 (en) | 2014-04-02 | 2016-08-16 | Cephalogics, LLC | Optical sensor array liner with optical sensor array pad |
| US10154815B2 (en) | 2014-10-07 | 2018-12-18 | Masimo Corporation | Modular physiological sensors |
| EP3417276B1 (en) * | 2016-02-16 | 2021-06-30 | Stichting Het Nederlands Kanker Instituut- Antoni van Leeuwenhoek Ziekenhuis | Method, apparatus and computer program for estimating a property of an optically diffuse medium |
| DE102017111957B4 (de) * | 2017-05-31 | 2019-05-16 | Bundesrepublik Deutschland, Vertreten Durch Das Bundesministerium Für Wirtschaft Und Energie, Dieses Vertreten Durch Den Präsidenten Der Physikalisch-Technischen Bundesanstalt | Phantom zum Prüfen eines Messgerätes zur zeitaufgelösten diffusen optischen Spektroskopie, insbesondere eines Gewebe-Oximeters und Verfahren zum Prüfen eines Messgerätes zur zeitaufgelösten diffusen optischen Spektroskopie an Gewebe |
| CN108648246A (zh) * | 2018-05-04 | 2018-10-12 | 中国科学院高能物理研究所 | 基于透射成像的多层结构体图像处理方法及装置 |
| US11654635B2 (en) | 2019-04-18 | 2023-05-23 | The Research Foundation For Suny | Enhanced non-destructive testing in directed energy material processing |
Family Cites Families (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3107914B2 (ja) * | 1992-07-20 | 2000-11-13 | 浜松ホトニクス株式会社 | 散乱吸収体内部の吸収情報計測装置及び方法 |
| JP3107927B2 (ja) * | 1992-10-06 | 2000-11-13 | 浜松ホトニクス株式会社 | 散乱吸収体の光学情報計測装置及び方法 |
| JP3433508B2 (ja) * | 1993-12-01 | 2003-08-04 | 浜松ホトニクス株式会社 | 散乱吸収体計測方法及び散乱吸収体計測装置 |
| JP3310782B2 (ja) * | 1994-07-14 | 2002-08-05 | 株式会社日立製作所 | 吸収物質濃度の空間分布画像化装置 |
| JP2780935B2 (ja) * | 1994-09-22 | 1998-07-30 | 浜松ホトニクス株式会社 | 散乱吸収体の吸収成分の濃度計測方法及び装置 |
| JP3433534B2 (ja) * | 1994-11-07 | 2003-08-04 | 浜松ホトニクス株式会社 | 散乱吸収体内の散乱特性・吸収特性の測定方法及び装置 |
| EP0758082B1 (en) * | 1995-08-08 | 2004-11-24 | Technology Research Association of Medical And Welfare Apparatus | Measurement apparatus for internal information in scattering medium |
| JPH0961359A (ja) * | 1995-08-29 | 1997-03-07 | Hamamatsu Photonics Kk | 濃度測定装置 |
| JP3662376B2 (ja) * | 1996-05-10 | 2005-06-22 | 浜松ホトニクス株式会社 | 内部特性分布の計測方法および装置 |
| JP3844815B2 (ja) * | 1996-08-30 | 2006-11-15 | 浜松ホトニクス株式会社 | 散乱体の吸収情報計測方法及び装置 |
-
1998
- 1998-05-26 JP JP14430098A patent/JP3887486B2/ja not_active Expired - Fee Related
-
1999
- 1999-05-24 DE DE69928392T patent/DE69928392T2/de not_active Expired - Lifetime
- 1999-05-24 AU AU38509/99A patent/AU3850999A/en not_active Abandoned
- 1999-05-24 WO PCT/JP1999/002696 patent/WO1999061893A1/ja active IP Right Grant
- 1999-05-24 EP EP99921231A patent/EP1081486B1/en not_active Expired - Lifetime
-
2000
- 2000-11-21 US US09/716,264 patent/US6335792B1/en not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| AU3850999A (en) | 1999-12-13 |
| EP1081486A1 (en) | 2001-03-07 |
| DE69928392D1 (de) | 2005-12-22 |
| WO1999061893A1 (fr) | 1999-12-02 |
| EP1081486A4 (en) | 2002-07-31 |
| DE69928392T2 (de) | 2006-08-03 |
| JPH11337476A (ja) | 1999-12-10 |
| US6335792B1 (en) | 2002-01-01 |
| EP1081486B1 (en) | 2005-11-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP3887486B2 (ja) | 散乱吸収体の内部特性分布の計測方法及び装置 | |
| JP3107927B2 (ja) | 散乱吸収体の光学情報計測装置及び方法 | |
| Delpy et al. | Quantification in tissue near–infrared spectroscopy | |
| JP3662376B2 (ja) | 内部特性分布の計測方法および装置 | |
| US9538926B2 (en) | Speckle contrast optical tomography | |
| US8416421B2 (en) | Optical coherence computed tomography | |
| US5983121A (en) | Absorption information measuring method and apparatus of scattering medium | |
| JPH06129984A (ja) | 散乱吸収体内部の吸収情報計測装置及び方法 | |
| Hebden et al. | The spatial resolution performance of a time‐resolved optical imaging system using temporal extrapolation | |
| JP3950243B2 (ja) | 散乱吸収体の内部情報の計測方法及び装置 | |
| US6975401B2 (en) | Method of calculating optical path distribution inside scattering absorber | |
| Maheswari et al. | Soft tissue optical property extraction for carcinoma cell detection in diffuse optical tomography system under boundary element condition | |
| Giacometti et al. | Diffuse optical tomography for brain imaging: continuous wave instrumentation and linear analysis methods | |
| Sato et al. | Extraction of depth-dependent signals from time-resolved reflectance in layered turbid media | |
| Okada et al. | Experimental measurements on phantoms and Monte Carlo simulation to evaluate the effect of inhomogeneity on optical pathlength | |
| Caffini et al. | Functional near infrared spectroscopy and diffuse optical tomography in neuroscience | |
| Lu et al. | Measurement of food optical properties | |
| Veesa et al. | Functional near infrared spectroscopy using spatially resolved data to account for tissue scattering | |
| Chance et al. | Photon dynamics in tissue imaging | |
| Yakimov et al. | Evaluating the Speckle-SFDI for the Quantification of Optical Properties of Biotissues: Modeling and Validation on Optical Phantoms | |
| Hennessy | Depth resolved diffuse reflectance spectroscopy | |
| Brooksby et al. | Internal refractive index changes affect light transport in tissue | |
| Willmann et al. | Small-volume frequency-domain oximetry: phantom experiments and first in vivo results | |
| Hunter et al. | Measurement of in-vivo hemoglobin concentration using diffuse reflectance | |
| Long et al. | Optical tomography in scattering media from photon time‐of‐flight diffuse reflectance measurements: a chemometric approach |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20041209 |
|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20041209 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20060725 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20060922 |
|
| 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: 20061121 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20061127 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| S531 | Written request for registration of change of domicile |
Free format text: JAPANESE INTERMEDIATE CODE: R313532 |
|
| R350 | Written notification of registration of transfer |
Free format text: JAPANESE INTERMEDIATE CODE: R350 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20091201 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20101201 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111201 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121201 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20131201 Year of fee payment: 7 |
|
| LAPS | Cancellation because of no payment of annual fees |