JP2014085188A - 蛍光像再構成方法及び装置 - Google Patents

Abstract

【課題】蛍光像再構成において、計算コストを下げながら精度のよい結果を得る。
【解決手段】モデル形状について適当な離散化距離のもとで吸収・散乱係数を異ならせて予め準備された複数のシステム行列を用い、以下の工程（Ａ）から（Ｄ）を備えて試料中の蛍光光源の分布を再構成する。
（Ａ）試料を実測してその三次元表面形状である実測形状を取得する実測形状取得工程、
（Ｂ）前記モデル形状と前記実測形状との間の形状誤差を算出する形状誤差算出工程、
（Ｃ）前記形状誤差に基づいて前記離散化距離を修正して修正離散化距離を得るとともに、前記複数のシステム行列のうちからの選定又は修正により最適システム行列を得る補正工程、及び、
（Ｄ）前記試料を実測して得た表面輝度分布と、前記補正工程で得られた前記修正離散化距離及び最適システム行列とに基づいて、前記試料中の蛍光光源の分布を再構成する蛍光像再構成工程。
【選択図】図２

Description

本発明は試料表面に現れた蛍光像（表面輝度分布）から試料内部の蛍光像（蛍光光源の分布）を再構成する蛍光像再構成方法とその方法を実現する装置に関するものである。試料は光透過性の物体内に蛍光体が存在するものであり、特に限定されるものではないが、例えば小動物などの生体試料である。生体試料を扱う蛍光像再構成方法に関する技術は蛍光バイオイメージング技術と称される。

蛍光画像取得装置を利用した蛍光像再構成技術は、理論的な物理モデルを用いて励起光と蛍光の伝搬解析を行う過程(順問題)と、測定した試料表面の蛍光画像の実測データに順問題の過程で得られた伝搬解析結果を適用して試料内部の蛍光体分布を推定する過程(逆問題)から成り立つ。

順問題を実施する上での課題の1つとして、境界条件として与える試料の形状が挙げられる。試料が円筒や直方体といった単純な形状でない場合、例えばマウスなどの生体試料で複雑形状をもつものの場合は、光伝搬解析を行う際にその試料形状を忠実に反映させなければ計算の精度が著しく低下する。一方、試料形状を取得する毎に複雑な形状の条件下で光伝搬解析を行うと、計算コストが膨大になる。

逆問題での解の安定性は順問題での計算精度に依存することから、順問題においていかに試料形状に即した光伝搬計算を短時間で実施するかが、蛍光像再構成技術の性能を決定付ける重要な鍵となる。

このような試料形状由来の問題点に対して、現在２つのシステムが適用されている。 1つは、試料表面の蛍光画像とともにＸ線ＣＴや光切断法などで試料の外形形状を測定するシステム（Ａ）であり、得られた形状を基に順問題の計算を行う（非特許文献１）。このシステムでは、試料形状を忠実に再現できるため順問題の計算精度は高いが、試料毎にその外形形状に基づいた順問題計算を実行する必要があるため、計算コストが膨大となる問題点がある。

もう１つは、試料形状を円筒形状や平板形状と仮定し、試料の形状(大きさを含む)を固定するシステム（Ｂ）であり、固定形状を基に順問題の計算を行う（特許文献１、非特許文献２参照。）。このシステムでは試料を型に固定して所定の形状に当てはめるため試料毎にその形状に基づく順問題計算を実行する必要がないため計算コストは軽減されるが、形状誤差によって計算精度が低下する問題や、試料測定時に試料を所定の型に固定することによる生理機能への影響は避けられない。

特開２０１０−１９７３８１号公報 ＷＯ２００８／０５９５７２ ＷＯ２００９／１３９０５８

FLECT: A True 360°FLuorescence Emission Computed Tomography System(2011 World Molecular Imaging Congress) Parkin Elmer社蛍光・定量トモグラフィ装置FMT、http://www.primetech.co.jp/search/detail.asp?pid=10086 Quantitative fluorescence tomography with functional and structural a priori information.(Lin Y, Yan H. Nalcioglu O, Gulsen G. Appl Opt. 2009 Mar l;48(7):1328-36.)

上述の順問題と逆問題計算から構成された蛍光像再構成技術を装置に組み込む場合、実用的な計算時間内で精度の良い結果を得ることが望ましい。

本発明は、既存のシステム（Ａ）と（Ｂ）の利点をうまく組み合わせて、計算コストを下げながら精度のよい結果を得ることのできる蛍光像再構成技術を提供することを目的とする。

本発明では、まず初めに、観測試料のモデル形状について予め順問題計算である光伝搬計算を実施する。モデル形状についての光伝搬計算では、適当な離散化距離のもとで吸収・散乱係数を異ならせて複数のシステム行列を求めておく。システム行列については後述する。対象試料は測定に用いる励起光と蛍光を透過させる物質の内部に蛍光物質が蛍光光源として存在するものであり、典型的には生体試料であるが、それに限らない。生体試料としては、典型的にはマウスやラットなどがあげられる。例えばマウスを測定対象とする場合には、モデル形状も標準試料として選んだマウスである。つまり、順問題計算を実施する段階では、システム（Ｂ）のように試料形状はモデル形状に固定されているものとみなして、モデル形状について順問題計算を行う。

次に、対象試料について、システム（Ａ）のような手法でその対象試料の実測形状を取得する。その実測形状により光伝搬計算を行うのではなく、実測形状とモデル形状との間の形状誤差を計算する。そして、先に求めたモデル形状による光伝搬計算結果をその形状誤差に基づいて補正処理を行ったものをその実測形状に対する光伝搬計算結果とする。補正処理では、形状誤差に合わせて吸収・散乱係数を変更することによって、実測形状に対応したシステム行列（最適システム行列という）を求める。言い換えれば、吸収・散乱係数を異ならせた複数のシステム行列から、形状誤差分を修正できるシステム行列を選択することで最適システム行列を求める。この際、形状誤差分に合わせて離散化距離を修正する。

このように、対象試料について、その都度、システム（Ａ）のような手法でその実測形状を取得して光伝搬計算を行っていると計算量が膨大になってコスト高を招くが、予め計算により求めておいたモデル形状についての光伝搬計算結果を補正するだけであればその計算量は少なくてすむため計算コストを削減することができる。

すなわち、本発明の蛍光像再構成方法は、モデル形状について適当な離散化距離のもとで光学定数を異ならせて予め準備された複数のシステム行列を用い、以下の工程（Ａ）から（Ｄ）を備えて試料中の蛍光光源の分布を再構成する。光学定数は吸収係数と散乱係数であり、光学定数に代えて吸収・散乱係数ということもある。
（Ａ）試料を実測してその三次元表面形状である実測形状を取得する実測形状取得工程、
（Ｂ）前記モデル形状と前記実測形状との間の形状誤差を算出する形状誤差算出工程、
（Ｃ）前記形状誤差に基づいて前記離散化距離を修正して修正離散化距離を得るとともに、前記複数のシステム行列のうちからの選定又は修正により最適システム行列を得る補正工程、及び、
（Ｄ）前記試料を実測して得た表面輝度分布と、前記補正工程で得られた前記修正離散化距離及び最適システム行列とに基づいて、前記試料中の蛍光光源の分布を再構成する蛍光像再構成工程。

好ましい一形態によれば、形状誤差算出工程（Ｂ）はモデル形状の断層画像と実測形状の断層画像について、それぞれの重心からそれぞれの断層画像輪郭までの同一方位の距離に基づいて形状誤差を求める。

好ましい他の形態によれば、補正工程（Ｃ）で得られる最適システム行列は、前記形状誤差算出工程（Ｂ）で算出された形状誤差により補正された吸収・散乱係数に一致するか又は最も近い吸収・散乱係数に対応したシステム行列として選定されたシステム行列である。

さらに他の好ましい他の形態によれば、補正工程（Ｃ）で得られる最適システム行列は、前記形状誤差算出工程（Ｂ）で算出された形状誤差により補正された吸収・散乱係数に近い２つの吸収・散乱係数に対応した２つのシステム行列から補完法により修正して得られたシステム行列である。

さらに他の好ましい他の形態によれば、モデル形状は形状の異なる複数個を含んでおり、モデル形状ごとに複数のシステム行列が予め準備されており、形状誤差算出工程（Ｂ）では全てのモデル形状と実測形状との間の形状誤差を算出し、補正工程（Ｃ）を形状誤差の最も小さいモデル形状について行う。

さらに他の好ましい他の形態によれば、モデル形状についての予め準備された複数のシステム行列は、モデル形状の内部形態情報を基にして組織毎に準備されている。

本発明の蛍光像再構成方法を実施する本発明の蛍光像再構成装置の好ましい一形態は、図１に示されるように、試料の少なくとも外表面の三次元形状を取得する三次元形状測定装置１２、試料表面に現れた蛍光像を取得する蛍光画像取得装置１４、及び三次元形状測定装置１２が取得した試料の三次元形状を基に励起光と蛍光の伝搬解析を行なうとともに、その伝搬解析と蛍光画像取得装置１４が取得した試料表面の蛍光像から試料内部の蛍光像を再構成する画像処理装置１６を備えている。そして、画像処理装置１６は、観測試料に関連する標準試料についての三次元形状測定装置１２による測定結果の三次元形状をモデル形状として保持するモデル形状保持部１８と、モデル形状保持部１８に保持されているモデル形状について光伝搬計算を行う光伝搬解析部２０と、光伝搬解析部２０による光伝搬計算結果としてモデル形状について適当な離散化距離のもとで吸収・散乱係数を異ならせて作成した複数のシステム行列を保持する光伝搬解析結果保持部２２と、観測試料についての三次元形状測定装置１２による測定結果の三次元形状とモデル形状保持部１８に保持されているモデル形状との間の形状誤差を算出する形状誤差算出部２４と、形状誤差に基づいて離散化距離を修正して修正離散化距離を得る離散化修正部２６と、光伝搬解析結果保持部２２に保持されている複数のシステム行列のうちから形状誤差に基づいて選定又は修正により最適システム行列を得る光伝搬解析結果補正部２８と、蛍光画像取得装置１４により取得された観測試料表面の輝度分布、離散化修正部２６で得られた修正離散化距離、及び光伝搬解析結果補正部２８で得られた最適システム行列から観測試料中の蛍光光源の分布を再構成する蛍光像再構成部３０と、を備えている。

本発明によれば、蛍光像再構成の順問題計算において、外形形状に対する補正処理を行うことで計算コストを軽減した状態で計算精度を向上させる。これによって、逆問題における解の向上を実現する。

蛍光像再構成装置の一実施例を示すブロック図である。 一実施例の動作を示すフローチャートである。 同実施例における形状誤差算出方法を示す断面画像である。 無限媒体における蛍光フルエンスレート分布を示す画像（ａ）とグラフ（ｂ）である。 図４の蛍光フルエンスレート分布における離散化と光学定数の補正を示すグラフである。 励起光順問題における形状誤差による補正処理を示す図である。 蛍光順問題における形状誤差による補正処理を示す図である。 励起光／蛍光順問題における補正処理の効果を検証するためのモデルを示す図である。 他の実施例の動作を示すフローチャートである。 同実施例における生体内部の組織毎の光学定数設定の例を示す端面画像である。 同実施例における形状誤差算出方法を示す断面画像である。

蛍光像再構成装置の一実施形態を示す図１において、三次元形状測定装置１２としては、特に限定されるものではないが、例えばＸ線ＣＴ、ＭＲＩ（磁気共鳴画像装置）、光切断法による画像取得装置、特許文献２に示される生体画像取得装置、特許文献３に示される生体画像撮影装置などを使用することができる。三次元形状として断面構造まで知る必要があるときはＸ線ＣＴやＭＲＩを使用する。

蛍光画像取得装置１４としては、これも特に限定されるものではないが、例えば特許文献２に示される生体画像取得装置、特許文献３に示される生体画像撮影装置などを使用することができる。

画像処理装置１６はコンピュータであり、この蛍光像再構成装置の専用コンピュータ又は汎用のパーソナルコンピュータである。モデル形状保持部１８、光伝搬解析部２０、光伝搬解析結果保持部２２、形状誤差算出部２４、離散化修正部２６、光伝搬解析結果補正部２８及び蛍光像再構成部３０は、そのコンピュータにより実現される機能である。

画像処理装置１６には蛍光像再構成部３０で作成した蛍光像を表示する液晶表示装置（ＬＣＤ）などの表示装置３２が接続されている。

まず、蛍光像再構成装置において、順問題解析と逆問題解析から断層画像を計算する処理手順を説明する。順問題解析では、任意の方向から励起光を試料へ照射した場合に励起光が試料内部を伝搬する過程、及び試料内部の任意の位置に蛍光剤が存在した場合に蛍光が試料内部を伝搬する過程から成り立つ。励起光が試料内部を伝搬する過程は次に示す光拡散方程式（１）によって計算することができる。
また、蛍光が試料内部を伝搬する過程は次に示す光拡散方程式（２）によって計算することができる。
ここで、
Ｄ：拡散定数、
ｒ：位置、
μａ：吸収係数、
φex(j)：ｊ番目の方向から励起光を与えた場合の位置ｒにおける励起光フルエンスレート（光密度）、
φem(k)：ｋ番目の位置に蛍光剤が存在した場合の位置ｒにおける蛍光フルエンスレート、
Ｓ：光源強度、
ε：モル吸光係数、
γ：量子収率、
Ｍ：モル濃度、を表す。

拡散定数Ｄは１／３（μｓ'＋μａ）と表わされ、μｓ'は等価散乱係数であり、μｓ'＝（μｓ（１−ｇ））と表わされる。μｓは散乱係数、ｇは非等方散乱係数と呼ばれる係数で、散乱の異方性を示しており、位相関数pの余弦平均として下記式で定義される。
p(θ)は散乱位相関数を表す。一般に、ｇ≒０だと等方散乱、ｇ＝０．９〜１だと前方散乱を示す。ｇは試料の種類によって定める。

（１）式と（２）式から励起・蛍光伝搬の理論値を求めることができ、この計算結果を元にシステム行列Ａを作成する。システム行列Ａとは、蛍光剤の空間分布ｆを与えた場合に試料表面で検出される蛍光分布の理論値を得るための行列である。具体的には、システム行列の列ベクトルａは、任意の位置における試料表面で検出される蛍光分布の理論値であり、この列ベクトルをｋ番目の位置までそれぞれ求めたのがシステム行列である。
ｇ＝Ａｆ （３）
Ａ＝［ａ（１），ａ（２），…，ａ（ｋ）］ （４）
ここで、
Ａ：システム行列、
ｆ：蛍光剤の空間分布ベクトル、
ｇ：試料表面での蛍光分布ベクトル、
を表す。例えば、蛍光剤の空間分布ｆ＝［０ ０ ０.５ １ ０.５ ０ … ０］^Tが与えられた場合、このベクトルにシステム行列を掛ける（Ａｆ）と、試料表面での蛍光分布ベクトルｇが得られる。ここで、［…］^Tは転置を示している。

逆問題解析での演算は、蛍光画像測定装置で得られた蛍光検出データと順問題解析で計算した理論値から蛍光剤の空間分布を求める過程から成り立つ。まず、順問題解析によって（４）式におけるシステム行列Ａを求めることができる。一方、蛍光画像測定装置によって（３）式における試料表面での蛍光分布ベクトルｇを得ることができる。逆問題解析では、順問題解析で求められたシステム行列Ａと試料表面での蛍光分布ベクトルｇから、未知の蛍光剤の空間分布ベクトルfを求める。ベクトルｆの求め方の一般的な方法は最小二乗法であり、例えば以下の（５）式の評価関数を最小化することで求めることができる。

図１のブロック図と図２のフローチャートを参照して、本発明を実施するための形態を説明する。ここでは、試料として生体試料を例にして説明する。

まず初めに、モデル形状を設定する（図２−ステップ１）。この操作は次のように行う。観測試料、例えばマウスやラット等、の標準試料を１又は複数用意し、その形状をＸ線ＣＴや光切断法などの三次元形状測定装置１２で精度よく測定する。この形状取得は好ましくは複数の標準試料について行う。例えば、マウスについて測定しようとする場合は、体型の異なる複数のマウスを標準試料として用意し、それらの複数のマウスについて三次元形状を測定する。そして、取得した複数の標準試料マウスの形状データをモデル形状（１〜Ｎ）として、モデル形状保持部１８に記憶させて保持しておく。これがモデル形状の設定である。

次に、各モデル形状について、適当な大きさのメッシュに離散化を行い、各メッシュの立方格子点毎に光伝搬解析部２０で光伝搬解析を行って光伝搬解析結果を求める（図２−ステップ２）。光伝搬解析結果はシステム行列である。システム行列の求め方は上述の通りである。ここでの離散化は、特に限定されるものではなく、モデル形状内の蛍光物質分布を適度に近似できる点数のデータが確保できるものであればよい。例えは、メッシュの大きさを１ｍｍとする。ここで、光学定数（吸収係数μａ及び散乱係数μs）を変えた条件での光伝搬解析も行う。例えば、あるモデル形状での光学定数をμａ＝０．０２［／ｍｍ］、μｓ＝１０［／ｍｍ］と設定した場合、両光学定数の値を適当な間隔で変化させた複数の光学定数の条件、例えば０．４、０．６、０．８、１．２、１．４、１．６倍した光学定数の条件、で（１），（２）式を解いてそれらの光学定数での光伝搬解析結果としてのシステム行列も求めておく。このようにして、各モデル形状について複数種類に光学定数を異ならせて求めたそれぞれのシステム行列を光伝搬解析結果保持部２２に記憶させて保持しておく。

ここまでの光伝搬解析は三次元形状に基づいて行うので計算に時間がかかるが、それはモデル形状に関して行うだけであり、観測試料に関してはこのような光伝搬解析計算は行わないので、この蛍光像再構成方法全体に占める光伝搬解析計算時間の占める割合は小さくなる。

次に、三次元形状測定装置１２に観測試料を設置し（図２−ステップ３）、実測形状として観測試料の三次元形状を取得する（図２−ステップ４）。ここまでの処理が終わった段階で、モデル形状と実測形状の二つの形状情報を得ることになる。

次に、形状誤差算出部２４は、モデル形状に対して実測形状がどの程度誤差を含むかを示す形状誤差を画像解析によって求める（図２−ステップ５）。

形状誤差の求め方の具体例を図３に示す。まず、１番目のモデル形状と実測形状について断層画像を取り出す。断層画像は、例えばモデル形状と実測形状を同じ向きに配置したときの断面画像とする。モデル形状と実測形状を配置する方向は特に限定されるものではなく、試料の方向を特定できて、モデル形状と実測形状の対応する部位の断面形状が比較できるようなものであればよい。例えば、生体試料を例とすると頭部から尾部に向かう生体軸の方向を挙げることができる。

そして、１番目のモデル形状の断層画像と実測形状の断層画像（図３−（ａ））から、それぞれの二値化画像（図３−（ｂ））を求める。二値化された画像を利用してそれぞれの画像の面積の重心点（図３−（ｃ））を求める。得られた重心情報から重心と外郭の距離を算出し、それぞれの距離をＤｍｉ（モデル）、Ｄｄｉ（実測）として形状誤差αｉ＝Ｄｄｉ／Ｄｍｉを求める。このとき、ＤｍｉとＤｄｉを求める外郭の位置ｉは、モデル形状と実測形状において重心からみて同じ方位になる位置である。形状誤差αｉは外郭の複数の位置（ｉ〜ｎ）で求め、生体表面の細かな凹凸を含めて、形状誤差αｉの平均値αを求める。これを１〜Ｎ番目のモデル形状に対して実施し、最も１に近いα値を求める。ここでは、モデル形状が複数個設定されている場合を示しているが、モデル形状が１個だけ設定されている場合はその１個のモデル形状について実測形状との形状誤差αｉの平均値αを求めるだけでよい。

「形状誤差」αは（Ｄｄ／Ｄｍ）として定義される関数であるので、αが１に近いほど形状誤差が小さいことを意味する。

本発明では、生体試料の個体差のみを対象とするため、形状誤差αはモデル形状に対する実測形状の大きさの比率であると考え、形状誤差αは拡大又は縮小で説明できるものとする。

次に、得られた形状誤差の情報を基に、離散化修正部２６ではモデル形状の離散化距離の修正（図２−ステップ６）を行い、光伝搬計算結果補正部２８では光伝搬計算結果の補正（図２−ステップ７）を行う。

具体例を図４と図５に示す。無限媒体に蛍光剤（強度：ｌｍＷ）が存在し、半径２０ｍｍまでの光伝搬解析を行った場合を考える。無限媒体における光拡散方程式の解析解は（６）式のようになり、半径２０ｍｍまでフルエンスレートは図４−（ｂ）のように求めることができる。

ここで、図２−ステップ６、７の処理の一例として、半径３０ｍｍ地点でのフルエンスレート（未知）を求めることを考える。図４−（ｂ）は離散化距離（ｌｍｍ）の下でのフルエンスレートの計算結果を表わしている。まず、その計算結果を離散化距離のみα倍（３０／２０＝１．５倍）する。これが離散化距離の修正である。図５−（ａ）は図４−（ｂ）と同じものを示しており、それの離散化距離の修正により、図５−（ｂ）のように半径３０ｍｍまでのフルエンスレートが求まる。しかし、この場合、離散化距離がα倍されただけなので、その求められたフルエンスレートは半径２０ｍｍまでのフルエンスレートが元の半径２０ｍｍまでのフルエンスレートとは異なる。

さらに光伝搬解析の補正を行う。散乱係数μｓ［／ｍｍ］、吸収係数μａ［／ｍｍ］はそれぞれ１ｍｍあたりのフォトンの散乱度合い、吸収度合いを示すパラメータである。そのため、離散化距離をα倍するのに伴って光学定数（吸収係数μａと散乱係数μｓ）もα倍すると、図５−（ｃ）のようにフルエンスレートが求まり、そのようにして求められたフルエンスレートは半径２０ｍｍまでのフルエンスレートが元の判例２０ｍｍまでのフルエンスレートと一致する。

図２−ステップ２において、異なる光学定数での光伝搬計算結果（システム行列）は予め求められて光伝搬解析結果保持部２２に保持されているので、光伝搬解析結果の補正の一例は、α倍された光学定数に対応するシステム行列を光伝搬解析結果保持部２２から呼び出してくることである。α倍された光学定数に一致する光学定数に対応するシステム行列がない場合は、α倍された光学定数に最も近い光学定数に対応するシステム行列を光伝搬解析結果保持部２２から呼び出してくるようにしてもよい。

さらに、α倍された光学定数に一致する光学定数に対応するシステム行列がない場合の補正の他の例は、形状誤差により補正された光学定数に近い２つの光学定数に基づいてそれぞれ算出された２つのシステム行列を光伝搬解析結果保持部２２から呼び出し、形状誤差により補正された光学定数に対応したシステム行列になるように補完法によりシステム行列を補正することである。そのような補完法の一例としてスプライン補完法を挙げることができる。補完法以外の計算方法により形状誤差により補正された光学定数に対応したシステム行列を導き出すようにしてもよい。

このように、図２−ステップ６、７の補正処理は、形状誤差により補正された光学定数に対応したシステム行列（最適システム行列）を光伝搬解析結果保持部２２から呼び出すか、又は呼び出したシステム行列から補完法等の計算により形状誤差により補正された光学定数に対応したシステム行列を導き出すものであるので、実測形状に基づいてシステム行列を実施するのに比べると順問題の計算簡略化が可能となる。

離散化修正部２６で形状誤差により修正された修正離散化距離が求められ、光伝搬解析結果補正部２８で形状誤差により補正された最適システム行列が求められるので、蛍光像再構成部３０は蛍光画像測定装置１４から観測試料表面での蛍光分布を得て逆問題解析を行い、観測試料内の蛍光剤の空間分布を計算により再構成する。蛍光像再構成部３０での逆問題解析方法は前述のとおりである。

次に、生体を模擬した円筒媒体に対して、一実施例を適用した例を図６〜７、表１に示す。

図６では励起光の順問題の計算を実施した。まず、直径２５ｍｍの円筒媒体において、離散化距離を例えば１ｍｍとして、有限要素法によって光拡散方程式の数値解を求めた（図６−（ａ））。ここで、実際の形状が直径３５ｍｍの円筒媒体だったと想定して、一実施例による補正処理を行った。直径２５ｍｍの計算結果と形状誤差α＝３５／２５＝１.４を利用して、離散化距離と光学定数（吸収係数μａと散乱係数μsの両方）を１.４倍にするように修正することで、図６−（ｂ）の結果を得た。

補正処理を行った結果の妥当性について検討するため、実際の直径３５ｍｍのサイズで計算を実施して光拡散方程式の数値解を求めた（図６−（ｃ））、図６−（ｃ）と図６−（ｂ）の結果を比較したところ、両者の結果が良く一致することを確認した。

なお、有限媒体の場合、φ→φ／α²のように値を更新することにより絶対値がよく一致することを確認している。φ→φ／α²に値を更新するとはシステム行列の各要素を１／α²倍することである。システム行列は（１）式と（２）式から理論的に求めることができ(φの記述)、光源強度を補正するために1/α²倍の処理を加えることである

図７では蛍光の順問題の計算を実施した。まず、直径２５ｍｍの円筒媒体において、離散化距離を例えば１ｍｍとして、有限要素法によって光拡散方程式の数値解を求めた（図７−（ａ））。ここで、実際の形状が直薩３５ｍｍの円筒媒体だったと想定して、一実施例による補正処理を行った。直径２５ｍｍの計算結果と形状誤差α＝３５／２５＝１．４を利用して、離散化距離と光学定数（吸収係数μａと散乱係数μsの両方）を１.４倍にするように修正することで、図７−（ｂ）の結果を得た。

補正処理を行った結果の妥当性について検討するため、実際の直径３５ｍｍのサイズで計算を実施して光拡散方程式の数値解を求めた（図７−（ｃ））、図７−（ｃ）と図７−（ｂ）の結果を比較したところ、両者の結果が良く一致することを確認した。

図６、７の結果を利用して、励起光／蛍光の順問題を実施し、補正処理の効果を検証した（図８）。

励起光源、検出器、蛍光剤ＩＣＧ（５(Indo-Cyanine Green）が図８のケース１、ケース２に示した位置にあると仮定し、本発明による補正処理を実施した結果を表１に示す。ケース１は蛍光剤が直径３５ｍｍの球の中心にある場合、ケース２は直径３５ｍｍの球の互いに近い位置に励起光入射位置と蛍光検出位置が設けられ、そこからほぼ反対側の位置（距離では３２.２ｍｍの位置）に蛍光剤がある場合を示している。

なお、蛍光剤ＩＣＧの濃度１［μＭ］、モル吸光係数１．２８×１０^-2［／ｍｍ・μＭ］、量子収率０．０１、円筒媒体の吸収係数０．０２［／ｍｍ］、散乱係数１０［／ｍｍ］、励起フルエンス（励起フルエンスレートのこと）は設定したＩＣＧの位置での値、蛍光フルエンス（蛍光フルエンスレートのこと）は設定した検出器付近での値とした。「Ｅ-10」という表示は「１０^-10」を示す。

表１の結果より、補正処理によって本来の解を近似できることが分かった。

以上により、モデル形状を用いた順問題での光伝搬解析結果を補正した結果が実試料の形状を用いて順問題を解いた結果とよく一致することを示した。したがって、そのように補正された離散化距離とシステム行列を用いて逆問題を解いて得られる実試料の蛍光光源分布の再構成結果も精度のよいものとなる。

本発明の他の実施例について説明する。順問題での計算において、特に試料が生体試料の場合、その形状と同様に内部組織の光学定数も計算精度を高める上で重要となる（非特許文献３参照。）。その際、ＭＲＩによる形態画像を利用すれば、組織に応じた光学定数を設定することができる。その一方、ＭＲＩ画像を蛍光測定毎に取得するのは多大な時間を要する。

そこで、この実施例は簡易に光学定数を補正する方法を提供するものである。具体例を図９〜図１１に示す。試料としては生体試料を採り上げる。

まず、生体試料をＭＲＩ撮像し（図９−ステップ１）、組織毎の光学定数（吸収係数μａと散乱係数μｓ）を設定した上で（図９−ステップ２、図１０）、光伝搬解析を行う（図９−ステップ３）。さらに吸収係数μａと散乱係数μｓを変えた条件でも光伝搬解析を行う。例えば、モデル形状で図１０のような光学定数を得た場合、各組織の光学定数の値を０．４、０．６、０．８、１．２、１．４、１．６倍した条件で光伝搬解析を行う。さらに、この結果を元に補間、例えばスプライン補間を利用することで、上記以外の光学定数についても結果を近似することができる。これらの結果は図９−ステップ７、ステップ８の補正計算で使用する。

一方、観察試料を設置して三次元形状を取得する（図９−ステップ４、ステップ５）。モデル形状と実測形状から形状誤差αを算出する（図９−ステップ６、図１１）。形状誤差αの算出は図３に示した方法と同様に行う。

この情報を元に、モデル形状の離散化及び組織毎の光学定数に対して補正処理を実施し（図９−ステップ７）、この情報を基に光伝搬計算の補正を行う。

以上の手法を利用することで、予め取得したモデル形状のＭＲＩ画像を利用して測定対象の生体試料における組織毎の光学定数の補正ができる。

外形に形状誤差があっても、内部組織同士には形状誤差がないといった場合や、その逆の場合も考えられる。そこで、図９−ステップ５の三次元測定をＸ線ＣＴやＭＲＩで撮像する場合は、次の補正手法を適用することで、内部形状の補正が実施できる。すなわち、形状全体の重心・形状誤差αとは別に、組織毎の重心を求めて、組織毎の形状誤差（β_1〜N）を算出する。また、組織が全体に占める割合をＷ_1〜N（＝０〜１）とおくと、求める修正形状誤差γは、
と求められる。図９−ステップ７、ステップ８での離散化修正、光学定数補正、光伝搬結果補正のための形状誤差αを修正形状誤差γに置き換えてそれらの処理を実行すれば、内部形状を考慮した補正が実施できる。

１２ 三次元形状測定装置
１４ 蛍光画像取得装置
１６ 画像処理装置
１８ モデル形状保持部
２０ 光伝搬解析部
２２ 光伝搬解析結果保持部
２４ 形状誤差算出部
２６ 離散化修正部
２８ 光伝搬解析結果補正部
３０ 蛍光像再構成部

Claims (7)

1. モデル形状について適当な離散化距離のもとで吸収・散乱係数を異ならせて予め準備された複数のシステム行列を用い、以下の工程（Ａ）から（Ｄ）を備えて試料中の蛍光光源の分布を再構成する蛍光像再構成方法。
   （Ａ）試料を実測してその三次元表面形状である実測形状を取得する実測形状取得工程、
   （Ｂ）前記モデル形状と前記実測形状との間の形状誤差を算出する形状誤差算出工程、
   （Ｃ）前記形状誤差に基づいて前記離散化距離を修正して修正離散化距離を得るとともに、前記複数のシステム行列のうちからの選定又は修正により最適システム行列を得る補正工程、及び、
   （Ｄ）前記試料を実測して得た表面輝度分布と、前記補正工程で得られた前記修正離散化距離及び最適システム行列とに基づいて、前記試料中の蛍光光源の分布を再構成する蛍光像再構成工程。
2. 前記形状誤差算出工程（Ｂ）は前記モデル形状の断層画像と前記実測形状の断層画像について、それぞれの重心からそれぞれの断層画像輪郭までの同一方位の距離に基づいて形状誤差を求める請求項１に記載の蛍光像再構成方法。
3. 前記補正工程（Ｃ）で得られる最適システム行列は、前記形状誤差算出工程（Ｂ）で算出された形状誤差により補正された吸収・散乱係数に一致するか又は最も近い吸収・散乱係数に対応したシステム行列として選定されたシステム行列である請求項１又は２に記載の蛍光像再構成方法。
4. 前記補正工程（Ｃ）で得られる最適システム行列は、前記形状誤差算出工程（Ｂ）で算出された形状誤差により補正された吸収・散乱係数に近い２つの吸収・散乱係数に対応した２つのシステム行列から補完法により修正して得られたシステム行列である請求項１又は２に記載の蛍光像再構成方法。
5. 前記モデル形状は形状の異なる複数個を含んでおり、モデル形状ごとに前記複数のシステム行列が予め準備されており、
   前記形状誤差算出工程（Ｂ）では全てのモデル形状と前記実測形状との間の形状誤差を算出し、
   前記補正工程（Ｃ）を形状誤差の最も小さいモデル形状について行う請求項１から４のいずれか一項に記載の蛍光像再構成方法。
6. モデル形状についての予め準備された前記複数のシステム行列は、モデル形状の内部形態情報を基にして組織毎に準備されている請求項１から５のいずれか一項に記載の蛍光像再構成方法。
7. 試料の少なくとも外表面の三次元形状を取得する三次元形状測定装置、試料表面に現れた蛍光像を取得する蛍光画像取得装置、及び前記三次元形状測定装置が取得した試料の三次元形状を基に励起光と蛍光の伝搬解析を行なうとともに、その伝搬解析と前記蛍光画像取得装置が取得した試料表面の蛍光像から試料内部の蛍光像を再構成する画像処理装置を備えた蛍光像再構成装置であって、
   前記画像処理装置は、
   観測試料に関連する標準試料についての前記三次元形状測定装置による測定結果の三次元形状をモデル形状として保持するモデル形状保持部と、
   前記モデル形状保持部に保持されているモデル形状について光伝搬計算を行う光伝搬解析部と、
   前記光伝搬解析部による光伝搬計算結果としてモデル形状について適当な離散化距離のもとで吸収・散乱係数を異ならせて作成した複数のシステム行列を保持する光伝搬解析結果保持部と、
   観測試料についての前記三次元形状測定装置による測定結果の三次元形状と前記モデル形状保持部に保持されているモデル形状との間の形状誤差を算出する形状誤差算出部と、
   前記形状誤差に基づいて離散化距離を修正して修正離散化距離を得る離散化修正部と、
   前記光伝搬解析結果保持部に保持されている前記複数のシステム行列のうちから前記形状誤差に基づいて選定又は修正により最適システム行列を得る光伝搬解析結果補正部と、
   前記蛍光画像取得装置により取得された前記観測試料表面の輝度分布、前記離散化修正部で得られた修正離散化距離、及び前記光伝搬解析結果補正部で得られた最適システム行列から前記観測試料中の蛍光光源の分布を再構成する蛍光像再構成部と、
   を備えている蛍光像再構成装置。