JP5794201B2 - 生体蛍光画像取得装置 - Google Patents

Description

本発明は小動物などの生体試料を対象とする蛍光バイオイメージング技術に関する。

生体蛍光画像取得装置には大きく分けて二つのシステムが存在する。一つは、短時間かつ効率的に分子プローブの定性解析を行うことを目的としたシステム（Ａ）である。システム（Ａ）では、生体試料の広範囲を照射して励起することが可能な光源とＣＣＤカメラなどの二次元検出器を有することによって、生体試料表面からの蛍光を検出する手法が一般的となっている。このシステムの特徴として、多検体の高速検査が可能となる一方、分子プローブの定量性や位置特定といった詳細な解析は困難な場合が多い。

これに対し、分子プローブの空間分布や濃度などを詳細に解析すること（蛍光断層画像化）を目的としたシステム（Ｂ）が存在する。システム（Ｂ）では、励起光源は生体試料に密着型の点光源とし、これを生体試料に対して多数配置して順次点灯していくことにより多点多方向照射として、生体内から発生する蛍光を観察することで断層化を行う手法が一般的となっている（特許文献１，２、非特許文献１参照。）。このシステムの特徴として、蛍光断層化によって蛍光を発する分子プローブの定量性や位置特定が可能となる一方、測定時間や断層化の計算時間が膨大となる場合が多い。

システム（Ｂ）における断層化の具体的な手法としては、生体に於ける励起光・蛍光の光伝搬解析（順問題）を行った上で、測定した実測データから内部の推定値を求める（逆問題）手法が用いられる。

光伝搬解析の具体的な解法としては、モンテカルロシミュレーションや光輸送方程式、光拡散方程式といった方法があるが、一般的には光拡散方程式が用いられる。任意の形状における光拡散方程式の解は、有限要素法等の数値解法によって求めることができる。

特開２００８−１４９１５４号公報 特開２０１０−１７５４６６号公報 国際公開第２００９／１３９０５８号

田嶋英朗他、ＤＯＩ−ＰＥＴ検出器を用いたＰＥＴ・蛍光ＣＴ同時イメージング装置における蛍光断層像画像再構成の検討、２００８年度日本医用画像工学会大会Ｐ１３，ｐ．４０，２００８ 島津評論 Ｖｏｌ．５７，Ｎｏ．１・２，ｐ．１４４（２０００．８）

生体機能の詳細な解析のためには、蛍光断層画像化が可能なシステム（Ｂ）を採用することが望ましいが、測定時間や断層化の計算時間が膨大となる問題点がある。

この課題解決として、システム（Ａ）とシステム（Ｂ）を組み合わせるアイデアが提案されている（特許文献３参照。）。即ち、一方向から広範囲に照射し励起して蛍光観察を行うとともに、その照射方向を複数の方向に切り換えるシステムである。このシステムが実現できれば、測定時間や断層化の問題点は解決するものと考えられる。

ただし、このアイデアには解決すべき課題がある。それは生体試料の形状に起因する問題であり、
（Ｉ）光拡散方程式を解く際の最適なメッシュサイズ（離散化）が分からない、
（II）広域照射の際に光源の照射ムラや生体試料の表面形状によって一様に励起照射されない、
といった課題が存在する。
ここで、メッシュサイについて説明する。光拡散方程式のような偏微分方程式の数値解法の一つに有限要素法がある。有限要素法とは、対象となる物理形状を要素(メッシュ)と呼ばれる小領域に分割(離散化)して、近似解を求める方法である。ここで、分割(離散化)の程度によってメッシュサイズが定義される。メッシュサイズを小さく設定すると、物理形状の離散化誤差が小さくなるため、解の精度は高くなることが一般的に知られている。

（Ｉ）に関しては、例えばモンテカルロシミュレーションでは、光が散乱するまでの距離（平均自由行程）を離散化の指標とし、平均自由行程以下を単位としてモデルを離散化する前例はあるが、光拡散方程式に関してはどの程度離散化すればよいかの指標は存在しない。例えば、メッシュのサイズを小さくすれば光拡散方程式を精度よく計算することができるが、それに伴いメッシュ数自体が増加することで計算時間が大幅に増加する問題点がある。この問題を解決する手法として、光源付近のメッシュのみ細かく設定するアダプティブメッシュ法が一般的に知られているが、メッシュの数だけ繰り返し演算を行う蛍光画像断層化においてはメッシュをその都度更新するアダプティブメッシュ法は逆に計算時間を増大させてしまう。したがって、計算に最適なメッシュサイズ（すなわちメッシュの数）を予め決定し、決定されたメッシュ構造の元で蛍光画像断層化を行うことが最も効率的な処理方法となる。

（II）に関しては、励起光源を点照射（試料密着型）から広範囲照射に変更するために、光源の照射ムラや生体試料の表面形状によって一様に励起照射されない可能性がある。また、一様照射の光源であっても、小動物などの複雑な形状を持つ生体試料においては、外形形状よっては励起光強度にかなりのバラツキが生じる。これによって、断層化計算による分子プローブの濃度の算出時に定量性が得られないか、あるいは照射ムラの影響がアーチファクトとして付加されてプローブの位置精度が劣化する、といった可能性がある。

本発明は光拡散方程式を用いた蛍光画像作成の計算精度を維持しつつ計算コストを下げることを目的とする。

本発明は生体試料の形状に起因する問題（Ｉ）と（II）を解決するために、生体試料をメッシュ構造で離散化する際に、条件式に基づいてメッシュサイズを最適化すること、及び励起照射強度の補正を行うことによって、上記目的を達成しようとするものである。

すなわち、本発明は、図1に示されるように、試料を離れた位置から、広がる光束の励起光で広範囲に照射し前記試料中の対象物を励起して蛍光を発生させる少なくとも１つの光源部と、試料からの蛍光を検出する二次元検出器４と、試料を複数の方向から同時に観測するとともに試料から放出される各方向の蛍光を二次元検出器４に導く検出光学系と、試料の三次元表面形状データ及び二次元検出器４が検出した蛍光実測データから試料中の蛍光画像を作成する画像処理装置２と、画像処理装置２が作成した蛍光画像を表示する表示装置６と、を備えた蛍光画像撮像装置である。

画像処理装置２は、図２に示されるように、試料における励起光と蛍光の光伝搬解析を光拡散方程式の有限要素解析によって理論計算する順問題解析部１０と、試料の三次元表面形状データ、二次元検出器４が検出した蛍光実測データ及び順問題解析部１０における解析結果に基づいて試料内部の蛍光画像を作成する逆問題解析部２０を備えている。そして、順問題解析部１０は試料について有限要素解析のための最適メッシュサイズを特定の条件式から導き出すメッシュサイズ導出部１２と、メッシュサイズ導出部１２により導出されたメッシュサイズの各メッシュに対して試料の三次元表面形状データ及び別途求められた励起光の照射分布測定値を用いて励起光強度を補正する励起光強度補正部１４とを備えており、逆問題解析部２０は順問題解析部１０で求められたメッシュサイズと補正された励起光強度を用いて蛍光画像を作成する。

メッシュサイズを導出するための条件式の好ましい一例は、メッシュサイズをmesh_size、光の平均自由行程をλ、形状関数の次数をε、及び比例係数をfactorとしたとき、
mesh_size＝λ×factor×ε
として表わされるものである。

本発明では、試料として生体試料を用いた場合でも最適なメッシュサイズを決定することができ、さらにその最適化されたメッシュに対して励起光照射の強度補正を行なうことができるので、光拡散方程式を用いた蛍光断層化の計算コストや計算精度を最適化することが可能となる。

本発明の主要部を示すブロック図である。 本発明における画像処理装置を示すブロック図である。 本発明の一実施形態における動作全体を示すフローチャートである。 同実施形態における順問題解析部の動作を示すフローチャートである。 光の平均自由行程λを説明する概略図である。 形状関数の次数εを説明するグラフである。 有限要素法による光伝搬解析のための簡単なモデルを示す概略斜視図である。 有限要素法におけるメッシュサイズと計算時間の関係を示す図である。 光拡散方程式における解析解と数値解の関係をメッシュサイズごとに示すグラフである。 光拡散方程式における解析解と数値解の関係を光の平均自由行程λごとに示すグラフである。 光拡散方程式における解析解と数値解の関係を形状関数の次数εごとに示すグラフである。 メッシュサイズ（メッシュ数）と得られる蛍光画像の関係を示す図表である。 拡散反射光と検出器への入射光との関係を示す概略図である。 生体試料とそれから得られる三次元表面形状データの一例を示す斜視図である。 生体試料表面への励起光照射分布を示す概略図である。 円筒媒体を離れた位置から広がる光束で広範囲に照射したときの励起光伝搬と励起光強度補正の関係を示す図である。 円筒媒体を離れた位置から広がる光束で広範囲に照射したときに得られる蛍光画像と励起光強度補正の関係を示す図である。 本発明の一実施形態を示す概略斜視図である。 同実施形態の概略正面断面図である。 概略正面断面図（図１８Ｂ）の導光路を示した図である。 同実施形態においてレンズ５から下をみたときの試料とその像を示す平面図である。 ５方向の像を撮影する実施例を示す概略正面断面図である。 同実施例において撮影される像を示す平面図である。 本発明の一実施形態で使用する光源部の一例を示す概略斜視図である。 同光源部における励起光と励起用干渉フィルタの特性の一例を示すスペクトル図である。

図３，４のフローチャートを参照して、本発明の一実施形態の動作を説明する。図３は一実施形態での処理の概要を示し、図４は同実施形態における具体的な手順を示している。

一実施形態での処理は下記２つの手順からなる。
手順１：条件式に従って、メッシュサイズの最適化を行う。
手順２：光源分布および立体形状の取得データと手順１で得た情報を元に、励起光の強度補正を行う。

（手順１）
本手順では、小動物などの複雑形状をメッシュ構造で有限要素解析する場合のメッシュサイズの最適化を行う。

一実施例形態では、次の条件式に基づいてメッシュサイズを規定する。
mesh_size＝λ×factor×ε・・・（１）
ここで、
λ：平均自由行程、
factor：定数（＝２）、
ε：形状関数の次数、
を表す。

有限要素法を用いた光伝搬解析では、メッシュサイズによって解の精度が決定されることが知られているが、精度よく短時間に計算できるメッシュサイズについては明確な指標がない。本発明の一実施例形態で提案する条件式（１）は、後で実証するように、光伝搬問題に対して簡易にメッシュサイズを提示できるものとなっている。

平均自由行程とは、図５に示されるように、フォトンが一回散乱するのに必要な距離であり、等価散乱係数μｓ’の逆数から、
λ＝１／μｓ’・・・（２）
と表すことができる。等価散乱係数は生体組織や光の波長によって値が異なることが知られている。

（１）式より、平均自由行程が長ければメッシュサイズは大きく設定することができ（図５中の（Ａ）参照。）、逆に平均自由行程が短ければメッシュサイズは小さく設定する必要がある（図５中の（Ｂ），（Ｃ）参照。）。

また、形状関数とは、図６に示されるように、有限要素法において要素内の任意点の変位を表す際に用いる近似関数である。（１）式より、形状関数の次数が高ければメッシュサイズは大きく設定することができ（図６中の（Ａ）参照。）、逆に形状関数の次数が低ければメッシュサイズは小さく設定する必要がある（図６中の（Ｂ）参照。）。

＜メッシュサイズの最適化１＞
ここで、条件式（１）で最適なメッシュサイズを導き出した例を述べる。

散乱係数μｓ＝１０［ｍｍ^-1］、吸収係数μａ＝０．０２［ｍｍ^-1］の光学特性を持った媒体に対して、光伝搬解析を行うことを考える。等価散乱係数μｓ’は、μｓ’＝（１−ｇ）μｓ表わされる。ｇは非等方散乱係数と呼ばれる係数で、散乱の異方性を示しており、位相関数pの余弦平均として下記式で定義される。
p(θ)は散乱位相関数を表す。一般に、ｇ≒０だと等方散乱、ｇ＝０．９〜１だと前方散乱を示す。ここでは、前方散乱のｇ＝０．９とする。

この場合、フォトンが生体媒体を伝搬する際に、フォトンが散乱するまでの平均自由行程λは（２）式を用いて「１ｍｍ」と求められる。

また、有限要素解析における形状関数の次数を２とした場合、（１）式より最適なメッシュサイズは「４ｍｍ」と求められる。

ここで求めた「メッシュサイズ：４ｍｍ」が妥当なものであるかどうか検証するために、有限要素法で求めた数値解を解析解と比較検証した結果について示す。具体的な検証方法は以下のとおりである。

まず始めに、生体内における光の伝搬過程を光拡散方程式によって近似する。
ここで、
Ｄ：拡散定数、
μａ：吸収係数、
φｉ：ｉ番目のメッシュに励起光を与えた場合の位置ｒにおける励起光フルエンスレート、
Ｉｔｉ：ｉ番目のメッシュに与えられた（補正）励起光強度、
を表す。

また、（３）式の境界条件は（４）式となる。
ここで、Ａは媒体表面での拡散光の全反射による反射率から決定される係数を表す。

（３）、（４）式から、解析解を求めると（５）式のようになる。
ここで、Ｒ：反射光強度、
ｒ１：点光源Ｓからの距離、
ｒ２：鏡像点光源（外挿境界を隔てて点光源と対称位置にある光源）からの距離、
を表す（非特許文献２の図５参考。）。

なお、解析解として、図７に示されるような単純な形状（例えば立方体など）において媒体の表面に光源を設置した場合の媒体表面からの反射光強度（Ｒ［ｍＷ／ｍｍ²］）を求めた。

次に、有限要素法を用いて（３）式の数値解（反射光Ｒ）を求めた。この際、図８に示されるように、メッシュサイズを変えてそれぞれ計算を行った。先に示したメッシュサイズの定義より、有限要素法における分割(離散化)の程度によって、メッシュサイズが決定される。ここでのメッシュサイズとは、３次元有限要素法における４面体メッシュの一辺の長さを示す。メッシュサイズを小さく設定すると、物理形状の離散化誤差が小さくなる。このため、解の精度は向上するが、計算時間は増大することが一般的に知られている。逆に、メッシュサイズを大きく設定すると、物理形状の離散化誤差が大きくなる。このため、計算時間は減少するが、解の精度は低下する。図８にはメッシュサイズとともに要した計算時間も示している。

このようにメッシュサイズを変えて光源からの距離（ρ［ｍｍ］）に対する反射光強度Ｒを求めた結果を図９に示す。この方法で求めた解析解（図９中の実線）と、有限要素法によって求めた数値解（図９中の×点を結ぶ線）を比較して、両者の解が誤差何％以内であるかを調べた。図８には、メッシュサイズの条件と１回の有限要素解析の計算時間を示している。なお、本計算の形状関数ではラグランジュ２次を用いた。また、計算には有限要素解析ツール「COMSOL Multiphysics Version3.5」を用いた。計算機として、Dell Precision T7500（CPU:Intel Xeon 3.19GHz 24GB）を使用した。

図８の（サンプル１）に示すように、メッシュ数を９７５０７、メッシュサイズを２ｍｍにした場合は１回の光伝搬解析に約２０秒程度要することが分かる。また、図９の比較結果より、解析解と数値解の誤差が２０％以内となるため、この計算結果の精度は十分であると結論付けられる。

一方、図８の（サンプル３）に示すように、メッシュ数を１１２０１、メッシュサイズを８ｍｍにした場合は１回の光伝搬解析が約０．９ｓｅｃ程度と短くなるが、図９の比較結果では解析解と数値解の誤差が極めて大きくなるため、この計算結果の精度は不十分であると結論付けられる。

このような方法で、計算結果が妥当かどうかの判定を実施し、最適なメッシュサイズを決定する。ここでの最適なメッシュとは、計算精度が十分確保されており、かつ計算時間が短いものを指す。図８、図９での例では、（サンプル２）（＝メッシュサイズ４ｍｍ程度）が最適なメッシュであると判断できる。

＜メッシュサイズの最適化２＞
また、平均自由行程λを（Ａ）０．５ｍｍ、（Ｂ）１．０ｍｍ、（Ｃ）２．０ｍｍと変えた場合に、（１）の条件式で最適なメッシュサイズが算出できるかどうかの検討を行った例を図１０に示す。

（１）式より、最適なメッシュサイズはそれぞれ、
（Ａ）の場合：２ｍｍ、
（Ｂ）の場合：４ｍｍ、
（Ｃ）の場合：８ｍｍ
となる。

一方、それぞれ実際のメッシュサイズは４ｍｍである。実際のメッシュサイズとは、有限要素法で計算した際のメッシュサイズを示す。 (Ａ)、(Ｂ)、(Ｃ)の平均自由行程λの場合について（１）式に従って求めた上の最適なメッシュサイズがそれぞれ２、４、６ｍｍであるのに対し、実際に有限要素法で計算した（データとして用いた）メッシュサイズは４ｍｍである。

（Ａ）（平均自由行程λが０．５ｍｍ）の場合、最適なメッシュサイズが２ｍｍであるのに対し、実際のメッシュサイズは４ｍｍなので条件を満たしていない。図１０（Ａ）の比較結果では解析解と数値解の誤差が極めて大きくなるため、この計算結果の精度は不十分であると結論付けられる。

（Ｂ）（平均自由行程λが１．０ｍｍ）の場合、最適なメッシュサイズが４ｍｍであるのに対し、実際のメッシュサイズも４ｍｍなので条件を満たしている。図１０（Ｂ）の比較結果では解析解と数値解の誤差が２０％以内となるため、この計算結果の精度は十分であると結論付けられる。

（Ｃ）（平均自由行程λが２．０ｍｍ）の場合、最適なメッシュサイズが８ｍｍであるのに対し、実際のメッシュサイズは４ｍｍなので、この場合も条件を満たしている。図１０（Ｃ）の比較結果でも解析解と数値解の誤差が２０％以内となるため、この計算結果の精度は十分であると結論付けられる。

以上の結果より、（１）の条件式で最適なメッシュサイズが算出できることを示した。

＜メッシュサイズの最適化３＞
また、形状関数の次数を（Ａ）２、（Ｂ）３と変えた場合に、（１）の条件式で最適なメッシュサイズが算出できるかどうかの検討を行った例を図１１に示す。

（１）式より、平均自由工程λ＝１ｍｍとして、最適なメッシュサイズはそれぞれ、
（Ａ）の場合：４ｍｍ、
（Ｂ）の場合：６ｍｍ
となる。

一方、それぞれ実際のメッシュサイズは６ｍｍである。

（Ａ）（形状関数の次数が２）の場合、最適なメッシュサイズが４ｍｍであるのに対し、実際のメッシュサイズは６ｍｍなので条件を満たしていない。図１１（Ａ）の比較結果では解析解と数値解の誤差が極めて大きくなるため、この計算結果の精度は不十分であると結論付けられる。

（Ｂ）（形状関数の次数が３）の場合、最適なメッシュサイズが６ｍｍであるのに対し、実際のメッシュサイズは４ｍｍなので条件を満たしている。図１１（Ｂ）の比較結果では解析解と数値解の誤差が２０％以内となるため、この計算結果の精度は十分であると結論付けられる。

＜メッシュサイズを最適化した上での断層化結果＞
メッシュサイズの最適化を行った場合に、断層化計算にどのような効果が期待できるか検討した実施例を図１２に示す。

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

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

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

一方、今回対象とする問題では、光の吸収・散乱によって測定ノイズの影響が強くなるため、（９）式の最小化問題では解が発散する可能性が高い。そこで、ノイズによる発散を防ぐ目的で以下の(１０)式の最小化問題を考える。この最小化問題（チホノフ正則化）は逆問題の研究分野で一般的に利用されており、ベクトルfのノルムを評価関数に組み入れることで解の発散を防ぐことができる。
ここで、
λ：正則化パラメータ、
を表す。

図１２では、直径５０ｍｍ、長さ２５ｍｍの円筒媒体を利用し、メッシュ数を変えた場合の計算時間と断層画像化計算結果を示している。図１２より、メッシュサイズを２ｍｍ（メッシュ数：８０１７３）にした場合（図８のサンプル１と同程度）は１回の光伝搬解析に１９秒程度かかっており、これを断層画像化計算に利用する場合は約１２４時間程度必要とするため、実用的ではない。

なお、順問題解析部でシステム行列を求める際、列ベクトルａをk番目の位置までそれぞれ求める必要がある。1回の光伝搬解析とは1つの列ベクトルを求める作業であり、これをk番目まで求めようとすると、メッシュサイズが２ｍｍの場合は約１２４時間程度必要になる。一方、ここでは「逆問題解析部」の計算時間は問題にしていない。つまり、断層画像化計算に必要な「順問題解析部」と「逆問題解析部」のうち、「順問題解析部」に膨大な時間が必要である。この計算時間がメッシュサイズ（メッシュ数）依存である。

一方、メッシュサイズを４ｍｍ（メッシュ数：２００００）にした場合（図８のサンプル２と同程度）は１回の光伝搬解析が２秒程度となり、断層画像化計算に利用するためには約１２時間程度で可能となる。したがって、実用上はメッシュサイズを４ｍｍ程度とすることが望ましい。

条件式（１）を利用すれば、実際の計算を実施する前に上述の結果（メッシュサイズ：４ｍｍ）を予想することができるため、（１）式による最適メッシュサイズの算出は大変有意義な方法といえる。

（手順２）
手順２では、光源分布及び立体形状の取得データと手順１で得た情報を元に、励起光の強度補正を行う。

第１に、励起光分布の取得を行う。
まず、表面に拡散反射を起こす材質（紙など）を用いた治具に対して励起光を照明し、その反射光データ取得する。この際、励起光強度が高い場合はＮＤフィルタ（減光フィルタ）やレンズの絞りを絞ることで検出器が飽和しないような機構を備える。

次に、光源と治具の位置関係と、得られた反射光強度データとから光源の入射光Ｉｉ（ｘ，ｙ，ｚ）を求める。入射光Ｉｉの具体的な計算手法を図１３に示す。紙などの表面の粗い物体に対して、反射光Ｉは等方拡散となり、物体の法線ベクトルと入射光のなす角度をθとした場合のｃｏｓθおよび拡散反射率Ｋｄに比例する。すなわち、
Ｉ＝Ｋｄ・Ｉｉ・ｃｏｓθ
である。したがって、光源と治具の位置関係が既知であれば観測データＩから入射光強度Ｉｉ（ｘ，ｙ，ｚ）を求めることができる。このような手順で最終的に各光源における励起光強度を求める。

第２に、立体形状の取得を行う。
立体形状データ取得の手法として多視点画像統合法、光切断法及びＸ線ＣＴなどが挙げられる。例として、Ｘ線ＣＴでマウスを撮像し、三次元表面形状データへ変換した結果を図１４に示す。三次元表面形状データは、任意の点で構成される平面データ（メッシュデータ）及び平面の法線ベクトルで構成される。そのメッシュサイズは手順１で決定することができる。

第３に、上述から得られた情報を元に、励起光の強度補正を行う。
具体的には、励起光照明分布のデータＩｉと、生体試料の三次元表面形状データから、実際に生体試料へ与えられる励起光Ｉｔを求める。

Ｉｔを求める計算手法を図１５に示す。入射光とメッシュの法線ベクトルのなす角をθ₁、屈折光とメッシュの法線ベクトルのなす角をθ₂、生体外の屈折率をｎ₁、生体内の屈折率をｎ₂とおくと、生体試料に入る励起光強度は、
と表すことができる。ここでＴは透過率である。ｉ＝１からＮまでのすべてのメッシュに対して同様の処理を行うことによって、実際に生体媒体に与えられる励起光の強度補正処理を実行することができる。

＜メッシュサイズの最適化／励起光強度補正を行った上での断層化結果＞
本実施形態を蛍光断層画像化に適用した例を図１６、１７に示す。

図１６（ａ）は、励起光がすべて一様と仮定した場合、円筒媒体での光伝搬を示したものである。一方、図１６（ｂ）は、本実施形態を利用して励起光の補正を行った場合の、円筒媒体での光伝搬を示したものである。図１６（ａ）と（ｂ）を比較すると、円筒媒体の両端点で形状による光伝搬の誤差が生じていることが分かる。つまり、照射方向と媒体の法線方向のなす角度が比較的大きな箇所で、補正の効果が顕著に表れることが分かる。また、光源分布を非均一にした場合でも、分布に対応して光伝搬の誤差が生じることを確認している。

本実施形態を蛍光断層画像化に組込んだ場合のシミュレーション結果を図１7に示す。シミュレーションでは、直径２５ｍｍ、長さ２５ｍｍの円筒媒体に、直径１ｍｍ、長さ４ｍｍの蛍光剤を埋め込んで得られた蛍光観察データを用いて計算を実施した。データには光ショットノイズ（光の粒子性に基づく雑音）を加えた。

また、断層画像化計算には一般的な逆問題手法であるチホノフ正則化を用いた。計算結果が逆問題手法に依存しないように、すべてのパラメータを固定して計算した。

図１７（ａ）は励起光の補正を行わずに断層画像化計算を行った結果であり、蛍光剤が元の位置に推定されていないことが分かる。図１７（ｂ）は本実施形態を適用して断層画像化計算を行った結果であり、蛍光剤が元の位置に推定されることが分かる。この比較結果より、データに付加されたノイズは、光ショットノイズよりも励起照射ムラによるノイズが支配的となっており、本発明（励起光照射の補正）によってその影響が効果的に緩和されることが分かった。

光学系の一実施形態を図１８Ａ−１８Ｄに示す。図１８Ａは斜視図、図１８Ｂは正面断面図、図１８Ｃはこの場合の導光路を示したものである。生体試料Ｓが載置される試料ホルダー２０と、試料ホルダー２０上の試料Ｓから放出される光の像を検出する二次元検出器４と、試料ホルダー２０上の試料Ｓを複数の方向から観測するとともに試料Ｓから放出される各方向の光の画像を二次元検出器４の方向に導く方向別の導光路と、二次元検出器４と導光路の間に配置され導光路により導かれた複数の画像を二次元検出器４上の方向毎の異なる場所に結像する主結像レンズ５とを備えている。

二次元検出器４は例えばＣＣＤである。以下、二次元検出器４を代表してＣＣＤ４と表示する。ＣＣＤ４には画像処理装置２が接続され、画像処理装置２には画像処理装置２が作成した蛍光画像を表示する液晶表示装置（ＬＣＤ）などの表示装置６が接続されている。画像処理装置２はコンピュータであり、この蛍光画像撮像装置の専用コンピュータ又は汎用のパーソナルコンピュータである。

この実施例では導光路に反射鏡Ｍ１，Ｍ２を含んでいるが、反射鏡、レンズ又は光ガイドなどの光学素子を含むものだけでなく、そのような光学素子を含まないものも総称して導光路と呼ぶ。

導光路として試料Ｓからの直接光を受けない導光路を含み、試料Ｓからの直接光を受けない導光路は試料Ｓの像Ｓ１，Ｓ２を主結像レンズ５の実質的合焦点範囲内に形成するとともに像Ｓ１，Ｓ２を形成した後の光線が主結像レンズ５の方向に進むように配置された光学素子Ｍ１，Ｍ２を含み、かつ少なくとも１つの導光路の光学素子Ｍ１，Ｍ２は像として実像を形成する光学素子である。

これにより、主結像レンズ５が導光路を経由した複数方向の画像を一括して二次元検出器４上に結像するようにしている。

導光路の光学素子として平面鏡だけを用いると収差が全くない利点はあるものの、試料Ｓの下部に虚像しかできず、かつ試料Ｓの裏面方向については平面鏡の虚像が遠方にできて、主結像レンズ５は異なる方向の像の総てに焦点を合わせることができない。またそれを緩和するために補助レンズを使うにしても結像レンズ５の口径に限界があるので、この実施例では以下の対策を講じた。

すなわち、この実施例では、平面鏡の像が遠くにできるような導光路に対して、平面鏡に代えて凹面鏡のような実像を形成する光学素子Ｍ１，Ｍ２を用い、虚像ではなく試料Ｓの実像を、主結像レンズ５から試料Ｓとほぼ等しい距離に作る。

ここでは実像を形成する光学素子として凹面鏡を取り上げて説明するが、凹面鏡に限定されるものではなく、平面鏡と凸レンズの組合せのような光学素子でも実現することができる。図１８Ａ−Ｃの方式では裏側像をＣＣＤ４側に導くのに大きな凹面鏡Ｍ１、Ｍ２を使い、試料Ｓの倒立実像Ｓ１、Ｓ２を試料Ｓの両隣りに作り、これら３つの像Ｓ、Ｓ１、Ｓ２を上から主結像レンズ５を通してＣＣＤ４で撮影する。試料Ｓと２つの像Ｓ１、Ｓ２の三者は主結像レンズからほぼ等距離にあるように凹面鏡Ｍ１、Ｍ２を設定する。

図１８Ａ、１８Ｂにおける像形成の重要なポイントは以下の２点である。
１）凹面鏡Ｍ１、Ｍ２のそれぞれの曲率中心の位置Ｃ１、Ｃ２が重要であって、試料Ｓとその像Ｓ１、Ｓ２の各対応点がＣ１又はＣ２に対して立体的な点対称になっていることである。結像レンズ５から例えばＱ点に焦点を合わせて見ているとし、結像レンズ５から試料Ｓの方に逆向きに進む主光線（結像レンズ５の中心を通る光）を考える。この主光線はＱ点の延長上のＧ点で凹面鏡Ｍ１にあたり、そこで反射され試料Ｓ上のＰ点に向かうが、反射の入射角と射出角が等しいことから、Ｐ点とＱ点の中点にミラーＭ１の曲率中心Ｃ１が存在する。すなわち、Ｇ点から進む光は、対称中心Ｃ１に対し、Ｑ点と点対称のＰ点に向かって反射することになる。試料Ｓとしてマウスを測定する場合を例として挙げると、Ｐ点が図１８Ａに示すような、マウスの頭部の下面であるとすると、Ｑ点はマウスの裏返った試料像Ｓ１の頭部の下部であって腹部がレンズの側を向いている。こうして主結像レンズ５から下を見ると、図１８Ｄのようなマウスの背中が写った正面像とその両隣にマウスの腹部が写った裏面像が観測され、画像においてもＰ点とＱ点がＣ１に対して点対称になっている。同様に試料Ｓの他の点、例えばＨ点は、試料像Ｓ２では、Ｃ２の反対側の点Ｉに写像されている。

なお、Ｃ１点、Ｃ２点の位置を試料Ｓに近づければ、像Ｓ１、Ｓ２が全体として試料Ｓの近くに寄るし、Ｃ１点、Ｃ２点の位置を試料Ｓから離せば像も離れる。Ｃ１点、Ｃ２点の位置を変化させることで像の位置を調節できるから、カメラの視野にバランスよく多方向の像を並べることができる。

２）ポイントの２つ目は実像を形成する光学素子の結像特性に関する。凹面鏡Ｍ１、Ｍ２のような大きな光学素子で試料の実像を作るときは、一般に収差が非常に大きくなる。例えば図１８Ａにおいて、試料像Ｓ１の位置に「紙」をおいて、試料像Ｓ１の結像を見るなら、凹面鏡Ｍ１の全面を使うため収差が大きくボケた像にしかならない。しかし実際には、主結像レンズ５を通して試料像Ｓ１を見るので、例えばＱ点に焦点を合わせるとき、主結像レンズ５に入射する光はＱ点の延長が凹面鏡Ｍ１に交わるＧ点を中心とするδＳで表した面積の部分の光線だけである。すなわち、結像に使われるのは、大きな凹面鏡Ｍ１の一部のδＳで表した小さい凹面鏡でしかない。したがって結像に使われる立体角ωは比較的小さい値（Ｆ値で表すとき、Ｆ６程度以上）となるから、このような大きな凹面鏡でも実用になる収差内に収まる。観測したいＱ点を動かすとき対応する凹面鏡Ｍ１の部分δＳは、Ｑ点につれて凹面鏡Ｍ１の全体を動くけれども、それぞれの場所で結像に使われるのは、上述の面積δＳの範囲に限られる。

さらに好都合な点は、球面の凹面鏡により等倍率の１対１の像を作るときは収差が少なくなることである。殊にＰ点とＱ点が近ければ、（角ＱＧＰが小さいとき）さらに綺麗な結像になる。極端な例は角ＱＧＰがゼロに近ければ、球面凹面鏡の曲率中心から出て同じ点に戻る光に極めて近い状態になって、収差が全くなくなることからも、球面凹面鏡による１対１の結像が収差の点で有利なことが理解できる。

さて、図１８Ａに見られるように像Ｓ１、Ｓ２が１対１の像ならば、試料Ｓを含め同じ大きさに３つの像を結像レンズ５がＣＣＤ４上に再度結像するのでサイズ上のバランスが良い上に、前述した１対１結像が収差的にも有利であるという幸運が重なり、球面凹面鏡を用いる光学素子を有利な方法としている。すなわち、非球面（楕円面など）の凹面鏡でなく、安価な球面凹面鏡で十分な効果が有ることを強調する。さらに特許文献３の「作用」の項には、この１対１の結像の収差及びβ（角ＱＧＰ）が少し大きくなるとどうなるかを、もう少し詳しい計算でも検討し、光学素子において凹面鏡の配置をどのようにすれば、凹面鏡の特性が生きるかのガイドラインについて述べられている。

以上述べたように、この実施例は裏面方向の像をＣＣＤの側に戻すための結像形光学素子を利用している。ここで、図の凹面鏡Ｍ１、Ｍ２が大きな凹面鏡であることも意味があることを追記する。主結像レンズ５から試料の像Ｓ１の１点、例えばＱを見るとき、向こう側の凹面鏡Ｍ１上の点Ｇが存在するが、Ｇが存在できる程度に凹面鏡Ｍ１のサイズが大きい必要がある。試料の像Ｓ１全体を見る必要があるとき、凹面鏡Ｍ１の全体のサイズは図の例では像Ｓ１のサイズの１．５倍程度である。レンズ５から見て試料の像Ｓ１の向こう側に凹面鏡が存在するだけのサイズを要するからである。試料の像Ｓ１の全体でなく、その一部（例えば頭部だけ）の観察でよいなら、凹面鏡Ｍ１、Ｍ２はもっと小さくてもよく、測定目的によってその大きさを選ぶことができる。

［５方向の同時観測法用の光学系］
５つの方向から同時観測する場合を例として、図１９Ａと図１９Ｂを用いて説明する。図１９Ａにおいて、中央に配置した試料Ｓに対して、真上（観測角度０度方向）に主結像レンズ５及びＣＣＤ４が配置され、４つの凹面鏡Ｍ３、Ｍ１、Ｍ２、Ｍ４がそれぞれ観測角度６０°、１３５°、２０５°、３００°の方向に配置されている。すなわち０度方向の導光路は試料Ｓから直接にレンズ５に進む光が占める領域であり、６０°方向の導光路は試料Ｓ→凹面鏡Ｍ３→平面鏡Ｍ５→レンズ５の順に進む光の占める領域であり、同様に１３５°方向は試料Ｓ→凹面鏡Ｍ１→レンズ５、２０５°方向は試料Ｓ→凹面鏡Ｍ２→レンズ５、３００°方向は試料Ｓ→凹面鏡Ｍ４→平面鏡Ｍ６→レンズ５にそれぞれ進む光の占める領域である。このうち１３５°、２０５°方向の導光路中、凹面鏡Ｍ１、Ｍ２による試料Ｓの実像がＳ１、Ｓ２の位置にできており、一旦結像した後主結像レンズ５に向かって進み、再びＣＣＤ４上で結像する。また、６０°方向及び３００°方向の導光路中では、斜め上の凹面鏡Ｍ３とＭ４の実像がＳ３、Ｓ４にできているが、これは折り曲げ鏡Ｍ５とＭ６の効果によって実像から進む光の方向と実像の場所を必要な条件に変換した結果による。

ここで、試料Ｓの下部の注目点Ｐから出る光は凹面鏡Ｍ１により実像Ｓ１上の点Ｑに一度結像し、主レンズ５により、二次元検出器４上にもう一度結像している。また、試料Ｓ上の上部の注目点Ｈから出る光は凹面鏡Ｍ３と折り曲げ鏡Ｍ５で反射し実像Ｓ３上のＩ点に一度結像し、さらにレンズ５によってもう一度二次元検出器４上に結像していることがわかる。

これら４つの実像Ｓ１−Ｓ４と試料Ｓは、レンズ５が焦点を結ぶことのできる同等の距離にあるから、試料Ｓと４つの像の合計５つの像を同時にＣＣＤ４上に結像し撮影することができる。この結果ＣＣＤ４上の試料の５方向像は、図１９Ｂのようになり、左からＳ３（観測方向６０度）、Ｓ１（観測方向１３５度）、Ｓ（観測方向０度）、Ｓ２（観測方向２０５度）及びＳ４（観測方向３００度）である。中央の像だけマウスの頭の像の向きが反対になっている。

この５方向測定光学系は、試料とカメラ（レンズ５とＣＣＤ４）の間に焦点調節用の補助レンズなどの光学部品を挿入する必要が全くないことが特徴で、５つの導光路が主結像レンズ５の近くで互いに重なっていても何ら問題なく、口径が大きい主結像レンズを用いることができる。可動部分もなく、多方向の観測像を一度に共通のＣＣＤ４上に結像することができる明るい観測系が実現できている。

照明用光源部ＩＬ１、ＩＬ２、ＩＬ３、ＩＬ４及びＩＬ５（ＩＬはilluminationを表す）が５箇所備えてあり、各光源部の中に含まれる蛍光励起光源が試料Ｓを照射し、試料Ｓから発生する蛍光を５方向の蛍光像としてＣＣＤ４に結像させて撮影する。なお主結像レンズ５の前に蛍光用（励起光除去）フィルタＦ_EMを挿入できるようになっており、蛍光測定時は励起波長成分をこの励起光除去フィルタＦ_EMで阻止することで、蛍光だけによる画像をＣＣＤ上に作るようになっている。

照明用光源部ＩＬ１、ＩＬ２、ＩＬ３、ＩＬ４及びＩＬ５の具体的な構成例を図２０Ａに示した。図２０Ａの例では照明用光源部は３つの異なる種類の光源Ａ、Ｂ及びＣから成り立っている。光源Ａ、Ｂ及びＣは試料を離れた位置から広がる光束で広範囲に照射することができる。光源Ａは後述するように、生体試料４を照明して外観画像を撮像するための白色ＬＥＤである。撮像された外観画像から本発明における順問題解析部１０及び逆問題解析部２０で使用する三次元表面形状データを得ることができる。光源Ｂ及びＣは蛍光励起光源であり、光源Ｂは波長λ１の励起光を発光する発光素子レーザダイオード（以降ＬＤと省略）ＬＤλ１と、発光素子ＬＤλ１の発光側に取り付けられた励起側干渉フィルタＦ_exλ１からなる。光源Ｃは波長λ２の励起光を発光するＬＤからなる発光素子ＬＤλ２と、発光素子ＬＤλ２の発光側に取り付けられた励起側干渉フィルタＦ_exλ２からなる。ただし、光源Ｂ及びＣは、必要に応じて干渉フィルタＦ_exλ１、Ｆ_exλ２の後に生体試料全体に光を照射するための発散レンズ（図示は省略）を追加することもできる。ここで、図２０Ａに示す３種類の光源Ａ、Ｂ及びＣの選択は、機械的な切替えなしで単に点灯の電源のオン・オフによって行うことができる。したがって、各照射方向ＩＬ１、ＩＬ２、ＩＬ３、ＩＬ４及びＩＬ５に配置されたそれぞれの光源Ａ、Ｂ及びＣの点灯のオン・オフによって、機械的な切替えなしで光源の方向ＩＬ１、ＩＬ２、ＩＬ３、ＩＬ４及びＩＬ５と光源の種類（Ａ、Ｂ及びＣ）の指定を行う。

次に、励起光源と励起用干渉フィルタの特性について図２０Ｂを用いて説明する。図２０Ｂは光源Ｂの発光素子ＬＤλ１とその発光側に取り付けられた干渉フィルタＦ_exλ１との関係を示したものである。発光素子ＬＤλ１は単一波長λ１の光だけを発光すると思われがちであるが、実際には弱い裾の発光を伴なっている。その裾の部分には二次元検出器ＣＣＤ４の入射側に設けられた励起光除去フィルタである蛍光側フィルタＦ_EMを透過する波長成分が含まれており、そのような波長の光が生体試料に向けて照射されると、生体試料で発生した蛍光成分とともに蛍光側フィルタＦ_EMを透過して二次元検出器４に入射する。その結果、撮影した蛍光撮像画像に試料の蛍光以外の漏れ光成分による迷光が重なり蛍光の検出感度を下げてしまう。

そこでこの実施例では、発光素子ＬＤλ１の発光側に発光波長の裾の部分の光を遮光する干渉フィルタＦ_exλ１を取り付けて、蛍光側フィルタＦ_EMの透過域にある発光素子ＬＤλ１からの励起光成分を除去している。これにより、蛍光側フィルタＦ_EMを透過する光は生体試料からの蛍光成分のみとなり、二次元検出器４上の撮像画像に迷光が混入にすることによる蛍光検出能力の低下を防止できる。

なお、上記では光源Ｂについて説明したが、光源Ｃについても光源Ｂと同様の構成である。

照明用光源部ＩＬ１、ＩＬ２、ＩＬ３、ＩＬ４及びＩＬ５に含まれる蛍光励起光源として、ＬＤやＬＥＤを使うなら電気回路のオン・オフにより必要な光源だけを自由に点灯消灯できるので、５方向の光源を総て点灯して、蛍光を測定することができるほか、必要な組合わせだけを選択的に点灯させて蛍光画像測定を行うこともできる。こうすると、試料Ｓから見るとき、前から励起光を照らすときの蛍光画像だけでなく、後ろからだけとか、横からだけ照らすときの画像が得られるから、試料Ｓからの５方向の画像それぞれに対して励起方向の異なる５つの画像を得るから合計２５枚の画像を５回の露光で得ることができる。２５枚の画像から動物の体内の浅い位置に発光源があるか、深い位置に発光原があるかを推定できる。すなわち、浅い位置の発光源なら、２５枚のどれかの画像中の被写体の小さい部分が強く光ることが推測されるのに対し、深い発光源ならば２５枚のどの画像にも拡散した発光分布になるからである。

こうして蛍光励起方法は可動部分が無くて単に励起光の点滅だけで、試料の前横後ろから励起する、励起方法を自在に設定できることである。こうして蛍光モードの場合でも、試料の全周にわたる多方向から励起・観測画像を簡単に得られることを示した。

上記において蛍光励起光は試料から「離れて」照射される。

図２０Ａに示した光源のうち、符号Ａで示した光源は、白色ＬＥＤであって試料の外観像を撮影するときに点灯させる。試料のどの部分が光っているかを知るために蛍光測定の前又は後に白色ＬＥＤ照明の画像を取得しておく。画像解析時に三次元表面形状データとして使用する他、外観照明のデータを蛍光画像に重ねた図で比較するとか、外観照明のデータと蛍光画像を並べて比較することで、該当の場所が外観図のどの部位に当たるかを容易に知ることができる。

光学系についての別の実施例は特許文献３の図６、図７、図８、図９、図１０、図１１及び図１４に記載されている。本発明の蛍光画像撮像装置の光学系としてはそれらの光学系を適用することもできる。

２ 画像処理装置
４ 二次元検出器
６ 表示装置
１０ 順問題解析部
１２ メッシュサイズ導出部
１４ 励起光強度補正部
２０ 逆問題解析部

Claims (2)

1. 試料を離れた位置から広がる光束の励起光で広範囲に照射し前記試料中の対象物を励起して蛍光を発生させる少なくとも１つの光源部と、前記試料からの蛍光を検出する二次元検出器と、前記試料を複数の方向から同時に観測するとともに前記試料から放出される各方向の蛍光を前記二次元検出器に導く検出光学系と、前記試料の三次元表面形状データ及び前記二次元検出器が検出した蛍光実測データから前記試料中の蛍光画像を作成する画像処理装置と、前記画像処理装置が作成した蛍光画像を表示する表示装置と、を備えた蛍光画像撮像装置において、
   前記画像処理装置は、前記試料における励起光と蛍光の光伝搬解析を光拡散方程式の有限要素解析によって理論計算する順問題解析部、並びに
   前記三次元表面形状データ、前記二次元検出器が検出した蛍光実測データ及び前記順問題解析部における解析結果に基づいて前記試料内部の蛍光画像を作成する逆問題解析部を備えており、
   前記順問題解析部は、前記試料について有限要素解析のための最適メッシュサイズを特定の条件式から導き出すメッシュサイズ導出部と、前記メッシュサイズ導出部により導出されたメッシュサイズの各メッシュに対して前記三次元表面形状データ及び別途求められた前記励起光の照射分布測定値を用いて励起光強度を補正する励起光強度補正部と、を備えており、
   前記逆問題解析部は前記順問題解析部で求められたメッシュサイズと補正された励起光強度を用いて前記蛍光画像を作成するものであることを特徴とする蛍光画像撮像装置。
2. 前記条件式は、メッシュサイズをmesh_size、光の平均自由行程をλ、形状関数の次数をε、及び比例係数をfactorとしたとき、
   mesh_size＝λ×factor×ε
   として表わされるものである請求項１に記載の蛍光画像撮像装置。