光断層情報生成装置、光強度分布算出方法および光強度分布算出プログラム
【課題】蛍光体の濃度分布を含む光断層画像を得るときの処理負荷の低減、装置の簡略化を可能とする。
【解決手段】計測部では、検体に励起光を照射して、これにより得られる蛍光の強度の計測データを取得する。また、画像処理部では、蛍光体の濃度分布に基づく蛍光体の吸収係数の初期値を設定すると、予め設定されている検体の吸収係数及び拡散係数に基づいて蛍光強度分布を演算し、計測データと演算結果を比較する（ＳＴ１０００〜１０８０）。このとき、所定値以内でなければ、最適化手法による光拡散方程式を用いた逆問題計算を行うことにより誤差を所定値以内とする蛍光体の吸収係数を推定し、この吸収係数に基づいた濃度分布から蛍光の強度分布の演算、誤差の評価を反復し、誤差が所定値以内となる濃度分布を取得する（ステップ１１００〜１１８０）。
【解決手段】計測部では、検体に励起光を照射して、これにより得られる蛍光の強度の計測データを取得する。また、画像処理部では、蛍光体の濃度分布に基づく蛍光体の吸収係数の初期値を設定すると、予め設定されている検体の吸収係数及び拡散係数に基づいて蛍光強度分布を演算し、計測データと演算結果を比較する（ＳＴ１０００〜１０８０）。このとき、所定値以内でなければ、最適化手法による光拡散方程式を用いた逆問題計算を行うことにより誤差を所定値以内とする蛍光体の吸収係数を推定し、この吸収係数に基づいた濃度分布から蛍光の強度分布の演算、誤差の評価を反復し、誤差が所定値以内となる濃度分布を取得する（ステップ１１００〜１１８０）。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、光を用いたトモグラフィ（Tomography）技術に係り、特に、励起光を照射して、この励起光によって発光する発光体の分布情報を用いて光断層画像を取得可能とする光断層情報生成装置、光強度分布算出方法および光強度分布算出プログラムに関する。
【背景技術】
【０００２】
生体などの断層画像を取得する方法としては、Ｘ線を用いたＸ線ＣＴ、超音波を用いた超音波ＣＴ、核磁気共鳴を適用したＮＭＲ−ＣＴ、プロトンなどの粒子線を用いた陽子線ＣＴなどがある。また、生体においては、光透過性を有することが知られており、小動物の断層画像に光を用いる光ＣＴ（Computed Tomography）が提案されている（例えば、特許文献１、特許文献２参照）。
【０００３】
生体に照射された光は、生体内で散乱し、散乱した光が生体の周囲から射出される。光ＣＴでは、生体内で乱反射して生体の周囲から射出される光を検出して、電気信号を取得し、それぞれの電気信号に対して所定の信号（画像信号）処理を施して得られる情報から画像再構成を行うことにより、生体の断層画像が得られるようにしたものである。
【０００４】
一方、病理学的な実験分野では、所定の波長の光で発光する蛍光体を含む薬剤などを生体に投与して、生体中での、その薬剤の移動、分布、特定部位へ集積するときの集積過程などを観察するときに光ＣＴ（以下、蛍光ＣＴとする）を用いることができる。すなわち、蛍光体に対する励起光を生体に照射し、この励起光に応じて生体から射出される発光（蛍光）を検出して、生体の二次元的な断層画像または三次元的な断層画像を再構成する。これにより、この断層画像から蛍光体、蛍光体を含む試薬や細胞などの位置、量などの情報を得ることができる。
【０００５】
このような蛍光ＣＴにおいても、生体の表面の一点へ励起光を照射し、これにより生体から射出される散乱蛍光を多点で検出する。これを、励起光の照射位置を変えながら繰り返すことで、照射点数×検出点数だけのデータを取得する。これらのデータの間には、蛍光物質の分布、体内での光の散乱、吸収特性に応じた関係が成立することから、この関係に基づいて断層画像を再構成することができる。
【０００６】
ところで、この蛍光ＣＴにおいて、断層画像の再構成を行うために蛍光体の濃度分布を演算するときには、励起光強度分布および蛍光強度分布の２つの光強度分布について逆問題計算を行うことが提案されている（例えば、特許文献３、非特許文献１参照）。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開平１１−１７３９７６号公報
【特許文献２】特開平１１−３３７４７６号公報
【特許文献３】米国特許出願公開第２００７／０２８６４６８号明細書
【非特許文献】
【０００８】
【非特許文献１】S. R. Arridge “Optical tomography in medical” Inverse Problems １５（１９９９）Ｒ４１−９３
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、光強度分布を得るための逆問題計算は、順問題計算と比較して計算負荷が大きく装置構成が大規模になりがちであり、計算に長時間を要するという問題がある。近似計算によりこれらの問題を軽減するアプローチも存在するが、ノイズ環境下では近似自体が精度良くできないことも多い。
【００１０】
本発明は上記事実に鑑みてなされたものであり、例えば蛍光ＣＴにおいて装置構成を簡略化するとともに、光強度分布を算出するにあたり計算処理負荷を軽減し、処理時間を短縮化することができる光断層情報生成装置、光強度分布算出方法および光強度分布算出プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１１】
上記目的を達成するために、本発明に係る光断層情報生成装置は、光拡散方程式に基づき、検体から射出される光の強度分布の計算値を求め、当該計算値と実測値との差の大小に応じて、ＡＲＴ（代数的再構成技術）を用いた再計算を行って前記計算値の更新値を得、当該更新値を用いて前記検体の光断層情報を生成する光断層情報生成装置であって、前記検体の光学特性値から、当該光学特性値の変化に対する前記光強度分布の変化の割合を示すヤコビ行列を算出するヤコビ行列算出手段と、前記ヤコビ行列を特異値分解して対角行列を得る特異値分解手段と、前記対角行列の閾値以下の対角成分を０で置換した単位対角行列を取得する単位対角行列取得手段と、前記ＡＲＴを用いた再計算において、ＡＲＴの逐次近似式に現れる前記計算値と実測値との差に対し前記単位対角行列を用いた置換を施して近似を行い、前記更新値を得る逐次近似手段と、を具備する構成を採る。
【００１２】
また、本発明に係る光強度分布算出方法は、光拡散方程式に基づき、検体から射出される光の強度分布の計算値を求め、当該計算値と実測値との差の大小に応じて、ＡＲＴを用いた再計算を行って前記計算値を更新する光強度分布算出方法であって、前記検体の光学特性値から、当該光学特性値の変化に対する前記光強度分布の変化の割合を示すヤコビ行列を算出するステップと、前記ヤコビ行列を特異値分解して対角行列を得るステップと、前記対角行列の閾値以下の対角成分を０で置換した単位対角行列を取得するステップと、前記ＡＲＴを用いた再計算において、ＡＲＴの逐次近似式に現れる前記計算値と実測値との差に対し前記単位対角行列を用いた置換を施して近似を行って、前記計算値を更新するステップと、を具備するようにした。
【発明の効果】
【００１３】
本発明によれば、光強度分布を算出する際の計算処理負荷が軽減されると共に、装置構成の簡略化、計測時間の短縮化を実現することができる。
【図面の簡単な説明】
【００１４】
【図１】実施の形態１に係る光断層情報生成装置の主要な外観構成を示した図である。
【図２】実施の形態１に係る計測部および画像処理部の詳細な構成を示した図である。
【図３】実施の形態１に係る画像処理部の入出力データおよびデータ処理を説明するための機能ブロック図である。
【図４】（Ａ）は検体内での励起光の伝播を示す概略図、（Ｂ）は検体内での蛍光の伝播を示す概略図である。
【図５】検体の光学的特性の一例を示す線図である。
【図６】実施の形態１に係る光強度分布算出方法の手順を示したフローチャートである。
【図７】実施の形態１に係る更新処理部のより詳細な構成を示した機能ブロック図である。
【図８】抗体投与マウスを２次元カメラで腹部方向から撮影して得られた画像である。
【図９】（Ａ）は従来のＡＲＴ法による再構成画像、（Ｂ）は実施の形態１に係る光断層情報生成装置によって得られた再構成画像の例である。
【図１０】ヤコビ行列Ｊのランク数と対角行列Ｄにおける特異値の値との関係を特性曲線として示した図である。
【図１１】（Ａ）〜（Ｆ）は閾値と再構成画像の精度との関係を実際の抗体投与マウスのデータによって検証した結果を示した図である。
【図１２】実施の形態２に係る光断層情報生成装置内の計測部および画像処理部の主要な構成を示すブロック図である。
【図１３】実施の形態２に係る各受光ヘッドの関係等を説明するための図である。
【図１４】実施の形態２に係る光強度分布算出方法の手順を示したフローチャートである。
【発明を実施するための形態】
【００１５】
以下、添付図面を参照しながら本発明の実施の形態を説明する。ここでは、光断層情報生成装置として蛍光トモグラフィ装置を用いる場合を例にとって説明する。
【００１６】
（実施の形態１）
図１は、本発明の実施の形態１に係る光断層情報生成装置１００の主要な外観構成を示した図である。この光断層情報生成装置１００は、計測部１２と画像処理部１４とを含んでいる。
【００１７】
画像処理部１４は、計測部１２から出力される電気信号に基づいた画像処理（信号処理）を行う。画像処理部１４には、表示手段としてＣＲＴ、ＬＣＤなどのモニタ１６が設けられており、このモニタ１６に、計測部１２の計測結果に基づいた画像が表示される。
【００１８】
光断層情報生成装置１００では、小動物（例えばヌードマウス）などの生体を観察対象とする検体１８として、この検体１８から得られる光断層情報に基づいた画像（以下、光断層画像とする）をモニタ１６に表示したり、表示画像の画像データ（光断層情報）を各種の記憶媒体に記憶可能となっている。
【００１９】
計測部１２は、計測ユニット２０を備えている。計測ユニット２０は、リング状の機枠２２を備え、この機枠２２の軸心部が計測位置となっている。計測部１２では、前記検体１８が機枠２２内の計測位置に配置される。
【００２０】
計測ユニット２０には、計測位置へ向けて所定波長の光を照射する光源ヘッド２４と、前記生体から射出される光を検出光Ｌ１として、この検出光Ｌ１を検出する複数の受光ヘッド２６とが、所定の角度間隔で（所定の角度θずつずらされて）機枠２２に取付けられている。本実施の形態に適用した光断層情報生成装置１００では、一例として光源ヘッド２４から３０°ずつ（θ＝３０°）ずらして１１台の受光ヘッド２６を配列している。
【００２１】
この構成により、計測部１２では、光源ヘッド２４から照射した光に対して、検体１８から射出された検出光Ｌ１を、１１台の受光ヘッド２６のそれぞれによって並行して受光可能となっている。
【００２２】
また、計測部１２において、機枠２２が軸心に対して所定の角度ずつ機械的に回転するよう構成されている。これにより、計測中は、計測部１２において、検体１８の周囲の複数点へ向けて光源ヘッド２４から光が照射され、受光ヘッド２６はそれぞれの位置での検出光Ｌ１の受光が可能となっている。ここでは、一例として、機枠２２を角度θ（θ＝３０°）ずつ回転するようにしている。また、計測部１２では、検体１８の周囲の１２点に対して、光を照射し、それぞれの照射点における１１箇所で検出光Ｌ１の検出が可能となっている。なお、光源ヘッド２４の数、受光ヘッド２６の数、これらの配列及び光源ヘッド２４と受光ヘッド２６の移動量（回転量）などは、これに限るものではなく、任意の数、配列及び移動量（回転量）を適用することができる。
【００２３】
計測部１２において、基台２８上に支柱３０が立設されており、この支柱３０に対して機枠２２が回転可能に支持されている。また、支柱３０は、基台２８上に、機枠２２の軸方向（図１の紙面表裏方向）に沿って移動可能に支持されている。即ち、機枠２２は、回転可能となっていると共に、その回転軸の軸方向に沿っても移動可能となっており、機枠２２の回転軸方向に沿った検体１８の任意の部位に対し計測可能となっている。この構成により、３次元の再構成画像を取得することが可能になっている。
【００２４】
なお、このような機枠２２の回転機構及び移動機構は、任意の構成を適用することができる。また、計測部１２では、機枠２２を回転するようにしているが、これに限らず、機枠２２内に配置される検体１８を回転する構成としても良く、また、検体１８と機枠２２のそれぞれが回転されるものであっても良い。
【００２５】
また、計測部１２には、制御ユニット３２が設けられている。
【００２６】
図２は、計測部１２および画像処理部１４の詳細な構成を示した図である。なお、図１と同一の構成要素には同一の符号を付している。
【００２７】
この図に示されるように、制御ユニット３２は、マイクロコンピュータを含んで形成されたコントローラ３４を備えている。また、制御ユニット３２には、例えば、光源ヘッド２４を駆動する発光駆動回路３６、受光ヘッド２６のそれぞれから出力される電気信号を増幅する増幅器（Ａｍｐ）３８、増幅器３８から出力される電気信号をデジタル信号に変換するＡ／Ｄ変換器４０を備え、光源ヘッド２４の発光、各受光ヘッド２６での検出光Ｌ１の受光、受光した検出光Ｌ１の強度を示す測定データの生成が、コントローラ３４によって制御される。
【００２８】
計測部１２は、計測ユニット２０の機枠２２を回転駆動する回転モータ４２、機枠２２を軸方向に移動する移動モータ４４及び、それぞれの駆動回路４６、４８を含み、これらがコントローラ３４に接続された構成とすることができる。
【００２９】
一方、画像処理部１４は、ＣＰＵ５０、ＲＯＭ５２、ＲＡＭ５４、記憶手段とされるＨＤＤ５６、キーボードやマウスなどの入力デバイス５８、モニタ１６等がバス６０に接続された一般的構成のコンピュータが形成されている。これにより、画像処理部１４では、ＲＯＭ５２やＨＤＤ５６に記憶されたプログラム、図示しないリムーバブルメモリなどに記憶されたプログラムに基づいた各種の制御、信号処理、画像処理などが可能となっている。なお、本願では、ＲＯＭ５２、ＲＡＭ５４、ＨＤＤ５６および図示しないリムーバブルメモリを総称して単にメモリ５１と呼ぶこととする。
【００３０】
また、画像処理部１４と計測部１２の制御ユニット３２との間では、制御信号の送受信及びデータの送受信が可能となっている。なお、このような構成は、任意の通信インターフェイスを用いて構成することができる。
【００３１】
図３は、画像処理部１４の入出力データおよび内部のデータ処理を説明するための機能ブロック図である。
【００３２】
この図に示されるように、画像処理部１４は、読取部７０を備え、計測部１２での検体１８の計測を制御しながら、計測部１２（制御ユニット３２）から出力される計測データをメモリ５１に読み込む。
【００３３】
また、画像処理部１４は、演算処理部７２、評価部７４、更新処理部７６、断層情報生成部７８及び断層画像生成部８０を備えている。
【００３４】
演算処理部７２は、蛍光体６２の吸収係数分布を含む予め設定されている検体１８の光学特性値に基づいて光拡散方程式を用いた順問題計算によって蛍光の強度分布を演算する。評価部７４は、演算された蛍光の強度分布と計測データから得られる蛍光の強度分布の差分を評価する。更新処理部７６は、光拡散方程式の逆問題計算を行うことにより、評価部７４の評価結果から得られる差分を減少させるように蛍光の強度分布から蛍光体の濃度分布に基づく吸収係数分布等を設定する。演算処理部７２は、更新処理部７６において蛍光体の濃度分布に基づく吸収係数分布等の更新が行われると、更新された蛍光体の濃度分布に基づく吸収係数分布等に基づいた蛍光強度分布の演算を行う。
【００３５】
このようにして蛍光の強度分布の更新と評価が繰り返され、例えば、演算された蛍光の強度分布と計測データとの差が許容範囲内と評価されると、断層情報生成部７８は、そのときの蛍光体の濃度分布に基づく吸収係数分布等から光断層情報である蛍光の強度分布を生成し、断層画像生成部８０は、この光断層情報に基づいて光断層画像を生成する。
【００３６】
このように、画像処理部１４は、計測部１２から読み込んだ計測データに対して、所定のデータ処理を行った後、この処理結果に基づいた画像処理を行うことにより、計測データに基づいた検体１８の光断層画像を再構成する。
【００３７】
ここで、本実施の形態に係る光断層情報生成装置１００で使用される光強度分布算出方法について、図面および数式を交えながら理論的な説明を行う。
【００３８】
本実施の形態では、図４（Ａ）および図４（Ｂ）に示すように、光断層情報生成装置１００の光源ヘッド２４から発する光を励起光として、この励起光が照射されることにより蛍光を発する蛍光体６２を含む物質ないし薬剤が検体１８に投与される。光断層情報生成装置１００は、検体１８の断層画像として、蛍光体６２の分布を含む画像を再構成し、検体１８内での各種臓器に対する蛍光体６２の分布が視認可能となるようにする。
【００３９】
具体的には、図４（Ａ）に示されるように、検体１８に励起光を照射すると、この励起光が検体１８内で散乱しながら蛍光体６２に達する。これにより、検体１８内の蛍光体６２が発光する。また、図４（Ｂ）に示されるように、蛍光体６２から発せられる蛍光は、検体１８内で散乱を繰り返しながら、検体１８から射出される。このとき、蛍光と共に励起光も射出されるが、受光ヘッド２６は、その前段に図示しない励起光カットフィルタを有しているため、検体１８から射出される蛍光のみを検出光として、この検出光（以下、蛍光とする）の強度を検出する。光断層情報生成装置１００は、この検出光の光強度分布から検体１８内での蛍光体６２の分布（濃度分布）を得る。
【００４０】
ここで、励起光などの光を検体１８に照射した場合、照射位置近傍の領域では、光に対する屈折率が方向によって異なるなどの異方散乱領域となっているが、検体１８内で所定距離以上に進行すると等方散乱領域となる。
【００４１】
検体１８内を散乱する光は、エネルギーを輸送する粒子と見なされることから、光強度の分布は、光強度のエネルギー保存式である光輸送方程式を用いて表すことができる。しかし、この光輸送方程式を解くことは現状では困難とされている。
【００４２】
一方、検体１８では、一般的に異方散乱領域が数ｍｍ程度であるために、数ｃｍ以上の大きさの検体１８では、その体内を実質的に等方散乱領域と見なすことができる。すなわち、検体１８での光の散乱は等方散乱として近似することができる。
【００４３】
ここから、光拡散方程式を用いることにより、光強度の分布が得られる。この光拡散方程式は、次式（１）で表される。ここで、Φ（ｒ、ｔ）は検体１８内の光密度、Ｄ（ｒ）は拡散係数分布、μａ（ｒ）は吸収係数分布、ｑ（ｒ、ｔ）は光源の光密度を表し、ｒは検体１８内の座標位置、ｔは時間を表す。
【００４４】
【数１】
【００４５】
ここで時間的に連続した光であれば、光拡散方程式は次式（２）で示すことができる。
【００４６】
【数２】
【００４７】
光学特性値である拡散係数分布Ｄ（ｒ）および吸収係数分布μａ（ｒ）が既知であるときに、この光拡散方程式（２）を用いて検体１８から射出される光強度分布を求める場合、順問題として計算することができる。しかし、光強度分布が既知であり、ここから、光拡散方程式を用いて検体１８の光学特性値である拡散係数分布Ｄ（ｒ）および吸収係数分布μａ（ｒ）を求めることは、逆問題計算となる。
【００４８】
ここで、検体１８での拡散係数分布Ｄ（ｒ）及び吸収係数分布μａ（ｒ）は、光の波長によって異なるため、励起光の波長に対する拡散係数分布をＤｓ（ｒ）、吸収係数分布をμａｓ（ｒ）とし、光源の光密度をｑｓ（ｒ）とすると、励起光に対する光拡散方程式は次式（３）で表され、一方、蛍光の波長に対する拡散係数分布をＤｍ（ｒ）、吸収係数分布をμａｍ（ｒ）とし、蛍光の光源をｑｍ（ｒ）とすると、蛍光に対する光拡散方程式は次式（４）で表される。また、蛍光の光源ｑｍ（ｒ）は、検体１８内の光密度Φｓ（ｒ）、蛍光体６２の量子効率γ、モル吸光係数εおよび濃度分布Ｎ（ｒ）を用いて、ｑｍ＝γ・ε・Ｎ（ｒ）・Φｓ（ｒ）と表すことができる。したがって、式（４）は次式（５）で置き換えられる。
【００４９】
【数３】
【００５０】
一方、図５に破線で示されるように、検体１８内のヘモグロビンは約７００ｎｍ以下の波長の光に対して吸収が強く、また、図５に二点鎖線で示されるように、検体１８内の水分は波長が約１μｍ以上の光に対して吸収が強い。換言すれば、７００ｎｍ〜１μｍの波長域は検体１８における吸収が弱いということになり、この帯域は所謂光学的窓と呼ばれている。したがって、検体１８は、光学的窓に相当したものとなっている。この波長域では、検体１８の吸収係数分布μａは、０．００２ｍｍ−１〜０．１ｍｍ−１の範囲となる。
【００５１】
また、検体１８内での光の散乱係数（散乱の強さ、図５に実線で示す）は、波長が長くなると小さくなるが、その変化は緩やかであり、光学的窓となっている７００ｎｍ〜１μｍの波長域では、散乱の強さを波長に依存せず略一定と見なすことができる。
【００５２】
そのため、本実施の形態に係る光断層情報生成装置１００は、光源ヘッド２４から発する励起光として、生体（検体１８）での光学的窓に該当する７００ｎｍ〜１μｍの波長の赤外線（近赤外線）を用いる。これにより、検体１８での光学的特性である吸収係数μａ、散乱係数（拡散係数Ｄ）を略一定の値（既知の値）とすることができ、式（３）及び式（５）において、Ｄｓ（ｒ）＝Ｄｍ（ｒ）＝Ｄ（ｒ）、μａｓ（ｒ）＝μａ（ｒ）＋ε・Ｎ（ｒ）、μａｍ（ｒ）＝μａ（ｒ）と簡略化することができる。ここで、εはモル吸光係数、Ｎ（ｒ）は蛍光体６２の濃度分布を表し、ε・Ｎ（ｒ）は蛍光体６２による励起光の吸収度合いを表す。したがって、式（３）は次式（６）で置き換えられ、式（５）は次式（７）で置き換えられる。
【００５３】
【数４】
【００５４】
また、本実施の形態で断層画像を観察する検体１８には、この近赤外線を励起光として発光する蛍光体６２を含む物質ないし薬剤が投与されている。このとき、上記式（７）に示されるように、蛍光体６２が光源となるときの蛍光の強度は、励起光の強度Φｓ（ｒ）に基づくもの、すなわちΦｓ（ｒ）の関数である。これは、拡散係数分布や吸収係数分布を予め設定すれば既知であり、光源の光強度ｑｓ（ｒ）も既知であることから、有限要素法などの数値解析手法により検体１８内の光強度Φｓ（ｒ）を順問題として求めることができる。
【００５５】
したがって、計測部１２で蛍光強度を計測し、この計測データに基づき蛍光強度分布Φｍ（ｒ）を求め、一つ（１系統）の逆問題計算によって検体１８内での蛍光体６２の濃度分布Ｎ（ｒ）を求めることができる。計測部１２で励起光強度を計測する必要はない。
【００５６】
得られた蛍光体６２の位置ｒにおける濃度分布Ｎ（ｒ）を検体１８の位置ｒにおける断層画像に合成することにより、検体１８内における蛍光体６２の濃度分布を光断層画像上で視認することができる。
【００５７】
次いで、上記理論に基づいて、本実施の形態に係る光断層情報生成装置１００が具体的にどのような演算処理を行っているかを説明する。
【００５８】
光断層情報生成装置１００では、計測部１２の計測ユニット２０に検体１８が配置されると、光源ヘッド２４から予め設定された波長の近赤外光を励起光として検体１８に照射する。この励起光は、検体１８内を拡散しながら伝播（透過）する。
【００５９】
ここで、検体１８に蛍光体６２が投与されていると、この蛍光体６２に励起光が照射され、これにより、蛍光体６２が発光する。この蛍光は、検体１８内を拡散しながら伝播して、検体１８から周囲へ射出される。
【００６０】
計測ユニット２０には、検体１８を囲うように受光ヘッド２６が所定の角度間隔で配列されており、計測部１２では、検体１８から射出される蛍光を検出光として、各受光ヘッド２６で受光する。
【００６１】
また、計測部１２では、機枠２２を回転することにより検体１８への励起光の照射位置及び検体１８から射出される蛍光の検出位置を相対的に変えて、励起光の照射、検出光の受光を繰り返す。これにより、検体１８の周囲に沿って照射した励起光に応じた蛍光の強度の測定データが得られる。
【００６２】
光断層情報生成装置１００の画像処理部１４では、この計測データに基づいて蛍光体６２の濃度分布Ｎ（ｒ）を演算する。
【００６３】
図６は、画像処理部１４が行う蛍光体６２の濃度分布計算の処理手順を示したフローチャートである。前述の通り、蛍光体の濃度分布を算出することは蛍光強度分布を算出することに相当するので、この図は、本実施の形態に係る光強度分布算出方法の手順を示していると言える。
【００６４】
なお、画像処理部１４では、検体１８の光学的特性に基づいて予め設定されている波長の近赤外線を励起光としており、ここから、吸収係数分布μａ（ｒ）および拡散係数分布Ｄ（ｒ）の値が予め設定されて記憶されている。なお、この吸収係数分布μａ（ｒ）および拡散係数分布Ｄ（ｒ）は、検体１８の個体差に合わせて適宜入力される値であっても良い。
【００６５】
このフローチャートでは、ステップ（以下ＳＴと略す）１０００で、計測部１２の複数の受光ヘッド２６によって得られた実測値、すなわち、検体１８から射出される蛍光強度分布Φｍ（ｒ）measがメモリ５１に読み込まれる。
【００６６】
一方、ＳＴ１０２０では、予め設定されている吸収係数分布μａ（ｒ）及び拡散係数分布Ｄ（ｒ）が用いられ、検体１８内に蛍光体が存在しなかった場合、すなわち蛍光体による励起光の吸収が生じていない場合の励起光強度分布Φｔ（ｒ）calcが下記光拡散方程式（８）に従って算出される。
【００６７】
【数５】
【００６８】
ＳＴ１０４０では、蛍光体６２の濃度分布Ｎ（ｒ）calcの初期値が設定される。ＳＴ１０６０では、ＳＴ１０４０において設定された蛍光体の濃度分布Ｎ（ｒ）calcと、予め設定されている吸収係数分布μａ（ｒ）及び拡散係数分布Ｄ（ｒ）とを用いて、既に説明した式（６）を適宜修正した下記光拡散方程式（９）に従って、励起光強度分布Φｓ（ｒ）calcが計算される。
【００６９】
【数６】
【００７０】
ＳＴ１０７０では、下記式（１０）に示すように蛍光体による励起光の吸収が生じていない場合の励起光強度分布Φｔ（ｒ）calcから蛍光体による励起光の吸収が生じている場合の励起光強度分布Φｓ（ｒ）calcを減じ、さらに励起光を蛍光に換算する係数γを乗ずることにより、検体１８から射出される蛍光強度の予測値である計算蛍光強度分布Φｍ（ｒ）calcが計算される。
【００７１】
【数７】
【００７２】
なお、ＳＴ１０２０における励起光強度分布Φｔ（ｒ）calcの算出やＳＴ１０６０における励起光強度分布Φｍ（ｒ）calcの算出は、数学的モデルである光拡散方程式を有限要素法などの数値解析手法を用いた公知の順問題計算として比較的容易に演算することができる。
【００７３】
ＳＴ１０８０では、計測データに基づいた蛍光強度分布Φｍ（ｒ）measと、演算結果に基づいた蛍光強度分布Φｍ（ｒ）calcとを比較し、ＳＴ１１００では、蛍光強度分布Φｍ（ｒ）measと蛍光強度分布Φｍ（ｒ）calcとの間の差分が許容範囲内か否か、具体的には予め設定した所定値以内か否かを判断する。
【００７４】
ＳＴ１０８０において、蛍光強度分布Φｍ（ｒ）measと蛍光強度分布Φｍ（ｒ）calcとの間の差分が所定値以内でないと判断されるときには、ＳＴ１１００で否定判定されてＳＴ１１２０へ移行する。
【００７５】
ＳＴ１１２０では、ヤコビ行列（Jacobian matrix）を用いた公知の手法で検体１８の光学特性値の変化に対する光強度分布の変化を演算する。
【００７６】
ＳＴ１１４０では、以下の光拡散方程式（１１）で表される逆問題に対して後述の最適化手法を適用することにより、蛍光体６２の濃度分布の推定値Ｎ（ｒ）estを求める。
【００７７】
【数８】
【００７８】
より詳細には、ＳＴ１１４０ではまず、蛍光強度分布Φｍ（ｒ）measと蛍光強度分布Φｍ（ｒ）calcとの間の誤差（例えば、二乗誤差ｙ）を評価する。すなわち、二乗誤差ｙは次式（１２）から得られ、この二乗誤差ｙを評価する。
【００７９】
【数９】
【００８０】
次に、このＳＴ１１４０では、この二乗誤差ｙを最小とする蛍光体６２での蛍光の吸収εＮ、すなわち、蛍光体６２の濃度分布Ｎ（ｒ）estを後述の最適化手法により推定する。
【００８１】
このようにして濃度分布の推定値Ｎ（ｒ）estが求まるが、これは、上記式（１１）のＮ（ｒ）calcに対する推定値を求めたに過ぎない。そこで、ＳＴ１１６０では、この濃度分布の推定値Ｎ（ｒ）estによりＮ（ｒ）calcを更新し、ＳＴ１０６０へ戻る。
【００８２】
画像処理部１４は、このＳＴ１０６０からＳＴ１１６０を繰返し行う。これにより、蛍光強度分布Φｍ（ｒ）measと蛍光強度分布Φｍ（ｒ）calcとの差が所定値以内と判断されると、ＳＴ１１００において肯定（ＹＥＳ）判定してＳＴ１１８０へ移行し、ＳＴ１１８０ではこのとき最終的に設定されている濃度分布Ｎ（ｒ）calcを計測データから得られた濃度分布Ｎ（ｒ）として出力する。
【００８３】
なお、ここでは、計測データに基づいた蛍光強度分布Φｍ（ｒ）measを取得した後、蛍光体６２の濃度分布Ｎ（ｒ）calcの初期値を設定するようにしているが、処理手順はこれに限らず、蛍光体６２の濃度分布Ｎ（ｒ）calcの初期値を設定した後に、蛍光強度分布Φｍ（ｒ）measを取得するようにしても良い。他の手順についても、当該フローの目的を逸脱しない範囲で適宜ステップの順番を変更することができる。
【００８４】
本発明では、ＳＴ１１４０で使用する逆問題に対する最適化手法としてＡＲＴ（Algebraic Reconstruction Technique：代数的再構成技術）法を採用する。しかし、ＡＲＴ法は、収束は早いがノイズに対する耐性（ノイズ耐性）に乏しいという特徴がある。そこで、本願発明者は、以下に示す改良型ＡＲＴ法を見い出した。なお、ノイズ耐性とは、ノイズが多い環境下においても所定時間内に演算を収束させることができる性質、またはノイズが多い環境下においても所定レベル以上の演算精度を維持できる性質のことである。
【００８５】
具体的に数式を交えて説明すると、上記逆問題を解くためには、式（１２）における二乗誤差ｙを最小にする必要がある。よって、この計算を行うために、まず式（１２）をεＮで微分して簡略化（Born近似）すると次式（１３）が得られる。
【００８６】
【数１０】
【００８７】
ここで、ΔｘはεＮの微小変化、具体的には、本実施の形態ではＡＲＴ法という繰り返し演算（反復法）を用いているためεＮの更新値を表している。また、ヤコビ行列Ｊは、光学特性値の変化に対する光強度分布の変化の割合を示している。なお、光強度分布は、光源の座標位置と検出地点の座標位置との関数として表される。
【００８８】
式（１３）においてΔｘが求めたい解であるため、式（１３）をΔｘについて書き直すと次式（１４）となる。
【００８９】
【数１１】
【００９０】
原理的には当該式（１４）を用いて計算を進めることができるはずであるが、実際には式（１４）をそのまま使用して計算機上で計算すると次のような問題が生じる。すなわち、ヤコビ行列Ｊ中に著しく小さい成分が存在する場合に逆行列Ｊ−１においては逆数となって現れるため、この成分がノイズのような小さな誤差を含んでいた場合にその誤差が拡大されて計算精度に著しく悪影響を与えてしまう。ＡＲＴ法がノイズ耐性に乏しい理由はこのためである。
【００９１】
そこで本願では、まず行列分解の一手法である特異値分解法を用いてヤコビ行列Ｊを分解して次式（１５）を得る。
【００９２】
【数１２】
【００９３】
ここで、Ｄはヤコビ行列の特異値をその対角成分とするｍ×ｎ行の対角行列、Ｕはｍ×ｍ行の正規直交行列、ＶＴはｎ×ｎ行の転置行列である。ＶＴはＶの転置行列である。ＵおよびＶには、ＵＴＵ＝ＶＴＶ＝Ｉ（単位行列）の関係がある。
【００９４】
式（１４）に代入できるように、式（１５）をＪ−１を求める形に書き換えると次式（１６）となる。
【００９５】
【数１３】
【００９６】
さらに、対角行列Ｄとは別個に新たな単位対角行列Ｈを導入し、Ｄにおいて特異値が閾値αth以下となっている対角成分を特定し、Ｈにおける同位置の成分を０で置換する（特異値が閾値αthよりも大きい成分は１を維持する）。このＨを用いて、式（１４）に対し次式（１７）で表されるΔｙの置換を行う。この置換において、０に置換した成分は閾値αth以下のもののみであるので、Δｙ全体の大きさ（絶対値）に大きな変化はないことが見込まれる。
【００９７】
【数１４】
【００９８】
そして、上記式（１４）、式（１６）および式（１７）から次式（１８）を得る。
【００９９】
【数１５】
【０１００】
式（１８）右辺のＤ−１Ｈの部分は、小さい値の成分については０との積となっており、また、他の対角成分についてはＤ−１の値をそのまま維持したものとなっている。すなわち、本願に係る改良型ＡＲＴ法により、ノイズの影響を減少させつつ蛍光体濃度分布の計算を行うことができる。
【０１０１】
本発明ではＡＲＴ法を用いているので、具体的には、式（１８）の代わりに次式（１９）のＡＲＴの逐次近似式を計算することになる。よって、式（１７）を用いたΔｙの置換は、式（１９）に対して行われる。ここで、λは緩和係数と呼ばれ、解の収束性を制御するため計算上使用されるパラメータである。
【０１０２】
【数１６】
【０１０３】
図７は、上記計算を実現する本実施の形態に係る更新処理部７６のより詳細な構成を示した機能ブロック図である。更新処理部７６は、ヤコビ行列算出部１５１、特異値分解部１５２、単位対角行列取得部１５３および逐次近似部１５４を備える。
【０１０４】
ヤコビ行列算出部１５１は、図３に示した評価部７４から出力される評価結果信号Ｓ１１が否定（ＮＯ）を示す場合、図３に示した演算処理部７２から出力されるデータ信号Ｓ１２からヤコビ行列Ｊを求め、特異値分解部１５２へこれを出力する。ここで、データ信号Ｓ１２は、検体１８の光学特性値、図１および図２に示した光源ヘッド２４の座標位置および受光ヘッド２６の位置情報を含み、ヤコビ行列算出部１５１は、このデータ信号Ｓ１２に基づいて上記式（１３）におけるヤコビ行列Ｊ、より具体的には上記式（１９）におけるヤコビ行列Ｊを求める。
【０１０５】
特異値分解部１５２は、ヤコビ行列算出部１５１で求められたヤコビ行列Ｊを上記式（１５）に従って特異値分解し、得られた対角行列Ｄを単位対角行列取得部１５３へ出力し、同様に得られた正規直交行列Ｕおよび転置行列ＶＴを逐次近似部１５４へ出力する。
【０１０６】
単位対角行列取得部１５３は、特異値分解部１５２から出力された対角行列Ｄに基づいて、既に説明した本願の手法によって単位対角行列Ｈを求め、逐次近似部１５４へ出力する。
【０１０７】
逐次近似部１５４は、特異値分解部１５２から出力された正規直交行列Ｕおよび転置行列ＶＴ、並びに単位対角行列取得部１５３から出力された単位対角行列Ｈを用いて上記式（１９）の逐次近似式を計算し、εＮの更新値であるΔｘを求めて図３に示した演算処理部７２へ出力する（信号Ｓ１３）。
【０１０８】
次いで、本実施の形態に係る光断層情報生成装置１００が奏する効果を実例を紹介しながら説明する。
【０１０９】
図８は、実際に腫瘍を有するヌードマウスに対し抗体を投与し、このマウスを２次元カメラで腹部方向から撮影して得られた画像である。この画像から、マウスの右腹部に腫瘍が存在することが確認できる。なお、本願における動物実験は、出願人である富士フイルム株式会社内の動物倫理規則および安全性評価センター動物実験規程に則って行われたものである。
【０１１０】
図９（Ｂ）は、光断層情報生成装置１００を用いて、この抗体投与マウスから得られた再構成画像の例である。図９（Ａ）に示した従来のＡＲＴ法による再構成画像と比較すると、従来のＡＲＴ法ではマウス体内に腫瘍の存在を認めることはできないが、本実施の形態に係る光断層情報生成装置１００で得られた再構成画像では、マウス右腹部に濃度分布の偏り、すなわち腫瘍の影が認められることがわかる。
【０１１１】
以上説明したように、本発明の実施の形態１に係る光断層情報生成装置１００によれば、励起光が照射された検体（生体）から射出される蛍光の強度を検出し、蛍光の強度（蛍光強度）の測定データを取得する。また、検体の光学特性値である散乱係数及び吸収係数の各々の分布を予め設定すると共に、蛍光体の濃度分布に基づく蛍光体の吸収係数分布を仮設定しておき、数学的モデルから蛍光の強度を取得する。この数学的モデルから取得した蛍光の強度と、測定データとして得られる蛍光の強度を比較し、この差分に対する評価を行う。このとき、差分が所定レベルより大きければ、逆問題計算を行うことにより、この差分を減少させるように吸収係数分布を推定し、推定した吸収係数分布に基づいた蛍光の強度と測定データを比較する処理を反復することにより、蛍光体の濃度分布に基づく吸収係数分布を決定する。このようにして得られる蛍光体の濃度分布に基づく吸収係数分布から検体内の蛍光体の濃度分布を求めて、検体の光断層画像を形成するための光断層情報を生成する。
【０１１２】
これにより、光断層情報を取得するための光拡散方程式を用いた逆問題計算は、蛍光の強度に対して行うのみでよいために、断層情報を生成するための処理負荷が軽減され、処理時間の短縮を図ることができる。また、断層画像を形成するための光の強度測定は、検体から射出される蛍光に対して行うのみでよいため、装置の簡略化、計測時間の短縮を図ることができる。例えば、受光ヘッド２６は蛍光のみを計測できる機能があれば充分で、励起光を計測できる必要はない。
【０１１３】
これに対して、従来構成では、散乱係数Ｄ（ｒ）及び吸収係数分布μａ（ｒ）を既知の値として設定しないため、測定データとして蛍光強度分布に加え励起光強度分布が必要となり、蛍光と励起光の双方に対応した受光ヘッドが必要となる。また、このときには、励起光に対する光拡散方程式、および蛍光に対する光拡散方程式のそれぞれに対して逆問題計算を行い、蛍光体の濃度分布を得る必要がある。従って、各々の逆問題計算を処理する装置構成が必要となり、装置が大規模化する。
【０１１４】
また、以上の構成において、前記光源を前記検体の周囲に沿って相対移動することにより前記励起光の前記検体への照射位置を変更しながら、それぞれの照射位置で前記蛍光の強度を検出する。このとき、前記蛍光を検出する複数の検出手段を前記検体の周囲に予め設定された間隔で設け、それぞれの検出手段を前記光源の前記検体に対する相対移動と対になって移動して前記蛍光の強度を検出する構成を取り得る。
【０１１５】
これにより、光断層画像生成装置における蛍光の強度を検出するための構成が簡略化できる。
【０１１６】
また、以上の構成において、本実施の形態１に係る光断層情報生成装置１００は、光強度分布を算出する際の逆問題においてＡＲＴ（代数的再構成技術）を使用し、このＡＲＴの逐次近似式中に現れるヤコビ行列Ｊに対して特異値分解を行って対角行列Ｄを得、このＤにおいて閾値αth以下となっている特異値を０で置換した単位対角行列ＨをＤに乗ずることにより、逆問題の計算負荷を減ずると共に、処理時間の短縮化を実現している。
【０１１７】
また、光断層情報生成装置１００は、ＡＲＴで逆問題を解く計算工程において、光学特性値の変化に対する光強度分布の変化の割合を複数の独立な要素に分解し、そのうちノイズの影響が支配的な要素を特定して除去してから逆問題を解くようにしている。これにより、ノイズが多く含まれている環境下においても、逆問題計算を迅速に収束させることができ、さらに再構成画像の精度を向上させることができる。なお、付言するならば、上記計算において、ノイズが支配的でない要素については値を維持したままとしているので、ノイズの影響を除去する過程において正味の信号成分まで除去してしまうことを防止しており、この点からも再構成画像の精度向上を実現できている。
【０１１８】
また、以上の構成において、このような、単位対角行列Ｈの導入は、ＡＲＴのような繰り返し近似計算法において計算を単純化簡易化することができ、実装置への実装化やプログラミングのし易さに繋がっている。
【０１１９】
また、以上の構成において、対角行列Ｄから単位対角行列Ｈを生成する際に、Ｄにおいて特異値が閾値αth以下となっている対角成分を０で置換することを説明した。この閾値αthとして妥当な値を設定することにより、本実施の形態に係る光強度分布の算出精度、ひいては光断層情報生成装置１００で得られる再構成画像の精度をより向上させることができる。この閾値αthの具体的な設定方法について以下説明する。
【０１２０】
図１０は、上記ヤコビ行列Ｊのランク数（横軸）と対角行列Ｄにおける特異値の値（縦軸）との関係を特性曲線として示した図である。この図からわかるように、Ｄにおける特異値は、ヤコビ行列のランク数が増加するのに伴いステップ状に減少する性質があることがわかった。
【０１２１】
そこで、本願発明者は、対角行列Ｄから単位対角行列Ｈを生成する際に使用する閾値αthとして、このステップ状に変化する箇所、即ち、その他の箇所と比較して急激に特異値が減少する箇所（図では破線で示す）の特異値の値を閾値αthとして用いることを思い付いた。閾値αth以下の特異値を０で置換するということは、ノイズの影響を除去する過程において正味の信号成分までも過大に除去してしまうおそれがある。すなわち、ノイズの影響除去と正味の信号成分の維持はトレードオフの関係にある。しかし、特異値が急激に変化する箇所を閾値αthに設定すれば、信号成分を過大に除去せずにノイズの影響を除去できる可能性が高まることに本願発明者は気付いたのである。具体的には、特異値が１０−３、１０−４、１０−５、１０−６の値を示す箇所の近傍において特異値がステップ状に変化している。そこで、本実施の形態に係る光断層情報生成装置では、１０−３、１０−４、１０−５、１０−６を閾値αthとして用いることにする。また、特異値が１０−８を示す箇所は厳密にはステップ状になっていないが、他の閾値との比較のためにこの値も閾値αthとして追加して用いることにする。
【０１２２】
なお、図１０に示した特性曲線は、光断層情報生成装置のシステム精度や計算環境によって変化するものである。従って、上記の閾値αthの具体的な値は実施の形態の一例にすぎない。
【０１２３】
図１１（Ａ）〜（Ｆ）は、実際の抗体投与マウスを用いて再構成画像の精度を検証した結果である。閾値αthとして上記の値を用いた場合に実際の再構成画像の精度がどのように変化するかをこの図から判断することができる。
【０１２４】
図１１（Ａ）は従来のＡＲＴ法の場合、図１１（Ｂ）は閾値αthとして１０−３を用いた場合、図１１（Ｃ）は閾値αthとして１０−４を用いた場合、図１１（Ｄ）は閾値αthとして１０−５を用いた場合、図１１（Ｅ）は閾値αthとして１０−６を用いた場合、図１１（Ｆ）は閾値αthとして１０−８を用いた場合である。この図からわかるように、閾値αthとして１０−４または１０−５を用いた場合に、腫瘍と思われる領域を他の領域と区別し得る精度の再構成画像を得ることができる。従って、本実施の形態に係る光断層情報生成装置が使用する閾値αthとしては１０−４または１０−５を用いる。また、閾値αthとして１０−５を用いた場合、腫瘍と思われる領域が他の領域と比べて識別可能な精度になるだけでなく、腫瘍の大きさについても一定の精度で再現できていることがわかる。よって、閾値αthとして１０−５を用いることが更に望ましい。
【０１２５】
このように、本実施の形態に係る光断層情報生成装置では、対角行列Ｄから単位対角行列Ｈを生成する際に使用する閾値αthの候補として、ヤコビ行列Ｊのランク数と対角行列Ｄにおける特異値の値との関係を示す特性曲線がステップ状に変化する複数の箇所の値を使用する。そして、複数の閾値αthの候補の中から最適な閾値αthを適宜選択する。これにより、本実施の形態に係る光強度分布の算出精度、ひいては光断層情報生成装置１００で得られる再構成画像の精度をより向上させることができる。
【０１２６】
なお、本実施の形態では、閾値αthとして、ヤコビ行列Ｊのランク数と対角行列Ｄにおける特異値の値との関係を示す特性曲線がステップ状に変化する箇所の値を使用する場合を例にとって説明したが、例えば、よりノイズの影響を除去したい場合には、上記設定方法で求めた閾値αthに対しオフセットを加えた値を閾値として用いても良い。
【０１２７】
また、本実施の形態では、光学特性値の変化に対する光強度分布の変化の割合を示すヤコビ行列の各成分のうち、ノイズの影響が大きく現れ得る成分を０に置換する場合、すなわち当該成分を完全に削除する場合を例にとって説明したが、この成分には正味の信号成分も含まれているはずであるから、完全に削除するのではなく、逆数をとったときに所定値以下の大きさに収まるような定数で置換するような構成としても良い。これにより、計算負荷等を軽減しつつ再構成画像の精度を向上させることができる。
【０１２８】
（実施の形態２）
実施の形態１では、検体から発せられる蛍光のみを実測して１系統の逆問題を解く態様を例にとって説明したが、本実施の形態２では、検体から発せられる蛍光だけでなく励起光も計測し、２系統の逆問題を解く態様について説明する。
【０１２９】
図１２は、本発明の実施の形態２に係る光断層情報生成装置２００内の計測部２１０および画像処理部１４ａの主要な構成を示すブロック図である。なお、実施の形態１で示した光断層情報生成装置１００と同一の構成要素には同一の符号を付し、その説明を省略する。また、基本的動作は同一であるが、詳細な点で違いがある構成要素には、同一の番号にアルファベットの小文字を付した符号を付して適宜説明を加える（以降の図の説明においても同様の表記方法を用いる）。例えば、本実施の形態２に係る画像処理部１４ａは、実施の形態１で示した画像処理部１４と基本的構成は同一であるが詳細な動作の点においては違いがあるため、同一の番号１４にアルファベットａを付している。
【０１３０】
光断層情報生成装置２００は、実施の形態１で示した光断層情報生成装置１００の構成に更に新規な構成である受光ヘッド２１１および増幅器（Ａｍｐ）２１２を有している。受光ヘッド２４は、実施の形態１で説明した通り蛍光Ｌ１用の受光ヘッドであるが、受光ヘッド２１１は、励起光Ｌ２用の受光ヘッドである。また、増幅器２１２は、受光ヘッド２１１用の増幅器であり、受光ヘッド２１１から出力される電気信号を増幅し、Ａ／Ｄ変換器４０へ出力する。
【０１３１】
ここで、受光ヘッド２１１〜Ａ／Ｄ変換器４０の経路は新たに設けられた経路であり、光断層情報生成装置１００における受光ヘッド２６〜Ａ／Ｄ変換器４０の経路は光断層情報生成装置２００においてもそのまま維持されている。即ち、光断層情報生成装置２００における受光ヘッド２６および増幅器３８の数は光断層情報生成装置１００のそれらと同数である。また、励起光用の受光ヘッド２１１および増幅器２１２の数も、蛍光用の受光ヘッド２６および増幅器３８のそれらと同数である。
【０１３２】
図１３は、受光ヘッド２６および受光ヘッド２１１の関係、およびこれらの構成の機能を説明するための図である。
【０１３３】
この図からわかるように、光断層情報生成装置２００の計測ユニット２０ａの基本構成は、実施の形態１に係る光断層情報生成装置１００の計測ユニット２０と類似している。しかし、検体１８から射出された光が受光ヘッド２６に到達するまでの経路においてダイクロイックミラー２１３が設けられ、検体１８から射出される光のうち所定の短波長成分のみを反射し、受光ヘッド２１１へ入射するように構成されている。検体１８から射出される光は、蛍光Ｌ１の成分のみならず励起光Ｌ２の成分も含む。そこで、ダイクロイックミラー２１３は、蛍光Ｌ１に比べ短波長である励起光Ｌ２の成分のみを反射するように構成され、受光ヘッド２１１には励起光成分のみが入射される。また、実施の形態１で説明した通り、受光ヘッド２６の前段には図示しない励起光カットフィルタが設けられているため、受光ヘッド２６には蛍光Ｌ１のみが入射される。
【０１３４】
図１４は、画像処理部１４ａが行う蛍光体６２の濃度分布計算の処理手順、すなわち本実施の形態に係る光強度分布算出方法の手順を示したフローチャートである。
【０１３５】
前述の通り、実施の形態１で測定されるのは蛍光のみであったが、本実施の形態２では蛍光に加え励起光も測定する。そのため、ＳＴ２０００では、検体１８から射出される蛍光強度分布の測定データΦｍ（ｒ）measがメモリ５１に読み込まれると共に、ＳＴ２０２０において、検体１８から射出される励起光強度分布の測定データΦｘ（ｒ）measもメモリ５１に読み込まれる。なお、ＳＴ２０００は、実施の形態１の図６で示したＳＴ１０００と同一の処理である。
【０１３６】
ＳＴ２０４０では、ＳＴ２０００で読み込まれた蛍光強度分布Φｍ（ｒ）measとＳＴ２０２０で読み込まれた励起光強度分布Φｘ（ｒ）measとから、次式（２０）に基づいて、蛍光体が存在しなかった場合の励起光強度分布Φｔ（ｒ）measが計算される。
【０１３７】
【数１７】
【０１３８】
ＳＴ２０６０では、ＳＴ２０８０以降のループ計算の準備として、蛍光体が存在する場合の励起光波長における吸収係数分布μａｘ（ｒ）の初期値、および蛍光体が存在しなかった場合の励起光波長における吸収係数分布μａｔ（ｒ）の初期値が設定される。なお、本実施の形態でも、実施の形態１と同様に拡散係数分布Ｄ（ｒ）は略一定の既知の値とする。
【０１３９】
ＳＴ２０８０では、下記光拡散方程式（２１）で示される順問題が解かれることにより、蛍光体が存在する場合に検体１８から射出される励起光強度分布Φｘ（ｒ）calcが求められ、同様に、下記光拡散方程式（２２）で示される順問題が解かれることより、蛍光体が存在しなかった場合に検体１８から射出される励起光強度分布Φｔ（ｒ）calcが求められる。
【０１４０】
【数１８】
【０１４１】
ＳＴ２１００では、蛍光体が存在する場合の励起光強度分布の計算結果（予想値）であるΦｘ（ｒ）calcと実測値であるΦｘ（ｒ）measとが比較され、また、蛍光体が存在しない場合の励起光強度分布の計算結果（予測値）であるΦｔ（ｒ）calcと実測値であるΦｔ（ｒ）measとが比較される。具体的には、比較結果として、各予測値と各実測値との差分が求められる。
【０１４２】
ＳＴ２１２０〜ＳＴ２１８０の計算は、実施の形態１と同様のループ計算の一部である。すなわち、ＳＴ２１２０において、ＳＴ２１００で求められた各予測値と各実測値との間の差分が許容範囲内か否か、具体的には予め設定した所定値以内か否かが判断され、差分が所定値以内でないと判断されるときには、ＳＴ２１４０〜ＳＴ２１８０の計算を経た後、ＳＴ２０８０に戻るループ計算である。ここで、実施の形態１と大きく異なる点は、ＳＴ２１６０において行われる逆問題計算である。実施の形態１では、検体から発せられる蛍光のみを実測して１系統の逆問題を解いていたが、本実施の形態２では、検体から発せられる蛍光だけでなく励起光も計測し、２系統の逆問題を解くからである。
【０１４３】
なお、ＳＴ２１００で求まる差の値は２個あるため、厳密には本実施の形態２に係るＳＴ２１２０〜ＳＴ２１８０の計算は、実施の形態１に係るＳＴ１１００〜１１６０の計算とは異なる。しかし差が２個あったとしても、差を所定値以内に収めようとする実施の形態１における繰り返し演算の趣旨は、本実施の形態でも同様であるため、具体的な繰り返し計算方法をどのように規定するかは当業者が適宜決定すれば良いことである。そこで、ＳＴ２１２０では説明を簡単にするため、あたかも差が１個であるが如く説明を行っている。また、ＳＴ２１６０およびＳＴ２１８０では、取り扱われるパラメータが実施の形態１の濃度分布と異なり吸収係数分布である。これは、ループ計算の初期値としてＳＴ２０６０において設定し、ＳＴ２１８０において更新するパラメータが、実施の形態１と異なり吸収係数分布だからである。
【０１４４】
ＳＴ２１２０において各予測値と各実測値との間の差分が所定値以内と判断されるとＳＴ２２００へ移行し、ＳＴ２２００ではこのとき最終的に設定されている吸収係数分布μａｘ（ｒ）および吸収係数分布μａｔ（ｒ）から次式（２３）に従って濃度分布Ｎ（ｒ）を算出し、本フローチャートは終了する。なお、吸収係数分布と濃度分布との間には次式（２４）および（２５）の関係があり、式（２３）はこれらの関係から求められたものである。
【０１４５】
【数１９】
【０１４６】
以上説明したように、本発明の実施の形態２に係る光断層情報生成装置２００によれば、検体から発せられる蛍光だけでなく励起光も計測し、２系統の逆問題を解くために、装置構成が実施の形態１よりは大規模化し計算負荷も増大するが、本発明に係る改良型ＡＲＴを用いているので、光強度分布の算出精度、および光断層情報生成装置で得られる再構成画像の精度を向上させることができる。
【０１４７】
なお、以上説明した本発明に係る各実施の形態は、本発明の一例を示すものであり、本発明の構成を限定するものではない。本発明に係る光断層情報生成装置は、上記各実施の形態に限定されず、本発明の目的を逸脱しない範囲で種々変更して実施することが可能である。
【０１４８】
例えば、光拡散方程式は蛍光以外の光であっても同様に適用できるので、光断層情報生成装置１００、２００に限らず、励起光を検体１８に照射し、検体１８から射出される蛍光以外の光を検出し、その光の強度に基づいて断層画像を再構成する任意の構成の光断層情報生成装置の構成としても良い。将来的には、励起光を照射することなく自然発光する光に対して本発明の光断層情報生成装置等を適用することも考えられる。
【０１４９】
また、実施の形態１では、光源ヘッド２４と受光ヘッド２６とが機枠２２と共に機械的に回転する場合を例にとって説明したが、逆問題を１系統で処理するという目的においてはこの構成に限定される必要はなく、例えば、光照射機能と受光機能との双方の機能を有した複数のヘッドを位置固定で配置し、測定中は所定の回転方向に光照射機能と受光機能とを順々に遷移させて蛍光等を測定するような構成としても良い。
【０１５０】
また、実施の形態２では、ダイクロイックミラー２１３を設け、検体１８から射出される蛍光Ｌ１の成分と励起光Ｌ２の成分とを切り分けて、蛍光用の受光ヘッド２６と励起光用の受光ヘッド２１１とにそれぞれ入射する態様を例にとって説明したが、例えば、ダイクロイックミラー２１３を設けずに、受光ヘッド２６と受光ヘッド２１１とをより密に交互に配置することによって受光位置（角度）の違いによる信号差の影響を少なくするような構成としても良い。
【０１５１】
また、本質的には、本発明の適用対象画像は断層画像でなくても良い。光拡散方程式に基づく逆問題を解いて再構成画像を取得するような態様であれば良い。
【０１５２】
また、本発明に係る光強度分布算出方法のアルゴリズムをプログラミング言語によって記述し、必要に応じてコンパイルし、この光強度分布算出プログラムをメモリ（記録媒体）に記憶して他の光断層情報生成装置の情報処理手段によって実行させれば、本発明に係る光断層情報生成装置と同様の機能を実現することができる。
【産業上の利用可能性】
【０１５３】
本発明に係る光断層情報生成装置、光強度分布算出方法および光強度分布算出プログラムは、例えば、励起光によって発光する発光体の分布情報を用いて光断層画像を取得可能とする装置、特に、蛍光トモグラフィ装置の用途に利用することができる。
【符号の説明】
【０１５４】
１００、２００ 光断層情報生成装置
１２、２１０ 計測部
１４、１４ａ 画像処理部
１６ モニタ
１８ 検体
２４ 光源ヘッド
２６、２１１ 受光ヘッド
６２ 蛍光体
７０ 読取部
７２ 演算処理部
７４ 評価部
７６ 更新処理部
７８ 断層情報生成部
８０ 断層画像生成部
【技術分野】
【０００１】
本発明は、光を用いたトモグラフィ（Tomography）技術に係り、特に、励起光を照射して、この励起光によって発光する発光体の分布情報を用いて光断層画像を取得可能とする光断層情報生成装置、光強度分布算出方法および光強度分布算出プログラムに関する。
【背景技術】
【０００２】
生体などの断層画像を取得する方法としては、Ｘ線を用いたＸ線ＣＴ、超音波を用いた超音波ＣＴ、核磁気共鳴を適用したＮＭＲ−ＣＴ、プロトンなどの粒子線を用いた陽子線ＣＴなどがある。また、生体においては、光透過性を有することが知られており、小動物の断層画像に光を用いる光ＣＴ（Computed Tomography）が提案されている（例えば、特許文献１、特許文献２参照）。
【０００３】
生体に照射された光は、生体内で散乱し、散乱した光が生体の周囲から射出される。光ＣＴでは、生体内で乱反射して生体の周囲から射出される光を検出して、電気信号を取得し、それぞれの電気信号に対して所定の信号（画像信号）処理を施して得られる情報から画像再構成を行うことにより、生体の断層画像が得られるようにしたものである。
【０００４】
一方、病理学的な実験分野では、所定の波長の光で発光する蛍光体を含む薬剤などを生体に投与して、生体中での、その薬剤の移動、分布、特定部位へ集積するときの集積過程などを観察するときに光ＣＴ（以下、蛍光ＣＴとする）を用いることができる。すなわち、蛍光体に対する励起光を生体に照射し、この励起光に応じて生体から射出される発光（蛍光）を検出して、生体の二次元的な断層画像または三次元的な断層画像を再構成する。これにより、この断層画像から蛍光体、蛍光体を含む試薬や細胞などの位置、量などの情報を得ることができる。
【０００５】
このような蛍光ＣＴにおいても、生体の表面の一点へ励起光を照射し、これにより生体から射出される散乱蛍光を多点で検出する。これを、励起光の照射位置を変えながら繰り返すことで、照射点数×検出点数だけのデータを取得する。これらのデータの間には、蛍光物質の分布、体内での光の散乱、吸収特性に応じた関係が成立することから、この関係に基づいて断層画像を再構成することができる。
【０００６】
ところで、この蛍光ＣＴにおいて、断層画像の再構成を行うために蛍光体の濃度分布を演算するときには、励起光強度分布および蛍光強度分布の２つの光強度分布について逆問題計算を行うことが提案されている（例えば、特許文献３、非特許文献１参照）。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開平１１−１７３９７６号公報
【特許文献２】特開平１１−３３７４７６号公報
【特許文献３】米国特許出願公開第２００７／０２８６４６８号明細書
【非特許文献】
【０００８】
【非特許文献１】S. R. Arridge “Optical tomography in medical” Inverse Problems １５（１９９９）Ｒ４１−９３
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、光強度分布を得るための逆問題計算は、順問題計算と比較して計算負荷が大きく装置構成が大規模になりがちであり、計算に長時間を要するという問題がある。近似計算によりこれらの問題を軽減するアプローチも存在するが、ノイズ環境下では近似自体が精度良くできないことも多い。
【００１０】
本発明は上記事実に鑑みてなされたものであり、例えば蛍光ＣＴにおいて装置構成を簡略化するとともに、光強度分布を算出するにあたり計算処理負荷を軽減し、処理時間を短縮化することができる光断層情報生成装置、光強度分布算出方法および光強度分布算出プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１１】
上記目的を達成するために、本発明に係る光断層情報生成装置は、光拡散方程式に基づき、検体から射出される光の強度分布の計算値を求め、当該計算値と実測値との差の大小に応じて、ＡＲＴ（代数的再構成技術）を用いた再計算を行って前記計算値の更新値を得、当該更新値を用いて前記検体の光断層情報を生成する光断層情報生成装置であって、前記検体の光学特性値から、当該光学特性値の変化に対する前記光強度分布の変化の割合を示すヤコビ行列を算出するヤコビ行列算出手段と、前記ヤコビ行列を特異値分解して対角行列を得る特異値分解手段と、前記対角行列の閾値以下の対角成分を０で置換した単位対角行列を取得する単位対角行列取得手段と、前記ＡＲＴを用いた再計算において、ＡＲＴの逐次近似式に現れる前記計算値と実測値との差に対し前記単位対角行列を用いた置換を施して近似を行い、前記更新値を得る逐次近似手段と、を具備する構成を採る。
【００１２】
また、本発明に係る光強度分布算出方法は、光拡散方程式に基づき、検体から射出される光の強度分布の計算値を求め、当該計算値と実測値との差の大小に応じて、ＡＲＴを用いた再計算を行って前記計算値を更新する光強度分布算出方法であって、前記検体の光学特性値から、当該光学特性値の変化に対する前記光強度分布の変化の割合を示すヤコビ行列を算出するステップと、前記ヤコビ行列を特異値分解して対角行列を得るステップと、前記対角行列の閾値以下の対角成分を０で置換した単位対角行列を取得するステップと、前記ＡＲＴを用いた再計算において、ＡＲＴの逐次近似式に現れる前記計算値と実測値との差に対し前記単位対角行列を用いた置換を施して近似を行って、前記計算値を更新するステップと、を具備するようにした。
【発明の効果】
【００１３】
本発明によれば、光強度分布を算出する際の計算処理負荷が軽減されると共に、装置構成の簡略化、計測時間の短縮化を実現することができる。
【図面の簡単な説明】
【００１４】
【図１】実施の形態１に係る光断層情報生成装置の主要な外観構成を示した図である。
【図２】実施の形態１に係る計測部および画像処理部の詳細な構成を示した図である。
【図３】実施の形態１に係る画像処理部の入出力データおよびデータ処理を説明するための機能ブロック図である。
【図４】（Ａ）は検体内での励起光の伝播を示す概略図、（Ｂ）は検体内での蛍光の伝播を示す概略図である。
【図５】検体の光学的特性の一例を示す線図である。
【図６】実施の形態１に係る光強度分布算出方法の手順を示したフローチャートである。
【図７】実施の形態１に係る更新処理部のより詳細な構成を示した機能ブロック図である。
【図８】抗体投与マウスを２次元カメラで腹部方向から撮影して得られた画像である。
【図９】（Ａ）は従来のＡＲＴ法による再構成画像、（Ｂ）は実施の形態１に係る光断層情報生成装置によって得られた再構成画像の例である。
【図１０】ヤコビ行列Ｊのランク数と対角行列Ｄにおける特異値の値との関係を特性曲線として示した図である。
【図１１】（Ａ）〜（Ｆ）は閾値と再構成画像の精度との関係を実際の抗体投与マウスのデータによって検証した結果を示した図である。
【図１２】実施の形態２に係る光断層情報生成装置内の計測部および画像処理部の主要な構成を示すブロック図である。
【図１３】実施の形態２に係る各受光ヘッドの関係等を説明するための図である。
【図１４】実施の形態２に係る光強度分布算出方法の手順を示したフローチャートである。
【発明を実施するための形態】
【００１５】
以下、添付図面を参照しながら本発明の実施の形態を説明する。ここでは、光断層情報生成装置として蛍光トモグラフィ装置を用いる場合を例にとって説明する。
【００１６】
（実施の形態１）
図１は、本発明の実施の形態１に係る光断層情報生成装置１００の主要な外観構成を示した図である。この光断層情報生成装置１００は、計測部１２と画像処理部１４とを含んでいる。
【００１７】
画像処理部１４は、計測部１２から出力される電気信号に基づいた画像処理（信号処理）を行う。画像処理部１４には、表示手段としてＣＲＴ、ＬＣＤなどのモニタ１６が設けられており、このモニタ１６に、計測部１２の計測結果に基づいた画像が表示される。
【００１８】
光断層情報生成装置１００では、小動物（例えばヌードマウス）などの生体を観察対象とする検体１８として、この検体１８から得られる光断層情報に基づいた画像（以下、光断層画像とする）をモニタ１６に表示したり、表示画像の画像データ（光断層情報）を各種の記憶媒体に記憶可能となっている。
【００１９】
計測部１２は、計測ユニット２０を備えている。計測ユニット２０は、リング状の機枠２２を備え、この機枠２２の軸心部が計測位置となっている。計測部１２では、前記検体１８が機枠２２内の計測位置に配置される。
【００２０】
計測ユニット２０には、計測位置へ向けて所定波長の光を照射する光源ヘッド２４と、前記生体から射出される光を検出光Ｌ１として、この検出光Ｌ１を検出する複数の受光ヘッド２６とが、所定の角度間隔で（所定の角度θずつずらされて）機枠２２に取付けられている。本実施の形態に適用した光断層情報生成装置１００では、一例として光源ヘッド２４から３０°ずつ（θ＝３０°）ずらして１１台の受光ヘッド２６を配列している。
【００２１】
この構成により、計測部１２では、光源ヘッド２４から照射した光に対して、検体１８から射出された検出光Ｌ１を、１１台の受光ヘッド２６のそれぞれによって並行して受光可能となっている。
【００２２】
また、計測部１２において、機枠２２が軸心に対して所定の角度ずつ機械的に回転するよう構成されている。これにより、計測中は、計測部１２において、検体１８の周囲の複数点へ向けて光源ヘッド２４から光が照射され、受光ヘッド２６はそれぞれの位置での検出光Ｌ１の受光が可能となっている。ここでは、一例として、機枠２２を角度θ（θ＝３０°）ずつ回転するようにしている。また、計測部１２では、検体１８の周囲の１２点に対して、光を照射し、それぞれの照射点における１１箇所で検出光Ｌ１の検出が可能となっている。なお、光源ヘッド２４の数、受光ヘッド２６の数、これらの配列及び光源ヘッド２４と受光ヘッド２６の移動量（回転量）などは、これに限るものではなく、任意の数、配列及び移動量（回転量）を適用することができる。
【００２３】
計測部１２において、基台２８上に支柱３０が立設されており、この支柱３０に対して機枠２２が回転可能に支持されている。また、支柱３０は、基台２８上に、機枠２２の軸方向（図１の紙面表裏方向）に沿って移動可能に支持されている。即ち、機枠２２は、回転可能となっていると共に、その回転軸の軸方向に沿っても移動可能となっており、機枠２２の回転軸方向に沿った検体１８の任意の部位に対し計測可能となっている。この構成により、３次元の再構成画像を取得することが可能になっている。
【００２４】
なお、このような機枠２２の回転機構及び移動機構は、任意の構成を適用することができる。また、計測部１２では、機枠２２を回転するようにしているが、これに限らず、機枠２２内に配置される検体１８を回転する構成としても良く、また、検体１８と機枠２２のそれぞれが回転されるものであっても良い。
【００２５】
また、計測部１２には、制御ユニット３２が設けられている。
【００２６】
図２は、計測部１２および画像処理部１４の詳細な構成を示した図である。なお、図１と同一の構成要素には同一の符号を付している。
【００２７】
この図に示されるように、制御ユニット３２は、マイクロコンピュータを含んで形成されたコントローラ３４を備えている。また、制御ユニット３２には、例えば、光源ヘッド２４を駆動する発光駆動回路３６、受光ヘッド２６のそれぞれから出力される電気信号を増幅する増幅器（Ａｍｐ）３８、増幅器３８から出力される電気信号をデジタル信号に変換するＡ／Ｄ変換器４０を備え、光源ヘッド２４の発光、各受光ヘッド２６での検出光Ｌ１の受光、受光した検出光Ｌ１の強度を示す測定データの生成が、コントローラ３４によって制御される。
【００２８】
計測部１２は、計測ユニット２０の機枠２２を回転駆動する回転モータ４２、機枠２２を軸方向に移動する移動モータ４４及び、それぞれの駆動回路４６、４８を含み、これらがコントローラ３４に接続された構成とすることができる。
【００２９】
一方、画像処理部１４は、ＣＰＵ５０、ＲＯＭ５２、ＲＡＭ５４、記憶手段とされるＨＤＤ５６、キーボードやマウスなどの入力デバイス５８、モニタ１６等がバス６０に接続された一般的構成のコンピュータが形成されている。これにより、画像処理部１４では、ＲＯＭ５２やＨＤＤ５６に記憶されたプログラム、図示しないリムーバブルメモリなどに記憶されたプログラムに基づいた各種の制御、信号処理、画像処理などが可能となっている。なお、本願では、ＲＯＭ５２、ＲＡＭ５４、ＨＤＤ５６および図示しないリムーバブルメモリを総称して単にメモリ５１と呼ぶこととする。
【００３０】
また、画像処理部１４と計測部１２の制御ユニット３２との間では、制御信号の送受信及びデータの送受信が可能となっている。なお、このような構成は、任意の通信インターフェイスを用いて構成することができる。
【００３１】
図３は、画像処理部１４の入出力データおよび内部のデータ処理を説明するための機能ブロック図である。
【００３２】
この図に示されるように、画像処理部１４は、読取部７０を備え、計測部１２での検体１８の計測を制御しながら、計測部１２（制御ユニット３２）から出力される計測データをメモリ５１に読み込む。
【００３３】
また、画像処理部１４は、演算処理部７２、評価部７４、更新処理部７６、断層情報生成部７８及び断層画像生成部８０を備えている。
【００３４】
演算処理部７２は、蛍光体６２の吸収係数分布を含む予め設定されている検体１８の光学特性値に基づいて光拡散方程式を用いた順問題計算によって蛍光の強度分布を演算する。評価部７４は、演算された蛍光の強度分布と計測データから得られる蛍光の強度分布の差分を評価する。更新処理部７６は、光拡散方程式の逆問題計算を行うことにより、評価部７４の評価結果から得られる差分を減少させるように蛍光の強度分布から蛍光体の濃度分布に基づく吸収係数分布等を設定する。演算処理部７２は、更新処理部７６において蛍光体の濃度分布に基づく吸収係数分布等の更新が行われると、更新された蛍光体の濃度分布に基づく吸収係数分布等に基づいた蛍光強度分布の演算を行う。
【００３５】
このようにして蛍光の強度分布の更新と評価が繰り返され、例えば、演算された蛍光の強度分布と計測データとの差が許容範囲内と評価されると、断層情報生成部７８は、そのときの蛍光体の濃度分布に基づく吸収係数分布等から光断層情報である蛍光の強度分布を生成し、断層画像生成部８０は、この光断層情報に基づいて光断層画像を生成する。
【００３６】
このように、画像処理部１４は、計測部１２から読み込んだ計測データに対して、所定のデータ処理を行った後、この処理結果に基づいた画像処理を行うことにより、計測データに基づいた検体１８の光断層画像を再構成する。
【００３７】
ここで、本実施の形態に係る光断層情報生成装置１００で使用される光強度分布算出方法について、図面および数式を交えながら理論的な説明を行う。
【００３８】
本実施の形態では、図４（Ａ）および図４（Ｂ）に示すように、光断層情報生成装置１００の光源ヘッド２４から発する光を励起光として、この励起光が照射されることにより蛍光を発する蛍光体６２を含む物質ないし薬剤が検体１８に投与される。光断層情報生成装置１００は、検体１８の断層画像として、蛍光体６２の分布を含む画像を再構成し、検体１８内での各種臓器に対する蛍光体６２の分布が視認可能となるようにする。
【００３９】
具体的には、図４（Ａ）に示されるように、検体１８に励起光を照射すると、この励起光が検体１８内で散乱しながら蛍光体６２に達する。これにより、検体１８内の蛍光体６２が発光する。また、図４（Ｂ）に示されるように、蛍光体６２から発せられる蛍光は、検体１８内で散乱を繰り返しながら、検体１８から射出される。このとき、蛍光と共に励起光も射出されるが、受光ヘッド２６は、その前段に図示しない励起光カットフィルタを有しているため、検体１８から射出される蛍光のみを検出光として、この検出光（以下、蛍光とする）の強度を検出する。光断層情報生成装置１００は、この検出光の光強度分布から検体１８内での蛍光体６２の分布（濃度分布）を得る。
【００４０】
ここで、励起光などの光を検体１８に照射した場合、照射位置近傍の領域では、光に対する屈折率が方向によって異なるなどの異方散乱領域となっているが、検体１８内で所定距離以上に進行すると等方散乱領域となる。
【００４１】
検体１８内を散乱する光は、エネルギーを輸送する粒子と見なされることから、光強度の分布は、光強度のエネルギー保存式である光輸送方程式を用いて表すことができる。しかし、この光輸送方程式を解くことは現状では困難とされている。
【００４２】
一方、検体１８では、一般的に異方散乱領域が数ｍｍ程度であるために、数ｃｍ以上の大きさの検体１８では、その体内を実質的に等方散乱領域と見なすことができる。すなわち、検体１８での光の散乱は等方散乱として近似することができる。
【００４３】
ここから、光拡散方程式を用いることにより、光強度の分布が得られる。この光拡散方程式は、次式（１）で表される。ここで、Φ（ｒ、ｔ）は検体１８内の光密度、Ｄ（ｒ）は拡散係数分布、μａ（ｒ）は吸収係数分布、ｑ（ｒ、ｔ）は光源の光密度を表し、ｒは検体１８内の座標位置、ｔは時間を表す。
【００４４】
【数１】
【００４５】
ここで時間的に連続した光であれば、光拡散方程式は次式（２）で示すことができる。
【００４６】
【数２】
【００４７】
光学特性値である拡散係数分布Ｄ（ｒ）および吸収係数分布μａ（ｒ）が既知であるときに、この光拡散方程式（２）を用いて検体１８から射出される光強度分布を求める場合、順問題として計算することができる。しかし、光強度分布が既知であり、ここから、光拡散方程式を用いて検体１８の光学特性値である拡散係数分布Ｄ（ｒ）および吸収係数分布μａ（ｒ）を求めることは、逆問題計算となる。
【００４８】
ここで、検体１８での拡散係数分布Ｄ（ｒ）及び吸収係数分布μａ（ｒ）は、光の波長によって異なるため、励起光の波長に対する拡散係数分布をＤｓ（ｒ）、吸収係数分布をμａｓ（ｒ）とし、光源の光密度をｑｓ（ｒ）とすると、励起光に対する光拡散方程式は次式（３）で表され、一方、蛍光の波長に対する拡散係数分布をＤｍ（ｒ）、吸収係数分布をμａｍ（ｒ）とし、蛍光の光源をｑｍ（ｒ）とすると、蛍光に対する光拡散方程式は次式（４）で表される。また、蛍光の光源ｑｍ（ｒ）は、検体１８内の光密度Φｓ（ｒ）、蛍光体６２の量子効率γ、モル吸光係数εおよび濃度分布Ｎ（ｒ）を用いて、ｑｍ＝γ・ε・Ｎ（ｒ）・Φｓ（ｒ）と表すことができる。したがって、式（４）は次式（５）で置き換えられる。
【００４９】
【数３】
【００５０】
一方、図５に破線で示されるように、検体１８内のヘモグロビンは約７００ｎｍ以下の波長の光に対して吸収が強く、また、図５に二点鎖線で示されるように、検体１８内の水分は波長が約１μｍ以上の光に対して吸収が強い。換言すれば、７００ｎｍ〜１μｍの波長域は検体１８における吸収が弱いということになり、この帯域は所謂光学的窓と呼ばれている。したがって、検体１８は、光学的窓に相当したものとなっている。この波長域では、検体１８の吸収係数分布μａは、０．００２ｍｍ−１〜０．１ｍｍ−１の範囲となる。
【００５１】
また、検体１８内での光の散乱係数（散乱の強さ、図５に実線で示す）は、波長が長くなると小さくなるが、その変化は緩やかであり、光学的窓となっている７００ｎｍ〜１μｍの波長域では、散乱の強さを波長に依存せず略一定と見なすことができる。
【００５２】
そのため、本実施の形態に係る光断層情報生成装置１００は、光源ヘッド２４から発する励起光として、生体（検体１８）での光学的窓に該当する７００ｎｍ〜１μｍの波長の赤外線（近赤外線）を用いる。これにより、検体１８での光学的特性である吸収係数μａ、散乱係数（拡散係数Ｄ）を略一定の値（既知の値）とすることができ、式（３）及び式（５）において、Ｄｓ（ｒ）＝Ｄｍ（ｒ）＝Ｄ（ｒ）、μａｓ（ｒ）＝μａ（ｒ）＋ε・Ｎ（ｒ）、μａｍ（ｒ）＝μａ（ｒ）と簡略化することができる。ここで、εはモル吸光係数、Ｎ（ｒ）は蛍光体６２の濃度分布を表し、ε・Ｎ（ｒ）は蛍光体６２による励起光の吸収度合いを表す。したがって、式（３）は次式（６）で置き換えられ、式（５）は次式（７）で置き換えられる。
【００５３】
【数４】
【００５４】
また、本実施の形態で断層画像を観察する検体１８には、この近赤外線を励起光として発光する蛍光体６２を含む物質ないし薬剤が投与されている。このとき、上記式（７）に示されるように、蛍光体６２が光源となるときの蛍光の強度は、励起光の強度Φｓ（ｒ）に基づくもの、すなわちΦｓ（ｒ）の関数である。これは、拡散係数分布や吸収係数分布を予め設定すれば既知であり、光源の光強度ｑｓ（ｒ）も既知であることから、有限要素法などの数値解析手法により検体１８内の光強度Φｓ（ｒ）を順問題として求めることができる。
【００５５】
したがって、計測部１２で蛍光強度を計測し、この計測データに基づき蛍光強度分布Φｍ（ｒ）を求め、一つ（１系統）の逆問題計算によって検体１８内での蛍光体６２の濃度分布Ｎ（ｒ）を求めることができる。計測部１２で励起光強度を計測する必要はない。
【００５６】
得られた蛍光体６２の位置ｒにおける濃度分布Ｎ（ｒ）を検体１８の位置ｒにおける断層画像に合成することにより、検体１８内における蛍光体６２の濃度分布を光断層画像上で視認することができる。
【００５７】
次いで、上記理論に基づいて、本実施の形態に係る光断層情報生成装置１００が具体的にどのような演算処理を行っているかを説明する。
【００５８】
光断層情報生成装置１００では、計測部１２の計測ユニット２０に検体１８が配置されると、光源ヘッド２４から予め設定された波長の近赤外光を励起光として検体１８に照射する。この励起光は、検体１８内を拡散しながら伝播（透過）する。
【００５９】
ここで、検体１８に蛍光体６２が投与されていると、この蛍光体６２に励起光が照射され、これにより、蛍光体６２が発光する。この蛍光は、検体１８内を拡散しながら伝播して、検体１８から周囲へ射出される。
【００６０】
計測ユニット２０には、検体１８を囲うように受光ヘッド２６が所定の角度間隔で配列されており、計測部１２では、検体１８から射出される蛍光を検出光として、各受光ヘッド２６で受光する。
【００６１】
また、計測部１２では、機枠２２を回転することにより検体１８への励起光の照射位置及び検体１８から射出される蛍光の検出位置を相対的に変えて、励起光の照射、検出光の受光を繰り返す。これにより、検体１８の周囲に沿って照射した励起光に応じた蛍光の強度の測定データが得られる。
【００６２】
光断層情報生成装置１００の画像処理部１４では、この計測データに基づいて蛍光体６２の濃度分布Ｎ（ｒ）を演算する。
【００６３】
図６は、画像処理部１４が行う蛍光体６２の濃度分布計算の処理手順を示したフローチャートである。前述の通り、蛍光体の濃度分布を算出することは蛍光強度分布を算出することに相当するので、この図は、本実施の形態に係る光強度分布算出方法の手順を示していると言える。
【００６４】
なお、画像処理部１４では、検体１８の光学的特性に基づいて予め設定されている波長の近赤外線を励起光としており、ここから、吸収係数分布μａ（ｒ）および拡散係数分布Ｄ（ｒ）の値が予め設定されて記憶されている。なお、この吸収係数分布μａ（ｒ）および拡散係数分布Ｄ（ｒ）は、検体１８の個体差に合わせて適宜入力される値であっても良い。
【００６５】
このフローチャートでは、ステップ（以下ＳＴと略す）１０００で、計測部１２の複数の受光ヘッド２６によって得られた実測値、すなわち、検体１８から射出される蛍光強度分布Φｍ（ｒ）measがメモリ５１に読み込まれる。
【００６６】
一方、ＳＴ１０２０では、予め設定されている吸収係数分布μａ（ｒ）及び拡散係数分布Ｄ（ｒ）が用いられ、検体１８内に蛍光体が存在しなかった場合、すなわち蛍光体による励起光の吸収が生じていない場合の励起光強度分布Φｔ（ｒ）calcが下記光拡散方程式（８）に従って算出される。
【００６７】
【数５】
【００６８】
ＳＴ１０４０では、蛍光体６２の濃度分布Ｎ（ｒ）calcの初期値が設定される。ＳＴ１０６０では、ＳＴ１０４０において設定された蛍光体の濃度分布Ｎ（ｒ）calcと、予め設定されている吸収係数分布μａ（ｒ）及び拡散係数分布Ｄ（ｒ）とを用いて、既に説明した式（６）を適宜修正した下記光拡散方程式（９）に従って、励起光強度分布Φｓ（ｒ）calcが計算される。
【００６９】
【数６】
【００７０】
ＳＴ１０７０では、下記式（１０）に示すように蛍光体による励起光の吸収が生じていない場合の励起光強度分布Φｔ（ｒ）calcから蛍光体による励起光の吸収が生じている場合の励起光強度分布Φｓ（ｒ）calcを減じ、さらに励起光を蛍光に換算する係数γを乗ずることにより、検体１８から射出される蛍光強度の予測値である計算蛍光強度分布Φｍ（ｒ）calcが計算される。
【００７１】
【数７】
【００７２】
なお、ＳＴ１０２０における励起光強度分布Φｔ（ｒ）calcの算出やＳＴ１０６０における励起光強度分布Φｍ（ｒ）calcの算出は、数学的モデルである光拡散方程式を有限要素法などの数値解析手法を用いた公知の順問題計算として比較的容易に演算することができる。
【００７３】
ＳＴ１０８０では、計測データに基づいた蛍光強度分布Φｍ（ｒ）measと、演算結果に基づいた蛍光強度分布Φｍ（ｒ）calcとを比較し、ＳＴ１１００では、蛍光強度分布Φｍ（ｒ）measと蛍光強度分布Φｍ（ｒ）calcとの間の差分が許容範囲内か否か、具体的には予め設定した所定値以内か否かを判断する。
【００７４】
ＳＴ１０８０において、蛍光強度分布Φｍ（ｒ）measと蛍光強度分布Φｍ（ｒ）calcとの間の差分が所定値以内でないと判断されるときには、ＳＴ１１００で否定判定されてＳＴ１１２０へ移行する。
【００７５】
ＳＴ１１２０では、ヤコビ行列（Jacobian matrix）を用いた公知の手法で検体１８の光学特性値の変化に対する光強度分布の変化を演算する。
【００７６】
ＳＴ１１４０では、以下の光拡散方程式（１１）で表される逆問題に対して後述の最適化手法を適用することにより、蛍光体６２の濃度分布の推定値Ｎ（ｒ）estを求める。
【００７７】
【数８】
【００７８】
より詳細には、ＳＴ１１４０ではまず、蛍光強度分布Φｍ（ｒ）measと蛍光強度分布Φｍ（ｒ）calcとの間の誤差（例えば、二乗誤差ｙ）を評価する。すなわち、二乗誤差ｙは次式（１２）から得られ、この二乗誤差ｙを評価する。
【００７９】
【数９】
【００８０】
次に、このＳＴ１１４０では、この二乗誤差ｙを最小とする蛍光体６２での蛍光の吸収εＮ、すなわち、蛍光体６２の濃度分布Ｎ（ｒ）estを後述の最適化手法により推定する。
【００８１】
このようにして濃度分布の推定値Ｎ（ｒ）estが求まるが、これは、上記式（１１）のＮ（ｒ）calcに対する推定値を求めたに過ぎない。そこで、ＳＴ１１６０では、この濃度分布の推定値Ｎ（ｒ）estによりＮ（ｒ）calcを更新し、ＳＴ１０６０へ戻る。
【００８２】
画像処理部１４は、このＳＴ１０６０からＳＴ１１６０を繰返し行う。これにより、蛍光強度分布Φｍ（ｒ）measと蛍光強度分布Φｍ（ｒ）calcとの差が所定値以内と判断されると、ＳＴ１１００において肯定（ＹＥＳ）判定してＳＴ１１８０へ移行し、ＳＴ１１８０ではこのとき最終的に設定されている濃度分布Ｎ（ｒ）calcを計測データから得られた濃度分布Ｎ（ｒ）として出力する。
【００８３】
なお、ここでは、計測データに基づいた蛍光強度分布Φｍ（ｒ）measを取得した後、蛍光体６２の濃度分布Ｎ（ｒ）calcの初期値を設定するようにしているが、処理手順はこれに限らず、蛍光体６２の濃度分布Ｎ（ｒ）calcの初期値を設定した後に、蛍光強度分布Φｍ（ｒ）measを取得するようにしても良い。他の手順についても、当該フローの目的を逸脱しない範囲で適宜ステップの順番を変更することができる。
【００８４】
本発明では、ＳＴ１１４０で使用する逆問題に対する最適化手法としてＡＲＴ（Algebraic Reconstruction Technique：代数的再構成技術）法を採用する。しかし、ＡＲＴ法は、収束は早いがノイズに対する耐性（ノイズ耐性）に乏しいという特徴がある。そこで、本願発明者は、以下に示す改良型ＡＲＴ法を見い出した。なお、ノイズ耐性とは、ノイズが多い環境下においても所定時間内に演算を収束させることができる性質、またはノイズが多い環境下においても所定レベル以上の演算精度を維持できる性質のことである。
【００８５】
具体的に数式を交えて説明すると、上記逆問題を解くためには、式（１２）における二乗誤差ｙを最小にする必要がある。よって、この計算を行うために、まず式（１２）をεＮで微分して簡略化（Born近似）すると次式（１３）が得られる。
【００８６】
【数１０】
【００８７】
ここで、ΔｘはεＮの微小変化、具体的には、本実施の形態ではＡＲＴ法という繰り返し演算（反復法）を用いているためεＮの更新値を表している。また、ヤコビ行列Ｊは、光学特性値の変化に対する光強度分布の変化の割合を示している。なお、光強度分布は、光源の座標位置と検出地点の座標位置との関数として表される。
【００８８】
式（１３）においてΔｘが求めたい解であるため、式（１３）をΔｘについて書き直すと次式（１４）となる。
【００８９】
【数１１】
【００９０】
原理的には当該式（１４）を用いて計算を進めることができるはずであるが、実際には式（１４）をそのまま使用して計算機上で計算すると次のような問題が生じる。すなわち、ヤコビ行列Ｊ中に著しく小さい成分が存在する場合に逆行列Ｊ−１においては逆数となって現れるため、この成分がノイズのような小さな誤差を含んでいた場合にその誤差が拡大されて計算精度に著しく悪影響を与えてしまう。ＡＲＴ法がノイズ耐性に乏しい理由はこのためである。
【００９１】
そこで本願では、まず行列分解の一手法である特異値分解法を用いてヤコビ行列Ｊを分解して次式（１５）を得る。
【００９２】
【数１２】
【００９３】
ここで、Ｄはヤコビ行列の特異値をその対角成分とするｍ×ｎ行の対角行列、Ｕはｍ×ｍ行の正規直交行列、ＶＴはｎ×ｎ行の転置行列である。ＶＴはＶの転置行列である。ＵおよびＶには、ＵＴＵ＝ＶＴＶ＝Ｉ（単位行列）の関係がある。
【００９４】
式（１４）に代入できるように、式（１５）をＪ−１を求める形に書き換えると次式（１６）となる。
【００９５】
【数１３】
【００９６】
さらに、対角行列Ｄとは別個に新たな単位対角行列Ｈを導入し、Ｄにおいて特異値が閾値αth以下となっている対角成分を特定し、Ｈにおける同位置の成分を０で置換する（特異値が閾値αthよりも大きい成分は１を維持する）。このＨを用いて、式（１４）に対し次式（１７）で表されるΔｙの置換を行う。この置換において、０に置換した成分は閾値αth以下のもののみであるので、Δｙ全体の大きさ（絶対値）に大きな変化はないことが見込まれる。
【００９７】
【数１４】
【００９８】
そして、上記式（１４）、式（１６）および式（１７）から次式（１８）を得る。
【００９９】
【数１５】
【０１００】
式（１８）右辺のＤ−１Ｈの部分は、小さい値の成分については０との積となっており、また、他の対角成分についてはＤ−１の値をそのまま維持したものとなっている。すなわち、本願に係る改良型ＡＲＴ法により、ノイズの影響を減少させつつ蛍光体濃度分布の計算を行うことができる。
【０１０１】
本発明ではＡＲＴ法を用いているので、具体的には、式（１８）の代わりに次式（１９）のＡＲＴの逐次近似式を計算することになる。よって、式（１７）を用いたΔｙの置換は、式（１９）に対して行われる。ここで、λは緩和係数と呼ばれ、解の収束性を制御するため計算上使用されるパラメータである。
【０１０２】
【数１６】
【０１０３】
図７は、上記計算を実現する本実施の形態に係る更新処理部７６のより詳細な構成を示した機能ブロック図である。更新処理部７６は、ヤコビ行列算出部１５１、特異値分解部１５２、単位対角行列取得部１５３および逐次近似部１５４を備える。
【０１０４】
ヤコビ行列算出部１５１は、図３に示した評価部７４から出力される評価結果信号Ｓ１１が否定（ＮＯ）を示す場合、図３に示した演算処理部７２から出力されるデータ信号Ｓ１２からヤコビ行列Ｊを求め、特異値分解部１５２へこれを出力する。ここで、データ信号Ｓ１２は、検体１８の光学特性値、図１および図２に示した光源ヘッド２４の座標位置および受光ヘッド２６の位置情報を含み、ヤコビ行列算出部１５１は、このデータ信号Ｓ１２に基づいて上記式（１３）におけるヤコビ行列Ｊ、より具体的には上記式（１９）におけるヤコビ行列Ｊを求める。
【０１０５】
特異値分解部１５２は、ヤコビ行列算出部１５１で求められたヤコビ行列Ｊを上記式（１５）に従って特異値分解し、得られた対角行列Ｄを単位対角行列取得部１５３へ出力し、同様に得られた正規直交行列Ｕおよび転置行列ＶＴを逐次近似部１５４へ出力する。
【０１０６】
単位対角行列取得部１５３は、特異値分解部１５２から出力された対角行列Ｄに基づいて、既に説明した本願の手法によって単位対角行列Ｈを求め、逐次近似部１５４へ出力する。
【０１０７】
逐次近似部１５４は、特異値分解部１５２から出力された正規直交行列Ｕおよび転置行列ＶＴ、並びに単位対角行列取得部１５３から出力された単位対角行列Ｈを用いて上記式（１９）の逐次近似式を計算し、εＮの更新値であるΔｘを求めて図３に示した演算処理部７２へ出力する（信号Ｓ１３）。
【０１０８】
次いで、本実施の形態に係る光断層情報生成装置１００が奏する効果を実例を紹介しながら説明する。
【０１０９】
図８は、実際に腫瘍を有するヌードマウスに対し抗体を投与し、このマウスを２次元カメラで腹部方向から撮影して得られた画像である。この画像から、マウスの右腹部に腫瘍が存在することが確認できる。なお、本願における動物実験は、出願人である富士フイルム株式会社内の動物倫理規則および安全性評価センター動物実験規程に則って行われたものである。
【０１１０】
図９（Ｂ）は、光断層情報生成装置１００を用いて、この抗体投与マウスから得られた再構成画像の例である。図９（Ａ）に示した従来のＡＲＴ法による再構成画像と比較すると、従来のＡＲＴ法ではマウス体内に腫瘍の存在を認めることはできないが、本実施の形態に係る光断層情報生成装置１００で得られた再構成画像では、マウス右腹部に濃度分布の偏り、すなわち腫瘍の影が認められることがわかる。
【０１１１】
以上説明したように、本発明の実施の形態１に係る光断層情報生成装置１００によれば、励起光が照射された検体（生体）から射出される蛍光の強度を検出し、蛍光の強度（蛍光強度）の測定データを取得する。また、検体の光学特性値である散乱係数及び吸収係数の各々の分布を予め設定すると共に、蛍光体の濃度分布に基づく蛍光体の吸収係数分布を仮設定しておき、数学的モデルから蛍光の強度を取得する。この数学的モデルから取得した蛍光の強度と、測定データとして得られる蛍光の強度を比較し、この差分に対する評価を行う。このとき、差分が所定レベルより大きければ、逆問題計算を行うことにより、この差分を減少させるように吸収係数分布を推定し、推定した吸収係数分布に基づいた蛍光の強度と測定データを比較する処理を反復することにより、蛍光体の濃度分布に基づく吸収係数分布を決定する。このようにして得られる蛍光体の濃度分布に基づく吸収係数分布から検体内の蛍光体の濃度分布を求めて、検体の光断層画像を形成するための光断層情報を生成する。
【０１１２】
これにより、光断層情報を取得するための光拡散方程式を用いた逆問題計算は、蛍光の強度に対して行うのみでよいために、断層情報を生成するための処理負荷が軽減され、処理時間の短縮を図ることができる。また、断層画像を形成するための光の強度測定は、検体から射出される蛍光に対して行うのみでよいため、装置の簡略化、計測時間の短縮を図ることができる。例えば、受光ヘッド２６は蛍光のみを計測できる機能があれば充分で、励起光を計測できる必要はない。
【０１１３】
これに対して、従来構成では、散乱係数Ｄ（ｒ）及び吸収係数分布μａ（ｒ）を既知の値として設定しないため、測定データとして蛍光強度分布に加え励起光強度分布が必要となり、蛍光と励起光の双方に対応した受光ヘッドが必要となる。また、このときには、励起光に対する光拡散方程式、および蛍光に対する光拡散方程式のそれぞれに対して逆問題計算を行い、蛍光体の濃度分布を得る必要がある。従って、各々の逆問題計算を処理する装置構成が必要となり、装置が大規模化する。
【０１１４】
また、以上の構成において、前記光源を前記検体の周囲に沿って相対移動することにより前記励起光の前記検体への照射位置を変更しながら、それぞれの照射位置で前記蛍光の強度を検出する。このとき、前記蛍光を検出する複数の検出手段を前記検体の周囲に予め設定された間隔で設け、それぞれの検出手段を前記光源の前記検体に対する相対移動と対になって移動して前記蛍光の強度を検出する構成を取り得る。
【０１１５】
これにより、光断層画像生成装置における蛍光の強度を検出するための構成が簡略化できる。
【０１１６】
また、以上の構成において、本実施の形態１に係る光断層情報生成装置１００は、光強度分布を算出する際の逆問題においてＡＲＴ（代数的再構成技術）を使用し、このＡＲＴの逐次近似式中に現れるヤコビ行列Ｊに対して特異値分解を行って対角行列Ｄを得、このＤにおいて閾値αth以下となっている特異値を０で置換した単位対角行列ＨをＤに乗ずることにより、逆問題の計算負荷を減ずると共に、処理時間の短縮化を実現している。
【０１１７】
また、光断層情報生成装置１００は、ＡＲＴで逆問題を解く計算工程において、光学特性値の変化に対する光強度分布の変化の割合を複数の独立な要素に分解し、そのうちノイズの影響が支配的な要素を特定して除去してから逆問題を解くようにしている。これにより、ノイズが多く含まれている環境下においても、逆問題計算を迅速に収束させることができ、さらに再構成画像の精度を向上させることができる。なお、付言するならば、上記計算において、ノイズが支配的でない要素については値を維持したままとしているので、ノイズの影響を除去する過程において正味の信号成分まで除去してしまうことを防止しており、この点からも再構成画像の精度向上を実現できている。
【０１１８】
また、以上の構成において、このような、単位対角行列Ｈの導入は、ＡＲＴのような繰り返し近似計算法において計算を単純化簡易化することができ、実装置への実装化やプログラミングのし易さに繋がっている。
【０１１９】
また、以上の構成において、対角行列Ｄから単位対角行列Ｈを生成する際に、Ｄにおいて特異値が閾値αth以下となっている対角成分を０で置換することを説明した。この閾値αthとして妥当な値を設定することにより、本実施の形態に係る光強度分布の算出精度、ひいては光断層情報生成装置１００で得られる再構成画像の精度をより向上させることができる。この閾値αthの具体的な設定方法について以下説明する。
【０１２０】
図１０は、上記ヤコビ行列Ｊのランク数（横軸）と対角行列Ｄにおける特異値の値（縦軸）との関係を特性曲線として示した図である。この図からわかるように、Ｄにおける特異値は、ヤコビ行列のランク数が増加するのに伴いステップ状に減少する性質があることがわかった。
【０１２１】
そこで、本願発明者は、対角行列Ｄから単位対角行列Ｈを生成する際に使用する閾値αthとして、このステップ状に変化する箇所、即ち、その他の箇所と比較して急激に特異値が減少する箇所（図では破線で示す）の特異値の値を閾値αthとして用いることを思い付いた。閾値αth以下の特異値を０で置換するということは、ノイズの影響を除去する過程において正味の信号成分までも過大に除去してしまうおそれがある。すなわち、ノイズの影響除去と正味の信号成分の維持はトレードオフの関係にある。しかし、特異値が急激に変化する箇所を閾値αthに設定すれば、信号成分を過大に除去せずにノイズの影響を除去できる可能性が高まることに本願発明者は気付いたのである。具体的には、特異値が１０−３、１０−４、１０−５、１０−６の値を示す箇所の近傍において特異値がステップ状に変化している。そこで、本実施の形態に係る光断層情報生成装置では、１０−３、１０−４、１０−５、１０−６を閾値αthとして用いることにする。また、特異値が１０−８を示す箇所は厳密にはステップ状になっていないが、他の閾値との比較のためにこの値も閾値αthとして追加して用いることにする。
【０１２２】
なお、図１０に示した特性曲線は、光断層情報生成装置のシステム精度や計算環境によって変化するものである。従って、上記の閾値αthの具体的な値は実施の形態の一例にすぎない。
【０１２３】
図１１（Ａ）〜（Ｆ）は、実際の抗体投与マウスを用いて再構成画像の精度を検証した結果である。閾値αthとして上記の値を用いた場合に実際の再構成画像の精度がどのように変化するかをこの図から判断することができる。
【０１２４】
図１１（Ａ）は従来のＡＲＴ法の場合、図１１（Ｂ）は閾値αthとして１０−３を用いた場合、図１１（Ｃ）は閾値αthとして１０−４を用いた場合、図１１（Ｄ）は閾値αthとして１０−５を用いた場合、図１１（Ｅ）は閾値αthとして１０−６を用いた場合、図１１（Ｆ）は閾値αthとして１０−８を用いた場合である。この図からわかるように、閾値αthとして１０−４または１０−５を用いた場合に、腫瘍と思われる領域を他の領域と区別し得る精度の再構成画像を得ることができる。従って、本実施の形態に係る光断層情報生成装置が使用する閾値αthとしては１０−４または１０−５を用いる。また、閾値αthとして１０−５を用いた場合、腫瘍と思われる領域が他の領域と比べて識別可能な精度になるだけでなく、腫瘍の大きさについても一定の精度で再現できていることがわかる。よって、閾値αthとして１０−５を用いることが更に望ましい。
【０１２５】
このように、本実施の形態に係る光断層情報生成装置では、対角行列Ｄから単位対角行列Ｈを生成する際に使用する閾値αthの候補として、ヤコビ行列Ｊのランク数と対角行列Ｄにおける特異値の値との関係を示す特性曲線がステップ状に変化する複数の箇所の値を使用する。そして、複数の閾値αthの候補の中から最適な閾値αthを適宜選択する。これにより、本実施の形態に係る光強度分布の算出精度、ひいては光断層情報生成装置１００で得られる再構成画像の精度をより向上させることができる。
【０１２６】
なお、本実施の形態では、閾値αthとして、ヤコビ行列Ｊのランク数と対角行列Ｄにおける特異値の値との関係を示す特性曲線がステップ状に変化する箇所の値を使用する場合を例にとって説明したが、例えば、よりノイズの影響を除去したい場合には、上記設定方法で求めた閾値αthに対しオフセットを加えた値を閾値として用いても良い。
【０１２７】
また、本実施の形態では、光学特性値の変化に対する光強度分布の変化の割合を示すヤコビ行列の各成分のうち、ノイズの影響が大きく現れ得る成分を０に置換する場合、すなわち当該成分を完全に削除する場合を例にとって説明したが、この成分には正味の信号成分も含まれているはずであるから、完全に削除するのではなく、逆数をとったときに所定値以下の大きさに収まるような定数で置換するような構成としても良い。これにより、計算負荷等を軽減しつつ再構成画像の精度を向上させることができる。
【０１２８】
（実施の形態２）
実施の形態１では、検体から発せられる蛍光のみを実測して１系統の逆問題を解く態様を例にとって説明したが、本実施の形態２では、検体から発せられる蛍光だけでなく励起光も計測し、２系統の逆問題を解く態様について説明する。
【０１２９】
図１２は、本発明の実施の形態２に係る光断層情報生成装置２００内の計測部２１０および画像処理部１４ａの主要な構成を示すブロック図である。なお、実施の形態１で示した光断層情報生成装置１００と同一の構成要素には同一の符号を付し、その説明を省略する。また、基本的動作は同一であるが、詳細な点で違いがある構成要素には、同一の番号にアルファベットの小文字を付した符号を付して適宜説明を加える（以降の図の説明においても同様の表記方法を用いる）。例えば、本実施の形態２に係る画像処理部１４ａは、実施の形態１で示した画像処理部１４と基本的構成は同一であるが詳細な動作の点においては違いがあるため、同一の番号１４にアルファベットａを付している。
【０１３０】
光断層情報生成装置２００は、実施の形態１で示した光断層情報生成装置１００の構成に更に新規な構成である受光ヘッド２１１および増幅器（Ａｍｐ）２１２を有している。受光ヘッド２４は、実施の形態１で説明した通り蛍光Ｌ１用の受光ヘッドであるが、受光ヘッド２１１は、励起光Ｌ２用の受光ヘッドである。また、増幅器２１２は、受光ヘッド２１１用の増幅器であり、受光ヘッド２１１から出力される電気信号を増幅し、Ａ／Ｄ変換器４０へ出力する。
【０１３１】
ここで、受光ヘッド２１１〜Ａ／Ｄ変換器４０の経路は新たに設けられた経路であり、光断層情報生成装置１００における受光ヘッド２６〜Ａ／Ｄ変換器４０の経路は光断層情報生成装置２００においてもそのまま維持されている。即ち、光断層情報生成装置２００における受光ヘッド２６および増幅器３８の数は光断層情報生成装置１００のそれらと同数である。また、励起光用の受光ヘッド２１１および増幅器２１２の数も、蛍光用の受光ヘッド２６および増幅器３８のそれらと同数である。
【０１３２】
図１３は、受光ヘッド２６および受光ヘッド２１１の関係、およびこれらの構成の機能を説明するための図である。
【０１３３】
この図からわかるように、光断層情報生成装置２００の計測ユニット２０ａの基本構成は、実施の形態１に係る光断層情報生成装置１００の計測ユニット２０と類似している。しかし、検体１８から射出された光が受光ヘッド２６に到達するまでの経路においてダイクロイックミラー２１３が設けられ、検体１８から射出される光のうち所定の短波長成分のみを反射し、受光ヘッド２１１へ入射するように構成されている。検体１８から射出される光は、蛍光Ｌ１の成分のみならず励起光Ｌ２の成分も含む。そこで、ダイクロイックミラー２１３は、蛍光Ｌ１に比べ短波長である励起光Ｌ２の成分のみを反射するように構成され、受光ヘッド２１１には励起光成分のみが入射される。また、実施の形態１で説明した通り、受光ヘッド２６の前段には図示しない励起光カットフィルタが設けられているため、受光ヘッド２６には蛍光Ｌ１のみが入射される。
【０１３４】
図１４は、画像処理部１４ａが行う蛍光体６２の濃度分布計算の処理手順、すなわち本実施の形態に係る光強度分布算出方法の手順を示したフローチャートである。
【０１３５】
前述の通り、実施の形態１で測定されるのは蛍光のみであったが、本実施の形態２では蛍光に加え励起光も測定する。そのため、ＳＴ２０００では、検体１８から射出される蛍光強度分布の測定データΦｍ（ｒ）measがメモリ５１に読み込まれると共に、ＳＴ２０２０において、検体１８から射出される励起光強度分布の測定データΦｘ（ｒ）measもメモリ５１に読み込まれる。なお、ＳＴ２０００は、実施の形態１の図６で示したＳＴ１０００と同一の処理である。
【０１３６】
ＳＴ２０４０では、ＳＴ２０００で読み込まれた蛍光強度分布Φｍ（ｒ）measとＳＴ２０２０で読み込まれた励起光強度分布Φｘ（ｒ）measとから、次式（２０）に基づいて、蛍光体が存在しなかった場合の励起光強度分布Φｔ（ｒ）measが計算される。
【０１３７】
【数１７】
【０１３８】
ＳＴ２０６０では、ＳＴ２０８０以降のループ計算の準備として、蛍光体が存在する場合の励起光波長における吸収係数分布μａｘ（ｒ）の初期値、および蛍光体が存在しなかった場合の励起光波長における吸収係数分布μａｔ（ｒ）の初期値が設定される。なお、本実施の形態でも、実施の形態１と同様に拡散係数分布Ｄ（ｒ）は略一定の既知の値とする。
【０１３９】
ＳＴ２０８０では、下記光拡散方程式（２１）で示される順問題が解かれることにより、蛍光体が存在する場合に検体１８から射出される励起光強度分布Φｘ（ｒ）calcが求められ、同様に、下記光拡散方程式（２２）で示される順問題が解かれることより、蛍光体が存在しなかった場合に検体１８から射出される励起光強度分布Φｔ（ｒ）calcが求められる。
【０１４０】
【数１８】
【０１４１】
ＳＴ２１００では、蛍光体が存在する場合の励起光強度分布の計算結果（予想値）であるΦｘ（ｒ）calcと実測値であるΦｘ（ｒ）measとが比較され、また、蛍光体が存在しない場合の励起光強度分布の計算結果（予測値）であるΦｔ（ｒ）calcと実測値であるΦｔ（ｒ）measとが比較される。具体的には、比較結果として、各予測値と各実測値との差分が求められる。
【０１４２】
ＳＴ２１２０〜ＳＴ２１８０の計算は、実施の形態１と同様のループ計算の一部である。すなわち、ＳＴ２１２０において、ＳＴ２１００で求められた各予測値と各実測値との間の差分が許容範囲内か否か、具体的には予め設定した所定値以内か否かが判断され、差分が所定値以内でないと判断されるときには、ＳＴ２１４０〜ＳＴ２１８０の計算を経た後、ＳＴ２０８０に戻るループ計算である。ここで、実施の形態１と大きく異なる点は、ＳＴ２１６０において行われる逆問題計算である。実施の形態１では、検体から発せられる蛍光のみを実測して１系統の逆問題を解いていたが、本実施の形態２では、検体から発せられる蛍光だけでなく励起光も計測し、２系統の逆問題を解くからである。
【０１４３】
なお、ＳＴ２１００で求まる差の値は２個あるため、厳密には本実施の形態２に係るＳＴ２１２０〜ＳＴ２１８０の計算は、実施の形態１に係るＳＴ１１００〜１１６０の計算とは異なる。しかし差が２個あったとしても、差を所定値以内に収めようとする実施の形態１における繰り返し演算の趣旨は、本実施の形態でも同様であるため、具体的な繰り返し計算方法をどのように規定するかは当業者が適宜決定すれば良いことである。そこで、ＳＴ２１２０では説明を簡単にするため、あたかも差が１個であるが如く説明を行っている。また、ＳＴ２１６０およびＳＴ２１８０では、取り扱われるパラメータが実施の形態１の濃度分布と異なり吸収係数分布である。これは、ループ計算の初期値としてＳＴ２０６０において設定し、ＳＴ２１８０において更新するパラメータが、実施の形態１と異なり吸収係数分布だからである。
【０１４４】
ＳＴ２１２０において各予測値と各実測値との間の差分が所定値以内と判断されるとＳＴ２２００へ移行し、ＳＴ２２００ではこのとき最終的に設定されている吸収係数分布μａｘ（ｒ）および吸収係数分布μａｔ（ｒ）から次式（２３）に従って濃度分布Ｎ（ｒ）を算出し、本フローチャートは終了する。なお、吸収係数分布と濃度分布との間には次式（２４）および（２５）の関係があり、式（２３）はこれらの関係から求められたものである。
【０１４５】
【数１９】
【０１４６】
以上説明したように、本発明の実施の形態２に係る光断層情報生成装置２００によれば、検体から発せられる蛍光だけでなく励起光も計測し、２系統の逆問題を解くために、装置構成が実施の形態１よりは大規模化し計算負荷も増大するが、本発明に係る改良型ＡＲＴを用いているので、光強度分布の算出精度、および光断層情報生成装置で得られる再構成画像の精度を向上させることができる。
【０１４７】
なお、以上説明した本発明に係る各実施の形態は、本発明の一例を示すものであり、本発明の構成を限定するものではない。本発明に係る光断層情報生成装置は、上記各実施の形態に限定されず、本発明の目的を逸脱しない範囲で種々変更して実施することが可能である。
【０１４８】
例えば、光拡散方程式は蛍光以外の光であっても同様に適用できるので、光断層情報生成装置１００、２００に限らず、励起光を検体１８に照射し、検体１８から射出される蛍光以外の光を検出し、その光の強度に基づいて断層画像を再構成する任意の構成の光断層情報生成装置の構成としても良い。将来的には、励起光を照射することなく自然発光する光に対して本発明の光断層情報生成装置等を適用することも考えられる。
【０１４９】
また、実施の形態１では、光源ヘッド２４と受光ヘッド２６とが機枠２２と共に機械的に回転する場合を例にとって説明したが、逆問題を１系統で処理するという目的においてはこの構成に限定される必要はなく、例えば、光照射機能と受光機能との双方の機能を有した複数のヘッドを位置固定で配置し、測定中は所定の回転方向に光照射機能と受光機能とを順々に遷移させて蛍光等を測定するような構成としても良い。
【０１５０】
また、実施の形態２では、ダイクロイックミラー２１３を設け、検体１８から射出される蛍光Ｌ１の成分と励起光Ｌ２の成分とを切り分けて、蛍光用の受光ヘッド２６と励起光用の受光ヘッド２１１とにそれぞれ入射する態様を例にとって説明したが、例えば、ダイクロイックミラー２１３を設けずに、受光ヘッド２６と受光ヘッド２１１とをより密に交互に配置することによって受光位置（角度）の違いによる信号差の影響を少なくするような構成としても良い。
【０１５１】
また、本質的には、本発明の適用対象画像は断層画像でなくても良い。光拡散方程式に基づく逆問題を解いて再構成画像を取得するような態様であれば良い。
【０１５２】
また、本発明に係る光強度分布算出方法のアルゴリズムをプログラミング言語によって記述し、必要に応じてコンパイルし、この光強度分布算出プログラムをメモリ（記録媒体）に記憶して他の光断層情報生成装置の情報処理手段によって実行させれば、本発明に係る光断層情報生成装置と同様の機能を実現することができる。
【産業上の利用可能性】
【０１５３】
本発明に係る光断層情報生成装置、光強度分布算出方法および光強度分布算出プログラムは、例えば、励起光によって発光する発光体の分布情報を用いて光断層画像を取得可能とする装置、特に、蛍光トモグラフィ装置の用途に利用することができる。
【符号の説明】
【０１５４】
１００、２００ 光断層情報生成装置
１２、２１０ 計測部
１４、１４ａ 画像処理部
１６ モニタ
１８ 検体
２４ 光源ヘッド
２６、２１１ 受光ヘッド
６２ 蛍光体
７０ 読取部
７２ 演算処理部
７４ 評価部
７６ 更新処理部
７８ 断層情報生成部
８０ 断層画像生成部
【特許請求の範囲】
【請求項１】
光拡散方程式に基づき、検体から射出される光の強度分布の計算値を求め、当該計算値と実測値との差の大小に応じて、ＡＲＴ（代数的再構成技術）を用いた再計算を行って前記計算値の更新値を得、当該更新値を用いて前記検体の光断層情報を生成する光断層情報生成装置であって、
前記検体の光学特性値から、当該光学特性値の変化に対する前記光強度分布の変化の割合を示すヤコビ行列を算出するヤコビ行列算出手段と、
前記ヤコビ行列を特異値分解して対角行列を得る特異値分解手段と、
前記対角行列の閾値以下の対角成分を０で置換した単位対角行列を取得する単位対角行列取得手段と、
前記ＡＲＴを用いた再計算において、ＡＲＴの逐次近似式に現れる前記計算値と実測値との差に対し前記単位対角行列を用いた置換を施して近似を行い、前記更新値を得る逐次近似手段と、
を具備する光断層情報生成装置。
【請求項２】
前記ヤコビ行列算出手段は、前記検体の光学特性値からヤコビ行列Ｊを算出し、
前記特異値分解手段は、前記ヤコビ行列ＪをＪ＝ＵＤＶＴと特異値分解して対角行列Ｄを得、
前記単位対角行列取得手段は、前記対角行列Ｄの閾値以下の対角成分を０で置換して単位対角行列Ｈを取得し、
前記逐次近似手段は、前記ＡＲＴの逐次近似式に現れる前記計算値と実測値との差ΔｙをＵＨＵＴΔｙとする置換を施して近似を行う、
請求項１記載の光断層情報生成装置。
【請求項３】
前記単位対角行列取得手段は、
前記ヤコビ行列のランク数の変化に対する特異値の変化を示す特性曲線において、ステップ状に変化する位置の特異値を前記閾値として用いる、
請求項２記載の光断層情報生成装置。
【請求項４】
光拡散方程式に基づき、検体から射出される光の強度分布の計算値を求め、当該計算値と実測値との差の大小に応じて、ＡＲＴを用いた再計算を行って前記計算値を更新する光強度分布算出方法であって、
前記検体の光学特性値から、当該光学特性値の変化に対する前記光強度分布の変化の割合を示すヤコビ行列を算出するステップと、
前記ヤコビ行列を特異値分解して対角行列を得るステップと、
前記対角行列の閾値以下の対角成分を０で置換した単位対角行列を取得するステップと、
前記ＡＲＴを用いた再計算において、ＡＲＴの逐次近似式に現れる前記計算値と実測値との差に対し前記単位対角行列を用いた置換を施して近似を行い、前記計算値の更新値を得るステップと、
を具備する光強度分布算出方法。
【請求項５】
光拡散方程式に基づき、検体から射出される光の強度分布の計算値を求め、当該計算値と実測値との差の大小に応じて、ＡＲＴを用いた再計算を行って前記計算値を更新する光強度分布算出方法であって、
前記検体の光学特性値から、当該光学特性値の変化に対する前記光強度分布の変化の割合を示すヤコビ行列を算出するステップと、
前記ヤコビ行列を特異値分解して対角行列を得るステップと、
前記対角行列の閾値以下の対角成分を０で置換した単位対角行列を取得するステップと、
前記ＡＲＴを用いた再計算において、ＡＲＴの逐次近似式に現れる前記計算値と実測値との差に対し前記単位対角行列を用いた置換を施して近似を行い、前記計算値の更新値を得るステップと、
をコンピュータに実行させる光強度分布算出プログラム。
【請求項１】
光拡散方程式に基づき、検体から射出される光の強度分布の計算値を求め、当該計算値と実測値との差の大小に応じて、ＡＲＴ（代数的再構成技術）を用いた再計算を行って前記計算値の更新値を得、当該更新値を用いて前記検体の光断層情報を生成する光断層情報生成装置であって、
前記検体の光学特性値から、当該光学特性値の変化に対する前記光強度分布の変化の割合を示すヤコビ行列を算出するヤコビ行列算出手段と、
前記ヤコビ行列を特異値分解して対角行列を得る特異値分解手段と、
前記対角行列の閾値以下の対角成分を０で置換した単位対角行列を取得する単位対角行列取得手段と、
前記ＡＲＴを用いた再計算において、ＡＲＴの逐次近似式に現れる前記計算値と実測値との差に対し前記単位対角行列を用いた置換を施して近似を行い、前記更新値を得る逐次近似手段と、
を具備する光断層情報生成装置。
【請求項２】
前記ヤコビ行列算出手段は、前記検体の光学特性値からヤコビ行列Ｊを算出し、
前記特異値分解手段は、前記ヤコビ行列ＪをＪ＝ＵＤＶＴと特異値分解して対角行列Ｄを得、
前記単位対角行列取得手段は、前記対角行列Ｄの閾値以下の対角成分を０で置換して単位対角行列Ｈを取得し、
前記逐次近似手段は、前記ＡＲＴの逐次近似式に現れる前記計算値と実測値との差ΔｙをＵＨＵＴΔｙとする置換を施して近似を行う、
請求項１記載の光断層情報生成装置。
【請求項３】
前記単位対角行列取得手段は、
前記ヤコビ行列のランク数の変化に対する特異値の変化を示す特性曲線において、ステップ状に変化する位置の特異値を前記閾値として用いる、
請求項２記載の光断層情報生成装置。
【請求項４】
光拡散方程式に基づき、検体から射出される光の強度分布の計算値を求め、当該計算値と実測値との差の大小に応じて、ＡＲＴを用いた再計算を行って前記計算値を更新する光強度分布算出方法であって、
前記検体の光学特性値から、当該光学特性値の変化に対する前記光強度分布の変化の割合を示すヤコビ行列を算出するステップと、
前記ヤコビ行列を特異値分解して対角行列を得るステップと、
前記対角行列の閾値以下の対角成分を０で置換した単位対角行列を取得するステップと、
前記ＡＲＴを用いた再計算において、ＡＲＴの逐次近似式に現れる前記計算値と実測値との差に対し前記単位対角行列を用いた置換を施して近似を行い、前記計算値の更新値を得るステップと、
を具備する光強度分布算出方法。
【請求項５】
光拡散方程式に基づき、検体から射出される光の強度分布の計算値を求め、当該計算値と実測値との差の大小に応じて、ＡＲＴを用いた再計算を行って前記計算値を更新する光強度分布算出方法であって、
前記検体の光学特性値から、当該光学特性値の変化に対する前記光強度分布の変化の割合を示すヤコビ行列を算出するステップと、
前記ヤコビ行列を特異値分解して対角行列を得るステップと、
前記対角行列の閾値以下の対角成分を０で置換した単位対角行列を取得するステップと、
前記ＡＲＴを用いた再計算において、ＡＲＴの逐次近似式に現れる前記計算値と実測値との差に対し前記単位対角行列を用いた置換を施して近似を行い、前記計算値の更新値を得るステップと、
をコンピュータに実行させる光強度分布算出プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図１０】
【図１２】
【図１３】
【図１４】
【図８】
【図９】
【図１１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図１０】
【図１２】
【図１３】
【図１４】
【図８】
【図９】
【図１１】
【公開番号】特開２０１１−２１５０３５（Ｐ２０１１−２１５０３５Ａ）
【公開日】平成２３年１０月２７日（２０１１．１０．２７）
【国際特許分類】
【出願番号】特願２０１０−８４３８０（Ｐ２０１０−８４３８０）
【出願日】平成２２年３月３１日（２０１０．３．３１）
【出願人】（３０６０３７３１１）富士フイルム株式会社 (25,513)
【Ｆターム（参考）】
【公開日】平成２３年１０月２７日（２０１１．１０．２７）
【国際特許分類】
【出願日】平成２２年３月３１日（２０１０．３．３１）
【出願人】（３０６０３７３１１）富士フイルム株式会社 (25,513)
【Ｆターム（参考）】
[ Back to top ]