JP3771339B2 - 光ct装置及び光ctによる画像再構成方法 - Google Patents

光ct装置及び光ctによる画像再構成方法 Download PDF

Info

Publication number
JP3771339B2
JP3771339B2 JP00781097A JP781097A JP3771339B2 JP 3771339 B2 JP3771339 B2 JP 3771339B2 JP 00781097 A JP00781097 A JP 00781097A JP 781097 A JP781097 A JP 781097A JP 3771339 B2 JP3771339 B2 JP 3771339B2
Authority
JP
Japan
Prior art keywords
light
subject
matrix
light source
computing
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
Application number
JP00781097A
Other languages
English (en)
Other versions
JPH09257694A (ja
Inventor
和義 太田
之雄 上田
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hamamatsu Photonics KK
Original Assignee
Hamamatsu Photonics KK
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hamamatsu Photonics KK filed Critical Hamamatsu Photonics KK
Priority to JP00781097A priority Critical patent/JP3771339B2/ja
Publication of JPH09257694A publication Critical patent/JPH09257694A/ja
Application granted granted Critical
Publication of JP3771339B2 publication Critical patent/JP3771339B2/ja
Anticipated expiration legal-status Critical
Expired - Lifetime legal-status Critical Current

Links

Images

Landscapes

  • Investigating Or Analysing Materials By Optical Means (AREA)

Description

【0001】
【発明の属する技術分野】
本発明は、生体等の被検体に光を投射し、その透過光の測定値から被検体の断層画像を再構成する光CT装置及び光CTによる画像再構成方法に関する。
【0002】
【従来の技術】
従来の光CT装置は、生体等の被検体に光を投射すると共に被検体内部を透過して来た光を検出器で測定し、その測定結果に基づいて、その光(光束)が生体内の特定領域を通過する際の通過頻度(平均光路長の空間分布)をモンテカルロシミュレーションにより求めて、断層画像を再構成するアルゴリズムを採用していた。つまり、被検体を複数の体積要素(VOXEL)に分割したモデルを想定し、各体積要素ごとに光束の衝突回数又は平均光路長(各体積要素を横切る飛行距離の総和または平均値)を計算することにより、画像再構成に必要な各体積要素の特徴抽出を行うこととしていた。
【0003】
【発明が解決しようとする課題】
しかし、従来の光CT装置にあっては、モンテカルロシミュレーションは、光束の散乱粒子での散乱事象を、乱数を用いて計算する光線追跡法であるため、実際の測定では検出器に到達しない光束があったとしても、想定され得る全ての光束について計算することとなる結果、無駄な情報(事象)についてまでも計算対象とするため、演算に要する時間が極めて長くなるという問題があった。
【0004】
特に、高品質の断層画像を再構成しようとすれば、当然に計算対象が増加し、無駄な情報に関する演算所要時間の占める割合が激増することになるので、実用に供し得る光CT装置を実現することが困難となっていた。
【0005】
本発明は、このような従来の技術に鑑みてなされたものであり、被検体の断層画像を高速に再構成することができる光CT装置を提供することを目的とする。
【0006】
【課題を解決するための手段】
このような目的を達成するために本発明は、被検体の複数の部位に個別排他的に投射光を投射する光源手段と、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、前記被検体を微少な複数個の体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第1の演算手段と、前記第1の演算手段で想定した前記集合体モデルの特定の体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第2の演算手段と、前記第1の演算手段で求められた前記光密度と前記第2の演算手段で求められた光密度との比を前記体積要素の影響度として演算して、前記複数の体積要素の配列に対応する影響度の行列を求める第3の演算手段と、前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第1の演算手段で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める第4の演算手段と、前記第3の演算手段により求められた影響度の行列と前記第4の演算手段により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める第5の演算手段と、を具備することとした。
【0007】
また、被検体の複数の部位に個別排他的に投射光を投射する光源手段と、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、前記被検体を複数個の要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第1の演算手段と、前記第1の演算手段で想定した前記集合体モデルの特定の要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第2の演算手段と、前記第1の演算手段で求められた光密度と前記第2の演算手段で求められた光密度との比を前記要素の影響度として演算して、前記複数の要素の配列に対応する影響度の行列を求める第3の演算手段と、前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第1の演算手段で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める第4の演算手段と、前記第3の演算手段により求められた影響度の行列と前記第4の演算手段により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める第5の演算手段と、を具備することとした。
【0008】
また、被検体の複数の部位に個別排他的に投射光を投射する光源手段と、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、前記被検体を複数の有限体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第1の演算手段と、前記第1の演算手段で想定した前記集合体モデルの特定の有限体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第2の演算手段と、前記第1の演算手段で求められた光密度と前記第2の演算手段で求められた光密度との比を前記有限体積要素の影響度として演算して、前記複数の有限体積要素の配列に対応する影響度の行列を求める第3の演算手段と、前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と、前記第1の演算手段で求められた光強度との相対比を演算し、前記複数の有限体積要素の配列に対応する相対比の行列を求める第4の演算手段と、前記第3の演算手段により求められた影響度の行列と前記第4の演算手段により求められた相対比の行列とのマトリクス演算により、前記各有限体積要素の特徴データを求める第5の演算手段と、を具備することとした。
【0009】
また、前記光源手段は、連続光を前記投射光として出射する構成とした。
【0010】
また、前記光源手段は、パルス光を前記投射光として出射し、前記光検出手段は、前記透過光を時間分解計測する構成とした。
【0011】
【発明の実施の形態】
(第1の実施の形態)
本発明の第1の実施の形態を図面と共に説明する。
【0012】
図1は、本実施の形態に係る光CT装置の構成を示す。同図において、この光CT装置は、被検体Mの光学特性を実測するための測定機構と、測定機構により得られた実測データを演算処理することにより被検体Mの断層画像を再構成する演算機構を備えている。
【0013】
演算機構は、本装置全体の動作を制御するマイクロプロセッサユニット(MPU)を有する制御部2と、予め決められたアルゴリズムに基づいて作成されたプログラムを実行することにより被検体Mの断層画像を再構成するための演算処理を行う演算部4と、演算部4が演算処理の際に種々のデータを格納等するための記憶部6と、上記演算処理により得られた被検体Mの断層画像等を表示するための表示部8とを具備している。
【0014】
測定機構は、被検体Mに所定波長の光を投射するための光投射部と、被検体Mの透過光を実測するための受光部とを備えている。
【0015】
即ち、上記光投射部にあっては、所定波長の光を出射する光源10の光出射端に、1入力n出力の光分岐路を有する光スイッチSWinの入力端が接続されると共に、光スイッチSWinの各出力端にn本の光ファイバF1〜Fnが接続されている。更に、n本の光ファイバF1〜Fnは、これらの光出射端が面一に揃えられて同一方向に沿って一列に配列され、機械的に安定ないわゆる光投射プローブとなっている。そして、光ファイバF1〜Fnの光出射端を被検体Mの一側端に対向又は接触させて使用される。
【0016】
上記光スイッチSWinは、制御部2からの切換制御に従って、n個中の1つの分岐路のみを選択的にオン、残余の分岐路をオフにする機能を有することにより、光源10からの光を個別排他的に各光ファイバF1〜Fnに導入し、その光の導入された光ファイバの光出射端のみから被検体Mに光を投射する。従って、被検体Mの位置の異なる部位へ順次に光を投射する等の走査が可能となっている。
【0017】
一方、上記受光部にあっては、n本の光ファイバB1〜Bnと、これらの光ファイバB1〜Bnの光出射端に接続されたn入力1出力の光分岐路を有する光スイッチSWout と、光スイッチSWout の光出力端に接続された光電変換デバイスなどを有する光検出部12を具備している。そして、n本の光ファイバB1〜Bnは、夫々の光入射端が面一に揃えられて同一方向に沿って一列に配列され且つ機械的に安定ないわゆる受光プローブとなっており、これらの光ファイバB1〜Bnの光入射端を被検体Mの一側端に対向又は接触させて使用される。
【0018】
上記光スイッチSWout は、制御部2からの切換制御に従って、n個中の1つの分岐路のみを選択的にオン、残余の分岐路をオフにする機能を有することにより、各光ファイバB1〜Bnからの光を個別排他的に光検出部12中の上記光電変換デバイスへ導出する。したがって、被検体Mの位置の異なる部位からの光を順次に検出する等の検出走査が可能となっている。
【0019】
そして、被検体Mを実測する際には、光ファイバF1〜Fnによる光投射プローブと光ファイバB1〜Bnによる受光プローブとを被検体Mを介して対向配置し、光源10からの光を各光ファイバF1〜Fnにより順次に被検体Mに投射し、被検体Mの内部を透過してきた光(光束)を各光ファイバB1〜Bnで順次に導入して光検出部12で光電変換及びA/D変換することにより、被検体Mの各部位の光学特性のデータを得るようになっている。
【0020】
また、図1に示す如く、光ファイバF1〜Fn(光投射プローブ)のn個の光出射端と光ファイバB1〜Bn(受光プローブ)のn個の光入射端によって画成される領域が、画像再構成されるべき被検体Mの一つの断層面に対応する。よって、光ファイバF1〜Fn(光投射プローブ)と光ファイバB1〜Bn(受光プローブ)を、被検体Mに対して、図1中の縦方向(Z座標軸方向)に相対移動させることにより、複数の断層面を画成することができ、これらの複数の断層面についての実測データに基づいて後述の画像処理を行うことにより、被検体M内部の三次元断層画像を再構成する。
【0021】
尚、説明の都合上、被検体MがXYZ直交座標系中に配置され、光ファイバF1〜Fn及び光ファイバB1〜BnはY座標軸方向に配列され、光ファイバF1〜Fnの光出射端と光ファイバB1〜Bnが光入射端がX座標軸方向に沿って対向配置されるものとする。更に、装置構成として、光スイッチの代わりに多数個の光源をオン・オフしたり、多数個の検出器を同時測定するようにしてもよい。
【0022】
次に、この光CT装置の動作及び被検体Mの断層画像の再構成原理を、実際の処理過程に対応させて説明する。尚、図2は、再構成原理を説明するための説明図、図3は、動作を処理過程に沿って示すフローチャート、図4は、本実施の形態の再構成原理の妥当性を示す実験結果を示す。
【0023】
(第1の処理)
図3のステップ100に示す如く、光ファイバF1〜Fn又は光ファイバB1〜Bnの各ピッチ間隔で、被検体MをXYZ座標上でn×n×n個の体積要素(VOXEL)に分割するものとし、かかる設定条件は制御部2及び演算部4に予め格納される。例えば、n=100として、被検体Mを1003 個の体積要素の集合体として扱う。
【0024】
更に、演算部4には、次に述べるアルゴリズムに基づく演算プログラムが予め格納されている。
【0025】
各体積要素に入射する連続光(光束)の定常光拡散方程式は、
【0026】
【数1】
Figure 0003771339
で表されるので、この定常拡散方程式(1) を差分形式で表した次式(2) が演算部4に予め格納されている。
【0027】
【数2】
Figure 0003771339
【0028】
ここで、添字i,j,kは、XYZ直交座標系の節点の座標を表し、φi,j,k は、XYZ座標系の位置(i ,j ,k )に存在する体積要素を通る光束を表し、1≦i≦n、1≦j≦n、1≦k≦nである。
【0029】
また、被検体Mの内部と外部との境界条件は、
【0030】
【数3】
Figure 0003771339
としている。尚、添字BLは、被検体M周辺の境界を表すものとする。また、式(3) は、被検体Mの表面で光が完全に吸収されてしまう様な条件、例えば、被検体Mの表面が真黒に塗られた状態と等価である。
【0031】
これらの式(2)(3)を用いて、被検体Mの断層画像を再構成するのであるが、説明を明確にするために、図2(a) に示す如く、一つの断層面に注目して述べることとする。即ち、XY座標面上の一つの断層面を、n本ずつの光ファイバF1〜FnとB1〜Bnに対応させて、n×n個の体積要素(VOXEL)に分割するものとする。そして、光ファイバF1〜Fnの各光出射端より投射される各投射光をH1〜Hnとすると共に、一般表現のHs (1≦s≦n)で総称するものとし、光ファイバB1〜Bnの各光入射端で検出される検出光をQ1〜Qnとすると共に、一般表現のQt(1≦t≦n)で総称するものとし、更に、各体積要素のXY座標上の位置を(i ,j )で表すものとする。
【0032】
(第2の処理)
次に、前記式(2)(3)を用いて第1のシミュレーションを行うことにより、かかる被検体Mのモデルに対して光H1〜Hnを投射したと仮定した場合の検出光Q1〜Qnの各光強度d11〜dnnを計算する。
【0033】
但し、このシミュレーションでは、被検体M中には何等の吸収体も存在しないと仮定し、更に、前記式(2) において完全拡散を想定することとし、被検体Mの大きさも1/μs ’より大きいとして、光ファイバF1〜Fnの出射端から被検体Mの対向面までの距離をLとする。例えば、散乱係数μs ’=1.0mm-1であれば、L≫1mmとする。また、被検体Mの吸収係数μa を定数とする。例えば、μa =0.0mm-1とする。
【0034】
そして、かかるシミュレーションを演算部4が行い、図3中のステップ110に示す如く、投射光Hs の位置s(1≦s≦n)を順次に切換えたときに検出される検出光Qt(1≦t≦n)の夫々の光密度dstを前記式(2) を用いてシミュレーションする。つまり、sとtを順次に換えたときに、前記式(2) により求まる各光束φi,j,k を光密度dstとする。
【0035】
尚、実際には、フィックの法則により、光流束密度ベクトルJは、
J=−D・grad{φ}
となり、その光流速密度はベクトルJの絶対値|J|で与えられる。しかし、光密度をdstとする場合と対比すると、相対的に拡散係数Dのみの相違であり、後述する式(5) 等のように割り算により拡散係数Dは約分されてしまうので、前記式(2) により求まる各光束φi,j,k を光密度dstとしても実質的に問題を生じることがない。即ち、各体積要素ΔvをΔv=ΔxΔyΔzとすれば、X座標における各体積要素の検出表面での光束φの傾きは、(φi −φi-1 )/Δxで表されることとなるが、Δx=1、φi =0とすると、光流速密度ベクトルJは、J=Dφi-1 となる。したがって、後述する式(5) 等の演算を行うときに、分母と分子の約分によって、拡散係数Dが消えることとなるので、各光束φi,j,k を光密度dstとしている。
【0036】
このシミュレーションによれば、投射光H1のみを投射したと仮定した場合の検出光Q1〜Qnの光密度d11〜d1n、投射光H2のみを投射したと仮定した場合の検出光Q1〜Qnの光密度d21〜d2n、以下同様にして、投射光Hnのみを投射したと仮定した場合の検出光Q1〜Qnの光密度dn1〜dnnが算出される。
これらの光密度dst(但し、1≦s≦n、1≦t≦n)のシミュレーション結果は記憶部6に格納される。
【0037】
(第3の処理)
次に、前記式(2)(3)を用いて第2のシミュレーションを行うことにより、各体積要素に吸収係数μa の吸収体が存在するものと仮定した場合の、検出光Q1〜Qnの各光密度stijを計算する。例えば、μa =0.01mm-1とする。
【0038】
即ち、このシミュレーションでは、被検体Mの一つの体積要素に吸収係数μa の吸収体が存在するものとし、その吸収体を順次に移動させた場合における、投射光H1〜Hnに対する検出光Q1〜Qnの光密度stijを計算する。
【0039】
ここで、stijの添字sは投射光Hs の位置、添字tは検出光Qt の位置、添字i,jは吸収係数μa の吸収体の存在する体積要素の位置(i,j)を表す。
【0040】
この第3の処理のシミュレーションを図3のステップ120において演算部4が行い、この第3の処理では吸収係数μa の吸収体の存在を考慮してシミュレーションし、上記第2の処理では吸収体が存在しないものとしてシミュレーションする点で相互に相違するが、前記式(2)(3)に基づいて同様の演算が行われる。
【0041】
第2のシミュレーションによると、位置(i,j)=(1,1)の体積要素に吸収体が存在するものとし且つ投射光H1のみを投射したと仮定した場合に、検出光Q1〜Qnの各光密度が111111011として求められ、位置(i,j)=(1,2)の体積要素に吸収体が存在するものとし且つ投射光H1のみを投射したと仮定した場合に、検出光Q1〜Qnの各光密度が111211012として求められ、以下同様にして、位置(i,j)=(n,n)に吸収体が存在するものとし且つ投射光Hnのみを投射したと仮定した場合に、検出光Q1〜Qnの各光密度がn1nnn10nnとして求められる。
【0042】
そして、この光密度stij(但し、1≦s≦n、1≦t≦n、1≦i≦n、1≦j≦n)は、記憶部6に格納される。
【0043】
(第4の処理)
次に、光が吸収係数aの媒質中を光路長Lだけ飛行した場合での光密度Iは、入射光密度をI0 とすると、ランバート・ベールの法則により、
【0044】
【数4】
Figure 0003771339
となる。そこで、第1のシミュレーションと第2のシミュレーションにより求められた各検出光の光密度dststijを上記式(4) に適用して変形することにより次式(5) を求め、位置(i,j)に存在する吸収係数μa の体積要素が光に対して影響を及ぼす程度を表す新たなパラメータとしてstijを計算する。尚、このパラメータを影響度stijと定義することとした。
【0045】
【数5】
Figure 0003771339
【0046】
ここで、影響度stijの符号sは投射光Hs の位置、符号tは検出光Qt の位置、符号i,jは被検体Mの各体積要素の位置(i,j)を示し、いずれも、1≦s≦n、1≦t≦n、1≦i≦n、1≦j≦nの関係にある変数である。
【0047】
そして、演算部4が図3中のステップ130において、上記式(5) の演算を行うことにより、全ての体積要素についての影響度の行列A={stij}を求め、この行列Aを記録部6に格納する。
【0048】
図2(b) 〜(k) は、n=10とし、102 個の体積要素に分割した一断層面のモデルを想定したシミュレーションの一例を示している。即ち、投射光H1を投射した場合の10個の各検出光Q1〜Q10について上記第1〜第4の処理(シミュレーション)を行うことにより得られた、影響度の計算結果stijを表示部8に再生表示した映像を示している。
【0049】
これらの図2(b) 〜(k) から明らかなように、位置(1,1)の体積要素に投射光H1を投射した場合に得られる検出光Q1〜Q10のシミュレーション結果から影響度stijを計算すると、投射光H1に対する各検出光Q1〜Q10の各体積要素での分散分布が得られることが分かる。
【0050】
このように、第1〜第4の処理(シミュレーション)は、予めモデル化した各体積要素についての影響度の行列Aを算出することを目的とし、この行列Aは、次に述べる被検体Mの実測データからその被検体Mの断層画像データを計算する際のいわゆる係数行列となる。
【0051】
(第5の処理)
次に、図1に示す測定機構を用いて、実際に被検体Mに光を投射し、そのときの検出光を測定する。即ち、光ファイバF1〜FnとB1〜Bnによる各光プローブを被検体Mに接触又は対向配置し、前記シミュレーションと同等の処理により、各投射光H1〜Hnを順次に切換えて投射したときの各検出光Q1〜Qnの光密度dst’(但し、添字sは投射光Hs の位置、添字tは検出光Qt の位置)を実測する。
【0052】
更に、演算部4が図3中のステップ140の処理、即ち、次式(6) の演算を行うことにより、実測の光密度dst’と前記シミュレーションにより予め求められていた光密度dstとの相対比stの行列B={st}を算出する。つまり、この行列B={st}は、投射光Hs に対応して検出された光の被検体M内での全影響度を表す。
【0053】
【数6】
Figure 0003771339
【0054】
(第6の処理)
被検体M内部の吸収係数の分布を行列X={xij}とすると、
行列B={st}と影響度の行列A={stij}を用いて、次式(7) が成り立つ。
【0055】
【数7】
Figure 0003771339
【0056】
そこで、演算部4は、図3中のステップ150において、行列Aの逆行列A-1を算出した後、ステップ160において、その逆行列A-1を前記式(7) に適用することによりマトリクス演算し、被検体M内部の特徴データである吸収係数の分布行列X={xij}を算出する。
【0057】
この吸収係数の分布行列X={xij}は、各体積要素の位置(i,j)に対応しているので、各体積要素の特徴データを分布行列X={xij}で特徴抽出することができる。
【0058】
(第7の処理)
次に、演算部4が、吸収係数の分布行列X={xij}の各値xijをXY座標系の各節点の座標に位置付け、更に、各値xijの大きさに応じて所定の配色処理を行うことによって、断層面のカラー画像データを形成する。そして、ステップ170において、制御部2が、このカラー画像データを表示部8へ転送することにより、被検体Mの断層面をカラー表示にて再構成する。
【0059】
図4は、図2に示したモデルに対応させて、散乱係数μs ’を1.0mm-1、所定の2領域中の体積要素に吸収係数μa が0.01mm-1の吸収体が存在し、残余の体積要素には吸収体が存在しない(μa =0.0mm-1)ものとして、上記第1〜第7の処理により、吸収係数の分布行列X={xij}を計算してグラフィック表示した結果を示している。
【0060】
図4から明らかなように、上記条件によるモデル(同図(a) )に対応して、鮮明な分布行列X={xij}が得られ(同図(b) )、本発明による画像再構成のアルゴリズムの妥当性が確認された。
【0061】
このように、この実施の形態によれば、モンテカルロ法を用いず、定常光拡散方程式について差分法(前記式(2) 参照)を用いて被検体Mの各体積要素の特徴抽出を行うようにしたので、短時間で被検体Mの断層画像を再構成することが可能となる。
【0062】
尚、説明の都合上、XY座標面上の一つの断層面の画像再構成について説明したが、実際にはZ座標軸上の位置kをもパラメータにして、各体積要素の位置(i,j,k)について上記第1〜第7の処理を行うことにより、被検体M内部の三次元断層画像を再構成している。
【0063】
また、前記光投射プローブ及び受光プローブの構成は、この実施の形態で説明したものに限らず、いわゆる走査が可能な構成であれば、他の構成であってもよい。
【0064】
更に又、前記(第2の処理)において、被検体Mの吸収係数μa を、μa =0.0mm-1として吸収が無いものとしたが、吸収がある場合(μa ≠0)でもよく、影響度stijを演算する場合には被検体Mと異なる吸収係数の体積要素を走査することにより、前記同様にstijを求めることができる。
【0065】
更に又、前記(第6の処理)における図3中のステップ150及び160では、前記式(7) に示すように、逆行列A-1から断層データX={xij}を直接計算する方法を採っているが、前記式(7) の代わりに、次式(8) ,(9) による逐次近似法により断層データXの解を求めるようにしてもよい。
【0066】
【数8】
Figure 0003771339
【0067】
【数9】
Figure 0003771339
ここで、1≦i≦n、1≦j≦n、1≦s≦n、1≦t≦nの範囲について上記式(8) ,(9) の演算を行う際に、n=1から開始するために、初期値として例えばs0 t=0、x0 ij=0として上記式(8) の演算を行い、次に、上記式(9) の演算を行う。そして、nを逐次増加していき、xn ij−xn-1 ij<εとなったときに演算を終了する。尚、εは必要とする演算精度を表す定数である。また、xn ijは、体積要素(i,j)のn番目の吸収係数の基準値からの変化(基準値が0のときは変化の絶対値)を表し、stijは、体積要素(i,j)の投射光Hs 対検出光Qt の対s,tによる影響度(重み係数)を表す。dstは、投射光Hs 対検出光Qt の関係についてシミュレーションにより予め求められていた検出光Qt の光密度を表し、dst’は、投射光Hs 対検出光Qt の関係について実際に計測された検出光Qt の光密度を表す。sn tは、投射光Hs 対検出光Qt の対s,tについての相対比(減衰値)であり、第n番目の値を示す。
【0068】
(第2の実施の形態)
次に、第2の実施の形態を説明する。
【0069】
第1の実施の形態では、被検体Mに対する投射光に連続光を用いた場合を示したが、第2の実施の形態では、パルス光を使用している。尚、装置の基本構成は、図1と同様であるが、光源10から所定のタイミングでパルス光を出射させ、そのタイミングに同期して光スイッチSWinを順次に切換えることにより、上記パルス光を投射光Hs としている。
【0070】
そして、第1の実施の形態で説明した第1〜第7の処理と同様の処理を行うことにより、最終的に吸収係数の分布行列X={xij}の各値xijをXY座標系の各節点の座標に位置付け、更に、各値xijの大きさに応じて所定の配色処理を行うことによって、断層面のカラー画像データを形成した後、被検体Mの断層面をカラー表示にて再構成する。
【0071】
この第2の実施の形態によれば、投射光Hs が被検体Mに投射された場合の検出光Qt をシミュレーションすると、図5に示す如く、時間分解波形の検出光Qt が求められる。更に、実際の被検体Mについて同じ条件で測定しても、同様の時間分解波形の検出光Qt が得られる。したがって、第1の実施の形態で説明した第1〜第7の処理をパルス光の投射光Hs を用いて行っても、被検体M内部を画像再構成することができる。
【0072】
更に注目すべき点は、シミュレーション及び実測の検出光Qt が、図5に示した如く、時間分解波形となるので、この波形の左側部分(時間的にτより前の部分)と、右側部分(時間的にτより後の部分)等に分割して、夫々の部分の情報に基づいて画像再構成の演算処理を行うことができる。そして、上記左側部分と右側部分とでは、被検体M内部での前記影響度が異なることから、第1の実施の形態と比較して情報量を増やすことができ、高精度の画像再構成を実現することができる。
【0073】
因みに、第1の実施の形態では、連続光を投射光Hs として用いるので、検出光Qt は、図5の上記左側部分と右側部分の合計に成るため、この第2の実施の形態と比較して情報量は少ない。
【0074】
また、図5の検出光の上記左側部分の情報を用いると、被検体M内部を早く伝搬した光束について処理することとなるので、位置による影響度の及ぼす領域が相対的に狭まったこととなり、鮮明な断層画像を再構成することができる。
【0075】
また、図5に示す時間分解波形の様々な時点での光密度値から画像再構成をすることにより、物理的特性の異なる種々の被検体Mについて、最も鮮明な断層画像を再構成するための調整処理を行うことが可能となる等の効果が得られる。
【0076】
更に、この第2の実施の形態では、上述した差分法を適用することによって、図5に示すような雑音成分の少ない検出光が得られるが、モンテカルロ法によれば、図6の様な雑音成分を持った検出光が得られることとなるので、本実施の形態による方が、モンテカルロ法よりも、短時間で鮮明な断層画像を再構成することができる。
【0077】
(第3の実施の形態)
次に、第3の実施の形態を説明する。
【0078】
これは、図1に示したシステム構成と同様の装置を用いるが、演算部4において有限要素法により被検体Mの断層画像を再構成するものである。即ち、被検体Mを有限個の小部分(要素)に分解し、各要素の特性を数学的モデルで近似し、それをまとめて全体的モデルを構成して、断層画像を解くものである。
【0079】
処理過程を示すフローチャート(図7参照)に基づいて原理を説明すると、解析対象を有限な要素に分割するものとしたときの、各要素内の光密度分布Pを次式(10)で表す。
【0080】
【数10】
Figure 0003771339
【0081】
ここで、[N(x,y,z) ]は、三次元直交座標における節点光密度と各要素内部の光密度とを結び付ける内挿関数マトリックスであり、{φ(t) }は時刻tにおける要素の節点光密度ベクトルである。
【0082】
内挿関数マトリックス[N]を重み関数としてガレルキン(Galerkin)法を適用することにより、次式(11)を得る。
【0083】
【数11】
Figure 0003771339
【0084】
ここで、上付き添字Tは転置を、ve は要素の領域を、μa は光の吸収係数を、Dは拡散係数を表す。
【0085】
更に、上記式を更に変形して整理すると、
【0086】
【数12】
Figure 0003771339
【0087】
ここで、
【0088】
【数13】
Figure 0003771339
【0089】
【数14】
Figure 0003771339
【0090】
【数15】
Figure 0003771339
【0091】
尚、qは光流束密度、sは領域vを囲む境界を表す。
【0092】
更に、解析対象全体の有限要素式を次式(16)で表す。
【0093】
【数16】
Figure 0003771339
【0094】
ここで、{Φ}は被検体M全体の節点光密度ベクトル、[K]は光伝導マトリックス、[C]は光容量マトリックス、{F}は光流束マトリックスと呼ぶこととする。また、各マトリクス[Φ]、[K]、[C]、{F}は、次式(17-a)〜(17-d)に示すように、被検体M全体の要素ve についてたし合わせることを表す。
【0095】
【数17】
Figure 0003771339
【0096】
【数18】
Figure 0003771339
リッツ法に基づく有限要素法でも全く同一の有限要素式が得られる。有限要素毎のマトリックス[k]とベクトル{f}を計算し、全体のマトリックス[K]とベクトル{F}を合成して、連立一次方程式を解くことにより、求めたい点の光密度φを求める。
【0097】
この計算方法を利用して被検体内の吸収係数をμa =μa1、影響度を計算するために移動させる吸収係数をμa =μa2として、第1の実施の形態において説明した差分法の場合と同様の走査方法により計算する。即ち、吸収係数をμa =μa1としたときの全要素の光密度(stφij)μa1と吸収係数をμa=μa2としたときの全要素の光密度(stφij)μa2とがシミュレーションによって求まるので、前記式(5) と同様の演算を行うことによって、影響度A={stij}を求め、更に実際に計測を行って得られる検出光の光強度(光密度)とシミュレーションによって求められた光強度(光密度)との相対比の行列B={sφt}を前記式(6) と同様の演算によって求め、更に、影響度Aと相対比の行列Bについて前記式(7) と同様の行列演算を行うことにより、被検体Mの全ての要素の吸収係数分布X={xij}を求める。そして、吸収係数分布Xに基づいて各要素を表示等することにより被検体Mの断層画像を再構成することができる。
【0098】
また、第3の実施の形態において、第2の実施の形態と同様に、投射光Hs にパルス光を用いても断層画像を再構成することができる。
【0099】
(第4の実施の形態)
次に、第4の実施の形態を説明する。これは、図1に示したシステム構成と同様の装置を用いるが、演算部4において有限体積(コントロール・ボリューム)法により被検体Mの断層画像を再構成するものである。即ち、被検体Mを有限個数の小部分(有限体積)に分割し、各有限体積の特性を数学的モデルで近似し、それをまとめて全体的モデルを構成して、断層画像を解くものである。
【0100】
解析対象である被検体Mを有限体積に分割するものとしたときの各有限体積に入射する光拡散方程式は前記式(1)から明らかな如く、次式(19)で表される。
【0101】
【数19】
Figure 0003771339
【0102】
そこで、上記式(19)を時間刻みΔt及び有限体積Δvで積分すると、次式(20)が得られ、更に、式(20)を離散化して整理することで式(21)が得られる。
【0103】
【数20】
Figure 0003771339
【0104】
【数21】
Figure 0003771339
ここで、φは単位体積当たりの光密度、μaは各有限体積の光吸収係数、Dは各有限体積の拡散係数(D=1/3μs’;但し、μs’は輸送散乱係数)、cは媒質中の光速である。更に、aP ,aE ,aW ,aN ,aS ,aT ,aB は、光密度φの定義格子点P及びその隣接格子点E,W,N,S,T,Bにおける係数を示す。即ち、図8に示す如く、格子点Pの位置を三次元のXYZ直交座標系における位置(i ,j ,k )とすると、格子点Wは(i-1 ,j ,k )、格子点Eは(i+1 ,j ,k )、格子点Sは(i ,j-1 ,k )、格子点Nは(i ,j+1 ,k )、格子点Bは(i ,j ,k-1 )、格子点Tは(i ,j ,k+1 )の各座標を表し、更に、格子点Pから各隣接格子点E,W,N,S,T,BまでのXYZ座標における夫々の距離は、図9に示す如く、(δx)w 、(δx)e 、(δy)s 、(δy)n 、(δz)t 、(δz)b であり、更に、各有限体積Δvは、Δv=ΔxΔyΔzの関係にあることを意味する。
【0105】
また、各係数aP ,aE ,aW ,aN ,aS ,aT ,aB は次式(22-a)〜(22-f),(22-i)で表され、更に、時刻t=0における定義格子点Pでの係数a0 P と初期値bは、次式(22-g)と(22-h)で表される。
【0106】
【数22】
Figure 0003771339
【0107】
そこで、以上に説明した各式(21),(22-a)〜(22-i)を用いて次に述べる処理を行うことにより、断層画像を構成する。
【0108】
(第1の処理)
まず、被検体Mの各有限体積の吸収係数μa をμa=μa1として順次に移動させ且つ、投射光Hs の位置s(1≦s≦n)を順次に切換え走査したときに検出される検出光Qt (1≦t≦n)の夫々の光密度(stφij)μa1 をシミュレーションにて求める。
【0109】
この処理内容を更に具体的に述べると、上記式(22-a)〜(22-i)に基づいて各係数aP ,aE ,aW ,aN ,aS ,aT ,aB ,bを求めた後、これらの係数を上記式(21)に代入することによって、各格子点における光密度φを求める。尚、この処理は図1中の演算部4が行い、演算にて求められた光密度{stφij}μa1 のデータが記憶部6に格納される。
【0110】
(第2の処理)
次に、被検体Mの各有限体積の吸収係数μa をμa=μa2として順次に移動させ且つ、投射光Hs の位置s(1≦s≦n)を順次に切換えたときに検出される検出光Qt (1≦t≦n)の夫々の光密度(stφij)μa2 をシミュレーションにて求める。
【0111】
この処理内容を更に具体的に述べると、上記式(22-a)〜(22-i)に基づいて各係数aP ,aE ,aW ,aN ,aS ,aT ,aB ,bを求めた後、これらの係数を上記式(21)に代入することによって、各格子点における光密度φを求める。尚、この処理は図1中の演算部4が行い、演算にて求められた光密度{stφij}μa2 のデータが記憶部6に格納される。
【0112】
(第3の処理)
次に、これらの光密度(stφij)μa1 と(stφij)μa2 の比を次式(23)に基づいて演算することにより、影響度マトリクスWG={stwgij}を求める。
【0113】
【数23】
Figure 0003771339
【0114】
尚、第1,第2の処理において、光密度(stφij)μa1 と(stφij)μa2 は各格子点における光密度φに夫々求められるので、影響度マトリクスWG={stwgij}も各格子点における光密度φに対応して求められる。
【0115】
そして、これらの影響度マトリクスWG={stwgij}も記憶部6に格納される。
【0116】
(第4の処理)
次に、実際に被検体Mに対して投射光Hs の位置sをずらしながら投光し、そのときの検出光Qt を実測することにより、投射光に対する検出光の対s,tの光密度(stφij)を得る。この実測された光密度(stφij)とシミュレーションによって求められている光密度(stφij)μa1 について次式(24)の演算を行うことにより、相対比の行列B={sφt}を求め、更に、前記のシミュレーションにより求められた影響マトリクスWGと行列Bについて次式(25)の演算を行うことにより、被検体M内の各有限体積の吸収係数の分布行列X={xij}を演算する。
【0117】
【数24】
Figure 0003771339
【0118】
【数25】
Figure 0003771339
このように、有限体積法によっても、断層画像を構築するための各有限体積の吸収係数を求めることができ、この分布行列Xに基づいて被検体Mの断層画像を再生することができる。
【0119】
尚、定常光(連続光)を用いる場合には、前記式(21),(22-a)〜(22-i)の代わりに、次式(26),(27-a)〜(27-g)を用いて上記第1の処理ないし第4の処理を行うことによって、断層画像を形成するための分布行列X={xij}を求めることができる。
【0120】
【数26】
Figure 0003771339
【0121】
【数27】
Figure 0003771339
又、以上に説明した処理方法では、全ての有限体積が所定の体積Δvであるとして断層画像を求める場合を説明したが、被検体Mと外部環境とが接する境界部分等の分布行列Xを求めるためには、図10(x座標を代表して示す)のように半有限体積を考え、前記式(21),(22-a)〜(22-i)の代わりに、次式(28-a)〜(28-d)を用いて、前記第1,第2の処理を行う。
【0122】
【数28】
Figure 0003771339
【0123】
尚、ここで、Iは定義格子点、BとEはx座標における隣接格子点を意味する。また、光流束qB の光伝達率をh、周囲光密度をφf 、境界面Bでの光密度をφB とすれば、qB=h(φf−φB)となり、b=hφf 、aB =aI +μa Δx+hとなる。
【0124】
また、第4の実施の形態において、第2の実施の形態と同様に、投射光Hs にパルス光を用いても断層画像を再構成することができる。
【0125】
【発明の効果】
以上に説明したように本発明によれば、従来のモンテカルロ法によらず、差分法等を適用した所定アルゴリズムに基づいて被検体の各体積要素の特徴抽出を行い、その特徴抽出データに基づいて断層画像を再構成するので、演算速度が速く高速の画像再構成が可能である。
【0126】
この結果、例えば生体の断層画像を再構成するための医療用光CT装置等の実用化等に貢献することができる優れた効果を発揮する。
【図面の簡単な説明】
【図1】実施の形態に係る光CT装置の構成を示す構成説明図である。
【図2】画像再構成の原理を説明するため、ディスプレイ上に表示した中間階調画像を示した図である。
【図3】光CT装置の動作を処理過程に沿って示すフローチャートである。
【図4】画像再構成の一例を説明するため、ディスプレイ上に表示した中間階調画像を示した図である。
【図5】パルス光を投射光に適用した場合に、差分法でシミュレーションして得られる検出光の波形を示すグラフである。
【図6】パルス光を投射光に適用した場合に、モンテカルロ法でシミュレーションして得られる検出光の波形を示すグラフである。
【図7】第3の実施の形態の処理行程を説明するためのフローチャートである。
【図8】第4の実施の形態の原理を説明するための説明図である。
【図9】第4の実施の形態の原理を更に説明するための説明図である。
【図10】第4の実施の形態の原理を更に説明するための説明図である。
【符号の説明】
2…制御部、4…演算部、6…記憶部、8…表示部、10…光源、12…光検出部、SWin,SWout …光スイッチ、F1〜Fn、B1〜Bn…光ファイバ、M…被検体、H1〜Hn…投射光、Q1〜Qn…検出光。

Claims (10)

  1. 被検体の複数の部位に個別排他的に投射光を投射する光源手段と、
    前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、
    前記被検体を微少な複数個の体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第1の演算手段と、
    前記第1の演算手段で想定した前記集合体モデルの特定の体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第2の演算手段と、
    前記第1の演算手段で求められた前記光密度と前記第2の演算手段で求められた光密度との比を前記体積要素の影響度として演算して、前記複数の体積要素の配列に対応する影響度の行列を求める第3の演算手段と、
    前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第1の演算手段で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める第4の演算手段と、
    前記第3の演算手段により求められた影響度の行列と前記第4の演算手段により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める第5の演算手段と、
    を具備することを特徴とする光CT装置。
  2. 被検体の複数の部位に個別排他的に投射光を投射する光源手段と、
    前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、
    前記被検体を複数個の要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第1の演算手段と、
    前記第1の演算手段で想定した前記集合体モデルの特定の要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第2の演算手段と、
    前記第1の演算手段で求められた光密度と前記第2の演算手段で求められた光密度との比を前記要素の影響度として演算して、前記複数の要素の配列に対応する影響度の行列を求める第3の演算手段と、
    前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第1の演算手段で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める第4の演算手段と、
    前記第3の演算手段により求められた影響度の行列と前記第4の演算手段により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める第5の演算手段と、
    を具備することを特徴とする光CT装置。
  3. 被検体の複数の部位に個別排他的に投射光を投射する光源手段と、
    前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する光検出手段と、
    前記被検体を複数の有限体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第1の演算手段と、
    前記第1の演算手段で想定した前記集合体モデルの特定の有限体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第2の演算手段と、
    前記第1の演算手段で求められた光密度と前記第2の演算手段で求められた光密度との比を前記有限体積要素の影響度として演算して、前記複数の有限体積要素の配列に対応する影響度の行列を求める第3の演算手段と、
    前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と、前記第1の演算手段で求められた光強度との相対比を演算し、前記複数の有限体積要素の配列に対応する相対比の行列を求める第4の演算手段と、
    前記第3の演算手段により求められた影響度の行列と前記第4の演算手段により求められた相対比の行列とのマトリクス演算により、前記各有限体積要素の特徴データを求める第5の演算手段と、
    を具備することを特徴とする光CT装置。
  4. 前記光源手段は、連続光を前記投射光として出射することを特徴とする請求項1乃至請求項3のいずれか1項に記載の光CT装置。
  5. 前記光源手段は、パルス光を前記投射光として出射し、前記光検出手段は、前記透過光を時間分解計測することを特徴とする請求項1乃至請求項3のいずれか1項にに記載の光CT装置。
  6. 光源手段により、被検体の複数の部位に個別排他的に投射光を投射すると共に、光検出手段により、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する第1の行程と、
    前記被検体を微少な複数個の体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第2の行程と、
    前記第2の行程で想定した前記集合体モデルの特定の体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく差分法により前記体積要素の光密度を演算する第3の行程と、
    前記第2の行程で求められた前記光密度と前記第3の行程で求められた光密度との比を前記体積要素の影響度として演算して、前記複数の体積要素の配列に対応する影響度の行列を求める演算をする第4の行程と、
    前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第1の演算手段で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める演算をする第5の行程と、
    前記第4の行程により求められた影響度の行列と前記第5の行程により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める演算をする第6の行程と、
    を備えることを特徴とする光CTによる画像再構成方法。
  7. 光源手段により、被検体の複数の部位に個別排他的に投射光を投射すると共に、光検出手段により、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する第1の行程と、
    前記被検体を複数個の要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第2の行程と、
    前記第2の行程で想定した前記集合体モデルの特定の要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限要素法により前記要素の光密度を演算する第3の行程と、
    前記第2の行程で求められた光密度と前記第3の行程で求められた光密度との比を前記要素の影響度として演算して、前記複数の要素の配列に対応する影響度の行列を求める第4の行程と、
    前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と前記第2の行程で求められた光強度との相対比を演算し、前記複数の体積要素の配列に対応する相対比の行列を求める第5の行程と、
    前記第4の行程により求められた影響度の行列と前記第5の行程により求められた相対比の行列とのマトリクス演算により、前記各体積要素の特徴データを求める第6の行程と、
    を備えることを特徴とする光CTによる画像再構成方法。
  8. 光源手段により、被検体の複数の部位に個別排他的に投射光を投射すると共に、光検出手段により、前記投射光が被検体内を透過して成る透過光を前記被検体の他の複数の部位から個別排他的に検出する第1の行程と、
    前記被検体を複数の有限体積要素に分割された集合体モデルとみなすと共に、被検体を一様な特定吸収係数又は吸収係数ゼロを有するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第2の行程と、
    前記第2の行程で想定した前記集合体モデルの特定の有限体積要素に他の光吸収体が存在するとみなして、前記光源手段及び前記光検出手段による測定条件と等価な条件下で前記集合体モデルを測定するとした場合の光拡散方程式に基づく有限体積法により前記有限体積要素の光密度を演算する第3の行程と、
    前記第2の行程で求められた光密度と前記第3の行程で求められた光密度との比を前記有限体積要素の影響度として演算して、前記複数の有限体積要素の配列に対応する影響度の行列を求める第4の行程と、
    前記光源手段及び前記光検出手段を用いて実測された前記被検体内を透過した透過光の光強度と、前記第2の演算手段で求められた光強度との相対比を演算し、前記複数の有限体積要素の配列に対応する相対比の行列を求める第5の行程と、
    前記第4の行程により求められた影響度の行列と前記第5の行程により求められた相対比の行列とのマトリクス演算により、前記各有限体積要素の特徴データを求める第6の行程と、
    を備えることを特徴とする光CTによる画像再構成方法。
  9. 前記光源手段により、連続光を前記投射光として出射することを特徴とする請求項6乃至請求項8のいずれか1項に記載の光CTによる画像再構成方法。
  10. 前記光源手段により、パルス光を前記投射光として出射し、前記光検出手段により、前記透過光を時間分解計測することを特徴とする請求項6乃至請求項8のいずれか1項に記載の光CTによる画像再構成方法。
JP00781097A 1996-01-18 1997-01-20 光ct装置及び光ctによる画像再構成方法 Expired - Lifetime JP3771339B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP00781097A JP3771339B2 (ja) 1996-01-18 1997-01-20 光ct装置及び光ctによる画像再構成方法

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
JP8-6619 1996-01-18
JP661996 1996-01-18
JP00781097A JP3771339B2 (ja) 1996-01-18 1997-01-20 光ct装置及び光ctによる画像再構成方法

Publications (2)

Publication Number Publication Date
JPH09257694A JPH09257694A (ja) 1997-10-03
JP3771339B2 true JP3771339B2 (ja) 2006-04-26

Family

ID=26340808

Family Applications (1)

Application Number Title Priority Date Filing Date
JP00781097A Expired - Lifetime JP3771339B2 (ja) 1996-01-18 1997-01-20 光ct装置及び光ctによる画像再構成方法

Country Status (1)

Country Link
JP (1) JP3771339B2 (ja)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2001020546A2 (en) * 1999-09-14 2001-03-22 The Research Foundation Of State University Of New York Imaging of scattering media using relative detector values
JP4364995B2 (ja) 2000-03-21 2009-11-18 浜松ホトニクス株式会社 散乱吸収体内部の光路分布計算方法
JP5644069B2 (ja) * 2009-07-07 2014-12-24 トヨタ自動車株式会社 塗膜反射率推定方法
CN113724351B (zh) * 2021-08-24 2023-12-01 南方医科大学 一种光声图像衰减校正方法

Also Published As

Publication number Publication date
JPH09257694A (ja) 1997-10-03

Similar Documents

Publication Publication Date Title
US5835617A (en) Optical computer tomographic apparatus and image reconstruction method using optical computer tomography
JP4271040B2 (ja) 実時間の光学的トモグラフィーの正規化された差方法の変更
Ren et al. Frequency domain optical tomography based on the equation of radiative transfer
US5070455A (en) Imaging system and method using scattered and diffused radiation
Gkioulekas et al. An evaluation of computational imaging techniques for heterogeneous inverse scattering
JP3310782B2 (ja) 吸収物質濃度の空間分布画像化装置
US5625458A (en) Method and system for imaging objects in turbid media using diffusive fermat photons
US8614707B2 (en) 3D and real time electrical capacitance volume-tomography sensor design and image reconstruction
Barbour et al. Imaging of subsurface regions of random media by remote sensing
JP3456721B2 (ja) 光音響コンピュータトモグラフィ装置
Rajendran et al. PIV/BOS synthetic image generation in variable density environments for error analysis and experiment design
JP2003532873A (ja) 不透明媒体における光学式コンピューター断層撮影
CN108627272A (zh) 一种基于四角度激光吸收光谱的二维温度分布重建方法
Wang et al. Imaging of scattering media by diffusion tomography: an iterative perturbation approach
JP3771339B2 (ja) 光ct装置及び光ctによる画像再構成方法
Allmaras et al. Reconstructions in ultrasound modulated optical tomography
US5758653A (en) Simultaneous absorption and diffusion imaging system and method using direct reconstruction of scattered radiation
Gong et al. Rapid simulation of X-ray scatter measurements for threat detection via GPU-based ray-tracing
KR20040094409A (ko) 산란 매체의 이미징에서 상호 파라미터 누화를 최소화하기위한 정규화 제한 알고리즘
US5746211A (en) Absorption imaging system and method using direct reconstruction of scattered radiation
WO2005110239A1 (ja) 画像再構成方法
Roggemann et al. Sensing three-dimensional index-of-refraction variations by means of optical wavefront sensor measurements and tomographic reconstruction
JP3836941B2 (ja) 光ct装置及び画像再構成方法
JP3730927B2 (ja) 吸収物質濃度空間分布画像化装置
Zalev et al. Fast 3-D opto-acoustic simulation for linear array with rectangular elements

Legal Events

Date Code Title Description
A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20060201

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: 20060207

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20060209

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: 20100217

Year of fee payment: 4

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20100217

Year of fee payment: 4

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20110217

Year of fee payment: 5

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20110217

Year of fee payment: 5

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120217

Year of fee payment: 6

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120217

Year of fee payment: 6

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130217

Year of fee payment: 7

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130217

Year of fee payment: 7

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20140217

Year of fee payment: 8

EXPY Cancellation because of completion of term