(57)【特許請求の範囲】
【請求項１】
被検体の複数の部位に個別排他的に投射光を投射する光源手段と、
前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、
前記被検体を微少な複数個の体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第１の演算手段と、
前記第１の演算手段で想定した前記集合体モデルの特定の体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第２の演算手段と、
前記第１の演算手段で求められた前記光密度と前記第２の演算手段で求められた光密度との比を前記体積要素の影響度として演算して、前記複数の体積要素の配列に対応する影響度の行列を求める第３の演算手段と、
前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第１の演算手段で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める第４の演算手段と、
前記第３の演算手段により求められた影響度の行列と前記第４の演算手段により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める第５の演算手段と、
を具備することを特徴とする光ＣＴ装置。
【請求項２】
被検体の複数の部位に個別排他的に投射光を投射する光源手段と、
前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、
前記被検体を複数個の要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第１の演算手段と、
前記第１の演算手段で想定した前記集合体モデルの特定の要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第２の演算手段と、
前記第１の演算手段で求められた光密度と前記第２の演算手段で求められた光密度との比を前記要素の影響度として演算して、前記複数の要素の配列に対応する影響度の行列を求める第３の演算手段と、
前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第１の演算手段で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める第４の演算手段と、
前記第３の演算手段により求められた影響度の行列と前記第４の演算手段により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める第５の演算手段と、
を具備することを特徴とする光ＣＴ装置。
【請求項３】
被検体の複数の部位に個別排他的に投射光を投射する光源手段と、
前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、
前記被検体を複数の有限体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第１の演算手段と、
前記第１の演算手段で想定した前記集合体モデルの特定の有限体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第２の演算手段と、
前記第１の演算手段で求められた光密度と前記第２の演算手段で求められた光密度との比を前記有限体積要素の影響度として演算して、前記複数の有限体積要素の配列に対応する影響度の行列を求める第３の演算手段と、
前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と、前記第１の演算手段で求められた光強度との相対比を演算し、前記複数の有限体積要素の配列に対応する相対比の行列を求める第４の演算手段と、
前記第３の演算手段により求められた影響度の行列と前記第４の演算手段により求められた相対比の行列とのマトリクス演算により、前記各有限体積要素の特徴データを求める第５の演算手段と、
を具備することを特徴とする光ＣＴ装置。
【請求項４】
前記光源手段は、連続光を前記投射光として出射することを特徴とする請求項１乃至請求項３のいずれか１項に記載の光ＣＴ装置。
【請求項５】
前記光源手段は、パルス光を前記投射光として出射し、前記光検出手段は、前記透過光を時間分解計測することを特徴とする請求項１乃至請求項３のいずれか１項にに記載の光ＣＴ装置。
【請求項６】
光源手段により、被検体の複数の部位に個別排他的に投射光を投射すると共に、光検出手段により、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する第１の行程と、
前記被検体を微少な複数個の体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第２の行程と、
前記第２の行程で想定した前記集合体モデルの特定の体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第３の行程と、
前記第２の行程で求められた前記光密度と前記第３の行程で求められた光密度との比を前記体積要素の影響度として演算して、前記複数の体積要素の配列に対応する影響度の行列を求める演算をする第４の行程と、
前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第１の演算手段で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める演算をする第５の行程と、
前記第４の行程により求められた影響度の行列と前記第５の行程により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める演算をする第６の行程と、
を備えることを特徴とする光ＣＴによる画像再構成方法。
【請求項７】
光源手段により、被検体の複数の部位に個別排他的に投射光を投射すると共に、光検出手段により、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する第１の行程と、
前記被検体を複数個の要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第２の行程と、
前記第２の行程で想定した前記集合体モデルの特定の要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第３の行程と、
前記第２の行程で求められた光密度と前記第３の行程で求められた光密度との比を前記要素の影響度として演算して、前記複数の要素の配列に対応する影響度の行列を求める第４の行程と、
前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第２の行程で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める第５の行程と、
前記第４の行程により求められた影響度の行列と前記第５の行程により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める第６の行程と、
を備えることを特徴とする光ＣＴによる画像再構成方法。
【請求項８】
光源手段により、被検体の複数の部位に個別排他的に投射光を投射すると共に、光検出手段により、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する第１の行程と、
前記被検体を複数の有限体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第２の行程と、
前記第２の行程で想定した前記集合体モデルの特定の有限体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第３の行程と、
前記第２の行程で求められた光密度と前記第３の行程で求められた光密度との比を前記有限体積要素の影響度として演算して、前記複数の有限体積要素の配列に対応する影響度の行列を求める第４の行程と、
前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と、前記第２の演算手段で求められた光強度との相対比を演算し、前記複数の有限体積要素の配列に対応する相対比の行列を求める第５の行程と、
前記第４の行程により求められた影響度の行列と前記第５の行程により求められた相対比の行列とのマトリクス演算により、前記各有限体積要素の特徴データを求める第６の行程と、
を備えることを特徴とする光ＣＴによる画像再構成方法。
【請求項９】
前記光源手段により、連続光を前記投射光として出射することを特徴とする請求項６乃至請求項８のいずれか１項に記載の光ＣＴによる画像再構成方法。
【請求項１０】
前記光源手段により、パルス光を前記投射光として出射し、前記光検出手段により、前記透過光を時間分解計測することを特徴とする請求項６乃至請求項８のいずれか１項に記載の光ＣＴによる画像再構成方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、生体等の被検体に光を投射し、その透過光の測定値から被検体の断層画像を再構成する光ＣＴ装置及び光ＣＴによる画像再構成方法に関する。
【０００２】
【従来の技術】
従来の光ＣＴ装置は、生体等の被検体に光を投射すると共に被検体内部を透過して来た光を検出器で測定し、その測定結果に基づいて、その光（光束）が生体内の特定領域を通過する際の通過頻度（平均光路長の空間分布）をモンテカルロシミュレーションにより求めて、断層画像を再構成するアルゴリズムを採用していた。つまり、被検体を複数の体積要素（ＶＯＸＥＬ）に分割したモデルを想定し、各体積要素ごとに光束の衝突回数又は平均光路長（各体積要素を横切る飛行距離の総和または平均値）を計算することにより、画像再構成に必要な各体積要素の特徴抽出を行うこととしていた。
【０００３】
【発明が解決しようとする課題】
しかし、従来の光ＣＴ装置にあっては、モンテカルロシミュレーションは、光束の散乱粒子での散乱事象を、乱数を用いて計算する光線追跡法であるため、実際の測定では検出器に到達しない光束があったとしても、想定され得る全ての光束について計算することとなる結果、無駄な情報（事象）についてまでも計算対象とするため、演算に要する時間が極めて長くなるという問題があった。
【０００４】
特に、高品質の断層画像を再構成しようとすれば、当然に計算対象が増加し、無駄な情報に関する演算所要時間の占める割合が激増することになるので、実用に供し得る光ＣＴ装置を実現することが困難となっていた。
【０００５】
本発明は、このような従来の技術に鑑みてなされたものであり、被検体の断層画像を高速に再構成することができる光ＣＴ装置を提供することを目的とする。
【０００６】
【課題を解決するための手段】
このような目的を達成するために本発明は、被検体の複数の部位に個別排他的に投射光を投射する光源手段と、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、前記被検体を微少な複数個の体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第１の演算手段と、前記第１の演算手段で想定した前記集合体モデルの特定の体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第２の演算手段と、前記第１の演算手段で求められた前記光密度と前記第２の演算手段で求められた光密度との比を前記体積要素の影響度として演算して、前記複数の体積要素の配列に対応する影響度の行列を求める第３の演算手段と、前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第１の演算手段で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める第４の演算手段と、前記第３の演算手段により求められた影響度の行列と前記第４の演算手段により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める第５の演算手段と、を具備することとした。
【０００７】
また、被検体の複数の部位に個別排他的に投射光を投射する光源手段と、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、前記被検体を複数個の要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第１の演算手段と、前記第１の演算手段で想定した前記集合体モデルの特定の要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第２の演算手段と、前記第１の演算手段で求められた光密度と前記第２の演算手段で求められた光密度との比を前記要素の影響度として演算して、前記複数の要素の配列に対応する影響度の行列を求める第３の演算手段と、前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第１の演算手段で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める第４の演算手段と、前記第３の演算手段により求められた影響度の行列と前記第４の演算手段により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める第５の演算手段と、を具備することとした。
【０００８】
また、被検体の複数の部位に個別排他的に投射光を投射する光源手段と、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、前記被検体を複数の有限体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第１の演算手段と、前記第１の演算手段で想定した前記集合体モデルの特定の有限体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第２の演算手段と、前記第１の演算手段で求められた光密度と前記第２の演算手段で求められた光密度との比を前記有限体積要素の影響度として演算して、前記複数の有限体積要素の配列に対応する影響度の行列を求める第３の演算手段と、前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と、前記第１の演算手段で求められた光強度との相対比を演算し、前記複数の有限体積要素の配列に対応する相対比の行列を求める第４の演算手段と、前記第３の演算手段により求められた影響度の行列と前記第４の演算手段により求められた相対比の行列とのマトリクス演算により、前記各有限体積要素の特徴データを求める第５の演算手段と、を具備することとした。
【０００９】
また、前記光源手段は、連続光を前記投射光として出射する構成とした。
【００１０】
また、前記光源手段は、パルス光を前記投射光として出射し、前記光検出手段は、前記透過光を時間分解計測する構成とした。
【００１１】
【発明の実施の形態】
（第１の実施の形態）
本発明の第１の実施の形態を図面と共に説明する。
【００１２】
図１は、本実施の形態に係る光ＣＴ装置の構成を示す。同図において、この光ＣＴ装置は、被検体Ｍの光学特性を実測するための測定機構と、測定機構により得られた実測データを演算処理することにより被検体Ｍの断層画像を再構成する演算機構を備えている。
【００１３】
演算機構は、本装置全体の動作を制御するマイクロプロセッサユニット（ＭＰＵ）を有する制御部２と、予め決められたアルゴリズムに基づいて作成されたプログラムを実行することにより被検体Ｍの断層画像を再構成するための演算処理を行う演算部４と、演算部４が演算処理の際に種々のデータを格納等するための記憶部６と、上記演算処理により得られた被検体Ｍの断層画像等を表示するための表示部８とを具備している。
【００１４】
測定機構は、被検体Ｍに所定波長の光を投射するための光投射部と、被検体Ｍの透過光を実測するための受光部とを備えている。
【００１５】
即ち、上記光投射部にあっては、所定波長の光を出射する光源１０の光出射端に、１入力ｎ出力の光分岐路を有する光スイッチＳＷinの入力端が接続されると共に、光スイッチＳＷinの各出力端にｎ本の光ファイバＦ１〜Ｆｎが接続されている。更に、ｎ本の光ファイバＦ１〜Ｆｎは、これらの光出射端が面一に揃えられて同一方向に沿って一列に配列され、機械的に安定ないわゆる光投射プローブとなっている。そして、光ファイバＦ１〜Ｆｎの光出射端を被検体Ｍの一側端に対向又は接触させて使用される。
【００１６】
上記光スイッチＳＷinは、制御部２からの切換制御に従って、ｎ個中の１つの分岐路のみを選択的にオン、残余の分岐路をオフにする機能を有することにより、光源１０からの光を個別排他的に各光ファイバＦ１〜Ｆｎに導入し、その光の導入された光ファイバの光出射端のみから被検体Ｍに光を投射する。従って、被検体Ｍの位置の異なる部位へ順次に光を投射する等の走査が可能となっている。
【００１７】
一方、上記受光部にあっては、ｎ本の光ファイバＢ１〜Ｂｎと、これらの光ファイバＢ１〜Ｂｎの光出射端に接続されたｎ入力１出力の光分岐路を有する光スイッチＳＷout と、光スイッチＳＷout の光出力端に接続された光電変換デバイスなどを有する光検出部１２を具備している。そして、ｎ本の光ファイバＢ１〜Ｂｎは、夫々の光入射端が面一に揃えられて同一方向に沿って一列に配列され且つ機械的に安定ないわゆる受光プローブとなっており、これらの光ファイバＢ１〜Ｂｎの光入射端を被検体Ｍの一側端に対向又は接触させて使用される。
【００１８】
上記光スイッチＳＷout は、制御部２からの切換制御に従って、ｎ個中の１つの分岐路のみを選択的にオン、残余の分岐路をオフにする機能を有することにより、各光ファイバＢ１〜Ｂｎからの光を個別排他的に光検出部１２中の上記光電変換デバイスへ導出する。したがって、被検体Ｍの位置の異なる部位からの光を順次に検出する等の検出走査が可能となっている。
【００１９】
そして、被検体Ｍを実測する際には、光ファイバＦ１〜Ｆｎによる光投射プローブと光ファイバＢ１〜Ｂｎによる受光プローブとを被検体Ｍを介して対向配置し、光源１０からの光を各光ファイバＦ１〜Ｆｎにより順次に被検体Ｍに投射し、被検体Ｍの内部を透過してきた光（光束）を各光ファイバＢ１〜Ｂｎで順次に導入して光検出部１２で光電変換及びＡ／Ｄ変換することにより、被検体Ｍの各部位の光学特性のデータを得るようになっている。
【００２０】
また、図１に示す如く、光ファイバＦ１〜Ｆｎ（光投射プローブ）のｎ個の光出射端と光ファイバＢ１〜Ｂｎ（受光プローブ）のｎ個の光入射端によって画成される領域が、画像再構成されるべき被検体Ｍの一つの断層面に対応する。よって、光ファイバＦ１〜Ｆｎ（光投射プローブ）と光ファイバＢ１〜Ｂｎ（受光プローブ）を、被検体Ｍに対して、図１中の縦方向（Ｚ座標軸方向）に相対移動させることにより、複数の断層面を画成することができ、これらの複数の断層面についての実測データに基づいて後述の画像処理を行うことにより、被検体Ｍ内部の三次元断層画像を再構成する。
【００２１】
尚、説明の都合上、被検体ＭがＸＹＺ直交座標系中に配置され、光ファイバＦ１〜Ｆｎ及び光ファイバＢ１〜ＢｎはＹ座標軸方向に配列され、光ファイバＦ１〜Ｆｎの光出射端と光ファイバＢ１〜Ｂｎが光入射端がＸ座標軸方向に沿って対向配置されるものとする。更に、装置構成として、光スイッチの代わりに多数個の光源をオン・オフしたり、多数個の検出器を同時測定するようにしてもよい。
【００２２】
次に、この光ＣＴ装置の動作及び被検体Ｍの断層画像の再構成原理を、実際の処理過程に対応させて説明する。尚、図２は、再構成原理を説明するための説明図、図３は、動作を処理過程に沿って示すフローチャート、図４は、本実施の形態の再構成原理の妥当性を示す実験結果を示す。
【００２３】
（第１の処理）
図３のステップ１００に示す如く、光ファイバＦ１〜Ｆｎ又は光ファイバＢ１〜Ｂｎの各ピッチ間隔で、被検体ＭをＸＹＺ座標上でｎ×ｎ×ｎ個の体積要素（ＶＯＸＥＬ）に分割するものとし、かかる設定条件は制御部２及び演算部４に予め格納される。例えば、ｎ＝１００として、被検体Ｍを１００3 個の体積要素の集合体として扱う。
【００２４】
更に、演算部４には、次に述べるアルゴリズムに基づく演算プログラムが予め格納されている。
【００２５】
各体積要素に入射する連続光（光束）の定常光拡散方程式は、
【００２６】
【数１】
0003771339.tif 000002
で表されるので、この定常拡散方程式(1) を差分形式で表した次式(2) が演算部４に予め格納されている。
【００２７】
【数２】
0003771339.tif 000003
【００２８】
ここで、添字ｉ，ｊ，ｋは、ＸＹＺ直交座標系の節点の座標を表し、φi,j,k は、ＸＹＺ座標系の位置（i ，j ，k ）に存在する体積要素を通る光束を表し、１≦ｉ≦ｎ、１≦ｊ≦ｎ、１≦ｋ≦ｎである。
【００２９】
また、被検体Ｍの内部と外部との境界条件は、
【００３０】
【数３】
0003771339.tif 000004
としている。尚、添字ＢＬは、被検体Ｍ周辺の境界を表すものとする。また、式(3) は、被検体Ｍの表面で光が完全に吸収されてしまう様な条件、例えば、被検体Ｍの表面が真黒に塗られた状態と等価である。
【００３１】
これらの式(2)(3)を用いて、被検体Ｍの断層画像を再構成するのであるが、説明を明確にするために、図２(a) に示す如く、一つの断層面に注目して述べることとする。即ち、ＸＹ座標面上の一つの断層面を、ｎ本ずつの光ファイバＦ１〜ＦｎとＢ１〜Ｂｎに対応させて、ｎ×ｎ個の体積要素（ＶＯＸＥＬ）に分割するものとする。そして、光ファイバＦ１〜Ｆｎの各光出射端より投射される各投射光をＨ１〜Ｈｎとすると共に、一般表現のＨs （１≦ｓ≦ｎ）で総称するものとし、光ファイバＢ１〜Ｂｎの各光入射端で検出される検出光をＱ１〜Ｑｎとすると共に、一般表現のＱt（１≦ｔ≦ｎ）で総称するものとし、更に、各体積要素のＸＹ座標上の位置を（i ，j ）で表すものとする。
【００３２】
（第２の処理）
次に、前記式(2)(3)を用いて第１のシミュレーションを行うことにより、かかる被検体Ｍのモデルに対して光Ｈ１〜Ｈｎを投射したと仮定した場合の検出光Ｑ１〜Ｑｎの各光強度ｄ11〜ｄnnを計算する。
【００３３】
但し、このシミュレーションでは、被検体Ｍ中には何等の吸収体も存在しないと仮定し、更に、前記式(2) において完全拡散を想定することとし、被検体Ｍの大きさも１／μs ’より大きいとして、光ファイバＦ１〜Ｆｎの出射端から被検体Ｍの対向面までの距離をＬとする。例えば、散乱係数μs ’＝１．０ｍｍ-1であれば、Ｌ≫１ｍｍとする。また、被検体Ｍの吸収係数μa を定数とする。例えば、μa ＝０．０ｍｍ-1とする。
【００３４】
そして、かかるシミュレーションを演算部４が行い、図３中のステップ１１０に示す如く、投射光Ｈs の位置ｓ（１≦ｓ≦ｎ）を順次に切換えたときに検出される検出光Ｑt（１≦ｔ≦ｎ）の夫々の光密度ｄstを前記式(2) を用いてシミュレーションする。つまり、ｓとｔを順次に換えたときに、前記式(2) により求まる各光束φi,j,k を光密度ｄstとする。
【００３５】
尚、実際には、フィックの法則により、光流束密度ベクトルＪは、
Ｊ＝−Ｄ・ｇｒａｄ｛φ｝
となり、その光流速密度はベクトルＪの絶対値｜Ｊ｜で与えられる。しかし、光密度をｄstとする場合と対比すると、相対的に拡散係数Ｄのみの相違であり、後述する式(5) 等のように割り算により拡散係数Ｄは約分されてしまうので、前記式(2) により求まる各光束φi,j,k を光密度ｄstとしても実質的に問題を生じることがない。即ち、各体積要素ΔｖをΔｖ＝ΔｘΔｙΔｚとすれば、Ｘ座標における各体積要素の検出表面での光束φの傾きは、（φi −φi-1 ）／Δｘで表されることとなるが、Δｘ＝１、φi ＝０とすると、光流速密度ベクトルＪは、Ｊ＝Ｄφi-1 となる。したがって、後述する式(5) 等の演算を行うときに、分母と分子の約分によって、拡散係数Ｄが消えることとなるので、各光束φi,j,k を光密度ｄstとしている。
【００３６】
このシミュレーションによれば、投射光Ｈ１のみを投射したと仮定した場合の検出光Ｑ１〜Ｑｎの光密度ｄ11〜ｄ1n、投射光Ｈ２のみを投射したと仮定した場合の検出光Ｑ１〜Ｑｎの光密度ｄ21〜ｄ2n、以下同様にして、投射光Ｈｎのみを投射したと仮定した場合の検出光Ｑ１〜Ｑｎの光密度ｄn1〜ｄnnが算出される。
これらの光密度ｄst（但し、１≦ｓ≦ｎ、１≦ｔ≦ｎ）のシミュレーション結果は記憶部６に格納される。
【００３７】
（第３の処理）
次に、前記式(2)(3)を用いて第２のシミュレーションを行うことにより、各体積要素に吸収係数μa の吸収体が存在するものと仮定した場合の、検出光Ｑ１〜Ｑｎの各光密度stｄijを計算する。例えば、μa ＝０．０１ｍｍ-1とする。
【００３８】
即ち、このシミュレーションでは、被検体Ｍの一つの体積要素に吸収係数μa の吸収体が存在するものとし、その吸収体を順次に移動させた場合における、投射光Ｈ１〜Ｈｎに対する検出光Ｑ１〜Ｑｎの光密度stｄijを計算する。
【００３９】
ここで、stｄijの添字ｓは投射光Ｈs の位置、添字ｔは検出光Ｑt の位置、添字ｉ，ｊは吸収係数μa の吸収体の存在する体積要素の位置（ｉ，ｊ）を表す。
【００４０】
この第３の処理のシミュレーションを図３のステップ１２０において演算部４が行い、この第３の処理では吸収係数μa の吸収体の存在を考慮してシミュレーションし、上記第２の処理では吸収体が存在しないものとしてシミュレーションする点で相互に相違するが、前記式(2)(3)に基づいて同様の演算が行われる。
【００４１】
第２のシミュレーションによると、位置（ｉ，ｊ）＝（１，１）の体積要素に吸収体が存在するものとし且つ投射光Ｈ１のみを投射したと仮定した場合に、検出光Ｑ１〜Ｑｎの各光密度が11ｄ11〜110ｄ11として求められ、位置（ｉ，ｊ）＝（１，２）の体積要素に吸収体が存在するものとし且つ投射光Ｈ１のみを投射したと仮定した場合に、検出光Ｑ１〜Ｑｎの各光密度が11ｄ12〜110ｄ12として求められ、以下同様にして、位置（ｉ，ｊ）＝（ｎ，ｎ）に吸収体が存在するものとし且つ投射光Ｈｎのみを投射したと仮定した場合に、検出光Ｑ１〜Ｑｎの各光密度がn1ｄnn〜n10ｄnnとして求められる。
【００４２】
そして、この光密度stｄij（但し、１≦ｓ≦ｎ、１≦ｔ≦ｎ、１≦ｉ≦ｎ、１≦ｊ≦ｎ）は、記憶部６に格納される。
【００４３】
（第４の処理）
次に、光が吸収係数ａの媒質中を光路長Ｌだけ飛行した場合での光密度Ｉは、入射光密度をＩ0 とすると、ランバート・ベールの法則により、
【００４４】
【数４】
0003771339.tif 000005
となる。そこで、第１のシミュレーションと第２のシミュレーションにより求められた各検出光の光密度ｄstとstｄijを上記式(4) に適用して変形することにより次式(5) を求め、位置（ｉ，ｊ）に存在する吸収係数μa の体積要素が光に対して影響を及ぼす程度を表す新たなパラメータとしてstａijを計算する。尚、このパラメータを影響度stａijと定義することとした。
【００４５】
【数５】
0003771339.tif 000006
【００４６】
ここで、影響度stａijの符号ｓは投射光Ｈs の位置、符号ｔは検出光Ｑt の位置、符号ｉ，ｊは被検体Ｍの各体積要素の位置（ｉ，ｊ）を示し、いずれも、１≦ｓ≦ｎ、１≦ｔ≦ｎ、１≦ｉ≦ｎ、１≦ｊ≦ｎの関係にある変数である。
【００４７】
そして、演算部４が図３中のステップ１３０において、上記式(5) の演算を行うことにより、全ての体積要素についての影響度の行列Ａ＝｛stａij｝を求め、この行列Ａを記録部６に格納する。
【００４８】
図２(b) 〜(k) は、ｎ＝１０とし、１０2 個の体積要素に分割した一断層面のモデルを想定したシミュレーションの一例を示している。即ち、投射光Ｈ１を投射した場合の１０個の各検出光Ｑ１〜Ｑ１０について上記第１〜第４の処理（シミュレーション）を行うことにより得られた、影響度の計算結果stａijを表示部８に再生表示した映像を示している。
【００４９】
これらの図２(b) 〜(k) から明らかなように、位置（１，１）の体積要素に投射光Ｈ１を投射した場合に得られる検出光Ｑ１〜Ｑ１０のシミュレーション結果から影響度stａijを計算すると、投射光Ｈ１に対する各検出光Ｑ１〜Ｑ１０の各体積要素での分散分布が得られることが分かる。
【００５０】
このように、第１〜第４の処理（シミュレーション）は、予めモデル化した各体積要素についての影響度の行列Ａを算出することを目的とし、この行列Ａは、次に述べる被検体Ｍの実測データからその被検体Ｍの断層画像データを計算する際のいわゆる係数行列となる。
【００５１】
（第５の処理）
次に、図１に示す測定機構を用いて、実際に被検体Ｍに光を投射し、そのときの検出光を測定する。即ち、光ファイバＦ１〜ＦｎとＢ１〜Ｂｎによる各光プローブを被検体Ｍに接触又は対向配置し、前記シミュレーションと同等の処理により、各投射光Ｈ１〜Ｈｎを順次に切換えて投射したときの各検出光Ｑ１〜Ｑｎの光密度ｄst’（但し、添字ｓは投射光Ｈs の位置、添字ｔは検出光Ｑt の位置）を実測する。
【００５２】
更に、演算部４が図３中のステップ１４０の処理、即ち、次式(6) の演算を行うことにより、実測の光密度ｄst’と前記シミュレーションにより予め求められていた光密度ｄstとの相対比sｂtの行列Ｂ＝｛sｂt｝を算出する。つまり、この行列Ｂ＝｛sｂt｝は、投射光Ｈs に対応して検出された光の被検体Ｍ内での全影響度を表す。
【００５３】
【数６】
0003771339.tif 000007
【００５４】
（第６の処理）
被検体Ｍ内部の吸収係数の分布を行列Ｘ＝｛ｘij｝とすると、
行列Ｂ＝｛sｂt｝と影響度の行列Ａ＝｛stａij｝を用いて、次式(7) が成り立つ。
【００５５】
【数７】
0003771339.tif 000008
【００５６】
そこで、演算部４は、図３中のステップ１５０において、行列Ａの逆行列Ａ-1を算出した後、ステップ１６０において、その逆行列Ａ-1を前記式(7) に適用することによりマトリクス演算し、被検体Ｍ内部の特徴データである吸収係数の分布行列Ｘ＝｛ｘij｝を算出する。
【００５７】
この吸収係数の分布行列Ｘ＝｛ｘij｝は、各体積要素の位置（ｉ，ｊ）に対応しているので、各体積要素の特徴データを分布行列Ｘ＝｛ｘij｝で特徴抽出することができる。
【００５８】
（第７の処理）
次に、演算部４が、吸収係数の分布行列Ｘ＝｛ｘij｝の各値ｘijをＸＹ座標系の各節点の座標に位置付け、更に、各値ｘijの大きさに応じて所定の配色処理を行うことによって、断層面のカラー画像データを形成する。そして、ステップ１７０において、制御部２が、このカラー画像データを表示部８へ転送することにより、被検体Ｍの断層面をカラー表示にて再構成する。
【００５９】
図４は、図２に示したモデルに対応させて、散乱係数μs ’を１．０ｍｍ-1、所定の２領域中の体積要素に吸収係数μa が０．０１ｍｍ-1の吸収体が存在し、残余の体積要素には吸収体が存在しない（μa ＝０．０ｍｍ-1）ものとして、上記第１〜第７の処理により、吸収係数の分布行列Ｘ＝｛ｘij｝を計算してグラフィック表示した結果を示している。
【００６０】
図４から明らかなように、上記条件によるモデル（同図(a) ）に対応して、鮮明な分布行列Ｘ＝｛ｘij｝が得られ（同図(b) ）、本発明による画像再構成のアルゴリズムの妥当性が確認された。
【００６１】
このように、この実施の形態によれば、モンテカルロ法を用いず、定常光拡散方程式について差分法（前記式(2) 参照）を用いて被検体Ｍの各体積要素の特徴抽出を行うようにしたので、短時間で被検体Ｍの断層画像を再構成することが可能となる。
【００６２】
尚、説明の都合上、ＸＹ座標面上の一つの断層面の画像再構成について説明したが、実際にはＺ座標軸上の位置ｋをもパラメータにして、各体積要素の位置（ｉ，ｊ，ｋ）について上記第１〜第７の処理を行うことにより、被検体Ｍ内部の三次元断層画像を再構成している。
【００６３】
また、前記光投射プローブ及び受光プローブの構成は、この実施の形態で説明したものに限らず、いわゆる走査が可能な構成であれば、他の構成であってもよい。
【００６４】
更に又、前記（第２の処理）において、被検体Ｍの吸収係数μa を、μa ＝０．０ｍｍ-1として吸収が無いものとしたが、吸収がある場合（μa ≠０）でもよく、影響度stａijを演算する場合には被検体Ｍと異なる吸収係数の体積要素を走査することにより、前記同様にstａijを求めることができる。
【００６５】
更に又、前記（第６の処理）における図３中のステップ１５０及び１６０では、前記式(7) に示すように、逆行列Ａ-1から断層データＸ＝｛ｘij｝を直接計算する方法を採っているが、前記式(7) の代わりに、次式(8) ，(9) による逐次近似法により断層データＸの解を求めるようにしてもよい。
【００６６】
【数８】
0003771339.tif 000009
【００６７】
【数９】
0003771339.tif 000010
ここで、１≦ｉ≦ｎ、１≦ｊ≦ｎ、１≦ｓ≦ｎ、１≦ｔ≦ｎの範囲について上記式(8) ，(9) の演算を行う際に、ｎ＝１から開始するために、初期値として例えばsｂ0t＝０、ｘ0ij＝０として上記式(8) の演算を行い、次に、上記式(9) の演算を行う。そして、ｎを逐次増加していき、ｘnij−ｘn-1ij＜εとなったときに演算を終了する。尚、εは必要とする演算精度を表す定数である。また、ｘn ijは、体積要素（ｉ，ｊ）のｎ番目の吸収係数の基準値からの変化（基準値が０のときは変化の絶対値）を表し、stａijは、体積要素（ｉ，ｊ）の投射光Ｈs 対検出光Ｑt の対ｓ，ｔによる影響度（重み係数）を表す。ｄstは、投射光Ｈs 対検出光Ｑt の関係についてシミュレーションにより予め求められていた検出光Ｑt の光密度を表し、ｄst’は、投射光Ｈs 対検出光Ｑt の関係について実際に計測された検出光Ｑt の光密度を表す。sｂntは、投射光Ｈs 対検出光Ｑt の対ｓ，ｔについての相対比（減衰値）であり、第ｎ番目の値を示す。
【００６８】
（第２の実施の形態）
次に、第２の実施の形態を説明する。
【００６９】
第１の実施の形態では、被検体Ｍに対する投射光に連続光を用いた場合を示したが、第２の実施の形態では、パルス光を使用している。尚、装置の基本構成は、図１と同様であるが、光源１０から所定のタイミングでパルス光を出射させ、そのタイミングに同期して光スイッチＳＷinを順次に切換えることにより、上記パルス光を投射光Ｈs としている。
【００７０】
そして、第１の実施の形態で説明した第１〜第７の処理と同様の処理を行うことにより、最終的に吸収係数の分布行列Ｘ＝｛ｘij｝の各値ｘijをＸＹ座標系の各節点の座標に位置付け、更に、各値ｘijの大きさに応じて所定の配色処理を行うことによって、断層面のカラー画像データを形成した後、被検体Ｍの断層面をカラー表示にて再構成する。
【００７１】
この第２の実施の形態によれば、投射光Ｈs が被検体Ｍに投射された場合の検出光Ｑt をシミュレーションすると、図５に示す如く、時間分解波形の検出光Ｑt が求められる。更に、実際の被検体Ｍについて同じ条件で測定しても、同様の時間分解波形の検出光Ｑt が得られる。したがって、第１の実施の形態で説明した第１〜第７の処理をパルス光の投射光Ｈs を用いて行っても、被検体Ｍ内部を画像再構成することができる。
【００７２】
更に注目すべき点は、シミュレーション及び実測の検出光Ｑt が、図５に示した如く、時間分解波形となるので、この波形の左側部分（時間的にτより前の部分）と、右側部分（時間的にτより後の部分）等に分割して、夫々の部分の情報に基づいて画像再構成の演算処理を行うことができる。そして、上記左側部分と右側部分とでは、被検体Ｍ内部での前記影響度が異なることから、第１の実施の形態と比較して情報量を増やすことができ、高精度の画像再構成を実現することができる。
【００７３】
因みに、第１の実施の形態では、連続光を投射光Ｈs として用いるので、検出光Ｑt は、図５の上記左側部分と右側部分の合計に成るため、この第２の実施の形態と比較して情報量は少ない。
【００７４】
また、図５の検出光の上記左側部分の情報を用いると、被検体Ｍ内部を早く伝搬した光束について処理することとなるので、位置による影響度の及ぼす領域が相対的に狭まったこととなり、鮮明な断層画像を再構成することができる。
【００７５】
また、図５に示す時間分解波形の様々な時点での光密度値から画像再構成をすることにより、物理的特性の異なる種々の被検体Ｍについて、最も鮮明な断層画像を再構成するための調整処理を行うことが可能となる等の効果が得られる。
【００７６】
更に、この第２の実施の形態では、上述した差分法を適用することによって、図５に示すような雑音成分の少ない検出光が得られるが、モンテカルロ法によれば、図６の様な雑音成分を持った検出光が得られることとなるので、本実施の形態による方が、モンテカルロ法よりも、短時間で鮮明な断層画像を再構成することができる。
【００７７】
（第３の実施の形態）
次に、第３の実施の形態を説明する。
【００７８】
これは、図１に示したシステム構成と同様の装置を用いるが、演算部４において有限要素法により被検体Ｍの断層画像を再構成するものである。即ち、被検体Ｍを有限個の小部分（要素）に分解し、各要素の特性を数学的モデルで近似し、それをまとめて全体的モデルを構成して、断層画像を解くものである。
【００７９】
処理過程を示すフローチャート（図７参照）に基づいて原理を説明すると、解析対象を有限な要素に分割するものとしたときの、各要素内の光密度分布Ｐを次式(10)で表す。
【００８０】
【数１０】
0003771339.tif 000011
【００８１】
ここで、［Ｎ(x,y,z) ］は、三次元直交座標における節点光密度と各要素内部の光密度とを結び付ける内挿関数マトリックスであり、｛φ(t) ｝は時刻ｔにおける要素の節点光密度ベクトルである。
【００８２】
内挿関数マトリックス［Ｎ］を重み関数としてガレルキン（Galerkin）法を適用することにより、次式(11)を得る。
【００８３】
【数１１】
0003771339.tif 000012
【００８４】
ここで、上付き添字Ｔは転置を、ｖe は要素の領域を、μa は光の吸収係数を、Ｄは拡散係数を表す。
【００８５】
更に、上記式を更に変形して整理すると、
【００８６】
【数１２】
0003771339.tif 000013
【００８７】
ここで、
【００８８】
【数１３】
0003771339.tif 000014
【００８９】
【数１４】
0003771339.tif 000015
【００９０】
【数１５】
0003771339.tif 000016
【００９１】
尚、ｑは光流束密度、ｓは領域ｖを囲む境界を表す。
【００９２】
更に、解析対象全体の有限要素式を次式(16)で表す。
【００９３】
【数１６】
0003771339.tif 000017
【００９４】
ここで、｛Φ｝は被検体Ｍ全体の節点光密度ベクトル、［Ｋ］は光伝導マトリックス、［Ｃ］は光容量マトリックス、｛Ｆ｝は光流束マトリックスと呼ぶこととする。また、各マトリクス［Φ］、［Ｋ］、［Ｃ］、｛Ｆ｝は、次式（17-a）〜（17-d）に示すように、被検体Ｍ全体の要素ｖe についてたし合わせることを表す。
【００９５】
【数１７】
0003771339.tif 000018
【００９６】
【数１８】
0003771339.tif 000019
リッツ法に基づく有限要素法でも全く同一の有限要素式が得られる。有限要素毎のマトリックス［ｋ］とベクトル｛ｆ｝を計算し、全体のマトリックス［Ｋ］とベクトル｛Ｆ｝を合成して、連立一次方程式を解くことにより、求めたい点の光密度φを求める。
【００９７】
この計算方法を利用して被検体内の吸収係数をμa ＝μa1、影響度を計算するために移動させる吸収係数をμa ＝μa2として、第１の実施の形態において説明した差分法の場合と同様の走査方法により計算する。即ち、吸収係数をμa ＝μa1としたときの全要素の光密度（stφij）μa1と吸収係数をμa＝μa2としたときの全要素の光密度（stφij）μa2とがシミュレーションによって求まるので、前記式(5) と同様の演算を行うことによって、影響度Ａ＝｛stａij｝を求め、更に実際に計測を行って得られる検出光の光強度（光密度）とシミュレーションによって求められた光強度（光密度）との相対比の行列Ｂ＝｛sφt｝を前記式(6) と同様の演算によって求め、更に、影響度Ａと相対比の行列Ｂについて前記式(7) と同様の行列演算を行うことにより、被検体Ｍの全ての要素の吸収係数分布Ｘ＝｛ｘij｝を求める。そして、吸収係数分布Ｘに基づいて各要素を表示等することにより被検体Ｍの断層画像を再構成することができる。
【００９８】
また、第３の実施の形態において、第２の実施の形態と同様に、投射光Ｈs にパルス光を用いても断層画像を再構成することができる。
【００９９】
（第４の実施の形態）
次に、第４の実施の形態を説明する。これは、図１に示したシステム構成と同様の装置を用いるが、演算部４において有限体積（コントロール・ボリューム）法により被検体Ｍの断層画像を再構成するものである。即ち、被検体Ｍを有限個数の小部分（有限体積）に分割し、各有限体積の特性を数学的モデルで近似し、それをまとめて全体的モデルを構成して、断層画像を解くものである。
【０１００】
解析対象である被検体Ｍを有限体積に分割するものとしたときの各有限体積に入射する光拡散方程式は前記式（1)から明らかな如く、次式(19)で表される。
【０１０１】
【数１９】
0003771339.tif 000020
【０１０２】
そこで、上記式(19)を時間刻みΔｔ及び有限体積Δｖで積分すると、次式(20)が得られ、更に、式(20)を離散化して整理することで式(21)が得られる。
【０１０３】
【数２０】
0003771339.tif 000021
【０１０４】
【数２１】
0003771339.tif 000022
ここで、φは単位体積当たりの光密度、μaは各有限体積の光吸収係数、Ｄは各有限体積の拡散係数（Ｄ＝１／３μs’；但し、μs’は輸送散乱係数）、ｃは媒質中の光速である。更に、ａP ，ａE ，ａW ，ａN ，ａS ，ａT ，ａB は、光密度φの定義格子点Ｐ及びその隣接格子点Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，Ｂにおける係数を示す。即ち、図８に示す如く、格子点Ｐの位置を三次元のＸＹＺ直交座標系における位置（i ，j ，k ）とすると、格子点Ｗは（i-1 ，j ，k ）、格子点Ｅは（i+1 ，j ，k ）、格子点Ｓは（i ，j-1 ，k ）、格子点Ｎは（i ，j+1 ，k ）、格子点Ｂは（i ，j ，k-1 ）、格子点Ｔは（i ，j ，k+1 ）の各座標を表し、更に、格子点Ｐから各隣接格子点Ｅ，Ｗ，Ｎ，Ｓ，Ｔ，ＢまでのＸＹＺ座標における夫々の距離は、図９に示す如く、（δｘ）w 、（δｘ）e 、（δｙ）s 、（δｙ）n 、（δｚ）t 、（δｚ）b であり、更に、各有限体積Δｖは、Δｖ＝ΔｘΔｙΔｚの関係にあることを意味する。
【０１０５】
また、各係数ａP ，ａE ，ａW ，ａN ，ａS ，ａT ，ａB は次式（22-a）〜（22-f），（22-i）で表され、更に、時刻ｔ＝０における定義格子点Ｐでの係数ａ0P と初期値ｂは、次式（22-g）と（22-h）で表される。
【０１０６】
【数２２】
0003771339.tif 000023
【０１０７】
そこで、以上に説明した各式（21），（22-a）〜（22-i）を用いて次に述べる処理を行うことにより、断層画像を構成する。
【０１０８】
（第１の処理）
まず、被検体Ｍの各有限体積の吸収係数μa をμa＝μa1として順次に移動させ且つ、投射光Ｈs の位置ｓ（１≦ｓ≦ｎ）を順次に切換え走査したときに検出される検出光Ｑt （１≦ｔ≦ｎ）の夫々の光密度（stφij）μa1 をシミュレーションにて求める。
【０１０９】
この処理内容を更に具体的に述べると、上記式（22-a）〜（22-i）に基づいて各係数ａP ，ａE ，ａW ，ａN ，ａS ，ａT ，ａB ，ｂを求めた後、これらの係数を上記式（21）に代入することによって、各格子点における光密度φを求める。尚、この処理は図１中の演算部４が行い、演算にて求められた光密度｛stφij｝μa1 のデータが記憶部６に格納される。
【０１１０】
（第２の処理）
次に、被検体Ｍの各有限体積の吸収係数μa をμa＝μa2として順次に移動させ且つ、投射光Ｈs の位置ｓ（１≦ｓ≦ｎ）を順次に切換えたときに検出される検出光Ｑt （１≦ｔ≦ｎ）の夫々の光密度（stφij）μa2 をシミュレーションにて求める。
【０１１１】
この処理内容を更に具体的に述べると、上記式（22-a）〜（22-i）に基づいて各係数ａP ，ａE ，ａW ，ａN ，ａS ，ａT ，ａB ，ｂを求めた後、これらの係数を上記式（21）に代入することによって、各格子点における光密度φを求める。尚、この処理は図１中の演算部４が行い、演算にて求められた光密度｛stφij｝μa2 のデータが記憶部６に格納される。
【０１１２】
（第３の処理）
次に、これらの光密度（stφij）μa1 と（stφij）μa2 の比を次式（23）に基づいて演算することにより、影響度マトリクスＷＧ＝｛stｗｇij｝を求める。
【０１１３】
【数２３】
0003771339.tif 000024
【０１１４】
尚、第１，第２の処理において、光密度（stφij）μa1 と（stφij）μa2 は各格子点における光密度φに夫々求められるので、影響度マトリクスＷＧ＝｛stｗｇij｝も各格子点における光密度φに対応して求められる。
【０１１５】
そして、これらの影響度マトリクスＷＧ＝｛stｗｇij｝も記憶部６に格納される。
【０１１６】
（第４の処理）
次に、実際に被検体Ｍに対して投射光Ｈs の位置ｓをずらしながら投光し、そのときの検出光Ｑt を実測することにより、投射光に対する検出光の対ｓ，ｔの光密度（stφij）を得る。この実測された光密度（stφij）とシミュレーションによって求められている光密度（stφij）μa1 について次式（24）の演算を行うことにより、相対比の行列Ｂ＝｛sφt｝を求め、更に、前記のシミュレーションにより求められた影響マトリクスＷＧと行列Ｂについて次式（25）の演算を行うことにより、被検体Ｍ内の各有限体積の吸収係数の分布行列Ｘ＝｛ｘij｝を演算する。
【０１１７】
【数２４】
0003771339.tif 000025
【０１１８】
【数２５】
0003771339.tif 000026
このように、有限体積法によっても、断層画像を構築するための各有限体積の吸収係数を求めることができ、この分布行列Ｘに基づいて被検体Ｍの断層画像を再生することができる。
【０１１９】
尚、定常光（連続光）を用いる場合には、前記式（21），（22-a）〜（22-i）の代わりに、次式（26），（27-a）〜（27-g）を用いて上記第１の処理ないし第４の処理を行うことによって、断層画像を形成するための分布行列Ｘ＝｛ｘij｝を求めることができる。
【０１２０】
【数２６】
0003771339.tif 000027
【０１２１】
【数２７】
0003771339.tif 000028
又、以上に説明した処理方法では、全ての有限体積が所定の体積Δｖであるとして断層画像を求める場合を説明したが、被検体Ｍと外部環境とが接する境界部分等の分布行列Ｘを求めるためには、図１０（ｘ座標を代表して示す）のように半有限体積を考え、前記式（21），（22-a）〜（22-i）の代わりに、次式（28-a）〜（28-d）を用いて、前記第１，第２の処理を行う。
【０１２２】
【数２８】
0003771339.tif 000029
【０１２３】
尚、ここで、Ｉは定義格子点、ＢとＥはｘ座標における隣接格子点を意味する。また、光流束ｑB の光伝達率をｈ、周囲光密度をφf 、境界面Ｂでの光密度をφB とすれば、ｑB＝ｈ（φf−φB）となり、ｂ＝ｈφf 、ａB ＝ａI ＋μa Δｘ＋ｈとなる。
【０１２４】
また、第４の実施の形態において、第２の実施の形態と同様に、投射光Ｈs にパルス光を用いても断層画像を再構成することができる。
【０１２５】
【発明の効果】
以上に説明したように本発明によれば、従来のモンテカルロ法によらず、差分法等を適用した所定アルゴリズムに基づいて被検体の各体積要素の特徴抽出を行い、その特徴抽出データに基づいて断層画像を再構成するので、演算速度が速く高速の画像再構成が可能である。
【０１２６】
この結果、例えば生体の断層画像を再構成するための医療用光ＣＴ装置等の実用化等に貢献することができる優れた効果を発揮する。
【図面の簡単な説明】
【図１】実施の形態に係る光ＣＴ装置の構成を示す構成説明図である。
【図２】画像再構成の原理を説明するため、ディスプレイ上に表示した中間階調画像を示した図である。
【図３】光ＣＴ装置の動作を処理過程に沿って示すフローチャートである。
【図４】画像再構成の一例を説明するため、ディスプレイ上に表示した中間階調画像を示した図である。
【図５】パルス光を投射光に適用した場合に、差分法でシミュレーションして得られる検出光の波形を示すグラフである。
【図６】パルス光を投射光に適用した場合に、モンテカルロ法でシミュレーションして得られる検出光の波形を示すグラフである。
【図７】第３の実施の形態の処理行程を説明するためのフローチャートである。
【図８】第４の実施の形態の原理を説明するための説明図である。
【図９】第４の実施の形態の原理を更に説明するための説明図である。
【図１０】第４の実施の形態の原理を更に説明するための説明図である。
【符号の説明】
２…制御部、４…演算部、６…記憶部、８…表示部、１０…光源、１２…光検出部、ＳＷin，ＳＷout …光スイッチ、Ｆ１〜Ｆｎ、Ｂ１〜Ｂｎ…光ファイバ、Ｍ…被検体、Ｈ１〜Ｈｎ…投射光、Ｑ１〜Ｑｎ…検出光。
【図１】
0003771339.tif 000030
【図２】
0003771339.tif 000031
【図３】
0003771339.tif 000032
【図４】
0003771339.tif 000033
【図５】
0003771339.tif 000034
【図６】
0003771339.tif 000035
【図７】
0003771339.tif 000036
【図８】
0003771339.tif 000037
【図９】
0003771339.tif 000038
【図１０】
0003771339.tif 000039