(57)【要約】【課題】  少ない光検出位置からの光ＣＴ計測の計測データから、精度良く光ＣＴ画像を再構成することが可能な画像再構成方法を提供する。【解決手段】  被検体に対して、所定の光照射位置からパルス光を照射し、この被検体内を透過して再出射されるパルス光を、光照射位置とは異なる複数の光検出位置で検出して、その検出光の時間分解波形から、被検体の内部情報を示す光ＣＴ画像を再構成する。光検出位置近傍では、被検体自身または被検体との間に配置される緩衝材により被検体から再出射されるパルス光の拡散を確保している。この拡散性を前提として、各検出位置における時間分解波形から時間的空間的補間により検出位置間における時間分解波形を求め、検出した時間分解波形と前記補間により求めた時間分解波形を基にして被検体の内部情報を算出する。【選択図】  図４
2005328916.tif 000002
【特許請求の範囲】
【請求項１】
被検体に対して、所定の光照射位置からパルス光を照射し、前記被検体内を透過して前記被検体から再出射されるパルス光を、前記光照射位置とは異なる複数の光検出位置で検出し、検出光の時間分解波形から、前記被検体の内部情報を示す光ＣＴ画像を再構成する画像再構成方法であって、
前記光検出位置近傍では、前記被検体自身または前記被検体との間に配置される緩衝材により被検体から再出射されるパルス光の拡散を確保しており、
この拡散性を前提として、各検出位置における時間分解波形から時間的空間的補間により検出位置間における時間分解波形を求め、
検出した時間分解波形と前記補間により求めた時間分解波形を基にして前記被検体の内部情報を算出することを特徴とする画像再構成方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、生体などの被検体に対してパルス光を照射し、被検体内をパルス光が透過された透過光を検出することによって、被検体の内部情報を示す光ＣＴ画像を取得して再構成する画像再構成方法に関するものである。
【背景技術】
【０００２】
現在、医療用の画像診断装置などにおいて、非侵襲で被検体の内部情報を計測する計測装置として、Ｘ線ＣＴ、超音波ＣＴ、ＭＲＩなどが用いられている。これに対して、近年、近赤外光が生体組織に対して高い透過性を有すること、生体組織中の酸素濃度を計測し得ること、Ｘ線などと比較して安全性が高いことなどから、非侵襲の計測方法として光ＣＴの開発と利用が進められている（例えば、非特許文献１参照）。
【０００３】
光ＣＴ計測においては、計測対象となる被検体に対してパルス光を照射する。そして、照射したパルス光が被検体内を透過して被検体の外部へと出力される透過光を検出し、その透過光の時間分解波形から被検体の内部情報を取得するとともに、取得された内部情報から被検体の光ＣＴ画像を再構成する（例えば、特許文献１参照）。また、連続光（ＣＷ光）を用いた光ＣＴ計測でも、同様に、被検体の内部情報を示す光ＣＴ画像を得ることができる。
【非特許文献１】上田之雄著「近赤外光を用いた生体計測」光技術コンタクトVol.35,No.7、1999年、397頁〜404頁
【特許文献１】特開平７−２３９２９９号公報
【発明の開示】
【発明が解決しようとする課題】
【０００４】
上述した光ＣＴ計測は、例えば、乳房内での腫瘍の有無あるいは乳房内組織の吸収分布などを画像化する非侵襲の乳がん診断装置（マンモグラフィー装置）などに用いることができる。このように、光ＣＴ計測を生体の診断等に用いる場合、光ＣＴ画像である生体の断層画像を、光ＣＴ計測の計測データから高精度かつ高速で再構成することが求められている。しかしながら、光ＣＴ装置において従来用いられている画像再構成方法では、充分な精度が得られていない。
【０００５】
すなわち、光ＣＴ計測では、通常、計測対象となる被検体が光の拡散距離などと比べて大きい。また、光ＣＴ画像での解像度を上げようとすれば、パルス光を照射する光源や、透過光を検出する光検出器などの個数を増やす必要がある。
【０００６】
しかしながら、光検出位置を増やすにも物理的な制限があり難しい。なぜなら、マンモグラフィー装置等では、各検出位置に到達する赤外線の光量は決して大きくなく、その時間分解波形を測定するには検出器を検出位置に直接配置するのではなく、光ファイバを用いて外部の光電子増倍管等へと導き、増幅した信号の時間変化を測定する手法が一般的である。このため、物理的な制約から検出位置同士の間隔を近づけることが困難であり、また、検出位置を増加させると、検出部全体の重量増大や大型化を招くため、測定の自由度が低下するという問題もある。
【０００７】
以上の問題点に鑑み、本発明は、少ない光検出位置からの光ＣＴ計測の計測データから、精度良く光ＣＴ画像を再構成することが可能な画像再構成方法を提供することを課題とする。
【課題を解決するための手段】
【０００８】
上記課題を解決するため、本発明にかかる画像再構成方法は、被検体に対して、所定の光照射位置からパルス光を照射し、この被検体内を透過して被検体から再出射されるパルス光を、光照射位置とは異なる複数の光検出位置で検出し、検出光の時間分解波形から、被検体の内部情報を示す光ＣＴ画像を再構成する画像再構成方法であって、光検出位置近傍では、被検体自身または被検体との間に配置される緩衝材により被検体から再出射されるパルス光の拡散を確保しており、この拡散性を前提として、各検出位置における時間分解波形から時間的空間的補間により検出位置間における時間分解波形を求め、検出した時間分解波形と補間により求めた時間分解波形を基にして被検体の内部情報を算出することを特徴とする。
【０００９】
検出位置周囲が拡散性の材質で囲まれている場合には、光検出位置へと到達する光においては、照射位置から直進して光検出位置へと到達した光（直進成分）の比率はほとんどなく、散乱を繰り返した結果、到達してきた光（散乱成分）がそのほとんどを占めることになる（図１２参照）。このように散乱が主体となる結果、各検出位置へ到達した光は、光検出位置と光照射位置を結ぶ直線を通過してきたものに限らず、この直線を包含する広い領域を通過・拡散してきた光成分を含むことになるから、結果的に、これらの広い領域の透過特性・散乱特性情報を含むことになる。このため、実際に光検出器が配置されている光検出位置間に到達する光の時間分解波形を時間的空間的補間により精度良く求めることが可能となる。
【発明の効果】
【００１０】
本発明によれば、物理的に検出点を増やすことなく、異なる位置の時間分解波形を取得することができるので、得られる光ＣＴ画像の解像度を向上させることが可能となる。そして、物理的な検出点を増やす必要がないため、検出部の質量増大や大型化を招くことがなく、測定の自由度を確保することができる。
【発明を実施するための最良の形態】
【００１１】
以下、添付図面を参照して本発明の好適な実施の形態について詳細に説明する。説明の理解を容易にするため、各図面において同一の構成要素に対しては可能な限り同一の参照番号を附し、重複する説明は省略する。なお、図面の寸法比率は、実際のものと必ずしも一致していない。
【００１２】
図１は、本発明による光ＣＴ装置１００の一実施形態の構成を概略的に示すブロック図である。図１に示す光ＣＴ装置１００は、生体などの被検体Ｍを計測対象として光ＣＴ計測を行うものである。本光ＣＴ装置１００は、被検体Ｍに対してパルス光を照射する光照射手段１と、被検体Ｍ内をパルス光が透過された透過光を検出する光検出手段２と、透過光を計測する時間分解波形計測装置３と、装置各部の動作制御等を行う制御装置４とを備える。
【００１３】
ここで、以下においては、被検体Ｍとして略正方形状の２次元ファントムを想定して、対応する光ＣＴ装置の構成、及び被検体Ｍの光ＣＴ画像の再構成方法について説明する。ただし、以下に示す光ＣＴ装置及び画像再構成方法を適用することが可能な被検体は、２次元ファントムには限られない。なお、この被検体は少なくとも周囲の領域（検出点付近だけでもよい。）が図１２に示されるように良拡散性の物質で覆われている必要がある。
【００１４】
光照射手段１は、パルス光供給部１０と、入力用光ガイド１２と、パルス光入力部１３とから構成されている。パルス光供給部１０には、４辺で囲まれた正方形状の２次元ファントムである被検体Ｍに対応して、所定波長及び所定時間幅のパルス光を出力するパルス光源１１が設けられている。パルス光源としては、例えばパルスレーザなどが用いられる。
【００１５】
このパルス光源１１に対してパルス光を供給するパルス光入力部１３が設けられている。パルス光源１１と、パルス光入力部１３とは、入力用光ガイド１２により光学的に接続されている。この入力用光ガイド１２としては、赤外光用光ファイバを用いることができる。
【００１６】
一方、光検出手段２は、透過光検出部２０と、出力用光ガイド２２と、透過光出力部２３とから構成されている。被検体Ｍをはさんで、上述した光照射位置に対向する位置には、所定の配列間隔で配列されたＮ個（Ｎは複数）の光検出位置が設定されている。そして、各光検出位置には、Ｎ個の透過光出力部２３１〜２３Ｎが設けられている。透過光出力部２３１〜２３Ｎのそれぞれは対応するＮ個の光検出器２１１〜２１ＮとＮ本の光ガイド２２１〜２２Ｎにより光学的に接続されている。これらの光ガイド２２１〜２２Ｎもまた、赤外光用光ファイバを用いることができる。また、光検出器２１１〜２１Ｎとしては、微弱な赤外線を精度良く検出することが可能な光電子増倍管等を用いることが好ましい。これにより、被検体Ｍを透過・拡散して光検出位置から外部へ出射された透過光は、光検出器
出力された２１１〜２１Ｎにより検出される。
【００１７】
光検出器２１１〜２１Ｎの出力信号（被検体Ｍからの透過光の検出信号）である透過光検出部２０の出力信号は、時間分解波形計測装置３に入力される。また、この時間分解波形計測装置３には、パルス光供給部１０から、パルス光源１１１〜１１４によるパルス光供給のタイミング信号を含む照射データが入力されている。時間分解波形計測装置３は、パルス光供給部１０からのパルス光の照射データを参照しつつ、透過光検出部２０からの透過光の検出データに基づいて、光検出器２１１〜２１Ｎのそれぞれで検出された透過光の時間分解波形を計測する。
【００１８】
光照射手段１、光検出手段２、及び時間分解波形計測装置３の動作制御及び必要なデータ処理等は、制御装置４によって行われる。本実施形態における制御装置４は、光ＣＴ計測制御装置４１と、画像再構成装置４２と、データ記憶装置４３とを有して構成されている。
【００１９】
光ＣＴ計測制御装置４１は、本光ＣＴ装置によって実行される光ＣＴ計測を制御する。具体的には、パルス光供給部１０のパルス光源１１によるパルス光の供給動作、透過光検出部２０の光検出器２１１〜２１Ｎによる透過光の検出動作、及び時間分解波形計測装置３による透過光の時間分解波形の計測動作など、装置各部の動作をそれぞれ制御する。
【００２０】
画像再構成装置４２は、時間分解波形計測装置３によって計測された透過光の時間分解波形から、被検体Ｍの内部情報を示す光ＣＴ画像を再構成する。また、データ記憶装置４３は、制御装置４においてデータ処理される各データ等を必要に応じて格納する。
【００２１】
また、この制御装置４には、表示装置５が接続されている。表示装置５には、例えば、画像再構成装置４２において再構成された被検体Ｍの光ＣＴ画像が表示される。あるいは、透過光の時間分解波形などの各データや、本装置における光ＣＴ計測の実行状況に関するデータなどを表示する構成としても良い。
【００２２】
図２は、図１に示した光ＣＴ装置について、そのパルス光供給部１０及び透過光検出部２０の構成の一例を示す平面図である。
【００２３】
図２に示す構成例においては、図１に関して上述したように、２次元ファントムである被検体Ｍの左側面上に光照射位置が設定され、パルス光源１１からのパルス光は、入力用光ガイド１２によりパルス光入力部１３に接続されている。
【００２４】
これに対して、被検体Ｍの右側面上には、８個の光検出位置が等間隔で配列されている。そして、これら８個の光検出位置（Ｎ＝８）のそれぞれに透過光出力部２３１〜２３８が設けられ、光ガイド２２１〜２２８により、透被検体Ｍからの透過光を光検出器２１１〜２１８へと導いている。
【００２５】
次に、測定手順を説明する。光照射手段１に設けられたパルス光源１１から所定のタイミングでパルス光を供給する。供給されたパルス光は、パルス光入力用の光ガイド１２及びパルス光入力部１３を介して、所定の光照射位置から被検体Ｍへと照射される。被検体Ｍ内に入力されたパルス光は、その各部において吸収や散乱を受けつつ被検体Ｍ内を伝搬する。そして、最終的に被検体Ｍ内を透過して光検出位置に到達したパルス光の光成分は、透過光出力部２３１〜２３８及び光ガイド２２１〜２２８を介して、それぞれの光検出器２１１〜２１８によって検出される。検出した光の光強度を表す信号は、光検出器２１１〜２１８から時間分解波形計測装置３へと送られ、各光強度の時間変化を示す時間分解波形が求められ、制御装置４へと送られる。
【００２６】
こうして、被検体Ｍの左側面側から光を照射して右側面側から出射する透過光の検出を行った後、パルス光入力部１３と、透過光出力部２３１〜２３８との対応関係を維持したまま、パルス光入力部１３と、透過光出力部２３１〜２３８の被検体Ｍに対する位置を順次、変更して同様に変更した位置における透過光をそれぞれ測定する。具体的には、まず、上述したように光照射位置を被検体Ｍの左側面側とし、光検出位置を被検体Ｍの右側面側として測定を行った後、光照射位置を被検体Ｍの上面側とし、光検出位置を被検体Ｍの下面側として測定を行い、次に、光照射位置を被検体Ｍの右側面側とし、光検出位置を被検体Ｍの左側面側として測定を行い、最後に、光照射位置を被検体Ｍの下面側とし、光検出位置を被検体Ｍの上面側として測定を行うことで、各方角からの測定を終了する。
【００２７】
次に、図１に示した光ＣＴ装置における被検体Ｍの内部情報の算出、及び光ＣＴ画像の再構成について説明する。図３は、図１に示した光ＣＴ装置における画像再構成装置の構成の一例を示すブロック図である。
【００２８】
図３に示す画像再構成装置４２は、波形データ取得部４２１と、補間情報算出部４２２と、内部情報算出部４２３とを有して構成されている。ここで、本構成例における内部情報算出部４２３では、光拡散方程式を離散化して有限要素法に基づいて被検体Ｍの内部情報である吸収係数や散乱係数等の分布を算出する。画像再構成装置４２は、ＣＰＵ、ＲＯＭ、ＲＡＭ、ハードディスク等の記憶装置、各種インターフェース等によって構成されており、各部が独立したハードウェアによって構成されていてもよいし、共通のハードウェア上で機能するソフトウェアとして構成されていてもよい。
【００２９】
図３に示した画像再構成装置４２の構成を参照しつつ、本発明による光ＣＴ画像の画像再構成方法について説明する。図４は、図３に示した画像再構成装置による画像再構成方法を概略的に示すフローチャートである。
【００３０】
時間分解波形計測装置３によって計測され、制御装置４の画像再構成装置４２へと送られた透過光の時間分解波形は、波形データ取得部４２１に入力される。波形データ取得部４２１は、透過光の時間分解波形データＹ２ｉ−１，j（ｔ）（ｉは検出点の位置であり、ｊはそのときの方角を示す。ここで、ｉ＝１〜８、ｊ＝１〜４）を取得する（ステップＳ１０１、波形データ取得ステップ）。
【００３１】
次に、補間情報算出部４２２は、得られた透過光の時間分解波形Ｙ２ｉ−１，j（ｔ）から時間的空間的補間により、各光検出位置の中間位置である仮想検出位置における透過光の時間分解波形Ｙ２ｉ，j（ｔ）（ｉ＝１〜７、ｊ＝１〜４）を求める（ステップＳ１０２、時間的空間的補間ステップ）。ここで、時間的空間的補間とは、同一の光出射位置（ｊ）に対する実光検出位置と仮想光検出位置における透過光の時間分解波形Ｙｋ，j（ｔ）（ｋ＝１〜１５）を、時間、光強度およびｋをそれぞれ座標軸として３次元空間に配置した場合（図５参照。）に、光強度によって形成される面が滑らかな曲面となるように補間を行うことを意味する。
【００３２】
図６は、ある特定の時刻ｔｘにおけるＹｋ，j（ｔｘ）の光強度の検出位置（仮想検出位置を含む）における分布を、図７は、特定のｋｘにおけるＹｋｘ，j（ｔ）の光強度の時間分布をそれぞれ示している。図６では、計測データを二重丸で示し、空間的時間的補間により求めた補間曲面上の点として求めた仮想検出位置におけるデータを白丸で示している。仮想光検出位置だけでなく、実光検出位置においても空間的時間的補間により実データを時間的に補間するデータを求めてもよい。
【００３３】
次に、内部情報算出部４２３において、求めた実検出位置と仮想検出位置双方の時間分解波形データに基づいて、被検体Ｍの内部情報の算出が行われる（ステップＳ１０３、内部情報算出ステップ）。この内部情報としては、物体内部における吸収係数μａ、散乱係数μｓがある。
【００３４】
物体内部における吸収係数μａ、散乱係数μｓが得られたら、これらをもとにして視覚情報として再構成し（ステップＳ１０４、画像再構成ステップ）。得られた吸収係数μａ、散乱係数μｓの分布を示す光ＣＴ画像を出力する（ステップＳ１０５）。
【００３５】
ここで、被検体Ｍの内部情報算出に用いられる計算モデルとしては、例えば、フィルタ法及び有限要素法を組み合わせた計算モデルを適用することが好ましい。このような計算モデルによれば、特に高精度かつ高速で、内部情報の算出及び光ＣＴ画像の再構成を行うことができる。また、これ以外の計算モデルを用いても良い。
【００３６】
図４のフローチャートに示した画像再構成方法について、その具体的な実施例を説明する。以下においては、被検体Ｍの内部情報を示し光ＣＴ計測による画像化の対象となる物理量μとして、主に吸収係数μａを想定している。
【００３７】
図２に示した２次元ファントムの被検体Ｍについて画像再構成を行う場合、図２の被検体Ｍ中に点線で示したように、被検体Ｍ内を２次元メッシュ状に区分して、被検体Ｍの各部に対応する計測要素を作成する。図２に示した区分例では、２次元の被検体Ｍを１６×１６の計測要素に区分している。
【００３８】
これらの区分された計測要素は、被検体Ｍの各部での物理量μの算出や、算出された物理量μの計測値の分布から再構成される光ＣＴ画像の表示などに用いられる。ただし、物理量μの算出を高精度で行うため、例えば、物理量μの算出においては、８０×８０にメッシュ状に区分した計測要素を用い、物理量μの計測値の表示においては、区分を粗くした１６×１６の計測要素を用いることとしても良い。
【００３９】
このような被検体Ｍの計測要素への区分については、物理量μに対して必要とされる算出精度や、物理量μを算出する演算処理の処理速度などを考慮して、その区分幅や区分数を好適に設定することが好ましい。なお、被検体Ｍの計測要素への分割幅は、例えば１ｍｍ程度である。
【００４０】
まず、計算モデルの作成について説明する。光拡散方程式及び有限要素法を用いて、光ＣＴの計測対象である被検体Ｍに対する順問題（forward problem）での計算モデルを作成する。
【００４１】
光照射手段にパルス光源として設けられているパルスレーザから、被検体Ｍに対してパルスレーザ光が照射されるとする。被検体Ｍの媒質中に入力されたパルス光は、その各部において吸収や散乱を受けつつ被検体Ｍ内を伝搬する。このような媒質中での光の伝搬は、以下の光拡散方程式
【００４２】
【数１】
2005328916.tif 000003
によって表される。ここで、Φ（ｒ，ｔ）は被検体Ｍ内の３次元位置ｒ（ｘ，ｙ，ｚ）及び時間ｔでの光密度、ｃは被検体Ｍの媒質中での光速、μａ（ｒ）は吸収係数、ｑ（ｒ，ｔ）は光源の光密度である。なお、被検体Ｍ内の位置ｒは、被検体Ｍが２次元ファントムである場合には、２次元位置ｒ（ｘ，ｙ）である。また、Ｄ（ｒ）は光拡散係数であり、次式
【００４３】
【数２】
2005328916.tif 000004
によって与えられる。ここで、μｓ’（ｒ）は換算散乱係数である。
【００４４】
次に、式（１）の偏微分方程式を状態空間システム（state-space system）に変換するため、ガラーキン近似法（Galerkin approximation）を適用し、位置ｒと時間ｔとに分離（離散化）する近似を行う。ここでの有限要素法（ＦＥＭ：Finite Element Method）を用いた計算モデルは、以下のガラーキン近似法に基づいている。
【００４５】
まず、ガラーキン近似のための基底として、位置ｒについての基底関数｛φ１，φ２，…，φｎ｝をとる。このとき、光密度の解Φ（ｒ，ｔ）は、次式
【００４６】
【数３】
2005328916.tif 000005
によって近似することができる。ここで、ｎは上記したガラーキン近似の基底の次元である。
【００４７】
式（３）に示す近似により、被検体Ｍ内の空間を表す元の空間における位置ｒは、位置ｒがｉ番目の部分空間（Subspace)にあることを示す番号ｉによって表されることとなる。すなわち、式（３）により、位置ｒによる３次元空間が、部分空間ｉ（ｉ＝１〜ｎ）による１次元アレイの状態空間システムに変換される。
【００４８】
また、光源の光密度ｑ（ｒ，ｔ）について、被検体Ｍに対してパルス光を照射するパルス光源がｋ個あるものとして、光密度Φ（ｒ，ｔ）と同様に近似する。このとき、光源の光密度ｑ（ｒ，ｔ）は、次の式（４）
【００４９】
【数４】
2005328916.tif 000006
によって近似される。ここで、番号ｎｉはｉ番目の光源に相当する部分空間を示している。
【００５０】
ガラーキン近似法によるこれらの近似式（３）及び（４）を、式（１）に示した光拡散方程式に適用することにより、計算モデルに用いる状態方程式を求めることができる。以下、被検体Ｍ内の媒質が一様（Homogeneous）な場合、及び非一様（Heterogeneous）な場合について、それぞれ説明する。
【００５１】
まず、被検体が一様な場合（Homogeneous case）について説明する。被検体内の媒質が一様であれば、吸収係数μａ（ｒ）及び散乱係数μｓ（ｒ）は、いずれも空間位置ｒによらず被検体内で一定である。このとき、光拡散方程式（１）に近似式（３）及び（４）を適用すると、次の式（５）
【００５２】
【数５】
2005328916.tif 000007
が得られる。ここで、式（５）中のＥ、Ｇ、Ｓ、及びＦは、それぞれ
【００５３】
【数６】
2005328916.tif 000008
である。
【００５４】
行列Ｅは正方行列であり、特異行列ではないので、式（５）は次式
【００５５】
【数７】
2005328916.tif 000009
のように書き直すことができる。また、この式（７）において、−（ＤＥ−１Ｇ＋ｃμａＩ）をＡとおき、Ｅ−１ＦをＢとおけば、式（８）
【００５６】
【数８】
2005328916.tif 000010
が得られる。さらに、光検出器の位置を表す応答行列Ｃを用いると、光検出器の応答ｙ（ｔ）は、次式
【００５７】
【数９】
2005328916.tif 000011
となり、線型の状態空間方程式が得られる。
【００５８】
次に、被検体が非一様な場合（Heterogeneous case）について説明する。被検体内の媒質が非一様であれば、吸収係数μａ（ｒ）は空間位置ｒに依存する。この吸収係数μａの位置による変化は、ｎ×ｎの対角行列Ｍａを用いて表すことができる。ここで、対角行列Ｍａにおいて、その対角成分Ｍａ（ｉ，ｉ）は、位置（部分空間）ｉでの吸収係数μａの値を示している。また、散乱係数μｓについては、被検体内で一定であると仮定しておく。このとき、光拡散方程式（１）に近似式（３）及び（４）を適用すると、次の式（１０）
【００５９】
【数１０】
2005328916.tif 000012
が得られる。ここで、式（１０）中のＳは、
【００６０】
【数１１】
2005328916.tif 000013
である。
【００６１】
行列Ｅは正方行列であり、特異行列ではないので、式（１０）は次式
【００６２】
【数１２】
2005328916.tif 000014
のように書き直すことができる。
【００６３】
吸収係数μａの位置による変化を表す行列Ｍａが、被検体が一様な場合における吸収係数μａからの摂動として、次式
【００６４】
【数１３】
2005328916.tif 000015
によって表されるとする。このとき、行列Ｍａは、次の式（１４）
【００６５】
【数１４】
2005328916.tif 000016
となる。ここで、μａＩは、被検体が一様媒質な場合での一定の吸収係数μａを表す行列であり、また、Ｐは、その対角成分に被検体が非一様媒質であることによる摂動情報を含む対角行列である。
【００６６】
非一様媒質での吸収係数μａについて、このような表現を用いることにより、被検体が非一様な場合での上述した状態空間システムは、以下のようになる。
【００６７】
【数１５】
2005328916.tif 000017
または、この式（１５）と同等な式として、次の式（１４）
【００６８】
【数１６】
2005328916.tif 000018
が得られる。
【００６９】
この式（１６）は、被検体が一様な場合での式（８）に対して、式中のＡをＡ＋ｃＰと入れ換えたものとなっている。また、非一様媒質での式（１６）における行列Ａ及びＢは、一様媒質での式（８）におけるものと等しくなる。このような非一様媒質の表現は、後述する逆問題の解法にとても便利なものとなる。
【００７０】
以上説明した被検体が一様媒質の場合、及び非一様媒質の場合をまとめると、被検体Ｍに対する計算モデルは、以下の式（１７．１）及び（１７．２）
【００７１】
【数１７】
2005328916.tif 000019
となる。ここで、Ａはシステム行列、Ｄは一様媒質と非一様媒質との吸収係数の差を表す対角行列、Ｂｕ（ｔ）は光源の位置情報を含むシステム入力値、Ｃは光検出器の位置情報を含む値、ｘ（ｔ）はシステムの状態（ＦＥＭにおける各ボクセルの時間ｔでの光密度）、ｙ（ｔ）は計測値である。
【００７２】
次に、状態空間システム（state-space system）を構築する。
【００７３】
被検体として２次元の媒質を考えると、空間位置ｒは２次元位置ｒ（ｘ，ｙ）となる。このような２次元媒質をｘ軸方向にＭｘ、ｙ軸方向にＭｙの区分数で区分すると、被検体がＭｘ×Ｍｙ個の計測要素である最小の正方形要素に区分されることになる。
【００７４】
このとき、上述したガラーキン近似の基底関数｛φ１，φ２，…，φｎ｝において、基底関数の個数がｎ＝ＭｘＭｙとなり、したがって、状態空間システムの次元がｎ＝ＭｘＭｙとなる。また、それぞれの空間位置ｒ（ｘ，ｙ）に対応する部分空間を示す番号ｉは、次式
【００７５】
【数１８】
2005328916.tif 000020
によって与えられる。
【００７６】
基底関数φｉの値は、次のように定義される。
【００７７】
【数１９】
2005328916.tif 000021
また、行列Ｅは、次元ｎの対角行列
【００７８】
【数２０】
2005328916.tif 000022
となる。
【００７９】
また、行列Ｇは、
【００８０】
【数２１】
2005328916.tif 000023
と表される。ただし、この式（２１）中のＨは、
【００８１】
【数２２】
2005328916.tif 000024
である。
【００８２】
ここで、計算モデルを作成する被検体が３次元媒質である場合や、物理量μの計測値の算出精度を高くするために被検体の区分を細かくして計測要素を小さくとった場合などには、ＦＥＭを用いた計算モデルでの行列は、非常に大きい行列となる。ただし、その行列は、スパース行列となる。
【００８３】
このように、計算モデルでの行列が大きくなった場合には、その演算に時間がかかることとなる。これに対して、直交射影によるモデル縮小で近似を行うことにより、その演算速度を飛躍的に向上することができる。具体的には、ランチョス法によってクリロフ部分空間（Krylov subspace）上に射影する方法を用いることができる。
【００８４】
まず、上述した行列Ａ、行列Ｂ、及び光検出器の位置を与える行列Ｃを、ランチョス法によってクリロフ部分空間Ｋｍ
【００８５】
【数２３】
2005328916.tif 000025
上に射影する。ここで、υ１は光源の単位ベクトルを表している。また、ランチョス法のアルゴリズムは、以下に示す通りである。
【００８６】
【数２４】
2005328916.tif 000026
ここで、ｑ１は光源を示しており、上記したυ１と同じである。
【００８７】
以上をまとめると、状態空間システムは、以下の式（２５．１）及び（２５．２）
【００８８】
【数２５】
2005328916.tif 000027
によって表現される。また、この解は、
【００８９】
【数２６】
2005328916.tif 000028
となる。
【００９０】
また、ランチョス法を用いて直交射影すると、射影されたシステムは、以下の式（２７．１）及び（２７．２）
【００９１】
【数２７】
2005328916.tif 000029
のようになる。これにより、状態空間システムの次元をより小さくすることができる（ｎ≪Ｎ）。また、この解は、
【００９２】
【数２８】
2005328916.tif 000030
となる。
【００９３】
ここで、連続時間ｔを離散時間とするため、次の式（２９）によってＳ及びＴを計算する。
【００９４】
【数２９】
2005328916.tif 000031
ここで、この式（２９）の左辺のＵ、Ｖ、及び右辺で０となっているところは、意味をなさない無関係な数値である。この演算によって、連続時間から離散時間へと変換することができる。そして、次の式（３０）
【００９５】
【数３０】
2005328916.tif 000032
によって、離散時間での応答を計算することができる。
【００９６】
以上により、基本的な時間応答ｘが得られることとなり、被検体Ｍに対する順問題の計算、すなわち被検体モデルに入力された光に対する応答波形の計算が可能となる。
【００９７】
次に、非線型フィルタを用いた被検体の光ＣＴ画像の再構成、すなわち、被検体Ｍに対する逆問題（inverse problem）について説明する。例えば、被検体の内部情報として算出対象となっている物理量μが吸収係数μａであるとすれば、光ＣＴ画像の再構成では、被検体Ｍの媒質中での吸収係数μａの分布を推定して画像化する。
【００９８】
逆問題を解いて画像再構成を行うため、上述した式（８）及び（９）を以下の式（３１）及び（３２）に変形する。
【００９９】
【数３１】
2005328916.tif 000033
【０１００】
【数３２】
2005328916.tif 000034
ここで、Ｄは一様媒質と非一様媒質との吸収係数の差を表す対角行列であり、未知のパラメータである。また、Ｄ＝０のときが一様（均質）媒質の場合を表し、Ｄ≠０のときが非一様（不均質）媒質の場合を表している。
【０１０１】
これらの式（３１）及び（３２）は、実際の光ＣＴ計測において発生するノイズを含んだ確率的な計算モデルを表している。また、Ｄは具体的には、吸収係数μａの基準値μＭａからの計測値差δμａの対角行列となっている。ここで、
【０１０２】
【数３３】
2005328916.tif 000035
である。また、ｗはシステムノイズ、Ｇはシステムにノイズがどのように入るかを記述するノイズ入力行列、ｖは光検出器による計測の不確定さを表すセンサーノイズをそれぞれ表している。
【０１０３】
ノイズが白色ガウス分布で互いに独立であり、また、モデルのパラメータが正確であれば、最適状態の推定にカルマンフィルタ（Kalman Filter）を用いることができる。また、セカンドオーダーフィルタ（Second Order Filter）の考え方は、非線型フィルタ理論からきている。
【０１０４】
まず、吸収の変化を表すｄ（ｔ）を用いて、状態ベクトルｚを次のように変形する。
【０１０５】
【数３４】
2005328916.tif 000036
この式（３４）を用いると、式（３１）、（３２）は以下の式（３５）、（３６）
【０１０６】
【数３５】
2005328916.tif 000037
【０１０７】
【数３６】
2005328916.tif 000038
のようになる。ここで、
【０１０８】
【数３７】
2005328916.tif 000039
である。
【０１０９】
また、セカンドオーダーフィルタを用いて計測値差δμａを予測する。連続的なモデルを記述している式（３５）、（３６）を離散化すると、
【０１１０】
【数３８】
2005328916.tif 000040
【０１１１】
【数３９】
2005328916.tif 000041
となる。セカンドオーダーフィルタでは、ｚについての２次導関数まで考慮したｆｄは、
【０１１２】
【数４０】
2005328916.tif 000042
となる。ここで、
【０１１３】
【数４１】
2005328916.tif 000043
とする。すなわち、ｆｄ（ｚｋ，ｕｋ）をｚｋの推測値のまわりでテイラー展開し、その２次の項まで考慮する。
【０１１４】
このとき、シグナルモデルは、以下のようになる。
【０１１５】
【数４２】
2005328916.tif 000044
また、フィルタ方程式は、以下のように与えられる。
【０１１６】
まず、初期化情報（Initialization, a priori information）は、
【０１１７】
【数４３】
2005328916.tif 000045
で与えられる。
【０１１８】
また、時間更新（Time Update）は、
【０１１９】
【数４４】
2005328916.tif 000046
で与えられる。ここで、
【０１２０】
【数４５】
2005328916.tif 000047
である。
【０１２１】
また、計測値更新（Measurement Update）は、
【０１２２】
【数４６】
2005328916.tif 000048
で与えられる。ここで、
【０１２３】
【数４７】
2005328916.tif 000049
である。
【０１２４】
これらの時間更新及び計測値更新を繰り返すことによって、計測値差δμａ及び計測値μａが求められる。以上の計算を、計測値との誤差が充分小さく最適な値となるまで繰り返し、終了した時点での吸収係数μａの値を出力して計算を完了する。
【０１２５】
ここで、本発明にかかる補正の有効性を確認するシミュレーション結果について説明する。図８は、画像再構成を行う被検体モデルにおける、２次元ファントム内での吸収体の配置を示す図である。ここでは、被検体を１６×１６の計測要素に区分している。以下、図の左下端の領域を（１，１）とし、右上端を（１６，１６）とし、左下端から右にｎ番目、上にｍ番目の領域を（ｎ，ｍ）で表すものとする。
【０１２６】
この被検体モデルは、全体の散乱係数が一律に１．０ｍｍ−１であり、周辺媒質の吸収係数は０．０１ｍｍ−１である。この周辺媒質中の（７，９）〜（９，１１）領域（３×３領域）に吸収係数０．０２ｍｍ−１の吸収体Ａが、（６，４）領域（１×１領域）に吸収係数０．０２ｍｍ−１の吸収体Ｂが、（１１，７）〜（１２，８）領域（２×２領域）に吸収係数０．０１５ｍｍ−１の吸収体Ｃがそれぞれ配置されている。このような被検体モデルに対し、８個の検出器で測定を行った生の時間分解波形データにより画像再構成を行った場合と、各検出器の中間位置で補間した７個の仮想位置データを含む１５個の時間分解波形データにより画像再構成を行った場合で再構成結果を比較した。
【０１２７】
図９は、補間を行うことなく測定データのみから行った画像再構成結果を示す図である。また、図１０は補間データを用いた場合の画像再構成結果を示す図である。図９に示される再構成結果では、吸収係数が強く、比較的大きな領域である吸収体Ａのみしか十分には把握ができておらず、これより領域の小さな吸収体Ｂ、吸収体Ｃについては把握ができていない。これに対して、図１０に示される本発明にかかる画像再構成結果によれば、吸収体Ｂ、吸収体Ｃについてもその存在を検出することができる。
【０１２８】
図１１（ａ）〜図１１（ｃ）は、図８〜図１０の画像結果から逆に各検出位置における時間分解波形データを構成し、これを連結して示したものである。補間データを用いた図１１（ｃ）の場合は、使わない図１１（ｂ）の場合に比較して、本来の状態である図１１（ａ）のピーク状態の再現性が高く、精度良く再現できていることが確認された。
【０１２９】
本発明にかかる画像再構成方法は、特に、生体等を計測する場合に有効である。生体においては、周囲を光拡散性の良好な媒体（頭なら頭蓋骨、体なら体脂肪等）により覆われていることもあり、照射位置から照射された光のうち検出位置まで直進してきて到達する成分はほとんどなく、大部分が両者を直線で結んだ進路上から離れた位置を経由して散乱により到達する光であるからである。
【０１３０】
被検体Ｍの検出位置近傍の拡散性が小さい場合には、被検体Ｍの外周を拡散性の良好な媒体（緩衝材）で覆い、この緩衝材の外側で測定を行うとよい。これにより、被検体Ｍの外周において時間分解波形データの位置依存性が大きく、補間に適さない場合（図５の曲面の凹凸が微細な場合に該当する。）であっても、この緩衝材の外側では、時間分解波形データの位置変化を滑らかにして補間に適した状況とすることができる。
【０１３１】
本発明による光ＣＴ装置及び画像再構成方法は、上記した実施形態及び実施例に限られるものではなく、様々な変形が可能である。例えば、光ＣＴ装置における光照射手段及び光検出手段の構成については、図１に示した構成に限らず、光ＣＴ計測の対象とする被検体Ｍの形状や大きさ等に応じて、様々な構成とすることが可能である。
【０１３２】
補間位置は検出箇所の中間位置に１つずつ設ける形態に限られるものではなく、検出位置を等間隔に配置しなくともよい。検出位置自体を不均等に配置したり、検出位置間に２個以上の補間位置を設けたり、検出位置間の一方の検出位置に偏った位置に補正点を配置してもよい。また、例えば、ｍ個の検出点から、補間によりこれより多いｎ個の仮想検出点におけるデータを求め、この仮想検出点のデータのみを用いて画像再構成を行うこともできる。これは、物理的制約で測定位置を再構成演算の際のメッシュに合致させて配置することが困難なときに有効であり、補間によりメッシュに合わせた測定位置を求めて境界条件を的確に設定することを可能とする。
【０１３３】
本発明によれば、実際の検出箇所を増やすことなく、仮想的に検出箇所を増やすことにより、検出箇所が少ない場合でも精度良く測定を行うことができる。また、検出箇所を増やす必要がないため、測定装置の軽量化、コンパクト化が可能となり、測定時の被検体Ｍへの負担を低減することができる。特に、生体測定の場合には被験者の負担を減らすことができ、好ましい。また、制御装置へ送るまでのデータ量を削減することができるので、データの取り扱い性が向上する。例えば、測定と画像再構成とを別システムで行うような場合に、システム間で転送するデータ量が少なくて済むという利点もある。
【図面の簡単な説明】
【０１３４】
【図１】本発明による光ＣＴ装置１００の一実施形態の構成を概略的に示すブロック図である。
【図２】図１に示した光ＣＴ装置について、そのパルス光供給部１０及び透過光検出部２０の構成の一例を示す平面図である。
【図３】図１に示した光ＣＴ装置における画像再構成装置の構成の一例を示すブロック図である。
【図４】本発明による画像再構成方法を概略的に示すフローチャートである。
【図５】時間的空間的補間について説明する図である。
【図６】ある特定の時刻ｔｘにおけるＹｋ，j（ｔｘ）の光強度の検出位置（仮想検出位置を含む）における分布を示す図である。
【図７】特定のｋｘにおけるＹｋｘ，j（ｔ）の光強度の時間分布を示す図である。
【図８】画像再構成を行う被検体モデルにおける、２次元ファントム内での吸収体の配置を示す図である。
【図９】補間を行うことなく測定データのみから行った画像再構成結果を示す図である。
【図１０】補間データを用いた場合の画像再構成結果を示す図である。
【図１１】画像再構成結果から求めた時間分解波形を比較したグラフである。
【図１２】拡散媒体で囲まれた被検体内の透過光を示す図である。
【符号の説明】
【０１３５】
１…光照射手段、２…光検出手段、３…時間分解波形計測装置、４…制御装置、５…表示装置、１０…パルス光供給部、１１…パルス光源、１２…入力用光ガイド、１３…パルス光入力部、２０…透過光検出部、２２…出力用光ガイド、２３…透過光出力部、４１…光ＣＴ計測制御装置、４２…画像再構成装置、４３…データ記憶装置、１００…光ＣＴ装置、４２１…波形データ取得部、４２２…補間情報算出部、４２３…内部情報算出部、Ｍ…被検体。
【図１】
2005328916.tif 000050
【図２】
2005328916.tif 000051
【図３】
2005328916.tif 000052
【図４】
2005328916.tif 000053
【図５】
2005328916.tif 000054
【図６】
2005328916.tif 000055
【図７】
2005328916.tif 000056
【図８】
2005328916.tif 000057
【図９】
2005328916.tif 000058
【図１０】
2005328916.tif 000059
【図１１】
2005328916.tif 000060
【図１２】
2005328916.tif 000061