JP3662376B2 - 内部特性分布の計測方法および装置 - Google Patents
内部特性分布の計測方法および装置 Download PDFInfo
- Publication number
- JP3662376B2 JP3662376B2 JP33467496A JP33467496A JP3662376B2 JP 3662376 B2 JP3662376 B2 JP 3662376B2 JP 33467496 A JP33467496 A JP 33467496A JP 33467496 A JP33467496 A JP 33467496A JP 3662376 B2 JP3662376 B2 JP 3662376B2
- Authority
- JP
- Japan
- Prior art keywords
- average
- coefficient
- light
- absorption
- absorption coefficient
- 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
- 238000009826 distribution Methods 0.000 title claims description 124
- 238000000034 method Methods 0.000 title claims description 117
- 238000010521 absorption reaction Methods 0.000 claims description 262
- 238000005259 measurement Methods 0.000 claims description 203
- 238000001514 detection method Methods 0.000 claims description 121
- 230000008859 change Effects 0.000 claims description 91
- 238000004364 calculation method Methods 0.000 claims description 28
- 230000002745 absorbent Effects 0.000 claims description 3
- 239000002250 absorbent Substances 0.000 claims description 3
- 230000006870 function Effects 0.000 description 69
- 239000000835 fiber Substances 0.000 description 52
- 230000003287 optical effect Effects 0.000 description 33
- 230000014509 gene expression Effects 0.000 description 22
- 238000012545 processing Methods 0.000 description 18
- 108010054147 Hemoglobins Proteins 0.000 description 15
- 102000001554 Hemoglobins Human genes 0.000 description 15
- 238000003384 imaging method Methods 0.000 description 14
- 238000010586 diagram Methods 0.000 description 13
- 239000013307 optical fiber Substances 0.000 description 11
- 238000002591 computed tomography Methods 0.000 description 10
- 238000004088 simulation Methods 0.000 description 9
- 238000007796 conventional method Methods 0.000 description 8
- 238000009792 diffusion process Methods 0.000 description 8
- 210000001519 tissue Anatomy 0.000 description 8
- 239000006096 absorbing agent Substances 0.000 description 6
- 238000004611 spectroscopical analysis Methods 0.000 description 5
- WQZGKKKJIJFFOK-GASJEMHNSA-N Glucose Natural products OC[C@H]1OC(O)[C@H](O)[C@@H](O)[C@@H]1O WQZGKKKJIJFFOK-GASJEMHNSA-N 0.000 description 4
- 238000000862 absorption spectrum Methods 0.000 description 4
- 230000002238 attenuated effect Effects 0.000 description 4
- 239000008103 glucose Substances 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 102000018832 Cytochromes Human genes 0.000 description 3
- 108010052832 Cytochromes Proteins 0.000 description 3
- 108010064719 Oxyhemoglobins Proteins 0.000 description 3
- 239000008280 blood Substances 0.000 description 3
- 210000004369 blood Anatomy 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000012937 correction Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 210000003722 extracellular fluid Anatomy 0.000 description 3
- 239000000463 material Substances 0.000 description 3
- 239000011159 matrix material Substances 0.000 description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 3
- 210000000481 breast Anatomy 0.000 description 2
- 210000004027 cell Anatomy 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 238000009543 diffuse optical tomography Methods 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 230000000644 propagated effect Effects 0.000 description 2
- 230000035945 sensitivity Effects 0.000 description 2
- 238000001179 sorption measurement Methods 0.000 description 2
- 230000007480 spreading Effects 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 238000003325 tomography Methods 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 102000036675 Myoglobin Human genes 0.000 description 1
- 108010062374 Myoglobin Proteins 0.000 description 1
- 238000002835 absorbance Methods 0.000 description 1
- 239000011358 absorbing material Substances 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 239000012472 biological sample Substances 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 210000000038 chest Anatomy 0.000 description 1
- 238000002939 conjugate gradient method Methods 0.000 description 1
- 239000000470 constituent Substances 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 230000031700 light absorption Effects 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 210000005171 mammalian brain Anatomy 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000005012 migration Effects 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 238000012634 optical imaging Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 210000003625 skull Anatomy 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 210000002784 stomach Anatomy 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
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/49—Scattering, i.e. diffuse reflection within a body or fluid
Landscapes
- Physics & Mathematics (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)
- Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
- Length Measuring Devices By Optical Means (AREA)
- Photometry And Measurement Of Optical Pulse Characteristics (AREA)
Description
【０００１】
【発明の属する技術分野】
本発明は、測定対象物の内部特性分布を計測するための方法およびそのための装置に関する。より詳細には、本発明は、光入射位置および光検出位置を測定対象物の表面に沿って移動させて断層像を得る光ＣＴ（コンピュータトモグラフィー）装置等に適用可能な内部特性分布計測方法並びにそのための装置に関する。
【０００２】
【従来の技術】
散乱吸収体である対象物の表面における１つの光入射位置から測定光を入射し、対象物によって散乱されつつ透過した測定光を前記対象物の表面における複数の光検出位置で受光するとともに、前記対象物の表面に添って光入射位置および光検出位置を移動させていき、散乱吸収体における内部特性の分布を求める光ＣＴ装置においては、例えばその内部の吸収係数の分布を求める方法として以下の方法が知られている。すなわち、“Imaging of Multiple Targets in Dense Scattering Media”（H.L.Graber, J.Chang, R.L.Barbour, SPIE vol.2570, p.219-p.234）、“Imaging diffusive media using time-independent and time-harmonic sources；dependence of image quality on imaging algorithms, target volume weight matrix, and view angles ”（Jenghwa Chang et.al., SPIE vol.2389）等に記載された方法がある。
【０００３】
このような従来の方法における基本的な画像化原理は、測定対象物内部を便宜的に複数のボクセルに区切り、対象物表面のある点から入射した光が測定対象物内部を透過して同表面の他の点で受光された際に、各ボクセル毎の吸収係数といった特定の内部特性に注目した場合の受光される光に対する寄与度を示す関数（ここでは便宜的に「広がり関数」という）と、受光された光との関係式を用いることにある。なお、ここでいうボクセルとは、測定対象物を複数の領域に分割した各領域（ボックスセル）のことをいう。
【０００４】
しかしながら、上記従来の方法においては、測定対象物とは別に吸収の無いファントムを用意し、それを用いて基準となる検出光量を測定し、さらにその状態の広がり関数を用いて目的とする散乱吸収体内部の吸収係数分布を求めていた。そして、このような方法を用いて画像化を行う場合、測定対象物と同一又は類似した形状を有しかつ既知の内部特性を持つように作られたファントムモデル（物理モデル）あるいはシミュレーションモデルを想定し、このようなモデルから得たデータを画像化の計算の際の基準値として用いる必要があった。そのため、これら従来の方法では、実際の測定対象物と物理モデルまたはシミュレーションモデルとの差や測定対象物の個体差等に起因する誤差を避けることができず、特に生体のように複雑な構造を有する測定対象への適用は非常に困難であった。
【０００５】
一方、ファントムを用いないで吸収物質濃度の空間分布を求める方法として、特開平８−２９３２９号公報に記載の方法がある。しかしながら、同公報記載の方法においては、測定対象物中の吸収成分が一種類の場合でも複数の波長を有する光を用いる必要があり、これら複数の波長間での平均光路長分布および減衰光量（散乱等の影響で減衰する光量）は同じであると仮定して吸収物質濃度の空間分布を求めていた。また、この方法においては、吸収物質を含まない仮想被検体を想定し、その仮想被検体内における平均光路長を用いて吸収物質濃度の空間分布を求めており、吸収による光路長の変化は考慮されていなかった。そのため、上記公報に記載の方法によっても、得られる内部特性分布の信頼性は未だ充分なものではなかった。
【０００６】
【発明が解決しようとする課題】
そこで本発明は、上記従来の課題に鑑みてなされたものであり、従来は必要であった物理モデルあるいはシミュレーションモデルから基準値を求めることなく、また測定対象物中の一成分に対して複数の波長を有する光を用いることなく、測定対象物についての測定値から直接的に基準値を求めることを可能とし、その基準値に基づいて測定対象物における内部特性分布の信頼性の高いすなわち高精度の測定が可能な方法および装置を提供することを目的とする。
【０００７】
【課題を解決するための手段】
本発明者らは、上記目的を達成すべく鋭意研究した結果、測定対象物の表面にありかつその対象物内の一点（例えば対象物の中心）に対して相対的に同じ位置関係にある光入射位置と光検出位置との複数の組み合わせによって求められた複数の測定値の平均値を内部特性分布を求めるための基準値として用いることによって上記課題が解決されることを見出し、本発明に到達した。
【０００８】
すなわち、本発明の内部特性分布の計測方法は、
測定対象物の表面における複数の光入射位置から順次該対象物中に測定光を入射するステップと、
該対象物中を透過した測定光を、前記対象物の表面における複数の光検出位置のうちの少なくとも１つの光検出位置であってかつ測定されるべき測定光が入射された光入射位置に対して所定の位置関係にある光検出位置で順次あるいは同時に検出するステップと、
各光検出位置で検出された各測定光に基づいて、該測定光の所定パラメータの測定値を求めるステップと、
前記位置関係が相対的に同じである前記光入射位置と前記光検出位置との複数の組み合わせによって求められた複数の前記測定値を抽出し、該測定値の平均値を算出して該位置関係における基準値を得るステップと、
前記複数の組み合わせによって求められた前記複数の測定値と、前記基準値とを用いて、複数の領域に分割された前記対象物の各領域における所定内部特性の変化量を算出して該対象物における内部特性変化量分布を求めるステップと、
を具備する方法である。
【０００９】
また、本発明の内部特性分布の計測装置は、
測定対象物の表面における複数の光入射位置から順次該対象物中に測定光を入射する光入射手段と、
該対象物中を透過した測定光を、前記対象物の表面における複数の光検出位置のうちの少なくとも１つの光検出位置であってかつ測定されるべき測定光が入射された光入射位置に対して所定の位置関係にある光検出位置で順次あるいは同時に検出する光検出手段と、
各光検出位置で検出された各測定光に基づいて、該測定光の所定パラメータの測定値を求める測定値取得手段と、
前記位置関係が相対的に同じである前記光入射位置と前記光検出位置との複数の組み合わせによって求められた複数の前記測定値を抽出し、該測定値の平均値を算出して該位置関係における基準値を得る基準値算出手段と、
前記複数の組み合わせによって求められた前記複数の測定値と、前記基準値とを用いて、複数の領域に分割された前記対象物の各領域における所定内部特性の変化量を算出して該対象物における内部特性変化量分布を求める内部特性変化量算出手段と、
を具備する装置である。
【００１０】
本発明の方法および装置では、測定対象物の表面にありかつその対象物内の一点（例えば対象物の中心）に対して相対的に同じ位置関係にある光入射位置と光検出位置との複数の組み合わせによって求められた複数の測定値の平均値を、内部特性分布を求めるための基準値として用いる。すなわち、上記基準値と各測定値とを用いて後述する方程式を解くことよって、複数の領域に分割された前記対象物の各領域における内部特性の変化量（差分）が求められる。
【００１１】
このように本発明では、実際の測定で得られる測定値の平均値から基準値を求め、この基準値に基づいて内部特性の変化量が算出される。したがって、本発明では物理モデルやシミュレーションモデルから予め求めた基準値を用いていないため、測定対象物の個体差や、実際の測定対象物と物理モデルやシミュレーションモデルとの間に生じる条件の差等に起因する誤差が発生する余地がなく、信頼性の高い高精度の測定が可能である。さらに、本発明によれば、物理モデル等を用いて基準値を予め求める手間が省かれ、測定時間が短縮化される。
【００１２】
基準値を実際の測定で得られる測定値の平均値から求め、その基準値に基づいて内部特性の変化量等を求めるということは、例えば図１で説明すると、Ａ、Ｂ、Ｃの各値は０からの値が分からなくても平均値との差分が分かれば求められるということと同様である。本発明の動作原理はこのような原理に基づいて内部特性の差分や絶対値を求めることである。
【００１３】
また、本発明では測定対象物中の一成分に対して複数の波長を有する光を用いることなく測定対象物中の内部特性分布が求められるため、複数の波長間での平均光路長分布および減衰光量（散乱等の影響で減衰する光量）は同じであるといった仮定に基づく誤差の発生が防止され、測定精度が向上する。さらに、本発明によれば、測定対象物中の多成分を複数の波長を有する光を用いて解析する場合においても誤差の発生は同様に防止される。すなわち、本発明においては、現実の物体には存在する散乱係数の波長依存性の変化にも対応して各波長毎に散乱係数を求めているためである。
【００１４】
なお、本発明にかかる光入射位置と光検出位置との位置関係は、例えば測定対象物の中心を基準に、すなわち光入射位置および該中心を結ぶ線と光検出位置および該中心を結ぶ線との間の角度で規定され、かかる位置関係が相対的に同じであるということは例えば上記角度が同じであることをいう。
【００１５】
また、本発明にかかる測定値としては、測定光の測定対象物内部での散乱および吸収に関係する所定パラメータの測定値が好ましく、例えば測定光の光量、位相差（または位相遅れ）、振幅、時間分解波形等のパラメータの測定値が好適に用いられる。
【００１６】
本発明の方法および装置によって測定可能な内部特性としては、吸収係数、等価散乱係数、屈折率が挙げられ、これらのうちのいずれかの特性を単独で、あるいは複数の特性を同時または順次に求めてもよい。
【００１７】
先ず、本発明の方法および装置によって測定されるべき内部特性が吸収係数である場合について説明する。
【００１８】
この場合、本発明の方法は、前記対象物の平均の吸収係数および平均の等価散乱係数を求める（好ましくは前記基準値に基づいて求める）ステップと、前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数（吸収係数用の広がり関数）を選択するステップとを更に具備していることが好ましく、前記内部特性変化量分布を求めるステップにおいて前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における吸収係数の変化量を算出することができる。
【００１９】
また、本発明の装置は、前記対象物の平均の吸収係数および平均の等価散乱係数を求める（好ましくは前記基準値に基づいて求める）平均吸収および散乱係数検出手段と、前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数（吸収係数用の広がり関数）を選択する広がり関数選択手段とを更に具備していることが好ましく、前記内部特性変化量算出手段において前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における吸収係数の変化量を算出することができる。
【００２０】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての吸収係数および等価散乱係数に対応して選択された広がり関数に基づいて、前記の各領域における吸収係数の変化量（差分）が求められる。したがって、このような本発明の方法および装置によれば、吸収係数および／または等価散乱係数をゼロとして計算する場合に比べて、それらによる実効光路長の変化に基づく誤差の発生が充分に防止され、測定精度が向上する。
【００２１】
さらに、上記本発明の方法は、前記吸収係数の変化量および前記平均の吸収係数を用いて、前記各領域における吸収係数の絶対値を算出して前記対象物における吸収係数絶対値分布を求めるステップ、および／または前記吸収係数の絶対値を用いて、前記各領域における吸収成分の濃度を算出して前記対象物における吸収成分濃度分布を求めるステップを更に具備してもよい。
【００２２】
また、上記本発明の装置は、前記吸収係数の変化量および前記平均の吸収係数を用いて、前記各領域における吸収係数の絶対値を算出して前記対象物における吸収係数絶対値分布を求める吸収係数絶対値算出手段、および／または前記吸収係数の絶対値を用いて、前記各領域における吸収成分の濃度を算出して前記対象物における吸収成分濃度分布を求める吸収成分濃度算出手段を更に具備してもよい。
【００２３】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての吸収係数に基づいて、前記の各領域における吸収係数の変化量（差分）から各領域における吸収係数の絶対値が求められる。このように本発明の方法および装置によれば、均一な吸収係数を有しかつ測定対象物と同一の外形を持つファントムから得た基準値を用いることなく各領域における吸収係数の絶対値が求められることとなる。そして、各領域における吸収係数の絶対値が求まれば、既知の吸収成分のモル吸収係数等を用いて各領域における吸収成分の濃度が求められる。本発明で求められる吸収係数の変化量分布の誤差は従来の方法による場合と比較して極めて少なくなるため、それに基づいて求められる吸収係数の絶対値分布並びに吸収成分の濃度分布の精度が高くなる。
【００２４】
上記本発明の方法を少なくとも２つの吸収成分を含有している対象物に適用する場合には、前記光入射ステップにおいて前記対象物中に入射される測定光が、該吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有していることが好ましい。この場合、前記光検出ステップにおいて前記少なくとも２つの波長を有する測定光をそれぞれ検出し、前記測定値を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記測定値を求め、前記基準値を得るステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記平均値を算出し、前記内部特性変化量分布を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の変化量を算出し、前記吸収係数絶対値分布を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の絶対値を算出し、前記吸収成分濃度分布を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収成分の濃度を算出することによって、前記対象物における前記各吸収成分の濃度分布を高精度に求めることが可能となる。
【００２５】
また、上記本発明の装置を用いて少なくとも２つの吸収成分を含有している対象物を測定する場合には、前記光入射手段において前記対象物中に入射される測定光が、該吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有していることが好ましい。この場合、前記光検出手段において前記少なくとも２つの波長を有する測定光をそれぞれ検出し、前記測定値取得手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記測定値を求め、前記基準値算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記平均値を算出し、前記内部特性変化量算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の変化量を算出し、前記吸収係数絶対値算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の絶対値を算出し、前記吸収成分濃度算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収成分の濃度を算出することによって、前記対象物における前記各吸収成分の濃度分布を高精度に求めることが可能となる。
【００２６】
次に、本発明の方法および装置によって測定されるべき内部特性が等価散乱係数である場合について説明する。
【００２７】
この場合、本発明の方法は、前記対象物の平均の吸収係数および平均の等価散乱係数を求める（好ましくは前記基準値に基づいて求める）ステップと、前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数（等価散乱係数用の広がり関数）を選択するステップとを更に具備していることが好ましく、前記内部特性変化量分布を求めるステップにおいて前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における等価散乱係数の変化量を算出することができる。
【００２８】
また、本発明の装置は、前記対象物の平均の吸収係数および平均の等価散乱係数を求める（好ましくは前記基準値に基づいて求める）平均吸収および散乱係数検出手段と、前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数（等価散乱係数用の広がり関数）を選択する広がり関数選択手段とを更に具備していることが好ましく、前記内部特性変化量算出手段において前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における等価散乱係数の変化量を算出することができる。
【００２９】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての吸収係数および等価散乱係数に対応して選択された広がり関数に基づいて、前記の各領域における等価散乱係数の変化量（差分）が求められる。したがって、このような本発明の方法および装置によれば、吸収係数および／または等価散乱係数をゼロとして計算する場合に比べて、それらによる実効光路長の変化に基づく誤差の発生が充分に防止され、測定精度が向上する。
【００３０】
さらに、上記本発明の方法は、前記等価散乱係数の変化量および前記平均の等価散乱係数を用いて、前記各領域における等価散乱係数の絶対値を算出して前記対象物における等価散乱係数絶対値分布を求めるステップを更に具備してもよい。
【００３１】
また、上記本発明の装置は、前記等価散乱係数の変化量および前記平均の等価散乱係数を用いて、前記各領域における等価散乱係数の絶対値を算出して前記対象物における等価散乱係数絶対値分布を求める等価散乱係数絶対値算出手段を更に具備してもよい。
【００３２】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての等価散乱係数に基づいて、前記の各領域における等価散乱係数の変化量（差分）から各領域における等価散乱係数の絶対値が求められる。このように本発明の方法および装置によれば、均一な等価散乱係数を有しかつ測定対象物と同一の外形を持つファントムから得た基準値を用いることなく各領域における等価散乱係数の絶対値が求められることとなる。本発明で求められる等価散乱係数の変化量分布の誤差は従来の方法による場合と比較して極めて少なくなるため、それに基づいて求められる等価散乱係数の絶対値分布の精度が高くなる。
【００３３】
次に、本発明の方法および装置によって測定されるべき内部特性が屈折率である場合について説明する。
【００３４】
この場合、本発明の方法は、前記対象物の平均の吸収係数、平均の等価散乱係数および平均の屈折率を求める（好ましくは前記基準値に基づいて求める）ステップと、前記平均の吸収係数、平均の等価散乱係数および平均の屈折率に対応する広がり関数（屈折率用の広がり関数）を選択するステップとを更に具備していることが好ましく、前記内部特性変化量分布を求めるステップにおいて前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における屈折率の変化量を算出することができる。
【００３５】
また、本発明の装置は、前記対象物の平均の吸収係数、平均の等価散乱係数および平均の屈折率を求める（好ましくは前記基準値に基づいて求める）平均吸収および散乱係数検出手段と、前記平均の吸収係数、平均の等価散乱係数および平均の屈折率に対応する広がり関数（屈折率用の広がり関数）を選択する広がり関数選択手段とを更に具備していることが好ましく、前記内部特性変化量算出手段において前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における屈折率の変化量を算出することができる。
【００３６】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての吸収係数、等価散乱係数および屈折率に対応して選択された広がり関数に基づいて、前記の各領域における屈折率の変化量（差分）が求められる。したがって、このような本発明の方法および装置によれば、吸収係数および／または等価散乱係数をゼロとして計算する場合に比べて、それらによる実効光路長の変化に基づく誤差の発生が充分に防止され、測定精度が向上する。
【００３７】
さらに、上記本発明の方法は、前記屈折率の変化量および前記平均の屈折率を用いて、前記各領域における屈折率の絶対値を算出して前記対象物における屈折率絶対値分布を求めるステップを更に具備してもよい。
【００３８】
また、上記本発明の装置は、前記屈折率の変化量および前記平均の屈折率を用いて、前記各領域における屈折率の絶対値を算出して前記対象物における屈折率絶対値分布を求める屈折率絶対値算出手段を更に具備してもよい。
【００３９】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての屈折率に基づいて、前記の各領域における屈折率の変化量（差分）から各領域における屈折率の絶対値が求められる。このように本発明の方法および装置によれば、均一な屈折率を有しかつ測定対象物と同一の外形を持つファントムから得た基準値を用いることなく各領域における屈折率の絶対値が求められることとなる。本発明で求められる屈折率の変化量分布の誤差は従来の方法による場合と比較して極めて少なくなるため、それに基づいて求められる屈折率の絶対値分布の精度が高くなる。
【００４０】
上述の本発明の方法は、前記の求められた分布に基づいて、前記対象物内部における分布を示す画像を表示するステップを更に具備してもよい。また、上述の本発明の装置は、前記の求められた分布に基づいて、前記対象物内部における分布を示す画像を表示する画像表示手段を更に具備してもよい。このような本発明の方法および装置によれば、高精度に求められた内部特性分布を画像化して表示することが可能である。
【００４１】
【発明の実施の形態】
以下、図面を参照しつつ本発明の好適な実施形態について詳細に説明する。尚、図面中、同一又は相当部分には同一符号を付することとする。
【００４２】
先ず、図２および図３を参照して本実施形態で使用した光ＣＴの画像化原理を説明する。なお、散乱されつつ進行する光に関しては３次元座標を用いて考える必要があるが、以下では説明を簡単にするために２次元座標を用いて説明する。
【００４３】
先ず、散乱吸収体内部をＮ個のボクセルに区切り、吸収係数が存在する条件下における散乱吸収体に対する入射光量と出射光量（検出光量）との関係を考える。等価散乱係数μ'sおよび吸収係数μaを均一とした散乱吸収体内部の模式図を図２に示す（Ｎ＝２５）。このとき入射光量をＩ0、検出光量をＩd0、散乱吸収体内部の等価散乱係数μ'sおよび吸収係数μaを均一とした場合の各ボクセル内での実効光路長をＷj、散乱および反射等によって入射光が散乱吸収体の外に出る割合を示す減衰定数をＤsrとすると、下記式（１）が成り立つ。
【００４４】
【数１】
【００４５】
次に、幾つかのボクセル内に等価散乱係数は同じであるが吸収係数が相違する媒質を挿入した以外は図２に示すものと同様の散乱吸収体内部の模式図を図３に示す。但し、図３に示す散乱吸収体に用いた各媒質の吸収係数μai（i=1,2,・・・,N）と図２に示す散乱吸収体に用いた媒質の吸収係数μaとの関係は以下の式（２）に示すような関係にある。
【００４６】
【数２】
【００４７】
このときの入射光量をＩ0、検出光量をＩd1とし、散乱および反射等によって入射光が散乱吸収体の外に出る割合を示す減衰定数Ｄsrは吸収係数を均一とした場合（図２）と変わらないと考えると、検出光量をＩd1は下記式（３）のように表わされる。
【００４８】
【数３】
【００４９】
したがって、式（３）より以下の式（４）が求められる。
【００５０】
【数４】
【００５１】
このように、基準の光量Ｉd0を用いれば、求めたい吸収係数μaと実際の実験系で測定可能な検出光量Ｉd1との関係より、実効光路長Ｗjが決まれば散乱吸収体内部の吸収係数μaの分布が求められることとなる。式（４）は１組の光入射位置−光検出位置で成り立つ関係を示す。したがって、例えばＮ個の吸収係数（未知数）を求める場合は、光入射位置−光検出位置の組み合わせをＮ組み揃え、それぞれの組み合わせについて成り立つＮ個の式（４）の連立方程式を解くことによってＮ個の吸収係数を求めることができる。
【００５２】
すなわち、光入射位置−光検出位置のＮ個の組み合わせについて成り立つＮ個の式（４）の連立方程式を行列式の形で表わすと、下記式（５）が得られる。
【００５３】
【数５】
【００５４】
但し、ΔＩはｌｎＩd0−ｌｎＩd1、Ｗは各ボクセルの実効光路長の分布を示す広がり関数を表わす。光入射位置Ｍ（Ｍ1〜ＭX）の数をＸ、光検出位置ｍ（ｍ1〜ｍx）の数をｘ、光入射位置Ｍ及び光検出位置ｍのときの光量の変化量をΔＩMm、光入射位置Ｍ及び光検出位置ｍのときの各ボクセルの広がり関数をＷMmとすると、［ΔＩMm］は（Ｘ×ｘ）×１、［ＬMm］は（Ｘ×ｘ）×Ｎ、［Δμan］はＮ×１のマトリックスとなる。したがって、吸収係数の変化量Δμanは下記式（６）の連立方程式を解くことにより求められることとなる。なお、式（６）の連立方程式は、Ｘ×ｘ＝ＮとなるようにＸおよびｘの値を選択して解くことが好ましい。
【００５５】
【数６】
【００５６】
かかる画像再構成方法で散乱吸収体内部の吸収係数の定量を行う場合、基本的に図２に示すような基準となる状態が必要となり、上記の形態では吸収係数が均一の状態を基準として想定していたため式（２）と式（４）から各ボクセルの吸収係数を求めていた。しかしながら、かかる画像化方法においては内部の吸収係数の値とその時の各光検出位置での光量さえ予め分かっていればよく、実際に用いるべき基準の状態は特に制限されない。すなわち、例えばある条件下での内部の吸収係数の値とその時の各光検出位置での光量を基準として内部の吸収係数を求めると、その吸収係数の値は基準の値との差という形で求められる。
【００５７】
従来はこのように基準となる内部の吸収係数の値とその時の各光検出位置での光量を、測定対象である散乱吸収体に対して内部の吸収係数のみが相違するファントムモデルあるいはシミュレーションモデルから得ていた。
【００５８】
しかしながら、本発明においては、測定対象物の表面にありかつその対象物内の一点（例えば対象物の中心）に対して相対的に同じ位置関係にある光入射位置と光検出位置との複数の組み合わせによって求められた複数の測定値の平均値を、内部特性分布を求めるための基準値として用いる。また、各ボクセルの実効光路長の作成方法は、本発明者らによる特願平８−６６１９「光ＣＴ装置及び光ＣＴによる画像再構成方法」等に記載されている。本実施形態においてはこの作成方法にしたがって、ある光入射位置−光検出位置の関係における各ボクセルの実効光路長の分布（すなわち広がり関数）を、吸収係数の平均値、等価散乱係数の平均値等に基づいて予め用意する。このように、本実施形態においては、従来は必要であった物理モデルやシミュレーションモデルから基準値を求めることなく、測定対象物についての測定値から直接的に基準値を求めることが可能となり、その基準値に基づいて測定対象物の内部特性の分布を高精度に測定することが可能となる。
【００５９】
以上、ＣＷ計測を用いる実施形態について説明したが、本発明において時間分解計測を適用することも可能である。すなわち、前記の式（３）においては検出光を時間０〜ｔ（ｓ）間に検出器で受光した検出光量の積分値として表わしているが、かかる関係式は光源にパルス光を用いた場合に検出器で得られる時間分解波形についても成り立つ。式（３）をある時間間隔ｔ1〜ｔ2について書き直すと以下の式（７）が得られる。
【００６０】
【数７】
【００６１】
但し、［Ｉd0］t1-t2、［Ｉd1］t1-t2は各検出光の時間ｔ1〜ｔ2間の時間分解波形の光量、［Ｗj］t1-t2は時間ｔ1〜ｔ2間の広がり関数を示す（ｊは各ボクセルの番号を示す）。また、０≦ｔ1＜ｔ2である。したがって、式（７）より以下の式（８）が求められる。
【００６２】
【数８】
【００６３】
このように、時間分解計測を用いる方法においては、測定時間を様々に区切って方程式の数を増やすことによって、ボクセル数と同じＮ個の方程式を解くことによって各ボクセルの吸収係数が求められる。
【００６４】
また、上記の実施形態においては測定されるべき内部特性が吸収係数の場合について説明したが、本発明は等価散乱係数、屈折率といった他の内部特性に適用することも可能である。すなわち、対象物の内部を透過して受光された光は、対象物の内部がもつ吸収係数、等価散乱係数のみならず、かかる対象物がもつ全ての内部特性の影響を受けており、さらにそれらはそれぞれ線形的かつ独立的に受光された光に作用を及ぼしている。なお、互いに影響を及ぼし合う内部特性を求める場合は、それらを１つの内部特性とみなすことにより独立性が保たれる。
【００６５】
これらの関係により、対象物がもつ全ての内部特性の値は、各ボクセル内における内部特性の受光される光に対する寄与度を示す関数（広がり関数）と、受光された光とを用いた式で表わされ、これらを用いれば例えば“Forward and Inverse Calculations for 3-D Frequency-Domain Diffuse Optical Tomography”（Brian W.Pogue et al., SPIE vol.2389, p.328-p.338）に記載されているように等価散乱係数、吸収係数の分布を求めることができる。このような場合においても、検出光の振幅や位相といったパラメータの基準値の導出法として本発明の方法は適用できる。したがって、吸収係数、等価散乱係数、屈折率といった内部特性に対しても、受光される光に対する寄与を示す関数（広がり関数）と受光された光とを用いた関係式に本発明にかかる基準値導出法を適用することによってこれらの内部特性の定量化が可能である。このような関係式としては、以下の諸式が挙げられる。すなわち、吸収係数および等価散乱係数を求める場合に適用可能な式として、前記式（４）および式（８）の変形式としてそれぞれ下記の式（４’）および式（８’）が挙げられる。
【００６６】
【数９】
【００６７】
【数１０】
【００６８】
また、吸収係数、等価散乱係数および屈折率を求める場合に適用可能な式として、前記式（４）および式（８）の変形式としてそれぞれ下記の式（４''）および式（８''）が挙げられる。
【００６９】
【数１１】
【００７０】
【数１２】
また、１つの内部特性を求める場合であっても、検出光に影響を及ぼす全ての内部特性に対応する広がり関数を用いる方が測定精度が向上する傾向にある。したがって、測定対象が生体のように内部特性として吸収係数、等価散乱係数、屈折率等をもつ場合は、吸収係数の分布についてのみ画像化する場合であっても少なくとも吸収係数の平均値および等価散乱係数の平均値（時間分解計測を行う場合は吸収係数の平均値、等価散乱係数の平均値および屈折率の平均値）に対応する広がり関数を用いることが好ましい場合がある。
【００７１】
次に、本発明の内部特性分布計測装置について説明する。図４に本発明の装置の一実施形態の模式図を示す。
【００７２】
図４に示す装置は、１２個の光ファイバーホルダー１〜１２（以下、場合により「光ファイバーホルダー群」と総称する）を備えており、各光ファイバーホルダー１〜１２は散乱吸収体ＳＭ(scattering medium) の一断面の周りに等間隔に配置されており（図４に示す装置においては、散乱吸収体ＳＭの中心から３０度間隔で放射線状に伸びる各線上に配置されている）、時計回りに１〜１２の番号が割り振られている。
【００７３】
光ファイバーホルダー１〜１２はそれぞれ光入射ファイバー１ａ〜１２ａと光検出ファイバー１ｂ〜１２ｂとを備える。光入射ファイバー１ａ〜１２ａと光検出ファイバー１ｂ〜１２ｂとは図４に示すように並列に束ねられている構造でもかまわないが、図５（ａ）に示すように光入射ファイバー１ａの回りを複数の光検出ファイバー１ｂ（バンドルファイバー）が包むように束ねられている構造、あるいは図５（ｂ）に示すように光入射ファイバー１ａと光検出ファイバー１ｂとを光ファイバーホルダー内で光カプラ１ｃでカップリングさせている構造であってもよい。図５（ａ）や図５（ｂ）に示すような構造を採用すると、散乱吸収体ＳＭの周囲にあたるファイバー端面は１つになり、両ファイバーを上下２段あるいは左右２列に並べる場合よりも光入射ファイバー端と光検出ファイバー端との位置ずれをなくすことができ、誤差を軽減できる傾向がある。
【００７４】
光入射ファイバー１ａ〜１２ａには波長選択器２０を介して光源３０が光学的に接続されている。そして、光源３０から発せられた光は波長選択器２０で波長選択され、光ファイバーホルダー１〜１２を通して計測対象である散乱吸収体ＳＭの表面に入射される。光源３０には、発光ダイオード、レーザーダイオード、Ｈｅ−Ｎｅレーザー等種々のものが使用できる。また、光源３０は、パルス光や方形波光、またはそれらの変調光を発生するものでもよい。
【００７５】
本実施形態において使用する光源３０は単波長の光（測定光）を発生するものであってもよいが、２波長以上の光（測定光）を発生可能なものであることが好ましい。計測に使用する光の波長は、計測対象に応じて適宜選択される。一般に生体では、ヘモグロビン等の吸収特性から７００ｎｍ以上の光、特に好ましくは可視光または近赤外光、を使用することが好ましい。例えば酸化ヘモグロビンと還元ヘモグロビンを対象とする場合は、図６に示すようにそれらの吸収係数が異なるから、適宜に波長を選択して使用することによってこれらを分離して計測することが可能となる。
【００７６】
光検出ファイバー１ｂ〜１２ｂには検出器４０が光学的に接続されている。そして、散乱吸収体ＳＭ中で散乱されつつ透過した光（測定光）は光ファイバーホルダー１〜１２の光検出ファイバー１ｂ〜１２ｂを介して検出器４０に導かれ、光検出器４０で受光信号を検出信号（電気信号）に変換して増幅し、それぞれに対応する検出信号が出力される。光検出器４０は光電子増倍管のほか、光電管、フォトダイオード、アバランシェフォトダイオード、ＰＩＮフォトダイオード等、あらゆる種類の光検出器を使用することができる。光検出器４０の選択に際しては、使用される測定光の波長の光が検出できる分光感度特性をもっていれば良い。また、光信号が微弱であるときは高感度あるいは高利得の光検出器を使用することが好ましい。なお、光検出ファイバー１ｂ〜１２ｂおよび光検出器４０の受光面以外の場所は、光を吸収あるいは遮光する構造にしておくことが望ましい。また、散乱吸収体ＳＭの内部を拡散伝搬した光が複数の波長の光を含む場合には、光検出器４０と散乱吸収体ＳＭとの間に波長選択フィルタ（図示せず）を適宜配置してもよい。
【００７７】
光源３０および検出器４０には制御部５０が接続されており、入射あるいは受光の際に使用される光ファイバーホルダー１〜１２の選択は制御部５０で行なわれる。すなわち、制御部５０によって光入射ファイバーから順次（例えば１ａ→２ａ→３ａ→・・・→１２ａ）一定時間間隔で散乱吸収体ＳＭ中に測定光が入射されるように制御され、それに同期して測定光が入射された光入射ファイバーに対して所定の位置関係にある光検出ファイバーから測定光が検出されるように制御される。本実施形態においては、測定光が入射された光入射ファイバーと異なる位置にある全ての光検出ファイバー（例えば光入射ファイバーが１ａの場合は光検出ファイバー２ｂ〜１２ｂ）から測定光がそれぞれ検出されるが、特にかかる組み合わせに限定されるものではない。
【００７８】
また、複数の波長を有する測定光を使用する場合は、入射される測定光の波長も制御部５０によって制御される。具体的な手法としては、異なる波長の光を時分割で入射させて使用する手法と、後述する異なる波長の光を同時に含む光を使用する手法とがある。具体的な波長選択手段としては、ミラーを用いた光ビーム切り換え器、フィルターを用いた波長切り換え器、光スイッチを用いた光切り換え器等がある（図７（ａ））。
【００７９】
上記の光入射ファイバー１ａ〜１２ａ、波長選択器２０、光源３０および制御部５０が本発明にかかる光入射手段を構成し、上記の光検出ファイバー１ｂ〜１２ｂ、検出器４０および制御部５０が本発明にかかる光検出手段を構成する。
【００８０】
制御部５０には演算処理部（例えばＣＰＵ）６０が電気的に接続されており、演算処理部６０には記憶部（例えばハードディスク、フレキシブルディスク）７０および表示部（例えばディスプレイ、プリンタ）８０が電気的に接続されている。そして、検出器４０から出力された検出信号は制御部５０を介して演算処理部６０へ導かれる。
【００８１】
上記の演算処理部６０および記憶部７０が、本発明にかかる測定値取得手段と、基準値算出手段と、内部特性変化量算出手段と、平均吸収および散乱係数検出手段と、広がり関数選択手段と、吸収係数絶対値算出手段と、吸収成分濃度算出手段と、等価散乱係数絶対値算出手段と、屈折率絶対値算出手段とを構成し、上記の表示部８０が画像表示手段を構成する。これらの本発明にかかる諸手段については、図８に示す本発明の方法の一実施形態のフローチャートに基づいて以下に詳細に説明する。
【００８２】
（１）図８に示す方法においては、先ず、以下に説明するようにして光ＣＴによる測定データ（Ｉd1｛Ｍ，ｍ｝））が取得される（Ｓ１００）。なお、Ｍは光入射ファイバーの番号を示し、ｍは光検出ファイバーの番号を示す。
【００８３】
すなわち、光入射ファイバー１ａ〜１２ａから順次散乱吸収体ＳＭ中に測定光を入射し、散乱吸収体ＳＭ中で散乱されつつ透過した各測定光をその測定光が入射された光入射ファイバーと異なる位置にある全ての光検出ファイバー（例えば光入射ファイバーが１ａの場合は光検出ファイバー２ｂ〜１２ｂ）から順次あるいは同時に検出する。なお、各測定光を同時に検出する場合は、光検出器４０が光検出ファイバー数分必要である。
【００８４】
そして、各光検出ファイバーで検出された各測定光に基づいて光検出器４０から検出信号が発せられる。そして、これらの検出信号はそれぞれ演算処理部６０において処理され、検出された測定光のそれぞれの検出光量に比例した測定値に変換されて、得られた測定値は記憶部７０に一次的に保存される。具体的には、演算処理部６０において、光源３０から発せられた光発生に同期した信号を利用して、検出信号に対する時間域での積分演算を行い、検出光量に比例した測定値を求める。ただし、パルス光等を利用する場合には、この同期信号を省略することができる。この種の演算処理は、演算処理手段に組み込んだマイクロコンピュータ等で高速に実行することができる。また、演算処理部６０において、測定値を平均化フィルタリングや最小二乗フィッティング等を利用して修正してもよい。
【００８５】
２）次に、演算処理部６０において位置関係が相対的に同じである光入射ファイバーと光検出ファイバーとの複数の組み合わせによって求められた複数の測定値が抽出され、それらの測定値の平均値である基準値（Ｉd0｛Ｍ，ｍ｝））が算出される（Ｓ１１０）。
【００８６】
すなわち、演算処理部６０では、吸収係数変化量等を求めるための基準値を得るために、光入射−光検出の位置関係が相対的に同じである光入射−光検出位置対毎に測定値の平均値を求める。
【００８７】
図４に基づいて具体的に説明すると、例えば光入射ファイバーと散乱吸収体ＳＭの中心と光検出ファイバーとのなす角度が１８０度の位置関係にある場合、相対的に同じ光入射−光検出位置の組み合わせは、（光入射ホルダー番号，光検出ホルダー番号）で表わすと(1,7)、(2,8)、(3,9)、(4,10)、(5,11)、(6,12)となる。なお、光の相反定理が成り立たない状態である場合は光入射−光検出位置の反対の組み合わせも考慮に入れる必要性がある。
【００８８】
各々の測定値をI(1,7)、I(2,8)、I(3,9)、I(4,10) 、I(5,11) 、I(6,12) とした場合、これらの平均値は下式となる。
【００８９】
I(ave_180)
＝｛I(1,7)＋I(2,8)＋I(3,9)＋I(4,10) ＋I(5,11) ＋I(6,12) ｝／６
このI(ave_180)を光入射−光検出の位置関係が１８０度の時の基準値とする。
【００９０】
同様にして、I(ave_150)、I(ave_120)、I(ave_90) 、I(ave_60) 、I(ave_30) も求め、これらの平均値を上記各位置関係における基準値として記憶部７０に一次的に保存する。
【００９１】
３）次に、本実施形態においては、上記各位置関係における基準値等に基づいて、平均吸収係数μa0および平均等価散乱係数μ'S0を光拡散理論等を利用して求めることができる（Ｓ１２０）。
【００９２】
すなわち、各角度毎の基準値から内部の吸収係数と等価散乱係数とを求め、さらにそれらの平均値を算出し、それらを散乱吸収体ＳＭ内部の平均吸収係数μa0と平均等価散乱係数μ'S0として記憶部７０に一次的に保存する。また、いずれか１つの角度の基準値、例えばI(ave_180)の値のみから散乱吸収体ＳＭ内部の平均吸収係数μa0と平均等価散乱係数μ'S0とを求めてもかまわない。
【００９３】
なお、上記基準値から散乱吸収体ＳＭ内部の吸収係数と等価散乱係数とを求める方法としては、例えば“Imaging diffusive media using time-independent and time-harmonic sources；dependence of image quality on imaging algorithms, target volume weight matrix, and view angles ”(Jenghwa Chang et.al., SPIE vol.2389）に記載されている方法がある。
【００９４】
４）次に、本実施形態においては、上記平均吸収係数μa0と平均等価散乱係数μ'S0に対応する広がり関数（Ｗθ）を選択する（Ｓ１３０）。
【００９５】
すなわち、予め用意されて記憶部７０に記憶してあった広がり関数の中から、上記で求めた平均吸収係数μa0および平均等価散乱係数μ'S0と合う広がり関数を選択する。この場合、実際の測定値から求めた吸収係数および等価散乱係数に基づいて広がり関数を選択するため、適当に仮定したそれらの値を用いるよりも誤差要因は除かれる。
【００９６】
なお、このような「広がり関数」とは、各ボクセルにおける光（測定光）の広がり方を示す関数を意味し、各ボクセルにおける実効光路長に関するいわゆる重み関数や各ボクセルにおける測定光に対する寄与度に関するいわゆる寄与関数を包含する概念である。本発明にかかる広がり関数は上記の重み関数あるいは寄与関数のいずれでもよい。このような広がり関数に関しては、例えば“A Perturbation Model for Imaging in Dense Scattering Media: Derivation and Evaluation of Imaging Operation”（H.L.Graber et.al., SPIE vol. IS11）、“Initial assessment of a simple system for frequency domein diffuse optical tomography”（B.W.Pogue et al.,Phys.Med.Biol.40(1995)p.1709-p.1729）、本発明者らによる特願平８−６６１９「光ＣＴ装置及び光ＣＴによる画像再構成方法」に記載されている。本実施形態においては特願平８−６６１９に記述された作成方法にしたがって、以下に示すような時間項の入っていない光拡散方程式を用いて広がり関数を予め用意する。
【００９７】
Δφ−μaＤ-1φ＝０
ここで、Ｄ＝1/{3(1-g)μs}＝1/3μ's、φ：光子密度、Ｄ：光拡散定数、μa：吸収係数、μ'S：等価散乱係数、ｇ：散乱吸収体による光子の散乱角の余弦の平均値である。また、屈折率分布を求める際に好適に用いられる時間項の入っている光拡散方程式を以下に示す。
【００９８】
【数１３】
【００９９】
ここで、Ｄ(r)＝1/{3(1-g)μs(r)}＝1/3μ's(r)、φ(r,t)：位置ｒ、時刻ｔにおける光子密度、Ｃ：媒質中の光速、Ｄ：光拡散定数、μa：吸収係数、Ｓ(r,t)：光源、μ'S：等価散乱係数、ｔ：時間、ｒ：位置、ｇ：散乱吸収体による光子の散乱角の余弦の平均値である。なお、真空中の光速をＣ’、被測定物の屈折率をｎとするとＣ＝Ｃ’ｎと表わすことができる。
【０１００】
また、平均吸収係数μa0及び平均等価散乱係数μ'S0と合う広がり関数とは、具体的には、測定対象物の平均吸収係数及び平均等価散乱係数と同じ平均吸収係数及び平均等価散乱係数を有しかつ測定対象物と同じ形状を有する対象物に対して光入射−光検出位置の相対的関係を実測の場合と同様にした場合に求められるであろう光の広がり方を示す関数のことであり、平均吸収係数μa0及び平均等価散乱係数μ'S0等に基づいて選択される。
【０１０１】
また、記憶部７０に対象物を複数のブロック（ボクセル）に分けた際に生じる歪みを補正するための補正項を記憶しておいてもよく、その場合は演算処理部６０で前記測定値および／または前記基準値を補正することが可能である（Ｓ１４０）。このようなボクセルに関する補正とは、例えば光入射位置ａと光検出位置ｂ，ｃとの間の距離が同じにも拘らずボクセルの切り方によってボクセル上でのトータルの距離が異なるような場合に、その相違分を測定値及び／又は基準値の補正項として利用することである。
【０１０２】
５）続いて、演算処理部６０において前記複数の組み合わせによって求められた複数の測定値と、前記基準値と、前記広がり関数とを用いて、複数の領域に分割された前記各領域における吸収係数の変化量Δμa を算出し（Ｓ１５０）、出力する（Ｓ１６０）。
【０１０３】
すなわち、前記各角度毎の基準値と前記測定値と前記広がり関数とを用いて吸収係数の変化量を求める。その際に成り立つ関係式は、前記式（４）に対応させると、例えば光入射−光検出の位置関係が１８０度の場合は基準値Ｉd0はI(ave_180)となり、測定値Ｉd1はI(1,7)、I(2,8)、I(3,9)、I(4,10) 、I(5,11) 、I(6,12) となる。この時、基準値Ｉd0の吸収係数は散乱吸収体ＳＭ内部の平均の吸収係数である。さらに、光入射−光検出の位置関係が１８０度の時の広がり関数をＷθとすると下記の式（4-1）〜式（4-6）が成り立ち、これらの連立方程式（前記式（８）の連立方程式）を各位置関係毎に立てて（つまり未知数分と同じ数の連立方程式を立てて）解くことによって各領域における吸収係数の変化量Δμa が算出される。
【０１０４】
【数１４】
【０１０５】
【数１５】
【０１０６】
【数１６】
【０１０７】
【数１７】
【０１０８】
【数１８】
【０１０９】
【数１９】
【０１１０】
なお、散乱吸収体ＳＭ内部の吸収係数の変化量の空間的分布や濃度変化の空間的分布を求める場合、このように成り立った関係式を散乱吸収体ＳＭの内部を分割したボクセル（ボックスセル）の個数分と同じ数の連立方程式を立てて解けばよい。本実施形態においては、共役勾配法を用いた。また、式数がボクセルの個数分より少ない場合や多い場合であっても、特異値分解法等を用いれば特異な問題を非特異な問題に変えられるため、内部特性の分布が求められる。
【０１１１】
このようにして求められた各領域における吸収係数の変化量Δμa に基づいて測定対象物内部における吸収係数変化量に関する吸収分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ１７０）。
【０１１２】
なお、上記に説明したような演算処理部６０におけるΔμa の算出方法から吸収分布を求めて画像を表示する方法としては他にも様々な方法が知られている。このような方法は、例えば“Optical Back Projection Tomography in Heterogeneous Diffusive Media ”（S.B.Cloak et al., in Advances in Optical Imaging and Photon Migration, 1996 Technical Digest; Optical Society of America, Washington DC, 1996, pp.147-149 ）、“Back-projection image reconstruction using photon density wave in tissues ”（S.A.Walker et al., SPIE vol.2389, pp.350, 1995）、“in Photon Propagation in Tissues”（S.Cloak et al., SPIE vol.2626, 1995 ）、“Optical tomography by the temporally extrapolated absorbance method ”（Ichiro Oda et al., APPLIED OPTICS, vol.35, No.01, 1996）等に記載されている。これらの方法はＸ線ＣＴで知られているバックプロジェクション（Back Projection ）法もしくはそれらを変形した方法であり、上記アルゴリズムの変わりに画像再構成をさせる方法である。
【０１１３】
また、吸収成分の既知のモル吸収係数を用いて上記各領域における吸収係数の変化量Δμa から各領域における吸収成分の濃度差の絶対値を算出することが可能であり（Ｓ１８０）、このようにして求められた各領域における吸収成分の濃度差の絶対値に基づいて測定対象物内部における吸収成分濃度差に関する分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ１９０）。
【０１１４】
更に、上記の各領域における吸収係数の変化量Δμa および前記平均吸収係数μa0を用いて上記各領域における吸収係数の絶対値μa を算出することが可能であり（Ｓ２００）、このようにして求められた各領域における吸収係数の絶対値μa に基づいて測定対象物内部における吸収係数絶対値に関する分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ２１０）。
【０１１５】
更にまた、吸収成分の既知のモル吸収係数を用いて上記の各領域における吸収係数の絶対値μa から各領域における吸収成分の濃度を算出することが可能であり（Ｓ２２０）、このようにして求められた各領域における吸収成分の濃度に基づいて測定対象物内部における吸収成分濃度に関する分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ２３０）。
【０１１６】
また、散乱吸収体ＳＭが少なくとも２つの吸収成分、例えば酸化および還元ヘモグロビン、を含有する場合は、それらの吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有する測定光を使用し、各波長を有する測定光のそれぞれに関して前記測定値、前記基準値を求め、それらに基づいて各波長を有する測定光のそれぞれに関して前記吸収係数変化量、前記吸収係数絶対値を求めることによって、各吸収成分の濃度分布が求められる。
【０１１７】
以下に、上記２波長分光法を利用したヘモグロビンの濃度の計測について説明する。
【０１１８】
哺乳類の脳における吸収成分の主なものは、水、チトクローム（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)を使用すれば上記の式はさらに簡単になる。
【０１２５】
以上、本発明の好適な実施形態について説明したが、本発明は勿論上記実施形態に限定されるものではない。
【０１２６】
すなわち、上記実施形態においては内部特性として吸収係数を求めているが、前述のように本発明は等価散乱係数の計測にも適用可能であり、吸収係数および等価散乱係数を求める一実施形態のフローチャートを図９に示す。
【０１２７】
図９に示す方法においては、吸収係数の測定に関しては図８に示す方法と同様であるが、平均吸収係数μa0と平均等価散乱係数μ'S0に対応する広がり関数（Ｗθ）を選択する際に、吸収係数用の広がり関数（Ｗμa,j）と等価散乱係数用の広がり関数（Ｗμ'S,j）とを選択しておくことが好ましい（Ｓ１３０）。
【０１２８】
そして、演算処理部６０において前記複数の組み合わせによって求められた複数の測定値と、前記基準値と、前記広がり関数とを用いて、複数の領域に分割された前記各領域における等価散乱係数の変化量Δμ'Sを算出し（Ｓ２４０）、出力する（Ｓ２５０）。すなわち、前述の式（４’）に基づいて、前記各角度毎の基準値と前記測定値と前記広がり関数と吸収係数の変化量とを用いて等価散乱係数の変化量を求める。より具体的には、前述の式（4-1）〜式（4-6）と同様に式（４’）に基づいて連立方程式が成り立ち、これらの連立方程式を各位置関係毎に立てて解くことによって各領域における等価散乱係数の変化量Δμ'Sが算出される。
【０１２９】
このようにして求められた各領域における等価散乱係数の変化量Δμ'Sに基づいて測定対象物内部における等価散乱係数変化量分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ２６０）。
【０１３０】
更に、上記の各領域における等価散乱係数の変化量Δμ'Sおよび前記平均等価散乱係数μ'S0を用いて上記各領域における等価散乱係数の絶対値μ'Sを算出することが可能であり（Ｓ２７０）、このようにして求められた各領域における等価散乱係数の絶対値μ'Sに基づいて測定対象物内部における等価散乱係数絶対値に関する分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ２８０）。
【０１３１】
また、本発明は屈折率の計測にも適用可能であり、吸収係数、等価散乱係数および屈折率を求める一実施形態のフローチャートを図１０に示す。
【０１３２】
図１０に示す方法においては、吸収係数および等価散乱係数の測定に関しては図８および図９に示す方法と同様であるが、各位置関係における基準値等に基づいて、平均吸収係数μa0および平均等価散乱係数μ'S0を求める際に、平均屈折率ｎ0も併せて求める（Ｓ１２０）。但し、平均屈折率ｎ0としては水の屈折率（1.33）を用いてもよい。また、広がり関数（Ｗθ）を選択する際に、吸収係数用の広がり関数（Ｗμa,j）と等価散乱係数用の広がり関数（Ｗμ'S,j）と屈折率用の広がり関数（Ｗn,j）とを選択しておくことが好ましい（Ｓ１３０）。
【０１３３】
そして、演算処理部６０において前記複数の組み合わせによって求められた複数の測定値と、前記基準値と、前記広がり関数とを用いて、複数の領域に分割された前記各領域における屈折率の変化量Δｎを算出し（Ｓ２９０）、出力する（Ｓ３００）。すなわち、前述の式（４''）に基づいて、前記各角度毎の基準値と前記測定値と前記広がり関数と吸収係数の変化量と等価散乱係数の変化量とを用いて屈折率の変化量を求める。より具体的には、前述の式（4-1）〜式（4-6）と同様に式（４''）に基づいて連立方程式が成り立ち、これらの連立方程式を各位置関係毎に立てて解くことによって各領域における屈折率の変化量Δｎが算出される。
【０１３４】
このようにして求められた各領域における屈折率の変化量Δｎに基づいて測定対象物内部における屈折率変化量分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ３１０）。
【０１３５】
更に、上記の各領域における屈折率の変化量Δｎおよび前記平均屈折率ｎ0を用いて上記各領域における屈折率の絶対値ｎを算出することが可能であり（Ｓ３２０）、このようにして求められた各領域における屈折率の絶対値ｎに基づいて測定対象物内部における屈折率絶対値に関する分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ３３０）。
【０１３６】
そして、このように屈折率に関する分布が求められると、血糖値の分布を求めることが可能となる。屈折率の変化により血糖値を検出する方法としては、例えば“Possible correlation between blood glucose concentration and the reduced scattering coefficient of tissues in the near infrared”（John S.Maier et al., OPTICS LETTERS vol.19, No.24, December 15, 1994）に記載の方法がある。生体組織のグルコース濃度は細胞外液の屈折率に大きな影響を与えており、組織中の等価散乱係数は細胞外液徒細胞との間の屈折率差に大きく依存するため、細胞外液の屈折率が変化して検出光に影響を及ぼす結果となる。したがって、検出光に基づいて屈折率分布を求めることによって組織内部の血糖値の分布を求めることが可能となるのである。
【０１３７】
また、上記実施形態においては複数の光入射及び光検出位置を散乱吸収体の一断面の周囲に配置しているが、図１１、図１２に示すように光入射及び／又は光検出位置（これらをＰで示す）を立体的に配置してもよい。すなわち、測定対象物として頭部や乳房を想定した場合は光入射及び／又は光検出位置（Ｐ）を図１１に示すように配置してもよく、測定対象物として腕部、脚部、胸部、乳房（圧迫時）を想定した場合は光入射及び／又は光検出位置（Ｐ）を図１２に示すように配置してもよい。
【０１３８】
また、上記実施形態においては時間積分法による光量測定値を測定値として用いているが、本発明に使用可能な測定値はこれに限定されるものではなく、例えば測定光の位相差（または位相遅れ）、振幅であってもよい。また、演算処理部６０において測定値を取得する具体的な手法も所望の測定値に応じて適宜選択され、例えば位相変調法による位相差及び／又は振幅測定、時間分解分光法による時間分解波形測定等の手段を採用してもよい。
【０１３９】
また、上記実施形態においては散乱吸収体ＳＭ内部の平均吸収係数および平均等価散乱係数を本発明にかかる光ＣＴ装置自身で得られたデータから求めているが、図１３〜図１５に示すように散乱吸収体ＳＭ内部の平均吸収係数および平均等価散乱係数を別の装置で求め（Ｓ１２０ａ）、それらに基づいて広がり関数を選択してもかまわない（Ｓ１３０ａ）。なお、図１３〜図１５における上記以外の諸工程は図８〜図１０における諸工程にそれぞれ対応する。この場合の利点は、例えば光ＣＴ装置で得られるデータはＣＷ（連続光）で測定し、平均吸収係数および平均等価散乱係数を求める装置でのみパルス光や変調光を用いればよくなるため、光ＣＴ装置のシステムの構成が簡単になる。なお、別の装置で平均吸収係数および平均等価散乱係数を求める手法は位相変調法や時間分解分光法であってもよい。
【０１４０】
このように測定対象物の内部の光学パラメータの分布を均一とみなして平均等価散乱係数μ's0と平均吸収係数μa0を測定する方法は、時間分解分光法については例えば“Devlopment of Time Resolved Spectroscopy System for Quantitative None-invasive Tissue Measurement ”M.Miwa et.al., SPIE vol. 2389 、位相変調法については例えば特開平６−２２１９１３号公報に記載されている。また、光ＣＴ装置に上記の手法を取り入れれば前記測定値の取得と同時に計算が可能である。
【０１４１】
また、光を生体等の散乱吸収体に入射する手段としては、図４および図１６（ｂ）に示した光ファイバーを使用する方法の他に、集光レンズ（図１６（ａ））、ピンホール（図１６（ｃ））を利用する方法、胃カメラのように体内から入射する方法（図１６（ｄ））等であってもよい。
【０１４２】
なお、生体試料等の散乱吸収体では平均拡散長が２ｍｍ程度であるため、入射光は約２ｍｍ直進するまでに散乱し、光の方向性がなくなる。したがって、数センチメートル以上の厚みの散乱吸収体では平均拡散長の影響を無視することができるのでスポット状に光を入射すれば良い。また、太いビーム状の光を散乱吸収体に入射してもよい。この場合には、複数のスポット状光源が並んでいると考えればよい。
【０１４３】
また、光入射ファイバーおよび光検出ファイバーと散乱吸収体ＳＭ表面との間の空間は、図４に示す実施形態では微小になっている。しかし実際には、これを大きくして、この空間に計測対象である散乱吸収体ＳＭとほぼ等しい屈折率および等価散乱係数をもつ液状体やゼリー状物体（以下、インターフェース材と呼ぶ）を満たしておいても良い。つまり、光はこのインターフェース材中を拡散伝搬して計測対象に入射するから何ら問題は生じない。また、散乱吸収体ＳＭの表面反射が問題になるときには、インターフェース材を適宜に選択することによって、表面反射等の影響を低減することができる。
【０１４４】
更に、上記実施形態においては異なる波長の光を時分割で入射させて使用する態様について説明したが、異なる波長の光を光混合器３５によって同軸状のビームにして、光入射点の直前に設けた波長選択フィルタ２０で波長選択して各波長の光を散乱吸収体に入射する方法、あるいはそのまま前記ビームを並行に散乱吸収体に入射する方法であってもよい（図７（ｂ））。但し、後者の場合は、図１７に示すように、光検出器４０ａ〜４０ｃの直前に波長選択フィルタ２５ａ〜２５ｃを設けて検出光を波長選択する必要がある。
【０１４５】
また、散乱吸収体の内部を拡散伝搬した光を受光して検出する手段としては、図４および図１８（ｂ）に示した光ファイバーを用いる方法以外に、直接検出する方法（図１８（ａ））、レンズを用いる方法（図１８（ｃ））等がある。
【０１４６】
また、光検出器４０で得られる信号を低雑音で増幅する必要がある場合には、狭帯域アンプ（図１９（ａ））、ロックインアンプ（図１９（ｂ））等を利用することができる。ロックインアンプを使用する場合には、参照信号として前記同期信号を利用する。この方法は、方形波光やパルス光を用いてダイナミックレンジの高い計測を行うときに有効である。
【０１４７】
さらに、上記実施形態においては複数の光入射ファイバーおよび光検出ファイバーを散乱吸収体の周囲に配置しておき、光入射および光検出に使用するファイバーをそれぞれ順次変えていくことによって光入射位置および光検出位置を移動させていたが、散乱吸収体に対する光入射位置と光検出位置とを同期して走査させてもよい。このようにすれば、散乱吸収体の各部の内部情報を求めて、フレームメモリーに蓄積し、これをテレビ方式で読み出すことによって内部情報の分布を示す画像が得られる。また、異なる時刻に計測すれば、内部情報の時間的変化を計測することができる。前記記憶部７０はこのようにして得た内部情報を記憶する機能をもち、表示部８０はこれらの途中経過や結果を表示するものである。この際、これらの演算処理は、メモリ７０、ディスプレイ８０等を備えるコンピュータ装置６０によって高速に実行することができる。
【０１４８】
さらに別な方法として、特開平６−２２１９１３号公報、特開平６−１２９９８４号公報に記載されているように、測定対象物の周囲に真円状のホルダーを取り付けて測定すれば、人の頭部を測定する際も個体差はそこで吸収され、実測の系をモデル系に近い状態で測定でき、精度が向上する傾向にある。
【０１４９】
また、上記実施形態においては既知の等価散乱係数より吸収係数の空間分布を算出したが、同一手法により、既知の吸収係数から等価散乱係数の空間的分布を求めることもできるため、したがって、今回提案する手法により吸収係数と等価散乱係数の両方の空間的分布を求めることが可能である。
【０１５０】
上記実施形態においては内部特性を求める式として光拡散方程式を用いたが、必ずしもこれに限定されるものではない。例えば、散乱吸収体内部の光吸収が伝播距離の関数として表わされることから導かれる式や、実験から得られる検出光と内部特性との関係式等を用いてもよい。
【０１５１】
また、病気や体の調子と検出光との相関が求まれば、本発明によって検出光から直接有用な情報を取得することが可能となる。例えば、検出光の変化が組織の構造変化と相関があれば、その相関を利用して検出光から構造変化を求めることができる。
【０１５２】
【実施例】
以下の実施例では、対象物がもつ全ての内部特性（吸収係数、等価散乱係数、屈折率等）の中で、吸収係数のみを変化させた場合における吸収係数の画像化の例を示す。
【０１５３】
本発明の有効性を実証するために、光入射ファイバーの配置間隔を２０度とし、光検出ファイバーの配置間隔を１０度とした以外は図４に示す装置と同様の装置を用いて以下の手順で実験を行なった。なお、本実施例で使用したファントムは高さ方向に同一形状であるため（ｚ軸方向に対称性を持つため）、３次元（立体）の問題を２次元（断面）の問題に次元を落として考えることができる。
【０１５４】
すなわち、波長８００ｎｍ、出力５０ｍＷのＣＷ光のレーザーを光入射ファイバーから図２０（ａ）、（ｂ）に示すファントムに入射し、ファントムを透過あるいは散乱反射した光を光検出ファイバーで検出し、検出器に導いた。なお、使用したファントムの仕様は以下の通りである。
【０１５５】
【０１５６】
検出された光信号は光電子増倍管にて検出信号に変換し、その検出信号をカウンターで読み取って１０秒間加算した。加算された値はＧＰＩＢ（計測インターフェース・バス）でコンピュータに送った。
【０１５７】
ここまでの測定を１回とし、光入射ファイバーを図２０に示すＡ〜Ｒまで２０度づつ、光検出ファイバーを１０度づつ左回りに移動させてデータをとった。なお、１つの光入射位置に対する光検出位置は、それぞれの光入射位置を基準にして０度から９０度までの１０箇所とした。
【０１５８】
実験結果を以下に示す。
【０１５９】
１）別に作製した吸収体の含まれていないファントムで基準値Ｉd0を取得し、画像再構成を行った結果を図２２（ａ）に示す（背景は黒で示し、図２２（ｂ）は吸収係数の指標を示す）。図２２（ａ）に示す画像では複数個の像が見られ、かつ位置ずれが生じている。なお、図２２（ａ）中に記載された円は、期待される像の位置を示す。
【０１６０】
２）本発明にしたがって光入射−光検出の位置関係が相対的に同じである光入射−光検出位置の組み合わせ毎に測定値の平均値を求め、それを基準値Ｉd0をとして画像再構成を行った結果を図２３（ａ）に示す（背景は黒で示し、図２３（ｂ）は吸収係数の指標を示す）。本発明によると、図２３（ａ）に示した画像に見られるように、従来の方法によって得た図２２（ａ）に示す画像比べて位置ずれと定量性が改善されることが明らかである。なお、図２３（ａ）中に記載された円は、期待される像の位置を示す。
【０１６１】
【発明の効果】
本発明によれば、従来は必要であった物理モデルあるいはシミュレーションモデルから基準値を求めることなく、また測定対象物中の一成分に対して複数の波長を有する光を用いることなく、測定対象物についての測定値から直接的に基準値を求めることが可能となり、その基準値に基づいて測定対象物における内部特性分布（例えば吸収係数変化量分布、吸収係数分布、吸収成分濃度分布、等価散乱係数変化量分布、等価散乱係数分布、屈折率変化量分布、屈折率分布）を求めて画像化することが可能となる。そのため、本発明によれば、実際の測定対象物と物理モデルまたはシミュレーションモデルとの差や測定対象物の個体差等に起因する誤差の発生を避けることでき、信頼性の高いすなわち高精度の測定が可能となる。また、本発明によれば、物理モデル等を用いて基準値を予め求める手間が省かれ、測定時間の短縮化が可能となる。さらに、本発明によれば、複数の波長間での平均光路長分布および減衰光量（散乱等の影響で減衰する光量）は同じであるといった仮定に基づく誤差の発生が防止され、測定精度が向上する。
【図面の簡単な説明】
【図１】本発明の動作原理の説明図である。
【図２】吸収が均一な散乱吸収体のモデルを示す模式図である。
【図３】吸収が不均一な散乱吸収体のモデルを示す模式図である。
【図４】本発明の内部特性分布計測装置の一例を示す模式図である。
【図５】（ａ）および（ｂ）はそれぞれ光入射ファイバーの一例を示す斜視図および模式図である。
【図６】ヘモグロビンおよびミオグロビンの吸収スペクトルを示すグラフである。
【図７】（ａ）および（ｂ）はそれぞれ光入射手段の一例を示す模式図である。
【図８】本発明の内部特性分布計測方法の一例を示すフローチャートである。
【図９】本発明の内部特性分布計測方法の他の例を示すフローチャートである。
【図１０】本発明の内部特性分布計測方法の更に他の例を示すフローチャートである。
【図１１】本発明にかかる光入射及び／又は光検出位置の配置の一例を示す模式図である。
【図１２】本発明にかかる光入射及び／又は光検出位置の配置の他の例を示す模式図である。
【図１３】本発明の内部特性分布計測方法の更に他の例を示すフローチャートである。
【図１４】本発明の内部特性分布計測方法の更に他の例を示すフローチャートである。
【図１５】本発明の内部特性分布計測方法の更に他の例を示すフローチャートである。
【図１６】（ａ）〜（ｄ）はそれぞれ散乱吸収体への光入射方法の一例を示す模式図である。
【図１７】光検出手段の一例を示す模式図である。
【図１８】（ａ）〜（ｃ）はそれぞれ光検出方法の一例を示す模式図である。
【図１９】（ａ）および（ｂ）はそれぞれ検出信号の低雑音増幅方法の一例を示す模式図である。
【図２０】（ａ）および（ｂ）はそれぞれ実施例で使用したファントムの斜視図および上面図である。
【図２１】実施例における光入射位置と光検出位置との関係を説明するための補助図である。
【図２２】（ａ）及び（ｂ）は従来の方法で画像再構成を行った結果、ディスプレー上に表示された中間調画像を示す写真である。
【図２３】（ａ）及び（ｂ）は本発明の方法で画像再構成を行った結果、ディスプレー上に表示された中間調画像を示す写真である。
【符号の説明】
１〜１２…光ファイバーホルダー、１ａ〜１２ａ…光入射ファイバー、１ｂ〜１２ｂ…光検出ファイバー、２０…波長選択器、３０…光源、４０…光検出器、５０…制御部、６０…演算処理部、７０…記憶部、８０…表示部、ＳＭ…散乱吸収体、Ｐ…光入射及び／又は光検出位置。
【発明の属する技術分野】
本発明は、測定対象物の内部特性分布を計測するための方法およびそのための装置に関する。より詳細には、本発明は、光入射位置および光検出位置を測定対象物の表面に沿って移動させて断層像を得る光ＣＴ（コンピュータトモグラフィー）装置等に適用可能な内部特性分布計測方法並びにそのための装置に関する。
【０００２】
【従来の技術】
散乱吸収体である対象物の表面における１つの光入射位置から測定光を入射し、対象物によって散乱されつつ透過した測定光を前記対象物の表面における複数の光検出位置で受光するとともに、前記対象物の表面に添って光入射位置および光検出位置を移動させていき、散乱吸収体における内部特性の分布を求める光ＣＴ装置においては、例えばその内部の吸収係数の分布を求める方法として以下の方法が知られている。すなわち、“Imaging of Multiple Targets in Dense Scattering Media”（H.L.Graber, J.Chang, R.L.Barbour, SPIE vol.2570, p.219-p.234）、“Imaging diffusive media using time-independent and time-harmonic sources；dependence of image quality on imaging algorithms, target volume weight matrix, and view angles ”（Jenghwa Chang et.al., SPIE vol.2389）等に記載された方法がある。
【０００３】
このような従来の方法における基本的な画像化原理は、測定対象物内部を便宜的に複数のボクセルに区切り、対象物表面のある点から入射した光が測定対象物内部を透過して同表面の他の点で受光された際に、各ボクセル毎の吸収係数といった特定の内部特性に注目した場合の受光される光に対する寄与度を示す関数（ここでは便宜的に「広がり関数」という）と、受光された光との関係式を用いることにある。なお、ここでいうボクセルとは、測定対象物を複数の領域に分割した各領域（ボックスセル）のことをいう。
【０００４】
しかしながら、上記従来の方法においては、測定対象物とは別に吸収の無いファントムを用意し、それを用いて基準となる検出光量を測定し、さらにその状態の広がり関数を用いて目的とする散乱吸収体内部の吸収係数分布を求めていた。そして、このような方法を用いて画像化を行う場合、測定対象物と同一又は類似した形状を有しかつ既知の内部特性を持つように作られたファントムモデル（物理モデル）あるいはシミュレーションモデルを想定し、このようなモデルから得たデータを画像化の計算の際の基準値として用いる必要があった。そのため、これら従来の方法では、実際の測定対象物と物理モデルまたはシミュレーションモデルとの差や測定対象物の個体差等に起因する誤差を避けることができず、特に生体のように複雑な構造を有する測定対象への適用は非常に困難であった。
【０００５】
一方、ファントムを用いないで吸収物質濃度の空間分布を求める方法として、特開平８−２９３２９号公報に記載の方法がある。しかしながら、同公報記載の方法においては、測定対象物中の吸収成分が一種類の場合でも複数の波長を有する光を用いる必要があり、これら複数の波長間での平均光路長分布および減衰光量（散乱等の影響で減衰する光量）は同じであると仮定して吸収物質濃度の空間分布を求めていた。また、この方法においては、吸収物質を含まない仮想被検体を想定し、その仮想被検体内における平均光路長を用いて吸収物質濃度の空間分布を求めており、吸収による光路長の変化は考慮されていなかった。そのため、上記公報に記載の方法によっても、得られる内部特性分布の信頼性は未だ充分なものではなかった。
【０００６】
【発明が解決しようとする課題】
そこで本発明は、上記従来の課題に鑑みてなされたものであり、従来は必要であった物理モデルあるいはシミュレーションモデルから基準値を求めることなく、また測定対象物中の一成分に対して複数の波長を有する光を用いることなく、測定対象物についての測定値から直接的に基準値を求めることを可能とし、その基準値に基づいて測定対象物における内部特性分布の信頼性の高いすなわち高精度の測定が可能な方法および装置を提供することを目的とする。
【０００７】
【課題を解決するための手段】
本発明者らは、上記目的を達成すべく鋭意研究した結果、測定対象物の表面にありかつその対象物内の一点（例えば対象物の中心）に対して相対的に同じ位置関係にある光入射位置と光検出位置との複数の組み合わせによって求められた複数の測定値の平均値を内部特性分布を求めるための基準値として用いることによって上記課題が解決されることを見出し、本発明に到達した。
【０００８】
すなわち、本発明の内部特性分布の計測方法は、
測定対象物の表面における複数の光入射位置から順次該対象物中に測定光を入射するステップと、
該対象物中を透過した測定光を、前記対象物の表面における複数の光検出位置のうちの少なくとも１つの光検出位置であってかつ測定されるべき測定光が入射された光入射位置に対して所定の位置関係にある光検出位置で順次あるいは同時に検出するステップと、
各光検出位置で検出された各測定光に基づいて、該測定光の所定パラメータの測定値を求めるステップと、
前記位置関係が相対的に同じである前記光入射位置と前記光検出位置との複数の組み合わせによって求められた複数の前記測定値を抽出し、該測定値の平均値を算出して該位置関係における基準値を得るステップと、
前記複数の組み合わせによって求められた前記複数の測定値と、前記基準値とを用いて、複数の領域に分割された前記対象物の各領域における所定内部特性の変化量を算出して該対象物における内部特性変化量分布を求めるステップと、
を具備する方法である。
【０００９】
また、本発明の内部特性分布の計測装置は、
測定対象物の表面における複数の光入射位置から順次該対象物中に測定光を入射する光入射手段と、
該対象物中を透過した測定光を、前記対象物の表面における複数の光検出位置のうちの少なくとも１つの光検出位置であってかつ測定されるべき測定光が入射された光入射位置に対して所定の位置関係にある光検出位置で順次あるいは同時に検出する光検出手段と、
各光検出位置で検出された各測定光に基づいて、該測定光の所定パラメータの測定値を求める測定値取得手段と、
前記位置関係が相対的に同じである前記光入射位置と前記光検出位置との複数の組み合わせによって求められた複数の前記測定値を抽出し、該測定値の平均値を算出して該位置関係における基準値を得る基準値算出手段と、
前記複数の組み合わせによって求められた前記複数の測定値と、前記基準値とを用いて、複数の領域に分割された前記対象物の各領域における所定内部特性の変化量を算出して該対象物における内部特性変化量分布を求める内部特性変化量算出手段と、
を具備する装置である。
【００１０】
本発明の方法および装置では、測定対象物の表面にありかつその対象物内の一点（例えば対象物の中心）に対して相対的に同じ位置関係にある光入射位置と光検出位置との複数の組み合わせによって求められた複数の測定値の平均値を、内部特性分布を求めるための基準値として用いる。すなわち、上記基準値と各測定値とを用いて後述する方程式を解くことよって、複数の領域に分割された前記対象物の各領域における内部特性の変化量（差分）が求められる。
【００１１】
このように本発明では、実際の測定で得られる測定値の平均値から基準値を求め、この基準値に基づいて内部特性の変化量が算出される。したがって、本発明では物理モデルやシミュレーションモデルから予め求めた基準値を用いていないため、測定対象物の個体差や、実際の測定対象物と物理モデルやシミュレーションモデルとの間に生じる条件の差等に起因する誤差が発生する余地がなく、信頼性の高い高精度の測定が可能である。さらに、本発明によれば、物理モデル等を用いて基準値を予め求める手間が省かれ、測定時間が短縮化される。
【００１２】
基準値を実際の測定で得られる測定値の平均値から求め、その基準値に基づいて内部特性の変化量等を求めるということは、例えば図１で説明すると、Ａ、Ｂ、Ｃの各値は０からの値が分からなくても平均値との差分が分かれば求められるということと同様である。本発明の動作原理はこのような原理に基づいて内部特性の差分や絶対値を求めることである。
【００１３】
また、本発明では測定対象物中の一成分に対して複数の波長を有する光を用いることなく測定対象物中の内部特性分布が求められるため、複数の波長間での平均光路長分布および減衰光量（散乱等の影響で減衰する光量）は同じであるといった仮定に基づく誤差の発生が防止され、測定精度が向上する。さらに、本発明によれば、測定対象物中の多成分を複数の波長を有する光を用いて解析する場合においても誤差の発生は同様に防止される。すなわち、本発明においては、現実の物体には存在する散乱係数の波長依存性の変化にも対応して各波長毎に散乱係数を求めているためである。
【００１４】
なお、本発明にかかる光入射位置と光検出位置との位置関係は、例えば測定対象物の中心を基準に、すなわち光入射位置および該中心を結ぶ線と光検出位置および該中心を結ぶ線との間の角度で規定され、かかる位置関係が相対的に同じであるということは例えば上記角度が同じであることをいう。
【００１５】
また、本発明にかかる測定値としては、測定光の測定対象物内部での散乱および吸収に関係する所定パラメータの測定値が好ましく、例えば測定光の光量、位相差（または位相遅れ）、振幅、時間分解波形等のパラメータの測定値が好適に用いられる。
【００１６】
本発明の方法および装置によって測定可能な内部特性としては、吸収係数、等価散乱係数、屈折率が挙げられ、これらのうちのいずれかの特性を単独で、あるいは複数の特性を同時または順次に求めてもよい。
【００１７】
先ず、本発明の方法および装置によって測定されるべき内部特性が吸収係数である場合について説明する。
【００１８】
この場合、本発明の方法は、前記対象物の平均の吸収係数および平均の等価散乱係数を求める（好ましくは前記基準値に基づいて求める）ステップと、前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数（吸収係数用の広がり関数）を選択するステップとを更に具備していることが好ましく、前記内部特性変化量分布を求めるステップにおいて前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における吸収係数の変化量を算出することができる。
【００１９】
また、本発明の装置は、前記対象物の平均の吸収係数および平均の等価散乱係数を求める（好ましくは前記基準値に基づいて求める）平均吸収および散乱係数検出手段と、前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数（吸収係数用の広がり関数）を選択する広がり関数選択手段とを更に具備していることが好ましく、前記内部特性変化量算出手段において前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における吸収係数の変化量を算出することができる。
【００２０】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての吸収係数および等価散乱係数に対応して選択された広がり関数に基づいて、前記の各領域における吸収係数の変化量（差分）が求められる。したがって、このような本発明の方法および装置によれば、吸収係数および／または等価散乱係数をゼロとして計算する場合に比べて、それらによる実効光路長の変化に基づく誤差の発生が充分に防止され、測定精度が向上する。
【００２１】
さらに、上記本発明の方法は、前記吸収係数の変化量および前記平均の吸収係数を用いて、前記各領域における吸収係数の絶対値を算出して前記対象物における吸収係数絶対値分布を求めるステップ、および／または前記吸収係数の絶対値を用いて、前記各領域における吸収成分の濃度を算出して前記対象物における吸収成分濃度分布を求めるステップを更に具備してもよい。
【００２２】
また、上記本発明の装置は、前記吸収係数の変化量および前記平均の吸収係数を用いて、前記各領域における吸収係数の絶対値を算出して前記対象物における吸収係数絶対値分布を求める吸収係数絶対値算出手段、および／または前記吸収係数の絶対値を用いて、前記各領域における吸収成分の濃度を算出して前記対象物における吸収成分濃度分布を求める吸収成分濃度算出手段を更に具備してもよい。
【００２３】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての吸収係数に基づいて、前記の各領域における吸収係数の変化量（差分）から各領域における吸収係数の絶対値が求められる。このように本発明の方法および装置によれば、均一な吸収係数を有しかつ測定対象物と同一の外形を持つファントムから得た基準値を用いることなく各領域における吸収係数の絶対値が求められることとなる。そして、各領域における吸収係数の絶対値が求まれば、既知の吸収成分のモル吸収係数等を用いて各領域における吸収成分の濃度が求められる。本発明で求められる吸収係数の変化量分布の誤差は従来の方法による場合と比較して極めて少なくなるため、それに基づいて求められる吸収係数の絶対値分布並びに吸収成分の濃度分布の精度が高くなる。
【００２４】
上記本発明の方法を少なくとも２つの吸収成分を含有している対象物に適用する場合には、前記光入射ステップにおいて前記対象物中に入射される測定光が、該吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有していることが好ましい。この場合、前記光検出ステップにおいて前記少なくとも２つの波長を有する測定光をそれぞれ検出し、前記測定値を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記測定値を求め、前記基準値を得るステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記平均値を算出し、前記内部特性変化量分布を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の変化量を算出し、前記吸収係数絶対値分布を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の絶対値を算出し、前記吸収成分濃度分布を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収成分の濃度を算出することによって、前記対象物における前記各吸収成分の濃度分布を高精度に求めることが可能となる。
【００２５】
また、上記本発明の装置を用いて少なくとも２つの吸収成分を含有している対象物を測定する場合には、前記光入射手段において前記対象物中に入射される測定光が、該吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有していることが好ましい。この場合、前記光検出手段において前記少なくとも２つの波長を有する測定光をそれぞれ検出し、前記測定値取得手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記測定値を求め、前記基準値算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記平均値を算出し、前記内部特性変化量算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の変化量を算出し、前記吸収係数絶対値算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の絶対値を算出し、前記吸収成分濃度算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収成分の濃度を算出することによって、前記対象物における前記各吸収成分の濃度分布を高精度に求めることが可能となる。
【００２６】
次に、本発明の方法および装置によって測定されるべき内部特性が等価散乱係数である場合について説明する。
【００２７】
この場合、本発明の方法は、前記対象物の平均の吸収係数および平均の等価散乱係数を求める（好ましくは前記基準値に基づいて求める）ステップと、前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数（等価散乱係数用の広がり関数）を選択するステップとを更に具備していることが好ましく、前記内部特性変化量分布を求めるステップにおいて前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における等価散乱係数の変化量を算出することができる。
【００２８】
また、本発明の装置は、前記対象物の平均の吸収係数および平均の等価散乱係数を求める（好ましくは前記基準値に基づいて求める）平均吸収および散乱係数検出手段と、前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数（等価散乱係数用の広がり関数）を選択する広がり関数選択手段とを更に具備していることが好ましく、前記内部特性変化量算出手段において前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における等価散乱係数の変化量を算出することができる。
【００２９】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての吸収係数および等価散乱係数に対応して選択された広がり関数に基づいて、前記の各領域における等価散乱係数の変化量（差分）が求められる。したがって、このような本発明の方法および装置によれば、吸収係数および／または等価散乱係数をゼロとして計算する場合に比べて、それらによる実効光路長の変化に基づく誤差の発生が充分に防止され、測定精度が向上する。
【００３０】
さらに、上記本発明の方法は、前記等価散乱係数の変化量および前記平均の等価散乱係数を用いて、前記各領域における等価散乱係数の絶対値を算出して前記対象物における等価散乱係数絶対値分布を求めるステップを更に具備してもよい。
【００３１】
また、上記本発明の装置は、前記等価散乱係数の変化量および前記平均の等価散乱係数を用いて、前記各領域における等価散乱係数の絶対値を算出して前記対象物における等価散乱係数絶対値分布を求める等価散乱係数絶対値算出手段を更に具備してもよい。
【００３２】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての等価散乱係数に基づいて、前記の各領域における等価散乱係数の変化量（差分）から各領域における等価散乱係数の絶対値が求められる。このように本発明の方法および装置によれば、均一な等価散乱係数を有しかつ測定対象物と同一の外形を持つファントムから得た基準値を用いることなく各領域における等価散乱係数の絶対値が求められることとなる。本発明で求められる等価散乱係数の変化量分布の誤差は従来の方法による場合と比較して極めて少なくなるため、それに基づいて求められる等価散乱係数の絶対値分布の精度が高くなる。
【００３３】
次に、本発明の方法および装置によって測定されるべき内部特性が屈折率である場合について説明する。
【００３４】
この場合、本発明の方法は、前記対象物の平均の吸収係数、平均の等価散乱係数および平均の屈折率を求める（好ましくは前記基準値に基づいて求める）ステップと、前記平均の吸収係数、平均の等価散乱係数および平均の屈折率に対応する広がり関数（屈折率用の広がり関数）を選択するステップとを更に具備していることが好ましく、前記内部特性変化量分布を求めるステップにおいて前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における屈折率の変化量を算出することができる。
【００３５】
また、本発明の装置は、前記対象物の平均の吸収係数、平均の等価散乱係数および平均の屈折率を求める（好ましくは前記基準値に基づいて求める）平均吸収および散乱係数検出手段と、前記平均の吸収係数、平均の等価散乱係数および平均の屈折率に対応する広がり関数（屈折率用の広がり関数）を選択する広がり関数選択手段とを更に具備していることが好ましく、前記内部特性変化量算出手段において前記複数の測定値と前記基準値と前記広がり関数とを用いて前記各領域における屈折率の変化量を算出することができる。
【００３６】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての吸収係数、等価散乱係数および屈折率に対応して選択された広がり関数に基づいて、前記の各領域における屈折率の変化量（差分）が求められる。したがって、このような本発明の方法および装置によれば、吸収係数および／または等価散乱係数をゼロとして計算する場合に比べて、それらによる実効光路長の変化に基づく誤差の発生が充分に防止され、測定精度が向上する。
【００３７】
さらに、上記本発明の方法は、前記屈折率の変化量および前記平均の屈折率を用いて、前記各領域における屈折率の絶対値を算出して前記対象物における屈折率絶対値分布を求めるステップを更に具備してもよい。
【００３８】
また、上記本発明の装置は、前記屈折率の変化量および前記平均の屈折率を用いて、前記各領域における屈折率の絶対値を算出して前記対象物における屈折率絶対値分布を求める屈折率絶対値算出手段を更に具備してもよい。
【００３９】
このような本発明の方法および装置によれば、内部が不均一な測定対象物に対して測定された平均値としての屈折率に基づいて、前記の各領域における屈折率の変化量（差分）から各領域における屈折率の絶対値が求められる。このように本発明の方法および装置によれば、均一な屈折率を有しかつ測定対象物と同一の外形を持つファントムから得た基準値を用いることなく各領域における屈折率の絶対値が求められることとなる。本発明で求められる屈折率の変化量分布の誤差は従来の方法による場合と比較して極めて少なくなるため、それに基づいて求められる屈折率の絶対値分布の精度が高くなる。
【００４０】
上述の本発明の方法は、前記の求められた分布に基づいて、前記対象物内部における分布を示す画像を表示するステップを更に具備してもよい。また、上述の本発明の装置は、前記の求められた分布に基づいて、前記対象物内部における分布を示す画像を表示する画像表示手段を更に具備してもよい。このような本発明の方法および装置によれば、高精度に求められた内部特性分布を画像化して表示することが可能である。
【００４１】
【発明の実施の形態】
以下、図面を参照しつつ本発明の好適な実施形態について詳細に説明する。尚、図面中、同一又は相当部分には同一符号を付することとする。
【００４２】
先ず、図２および図３を参照して本実施形態で使用した光ＣＴの画像化原理を説明する。なお、散乱されつつ進行する光に関しては３次元座標を用いて考える必要があるが、以下では説明を簡単にするために２次元座標を用いて説明する。
【００４３】
先ず、散乱吸収体内部をＮ個のボクセルに区切り、吸収係数が存在する条件下における散乱吸収体に対する入射光量と出射光量（検出光量）との関係を考える。等価散乱係数μ'sおよび吸収係数μaを均一とした散乱吸収体内部の模式図を図２に示す（Ｎ＝２５）。このとき入射光量をＩ0、検出光量をＩd0、散乱吸収体内部の等価散乱係数μ'sおよび吸収係数μaを均一とした場合の各ボクセル内での実効光路長をＷj、散乱および反射等によって入射光が散乱吸収体の外に出る割合を示す減衰定数をＤsrとすると、下記式（１）が成り立つ。
【００４４】
【数１】
次に、幾つかのボクセル内に等価散乱係数は同じであるが吸収係数が相違する媒質を挿入した以外は図２に示すものと同様の散乱吸収体内部の模式図を図３に示す。但し、図３に示す散乱吸収体に用いた各媒質の吸収係数μai（i=1,2,・・・,N）と図２に示す散乱吸収体に用いた媒質の吸収係数μaとの関係は以下の式（２）に示すような関係にある。
【００４６】
【数２】
このときの入射光量をＩ0、検出光量をＩd1とし、散乱および反射等によって入射光が散乱吸収体の外に出る割合を示す減衰定数Ｄsrは吸収係数を均一とした場合（図２）と変わらないと考えると、検出光量をＩd1は下記式（３）のように表わされる。
【００４８】
【数３】
したがって、式（３）より以下の式（４）が求められる。
【００５０】
【数４】
このように、基準の光量Ｉd0を用いれば、求めたい吸収係数μaと実際の実験系で測定可能な検出光量Ｉd1との関係より、実効光路長Ｗjが決まれば散乱吸収体内部の吸収係数μaの分布が求められることとなる。式（４）は１組の光入射位置−光検出位置で成り立つ関係を示す。したがって、例えばＮ個の吸収係数（未知数）を求める場合は、光入射位置−光検出位置の組み合わせをＮ組み揃え、それぞれの組み合わせについて成り立つＮ個の式（４）の連立方程式を解くことによってＮ個の吸収係数を求めることができる。
【００５２】
すなわち、光入射位置−光検出位置のＮ個の組み合わせについて成り立つＮ個の式（４）の連立方程式を行列式の形で表わすと、下記式（５）が得られる。
【００５３】
【数５】
但し、ΔＩはｌｎＩd0−ｌｎＩd1、Ｗは各ボクセルの実効光路長の分布を示す広がり関数を表わす。光入射位置Ｍ（Ｍ1〜ＭX）の数をＸ、光検出位置ｍ（ｍ1〜ｍx）の数をｘ、光入射位置Ｍ及び光検出位置ｍのときの光量の変化量をΔＩMm、光入射位置Ｍ及び光検出位置ｍのときの各ボクセルの広がり関数をＷMmとすると、［ΔＩMm］は（Ｘ×ｘ）×１、［ＬMm］は（Ｘ×ｘ）×Ｎ、［Δμan］はＮ×１のマトリックスとなる。したがって、吸収係数の変化量Δμanは下記式（６）の連立方程式を解くことにより求められることとなる。なお、式（６）の連立方程式は、Ｘ×ｘ＝ＮとなるようにＸおよびｘの値を選択して解くことが好ましい。
【００５５】
【数６】
かかる画像再構成方法で散乱吸収体内部の吸収係数の定量を行う場合、基本的に図２に示すような基準となる状態が必要となり、上記の形態では吸収係数が均一の状態を基準として想定していたため式（２）と式（４）から各ボクセルの吸収係数を求めていた。しかしながら、かかる画像化方法においては内部の吸収係数の値とその時の各光検出位置での光量さえ予め分かっていればよく、実際に用いるべき基準の状態は特に制限されない。すなわち、例えばある条件下での内部の吸収係数の値とその時の各光検出位置での光量を基準として内部の吸収係数を求めると、その吸収係数の値は基準の値との差という形で求められる。
【００５７】
従来はこのように基準となる内部の吸収係数の値とその時の各光検出位置での光量を、測定対象である散乱吸収体に対して内部の吸収係数のみが相違するファントムモデルあるいはシミュレーションモデルから得ていた。
【００５８】
しかしながら、本発明においては、測定対象物の表面にありかつその対象物内の一点（例えば対象物の中心）に対して相対的に同じ位置関係にある光入射位置と光検出位置との複数の組み合わせによって求められた複数の測定値の平均値を、内部特性分布を求めるための基準値として用いる。また、各ボクセルの実効光路長の作成方法は、本発明者らによる特願平８−６６１９「光ＣＴ装置及び光ＣＴによる画像再構成方法」等に記載されている。本実施形態においてはこの作成方法にしたがって、ある光入射位置−光検出位置の関係における各ボクセルの実効光路長の分布（すなわち広がり関数）を、吸収係数の平均値、等価散乱係数の平均値等に基づいて予め用意する。このように、本実施形態においては、従来は必要であった物理モデルやシミュレーションモデルから基準値を求めることなく、測定対象物についての測定値から直接的に基準値を求めることが可能となり、その基準値に基づいて測定対象物の内部特性の分布を高精度に測定することが可能となる。
【００５９】
以上、ＣＷ計測を用いる実施形態について説明したが、本発明において時間分解計測を適用することも可能である。すなわち、前記の式（３）においては検出光を時間０〜ｔ（ｓ）間に検出器で受光した検出光量の積分値として表わしているが、かかる関係式は光源にパルス光を用いた場合に検出器で得られる時間分解波形についても成り立つ。式（３）をある時間間隔ｔ1〜ｔ2について書き直すと以下の式（７）が得られる。
【００６０】
【数７】
但し、［Ｉd0］t1-t2、［Ｉd1］t1-t2は各検出光の時間ｔ1〜ｔ2間の時間分解波形の光量、［Ｗj］t1-t2は時間ｔ1〜ｔ2間の広がり関数を示す（ｊは各ボクセルの番号を示す）。また、０≦ｔ1＜ｔ2である。したがって、式（７）より以下の式（８）が求められる。
【００６２】
【数８】
このように、時間分解計測を用いる方法においては、測定時間を様々に区切って方程式の数を増やすことによって、ボクセル数と同じＮ個の方程式を解くことによって各ボクセルの吸収係数が求められる。
【００６４】
また、上記の実施形態においては測定されるべき内部特性が吸収係数の場合について説明したが、本発明は等価散乱係数、屈折率といった他の内部特性に適用することも可能である。すなわち、対象物の内部を透過して受光された光は、対象物の内部がもつ吸収係数、等価散乱係数のみならず、かかる対象物がもつ全ての内部特性の影響を受けており、さらにそれらはそれぞれ線形的かつ独立的に受光された光に作用を及ぼしている。なお、互いに影響を及ぼし合う内部特性を求める場合は、それらを１つの内部特性とみなすことにより独立性が保たれる。
【００６５】
これらの関係により、対象物がもつ全ての内部特性の値は、各ボクセル内における内部特性の受光される光に対する寄与度を示す関数（広がり関数）と、受光された光とを用いた式で表わされ、これらを用いれば例えば“Forward and Inverse Calculations for 3-D Frequency-Domain Diffuse Optical Tomography”（Brian W.Pogue et al., SPIE vol.2389, p.328-p.338）に記載されているように等価散乱係数、吸収係数の分布を求めることができる。このような場合においても、検出光の振幅や位相といったパラメータの基準値の導出法として本発明の方法は適用できる。したがって、吸収係数、等価散乱係数、屈折率といった内部特性に対しても、受光される光に対する寄与を示す関数（広がり関数）と受光された光とを用いた関係式に本発明にかかる基準値導出法を適用することによってこれらの内部特性の定量化が可能である。このような関係式としては、以下の諸式が挙げられる。すなわち、吸収係数および等価散乱係数を求める場合に適用可能な式として、前記式（４）および式（８）の変形式としてそれぞれ下記の式（４’）および式（８’）が挙げられる。
【００６６】
【数９】
【数１０】
また、吸収係数、等価散乱係数および屈折率を求める場合に適用可能な式として、前記式（４）および式（８）の変形式としてそれぞれ下記の式（４''）および式（８''）が挙げられる。
【００６９】
【数１１】
【数１２】
【００７１】
次に、本発明の内部特性分布計測装置について説明する。図４に本発明の装置の一実施形態の模式図を示す。
【００７２】
図４に示す装置は、１２個の光ファイバーホルダー１〜１２（以下、場合により「光ファイバーホルダー群」と総称する）を備えており、各光ファイバーホルダー１〜１２は散乱吸収体ＳＭ(scattering medium) の一断面の周りに等間隔に配置されており（図４に示す装置においては、散乱吸収体ＳＭの中心から３０度間隔で放射線状に伸びる各線上に配置されている）、時計回りに１〜１２の番号が割り振られている。
【００７３】
光ファイバーホルダー１〜１２はそれぞれ光入射ファイバー１ａ〜１２ａと光検出ファイバー１ｂ〜１２ｂとを備える。光入射ファイバー１ａ〜１２ａと光検出ファイバー１ｂ〜１２ｂとは図４に示すように並列に束ねられている構造でもかまわないが、図５（ａ）に示すように光入射ファイバー１ａの回りを複数の光検出ファイバー１ｂ（バンドルファイバー）が包むように束ねられている構造、あるいは図５（ｂ）に示すように光入射ファイバー１ａと光検出ファイバー１ｂとを光ファイバーホルダー内で光カプラ１ｃでカップリングさせている構造であってもよい。図５（ａ）や図５（ｂ）に示すような構造を採用すると、散乱吸収体ＳＭの周囲にあたるファイバー端面は１つになり、両ファイバーを上下２段あるいは左右２列に並べる場合よりも光入射ファイバー端と光検出ファイバー端との位置ずれをなくすことができ、誤差を軽減できる傾向がある。
【００７４】
光入射ファイバー１ａ〜１２ａには波長選択器２０を介して光源３０が光学的に接続されている。そして、光源３０から発せられた光は波長選択器２０で波長選択され、光ファイバーホルダー１〜１２を通して計測対象である散乱吸収体ＳＭの表面に入射される。光源３０には、発光ダイオード、レーザーダイオード、Ｈｅ−Ｎｅレーザー等種々のものが使用できる。また、光源３０は、パルス光や方形波光、またはそれらの変調光を発生するものでもよい。
【００７５】
本実施形態において使用する光源３０は単波長の光（測定光）を発生するものであってもよいが、２波長以上の光（測定光）を発生可能なものであることが好ましい。計測に使用する光の波長は、計測対象に応じて適宜選択される。一般に生体では、ヘモグロビン等の吸収特性から７００ｎｍ以上の光、特に好ましくは可視光または近赤外光、を使用することが好ましい。例えば酸化ヘモグロビンと還元ヘモグロビンを対象とする場合は、図６に示すようにそれらの吸収係数が異なるから、適宜に波長を選択して使用することによってこれらを分離して計測することが可能となる。
【００７６】
光検出ファイバー１ｂ〜１２ｂには検出器４０が光学的に接続されている。そして、散乱吸収体ＳＭ中で散乱されつつ透過した光（測定光）は光ファイバーホルダー１〜１２の光検出ファイバー１ｂ〜１２ｂを介して検出器４０に導かれ、光検出器４０で受光信号を検出信号（電気信号）に変換して増幅し、それぞれに対応する検出信号が出力される。光検出器４０は光電子増倍管のほか、光電管、フォトダイオード、アバランシェフォトダイオード、ＰＩＮフォトダイオード等、あらゆる種類の光検出器を使用することができる。光検出器４０の選択に際しては、使用される測定光の波長の光が検出できる分光感度特性をもっていれば良い。また、光信号が微弱であるときは高感度あるいは高利得の光検出器を使用することが好ましい。なお、光検出ファイバー１ｂ〜１２ｂおよび光検出器４０の受光面以外の場所は、光を吸収あるいは遮光する構造にしておくことが望ましい。また、散乱吸収体ＳＭの内部を拡散伝搬した光が複数の波長の光を含む場合には、光検出器４０と散乱吸収体ＳＭとの間に波長選択フィルタ（図示せず）を適宜配置してもよい。
【００７７】
光源３０および検出器４０には制御部５０が接続されており、入射あるいは受光の際に使用される光ファイバーホルダー１〜１２の選択は制御部５０で行なわれる。すなわち、制御部５０によって光入射ファイバーから順次（例えば１ａ→２ａ→３ａ→・・・→１２ａ）一定時間間隔で散乱吸収体ＳＭ中に測定光が入射されるように制御され、それに同期して測定光が入射された光入射ファイバーに対して所定の位置関係にある光検出ファイバーから測定光が検出されるように制御される。本実施形態においては、測定光が入射された光入射ファイバーと異なる位置にある全ての光検出ファイバー（例えば光入射ファイバーが１ａの場合は光検出ファイバー２ｂ〜１２ｂ）から測定光がそれぞれ検出されるが、特にかかる組み合わせに限定されるものではない。
【００７８】
また、複数の波長を有する測定光を使用する場合は、入射される測定光の波長も制御部５０によって制御される。具体的な手法としては、異なる波長の光を時分割で入射させて使用する手法と、後述する異なる波長の光を同時に含む光を使用する手法とがある。具体的な波長選択手段としては、ミラーを用いた光ビーム切り換え器、フィルターを用いた波長切り換え器、光スイッチを用いた光切り換え器等がある（図７（ａ））。
【００７９】
上記の光入射ファイバー１ａ〜１２ａ、波長選択器２０、光源３０および制御部５０が本発明にかかる光入射手段を構成し、上記の光検出ファイバー１ｂ〜１２ｂ、検出器４０および制御部５０が本発明にかかる光検出手段を構成する。
【００８０】
制御部５０には演算処理部（例えばＣＰＵ）６０が電気的に接続されており、演算処理部６０には記憶部（例えばハードディスク、フレキシブルディスク）７０および表示部（例えばディスプレイ、プリンタ）８０が電気的に接続されている。そして、検出器４０から出力された検出信号は制御部５０を介して演算処理部６０へ導かれる。
【００８１】
上記の演算処理部６０および記憶部７０が、本発明にかかる測定値取得手段と、基準値算出手段と、内部特性変化量算出手段と、平均吸収および散乱係数検出手段と、広がり関数選択手段と、吸収係数絶対値算出手段と、吸収成分濃度算出手段と、等価散乱係数絶対値算出手段と、屈折率絶対値算出手段とを構成し、上記の表示部８０が画像表示手段を構成する。これらの本発明にかかる諸手段については、図８に示す本発明の方法の一実施形態のフローチャートに基づいて以下に詳細に説明する。
【００８２】
（１）図８に示す方法においては、先ず、以下に説明するようにして光ＣＴによる測定データ（Ｉd1｛Ｍ，ｍ｝））が取得される（Ｓ１００）。なお、Ｍは光入射ファイバーの番号を示し、ｍは光検出ファイバーの番号を示す。
【００８３】
すなわち、光入射ファイバー１ａ〜１２ａから順次散乱吸収体ＳＭ中に測定光を入射し、散乱吸収体ＳＭ中で散乱されつつ透過した各測定光をその測定光が入射された光入射ファイバーと異なる位置にある全ての光検出ファイバー（例えば光入射ファイバーが１ａの場合は光検出ファイバー２ｂ〜１２ｂ）から順次あるいは同時に検出する。なお、各測定光を同時に検出する場合は、光検出器４０が光検出ファイバー数分必要である。
【００８４】
そして、各光検出ファイバーで検出された各測定光に基づいて光検出器４０から検出信号が発せられる。そして、これらの検出信号はそれぞれ演算処理部６０において処理され、検出された測定光のそれぞれの検出光量に比例した測定値に変換されて、得られた測定値は記憶部７０に一次的に保存される。具体的には、演算処理部６０において、光源３０から発せられた光発生に同期した信号を利用して、検出信号に対する時間域での積分演算を行い、検出光量に比例した測定値を求める。ただし、パルス光等を利用する場合には、この同期信号を省略することができる。この種の演算処理は、演算処理手段に組み込んだマイクロコンピュータ等で高速に実行することができる。また、演算処理部６０において、測定値を平均化フィルタリングや最小二乗フィッティング等を利用して修正してもよい。
【００８５】
２）次に、演算処理部６０において位置関係が相対的に同じである光入射ファイバーと光検出ファイバーとの複数の組み合わせによって求められた複数の測定値が抽出され、それらの測定値の平均値である基準値（Ｉd0｛Ｍ，ｍ｝））が算出される（Ｓ１１０）。
【００８６】
すなわち、演算処理部６０では、吸収係数変化量等を求めるための基準値を得るために、光入射−光検出の位置関係が相対的に同じである光入射−光検出位置対毎に測定値の平均値を求める。
【００８７】
図４に基づいて具体的に説明すると、例えば光入射ファイバーと散乱吸収体ＳＭの中心と光検出ファイバーとのなす角度が１８０度の位置関係にある場合、相対的に同じ光入射−光検出位置の組み合わせは、（光入射ホルダー番号，光検出ホルダー番号）で表わすと(1,7)、(2,8)、(3,9)、(4,10)、(5,11)、(6,12)となる。なお、光の相反定理が成り立たない状態である場合は光入射−光検出位置の反対の組み合わせも考慮に入れる必要性がある。
【００８８】
各々の測定値をI(1,7)、I(2,8)、I(3,9)、I(4,10) 、I(5,11) 、I(6,12) とした場合、これらの平均値は下式となる。
【００８９】
I(ave_180)
＝｛I(1,7)＋I(2,8)＋I(3,9)＋I(4,10) ＋I(5,11) ＋I(6,12) ｝／６
このI(ave_180)を光入射−光検出の位置関係が１８０度の時の基準値とする。
【００９０】
同様にして、I(ave_150)、I(ave_120)、I(ave_90) 、I(ave_60) 、I(ave_30) も求め、これらの平均値を上記各位置関係における基準値として記憶部７０に一次的に保存する。
【００９１】
３）次に、本実施形態においては、上記各位置関係における基準値等に基づいて、平均吸収係数μa0および平均等価散乱係数μ'S0を光拡散理論等を利用して求めることができる（Ｓ１２０）。
【００９２】
すなわち、各角度毎の基準値から内部の吸収係数と等価散乱係数とを求め、さらにそれらの平均値を算出し、それらを散乱吸収体ＳＭ内部の平均吸収係数μa0と平均等価散乱係数μ'S0として記憶部７０に一次的に保存する。また、いずれか１つの角度の基準値、例えばI(ave_180)の値のみから散乱吸収体ＳＭ内部の平均吸収係数μa0と平均等価散乱係数μ'S0とを求めてもかまわない。
【００９３】
なお、上記基準値から散乱吸収体ＳＭ内部の吸収係数と等価散乱係数とを求める方法としては、例えば“Imaging diffusive media using time-independent and time-harmonic sources；dependence of image quality on imaging algorithms, target volume weight matrix, and view angles ”(Jenghwa Chang et.al., SPIE vol.2389）に記載されている方法がある。
【００９４】
４）次に、本実施形態においては、上記平均吸収係数μa0と平均等価散乱係数μ'S0に対応する広がり関数（Ｗθ）を選択する（Ｓ１３０）。
【００９５】
すなわち、予め用意されて記憶部７０に記憶してあった広がり関数の中から、上記で求めた平均吸収係数μa0および平均等価散乱係数μ'S0と合う広がり関数を選択する。この場合、実際の測定値から求めた吸収係数および等価散乱係数に基づいて広がり関数を選択するため、適当に仮定したそれらの値を用いるよりも誤差要因は除かれる。
【００９６】
なお、このような「広がり関数」とは、各ボクセルにおける光（測定光）の広がり方を示す関数を意味し、各ボクセルにおける実効光路長に関するいわゆる重み関数や各ボクセルにおける測定光に対する寄与度に関するいわゆる寄与関数を包含する概念である。本発明にかかる広がり関数は上記の重み関数あるいは寄与関数のいずれでもよい。このような広がり関数に関しては、例えば“A Perturbation Model for Imaging in Dense Scattering Media: Derivation and Evaluation of Imaging Operation”（H.L.Graber et.al., SPIE vol. IS11）、“Initial assessment of a simple system for frequency domein diffuse optical tomography”（B.W.Pogue et al.,Phys.Med.Biol.40(1995)p.1709-p.1729）、本発明者らによる特願平８−６６１９「光ＣＴ装置及び光ＣＴによる画像再構成方法」に記載されている。本実施形態においては特願平８−６６１９に記述された作成方法にしたがって、以下に示すような時間項の入っていない光拡散方程式を用いて広がり関数を予め用意する。
【００９７】
Δφ−μaＤ-1φ＝０
ここで、Ｄ＝1/{3(1-g)μs}＝1/3μ's、φ：光子密度、Ｄ：光拡散定数、μa：吸収係数、μ'S：等価散乱係数、ｇ：散乱吸収体による光子の散乱角の余弦の平均値である。また、屈折率分布を求める際に好適に用いられる時間項の入っている光拡散方程式を以下に示す。
【００９８】
【数１３】
ここで、Ｄ(r)＝1/{3(1-g)μs(r)}＝1/3μ's(r)、φ(r,t)：位置ｒ、時刻ｔにおける光子密度、Ｃ：媒質中の光速、Ｄ：光拡散定数、μa：吸収係数、Ｓ(r,t)：光源、μ'S：等価散乱係数、ｔ：時間、ｒ：位置、ｇ：散乱吸収体による光子の散乱角の余弦の平均値である。なお、真空中の光速をＣ’、被測定物の屈折率をｎとするとＣ＝Ｃ’ｎと表わすことができる。
【０１００】
また、平均吸収係数μa0及び平均等価散乱係数μ'S0と合う広がり関数とは、具体的には、測定対象物の平均吸収係数及び平均等価散乱係数と同じ平均吸収係数及び平均等価散乱係数を有しかつ測定対象物と同じ形状を有する対象物に対して光入射−光検出位置の相対的関係を実測の場合と同様にした場合に求められるであろう光の広がり方を示す関数のことであり、平均吸収係数μa0及び平均等価散乱係数μ'S0等に基づいて選択される。
【０１０１】
また、記憶部７０に対象物を複数のブロック（ボクセル）に分けた際に生じる歪みを補正するための補正項を記憶しておいてもよく、その場合は演算処理部６０で前記測定値および／または前記基準値を補正することが可能である（Ｓ１４０）。このようなボクセルに関する補正とは、例えば光入射位置ａと光検出位置ｂ，ｃとの間の距離が同じにも拘らずボクセルの切り方によってボクセル上でのトータルの距離が異なるような場合に、その相違分を測定値及び／又は基準値の補正項として利用することである。
【０１０２】
５）続いて、演算処理部６０において前記複数の組み合わせによって求められた複数の測定値と、前記基準値と、前記広がり関数とを用いて、複数の領域に分割された前記各領域における吸収係数の変化量Δμa を算出し（Ｓ１５０）、出力する（Ｓ１６０）。
【０１０３】
すなわち、前記各角度毎の基準値と前記測定値と前記広がり関数とを用いて吸収係数の変化量を求める。その際に成り立つ関係式は、前記式（４）に対応させると、例えば光入射−光検出の位置関係が１８０度の場合は基準値Ｉd0はI(ave_180)となり、測定値Ｉd1はI(1,7)、I(2,8)、I(3,9)、I(4,10) 、I(5,11) 、I(6,12) となる。この時、基準値Ｉd0の吸収係数は散乱吸収体ＳＭ内部の平均の吸収係数である。さらに、光入射−光検出の位置関係が１８０度の時の広がり関数をＷθとすると下記の式（4-1）〜式（4-6）が成り立ち、これらの連立方程式（前記式（８）の連立方程式）を各位置関係毎に立てて（つまり未知数分と同じ数の連立方程式を立てて）解くことによって各領域における吸収係数の変化量Δμa が算出される。
【０１０４】
【数１４】
【数１５】
【数１６】
【数１７】
【数１８】
【数１９】
なお、散乱吸収体ＳＭ内部の吸収係数の変化量の空間的分布や濃度変化の空間的分布を求める場合、このように成り立った関係式を散乱吸収体ＳＭの内部を分割したボクセル（ボックスセル）の個数分と同じ数の連立方程式を立てて解けばよい。本実施形態においては、共役勾配法を用いた。また、式数がボクセルの個数分より少ない場合や多い場合であっても、特異値分解法等を用いれば特異な問題を非特異な問題に変えられるため、内部特性の分布が求められる。
【０１１１】
このようにして求められた各領域における吸収係数の変化量Δμa に基づいて測定対象物内部における吸収係数変化量に関する吸収分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ１７０）。
【０１１２】
なお、上記に説明したような演算処理部６０におけるΔμa の算出方法から吸収分布を求めて画像を表示する方法としては他にも様々な方法が知られている。このような方法は、例えば“Optical Back Projection Tomography in Heterogeneous Diffusive Media ”（S.B.Cloak et al., in Advances in Optical Imaging and Photon Migration, 1996 Technical Digest; Optical Society of America, Washington DC, 1996, pp.147-149 ）、“Back-projection image reconstruction using photon density wave in tissues ”（S.A.Walker et al., SPIE vol.2389, pp.350, 1995）、“in Photon Propagation in Tissues”（S.Cloak et al., SPIE vol.2626, 1995 ）、“Optical tomography by the temporally extrapolated absorbance method ”（Ichiro Oda et al., APPLIED OPTICS, vol.35, No.01, 1996）等に記載されている。これらの方法はＸ線ＣＴで知られているバックプロジェクション（Back Projection ）法もしくはそれらを変形した方法であり、上記アルゴリズムの変わりに画像再構成をさせる方法である。
【０１１３】
また、吸収成分の既知のモル吸収係数を用いて上記各領域における吸収係数の変化量Δμa から各領域における吸収成分の濃度差の絶対値を算出することが可能であり（Ｓ１８０）、このようにして求められた各領域における吸収成分の濃度差の絶対値に基づいて測定対象物内部における吸収成分濃度差に関する分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ１９０）。
【０１１４】
更に、上記の各領域における吸収係数の変化量Δμa および前記平均吸収係数μa0を用いて上記各領域における吸収係数の絶対値μa を算出することが可能であり（Ｓ２００）、このようにして求められた各領域における吸収係数の絶対値μa に基づいて測定対象物内部における吸収係数絶対値に関する分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ２１０）。
【０１１５】
更にまた、吸収成分の既知のモル吸収係数を用いて上記の各領域における吸収係数の絶対値μa から各領域における吸収成分の濃度を算出することが可能であり（Ｓ２２０）、このようにして求められた各領域における吸収成分の濃度に基づいて測定対象物内部における吸収成分濃度に関する分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ２３０）。
【０１１６】
また、散乱吸収体ＳＭが少なくとも２つの吸収成分、例えば酸化および還元ヘモグロビン、を含有する場合は、それらの吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有する測定光を使用し、各波長を有する測定光のそれぞれに関して前記測定値、前記基準値を求め、それらに基づいて各波長を有する測定光のそれぞれに関して前記吸収係数変化量、前記吸収係数絶対値を求めることによって、各吸収成分の濃度分布が求められる。
【０１１７】
以下に、上記２波長分光法を利用したヘモグロビンの濃度の計測について説明する。
【０１１８】
哺乳類の脳における吸収成分の主なものは、水、チトクローム（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)を使用すれば上記の式はさらに簡単になる。
【０１２５】
以上、本発明の好適な実施形態について説明したが、本発明は勿論上記実施形態に限定されるものではない。
【０１２６】
すなわち、上記実施形態においては内部特性として吸収係数を求めているが、前述のように本発明は等価散乱係数の計測にも適用可能であり、吸収係数および等価散乱係数を求める一実施形態のフローチャートを図９に示す。
【０１２７】
図９に示す方法においては、吸収係数の測定に関しては図８に示す方法と同様であるが、平均吸収係数μa0と平均等価散乱係数μ'S0に対応する広がり関数（Ｗθ）を選択する際に、吸収係数用の広がり関数（Ｗμa,j）と等価散乱係数用の広がり関数（Ｗμ'S,j）とを選択しておくことが好ましい（Ｓ１３０）。
【０１２８】
そして、演算処理部６０において前記複数の組み合わせによって求められた複数の測定値と、前記基準値と、前記広がり関数とを用いて、複数の領域に分割された前記各領域における等価散乱係数の変化量Δμ'Sを算出し（Ｓ２４０）、出力する（Ｓ２５０）。すなわち、前述の式（４’）に基づいて、前記各角度毎の基準値と前記測定値と前記広がり関数と吸収係数の変化量とを用いて等価散乱係数の変化量を求める。より具体的には、前述の式（4-1）〜式（4-6）と同様に式（４’）に基づいて連立方程式が成り立ち、これらの連立方程式を各位置関係毎に立てて解くことによって各領域における等価散乱係数の変化量Δμ'Sが算出される。
【０１２９】
このようにして求められた各領域における等価散乱係数の変化量Δμ'Sに基づいて測定対象物内部における等価散乱係数変化量分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ２６０）。
【０１３０】
更に、上記の各領域における等価散乱係数の変化量Δμ'Sおよび前記平均等価散乱係数μ'S0を用いて上記各領域における等価散乱係数の絶対値μ'Sを算出することが可能であり（Ｓ２７０）、このようにして求められた各領域における等価散乱係数の絶対値μ'Sに基づいて測定対象物内部における等価散乱係数絶対値に関する分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ２８０）。
【０１３１】
また、本発明は屈折率の計測にも適用可能であり、吸収係数、等価散乱係数および屈折率を求める一実施形態のフローチャートを図１０に示す。
【０１３２】
図１０に示す方法においては、吸収係数および等価散乱係数の測定に関しては図８および図９に示す方法と同様であるが、各位置関係における基準値等に基づいて、平均吸収係数μa0および平均等価散乱係数μ'S0を求める際に、平均屈折率ｎ0も併せて求める（Ｓ１２０）。但し、平均屈折率ｎ0としては水の屈折率（1.33）を用いてもよい。また、広がり関数（Ｗθ）を選択する際に、吸収係数用の広がり関数（Ｗμa,j）と等価散乱係数用の広がり関数（Ｗμ'S,j）と屈折率用の広がり関数（Ｗn,j）とを選択しておくことが好ましい（Ｓ１３０）。
【０１３３】
そして、演算処理部６０において前記複数の組み合わせによって求められた複数の測定値と、前記基準値と、前記広がり関数とを用いて、複数の領域に分割された前記各領域における屈折率の変化量Δｎを算出し（Ｓ２９０）、出力する（Ｓ３００）。すなわち、前述の式（４''）に基づいて、前記各角度毎の基準値と前記測定値と前記広がり関数と吸収係数の変化量と等価散乱係数の変化量とを用いて屈折率の変化量を求める。より具体的には、前述の式（4-1）〜式（4-6）と同様に式（４''）に基づいて連立方程式が成り立ち、これらの連立方程式を各位置関係毎に立てて解くことによって各領域における屈折率の変化量Δｎが算出される。
【０１３４】
このようにして求められた各領域における屈折率の変化量Δｎに基づいて測定対象物内部における屈折率変化量分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ３１０）。
【０１３５】
更に、上記の各領域における屈折率の変化量Δｎおよび前記平均屈折率ｎ0を用いて上記各領域における屈折率の絶対値ｎを算出することが可能であり（Ｓ３２０）、このようにして求められた各領域における屈折率の絶対値ｎに基づいて測定対象物内部における屈折率絶対値に関する分布が求められ、その分布を示す画像が表示部８０に表示される（Ｓ３３０）。
【０１３６】
そして、このように屈折率に関する分布が求められると、血糖値の分布を求めることが可能となる。屈折率の変化により血糖値を検出する方法としては、例えば“Possible correlation between blood glucose concentration and the reduced scattering coefficient of tissues in the near infrared”（John S.Maier et al., OPTICS LETTERS vol.19, No.24, December 15, 1994）に記載の方法がある。生体組織のグルコース濃度は細胞外液の屈折率に大きな影響を与えており、組織中の等価散乱係数は細胞外液徒細胞との間の屈折率差に大きく依存するため、細胞外液の屈折率が変化して検出光に影響を及ぼす結果となる。したがって、検出光に基づいて屈折率分布を求めることによって組織内部の血糖値の分布を求めることが可能となるのである。
【０１３７】
また、上記実施形態においては複数の光入射及び光検出位置を散乱吸収体の一断面の周囲に配置しているが、図１１、図１２に示すように光入射及び／又は光検出位置（これらをＰで示す）を立体的に配置してもよい。すなわち、測定対象物として頭部や乳房を想定した場合は光入射及び／又は光検出位置（Ｐ）を図１１に示すように配置してもよく、測定対象物として腕部、脚部、胸部、乳房（圧迫時）を想定した場合は光入射及び／又は光検出位置（Ｐ）を図１２に示すように配置してもよい。
【０１３８】
また、上記実施形態においては時間積分法による光量測定値を測定値として用いているが、本発明に使用可能な測定値はこれに限定されるものではなく、例えば測定光の位相差（または位相遅れ）、振幅であってもよい。また、演算処理部６０において測定値を取得する具体的な手法も所望の測定値に応じて適宜選択され、例えば位相変調法による位相差及び／又は振幅測定、時間分解分光法による時間分解波形測定等の手段を採用してもよい。
【０１３９】
また、上記実施形態においては散乱吸収体ＳＭ内部の平均吸収係数および平均等価散乱係数を本発明にかかる光ＣＴ装置自身で得られたデータから求めているが、図１３〜図１５に示すように散乱吸収体ＳＭ内部の平均吸収係数および平均等価散乱係数を別の装置で求め（Ｓ１２０ａ）、それらに基づいて広がり関数を選択してもかまわない（Ｓ１３０ａ）。なお、図１３〜図１５における上記以外の諸工程は図８〜図１０における諸工程にそれぞれ対応する。この場合の利点は、例えば光ＣＴ装置で得られるデータはＣＷ（連続光）で測定し、平均吸収係数および平均等価散乱係数を求める装置でのみパルス光や変調光を用いればよくなるため、光ＣＴ装置のシステムの構成が簡単になる。なお、別の装置で平均吸収係数および平均等価散乱係数を求める手法は位相変調法や時間分解分光法であってもよい。
【０１４０】
このように測定対象物の内部の光学パラメータの分布を均一とみなして平均等価散乱係数μ's0と平均吸収係数μa0を測定する方法は、時間分解分光法については例えば“Devlopment of Time Resolved Spectroscopy System for Quantitative None-invasive Tissue Measurement ”M.Miwa et.al., SPIE vol. 2389 、位相変調法については例えば特開平６−２２１９１３号公報に記載されている。また、光ＣＴ装置に上記の手法を取り入れれば前記測定値の取得と同時に計算が可能である。
【０１４１】
また、光を生体等の散乱吸収体に入射する手段としては、図４および図１６（ｂ）に示した光ファイバーを使用する方法の他に、集光レンズ（図１６（ａ））、ピンホール（図１６（ｃ））を利用する方法、胃カメラのように体内から入射する方法（図１６（ｄ））等であってもよい。
【０１４２】
なお、生体試料等の散乱吸収体では平均拡散長が２ｍｍ程度であるため、入射光は約２ｍｍ直進するまでに散乱し、光の方向性がなくなる。したがって、数センチメートル以上の厚みの散乱吸収体では平均拡散長の影響を無視することができるのでスポット状に光を入射すれば良い。また、太いビーム状の光を散乱吸収体に入射してもよい。この場合には、複数のスポット状光源が並んでいると考えればよい。
【０１４３】
また、光入射ファイバーおよび光検出ファイバーと散乱吸収体ＳＭ表面との間の空間は、図４に示す実施形態では微小になっている。しかし実際には、これを大きくして、この空間に計測対象である散乱吸収体ＳＭとほぼ等しい屈折率および等価散乱係数をもつ液状体やゼリー状物体（以下、インターフェース材と呼ぶ）を満たしておいても良い。つまり、光はこのインターフェース材中を拡散伝搬して計測対象に入射するから何ら問題は生じない。また、散乱吸収体ＳＭの表面反射が問題になるときには、インターフェース材を適宜に選択することによって、表面反射等の影響を低減することができる。
【０１４４】
更に、上記実施形態においては異なる波長の光を時分割で入射させて使用する態様について説明したが、異なる波長の光を光混合器３５によって同軸状のビームにして、光入射点の直前に設けた波長選択フィルタ２０で波長選択して各波長の光を散乱吸収体に入射する方法、あるいはそのまま前記ビームを並行に散乱吸収体に入射する方法であってもよい（図７（ｂ））。但し、後者の場合は、図１７に示すように、光検出器４０ａ〜４０ｃの直前に波長選択フィルタ２５ａ〜２５ｃを設けて検出光を波長選択する必要がある。
【０１４５】
また、散乱吸収体の内部を拡散伝搬した光を受光して検出する手段としては、図４および図１８（ｂ）に示した光ファイバーを用いる方法以外に、直接検出する方法（図１８（ａ））、レンズを用いる方法（図１８（ｃ））等がある。
【０１４６】
また、光検出器４０で得られる信号を低雑音で増幅する必要がある場合には、狭帯域アンプ（図１９（ａ））、ロックインアンプ（図１９（ｂ））等を利用することができる。ロックインアンプを使用する場合には、参照信号として前記同期信号を利用する。この方法は、方形波光やパルス光を用いてダイナミックレンジの高い計測を行うときに有効である。
【０１４７】
さらに、上記実施形態においては複数の光入射ファイバーおよび光検出ファイバーを散乱吸収体の周囲に配置しておき、光入射および光検出に使用するファイバーをそれぞれ順次変えていくことによって光入射位置および光検出位置を移動させていたが、散乱吸収体に対する光入射位置と光検出位置とを同期して走査させてもよい。このようにすれば、散乱吸収体の各部の内部情報を求めて、フレームメモリーに蓄積し、これをテレビ方式で読み出すことによって内部情報の分布を示す画像が得られる。また、異なる時刻に計測すれば、内部情報の時間的変化を計測することができる。前記記憶部７０はこのようにして得た内部情報を記憶する機能をもち、表示部８０はこれらの途中経過や結果を表示するものである。この際、これらの演算処理は、メモリ７０、ディスプレイ８０等を備えるコンピュータ装置６０によって高速に実行することができる。
【０１４８】
さらに別な方法として、特開平６−２２１９１３号公報、特開平６−１２９９８４号公報に記載されているように、測定対象物の周囲に真円状のホルダーを取り付けて測定すれば、人の頭部を測定する際も個体差はそこで吸収され、実測の系をモデル系に近い状態で測定でき、精度が向上する傾向にある。
【０１４９】
また、上記実施形態においては既知の等価散乱係数より吸収係数の空間分布を算出したが、同一手法により、既知の吸収係数から等価散乱係数の空間的分布を求めることもできるため、したがって、今回提案する手法により吸収係数と等価散乱係数の両方の空間的分布を求めることが可能である。
【０１５０】
上記実施形態においては内部特性を求める式として光拡散方程式を用いたが、必ずしもこれに限定されるものではない。例えば、散乱吸収体内部の光吸収が伝播距離の関数として表わされることから導かれる式や、実験から得られる検出光と内部特性との関係式等を用いてもよい。
【０１５１】
また、病気や体の調子と検出光との相関が求まれば、本発明によって検出光から直接有用な情報を取得することが可能となる。例えば、検出光の変化が組織の構造変化と相関があれば、その相関を利用して検出光から構造変化を求めることができる。
【０１５２】
【実施例】
以下の実施例では、対象物がもつ全ての内部特性（吸収係数、等価散乱係数、屈折率等）の中で、吸収係数のみを変化させた場合における吸収係数の画像化の例を示す。
【０１５３】
本発明の有効性を実証するために、光入射ファイバーの配置間隔を２０度とし、光検出ファイバーの配置間隔を１０度とした以外は図４に示す装置と同様の装置を用いて以下の手順で実験を行なった。なお、本実施例で使用したファントムは高さ方向に同一形状であるため（ｚ軸方向に対称性を持つため）、３次元（立体）の問題を２次元（断面）の問題に次元を落として考えることができる。
【０１５４】
すなわち、波長８００ｎｍ、出力５０ｍＷのＣＷ光のレーザーを光入射ファイバーから図２０（ａ）、（ｂ）に示すファントムに入射し、ファントムを透過あるいは散乱反射した光を光検出ファイバーで検出し、検出器に導いた。なお、使用したファントムの仕様は以下の通りである。
【０１５５】
検出された光信号は光電子増倍管にて検出信号に変換し、その検出信号をカウンターで読み取って１０秒間加算した。加算された値はＧＰＩＢ（計測インターフェース・バス）でコンピュータに送った。
【０１５７】
ここまでの測定を１回とし、光入射ファイバーを図２０に示すＡ〜Ｒまで２０度づつ、光検出ファイバーを１０度づつ左回りに移動させてデータをとった。なお、１つの光入射位置に対する光検出位置は、それぞれの光入射位置を基準にして０度から９０度までの１０箇所とした。
【０１５８】
実験結果を以下に示す。
【０１５９】
１）別に作製した吸収体の含まれていないファントムで基準値Ｉd0を取得し、画像再構成を行った結果を図２２（ａ）に示す（背景は黒で示し、図２２（ｂ）は吸収係数の指標を示す）。図２２（ａ）に示す画像では複数個の像が見られ、かつ位置ずれが生じている。なお、図２２（ａ）中に記載された円は、期待される像の位置を示す。
【０１６０】
２）本発明にしたがって光入射−光検出の位置関係が相対的に同じである光入射−光検出位置の組み合わせ毎に測定値の平均値を求め、それを基準値Ｉd0をとして画像再構成を行った結果を図２３（ａ）に示す（背景は黒で示し、図２３（ｂ）は吸収係数の指標を示す）。本発明によると、図２３（ａ）に示した画像に見られるように、従来の方法によって得た図２２（ａ）に示す画像比べて位置ずれと定量性が改善されることが明らかである。なお、図２３（ａ）中に記載された円は、期待される像の位置を示す。
【０１６１】
【発明の効果】
本発明によれば、従来は必要であった物理モデルあるいはシミュレーションモデルから基準値を求めることなく、また測定対象物中の一成分に対して複数の波長を有する光を用いることなく、測定対象物についての測定値から直接的に基準値を求めることが可能となり、その基準値に基づいて測定対象物における内部特性分布（例えば吸収係数変化量分布、吸収係数分布、吸収成分濃度分布、等価散乱係数変化量分布、等価散乱係数分布、屈折率変化量分布、屈折率分布）を求めて画像化することが可能となる。そのため、本発明によれば、実際の測定対象物と物理モデルまたはシミュレーションモデルとの差や測定対象物の個体差等に起因する誤差の発生を避けることでき、信頼性の高いすなわち高精度の測定が可能となる。また、本発明によれば、物理モデル等を用いて基準値を予め求める手間が省かれ、測定時間の短縮化が可能となる。さらに、本発明によれば、複数の波長間での平均光路長分布および減衰光量（散乱等の影響で減衰する光量）は同じであるといった仮定に基づく誤差の発生が防止され、測定精度が向上する。
【図面の簡単な説明】
【図１】本発明の動作原理の説明図である。
【図２】吸収が均一な散乱吸収体のモデルを示す模式図である。
【図３】吸収が不均一な散乱吸収体のモデルを示す模式図である。
【図４】本発明の内部特性分布計測装置の一例を示す模式図である。
【図５】（ａ）および（ｂ）はそれぞれ光入射ファイバーの一例を示す斜視図および模式図である。
【図６】ヘモグロビンおよびミオグロビンの吸収スペクトルを示すグラフである。
【図７】（ａ）および（ｂ）はそれぞれ光入射手段の一例を示す模式図である。
【図８】本発明の内部特性分布計測方法の一例を示すフローチャートである。
【図９】本発明の内部特性分布計測方法の他の例を示すフローチャートである。
【図１０】本発明の内部特性分布計測方法の更に他の例を示すフローチャートである。
【図１１】本発明にかかる光入射及び／又は光検出位置の配置の一例を示す模式図である。
【図１２】本発明にかかる光入射及び／又は光検出位置の配置の他の例を示す模式図である。
【図１３】本発明の内部特性分布計測方法の更に他の例を示すフローチャートである。
【図１４】本発明の内部特性分布計測方法の更に他の例を示すフローチャートである。
【図１５】本発明の内部特性分布計測方法の更に他の例を示すフローチャートである。
【図１６】（ａ）〜（ｄ）はそれぞれ散乱吸収体への光入射方法の一例を示す模式図である。
【図１７】光検出手段の一例を示す模式図である。
【図１８】（ａ）〜（ｃ）はそれぞれ光検出方法の一例を示す模式図である。
【図１９】（ａ）および（ｂ）はそれぞれ検出信号の低雑音増幅方法の一例を示す模式図である。
【図２０】（ａ）および（ｂ）はそれぞれ実施例で使用したファントムの斜視図および上面図である。
【図２１】実施例における光入射位置と光検出位置との関係を説明するための補助図である。
【図２２】（ａ）及び（ｂ）は従来の方法で画像再構成を行った結果、ディスプレー上に表示された中間調画像を示す写真である。
【図２３】（ａ）及び（ｂ）は本発明の方法で画像再構成を行った結果、ディスプレー上に表示された中間調画像を示す写真である。
【符号の説明】
１〜１２…光ファイバーホルダー、１ａ〜１２ａ…光入射ファイバー、１ｂ〜１２ｂ…光検出ファイバー、２０…波長選択器、３０…光源、４０…光検出器、５０…制御部、６０…演算処理部、７０…記憶部、８０…表示部、ＳＭ…散乱吸収体、Ｐ…光入射及び／又は光検出位置。
Claims (26)
- 測定対象物の表面における複数の光入射位置から順次該対象物中に測定光を入射するステップと、
該対象物中を透過した測定光を、前記対象物の表面における複数の光検出位置のうちの少なくとも１つの光検出位置であってかつ測定されるべき測定光が入射された光入射位置に対して所定の位置関係にある光検出位置で順次あるいは同時に検出するステップと、
各光検出位置で検出された各測定光に基づいて、該測定光の所定パラメータの測定値を求めるステップと、
前記位置関係が相対的に同じである前記光入射位置と前記光検出位置との複数の組み合わせによって求められた複数の前記測定値を抽出し、該測定値の平均値を算出して該位置関係における基準値を得るステップと、
前記複数の組み合わせによって求められた前記複数の測定値と、前記基準値とを用いて、複数の領域に分割された前記対象物の各領域における所定内部特性の変化量を算出して該対象物における内部特性変化量分布を求めるステップと、
を具備する、内部特性分布の計測方法。 - 前記内部特性が吸収係数であり、
前記対象物の平均の吸収係数および平均の等価散乱係数を求めるステップと、
前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数を選択するステップと、
を更に具備しており、
前記内部特性変化量分布を求めるステップにおいて、前記複数の測定値と前記基準値と前記広がり関数とを用いて、前記各領域における吸収係数の変化量を算出する、請求項１記載の方法。 - 前記平均の吸収係数および平均の等価散乱係数を求めるステップにおいて、前記基準値に基づいて該平均の吸収係数および平均の等価散乱係数を求める、請求項２記載の方法。
- 前記吸収係数の変化量および前記平均の吸収係数を用いて、前記各領域における吸収係数の絶対値を算出して前記対象物における吸収係数絶対値分布を求めるステップを更に具備する、請求項２または３記載の方法。
- 前記吸収係数の絶対値を用いて、前記各領域における吸収成分の濃度を算出して前記対象物における吸収成分濃度分布を求めるステップを更に具備する、請求項４記載の方法。
- 前記対象物が少なくとも２つの吸収成分を含有しており、
前記光入射ステップにおいて前記対象物中に入射される測定光が、該吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有しており、
前記光検出ステップにおいて前記少なくとも２つの波長を有する測定光をそれぞれ検出し、
前記測定値を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記測定値を求め、
前記基準値を得るステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記平均値を算出し、
前記内部特性変化量分布を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の変化量を算出し、
前記吸収係数絶対値分布を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の絶対値を算出し、
前記吸収成分濃度分布を求めるステップにおいて前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収成分の濃度を算出し、前記対象物における前記各吸収成分の濃度分布を求める、請求項５記載の方法。 - 前記内部特性が等価散乱係数であり、
前記対象物の平均の吸収係数および平均の等価散乱係数を求めるステップと、
前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数を選択するステップと、
を更に具備しており、
前記内部特性変化量分布を求めるステップにおいて、前記複数の測定値と前記基準値と前記広がり関数とを用いて、前記各領域における等価散乱係数の変化量を算出する、請求項１記載の方法。 - 前記平均の吸収係数および平均の等価散乱係数を求めるステップにおいて、前記基準値に基づいて該平均の吸収係数および平均の等価散乱係数を求める、請求項７記載の方法。
- 前記等価散乱係数の変化量および前記平均の等価散乱係数を用いて、前記各領域における等価散乱係数の絶対値を算出して前記対象物における等価散乱係数絶対値分布を求めるステップを更に具備する、請求項７または８記載の方法。
- 前記内部特性が屈折率であり、
前記対象物の平均の吸収係数、平均の等価散乱係数および平均の屈折率を求めるステップと、
前記平均の吸収係数、平均の等価散乱係数および平均の屈折率に対応する広がり関数を選択するステップと、
を更に具備しており、
前記内部特性変化量分布を求めるステップにおいて、前記複数の測定値と前記基準値と前記広がり関数とを用いて、前記各領域における屈折率の変化量を算出する、請求項１記載の方法。 - 前記平均の吸収係数、平均の等価散乱係数および平均の屈折率を求めるステップにおいて、前記基準値に基づいて少なくとも該平均の吸収係数および平均の等価散乱係数を求める、請求項１０記載の方法。
- 前記屈折率の変化量および前記平均の屈折率を用いて、前記各領域における屈折率の絶対値を算出して前記対象物における屈折率絶対値分布を求めるステップを更に具備する、請求項１０または１１記載の方法。
- 前記の求められた分布に基づいて、前記対象物内部における該分布を示す画像を表示するステップを更に具備する、請求項１〜１２のうちのいずれか一項記載の方法。
- 測定対象物の表面における複数の光入射位置から順次該対象物中に測定光を入射する光入射手段と、
該対象物中を透過した測定光を、前記対象物の表面における複数の光検出位置のうちの少なくとも１つの光検出位置であってかつ測定されるべき測定光が入射された光入射位置に対して所定の位置関係にある光検出位置で順次あるいは同時に検出する光検出手段と、
各光検出位置で検出された各測定光に基づいて、該測定光の所定パラメータの測定値を求める測定値取得手段と、
前記位置関係が相対的に同じである前記光入射位置と前記光検出位置との複数の組み合わせによって求められた複数の前記測定値を抽出し、該測定値の平均値を算出して該位置関係における基準値を得る基準値算出手段と、
前記複数の組み合わせによって求められた前記複数の測定値と、前記基準値とを用いて、複数の領域に分割された前記対象物の各領域における所定内部特性の変化量を算出して該対象物における内部特性変化量分布を求める内部特性変化量算出手段と、
を具備する、内部特性分布の計測装置。 - 前記内部特性が吸収係数であり、
前記対象物の平均の吸収係数および平均の等価散乱係数を求める平均吸収および散乱係数検出手段と、
前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数を選択する広がり関数選択手段と、
を更に具備しており、
前記内部特性変化量算出手段において、前記複数の測定値と前記基準値と前記広がり関数とを用いて、前記各領域における吸収係数の変化量を算出する、請求項１４記載の装置。 - 前記平均吸収および散乱係数検出手段において、前記基準値に基づいて前記平均の吸収係数および平均の等価散乱係数を求める、請求項１５記載の装置。
- 前記吸収係数の変化量および前記平均の吸収係数を用いて、前記各領域における吸収係数の絶対値を算出して前記対象物における吸収係数絶対値分布を求める吸収係数絶対値算出手段を更に具備する、請求項１５または１６記載の装置。
- 前記吸収係数の絶対値を用いて、前記各領域における吸収成分の濃度を算出して前記対象物における吸収成分濃度分布を求める吸収成分濃度算出手段を更に具備する、請求項１７記載の装置。
- 前記対象物が少なくとも２つの吸収成分を含有しており、
前記光入射手段において前記対象物中に入射される測定光が、該吸収成分に対する吸収係数が互いに相違する少なくとも２つの波長を有しており、
前記光検出手段において前記少なくとも２つの波長を有する測定光をそれぞれ検出し、
前記測定値取得手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記測定値を求め、
前記基準値算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記平均値を算出し、
前記内部特性変化量算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の変化量を算出し、
前記吸収係数絶対値算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収係数の絶対値を算出し、
前記吸収成分濃度算出手段において前記少なくとも２つの波長を有する測定光に関してそれぞれ前記吸収成分の濃度を算出し、前記対象物における前記各吸収成分の濃度分布を求める、請求項１８記載の装置。 - 前記内部特性が等価散乱係数であり、
前記対象物の平均の吸収係数および平均の等価散乱係数を求める平均吸収および散乱係数検出手段と、
前記平均の吸収係数および平均の等価散乱係数に対応する広がり関数を選択する広がり関数選択手段と、
を更に具備しており、
前記内部特性変化量算出手段において、前記複数の測定値と前記基準値と前記広がり関数とを用いて、前記各領域における等価散乱係数の変化量を算出する、請求項１４記載の装置。 - 前記平均吸収および散乱係数検出手段において、前記基準値に基づいて前記平均の吸収係数および平均の等価散乱係数を求める、請求項２０記載の装置。
- 前記等価散乱係数の変化量および前記平均の等価散乱係数を用いて、前記各領域における等価散乱係数の絶対値を算出して前記対象物における等価散乱係数絶対値分布を求める等価散乱係数絶対値算出手段を更に具備する、請求項２０または２１記載の装置。
- 前記内部特性が屈折率であり、
前記対象物の平均の吸収係数、平均の等価散乱係数および平均の屈折率を求める平均吸収および散乱係数検出手段と、
前記平均の吸収係数、平均の等価散乱係数および平均の屈折率に対応する広がり関数を選択する広がり関数選択手段と、
を更に具備しており、
前記内部特性変化量算出手段において、前記複数の測定値と前記基準値と前記広がり関数とを用いて、前記各領域における屈折率の変化量を算出する、請求項１４記載の装置。 - 前記平均吸収および散乱係数検出手段において、前記基準値に基づいて少なくとも前記平均の吸収係数および平均の等価散乱係数を求める、請求項２３記載の装置。
- 前記屈折率の変化量および前記平均の屈折率を用いて、前記各領域における屈折率の絶対値を算出して前記対象物における屈折率絶対値分布を求める屈折率絶対値算出手段を更に具備する、請求項２３または２４記載の装置。
- 前記の求められた分布に基づいて、前記対象物内部における該分布を示す画像を表示する画像表示手段を更に具備する、請求項１４〜２５のうちのいずれか一項記載の装置。
Priority Applications (4)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP33467496A JP3662376B2 (ja) | 1996-05-10 | 1996-11-29 | 内部特性分布の計測方法および装置 |
| EP97303149A EP0806650B1 (en) | 1996-05-10 | 1997-05-08 | Method and apparatus for measuring internal property distribution |
| DE69733100T DE69733100T2 (de) | 1996-05-10 | 1997-05-08 | Vorrichtung und Verfahren zur Bestimmung der Verteilung einer inneren Eigenschaft |
| US08/853,879 US6075610A (en) | 1996-05-10 | 1997-05-09 | Method and apparatus for measuring internal property distribution |
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP8-140711 | 1996-05-10 | ||
| JP14071196 | 1996-05-10 | ||
| JP33467496A JP3662376B2 (ja) | 1996-05-10 | 1996-11-29 | 内部特性分布の計測方法および装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH1026585A JPH1026585A (ja) | 1998-01-27 |
| JP3662376B2 true JP3662376B2 (ja) | 2005-06-22 |
Family
ID=26473149
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP33467496A Expired - Lifetime JP3662376B2 (ja) | 1996-05-10 | 1996-11-29 | 内部特性分布の計測方法および装置 |
Country Status (4)
| Country | Link |
|---|---|
| US (1) | US6075610A (ja) |
| EP (1) | EP0806650B1 (ja) |
| JP (1) | JP3662376B2 (ja) |
| DE (1) | DE69733100T2 (ja) |
Families Citing this family (52)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5952664A (en) * | 1997-01-17 | 1999-09-14 | Imaging Diagnostic Systems, Inc. | Laser imaging apparatus using biomedical markers that bind to cancer cells |
| WO2000075633A1 (fr) * | 1997-12-12 | 2000-12-14 | Hamamatsu Photonics K.K. | Dispositif optique ct et procede de reformation d'images |
| JP4018799B2 (ja) * | 1998-04-02 | 2007-12-05 | 浜松ホトニクス株式会社 | 散乱吸収体の吸収成分の濃度計測方法及び装置 |
| JP3887486B2 (ja) * | 1998-05-26 | 2007-02-28 | 浜松ホトニクス株式会社 | 散乱吸収体の内部特性分布の計測方法及び装置 |
| ES2402233T7 (es) * | 1998-10-13 | 2015-04-23 | Covidien Lp | Oxímetro tisular multicanal no invasivo |
| DE69929380T2 (de) * | 1998-12-07 | 2006-08-31 | Hitachi, Ltd. | Vorrichtung zur Steuerung von Anlagen mit Hilfe von Körpersignalen |
| US7107116B2 (en) * | 1999-03-29 | 2006-09-12 | Genex Technologies, Inc. | Diffuse optical tomography system and method of use |
| US7142304B1 (en) * | 1999-09-14 | 2006-11-28 | The Research Foundation Of State University Of New York | Method and system for enhanced imaging of a scattering medium |
| CA2384839C (en) * | 1999-09-14 | 2013-12-24 | Randall L. Barbour | Imaging of scattering media using relative detector values |
| WO2001029534A1 (en) * | 1999-10-15 | 2001-04-26 | The Administrators Of The Tulane Educational Fund | Device for and method of simultaneously measuring light scatter from multiple liquid samples |
| JP4364995B2 (ja) * | 2000-03-21 | 2009-11-18 | 浜松ホトニクス株式会社 | 散乱吸収体内部の光路分布計算方法 |
| US20020033454A1 (en) * | 2000-08-04 | 2002-03-21 | Xuefeng Cheng | Optical imaging system with direct image construction |
| WO2002040971A1 (en) * | 2000-11-15 | 2002-05-23 | Bang & Olufsen Medicom A/S | Method and apparatus for determination of optical parameters of turbid media and use thereof |
| GB2371358A (en) * | 2001-01-22 | 2002-07-24 | Optokem Ltd | Light scattering particle characterisation apparatus and detection means |
| NO325061B1 (no) * | 2001-03-06 | 2008-01-28 | Photosense As | Fremgangsmate og arrangement for bestemmelse av den optiske egenskap av et multisjiktvev |
| US8506550B2 (en) * | 2001-09-07 | 2013-08-13 | Medtronic Minimed, Inc. | Method and system for non-vascular sensor implantation |
| US8465466B2 (en) * | 2001-10-23 | 2013-06-18 | Medtronic Minimed, Inc | Method and system for non-vascular sensor implantation |
| EP1514093B1 (en) | 2002-06-04 | 2021-04-28 | Visen Medical, Inc. | Imaging volumes with arbitrary geometries in non-contact tomography |
| US7349731B2 (en) * | 2002-06-07 | 2008-03-25 | Clemson University Research Foundation | Reconstructed refractive index spatial maps and method with algorithm |
| US6850314B2 (en) * | 2002-08-08 | 2005-02-01 | Board Of Reagents University Of Houston | Method for optical sensing |
| US6798508B2 (en) * | 2002-08-23 | 2004-09-28 | Coulter International Corp. | Fiber optic apparatus for detecting light scatter to differentiate blood cells and the like |
| US7736309B2 (en) * | 2002-09-27 | 2010-06-15 | Medtronic Minimed, Inc. | Implantable sensor method and system |
| US6954663B2 (en) * | 2003-01-07 | 2005-10-11 | Art Advanced Research Technologies Inc. | Continuous wave optical imaging assuming a scatter-law |
| WO2004072906A1 (en) * | 2003-02-05 | 2004-08-26 | The General Hospital Corporation | Method and system for free space optical tomography of diffuse media |
| JP2006521869A (ja) * | 2003-04-01 | 2006-09-28 | グルコン インク | 光音響分析評価方法と装置 |
| CA2520602A1 (en) * | 2004-09-21 | 2006-03-21 | Art Recherches Et Technologies Avancees Inc./Art Advanced Research Techn | Method for selecting wavelengths for optical data acquisition |
| US20070060808A1 (en) | 2005-09-12 | 2007-03-15 | Carine Hoarau | Medical sensor for reducing motion artifacts and technique for using the same |
| JP2008039539A (ja) * | 2006-08-04 | 2008-02-21 | Shimadzu Corp | 光散乱検出装置 |
| US10775308B2 (en) * | 2006-08-24 | 2020-09-15 | Xenogen Corporation | Apparatus and methods for determining optical tissue properties |
| 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 |
| JP4263767B2 (ja) * | 2007-09-20 | 2009-05-13 | 日東電工株式会社 | 充填用発泡組成物、充填発泡部材および充填用発泡体 |
| US8442608B2 (en) | 2007-12-28 | 2013-05-14 | Covidien Lp | System and method for estimating physiological parameters by deconvolving artifacts |
| US20090240138A1 (en) * | 2008-03-18 | 2009-09-24 | Steven Yi | Diffuse Optical Tomography System and Method of Use |
| US20090240139A1 (en) * | 2008-03-18 | 2009-09-24 | Steven Yi | Diffuse Optical Tomography System and Method of Use |
| US8292809B2 (en) | 2008-03-31 | 2012-10-23 | Nellcor Puritan Bennett Llc | Detecting chemical components from spectroscopic observations |
| US20100220315A1 (en) * | 2009-02-27 | 2010-09-02 | Beckman Coulter, Inc. | Stabilized Optical System for Flow Cytometry |
| JP5463545B2 (ja) * | 2009-03-31 | 2014-04-09 | セイコーエプソン株式会社 | 濃度定量装置、濃度定量方法及びプログラム |
| JP5465924B2 (ja) * | 2009-05-15 | 2014-04-09 | 宜彦 水島 | 不均一物質の光分析方法 |
| US8401608B2 (en) * | 2009-09-30 | 2013-03-19 | Covidien Lp | Method of analyzing photon density waves in a medical monitor |
| JP4892624B2 (ja) * | 2010-04-22 | 2012-03-07 | 株式会社日立製作所 | 生体光計測装置およびそれに用いる被検体装着具 |
| US9775545B2 (en) | 2010-09-28 | 2017-10-03 | Masimo Corporation | Magnetic electrical connector for patient monitors |
| EP2621333B1 (en) | 2010-09-28 | 2015-07-29 | Masimo Corporation | Depth of consciousness monitor including oximeter |
| JP2012237595A (ja) * | 2011-05-10 | 2012-12-06 | Sumitomo Electric Ind Ltd | 光トモグラフィ装置 |
| US20150253428A1 (en) * | 2013-03-15 | 2015-09-10 | Leap Motion, Inc. | Determining positional information for an object in space |
| US11493998B2 (en) | 2012-01-17 | 2022-11-08 | Ultrahaptics IP Two Limited | Systems and methods for machine control |
| USD763939S1 (en) | 2014-04-02 | 2016-08-16 | Cephalogics, LLC | Optical sensor array liner with optical sensor array pad |
| USD763938S1 (en) | 2014-04-02 | 2016-08-16 | Cephalogics, LLC | Optical sensor array |
| US10154815B2 (en) | 2014-10-07 | 2018-12-18 | Masimo Corporation | Modular physiological sensors |
| JP6570852B2 (ja) * | 2015-03-20 | 2019-09-04 | 株式会社東芝 | 生体成分推定装置、生体成分推定方法、およびプログラム |
| US10753862B2 (en) * | 2016-02-16 | 2020-08-25 | Stichting Het Nederlands Kanker Instituut-Antoni van Leeuwenhoek Ziekenhuis | Method, apparatus and software for detection and localization of hidden defects in optically diffuse media |
| WO2018089478A1 (en) * | 2016-11-11 | 2018-05-17 | Carrier Corporation | High sensitivity fiber optic based detection |
| EP3823525A4 (en) | 2018-07-16 | 2022-04-20 | BBI Medical Innovations, LLC | PERFUSION AND OXYGENATION MEASUREMENT |
Family Cites Families (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5137355A (en) * | 1988-06-08 | 1992-08-11 | The Research Foundation Of State University Of New York | Method of imaging a random medium |
| WO1989012223A1 (en) * | 1988-06-08 | 1989-12-14 | Barbour Randall L | A method of imaging a random medium |
| US5353799A (en) * | 1991-01-22 | 1994-10-11 | Non Invasive Technology, Inc. | Examination of subjects using photon migration with high directionality techniques |
| GB2228314B (en) * | 1989-02-16 | 1992-11-18 | Hamamatsu Photonics Kk | Examination apparatus |
| US5070455A (en) * | 1989-11-22 | 1991-12-03 | Singer Imaging, Inc. | Imaging system and method using scattered and diffused radiation |
| DE4341063B4 (de) * | 1992-12-09 | 2005-04-21 | Carl Zeiss | Vorrichtung und Verfahren zur optischen, ortsauflösenden Bestimmung von Dichteverteilungen in biologischem Gewebe |
| JP3433498B2 (ja) * | 1993-06-02 | 2003-08-04 | 浜松ホトニクス株式会社 | 散乱吸収体の内部情報計測方法及び装置 |
| JP3310782B2 (ja) * | 1994-07-14 | 2002-08-05 | 株式会社日立製作所 | 吸収物質濃度の空間分布画像化装置 |
| JP3433534B2 (ja) * | 1994-11-07 | 2003-08-04 | 浜松ホトニクス株式会社 | 散乱吸収体内の散乱特性・吸収特性の測定方法及び装置 |
-
1996
- 1996-11-29 JP JP33467496A patent/JP3662376B2/ja not_active Expired - Lifetime
-
1997
- 1997-05-08 EP EP97303149A patent/EP0806650B1/en not_active Expired - Lifetime
- 1997-05-08 DE DE69733100T patent/DE69733100T2/de not_active Expired - Lifetime
- 1997-05-09 US US08/853,879 patent/US6075610A/en not_active Expired - Lifetime
Also Published As
| Publication number | Publication date |
|---|---|
| DE69733100T2 (de) | 2006-03-09 |
| JPH1026585A (ja) | 1998-01-27 |
| EP0806650A3 (en) | 1998-03-18 |
| EP0806650A2 (en) | 1997-11-12 |
| EP0806650B1 (en) | 2005-04-27 |
| US6075610A (en) | 2000-06-13 |
| DE69733100D1 (de) | 2005-06-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP3662376B2 (ja) | 内部特性分布の計測方法および装置 | |
| US7356365B2 (en) | Method and apparatus for tissue oximetry | |
| JP3577335B2 (ja) | 散乱吸収体計測方法及び装置 | |
| JP3433534B2 (ja) | 散乱吸収体内の散乱特性・吸収特性の測定方法及び装置 | |
| JP3887486B2 (ja) | 散乱吸収体の内部特性分布の計測方法及び装置 | |
| JP3433498B2 (ja) | 散乱吸収体の内部情報計測方法及び装置 | |
| JP3107914B2 (ja) | 散乱吸収体内部の吸収情報計測装置及び方法 | |
| US8788003B2 (en) | Monitoring blood constituent levels in biological tissue | |
| US8050744B2 (en) | Optical system for measuring metabolism in a body and imaging method | |
| JP3107927B2 (ja) | 散乱吸収体の光学情報計測装置及び方法 | |
| US5596987A (en) | Optical coupler for in vivo examination of biological tissue | |
| US8280494B2 (en) | Apparatus and method to measure a spectroscopic characteristic in an object | |
| Kacprzak et al. | Time-resolved optical imager for assessment of cerebral oxygenation | |
| JPH07209177A (ja) | 散乱吸収体計測方法及び散乱吸収体計測装置 | |
| JPH1073481A (ja) | 散乱体の吸収情報計測方法及び装置 | |
| US20090240138A1 (en) | Diffuse Optical Tomography System and Method of Use | |
| WO2001019241A1 (en) | Calibration methods and systems for diffuse optical tomography and spectroscopy | |
| JP2004167080A (ja) | 酸素飽和度測定装置 | |
| CN101262812B (zh) | 光学层析成像装置 | |
| JPH11311569A (ja) | 内部特性分布の計測方法及び装置 | |
| US10682053B2 (en) | Ophthalmic photothermal optical coherence tomography apparatus | |
| JP4077476B2 (ja) | 散乱体の吸収情報計測方法及び装置 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20050304 |
|
| 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: 20050314 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20050323 |
|
| 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: 20090401 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090401 Year of fee payment: 4 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100401 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100401 Year of fee payment: 5 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110401 Year of fee payment: 6 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120401 Year of fee payment: 7 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130401 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130401 Year of fee payment: 8 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20140401 Year of fee payment: 9 |
|
| EXPY | Cancellation because of completion of term |