JP4304710B2 - 光測定装置 - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明は光測定装置に関し、主に、生体の脳機能診断や循環器系障害診断等の医療分野に適用して光診断画像装置とすることができる。
【０００２】
【従来の技術】
生体内部を生体に影響を与えずに測定する装置及び方法が医療分野で望まれており、可視から近赤外領域の波長の光を生体に照射し、照射位置から離れた位置から放出される光によって生体内部を測定する装置が提案されている。
光測定装置を医療分野に適用した光画像診断装置では、頭部や循環器等の各測定対象の酸素代謝の活性化状態を画像表示することが有効的である。
光測定装置は主に散乱吸収体である生体を被測定対象とし、被測定対象の被測定表面に複数の送光点・受光点を配置する。送光点・受光点の配置は、１つの送光点と複数の受光点および複数の送光点と１つの検出点の場合も含み、様々な送光点−受光点の組み合わせについてその散乱透過光を検出する。
【０００３】
図２４は従来の光測定装置の概要を説明する図である。図２４において、被検体には、レーザ光源（λ１〜λ３）からの光を照射する送光点と、被検体を散乱透過してきた光を検出する受光点とが６点ずつ配置され、受光点からの検出光を検出する検出器と検出信号を増幅する増幅器と検出信号をデジタル信号に変換するＡ／Ｄ器とデータ処理部と画像表示部を備える。図２５は、この送受光点の配置関係を示している。
一方、生体内部には酸素代謝機能により光スペクトル変化を起こす生体色素（例えば、ヘモグロビン、ミオグロビン、チトクロームａａ₃など）が存在し、これら生体色素に由来する光吸収（あるいは光吸収変化）を測定することで生体内部の酸素代謝に関する情報を得ることができる。例えば、７００ｎｍ〜９００ｎｍの複数波長における吸光度（あるいは吸収度変化）から酸素化／脱酸素化ヘモグロビン量(あるいは変化)に関する情報を得る装置（例えば酸素モニタ）が開発されている。
【０００４】
また、生体内に色素等の光造影剤を注入し、その濃度変化を光で追跡して生体内部の情報を得ることも提案されている。光測定装置及び光測定装置を用いた光画像診断装置は、複数の送受光を生体表面に配置して測定した吸光度情報から、生体内部の情報を二次元分布・三次元分布として画像表示し、該画像を用いて診断を行う。現在、実現されている光測定装置あるいは光画像診断装置は、例えば酸素化／脱酸素化ヘモグロビン変化情報等の生体内部情報を、図２５に示した送受光点平面に表示し、内部情報を生体表面に投影させた表面マッピング画像（トポグラフィ）として表している。
【０００５】
光測定装置を用いて得た光画像診断装置の臨床応用としては、人頭部に適用した場合の脳機能診断や局所的脳内出血診断など、腹部測定時の臓器機能診断などを挙げることができる。特に、例えば四肢を運動させた時に起こる脳運動野の賦活状態といった測定では、画像として捉えることが非常に有効である。こうした画像診断においては、変化領域の位置精度あるいは空間分解能や定量的な精度も求められる。
変化領域の位置精度あるいは空間分解能や定量的な精度に対する要求に対して、従来の光測定装置の画像表示、例えば、特開平９−１９４０８号公報では、その測定吸光度（あるいは吸収度変化）、又は複数の波長により測定された吸光度（あるいは吸収度変化）から演算される酸素化／脱酸素化ヘモグロビン量を送受光点の中点にある画素に代入し、送受光点の中点でない画素については複数の測定吸光度（あるいは吸収度変化）を補間した値を代入することで二次元画像を表示している。
【０００６】
【発明が解決しようとする課題】
測定信号を単に画素に代入・補間する方法で得られる画像は、大まかな変化分布の傾向をみることができるため現在臨床分野で利用されつつあるが、測定上不便なところや得られる画像に対する不満が指摘されている。
従来技術のように、測定信号を単に画素に代入・補間することによって画像表示を行う場合には、以下のような問題がある。
【０００７】
（１）送受光間距離の異なる測定では、同一画像を構成するデータとして扱うことができないため、従来の光測定装置では、全ての送受光点の組み合わせが同じ距離になるようにプローブを配置する必要があり、送受光点の設置位置に制限があるという問題がある。
（２）被検体内部にある吸収（あるいは吸収変化）部分と送受光点との位置関係によって、得られる画像の値が異なって表示されるという問題がある。
（３）被検体の表面からの深さ方向に対する情報に拡張することができず、三次元画像への拡張が困難であるという問題がある。
【０００８】
第１の問題となるのは、送受光プローブの配置が不自由であるという問題である。従来技術による画像方法では、１つの画像を再構成するのに用いることができるデータは全て送受光点距離が等しいもののみに限られる。例えば、６照射点・６検出点の組み合わせで測定する場合、図２５に示すように各送光点・受光点が等間隔の千鳥格子になるように配列すれば、送光点（Ａ）-受光点[ａ]，送光点（Ａ）-受光点[ｃ]，送光点（Ｂ）-受光点[ａ]，‥‥，送光点（Ｆ）-受光点[ｆ]の１６通りの送受光点の組み合わせで得られるデータによって画像を得ることができる。この場合、仮に送光点（Ａ）-受光点[ｂ]，送光点（Ａ）-受光点[ｄ]，‥‥といった組み合わせで得られる測定データは、前記した１６データとは送受光距離が異なるため、組み合わせて画像化することはできない。
【０００９】
更に、図２６において、測定部位の一部に傷があるとか他の装置のプローブと干渉するなどの理由で送光点（Ｂ）が等間隔の位置に置けない場合、例え僅か図中の上側に（Ｂ)をずらすことによって配置できるとしても、送光点（Ｂ）−受光点[ａ]，送光点（Ｂ）−受光点[ｂ]，送光点（Ｂ）−受光点[ｃ]，送光点（Ｂ）−受光点[ｄ]の４つのデータは他の１２データと合せて画像化することができず、補間に対する影響もあって結果的に斜線部分は欠損した画像になる。
【００１０】
また、測定する部位が横長な形をしているとき、同じ６送光点・６受光点を図２７のように配置することもできるが、長さaとbは等しくしなければならず、送受光点の組み合わせは１７通りとなる。結局、この配置によって測定できる測定部位の形は横３:縦２（あるいは縦３:横２）のみに限定される。
こうした送受光点の配置問題は、生体を対象とした装置であることを考えると大きな制限になるばかりでなく、特に脳機能計測分野では他の装置との組み合わせた計測が望まれていることからも好ましいことではない。
【００１１】
第２に問題になるのが、生体内部にある吸収（あるいは吸収変化）部分が送受光点に対してどこにあるかによって、得られる吸収（あるいは吸収変化)部分の画像値が異なって表示されるという問題である。生体内部において、送受光点に対する吸収（あるいは吸収変化）部分の位置によって、画像表示される変化部分の大きさや変化量(画像値)が異なって表示され、生体内部で同じ変化が生じたとしても、送受光点プローブの配置位置が少しずれるだけで表示画像が異なり、観測者に異なったイメージを与えてしまっていることになる。
【００１２】
こうした再構成画像の定量的な問題は、異なる画像間の比較はもとより、同一画像内の異なる部分の定量的な比較すらも確実ではなく、例えば、画像の左上側の変化は中央部分の変化より大きいといった知見も不確かなものとなる。これは、内部変化を画像として表示しているメリットを著しく損ねてしまうことにもなる。
第３に問題になるのは、被検体の深さ方向に対する拡張性の問題である。これまでの表面マッピング画像は、表面からの奥行方向に関して変化を重ね合わせた画像として表示している。生体内部の変化が表面からどの程度の深さで起こっているのかは、臨床的にも有用なデータである。特に脳機能測定において、ヘモグロビンの酸素化／脱酸素化変化が皮膚表面に近いところで起こっているのか、それより深い位置にある脳表面で起こっているのかは重要な情報である。しかしながら、従来技術による画像方法では、例え送受光点を増やして得られるデータを多くしてもこうした要求に原理的に応えられない。
【００１３】
そこで、本発明は従来の課題を解決して、送受光点の配置位置を制限することなく任意の配置位置で画像表示を行うことを目的とし、また、被検体内部の吸収（あるいは吸収変化）部分と送受光点との位置関係に影響されない画像表示を行うことを目的とし、また、被検体の表面からの深さ方向に対する情報に拡張して三次元の画像表示を行うことを目的とする。
【００１４】
【課題を解決するための手段】
前記した従来技術の各問題は、従来技術、例えば特開平９−１９４０８号公報では、
（ａ）光が散乱吸収体である被検体内部を透過する際の感度分布（飛行距離分布）は、散乱により送受光点間で拡がりをもった形となっているのにも関わらず、送受光点の中点のみを飛行するとして近似し、画素値を算出する点。
（ｂ）送受光点の中点でない画素は、上記の近似では不定解になってしまうため、中点の画素値を補間することによって画像値を求める。そのため、中点である画素とそうでない画素との扱いが全く異なる点、に起因している。
【００１５】
そこで、本発明では
（Ａ）光の感度分布（飛行距離分布）を求め、この感度分布を利用して各画素値を求め画像を構成する。
（Ｂ）画像の各画素値を求める演算式において、求める画像の画素数に対して送受光点の組み合わせによる測定値が少ないことによる不定要素を補うために、求める画像の空間的特性を制約する条件を付加することで、送受光点の配置位置や、また、画素位置と送受光点との位置関係に影響されない画素値を求め、二次元または三次元の画像表示を行うことができる。
【００１６】
なお、感度分布は、拡散理論やMonte・Carloシミュレーションなどによって求めることができる。また、測定される吸光度（吸収度変化）は、内部の吸収（吸収変化）と内部を飛行した光の距離の1次結合で近似されることから、内部の平均飛行距離分布を光拡散理論などで求め、更に求める画素に空間的な特性を与えることで二次元画像化する。なお、前記1次結合による近似は、散乱系において内部吸収変化による光の飛行分布変化が無視できるような吸光度変化であるとしている。
【００１７】
本発明の第１の態様は、複数の送光点から被検体に対して光を照射し、複数の受光点で被検体を透過してきた光を受光する光送受光手段と、光送受光手段で検出された光信号の光強度信号を測定する光強度信号測定手段と、光強度信号に対する重みを被検体内の各位置毎に表す感度分布を、複数の送受光点中の送受光点の一組み合わせ毎に形成する感度分布形成手段と、光強度信号と感度分布を用いて被検体の内部状態を表す画像信号を演算する画像形成手段とを備えた構成とし、画像形成手段で求めた画像信号を用いて被検体の内部状態を画像表示する。
【００１８】
本発明の第２の態様は、第１の態様に加えて、画像形成手段が形成する画像領域の空間特性を制限する空間的特性制約条件を形成する空間的特性制約条件形成手段を備えた構成とする。
画像形成手段は、求める画像の画素数に対して送受光点の組み合わせによる測定値が少ない場合、この空間的特性制約条件を加えることによって不定要素を補って画像信号を求める。空間的特性制約条件は、画像領域中の隣接する画像信号間は滑らかに変化する条件とすることができる。
【００１９】
また、本発明の光強度信号測定手段が形成する光強度信号は種々の態様とすることができる。光強度信号の第１の態様は、各時点における検出信号の信号強度、あるいは吸光度又は吸光度変化とすることができ、第２の態様は、被検体にパルス光を照射した時間を基準として定まる所定時間範囲のゲート時間内における光信号を積算して得られるゲート時間内光強度信号、あるいは吸光度又は吸光度変化とすることができる。
【００２０】
また、光強度信号の第３の態様は、被検体内の平均的な吸収係数あるいは吸収係数変化とすることができる。光強度信号測定手段は、光信号から被検体内の平均的飛行距離を算出して平均光路長を求める平均光路長算出部と、光信号から吸光度を算出する吸光度算出部と、吸光度を平均光路長で除算して平均吸収係数を算出する除算部を備える。また、この態様では、感度分布形成手段は、送受光点における総飛行距離を算出する総飛行距離算出部と、前記感度分布を総飛行距離で除算する除算部を含み、画像演算手段は、平均的な吸収係数あるいは吸収係数変化と、総飛行距離で除算した感度分布を用いて画像信号を演算する。
【００２１】
また、感度分布形成手段は、被検体内の光学的特性と画像信号から求めた光学的特性とを比較して感度分布を補正する処理を繰り返すことによって、感度分布の精度を高めることができる。
本発明の光測定装置は、光強度信号測定を複数波長で行い、画像形成手段で複数波長について画像信号を求め、この画像信号を用いて波長間演算を行う構成、あるいは、光強度信号測定を複数波長で行い、複数波長の光強度信号を用いて波長間演算を行って得られた信号を光信号として画像信号を形成する構成とすることができ、これによって、例えば酸素化／脱酸素化ヘモグロビン等の生体に特有の成分に関する画像を表示することができる。
【００２２】
本発明の光測定装置は、二次元の感度分布を用いて二次元画像を表示したり、あるいは三次元の感度分布を用いて二次元画像若しくは三次元画像を表示することができる。
また、本発明の光測定装置の画像形成手段は、任意の画素における画像信号の時間変化、あるいは任意の複数画素における画像信号の平均値の時間変化を基準とし、他画素における画像信号の時間変化との相関を求め、この相関によって被検体の内部状態を画像表示することができる。
さらに、本発明の光測定装置は、被検体の内部状態の画像表示、あるいは被検体の内部状態の画像の画素時間変化に加えて、磁気共鳴画像、Ｘ線画像、超音波画像、脳磁気画像、脳発生画像の少なくとも一画像を同一画面上に重ねて表示することができる。
【００２３】
本発明の光測定装置によれば、
（１）送受光間距離の異なる測定データを用いて画像を構成することができ、被測定部位に対してプローブを自由に制限なく配置することができる。また、再構成画像に用いるメッシュあるいはボクセルをプローブ配置に関わりなく、自由に設定できる。
（２）被検体内部にある吸収部分（あるいは吸収変化部分）と送受光点との位置関係に依存しない画素値を得ることができ、位置関係の相違による画像の変化を低減して表示することができる。
（３）被検体の表面情報に深さ方向に対する情報を加えて、三次元画像に拡張することができる。
（４）画像演算で計算される値は、平均吸収係数等の、物理的に意味のある内部吸収のみに依存する値であるため、ターゲットとしている物質のモル吸光係数から平均濃度に換算することができる。従って、異なるプローブ配置で求めた２つの画像間で比較検討を行なうことができる。
（５）物理的に意味のある感度分布を用いていることによって、感度分布の選択的適用や繰り返し演算を行なうことで、再構成される画像の空間分解能や濃度の定量性を向上させることができる。
（６）物理的に意味のある画像の空間的制限を用いることによって、測定データ数などの状態により空間分解能を上げて再構成することができる。
という効果を奏することができる。
【００２４】
【発明の実施の形態】
以下、本発明の実施の形態を図を参照しながら詳細に説明する。
はじめに、本発明の光測定装置による画像構成の原理について説明する。本発明では、画像の再構成を行うために被検体を適当なメッシュあるいはボクセル（微小立体）に分割する。メッシュあるいはボクセルは、画像を構成する画像領域を形成する。なお、メッシュは二次元で分割する場合であり、ボクセル三次元で分割する場合である。なお、以下、二次元の場合について説明する。
図１は、被検体の送受光点及び各メッシュ点の関係を説明する図である。図１に示すように、被検体の測定範囲を適当なメッシュで区切り、該メッシュに必要数の送光点（白丸印）と受光点（黒丸印）を配する。送光点と受光点（黒丸印）の組み合わせは、各個数を掛け合わせた個数存在するが、図１ではｉ番目の送光点とｊ番目の受光点の組み合わせのみを示している。この送光点ｉと受光点ｊの組み合わせをｍ番目の測定データ（吸光度変化ΔＡｂｓ（ｍ））とする。
【００２５】
一方、各メッシュを各画素と対応付けると共に番号を付し、ｋ番目のメッシュ（画素）における吸収係数の変化をΔμ（ｋ）とすると、ｍ番目の測定データは、（１）式で表されるように、各メッシュにおける吸収係数変化Δμ（ｋ）と重みとの一次結合式で表すことができる。

ここでＷ(ｍ，ｋ)は、ｋ番目のメッシュがｍ番目の測定データに対して持つ重み(感度分布)である。また、Ｗ(ｋ，ｍ)は単位が距離になり、ｋ番目のメッシュがｍ番目の測定データに寄与する平均光路長であると言うことができる。
【００２６】
送光点数をI点，受光点数をJ点とすれば、測定数の組は最大I×J個まで可能であるが、実際上、送受光点の距離が遠すぎるものはＳ／Ｎ比が悪く計算に使えない。そのため、最大値I×J個より少ない適当な数Ｍ個だけの式（１）が得られる。この式の個数Ｍがメッシュの総数Ｋよりも多ければそのまま連立方程式として解けるが、通常ＭはＫより小さいため不定解となって解くことができない。
そこで、画像領域の空間特性を制限する空間的特性制約条件を付加することによって連立方程式の個数を増やし、連立方程式を解ける状態とすることができる。
【００２７】
空間的特性制約条件は、各メッシュが持つの値Δμ（ｋ）に対して制約条件の式をｎ個形成して式（１）に追加することで、式の個数を（ｍ＋ｎ）をＭ以上にすることができる。これを行列表示した連立方程式が式（２）である。
式の数を十分大きくした後、Ｍ＝ｍ＋ｎならばＷの逆行列を、Ｍ＞ｍ＋ｎならば最小２乗法で連立方程式を解き、Δμ（ｋ）の画像の式（３）を得る。
【００２８】
【数１】

【００２９】
なお、ａは測定データを示し、ｘはΔμ（ｋ）を示している。
被検体のヘモグロビン濃度について画像表示する場合には、複数波長で吸収係数Δμ（ｋ）を求め、この吸収係数Δμ（ｋ）に所定係数を掛けて加え合わせることによって算出する。
なお、式（３）を非線形な拘束条件のもとで解かない限り、線形関係が保たれる。従って、最初に複数波長の測定データからヘモグロビン濃度に関する値を計算した後、このヘモグロビン濃度に対して前記演算を施して画像化することもできる。
【００３０】
図２の簡単な例を用いて上述の拘束条件を具体的に説明する。ここでは、３×３の合計９個のメッシュについて２つの送光点Ａ，Ｂと２つの受光点ａ，ｂを配置した例を示している。Ａ−ａ、Ａ−ｂ、Ｂ−ａ、Ｂ−ｂの４つの送受光点の組み合わせから図２（ｂ）〜図２（ｅ）の１番目〜４番目の測定データΔＡ₁〜ΔＡ₄と９個のメッシュの吸収係数の変化Δμ₁〜Δμ₉に対して４つの式が成立する。
ΔＡ₁＝Ｗ_1,1Δμ₁＋Ｗ_1,2Δμ₂＋…＋Ｗ_1,9Δμ₉
ΔＡ₂＝Ｗ_2,1Δμ₁＋Ｗ_2,2Δμ₂＋…＋Ｗ_2,9Δμ₉
ΔＡ₃＝Ｗ_3,1Δμ₁＋Ｗ_3,2Δμ₂＋…＋Ｗ_3,9Δμ₉
ΔＡ₄＝Ｗ_4,1Δμ₁＋Ｗ_4,2Δμ₂＋…＋Ｗ_4,9Δμ₉
しかし、未知数である吸収係数Δμ（ｋ）の個数はメッシュ数と同数の９個であるため、少なくとも５個の等式を必要とする。このため、感度分布を用いて以下の式を形成する。
０＝Ｑ_1,1Δμ₁＋Ｑ_1,2Δμ₂＋…＋Ｑ_1,9Δμ₉
０＝Ｑ_2,1Δμ₁＋Ｑ_2,2Δμ₂＋…＋Ｑ_2,9Δμ₉
:
０＝Ｑ_5,1Δμ₁＋Ｑ_5,2Δμ₂＋…＋Ｑ_5,9Δμ₉
これら５つの式が空間的特性制約条件であり、測定データから求めた４つの式及び感度分布から求めた５つの式を組み合わせることによって、解を得るに必要数の９つの連立方程式を形成することができる。
【００３１】
次に、本発明の光測定装置の構成例について説明する。
図３は、本発明の光測定装置の第１の構成例を説明するための概略図である。
光測定装置１の第１の構成は、光送受光手段２と光強度信号測定手段３と感度分布形成手段４と画像形成手段５とを備える。
光送受光手段２は、複数の送光点から被検体に対して光を照射し、複数の受光点で被検体を透過してきた光を受光する。光強度信号測定手段３は、光送受光手段２で検出された光信号の光強度信号を測定し、光強度信号を画像演算手段５に送る。
【００３２】
感度分布形成手段４において、感度は被検体が光強度信号に対して寄与する重みの程度を表し、感度分布は被検体内の各位置毎の感度であり、複数の送受光点中の送受光点の一組み合わせに対して形成される。
画像演算手段５は、光強度信号測定手段３で光強度信号と感度分布形成手段４で形成した感度分布を用いて、被検体の内部状態を表す画像信号を演算する。画像信号を表示することによって、被検体の内部状態を画像表示することができる。
【００３３】
図４は、光測定装置の第１の構成による画像表示の手順を示すフローチャートである。はじめに、光送受光手段２によって測定データを求める。この測定データは、送光点と受光点の組み合わせの中から有意な測定データを選択する（ステップＳ１）。一方、例えば被検体の形状や被検体内部の光学的パラメータ分布を用いて、拡散方程式に基づく理論解や有限要素法及びモンテカルロシミュレーションにより感度分布を求める（ステップＳ２）。なお、このステップＳ２はステップＳ１の測定工程の前にあらかじめ計算しておいてもよい。
次に、感度分布の値と被検体の内部状態を示す値（例えば吸収係数変化μａ）との一次結合式を各測定データ毎に形成して連立方程式を形成し、この連立方程式を解いて内部状態を示す値を求める（ステップＳ３）。求めた内部状態を示す値あるいは該値から派生した値を各画素値として、画像表示する（ステップＳ４）。
【００３４】
図５は、本発明の光測定装置の第２の構成例を説明するための概略図である。
光測定装置１の第２の構成は、第１の構成例に空間的制約条件検出手段６を付加したものである。空間的制約条件検出手段６は、光強度信号と感度分布を用いた画像演算において連立方程式の式数を増やすために、画像領域の空間特性を制限する空間的特性制約条件を形成する。画像演算手段５は、光強度信号と感度分布と空間的特性制約条件を用いて、被検体の内部状態を表す画像信号を演算し、画像表示する。
【００３５】
図６は、光測定装置の第２の構成による画像表示の手順を示すフローチャートであり、図４のフローチャートのステップＳ１と同様の測定データを求めて選択する工程（ステップＳ１１）と、ステップＳ２と同様の感度分布を求める工程（ステップＳ１２）に加えて、空間的制約条件を求める工程を付加する（ステップＳ１３）。
次に、感度分布の値と被検体の内部状態を示す値（例えば吸収係数変化μａ）との一次結合式を各測定データ毎に形成して連立方程式を形成し、更に、内部状態を示す値について空間的制約条件を用いた方程式を形成し、連立方程式に付加する。この式を付加した連立方程式を解いて内部状態を示す値を求める（ステップＳ１４）。求めた内部状態を示す値あるいは該値から派生した値を各画素値として、画像表示する（ステップＳ１５）。
【００３６】
次に、本発明で用いる感度分布について説明する。
本発明による光測定によって被検体の内部状態を表す画像を再構成するためには、各画素が測定データに与える影響の度合を示す重み関数を感度分布として求めておく必要がある。
感度分布は、対象とする被検体（例えば生体）の形状および被検体内部の光学パラメータ（散乱係数μｓ，吸収係数μａ，非等方パラメータｇ，屈折率ｎ等)の分布を用い、拡散方程式に基づく理論解や有限要素法（ＦＥＭ）・Ｍｏｎｔｅ・Ｃａｒｌｏシミュレーションなどを用いて計算することができる。
被検体によっては簡易なモデルを仮定して求めることができる。例えば、無限平板でかつ光学パラメータが均一であるとして感度分布を求めることができる。
図７は点Ａを送光点とし点Ｂを受光点とするとき、点Ｐ（ｘ₁，ｙ₁，ｚ₁）の感度を求めるための概略図であり、送光点Ａを原点、送受光点Ａ−Ｂ間の距離をρとしたとき、感度δＪ（ρ）は拡散方程式に基づいて以下の式で表すことができる。
δＪ（ρ）＝（ｚ₁dＳ／８π²Ｄ）・（１／ｂ³）・（１＋（μａ／Ｄ）^1/2・ｂ）・
exp（−（μａ／Ｄ）^1/2・ｂ）[（１／ａ₁）・exp（−（μａ／Ｄ）^1/2・ａ₁）
−（１／ａ₂・exp（−（μａ／Ｄ）^1/2・ａ₂）] …（４）
なお、Ｄは拡散係数（＝１／３（μ＋μｓ´））、ｄＳは点Ｐの断面積、ａ₁はＡＰ間距離、ａ₂はＡ´Ｐ間距離、ｂはＰＢ間距離である。
【００３７】
図８に示すようなモデル（送受光間距離２５ｍｍ）を考えたとき、送受光点を含む深さ方向の平面（ｘ−ｚ平面，ｙ＝０ｍｍ）の測定データに対する感度分布は図９のようになる。なお、図９の縦軸は図８で設定した構成微小立体（ボクセル）内を光が飛行する平均的な距離（ｍｍ）である。これをz軸方向に対して積算することによって、深さ方向（z軸方向）を積算したｘ軸上の感度分布となる。更にｙ軸方向に走査すると、図１０（ａ）に示すような二次元の感度分布（重み関数）を求めることができる。ただし、図１０（ａ）中のバースケールは平均飛行距離の常与対数表示としている。なお、図１０（ｂ）は、例えば測定データを送受光点の中点に置く従来方式に従った場合を、本発明の感度分布に対応させて示している。
【００３８】
したがって、対象とする生体の形状および生体内部の光学パラメータ分布などの情報を用いて厳密に求めた感度分布に代えて、経験的に知られている生理学的情報によって被検体に近似したモデルを設定して求めた感度分布によっても、従来の方式より充分良い感度分布を与えることができる。
【００３９】
また、繰り返し演算によって、より精度の高い感度分布を求めることができる。図１１，１２は繰り返し演算による感度分布の算出構成図及びフローチャートである。図１１に示す構成は、第１，２の構成例に光学的特性測定（算出）手段７及びパラメータ比較手段８を付加したものである。
測定データを求めておくと共に（ステップＳ２１）、光学的特性測定（算出）手段７によって被検体の光学パラメータ分布を仮定して初期パラメータを設定しておく（ステップＳ２２）。感度分布形成手段４は、仮定した初期パラメータを使い、前記した拡散方程式の解や有限要素法によるシミュレーションによって感度分布（重み関数）を求める（ステップＳ２３）。画像演算手段５は、測定データ、感度分布、及び必要な場合には空間的制約条件形成手段６で形成した空間的特性制約を用いて、画像信号を求めて画像を再構成する（ステップＳ２４）。
【００４０】
次に、パラメータ比較手段８によって、得られた画像のパラメータと感度分布を計算したパラメータとを比較して、その差が充分に小さければ終了し、差が大きければパラメータを修正してステップＳ２３，２４により再度感度分布の計算をやり直し、収束するまで行なう（ステップＳ２５）。
この繰り返し演算によって、より精度の高い感度分布を求め、内部の光学パラメータ分布を反映させた感度分布による画像の再構成を行なう（ステップＳ２６）。
【００４１】
次に、二次元画像の空間的特性を制限する制約条件について、具体的な設定方法を説明する。
空間的特性制約条件は、画像中の隣接位置の解の連続性を拘束条件とするものであり、被検体の空間的な吸収係数変化が滑らかであり、各メッシュが持つ値がその近傍にあるメッシュの値に対して緩やかな連続している関係であることを利用するものである。この拘束条件によれば、得られる画像は滑らかな像となる。
【００４２】
一般に、滑らかな像は、元画像に高周波成分をカットする処理を施すことによって得ることができる。したがって、画像ｇ（ｘ，ｙ）が滑らかである場合には、二次元フーリエ面における画像Ｇ（ｕ，ｖ）が高周波成分を含まない。また、二次元画像のフィルタ処理において高周波成分をカットするには、高周波成分をカットするフィルタ関数Ｆ（ｕ，ｖ）を元画像Ｒ（ｕ，ｖ）に掛合わせる信号処理を行う。なお、同様な処理を二次元実空間で行なうには、二次元フーリエ面のフィルタ関数Ｆ（ｕ，ｖ）を逆二次元フーリエ変換したｆ（ｘ、ｙ）を元画像ｒ（ｘ，ｙ）に二次元コンボーリューションする。図１３（ａ）は高周波成分をカットするフィルタ関数の例を示し、図１３（ｂ）は二次元フーリエ面で表示したフィルタ関数であり、逆二次元フーリエ変換することによって図１３（ａ）の二次元実空間のフィルタ関数を得ることができる。
【００４３】
したがって、メッシュが持つ値が滑らかであるという拘束条件を満たす場合には、画像に対して高周波成分をカットする制約を付けて再構成を行なった場合、元画像及び処理後の画像は共に高周波成分がカットされたものになっているはずである。言い換えれば、再構成で得られる画像は、二次元フーリエ面上において図１３（ｂ）のフィルタ関数が掛かった状態にあり、その像に図１３（ｂ）のフィルタ関数を掛合わせても同じ画像になる。
【００４４】
このことは、二次元実空間上では、図１３（ａ）のフィルタ関数を二次元コンボリューションした画像は元の画像と同一とならなければならない。この拘束条件を画像に対する空間的制約条件とすることができ、離散的フィルタ関数をf(ｉ,j)(ただし、ｉ，ｊ＝‥‥,−２,−１,０,１,２,‥‥)とし、離散的画像をｒ(ｉ,j)(ただし、ｉ，ｊ＝‥‥,−２,−１,０,１,２,‥‥)とするとき、ｆとｒを二次元コンボリューションしたときの画像が元画像ｒ（ｉ，ｊ）になるという関係が成り立ち、全ての画素(ｉ,j)に対して次式（５）が成り立つことになる。
ｒ(ｉ,j)＝Σ_mΣ_nｒ(ｉ−ｍ,j−ｎ)・ｆ（ｍ,ｎ) …(５）
ただし、ｍ＝…，−２，−１，０，１，２，…、ｎ＝…，−２，−１，０，１，２，…とする。
【００４５】
画像に対して滑らかになるという空間的な制約条件は、上記した高周波成分をカットしたフィルタを掛けても同じ画像になる制約条件の他に、画像の高次微分が０になるという制約条件とすることもできる。
この空間的制約条件によって得た式を感度分布によって形成した式に付加することによって、前記式(２)を形成することができる。
【００４６】
前記説明では二次元の場合について示したが、本発明の光測定装置は三次元についても適用することができる。以下、図１４〜図１６を用いて三次元への適用について説明する。図１で示した平面メッシュを図１４のように立体的なメッシュ（ボクセル）とすることによって、三次元画像に拡張することができる。このとき、用いる感度分布は図９に示した３次元分布をｘ，ｙ，ｚの立体で算出して使用する。また、三次元空間の特性を制限するためのフィルタ関数としては、例えば図１５に示すような三次元フーリエ面における高周波をカットするような関数があり、これを三次元逆フーリエ変換した関数から式(５)と同様に設定し、三次元コンボリューション演算を行う。
【００４７】
得られた三次元画像データは、図１６（ａ）に示すように立体的な形状のある切り口を濃淡表示したり、図１６（ｂ）に示すように表面から深さ方向に層状に見立てて平面画面表示することもできる。特に、脳機能計測では大脳皮質表面の活性分布が求められるので、例えば表面から１ｃｍ〜１．５ｃｍにあたる層を表示することで、疑似的に大脳皮質表面マッピングすることができる。更に、図１６（ｃ）に示すように各深層マッピング画像を深さ方向に足し合わせて、二次元再構成画像と同様な表面マッピング画像を求めることができる。また、ある深さ以上の深層マッピング画像のみを積算対象としたマッピング画像を表示することも可能である。
以上のように、従来技術では求められなかった深さ方向に対する情報抽出が可能になるため、特に光脳機能診断で問題になっている皮膚下血液変化と大脳血液変化とを分離して議論できる点で新たなる臨床的見地が得られる可能性がある。
【００４８】
次に、本発明の光測定装置の他の実施形態について説明する。
他の実施形態の第１は、測定データをゲート時間内吸光度（吸光度変化）とする。
図１７に示すような光時間分解計測システムを用いることで、極短パルス光を被検体に照射したときに再放射される光を時間分解した波形（時間プロファイル）として計測することができる。光時間分解計測システムは、極短パルス光源と時間分解検出器、及び制御部、演算部、画像表示部を備え、極短パルス光源から発した極短パルス光を被検体に照射し、被検体から放出される光を時間変化を検出する。
【００４９】
観測される時間分解波形は、図１８に示すように極短パルスで入射された光が被検体内部の散乱により拡がりを生じている。そこで、この時間分解波形に対して特定の時間範囲（飛行距離）のみを使って吸光度を求めることによりゲート時間内吸光度を得ることができ、ゲート時間内吸光度の変化からゲート時間内吸光度変化を得ることができる。
このゲート時間内吸光度（吸光度変化）を測定データとして本考案に適用するには、前記式（２）における感度分布を対応するゲート時間内感度分布とする。このゲート時間内吸光度（吸光度変化）を用いることによって、１つの時間分解波形に対して設定する時間範囲（飛行距離）を変えることで、複数の測定データを得ることができ、異なる感度分布の測定データを数多く使用して演算することによって高分解・高精度の再構成画像を得ることができる。
【００５０】
他の実施形態の第２は、前記他の実施形態の第１と同様に光時間分解計測システムから得られる測定データを用いるものであり、吸収係数として平均的な吸収係数を算出して画像構成を行う。図１９に示す構成例は、平均的な吸収係数を算出する手段として、平均光路長算出手段３ａ、吸光度算出手段３ｂ、除算手段３ｃを備え、平均的な吸収係数に対応する感度分布を求める手段として、感度分布形成手段４ａ、総飛行距離算出手段４ｂ、除算手段４ｃを備える。
平均光路長算出手段３ａは、光時間分解計測システムで得られる時間分解波形を用いて、各送受光点における被検体内の平均的な飛行距離（平均光路長）を計算する。平均光路長は、時間分解波形の時間軸に対する1次モーメントとも呼ばれ、図２０において斜線部面積の重心点となる。除算手段３ｃは、吸光度算出手段３ｂで求めた吸光度（吸光度変化）を平均光路長で除算して、その送受光点における被検体内の平均的な吸収係数（吸光度変化）を算出する。
【００５１】
また、除算手段４ｃは、感度分布形成手段４ａで形成した感度分布を総飛行距離算出手段４ｂで求めた総飛行距離で除算して、平均的な吸収係数に対応する感度分布を求める。なお、総飛行距離は全メッシュにおける飛行距離の積算値によって求めることができる。平均吸収係数を用いた画像の再構成は、得られる画像の定量性が向上するという効果がある。
平均吸収係数から表面マッピング画像を得るための方法としては、演算で得られた平均吸収係数を送受光点の線上に配置する方法も考えられるが、本発明は、被検体内部の吸収（吸収変化）部分と送受光点との位置関係から影響を受けにくいという利点を有している。
【００５２】
なお、平均吸収係数を求める他の方法として、拡散方程式を時間分解波形にフィティングする方法（例えば特開平７−３２３０３３号公報）を適用することもできる。
また、平均吸収係数を時間ゲートを用いて算出することもできる。この場合には、図２１のようにゲート時間内の検出信号から得られるゲート時間内吸光度（吸光度変化）と、ゲート時間内平均光路長から算出されるゲート時間内平均吸収係数を用いて、画像を再構成する。
【００５３】
他の実施形態の第３は、感度分布を固定的ではなく選択的に求めるものであり、被検体からの測定データに基づいて感度分布を求める。
感度分布は、あらかじめ推定した生体内部の光学パラメータ分布を用いて算出の他に、拡散方程式を時間分解波形にフィティングすることによる送受光間毎の生体内部の平均的な光学パラメータの算出、空間的な光の減衰量を用いた光学パラメータの算出等、被検体からの測定データに基づいて感度分布を求める手法で考えられる。
【００５４】
選択的適用としては、画像再構成アルゴリズム内で得られた光学パラメータを用いて感度分布を計算する方法、あらかじめ複数の光学パラメータに対する感度分布を計算しておいて測定された値から選択及び補間を行う方法、推定した光学パラメータで得られた感度分布を測定された値により修正する方法等がある。 例えば、拡散方程式を時間分解波形にフィティングする方法では、送受光間における散乱係数μｓ´および吸収係数μａを求め、それらの値を用いて感度分布を算出する。
【００５５】
また、図１１，１２に示したような繰り返し演算により画像を求める方法では、演算で得られた内部の光学パラメータ分布を用いて感度分布を再計算あるいは選択して感度分布を算出する。
また、あらかじめ求めた感度分布を総飛行距離で除算して平均吸収係数に対する感度分布を求めることによって、感度分布を修正することができる。
内部の光学パラメータに関する情報を用いて感度分布を選択的に適用することによって、実際の被検体に対応した感度分布を求めることができ、再構成される画像の空間分解能や濃度の定量性を向上させることができる。
【００５６】
他の実施形態の第４は、空間的特性制約の選択的に適用するものである。本発明の空間的特性制約条件は、求める画像の画素数に対して送受光点の組み合わせによる測定値が少ない不定要素を補うものであり、送受光点の組み合わせ数（測定データ数）によって、その特性を選択的に適用することができる。
例えば、測定データ数が少ない程、空間的制約によってカットする高周波成分を増やした条件で再構成し、逆に測定データ数が多い程、空間的制約がカットする高周波成分を減らした条件で再構成することによって、測定データ数に応じた最適な空間分解能を得ることができる。
【００５７】
この空間的特性制約の選択的適用は、１画像における各画素に対して適用することも、あるいは、複数の画素を含む画像領域に対して適用することもできる。例えば、ある画素の近辺において、測定データ密度が低い領域と高い領域で異なる空間的特性制約を設定することもでき、実際の被検体に対応した画像を得ることができる。
他の実施形態の第５は、送受光間距離の選択的に適用するものである。
本発明の感度分布は送受光間距離をパラメータとして計算される分布であり、測定した送受光間距離を用いて計算される感度分布を用いることによって、異なる送受光間距離の測定データであっても、画像を再構成するためのデータとして同等に用いることができる。
【００５８】
次に、本発明によるコンピュータシミュレーション結果例を示す。
このシミュレーションでは送受光プローブを各８本ずつ、２５ｍｍピッチで配置した試料において内部の一部（１８．７５ｍｍ×１８．７５ｍｍ×６．００ｍｍ)が吸収変化を起こしたときを考える。図２２はモデル図を示し、図２３は本発明による構成画像と従来技術による構成画像を比較して示している。
図２３において、（２−ａ）及び（３−ａ）は、（１−ａ）に示すように送光点の真下において吸収変化した場合について、従来技術による画像および本発明による画像を示している。同様に、吸収変化の位置が一組の送受光点の中間にある場合（図２３（１−ｂ），（２−ｂ），（３−ｂ））と二組の送受光点の中間にある場合（図２３（１−ｃ），（２−ｃ），（３−ｃ））について示している。
【００５９】
画像の比較から、従来技術による画像は同じ吸収変化が生じた場合であっても、吸収変化の位置と送受光点との相対的な位置関係が変わるだけで再構成される吸収変化の画像が大きさ・濃淡の違いとして現れる。これに対して、本発明による画像は、従来技術による画像と比較して、比較的同じような濃淡で再構成される。
図２３（４−ａ），（４−ｂ），（４−ｃ）は各画像における吸収変化画像の中心における断面プロット図である。この図から、吸収変化の中心における値が従来技術による画像（トポグラフィで示す）では吸収変化の位置によって約３倍程度になるのに対して、本発明による画像（新アルゴリズムで示す）では約１．５倍程度と改善されていることが分かる。
【００６０】
更に本発明において吸収変化が非負の条件を用いて画像を再構成すると、図２３（５−ａ），（５−ｂ），（５−ｃ）となり、ほぼ与えた吸収変化の形と同じになり、吸収変化の中心における値も更に改善される。
この吸収変化が非負という条件は本発明のみで有効であり、従来技術に適用しても改善されない。
なお、図２３の画像は、表層から３ｍｍまでは吸収変化が起こらないとし、式（３）は重み付きの最小２乗法で計算している。
【００６１】
【発明の効果】
以上説明したように、本発明の光測定装置によれば、送受光点の配置位置を制限することなく任意の配置位置で画像表示を行うことができる。また、被検体内部の吸収（あるいは吸収変化）部分と送受光点との位置関係に影響されない画像表示を行うことができる。また、被検体の表面からの深さ方向に対する情報に拡張して三次元の画像表示を行うことができる。
【図面の簡単な説明】
【図１】本発明の被検体の送受光点及び各メッシュ点の関係を説明する図である。
【図２】本発明の被検体のメッシュ点及び送受光点例を説明する図である。
【図３】本発明の光測定装置の第１の構成例を説明するための概略図である。
【図４】本発明の光測定装置の第１の構成による画像表示の手順を示すフローチャートである。
【図５】本発明の光測定装置の第２の構成例を説明するための概略図である。
【図６】本発明の光測定装置の第２の構成による画像表示の手順を示すフローチャートである。
【図７】感度を求めるための概略図である。
【図８】感度分布を求めるモデル図である。
【図９】感度分布の一例を示す図である。
【図１０】二次元の感度分布の一例を示す図である。
【図１１】繰り返し演算による感度分布の算出構成図である。
【図１２】繰り返し演算による感度分布のフローチャートである。
【図１３】高周波成分をカットするフィルタ関数の例である。
【図１４】本発明による画像の立体的なメッシュ（ボクセル）を説明する図である。
【図１５】三次元フーリエ面における高周波をカットする関数例である。
【図１６】三次元画像の表示例である。
【図１７】光時間分解計測システムを説明する図である。
【図１８】ゲート時間内検出光信号の一例である。
【図１９】平均的な吸収係数を算出する構成例である。
【図２０】平均光路長による平均吸収係数を説明する検出信号図である。
【図２１】時間ゲートによる平均吸収係数を説明する検出信号図である。
【図２２】本発明によるコンピュータシミュレーションのモデル図である。
【図２３】本発明によるコンピュータシミュレーションによる構成画像と従来技術による構成画像を比較する図である。
【図２４】従来の光測定装置の概要を説明する図である。
【図２５】送受光点の配置関係を示す図である。
【図２６】送受光点の配置関係を示す図である。
【図２７】送受光点の配置関係を示す図である。
【符号の説明】
１…光測定装置、２…送受光手段、３…光強度信号測定手段、３ａ…平均光路長算出手段、３ｂ…吸光度算出手段、３ｃ，４ｃ…除算手段、４，４ａ…感度分布形成手段、４ｂ…総飛行距離算出手段、５…画像演算手段、６…空間的制約条件形成手段、７…光学的特性測定手段、８…パラメータ比較手段。

Claims (11)

1. 複数の送光点から被検体に対して光を照射し、被検体中を通った後再放出される光を複数の受光点で受光する光送受光手段と、
   前記光送受光手段で検出された光信号の光強度信号を測定する光強度信号測定手段と、
   光強度信号に対する重みを被検体内の各位置毎に表す感度分布を、前記複数の送受光点中の送受光点の一組み合わせ毎に形成する感度分布形成手段と、
   前記光強度信号と感度分布を用いて被検体の内部状態を表す画像信号を演算する画像形成手段と、
   前記画像形成手段が形成する画像において、送光点と受光点との組み合わせによる測定数の組が画像のメッシュの総数より小さいため不定解となって解くことができない条件の下で、各メッシュが持つ値がその近傍にあるメッシュの値に対して緩やかな連続している関係であることを利用して空間的特性制約条件を形成する空間的特性制約条件形成手段とを備え、
   前記画像形成手段は、光強度信号と感度分布と空間的特性制約条件を用いて、被検体の内部状態を表す画像信号を演算することを特徴とする、光測定装置。
2. 前記空間的特性制約条件は、画像領域中の隣接する画像信号間は滑らかに変化する条件であることを特徴とする、請求項１記載の光測定装置。
3. 前記光強度信号測定手段は、被検体にパルス光を照射した時間を基準として定まる所定時間範囲のゲート時間内における光信号を積算して得られるゲート時間内光強度信号を光強度信号とすることを特徴とする、請求項１又は２記載の光測定装置。
4. 前記光強度信号測定手段は、光信号から被検体内の平均的飛行距離を算出して平均光路長を求める平均光路長算出部と、光信号から吸光度を算出する吸光度算出部と、吸光度を平均光路長で除算して平均吸収係数を算出する除算部を備え、
   前記感度分布形成手段は、送受光点における総飛行距離を算出する総飛行距離算出部と、前記感度分布を総飛行距離で除算する除算部を含むことを特徴とする、請求項１記載の光測定装置。
5. 前記感度分布形成手段は、被検体内の光学的特性と画像信号から求めた光学的特性とを比較して感度分布を補正することを特徴とする、請求項１記載の光測定装置。
6. 光強度信号測定を複数波長で行い、
   画像形成手段で得られた複数波長の画像信号を用いて波長間演算を行い、
   得られた信号を用いて画像表示することを特徴とする、請求項１記載の光測定装置。
7. 光強度信号測定を複数波長で行い、該複数波長の光強度信号を用いて波長間演算を行い、得られた信号を前記光信号とすることを特徴とする、請求項１記載の光測定装置。
8. 二次元感度分布による二次元画像の表示、あるいは三次元感度分布による二次元画像若しくは三次元画像の表示を行うことを特徴とする、請求項１から請求項７のいずれかに記載の光測定装置。
9. 前記画像形成手段は、任意の画素における画像信号の時間変化、あるいは任意の複数画素における画像信号の平均値の時間変化を基準とし、他画素における画像信号の時間変化との相関を求め、該相関によって被検体の内部状態を画像表示することを特徴とする、請求項１から請求項８のいずれかに記載の光測定装置。
10. 前記表示は、前記被検体の内部状態の画像表示に、磁気共鳴画像、Ｘ線画像、超音波画像、脳磁気画像、脳発生画像の少なくとも一画像を同一画面上に重ねて表示することを特徴とする、請求項１から請求項９のいずれかに記載の光測定装置。
11. 前記表示は、前記被検体の内部状態の画像の画素時間変化と、磁気共鳴画像、Ｘ線画像、超音波画像、脳磁気画像、脳発生画像の少なくとも一画像の画素時間変化との相関を求め、該相関によって被検体の内部状態を画像表示することを特徴とする、請求項１から請求項１０のいずれかに記載の光測定装置。

Images