JP3730927B2 - 吸収物質濃度空間分布画像化装置 - Google Patents
吸収物質濃度空間分布画像化装置 Download PDFInfo
- Publication number
- JP3730927B2 JP3730927B2 JP2002078507A JP2002078507A JP3730927B2 JP 3730927 B2 JP3730927 B2 JP 3730927B2 JP 2002078507 A JP2002078507 A JP 2002078507A JP 2002078507 A JP2002078507 A JP 2002078507A JP 3730927 B2 JP3730927 B2 JP 3730927B2
- Authority
- JP
- Japan
- Prior art keywords
- light
- subject
- optical path
- spatial distribution
- equation
- 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 108
- 239000000126 substance Substances 0.000 title claims description 90
- 238000003384 imaging method Methods 0.000 title claims description 23
- 230000003287 optical effect Effects 0.000 claims description 220
- 238000001514 detection method Methods 0.000 claims description 123
- 238000000034 method Methods 0.000 claims description 99
- 238000005259 measurement Methods 0.000 claims description 80
- 239000011159 matrix material Substances 0.000 claims description 78
- 238000000342 Monte Carlo simulation Methods 0.000 claims description 47
- 239000011358 absorbing material Substances 0.000 claims description 25
- 238000004088 simulation Methods 0.000 claims description 8
- 230000001678 irradiating effect Effects 0.000 claims description 3
- 239000013307 optical fiber Substances 0.000 description 59
- 238000010521 absorption reaction Methods 0.000 description 36
- 238000010586 diagram Methods 0.000 description 29
- 230000008569 process Effects 0.000 description 23
- 238000004364 calculation method Methods 0.000 description 18
- 230000006870 function Effects 0.000 description 18
- 230000002745 absorbent Effects 0.000 description 16
- 239000002250 absorbent Substances 0.000 description 16
- 238000001228 spectrum Methods 0.000 description 16
- 238000012545 processing Methods 0.000 description 15
- 239000006096 absorbing agent Substances 0.000 description 12
- 230000008033 biological extinction Effects 0.000 description 12
- 108010054147 Hemoglobins Proteins 0.000 description 6
- 102000001554 Hemoglobins Human genes 0.000 description 6
- 230000008827 biological function Effects 0.000 description 6
- 239000008186 active pharmaceutical agent Substances 0.000 description 5
- 238000000149 argon plasma sintering Methods 0.000 description 5
- 238000007796 conventional method Methods 0.000 description 5
- 239000000975 dye Substances 0.000 description 5
- 210000001519 tissue Anatomy 0.000 description 5
- 210000000481 breast Anatomy 0.000 description 4
- 238000002059 diagnostic imaging Methods 0.000 description 4
- 239000000463 material Substances 0.000 description 4
- 230000008054 signal transmission Effects 0.000 description 4
- 230000005540 biological transmission Effects 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 3
- 238000003745 diagnosis Methods 0.000 description 3
- 239000000835 fiber Substances 0.000 description 3
- 230000017105 transposition Effects 0.000 description 3
- 102000036675 Myoglobin Human genes 0.000 description 2
- 108010062374 Myoglobin Proteins 0.000 description 2
- 206010028980 Neoplasm Diseases 0.000 description 2
- 239000012620 biological material Substances 0.000 description 2
- 230000001186 cumulative effect Effects 0.000 description 2
- 238000013016 damping Methods 0.000 description 2
- 239000003814 drug Substances 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000004907 flux Effects 0.000 description 2
- 230000005484 gravity Effects 0.000 description 2
- 238000001727 in vivo Methods 0.000 description 2
- 230000031700 light absorption Effects 0.000 description 2
- 210000003205 muscle Anatomy 0.000 description 2
- 230000003449 preventive effect Effects 0.000 description 2
- 230000005855 radiation Effects 0.000 description 2
- 230000004044 response Effects 0.000 description 2
- 239000004576 sand Substances 0.000 description 2
- 208000003174 Brain Neoplasms Diseases 0.000 description 1
- 208000014644 Brain disease Diseases 0.000 description 1
- 102000000634 Cytochrome c oxidase subunit IV Human genes 0.000 description 1
- 108050008072 Cytochrome c oxidase subunit IV Proteins 0.000 description 1
- 102000018832 Cytochromes Human genes 0.000 description 1
- 108010052832 Cytochromes Proteins 0.000 description 1
- 108010064719 Oxyhemoglobins Proteins 0.000 description 1
- RTAQQCXQSZGOHL-UHFFFAOYSA-N Titanium Chemical compound [Ti] RTAQQCXQSZGOHL-UHFFFAOYSA-N 0.000 description 1
- 238000000862 absorption spectrum Methods 0.000 description 1
- QVGXLLKOCUKJST-UHFFFAOYSA-N atomic oxygen Chemical compound [O] QVGXLLKOCUKJST-UHFFFAOYSA-N 0.000 description 1
- 239000008280 blood Substances 0.000 description 1
- 210000004369 blood Anatomy 0.000 description 1
- 230000003925 brain function Effects 0.000 description 1
- 201000011510 cancer Diseases 0.000 description 1
- 230000002490 cerebral effect Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000009792 diffusion process 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
- 230000010354 integration Effects 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 230000002503 metabolic effect Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000013421 nuclear magnetic resonance imaging Methods 0.000 description 1
- 229910052760 oxygen Inorganic materials 0.000 description 1
- 239000001301 oxygen Substances 0.000 description 1
- 238000006213 oxygenation reaction Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 239000010453 quartz Substances 0.000 description 1
- 230000004043 responsiveness Effects 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- VYPSYNLAJGMNEJ-UHFFFAOYSA-N silicon dioxide Inorganic materials O=[Si]=O VYPSYNLAJGMNEJ-UHFFFAOYSA-N 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 239000000725 suspension Substances 0.000 description 1
- 229910052719 titanium Inorganic materials 0.000 description 1
- 239000010936 titanium Substances 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Landscapes
- Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Description
【発明の属する技術分野】
本発明は、散乱物質中に吸収物質が存在してなる被検体において、吸収物質濃度の空間分布を求める方法に関する。特に、光を用いて生体内部あるいは大気中の吸収物質濃度の空間分布を計測する方法、装置に関する。
【0002】
【従来の技術】
現在、医療用の画像診断装置として、MRI(核磁気共鳴イメージング装置)X線CT、超音波診断装置等がある。それぞれ、臨床医療の現場で多用され有用性が認められている。現在、臨床医療に使用されている画像診断装置は、主に生体内の組織、構造を断層像として得るものが多く、生体内に分布する化学物質の代謝活動に起因する生体機能を計測するための簡便で実用的な画像診断装置は少ない。生体機能の計測により得られる情報は、脳疾患やガンの予防診断、脳機能の解明に有用と期待されており、簡便、実用的な生体機能画像診断装置の開発が医療の分野では待たれている。散乱物質中に、光の波長により吸収特性が異なる吸収物質が存在してなる被検体中の吸収物質の濃度の空間分布を画像化する光計測装置として、生体内の特定の色素濃度の空間分布を画像化する生体光計測装置がある。生体中には、生体機能に関与する吸収物質が多く存在し、吸収物質の濃度の空間分布を計測できる生体光計測装置は、生体機能を画像化して診断に有効となる装置の一つとして開発が急がれている。
【0003】
生体機能の中でも、生体内酸素飽和度は、生体中に多く含まれる特定色素の濃度に対応しており、この色素濃度は可視から近赤外領域における光吸収の量から計測することが可能である。特定色素として、酸化型あるいは還元型の別により、異なる光吸収の分光特性を持つことが知られている、血液中に存在するヘモグロビン(以下では、Hbと略記する)、細胞中に存在するチトクロームaa3、筋肉中に存在するミオグロビン等がある。さらに、光による計測は、光ファイバーや小型光源を用いることにより、簡便、小型な装置として開発される可能性を持っている。
【0004】
可視から近赤外の光を用いた生体機能を計測する装置が、例えば、特開昭57−115232号公報、特開昭63−275323号公報に記載されている。生体光計測装置で重要な課題の1つとして、生体内に存在する吸収物質の空間分布を画像化する方法の確立が挙げられる。この課題は吸収物質の濃度と位置を求める問題である。従来技術での画像化方法として、(1)X線CTのアルゴリズムを用いる方法、(2)輸送方程式を用いる方法、(3)モンテカルロシミュレーションの結果を用いる方法、に大別される。
【0005】
(1)の方法の具体的例として、時間ゲート法を用いて得た計測データから、X線CTにおいて使用されてきた画像再構成のアルゴリズムにより画像を得る方法(特開平01−209342号公報に記載)が挙げられる。
(2)の方法の具体的例として、アリッジ( Aridge )等の方法(" Reconstruction methods for infrared absorption imaging ", Proc. SPIE 1431, 204-215(1991)に記載)や、シンガー( Singer )等の方法(" Image reconstruction of the interior of bodies that diffuse radiation ", Science, 248, 990-993(1990)に記載)が挙げられる。両者の方法共に、散乱と吸収現象を表わす数学モデルに基づき、照射光強度と散乱を表わす変数(以下では、散乱変数と略記する)と吸収物質の吸収を表わす変数(以下では、吸収変数と略記する)とから、輸送方程式を立て検出光強度を計算する。ここで、散乱変数と吸収変数は散乱体内部の位置の関数であり、また、吸収変数は吸収物質の濃度を反映する値である。上記の輸送方程式の散乱変数と吸収変数の値を散乱体内部の任意の位置で独立に変化させ、計算された出光強度と実際に計測された検出光強度との差を評価し、その差が最小となるまで計算を繰り返す。この繰り返しの結果得られた、吸収変数の空間分布が吸収物質濃度の空間分布を表わしており、その分布を画像化して表示する。アリッジ( Aridge )等の方法では、散乱変数として散乱係数を、吸収変数として吸収係数を用いており、シンガー( Singer )等の方法では、散乱変数として散乱体内に設定された面積要素もしくは体積要素から出て行く光子の確率を、吸収変数として要素内で吸収される光子の確率を用いている。
【0006】
(3)の方法の具体的例として、バーバー ランダル エル等によるモンテカルロ法を利用した画像化方法がある(特表平3−505922号公報に記載)。この方法の概要は以下の通りである。媒体中の全ての点(位置及び方向)に対し、フラックスによって与えられる重みを割当て、その点で検出器の吸収体の存在しない媒体に対する応答に対する光子の期待される寄与を割り当てることによって、複雑な媒体を画像化する技術に関する。照射される輻射線の測定された強度が、予定される強度と比較されることにより、対応する散乱体積の透過係数が決定される。予定される強度は、予め規定されたモデル媒体に対するモンテカルロシュミレーションを用いた計算により得られ、透過強度係数が重み関数に適用される。媒体中における微小な吸収の検出器応答特性に対する影響は、影響を受ける体積要素に対する吸収散乱断面積と対応する重み関数との積となるものと仮定される。重み関数は、モンテカルロ法を用いて計算され、領域内の全ての点の、与えられた光源−検出器の配置に対する検出器応答性への相対寄与の3次元マップを与える。画像化は、透過係数と、光源−検出器の各配置ごとに、各体積要素に対応する重み関数との積を計算することにより達成される。これらの積の総和は、光源−検出器の配置のすべてに対するすべての体積要素の重ね合わせに対応する。この計算は、全ての光源の位置の関数として、全ての領域にわたる相対吸収断面積の3次元マップを与える。
【0007】
【発明が解決しようとする課題】
生体光計測装置は、Hb(ヘモグロビン)のように、生体内で機能する吸収物質濃度の空間分布を、画像として表示して、例えば癌の様な疾病の予防診断への応用が期待されている。しかし、生体は光に対して強い散乱特性を有しているため(散乱係数は1.0〔mm- 1〕程度である)、生体内部の吸収物質濃度の空間分布の精度良い画像化は非常に困難である。生体光計測装置を臨床医療で実用化するためには、生体の様な散乱体の内部に存在する吸収物質濃度の空間分布を精度良く画像化する方法が重要である。現在、このような生体内部の吸収物質濃度の画像化方法として、実用化レベルに達していると考えられる方法はまだ存在しない。従来技術における方法では、様々な問題があり、上記(1)の方法では、3次元的な散乱現象を2次元空間上で扱うため再構成された画像の空間分解能が低い。上記(2)の方法では、繰り返し計算をするため長時間の演算が必要で、現状の計算機速度では高分解能の画像を得ることは実質的に困難であり、また局所最小値( local minimun )に落ちこみ真の最適解が得られにくい場合がある。上記(3)の方法では、不透明な媒体中にある標的物体を3次元的に画像化することの記載はあるが、画像化する具体的な手順は数式的に明確には示されておらず、さらに媒体中の吸収物質の濃度を計算することの記載はない。
【0008】
本願発明の目的は、3次元的な散乱現象を考慮して、特に生体内で機能する吸収物質濃度を精度良く、短時間に空間分解能で求め、吸収物質の濃度の3次元空間分布を求める方法、及び装置を提供することにある。
【0009】
【課題を解決するための手段】
本願発明は、被検体を通過して検出された光強度が、被検体中に存在する吸収物質の濃度と、散乱に基づく光子の平均光路長との線形結合により表わされると仮定し、平均光路長の空間分布を使用して、被検体中に存在する吸収物質の濃度の3次元分布を求めることに特徴を有する。即ち、光吸収物質を含んでなる光散乱体からなる被検体に、所定の光照射位置から所定の波長を有するパルス光又は連続光を照射する照射ステップと、被検体を通過する光の強度を所定の光検出位置で検出する検出ステップと、を有し被検体内に分布する吸収物質濃度の空間分布画像化方法において、光照射位置から光検出位置までの複数の光路を求めるステップと、複数の光路から平均光路長を求めるステップと、平均光路長と、被検体に照射される光の強度と、光検出位置で検出される検出光強度と、被検体に照射される光に対する吸収物質の光学定数とから、各体積要素毎の吸収物質の濃度を求めることを特徴とする。平均光路長の空間分布は、モンテカルロ法により模擬して計算により求める。この計算では、被検体と実質的に同一の散乱係数を有し、仮想被検体を包含する大きさの仮想散乱体内の一点から発する複数光子のそれぞれのたどる複数光路を求める。仮想被検体は被検体と同一又は類似した形状を有し吸収物質を含まないものとする。
【0010】
より詳細に説明すると、本願発明の吸収物質濃度の空間分布画像化方法では、被検体と同一又は類似した形状を有し、吸収物質を含まない仮想被検体を想定して、仮想被検体において設定される光照射位置と光検出位置の複数の対のそれぞれで、光子がたどる光照射位置から光検出位置に至るまでの複数の光路を求めるステップと、複数の光路から光照射位置から光検出位置までの仮想被検体内における平均光路長を求めるステップと、平均光路長と、被検体に照射される照射光強度と、光検出位置で検出される検出光強度と、被検体に照射される光に対する吸収物質の光学定数(吸光係数)と、から吸収物質の濃度の空間分布を求めるステップと、吸収物質の濃度の空間分布を表示するステップと、を有することに特徴を有する。被検体の形状を表す空間は複数の体積要素に分割され、各体積要素ごとにに平均光路長が求められ、各体積要素内の吸収物質濃度が一様であると仮定し、各体積要素毎の吸収物質濃度が求められる。
【0011】
より具体的な吸収物質濃度の空間分布画像化のステップは、上記対における照射光強度と検出光強度との比に基づいて計測値行列を生成するステップと、上記対における平均光路長に基づいて平均光路長の空間分布行列を生成するステップと、計測値行列と平均光路長の空間分布行列から吸収物質濃度を求めるステップ、とからなる。予め与えられている所定の形状と所定の散乱係数を有する被検体に本発明の手順を適用する場合は、平均光路長の空間分布行列を生成するステップで求められ記憶手段に予め記憶されている平均光路長の空間分布行列を使用する。
【0012】
光路を求めるステップは、被検体と実質的に同一の散乱係数を有し、仮想被検体を包含する大きさの仮想散乱体内の一点から発する光子のたどる光路を模擬するモンテカルロ法により模擬して求められ、モンテカルロ法では、光子が順次散乱される複数の散乱点の内の連続する2つの散乱点のそれぞれの空間座標位置を関連づける散乱角度特性を表す関数、及び乱数を使用する。モンテカルロ法により求められた複数の光路のそれぞれを表す複数の散乱点の空間座標位置が記憶され、複数の光路データが記憶手段に記憶される。光路データは、仮想散乱体内の一点から発する単一光子のたどる有限光路長を有する光路データを求めるか、あるいは仮想被検体の2点を結ぶ最大距離以上の半径を有する球を仮想散乱体として、球の中心から複数光子を発して、複数光子のそれぞれの光子のたどる光路のうち、球の内部に存在する光路データを求める。
【0013】
仮想被検体での光照射位置から光検出位置までの光路は、仮想被検体での光照射位置と仮想散乱体内の1点とを一致させ、光検出位置とモンテカルロ法により求められた光路とが交差する点を、光路を表す散乱点の空間座標、又は仮想被検体の形状を表す空間座標のいずれか一方の空間座標を、他方の空間座標に対して座標変換(回転移動及び/又は並進移動)して求める。求められた仮想被検体での光照射位置から仮想被検体での光検出位置までの複数光路のそれぞれのを、被検体での対応する光照射位置から光検出位置までの光路とする。なお、仮想被検体の形状を表す空間は、被検体の形状を表す空間を複数の体積要素に分割したと同様に体積要素に分割される。また、光照射位置と光検出位置の対の数は、被検体、仮想被検体を分割した体積要素の総数以上の値に設定される。上記の複数光路データから仮想被検体の各体積要素を通過する平均光路長を体積要素ごとに求め、被検体の各体積要素を通過する平均光路長とする。光照射位置(i)と光検出位置(j)の対(i−j)に対して得られた光路に関するデータに基づいて、光照射位置と光検出位置の位置関係が入れ替わった光照射位置(j)と光検出位置(i)の対(j−i)に対する光路に関するデータを、対(i−j)に対して得られた光路に関するデータから、座標変換により求めることができる。
【0014】
照射ステップにおいて被検体に照射される光は、400nmから2000nmの範囲にある複数の波長を有する光であり、複数の波長の差が100nm以下とする。被検体を通過する光の強度が時間分解計測される場合には、検出光強度は時間分解計測された光強度の所定の時間内での光強度の積分値から求められ、所定の時間内に光照射位置から光検出位置に到達する光子に関する光路から求められる平均光路長に基づいて平均光路長の空間分布行列が生成される。
【0015】
本発明では、計測された検出光強度から計測値行列Yを生成し、散乱体内部の吸収物質の吸光係数と、吸収物質の存在しない散乱体を模擬して得る散乱体内における光が通過する光路データをモンテカルロ法により求め、光路データから、平均光路長の空間分布行列Aを生成し、行列Yと行列Aとから、未知数行列Xを導出して、求められた未知数行列Xから被検体内部の吸収物質濃度を求めこの空間分布を表示する。さらに、異なる複数位置に配置される光照射位置と光検出位置の対の数は、被検体、仮想被検体を分割した体積要素の総数以上の値に設定するので、ある体積要素に注目するとき、複数方向からの散乱の影響が考慮されるので、従来方法に比較して、約3倍以上の空間分解能を有する3次元的な吸収物質の空間分布を得ることができ、10〜100倍程度計算時間が短縮され、さらに定量性の改善が可能である。本方法では、局所の最適解に落ち込むことなく解を得ることができる。さらに、モンテカルロ法による光路の計算が、予め計算された光路データが記憶保存されている場合には、記憶保存されている光路データを使用するので、従来のモンテカルロ法を適用する方法と比較して、平均光路長の空間分布を約100倍以上の速度で求めることができる。
【0016】
【発明の実施の形態】
従来技術として、光を用いて生体内部に分布する特定色素濃度の空間分布を計測する生体光計測装置があり、また大気中に存在する吸収物質濃度の空間分布の計測として、レーザーレーダがある。この様に、散乱物質中に存在する吸収物質の濃度の空間分布を求め画像化する方法を必要とする分野は広い。特に以下では、生体光計測装置を代表的な例として説明をする。
【0017】
本願発明は、光計測装置で計測された検出光強度から、被検体内部に存在する吸収物質の濃度を画像化する方法に係り、本願発明における手順は要約すると、計測された検出光強度から計測値行列Yを生成するstep1、散乱体内部の吸収物質の吸光係数と、吸収物質の存在しない散乱体を模擬して得る散乱体内における光が通過する光路データをモンテカルロ法により求め、光路データから、平均光路長の空間分布行列Aを生成するstep2、行列Yと行列Aとから未知数行列Xを導出するstep3、求められた未知数行列Xから被検体内部の吸収物質濃度を求めこの空間分布を表示するstep4からなる。モンテカルロ法による光路の計算が、step2の処理開始以前に既になされ、光路データが記憶保存されている場合には、記憶保存されている光路データを使用する。
【0018】
まず最初に、被検体によって照射光が散乱されない場合(即ち、被検体が非散乱体である場合)に関して説明し、次いで、被検体によって照射光が散乱される場合(即ち、被検体が散乱体である場合)に関して説明する。図1に示す様な、光源1から被検体2へ所定の波長の光を照射し、光路3を通過してきた光の強度を検出器4で検出する光計測系を考える。図1は、光路3を含む被検体2の断面を示す。被検体2中に1成分の吸収物質が一様濃度cで存在する場合に、光源1からの照射光強度をI0、検出器4で検出される検出光強度をId、吸収物質の吸光係数をε、光照射位置と光検出位置間の幾何学的距離で決まる光路長をL’とすると、(数1)で定義されるベール・ランバート則が成り立つ。
Id/I0=exp(−εcL’) …(数1)
図1と同様の光計測系を考え、図2に示すように、被検体2が、吸収物質濃度の異なる2個の体積要素5−1、5−2に分けられる場合に、体積要素5−1中での吸収物質濃度をc1、光路長をL1’、体積要素5−2中での吸収物質濃度をc2、光路長をL2’とし、光源1からの照射光強度をI0、体積要素5−1から5−2に入射する光の強度をIC、検出器4で検出される検出光強度Idとすると、(数1)より、体積要素5−1、体積要素5−2に関し、それぞれ(数2)、(数3)が成り立つ。
IC/I0=exp(−εc1L1’) …(数2)
Id/IC=exp(−εc2L2’) …(数3)
体積要素5−1から5−2に入射する光の強度ICは計測されないので、(数2)、(数3)を用いてICを消去すると、(数4)が成り立つ。
Id/I0=exp(−ε(c1L1’+c2L2’)) …(数4)
次に、被検体が散乱体である場合について記述する(以下の説明では、特に明記しない場合は、被検体は散乱体とする)。
【0019】
図3に示すような光計測系で、被検体2中に1成分の吸収物質が一様な濃度cで存在する場合に、光源1からの照射光強度をI0、検出器4で検出される検出光強度をId、吸収物質の吸光係数をε、被検体内部に吸収物質が存在しないときの光路長(以下光路長と略す)をL、散乱によって被検体2の外に出て行く成分を減衰定数DSとすると、(数5)が成り立つ( Victor Twersky,"Absorption and Multiple Scattering by Biological Suspensions ", J. Opt. Soc. Am., 60, 1089(1970)に記載)。
Id/I0=DS・exp(−εcL) …(数5)
ここで、光を光子として扱うと、光源1から照射された光子は被検体2中で散乱を繰り返して検出器4に到達する。光子は散乱の度に進行方向が変わるので、その光路は、図3中に示す光路3の様にジグザグを描く。なお、図3では、図1、図2と同様に被検体2は立体であり、簡単のために、2本の光路3が1つの断面内にある場合を例示しており、一般に、光路3は3次元的にジグザグを描く。散乱による進行方向の変化はランダムであり光路3は各光子によって異なり、(数5)中の光路長Lは複数個存在し一義的に決定できない。そこで、複数個存在する光路長の平均値(以下では、平均光路長と略記する)を近似値として用いる。検出器に到達する光子数をm本とすると、光路の数はm本存在し、m本ある光路のうちk番目の光路長をLKとして、平均光路長LAを(数6)で定義し、さらに、平均光路長LAを用いて(数5)を書き換えると(数7)になる。
LA=(ΣLK)/m …(数6)
(数6)において、加算Σはk=1〜mについて行なう。
Id/I0=DS・exp(−εcLA) …(数7)
図3に示す光計測系を、図4に示すように、被検体2を、n個の体積要素に分割して考え、光源から発した光子は、全ての体積要素を通過して検出器4に至るものと仮定する。さらに各体積要素における吸収物質濃度は一般に異なると考える。図4は、簡単のために、図3で例示したのと同様に、光源1から照射された光子が被検体2中で散乱を繰り返して検出器4に到達する光路が全て、厚さが1体積要素からなる1断面内にある場合を示し、各体積要素の中心部分に示される太い線分は、この1断面内に存在する複数個の光路を各体積要素に分割して考え、各体積要素において複数個の光路のそれぞれ進行方向での距離を加算し一定値で割算して得た光路長を示し(得られる光路長の最大長を体積要素の辺長さ以下とするために一定値で割算した)、光源1と検出器4を結ぶ方向に平行に示したものである。言いかえると、図4は、被検体2として2次元物体を想定し、被検体2を空間(面積)要素に分割して考え、光源1からの複数光子のそれぞれが、2次元平面内で散乱の度に進行方向を変えて光路を形成し進む場合に、各面積要素において複数光子の描くそれぞれの光路を、それぞれの光子の進行方向での距離を加算し一定値で割算して得た光路長として太い線分により、光源1と検出器4を結ぶ方向に平行に示したものである。吸収物質が存在せず、かつ一様な散乱体では、図4に示すように、この太い線分の長さの分布は光源1と検出器4を結ぶ直線に対して対称となり、この直線から距離が離れる程太い線分の長さは短くなる。しかし、一般に被検体2は一様な散乱体ではなくしかも、被検体2の各部位は異なる濃度の吸収物質が存在するため、このような場合には図4に示す太い線分の長さは、各要素で異なることになる。
【0020】
本願発明では、光源1からの複数光子のそれぞれは、被検体2を構成する全ての体積要素を一般に通過していくと想定し、各体積要素が図3に示す被検体であると見做して(数7)を適用し、さらに複数の異なる濃度の吸収物質が存在する場合に適用される(数4)を適用する。即ち、光源1からの複数光子のそれぞれは、散乱のためジグザグを描いて進行するが、このとき各光子は、全ての体積要素を通過して検出器4に至るものと仮定すると、(数4)、(数7)を考慮して(数8)が成立するものとする。但し、(数8)において、各体積要素中での、吸収物質濃度をcj、平均光路長をLAjとし、添字jは、j番目の体積要素内における吸収物質濃度、平均光路長の値であることを示し、加算Σはj=1〜nについて行なう。nは、体積要素の総数であり、図4に示す例では、5×5=25である。
なお、光源からの光子は、被検体を構成する全ての体積要素を一般に通過していくと想定するが、光子が通過しない体積要素では、LAj=0であり、体積要素毎の平均光路長LAj(j=1〜n)を求める方法に関しては後述するので、LAj(j=1〜n)は既知として扱う。(数8)において、減衰定数DSを直接求めることは困難であるため、照射光として波長の異なる2種類の光を用い、各波長毎に計測を行ないDSを消去する。波長λ1の照射光、波長λ2の照射光を用いた計測に関して、それぞれ(数9)、(数10)が成り立つ。ここで、平均光路長LAjは波長に依存しないものとする。DS、εは波長に依存する。
(数9)、(数10)において、加算Σはj=1〜nについて行なう。
【0021】
(数9)及び(数10)を(数11)及び(数12)を用いて変数変換をすると、(数13)及び(数14)が得られる。
例えば、被検体が生体組織の場合には、波長λ1、λ2の波長差の絶対値(│λ1−λ2│)が、少なくとも100nm、特に30nm以内であれば、近似的に次の(数15)が成り立つ。
と置いて、(数16)の両辺の自然対数(Loge)とり(数18)を得て、
さらに、
と置いて、
y=Σ(χ(j,λ1,λ2)LAj) …(数20)
を得る。(数20)において、加算Σはj=1〜nについて行なう。
【0022】
(数20)は、線形方程式であり、少なくともn組の同様の方程式を用意すればn個ある未知数χ(j,λ1,λ2)(j=1〜n)を求めることができる。1つの方程式は被検体に関する1つの計測位置(光源1と検出器4の位置関係)毎に成立するので、図5に示す様に、1つの光源1に対して検出器4を被検体2に複数配置して複数位置で検出するか、あるいは複数の照射位置から光を照射し、n個の検出位置で光検出を行うと、n組の(数20)の方程式が得られ、未知数χ(j,λ1,λ2)(j=1〜n)と同数の連立方程式(数21)を得ることができる。但し、iはi番目の光検出位置を示す添え字である。なお図5は、被検体2が2次元体である例を示す。
(数22)、(数23)、(数24)により、行列Y、A、Xを定義すると、(数21)で表される連立方程式は(数25)となる。以下の説明においては、Yを計測値行列、Aを光路長の空間分布行列、Xを未知数行列と呼ぶ。
Y=│y1…yi…yn│T …(数22)
(数22)において、Tは転置を表し、Yはn行1列の行列である。
A=│LAij│ …(数23)
(数23)において、Aはi行j列の要素の値をLAijとするn行n列の行列である。
(数24)において、Tは転置を表し、Xはn行1列の行列である。
Y=AX …(数25)
(数25)から未知数行列Xを求めるためには、(数26)で示すように、Aの逆行列を両辺にかければよい。
A-1Y=A-1AX=X …(数26)
ここで、求められた未知数行列Xの各成分は、(数27)で表される。
(数27)に照射光として用いた各波長λ1、λ2に対する吸光係数を代入して、各体積要素中の吸収物質濃度cjを得る。
【0023】
以上の説明では、被検体2をn個の体積要素に分割して、n個の検出位置(光源1(光照射位置)と検出器4の対)で光検出を行ない、n組の(数20)の方程式が得たが、nをn’≧nなるn’に置き換え、即ちn’個の検出位置(光源1(光照射位置)と検出器4の対)で光検出を行ない、n’組の(数20)の方程式が得ても、以下のようにして未知数行列Xを得ることができる。以上の説明において(数22)、(数23)を以下のように置き換えればよい。即ち、(数22)は、
とする。(数28)においてTは転置を表し、Yはn’行1列の行列である。
【0024】
(数23)は、
とする。(数29)において、Aはi行j列の要素の値をLAijとするn’行n列の行列である。このような置き換えを行なった(数25)の解は一般に、
X=((AT)WA)-1((AT)W)Y …(数30)
として得られる。なお(数30)において、Tは反転を表し、ATはAの反転行列であり、Wは重み行列を表すn’行n’列の行列である。重み行列Wをn’行n’列の単位行列とするとき(数30)は、
X=((AT)A)-1(AT)Y …(数31)
となる。
【0025】
このように、被検体をn個の体積要素に分割して、n’≧nなるn’個の各種方向での検出位置(光源1(光照射位置)と検出器4の対)で光検出を行なうことにより、精度よく未知行列X、即ち被検体の各体積要素中の吸収物質濃度cjを得ることができる。
【0026】
さらに以上の説明において、計測対象である被検体が多成分の吸収物質を含む場合には、照射光に用いる波長の種類を増やして、多成分の吸収物質の濃度を求めることができる。例えば、被検体がαとβの2種類の吸収物質を含む場合、(数5)を拡張して(数32)が成り立つ。但し、吸収物質α、βの分子吸光係数と濃度をそれぞれ、εα、εβ、cα、cβで表す。
Id/I0=DS・exp(−(εαcα+εβcβ)LA) …(数32)
図4に示すように、被検体2をn個の体積要素に分割する場合には、吸収物質が1種類ある場合と同様にして、(数32)から(数33)が得られる。
(数33)において、加算Σはj=1…nについて行なう。
【0027】
減衰定数DSを2種類の波長λ1、λ2の照射光を用いて消去すると(数34)が得られ、(数35)の置換を行なうと、(数34)は(数36)で表される。
同様にして2種類の波長λ1、λ3の照射光を用いて、(数36)と同様の式を得る。
【0028】
(数34)、(数36)において、加算Σはj=1…nについて行なう。(数36)の両辺の自然対数(Loge)をとると線形な式となるので、以下、先に説明した1種類の吸収物質がある場合と同様に、n’≧nなるn’個の検出位置(光源1(光照射位置)と検出器4の対)で光検出を行ない、(数36)を解くことができ、各体積要素内における変数χ(j,λ1,λ2)が求まる。ここで、波長λ1と波長λ3から得られる各体積要素内における変数χ(j,λ1,λ3)と前記の変数χ(j,λ1,λ2)により連立方程式(数37)が成り立つ。
ここで、吸収物質α、βの吸光係数が既知のとき、連立方程式(数37)に吸光係数の値を代入して、各体積要素内の2種類の吸収物質の濃度cα、cβを求めることができる。従って、被検体内部にγ種類の吸収物質が存在する場合、(γ+1)種類の波長の照射光を用いて、各吸収物質の濃度分布を求めることができる。
【0029】
以上の説明した被検体を複数の体積要素に分割した場合において、特に被検体を分割しない場合(即ち、分割の数をn=1とする)にも、以上の方法は適用できることは言うまでもない。この場合、(数29)はn’行1列の行列となり、異なる各種の方向位置に配置される光照射位置と光検出位置の対により計測されるn’個の計測値からなる計測値行列と、上記の対におけるそれぞれの平均光路長とから吸収物質の濃度が求められる。もちろん、被検体が複数種類の吸収物質を含む場合にも、被検体を分割せず分割の数をn=1として、(数37)を適用できる。このようにして求められる吸収物質の濃度は、被検体の2次元的空間(被検体の特定の断面の周囲に2次元的に上記の対が配置されるとき)、あるいは3次元的空間(被検体の外周に3次元的に上記の対が配置されるとき)での平均的な値が求められることになる。
【0030】
以上の説明において、行列Xを求めるためには、行列Aが既知でなければならない。行列Aを構成する成分LAijは、i番目の計測位置(光源と検出器の位置関係)におけるj番目の体積要素内の平均光路長を表している。即ち、LAijは、各計測位置毎に検出器に到達する光子の光路が既知であれば求めることができる。図6は、光源1から被検体2を通過し、検出器4に到達する光子の光路3を、被検体2が2次元体である例をとって模式的に表した図である。
【0031】
本発明では、光源から検出器に至る複数の光路を求める方法として、被検体内の光散乱を模擬するモンテカルロ法シミュレーションを用いる。被検体内の光散乱を模擬するモンテカルロ法シミュレーションの例として、B.C.Wilson等の方法(" A Monte Carlo model for the absorption and flux distributions of light in tissue ", Med. Phys., 10(6), 824ー830(1983)に記載)が知られている。本願発明で使用するモンテカルロ法は、よく知られた方法であるが、この方法の概要を以下に説明する。光源から発した1光子のたどる経路は各散乱点を表す、極座標(r,θ,φ)で示される。各散乱点では、1光子が次に達する散乱点を与える、r、θ、φの各座標値が確率的に定められ、散乱点間の距離は散乱距離LSから、光子の進行方向を定めるための角度方向は、θ方向の散乱角度θSとφ方向の散乱角度φSとから決定される。なお、光子がたどる物体中には吸収物質は含まれないものとする。各散乱毎の散乱距離LSは散乱係数μSを用いて(数38)で表わされ、
LS=−(loge(R))/μS …(数38)
である。Rは一様乱数(0<R<1)である。θ方向の散乱角度θSは、例えばMie散乱理論から近似して得られる位相関数の近似式である、Heney−Greensteinの式を用いて、一様乱数R(0<R<1)と位相関数( phase function )(即ち1回散乱における散乱強度の角度分布を示す関数)のコサイン平均gとを使用して(数39)で表わされ、
である。なお、gは、位相関数をPfとし、積分∫を立体角ω=0〜4πの範囲について行なうものとして、g=(∫Pf・cosθdω)/(∫Pfdω)により与えれる。一般の生体関連物質では、g≒0.9であることが知られている。特に、g=0の場合には等方的散乱を表し、(数39)は、
θS=cos-1(2R−1) …(数40)
となる。φ方向の散乱角度φSは一様な確率分布とみなすことができ、一様乱数Rを用いた(数41)で表わされる。
φS=2πR …(数41)
光子のたどる各散乱点は、このように、一様乱数(0<R<1)Rにより、確率的に与えられる(LS,θS,φS)により、順次与えられる。以上が、本願発明で使用するモンテカルロ法の概要である。
【0032】
このモンテカルロ法による光子のたどる経路シミュレーションを用いて、i番目の計測位置(光照射位置(光源)と光検出位置の組み合わせ)における1光路を求めた時に、j番目の体積要素内での光路長をp(i,j)とする。実際の被検体内での光の伝播を実質的に再現するに十分な光路数m本(被検体形状に依存するが、少なくとも1000本程度とする)を、上記で説明したモンテカルロ法シミュレーションを用いて計算し記憶しておき、各光路毎にp(i,j)を求める。k(k=1〜m)番目の光路におけるp(i,j)をp(i.j.k)とすると、LAijは(数42)で計算される。
LAij=(Σp(i,j,k))/m …(数42)
(数42)において、加算Σはk=1〜mについて行なう。全ての計測位置(i)と、全ての体積要素(j)に関してLAijを計算すれば、平均光路長の空間分布行列Aが求められる。
【0033】
以上説明した本発明の手順は、照射光として連続光を用い被検体を通過する光強度を計測信号とする場合、あるいは、照射光としてパルス光を用い、被検体を通過するパルス光に対応する全光強度を計測信号とする場合であるが、被検体を通過する光強度を時間分解計測し、計測された光強度の1部を計測信号とする場合にも、同様の手順で吸収物質の濃度分布を求めることができる。
【0034】
図7に示す様な時間スペクトルを有するδ関数状のパルス光を、被検体、例えば生体の所定部分に照射し、照射位置と対向する側の位置で通過光強度を、ストリークカメラを用いて時間分解計測すると、図8に示す様に時間的に拡がった通過光強度の時間スペクトルId(t)が得られる。ただし、横軸はパルス光が照射されてから通過光が検出されるまでの時間(検出時間)を表す。検出光強度の時間スペクトルId(t)が、時間的に広がる理由は、被検体が散乱体である場合には、各光子の光路長に応じて検出される時間が変化するためである。図6の様な系で、i番目の各計測位置毎に検出光強度の時間スペクトルIdi(t)を計測し、検出時間t1からt2までIdi(t)を積分した値を検出光強度Idiとして用い計測値行列Yを求め、検出時間t1からt2までの平均光路長の空間分布行列Aを求めれば、同様の手順で未知数行列Xを求めることができる。検出時間t1からt2までの平均光路長の空間分布行列Aは、検出時間t1からt2までの間に検出器に到達する光子の光路のみについて(数42)を用いて計算すれば求めることができる。このように検出時間t1からt2まで、時間ゲートを設定してIdi(t)を積分した値を求める際に、時間ゲート幅、(t2−t1)の大きさを、小さくすればする程、最終的に再構成される被検体中の吸収体の3次元空間分布の空間分解能、及び定量性を向上させることが可能となる。
【0035】
以下、以上説明した内容を実際に適用した実施例について詳細に説明する。まず簡単のため、本発明を、生体の散乱特性を模擬した立方体の被検体に適用した例を示す。なお、本願発明は被検体が立方体の形状に限定されるものではなく、任意の形状に適用できることは言うまでもない。任意形状の被検体では、被検体の形状を小さい立方体の集合として近似して表してもよいし、被検体の内部を小さい立方体の集合とし、被検体の外面を含む部分を、一部を欠損する小さい立方体で近似して表してもよい。さらに、検査対象物体の特定部分のみの空間、例えば、特定組織部分を含む空間部分、生体表面もしくは表面近傍にある腫瘍を含む空間部分等を、被検体部位として上記と同様の方法で表すことも可能である。
【0036】
図9は、本実施例に用いた被検体の形状を示す図である。被検体2は1辺25mmの立方体で、被検体内部の位置を表す座標系(X0,Y0,Z0)は図9に示す通りである。被検体2の各面を図9に示す様に、P1〜P6とする。ここで、P1はZ0=0(mm)の面、P2はX0=25(mm)の面、P3はZ0=25(mm)の面、P4はX0=0(mm)の面、P5はY0=25(mm)の面、P6はY0=0(mm)の面である。被検体2の散乱特性として、(数38)で使用する散乱係数を生体に近い値であるμS=1.0〔mm-1〕とし、(数39)で使用する先に説明した位相関数のコサイン平均は、g=0とした。即ち、1回散乱における散乱光強度の角度分布は一様であり等方的であるとした。
【0037】
また、図10に示す様に、被検体2内部を5×5×5の仮想の体積要素に分割し、それぞれの体積要素の位置を(Vx,Vy,Vz)で表し、各体積要素番号jは(数43)で定義する。例えば、図10の体積要素5の位置は(5,5,1)であり、体積要素番号j=25となる。
j=Vx+5(Vy−1)+25(Vz−1) …(数43)
図11は、照射光7を図9のP1の面上にある光照射位置6−1から照射していることを示している。本実施例では、全ての光照射位置、及び光検出位置は、被検体表面に接する任意の体積要素の面の重心位置に設定されている。図9の様に、光を照射するi番目の面Piに固有の座標(X,Y)Piを定義すると、光照射位置6−1は(2,2)P1、光照射位置6−2は(2,2)P2、光照射位置6−3は(3,3)P5で表される。また、本実施例では、光検出位置を光照射位置が対向する側の被検体外面にある全ての体積要素面の重心位置に設けた。例えば、図11において、P1面上の光照射位置6−1から照射光7を照射した時には、対向側のP3面上に25点の光検出位置を設けており、それぞれの光検出位置で検出光強度が計測される(面Piの定義については図9を参照)。本実施例では、被検体2を合計125個の体積要素で分割しているので、それぞれの体積要素内で一様な吸収物質濃度を未知数とすると、求めるべき125個の未知数があることになる。未知数を決定するためには、光計測位置(光源と光検出器の位置関係)も125組だけ必要となる。そこで、本実施例では、照射点を(2,2)P1、(2,2)P2、(4,4)P3、(4,2)P4、(3,3)P5の5点とし、光検出位置は各光照射位置の対向側の面上に25点設定し、合計125(5×25)組の光計測位置を用意した(面Piの定義については図9を参照)。
【0038】
本実施例では、被検体として、吸収物質の位置の異なる2通りを用意した。これら被検体内部の吸収物質の位置を図12と図13に示す。図12と図13における吸収物質8の位置を(Vx,Vy,Vz)で表すと、図12の被検体では、吸収物質が(3,3,3)の体積要素内に均一濃度で存在し、図13の被検体では、吸収物質が(3,3,3)と(1,3,1)の体積要素内に均一濃度でそれぞれ存在している。照射光として2波長λ1、λ2のパルス光を用い、各検出位置において通過光強度を時間分解計測した。波長λ1、λ2のそれぞれに対する吸収物質の吸光係数は、ε(λ1)=1.0〔(mm・mM)-1〕、ε(λ2)=0.0〔 (mm・mM)-1)〕とし、吸収物質濃度は全て0.01〔mM〕とした。図14に、先に説明したモンテカルロシュミレーションにより得られた、透過光強度(検出光強度)の時間スペクトルIdi(λ1,t)、Idi(λ2,t)の1例を示す。横軸は検出時間を表しており、光子が散乱されずに直接検出位置に到達する時刻を原点としている。なお、図7、図8と同様に、図14の横軸は、それぞれの横軸の値に、光子の進行する媒体での光速度を乗算して、光子の進行距離として示すこともできる。本実施例では、i番目の計測位置における検出光強度Idi(λ1)、Idi(λ2)として、各計測位置(i)毎に得られる検出光強度の時間スペクトルを検出時間t1=0〔ps〕から検出時間t2=1500〔ps〕まで積分した値を用いた。なお、t2を1500〔ps〕より、小さい値に設定することにより、最終的に得られる吸収体の3次元分布の空間分解能、及び定量精度は、更に向上する。
【0039】
計測された検出光強度から、被検体内部に存在する吸収物質の濃度を画像化する手順の概要を図15に示す。この画像化の手順は大きく分けると、計測値行列Yを生成するstep1、平均光路長の空間分布行列Aを生成するstep2、行列Yと行列A逆行列から未知数行列Xを導出するstep3、求められた未知数行列Xから被検体内部の吸収物質濃度を求めこの空間分布を表示するstep4からなる。なお、モンテカルロ法による光路の計算が、step2の処理開始以前に既になされ、光路データが記憶保存されている場合には、記憶保存されている光路データを使用することは言うまでもない。以下、各stepにおける詳細な手順を、図16と図17にフローチャートにより説明する。
【0040】
〔step1〕: 計測値行列Yの生成。
【0041】
step1−1: 各計測位置(i)毎に通過光強度を計測して、それぞれ検出光強度Idi(λ1)、Idi(λ2)を求める。
【0042】
step1−2: 各計測位置(i)毎に照射光に用いた2波長の照射光強度Id0(λ1)、Id0(λ2)を計測する。
【0043】
step1−3: 全計測位置(i)に関して計測終了後、step1−1、及びstep1−2で求められた検出光強度Idi(λ1)、Idi(λ2)、及びId0(λ1)、Id0(λ2)を(数19)に代入してyiを得る。
【0044】
step1−4: step1−3で得たyi(i=1〜125)から計測値行列Y(数22)を生成する。
【0045】
〔step2〕: 平均光路長の空間分布行列Aの生成。
【0046】
step2−1: 計測対象である被検体の3次元的な外形形状を計測して、外形形状を表す座標データを演算処理装置に入力し、仮想被検体データを生成する。外形形状を計測する方法としては、光切断法、MRI、X線CT等その他方法を用いる。本実施例では被検体は1辺25mmの立方体であり、この外形形状データを入力する。
【0047】
step2−2: 計測時に設定される計測位置(光照射位置と光検出位置)の座標を入力する。本実施例では先に説明した125組の計測位置を演算処理装置に入力する。
【0048】
step2−3: 計測対象となる被検体の散乱特性として、被検体の散乱係数μS、及び位相関数のコサイン平均gを演算処理装置に入力する。本実施例では、散乱係数μS=1.0〔mm-1〕、指向性を表わす変数としてg=0を用いる。
【0049】
step2−4: step2−1〜step2−3で入力したパラメータを基に、モンテカルロ法を用いて吸収物質の存在しない被検体内の光散乱を模擬する計算を行なう。設定された各計測位置(i)毎に複数の光路を計算し、各計測位置(i)毎の各体積要素(j)毎における平均光路長LAij(平均光路長の空間分布)をメモリに記憶しておく。この時、通過光強度を時間分解計測する場合には、計測信号として用いる検出時間t1からt2の間に到達する光子の光路のみから、平均光路長の空間分布LAijを求めて、メモリに記憶しておく。本実施例では、各計測位置毎に30000本の光路を計算し、検出時間は0〔psec〕から1500〔psec〕とした。本ステップで用いるモンテカルロ法に関しては先に概要を説明した。なお、モンテカルロ法による光路の計算が、step2−4の処理開始以前に既になされ、光路データが記憶保存されている場合には、記憶保存されている光路データを使用することは言うまでもない。光路データが記憶保存されている場合の、各計測位置(i)毎の各体積要素(j)毎における平均光路長LAijを求めるための処理方法に関しては、後で詳細に説明する。
【0050】
step2−5: step2−4で全計測位置(i)と全体積要素内(j)に関して平均光路長LAij(i=1〜125、j=1〜125)を求めた後、平均光路長の空間分布行列A(数23)を生成する。ただし、計測対象である被検体の形状が予想される場合、あるいは被検体の形状と類似した形状を用いる場合には、あらかじめ平均光路長の空間分布行列Aを計算しておく。例えば、生体頭部や乳房が被検体である場合には、複数種類の頭部形状や乳房形状について、平均光路長の空間分布行列Aを計算しておき、実際に計測対象となる生体頭部や乳房と最も近い形状に関して計算した平均光路長の空間分布行列Aを選択する。また、生体頭部等は、球体あるいは楕円体に近似して計算する。
【0051】
〔step3〕: 未知数行列Xの導出。
【0052】
step3−1: step1−4で生成された計測値行列Yとstep2−5で生成された平均光路長の空間分布行列Aの逆行列A-1を求め(数26)に代入して、未知数行列Xを求める。もちろん(数30)により未知数行列Xを求めてもよい。
【0053】
step3−2: 未知数行列Xの各成分(数27)より、照射光の各波長に対する吸収物質の吸光係数を用いて、各体積要素毎の吸収物質濃度を求める。
【0054】
〔step4〕: 画像表示。
【0055】
step3で得られた吸収物質濃度の空間分布から、任意断層面の吸収物質濃度分布を画像表示する。
【0056】
図18から図23に、上記で説明したフローチャートに従って求められた吸収物質濃度の空間分布の結果を表示する。図18から図20は、図12に示した被検体における吸収物質濃度の空間分布を求めた結果である。図18はVy=2の面、図19はVy=3の面、図20はVy=4の面を表示している。図21から図23は、図13に示した被検体における吸収物質濃度の空間分布を求めた結果である。図21はVy=2の面、図22はVy=3の面、図23はVy=4の面を表示している。図18から図23において、x=、y=、z=のそれぞれに続いて記される数は対応する座標の値を示す。以上の計算では、被検体の散乱係数として、生体に近い散乱係数、μ=1.0〔mm-1〕を設定し、(数39)中の位相関数のコサイン平均をg=0(等方散乱)と仮定している。一般の生体関連物質では、g≒0.9であることが知られているが、上記の計算ではg=0の仮定にもかかわらず、得られた結果は、明確に予め設定された被検体における吸収物質濃度の空間分布を再現しており、吸収物質濃度が異なる5×5×5のサイズの空間が明確に識別されており、が空間分解能も良好である。
【0057】
図24に、X線CTのアルゴリズム( Filtered Back Projection 法)を用いた結果例を示す。図24は、外径25mm、高さ100mmを有する円柱(散乱係数1.0mm-1である)の中心軸位置に、外径5.mm、高さ100mmを有する円柱(空間分解能の評価において広く仮定される完全吸収体とし、吸収系数は無限大とする)の吸収体を配置した被検体に関する再構成画像である。画像の再構成に使用した計測データは、被検体に照射する光をδ関数と見做して、先に概要を説明したモンテカルロ法を使用したシミュレーションにより、図14に示すような検出光強度の時間スペクトルを生成した。画像の空間分解能を上げるために、時間ゲート幅を小さく設定して、検出光強度の時間スペクトルの0〔psec〕から300〔psec〕までの積分信号を計測信号とした。断層像は、計算により円柱の高さ方向の中心位置で Filtered Back Projection 法で求めた。
【0058】
図24では、再構成画像の結果を立体図で示しており、高さ方向の値の最大値を1.0、高さ方向の表示の最小間隔を0.05として示している。図18〜図23と図24との比較から、本発明による画像化方法は従来の方法と比較して明らかに空間分解能が向上していることがわかる。 Filtered Back Projection 法では、散乱の影響を十分考慮することができず、検出光強度の時間スペクトルの積分時間間隔を仮に短くした計測信号を使用しても、S/N比が劣化してしまい高い分解能を得ることは困難である。図24の結果では、外径5mmの吸収体は、空間分解能は半値幅で評価すると15mmであり、5mm×5mm×5mmの大きさの吸収体を再構成した、図18〜図23に示す結果では、吸収体は明確に鋭く(シャープ)に検出されている。このように、本発明での空間分解能は、 Filtered Back Projection 法に比較して約3倍以上改善されている。
【0059】
図25は、本発明の方法が適用される装置構成を説明する図である。光源9は少なくとも2波長の光を発する光源で、例えばチタンサファイアレーザ、もしくは複数の半導体レーザで構成される。波長の切り替えは、計算機24で制御する。まず、光源9から1波長(λ1)の光が選択され、光ファイバー、レンズ等の光学系で構成される光源用導波路10(以下、特に断りの無い場合、導波路は光ファイバーあるいは光学系で構成される)を通り、光は分波器11で参照光用導波路12と光源用光ファイバー13の2系統に分岐される。ストリークカメラ等の光検出器20により、参照光用導波路12に分岐した光は照射光強度I0をもつ参照光として計測される。光源用光ファイバー13に分岐した光は光スイッチ14に入射する。光スイッチ14は、被検体2の表面上に配置された複数の光照射用光ファイバー15の1本を選択し、光源用光ファイバー13と接続され、光源用光ファイバー13からの光は被検体2の任意位置から被検体2内部に照射される。
【0060】
ここで、光照射用光ファイバー15は、光ファイバー固定ホルダー16で位置を固定維持されている。被検体2を通過した光は、被検体2の表面上に配置された複数の光検出用光ファイバー17に入射する。光スイッチ18は、複数の光検出用光ファイバー17の中から必要な光ファイバーのみを選択し、検出器用光ファイバー19と接続する。選択された光検出用光ファイバー17を通過する検出光は光検出器20に導かれ、各検出器用光ファイバー19毎に検出光強度の時間分解信号が計測される。光源の強度信号と、検出光強度の時間分解信号は信号伝送ケーブル21でA/D変換器22に転送されデジタル信号に変換される。デジタル信号に変換された光源の強度信号と検出光強度の時間分解信号は、信号伝送ケーブル23で計算機24に転送され、計算機内の記憶装置に保存される。
【0061】
引続き、上記で用いた波長(λ1)と異なる波長(λ2)の光を選択し、上記と同様の計測を行なう。1照射位置に関して計測終了後、前回の照射位置あるいは検出位置と異なる位置における光計測をするため、光スイッチ14及び18を光スイッチ制御ケーブル26及び27を介して計算機24により制御して、光源用光ファイバー13と光照射用光ファイバー15の接続位置及び検出器用光ファイバー19と光検出用光ファイバー17の接続位置を切り替える。必要な計測が終了した後、先に説明した図15、図16、図17のフローチャートの手順に従い、吸収物質濃度の空間分布を求め表示装置25に表示する。
【0062】
図26に、光照射用光ファイバー15と光検出用光ファイバー17の配置例を示す。被検体の面の番号Piを図9と同様に定義し、照射位置と検出位置を図9と同様にする場合には、照射用光ファイバー13を(2,2)P1、(2,2)P2、(4,4)P3、(4,2)P4、(3,3)P5の5点に配置し、検出用光ファイバー17はP1〜P5の各面の各体積要素面の重心位置に配置する。ただし、前記の5点の光照射位置には、光照射用光ファイバー15と光検出用光ファイバー17の2本の光ファイバーが配置されるので、それぞれの光ファイバーの先端位置を多少ずらして配置するものとする。光照射用光ファイバー15、及び光検出用光ファイバー17には、例えば径1mmの石英ファイバーを用いる。
【0063】
図27に、被検体2を生体頭部とする場合の、光照射用光ファイバー15と光検出用光ファイバー17の配置例を示す。被検体が生体頭部の場合には、ヘルメット状の光ファイバー固定具28を用いて光照射用光ファイバー15、あるいは光検出用光ファイバー17を固定維持する。光ファイバー固定具28には、光照射用光ファイバー15、あるいは光検出用光ファイバー17を挿入する穴が設けられ、各光ファイバーを任意の穴位置に挿入する。被検体が生体の様な場合にも、仮想の体積要素5を立方体とするが、図27に示す様に、表面に接する体積要素5は立方体とはならない。
【0064】
ここで、被検体を1000個の体積要素に分割する場合には、計測位置(照射位置と検出位置の関係)を少なくとも1000組用意する必要があり、計測位置を1000組設定するために、図25に示す光照射用ファイバー15と光検出用光ファイバー17の配置方法としては、様々の方法がある。例えば、光照射用光ファイバー15を少なくとも20本、光検出用光ファイバー17を少なくとも1000本、それぞれ配置し、少なくとも20本の光照射用光ファイバー15から光を照射した時に、少なくとも1000本の光検出用光ファイバー17の中から少なくとも50本を選択して、検出光強度を計測すればよい。この計測は、図25と同様な装置構成で実現できる。
【0065】
生体に照射光として近赤外光を用いると、吸収物質としては酸化ヘモグロビン(以下HbO2と略す)と脱酸素化ヘモグロビン(以下Hbと略す)が支配的であり、これら2成分の吸収物質濃度の空間分布を求めることができる。この場合は、照射光として3種類の波長(λ1,λ2,λ3)の光を用いて、各波長毎に光計測を行い、(λ1,λ2)の波長を照射光として用いた計測から図15、図16、図17に記載の、step1からstep3までと同様のフローチャートの手順に従い、未知行列X、即ちχ(j,λ1,λ2)を求め、次に(λ1,λ3)の波長を照射光として用いた計測から同様にしてχ(j,λ1,λ3)を求め、得られたχ(j,λ1,λ2)、χ(j,λ1,λ3)から連立方程式(数37)を立て、連立方程式(数37)に各波長の光に対する吸収物質の吸光係数εを代入して解き、各体積要素毎のHbの濃度c(Hb,j)と、HbO2の濃度c(HbO2,j)を求める。
【0066】
図28に、連立方程式を解く際に用いるHbとHbO2の吸光係数の波長スペクトルを示す( S.Wray 等, " Characterization of the near infrared absorption spectra of cytochrome aa3 and hemoglobin for the non-invasive monitoring of cerebral oxygenation "; Biochim. Biophys. Acta 933, 184-192(1988) より数値を抜粋してグラフ化した)。図28中で、実線がHbの吸光係数の波長スペクトルを、点線がHbO2の吸光係数の波長スペクトルを、それぞれ示す。
【0067】
図28に示すような吸収物質に関する吸光係数の波長スペクトルをあらかじめ計算機内の記憶装置に記憶しておき、照射光として用いる波長に対する吸光係数を求める。頭部以外の生体組織、例えば乳房等の様な部位を被検体とする場合も、専用の光ファイバー固定具を使用し、照射光として近赤外光を用いて、被検体内部のHbとHbO2の空間分布を求めることができる。
【0068】
また、被検体中にγ種類の吸収物質が存在する場合には、照射光として(γ+1)種類の波長の光を用いて、被検体中のγ種類の吸収物質濃度の空間分布を得ることができる。例えば、近赤外光を用いて筋肉中のHbとHbO2とミオグロビンの空間分布を求める場合には、4種類の波長の照射光を用いればよい。
【0069】
被検体の一部の領域についての吸収物質濃度の空間分布を求める方法を、図29を用いて説明する。図29では、生体頭部を被検体2とし、被検体の一部である注目領域31内の吸収物質濃度の空間分布を求めることを示している。図29において、注目領域31の形状は25mm角の立方体で、5×5×5の合計125個の体積要素に分割し、各体積要素内の吸収物質濃度を求める。注目領域31の近傍に配置した光照射用光ファイバー15から照射光29を被検体2に照射し、被検体2内で散乱された光が、注目領域31の近傍に配置された光検出用光ファイバー17に到達し検出光30として検出される(この時検出される光は、光照射位置と同一方向に配置された検出器により検出される、反射光である)。但し、光照射用光ファイバー15から光検出用光ファイバー17までの光路3は、注目領域31外の領域を通過する光路もあるので、光が到達する領域32(以下では、光到達領域と略記する)の内部全体についての吸収物質濃度の空間分布を求めなければならない。光到達領域32の形状は、あらかじめモンテカルロ法によるシミュレーションや、拡散方程式を用いたシミュレーション等で決定する。
【0070】
例えば、光照射用光ファイバー15から光検出用光ファイバー17に到達する光の95%が通過する領域を包絡するように、光到達領域32を定める。
【0071】
光到達領域32内の吸収物質濃度の空間分布を求めるには、光到達領域32を注目領域31と同様に体積要素に分割して求める方法と、注目領域31外の光到達領域32に含まれる領域を1体積要素と扱って求める方法がある。両方法とも体積要素数が異なるだけであり、本質的には同じであるので、本実施例では注目領域31外の領域を1体積要素として説明する。従って、光吸収物質濃度の空間分布を求めるべき光到達領域32内には、体積要素数が126個存在することになり、異なる計測位置(照射位置と検出位置の位置関係)を少なくとも126組用意しなければならない。計測位置の配置方法は無数にあるが、図29において、例えば光照射用光ファイバー15を注目領域31付近に少なくとも9本配置し、光検出用光ファイバー17を注目領域付近に少なくとも14本配置すれば、少なくとも126(9×14)組の計測位置を用意できる。図25と同様な装置構成で、各計測位置毎に被検体2を通過してくる検出光30の強度を計測し、図15、図16、図17と同様のフローチャートの手順で、光到達領域32内の吸収物質濃度の空間分布を求めることができる。ただし、照射光として近赤外光を用い、各体積要素中のHbO2とHbの濃度を求める場合には、照射光として3種類の波長を用いる。また、平均光路長の空間分布行列Aを計算する際には、光到達領域32の外部形状が入力される。光到達領域32の外部形状は、あらかじめ求めることができるので、平均光路長の空間分布行列Aも、あらかじめ計算してメモリに記憶しておいても構わない。
【0072】
以下、被検体の任意の光照射位置から任意の光検出位置までの、複数の光路を計算して一旦メモリに光路データとして記憶しておき、記憶された光路データに基づいて、平均光路長の空間分布行列Aをもとめるための手順(以下では、簡単のため回転移動型モンテカルロ法と略記する)について図を用いて、詳細に説明する。任意の光照射位置から任意の光検出位置までの複数光路は、被検体内の光散乱を被検体の散乱特性を用いて模擬するモンテカルロ法で計算して求めた。従来の光散乱を模擬し光路を計算する従来のモンテカルロ法は、ウィルソン等の行った方法が一般的であり、先に説明した。被検体の任意の光照射位置から任意の光検出位置までの光路を、従来のモンテカルロ法を用いて、必要に応じてその都度計算してもよいが、その計算のために膨大な計算時間を要するため実用的ではない。本発明では、光路データを記憶しておき、これを使用することより、その都度光路を計算する従来方法より、約100倍以上高速化できる。
【0073】
以下では被検体形状として、図12あるいは図13と同様な立方体を用いて説明する。図30は、回転移動型モンテカルロ法を説明するために想定した被検体形状と、照射光29を照射する光照射位置33と検出光30を検出する光検出位置34の相対位置関係を示す図である。図30において光照射位置33から光検出位置34に至るまでの複数の光路を回転移動型モンテカルロ法で求める手順を、代表例として以下説明する。図30において、立方体の被検体2は、仮想的に5×5×5の体積要素に分割され、各体積要素5は立方体である。回転移動型モンテカルロ法の概念を図31に示す。図31は、簡単のため、光照射位置33と光検出位置34の位置を含む断面を示しており、この断面内に存在する光路37が示される。最初に、半径が被検体2の任意の2点をとりその間の距離の最大値(図30の被検体形状では、頂点35−1と頂点35−2の距離)と等しい、球形の仮想散乱体36を計算機内で発生する。球形の仮想散乱体36(図30では図示せず)の中心Oから発射される1光子は、被検体2と等価な散乱特性で多重散乱を繰り返し、仮想散乱体36の境界面まで到達する。ここで与えている散乱特性とは、先に説明した被検体の散乱係数と位相関数で与えられる。
【0074】
仮想散乱体36の境界面に到達するまでの光子の散乱点を、光路データとして記憶保存しておくことにより、1光子の仮想散乱体36中における光路37を再現することができる。発射された1光子が境界面に到達した後、同様の過程を繰り返して複数の光子に関する光路を得る。統計的な精度を実質的に補償するに十分な光路数が得られた後、図31に示す様に、被検体2の外形形状を有する仮想被検体38を発生する。被検体2上の光照射位置33と光検出位置34の位置に対応して、仮想被検体38上に光照射位置33と光検出位置34を設定する。続いて、仮想被検体38上の光照射位置33を、球形の仮想散乱体36の中心位置に配置し、球形の仮想散乱体36の中心Oを回転中心として仮想被検体38を回転させて、仮想被検体38と光路の相対位置を変化させ、仮想被検体38上の光検出位置34と、上記の過程で記憶保存されている光路とが交差する点を求め、中心Oからこの交差点に至るまでの光路を再現する散乱点の座標を仮想被検体38に設定されている座標系に変換する。即ち、球形の仮想散乱体36の中心Oを回転中心として、仮想被検体38上の光照射位置33と光検出位置34を結ぶ線分を回転させて、この線分の先端点である光検出位置34と記憶保存されている光路とが交差する点を求める。上記の保存されている光路データのそれぞれに関して、同様の操作を繰り返し、仮想被検体38上の光照射位置33から光検出位置34に至るまでの複数の光路が求められる。例えば、図31中に示す様に、1光路に関して処理が終了した後、仮想散乱体38を中心Oの周りに矢印の方向に回転し、次の光路データに関する処理を行う。求められた仮想被検体36上の光照射位置33から光検出位置34に至るまでの光路データ(各光路に属する散乱点のデータ)から、各体積要素毎の光路長を計算することができ、各体積要素毎の光路長を記憶しておき、被検体2上の光照射位置33と光検出位置34の計測位置における、平均光路長の空間分布を(数42)から計算することができる。以上が、回転移動型モンテカルロ法の2次元での説明であるが、3次元においても同様の処理ができることは言うまでもない。被検体の形状が一般の場合にも、被検体を複数の空間要素に分割して、被検体の外形形状を有する仮想被検体を処理装置において生成して、以上の説明と同様の処理を行なって、任意形状の被検体に関しての、光照射位置と光検出位置とが設定されたとき、この光照射位置と光検出位置とを通過する平均光路長の空間分布を(数42)から計算することができる。
【0075】
以上説明した処理の手順をフローチャートとして図32、及び図33に示し、以下にその詳細を各ステップ毎に説明する。ただし、ここでは図17に示されている様に、既に被検体の外形形状を表すデータが入力されており、さらに、被検体の散乱特性も入力されているものとする。
【0076】
step1: 被検体の表面上に任意の2点を取り、その距離の最大値が半径となるような球の仮想散乱体36を発生する。仮想散乱体36の中心から光子を出射し、仮想散乱体境界面から光子が出ていくまでの光路データを保存する。統計精度を実質的に補償するために、少なくとも1000本の光路データを保存する。仮想散乱体36の散乱特性は、被検体の散乱特性を反映する。
【0077】
step2: 体積要素に分割された仮想被検体を発生し、被検体上の光照射位置に対応して、光照射位置が仮想被検体に設定される。仮想被検体上の光照射位置を球形の仮想散乱体中心に配置する。
【0078】
step3: step1で保存されている任意の光路データから、光路の1本を選択し、被検体上の光照射位置と同位置に設定される仮想被検体上の光検出位置と、光路が交差するように被検体を回転する。この時、回転中心は仮想散乱体の中心(即ち光照射位置)とする。本ステップで任意光路を選ぶ際に、既に選択された光路は重複して選択しない。
【0079】
step4: step3で光照射位置と光検出位置を通過した光路が、光路中で仮想被検体の境界面と交差するか否かを判断する。交差する場合はstep3に戻り、次の光路に関して処理を行ない、交差しない場合はstep5に進む。
【0080】
step5: 光照射位置から光検出位置への到達時間tを総光路長から計算し、検出時間条件t1、t2に対して、t1≦t≦t2の条件を満たしているか否かを判断する。条件を満たしていない場合は、step3に戻り次の光路データに関して処理を行ない、条件を満たしている場合は、step6に進む。
step6: 光照射位置から光検出位置に至るまでの光路の散乱点の座標を、仮想散乱体の座標系から仮想被検体の座標系に変換した後、各体積要素毎の光路長p(i,j,k)を計算し、記憶保存する。
step7: 保存されている全光路データに関する以上の各stepの処理が終了したか否かを判断する。計算終了でない場合はstep3に戻り次の光路に関して処理を行ない、処理終了の場合はstep9へ進む。
【0081】
step8: 本ステップまでに処理された光路の累積数mは、(数38)の右辺の分母に相当するmであり、検出器への到達時間〔t1〜t2〕の時間区間に到達する光子数mと等しく、step6で保存されている各体積要素毎の光路長p(i,j,k)から、各体積要素毎の累積光路長((数38)の右辺の分子)を計算し、光路の累積数mで除算する。本stepでi番目の光計測位置(光照射位置と光検出位置の組)における、j番目の体積要素毎の平均光路長LAij(平均光路長の空間分布)を求めたことになる。異なる計測位置に関して平均光路長の空間分布を計算する場合には、step9に進む。
【0082】
step9: 全計測位置(光照射位置と光検出位置の組)に関して、処理を終了したか否かを判断する。処理終了していない場合は、step2に戻り次の計測位置に関して処理を行なう。処理終了の場合には、step10に進む。
【0083】
step10: 全計測位置に関する平均光路長の空間分布を記憶保存する。
【0084】
以上のフローチャートに示す手順では、全計測位置に関して必ずしも処理する必要はない。幾何学的に等価な関係にある異なる計測位置(光照射位置と光検出位置の組)が複数ある場合に関しては、その中から1つの計測位置を選択して、光路を求め、得られた光路が、他の幾何学的に等価な計測位置を通過するように座標変換して、処理時間を短縮できる。この処理は、前記のstep6で行えばよい。例えば、図30において、光照射位置33が(3,3)P1であり、光検出位置34が(4,3)P2である計測位置は、光照射位置が(3,3)P1であり、光検出位置が面P2と平行な面P4上の点である(4,3)P4、等の計測位置と幾何学的に等価と考える(面Piの定義は図9参照)。
【0085】
以上が、回転移動型モンテカルロ法の説明であるが、同様の考え方により、図34に示すように、仮想散乱体36(十分大きいサイズを有する空間であれば任意の形状でよい)の1点(O0)から、発する1光子の複数の散乱点を、光路データとして記憶保存しておき、光路37を再現することができる(この方法を、移動回転型モンテカルロ法と呼ぶ)。なお、図34では、図31と同様に、簡単のため、図30に示される光照射位置33と光検出位置34の位置を含む断面を示しており、この断面内に存在する光路37が示される。仮想散乱体36中における仮想被検体38上の光照射位置33を、仮想散乱体36の任意の位置Oに配置し、仮想散乱体36の中心Oを回転中心として仮想被検体38を回転させて、仮想被検体38と光路37の相対位置を変化させ、仮想被検体38上の光検出位置34と、保存された光路37とが交差する点を求め、点Oからこの交差点に至るまでの光路を再現する散乱点の座標を仮想被検体38に設定されている座標系に変換し、各体積要素毎の光路長p(i,j,k)を計算し、記憶保存する。次いで、光照射位置33の位置を点Oから、光路の沿って点O’に移動させて同様の処理を行なう。点Oと点O’との間の距離は、仮想被検体の光照射位置33と光検出位置34との間の距離より大きくすることが望ましい。その他の処理に関しては、先に説明した回転移動型モンテカルロ法の場合と同様である。以上が、移動回転型モンテカルロ法の2次元での説明であるが、3次元においても同様の処理ができることは言うまでもない。同様にして、被検体の形状が一般の場合にも、被検体に光照射位置と光検出位置とが設定されたとき、この光照射位置と光検出位置とを通過する平均光路長の空間分布を(数38)から計算することができる。
【0086】
【発明の効果】
従来の画像化方法では、空間分解能、定量性、計算時間の点で問題点があったが、本発明の方法では、従来方法に比較して、空間分解能で約3倍以上の改善がなされた3次元的な吸収物質の空間分布を得ることができ、10〜100倍程度計算時間が短縮され、局所の最適解に落ち込むことなく解を得ることができ、さらに定量性の改善が可能である。また、回転移動型、もしくは移動回転型モンテカルロ法を用いることにより、従来のモンテカルロ法を適用する方法と比較して、平均光路長の空間分布を約100倍以上の速度で求めることができる。
【図面の簡単な説明】
【図1】吸収物質を一様に含む被検体中での光路を含む断面図。
【図2】異なる層状の吸収物質からなる被検体中での光路を含む断面図。
【図3】散乱体に吸収物質を含む被検体中での光路を含む断面図。
【図4】散乱体に吸収物質を含む被検体中を複数体積要素に分割したときの各体積要素での光路を模式的に示す図。
【図5】被検体に光源と検出器を配置する例を示す図。
【図6】検出器に到達する光子の光路例を示す図。
【図7】照射光の時間スペクトルの例を示す図。
【図8】被検体を通過し検出される光の時間スペクトルの例を示す図。
【図9】実施例に用いた被検体の形状を示す図。
【図10】実施例に用いた被検体を体積要素に分割する例を示す図。
【図11】実施例に用いた光照射位置の配置例を示す図。
【図12】実施例に用いた被検体内での吸収体の位置を示す図。
【図13】実施例に用いた被検体内での吸収体の位置を示す図。
【図14】実施例に用いた計算により求めた検出光の時間スペクトル例を示す図。
【図15】被検体内の吸収物質濃度の空間分布を求める処理を示すフローチャート。
【図16】被検体内の吸収物質濃度の空間分布を求める処理を示すフローチャート。
【図17】被検体内の吸収物質濃度の空間分布を求める処理を示すフローチャート。
【図18】本発明により得られた吸収物質濃度の空間分布の例を示す図。
【図19】本発明により得られた吸収物質濃度の空間分布の例を示す図。
【図20】本発明により得らえた吸収物質濃度の空間分布の例を示す図。
【図21】本発明により得らえた吸収物質濃度の空間分布の例を示す図。
【図22】本発明により得られた吸収物質濃度の空間分布の例を示す図。
【図23】本発明により得られた吸収物質濃度の空間分布の例を示す図。
【図24】従来方法により得られた吸収物質濃度の空間分布の例を示す図。
【図25】本発明が適用される装置の構成例を示す図。
【図26】被検体に光照射用、及び光検出用光ファイバーを配置する例を示す図。
【図27】生体頭部に光照射用光、及び光検出用ファイバーを配置する例を示す図。
【図28】酸素化ヘモグロビンと脱酸素化ヘモグロビンの吸光係数の波長スペクトル。
【図29】本発明を生体の1部の領域に適用する例を示す図。
【図30】本発明による回転移動型モンテカルロ法を適用する被検体の例を示す図。
【図31】本発明による回転移動型モンテカルロ法の概念を説明するための図。
【図32】本発明による回転移動型モンテカルロ法により平均光路長の空間分布を求める手順を示すフローチャート。
【図33】本発明による回転移動型モンテカルロ法により平均光路長の空間分布を求める手順を示すフローチャート。
【図34】本発明による移動回転型モンテカルロ法の概念を説明するための図。
【符号の説明】
P1…被検体の第1面、P2…被検体の第2面、P3…被検体の第3面、P4…被検体の第4面、P5…被検体の第5面、P6…被検体の第6面、1…光源、2…被検体、3…光路、4…検出器、5…体積要素、6…光照射位置、7…照射光、8…吸収物質、9…光源、10…光源用導波路、11…分波器、12…参照光用導波路、13…光源用光ファイバー、14…光スイッチ、15…光照射用光ファイバー、16…光ファイバー固定具、17…光検出用光ファイバー、18…光スイッチ、19…検出器用光ファイバー、20…光検出器、21…信号伝送ケーブル、22…A/D変換器、23…信号伝送ケーブル、24…計算機、25…表示装置、26…光スイッチ制御ケーブル、27…光スイッチ制御ケーブル、28…光ファイバー固定具、29…照射光、30…検出光、31…注目領域、32…注目領域外の領域、33…光照射位置、34…光検出位置、35…被検体の頂点、36…仮想散乱体、37…光路、38…仮想被検体。
Claims (4)
- 被検体に光を照射し、前記被検体内の光吸収物質の濃度の空間分布を画像化する方法であって、
前記被検体上の光照射位置と光検出位置の複数の対に対してそれぞれ、前記光照射位置から前記光検出位置に至るまでの散乱光の光路をシミュレーションにより求めるステップと、
前記光路から前記複数の対ごとに平均光路長を算出するステップと、
複数の前記平均光路長の空間分布を示す空間分布行列を作るステップと、
前記空間分布行列と、前記光吸収物質の光学定数と、前記光検出位置で検出された光の強度とから、前記光吸収物質の濃度の空間分布を算出するステップと、
前記光吸収物質の濃度の空間分布を表示するステップとを有することを特徴とする吸収物質濃度空間分布画像化方法。 - 前記光路は、モンテカルロ法により求められることを特徴とする請求項1記載の吸収物質濃度空間分布画像化方法。
- 複数の前記光検出位置で検出された光の強度から計測値行列を作るステップを有し、
前記空間分布行列と、前記光吸収物質の光学定数と、前記計測値行列とから前記光吸収物質の濃度の空間分布を算出することを特徴とする請求項1記載の吸収物質濃度空間分布画像化方法。 - 前記空間分布行列を記憶手段に保存するステップを有することを特徴とする請求項1記載の吸収物質濃度空間分布画像化方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2002078507A JP3730927B2 (ja) | 2002-03-20 | 2002-03-20 | 吸収物質濃度空間分布画像化装置 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2002078507A JP3730927B2 (ja) | 2002-03-20 | 2002-03-20 | 吸収物質濃度空間分布画像化装置 |
Related Parent Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP16182394A Division JP3310782B2 (ja) | 1994-07-14 | 1994-07-14 | 吸収物質濃度の空間分布画像化装置 |
Publications (2)
Publication Number | Publication Date |
---|---|
JP2002333400A JP2002333400A (ja) | 2002-11-22 |
JP3730927B2 true JP3730927B2 (ja) | 2006-01-05 |
Family
ID=19193320
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2002078507A Expired - Lifetime JP3730927B2 (ja) | 2002-03-20 | 2002-03-20 | 吸収物質濃度空間分布画像化装置 |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP3730927B2 (ja) |
Families Citing this family (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6944322B2 (en) * | 2001-03-28 | 2005-09-13 | Visiongate, Inc. | Optical tomography of small objects using parallel ray illumination and post-specimen optical magnification |
CN100432699C (zh) * | 2006-12-29 | 2008-11-12 | 成都川大奇林科技有限责任公司 | 一种测量医用加速器光子束能谱的方法 |
JP5465924B2 (ja) * | 2009-05-15 | 2014-04-09 | 宜彦 水島 | 不均一物質の光分析方法 |
JP6212312B2 (ja) | 2012-08-13 | 2017-10-11 | パナソニック株式会社 | 物体内部推定装置およびその方法 |
JP6154613B2 (ja) * | 2013-01-18 | 2017-06-28 | 浜松ホトニクス株式会社 | 断面画像計測装置及び計測方法 |
-
2002
- 2002-03-20 JP JP2002078507A patent/JP3730927B2/ja not_active Expired - Lifetime
Also Published As
Publication number | Publication date |
---|---|
JP2002333400A (ja) | 2002-11-22 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP3310782B2 (ja) | 吸収物質濃度の空間分布画像化装置 | |
US6956650B2 (en) | System and method for enabling simultaneous calibration and imaging of a medium | |
Wilson et al. | Time-dependent optical spectroscopy and imaging for biomedical applications | |
US5931789A (en) | Time-resolved diffusion tomographic 2D and 3D imaging in highly scattering turbid media | |
US6108576A (en) | Time-resolved diffusion tomographic 2D and 3D imaging in highly scattering turbid media | |
US5813988A (en) | Time-resolved diffusion tomographic imaging in highly scattering turbid media | |
JP3112025B2 (ja) | 生体計測装置 | |
JPH06221913A (ja) | 散乱吸収体の光学情報計測装置及び方法 | |
US5758653A (en) | Simultaneous absorption and diffusion imaging system and method using direct reconstruction of scattered radiation | |
Murad et al. | Reconstruction and localization of tumors in breast optical imaging via convolution neural network based on batch normalization layers | |
JP3730927B2 (ja) | 吸収物質濃度空間分布画像化装置 | |
JPH03505922A (ja) | ランダムな媒体を画像化するためのシステム | |
JP3164112B2 (ja) | 生体計測装置 | |
Rasmussen et al. | Radiative transport in fluorescence‐enhanced frequency domain photon migration | |
JP3276947B2 (ja) | 吸収物質濃度空間分布画像装置 | |
US5832922A (en) | Diffusion Tomography system and method using direct reconstruction of scattered radiation | |
US5746211A (en) | Absorption imaging system and method using direct reconstruction of scattered radiation | |
US5787888A (en) | Absorption tomography system and method using direct reconstruction of scattered radiation | |
Vidal-Rosas | Advanced tomographic image reconstruction algorithms for Diffuse Optical Imaging | |
Hashmi et al. | Optical tomography in medical imaging and diagnostic engineering | |
JP2005331292A (ja) | 画面再構成方法 | |
Chance et al. | Photon dynamics in tissue imaging | |
Hideo | 3-7 Measurement of Brain Activity by Near Infrared Light | |
Böcklin et al. | Determining light distribution in human head using 3D Monte Carlo simulations | |
Haller | Time resolved breast transillumination: analytical, numerical and experimental study |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20040608 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20040722 |
|
A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20050405 |
|
A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20050511 |
|
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: 20050927 |
|
A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20051007 |
|
R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20091014 Year of fee payment: 4 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20091014 Year of fee payment: 4 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20101014 Year of fee payment: 5 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20111014 Year of fee payment: 6 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121014 Year of fee payment: 7 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20121014 Year of fee payment: 7 |
|
FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20131014 Year of fee payment: 8 |
|
EXPY | Cancellation because of completion of term |