JP2018151334A - 計測システム - Google Patents

Abstract

【課題】膜電位分布を従来よりも高速に計測し得る計測システムを提供する。
【解決手段】計測システム１では、パターン光Ｌ１を組織表面６ａに照射することで、蛍光を多点で同時計測できる蛍光分布画像を得ることができ、１点１点蛍光を計測していた従来の計測システムに比べて、膜電位分布Ｖを高速に計測し得る。また、計測システム１では、非ゼロ要素の数が少ないスパース行列αを再構成し、逆変換によってスパースでない膜電位分布Ｖを求めるようにしたことで、圧縮センシングの理論の下、少ない枚数の蛍光分布画像からでも膜電位分布Ｖを求めることができる。よって、計測システム１では、蛍光分布画像の枚数を減らせる分、処理負担が軽減し、膜電位分布を高速に計測し得る。
【選択図】図１

Description

本発明は、計測システムに関するものである。

近年、メゾスコピックな領域での拡散光トモグラフィを試みている研究が報告されている（例えば、非特許文献１〜３参照）。この手法は、Ｌａｍｉｎａｒ Ｏｐｔｉｃａｌ Ｔｏｍｏｇｒａｐｈｙ（ＬＯＴ）と呼ばれている。この種の計測システムでは、例えば１本の励起レーザ光と、２つのガルバノミラーと、を用いて組織表面上をラスタースキャンする。組織表面で励起された蛍光はガルバノミラーを通り、ダイクロイックミラーによって反射され、例えば７個の光電子増倍管で検出される。レーザ照射点と光電子増倍管との位置関係は、予め固定されており、レーザ照射点から近い点から遠い点まで蛍光を計測している。そして、用意した生体組織（以下、単に組織と呼ぶ）の光拡散モデルを用いて、計測された蛍光のデータから３次元情報を再構成している。

Hillman, E. M., Devor, A., Bouchard, M. B., Dunn, A. K., Krauss, G. W., Skoch, J., ... & Boas, D. A. (2007). Depth-resolved optical imaging and microscopy of vascular compartment dynamics during somatosensory stimulation. Neuroimage, 35(1), 89-104. Hillman, E. M., Dunn, A. K., & Boas, D. A. (2004, April). Laminar optical tomography of rat cortical activation: resolving depth-dependent hemodynamics from 0 to 2mm. In Biomedical Topical Meeting (p. FC5). Optical Society of America. Hillman, E. M., Bernus, O., Pease, E., Bouchard, M. B., & Pertsov, A. (2007). Depth-resolved optical imaging of transmural electrical propagation in perfused heart. Optics express, 15(26), 17827-17841.

しかしながら、かかる計測システムでは、１本の励起レーザ光と、複数個の光電子増倍管（ディテクタ）との位置関係を固定した状態で２枚のガルバノミラーによって１点１点、組織表面上をラスタースキャン計測してゆくことから、計測速度はガルバノミラーの動作速度に依存する。励起レーザ光によるスキャンスピードは、組織表面とガルバノミラーとの距離を開けることで、速めることが可能であるが、励起レーザ光の広がりや、組織表面での蛍光輝度値がさらに微弱になるという問題点があるため、膜電位分布を高速に計測し難いという問題があった。

本発明は、上記のような問題に鑑みてなされたものであり、膜電位分布を従来よりも高速に計測し得る計測システムを供することを目的とする。

本発明の計測システムは、光源から照射された励起光を基に、所定の照射点分布パターンのパターン光を生成し、前記パターン光を組織表面に照射する光学機器と、前記組織表面における蛍光分布を撮像して該組織表面の蛍光分布画像を取得する撮像機器と、前記蛍光分布画像に基づいて組織の膜電位分布Ｖを求める演算処理装置と、を備え、前記演算処理装置は、前記励起光および蛍光の強度分布関係をシミュレーションから求めた計測感度行列Ｊと、前記計測感度行列Ｊをスパース行列αにするための空間変換行列Ψと、を記憶する記憶部を備える。また、前記演算処理装置は、前記計測感度行列Ｊおよび前記空間変換行列Ψを基底行列として扱い、前記蛍光分布画像から求めた蛍光輝度値から前記スパース行列αを再構成する再構成演算部と、前記再構成演算部で再構成した前記スパース行列αを、前記空間変換行列Ψの逆変換によって前記膜電位分布Ｖに変換する膜電位分布演算部と、を備える。

本発明によれば、照射点分布パターンをもつパターン光を組織表面に照射することで、蛍光を多点で同時計測できる蛍光分布画像を得ることができ、１点１点蛍光を計測していた従来の計測システムに比べて、膜電位分布を高速に計測し得る。また、本発明では、非ゼロ要素の数が少ないスパース行列を再構成し、逆変換によってスパースでない膜電位分布を求めるようにしたことで、少ない枚数の蛍光分布画像からでも膜電位分布を求めることができ、その分、処理負担が軽減し得、膜電位分布を高速に計測し得る。

本発明による計測システムの全体構成を示す概略図である。 ＤＭＤのマイクロミラーのオンオフ動作を示す概略図である。 組織表面に照射されたパターン光の照射点分布パターンの一例を示す概略図である。 ＬＯＴにおける計測原理の説明に供する概略図である。 図５Ａは組織内における励起光強度分布Ｈの説明に供する概略図であり、図５Ｂは計測感度行列Ｊの説明に供する概略図である。 圧縮センシングの説明に供する概略図である。 圧縮センシング適用時のデータ構造を示す概略図である。 不応期と興奮領域の説明に供する概略図である。 ＤＷＴにおけるマルチレベル分解構造の説明に供する概略図である。 ３次元データにおけるＤＷＴマルチレベル分解構造の説明に供する概略図である。 演算処理装置の回路構成を示すブロック図である。 膜電位分布計測手順を示すフローチャートである。 画像取得処理手順を示すフローチャートである。 図１４Ａは、上から下に興奮伝播している３次元の膜電位分布モデル「Ｖｅｒｔｉｃａｌ１」を示す概略図であり、図１４Ｂは、下から上に興奮伝播している３次元の膜電位分布モデル「Ｖｅｒｔｉｃａｌ２」を示す概略図であり、図１４Ｃは、鉛直波面で興奮伝播している３次元の膜電位分布モデル「Ｈｏｒｉｚｏｎｔａｌ１」を示す概略図であり、図１４Ｄは、＋４５度に傾いた波面で興奮伝播している３次元の膜電位分布モデル「Ｈｏｒｉｚｏｎｔａｌ２」を示す概略図であり、図１４Ｅは、−４５度に傾いた波面で興奮伝播している３次元の膜電位分布モデル「Ｈｏｒｉｚｏｎｔａｌ３」を示す概略図であり、図１４Ｆは、中心から球面波で興奮伝播している３次元の膜電位分布モデル「Ｓｐｈｅｒｅ」を示す概略図である。 図１５Ａは、ＤＷＴによる非ゼロ要素数の平均（カッコ内は標準偏差）をまとめた表であり、図１５Ｂは、上位５０個の係数の累積エネルギー平均（カッコ内は標準偏差）をまとめた表である。 シミュレーション実験におけるモンテカルロシミュレーションのパラメータをまとめた表である。 照射点分布パターンの説明に供する概略図である。 各ボクセルの２値化処理に用いる閾値の説明に供するグラフである。 図１９Ａは、ＤＣＴの計測回数および照射点数を変えたときの再構成誤差ボクセル数を示すグラフであり、図１９Ｂは、Ｈａａｒ １．３の計測回数および照射点数を変えたときの再構成誤差ボクセル数を示すグラフであり、図１９Ｃは、Ｈａａｒ ３．３の計測回数および照射点数を変えたときの再構成誤差ボクセル数を示すグラフである。 図２０Ａは、膜電位分布モデル「Ｖｅｒｔｉｃａｌ１」での再構成誤差ボクセル数を示すグラフであり、図２０Ｂは、膜電位分布モデル「Ｖｅｒｔｉｃａｌ２」での再構成誤差ボクセル数を示すグラフであり、図２０Ｃは、膜電位分布モデル「Ｈｏｒｉｚｏｎｔａｌ１」での再構成誤差ボクセル数を示すグラフであり、図２０Ｄは、膜電位分布モデル「Ｈｏｒｉｚｏｎｔａｌ２」での再構成誤差ボクセル数を示すグラフであり、図２０Ｅは、膜電位分布モデル「Ｈｏｒｉｚｏｎｔａｌ３」での再構成誤差ボクセル数を示すグラフであり、図２０Ｆは、膜電位分布モデル「Ｓｐｈｅｒｅ」での再構成誤差ボクセル数を示すグラフである。 図２１Ａは、図１４Ａの膜電位分布モデル「Ｖｅｒｔｉｃａｌ１」の再構成結果を示す概略図であり、図２１Ｂは、図１４Ｂの３次元膜電位分布モデル「Ｖｅｒｔｉｃａｌ２」の再構成結果を示す概略図であり、図２１Ｃは、図１４Ｃの３次元膜電位分布モデル「Ｈｏｒｉｚｏｎｔａｌ１」の再構成結果を示す概略図であり、図２１Ｄは、図１４Ｄの３次元膜電位分布モデル「Ｈｏｒｉｚｏｎｔａｌ２」の再構成結果を示す概略図であり、図２１Ｅは、図１４Ｅの３次元膜電位分布モデル「Ｈｏｒｉｚｏｎｔａｌ３」の再構成結果を示す概略図であり、図２１Ｆは、図１４Ｆの３次元膜電位分布モデル「Ｓｐｈｅｒｅ」の再構成結果を示す概略図である。

（１）本発明の計測システムについて
図１は、本発明の計測システムの全体構成を示す概略図である。計測システム１は、光学機器２と撮像機器３と演算処理装置４とで構成されている。計測システム１は、光学機器２で生成されたパターン光Ｌ１を、計測対象となる組織６に照射し、組織表面６ａにおける蛍光分布を撮像機器３により撮像する。演算処理装置４は、撮像機器３で取得した蛍光分布画像を基に膜電位分布を求めることができる。

光学機器２は、光源９から出射した励起光をエキスパンダー１０で拡大し、デジタルマイクロデバイス（Digital Micro-mirror Device：ＤＭＤ）１２に入射する。ＤＭＤ１２は、照射点分布パターンがランダムに変わるパターン光Ｌ１を生成し、当該パターン光を対物レンズ１３およびダイクロイックミラー１４を介して組織表面６ａに照射する。対物レンズ１３は、組織表面６ａに対してパターン光Ｌ１を拡大および縮小し、当該パターン光Ｌ１の照射点分布パターンを組織表面６ａの最適な領域に投影させる。

組織表面６ａに照射されたパターン光Ｌ１により励起され、組織表面６ａ上に現れた色素の蛍光Ｌ２は、ダイクロイックミラー１４で分離され、イメージングレンズ１５を介して撮像機器３により撮像される。撮像機器３は、照射点分布パターン毎に得られる組織表面６ａの蛍光分布を撮像してゆき、取得した複数の蛍光分布画像を演算処理装置４に送出する。演算処理装置４は、照射点分布パターン毎に得られた複数の蛍光分布画像を用いて、圧縮センシングの理論を利用した演算処理を実行し、組織６における膜電位分布を求める。

ここで、ＤＭＤ１２は、図２に示すように、アレイ状に配置された複数のマイクロミラー１６を有しており、各マイクロミラー１６をそれぞれ独立にオンオフ動作することにより所定のマイクロミラー１６だけを傾かせる。オン動作したマイクロミラー１６は、励起光Ｌ３を反射し、対物レンズ１３を透過させて組織表面６ａに励起光Ｌ３を照射し得る。ＤＭＤ１２は、ミラーによる光反射のため光効率が高く、かつ各マイクロミラー１６間の隙間も狭いため高い光利用効率を有する。

また、ＤＭＤ１２は、所望の動作周波数に変更でき、各マイクロミラー１６を同期させてオンオフ動作することで、大きな１つのミラーとして扱え、空間解像度についても自由に変更できる。この実施形態の場合、動作周波数が９．５［ｋＨｚ］以下、ピクセルサイズが２５６０×１６００ピクセルでなるＤＭＤ１２を用いることができる。図３に示すように、組織表面６ａには、ＤＭＤ１２によって、励起光Ｌ３の照射点Ｐ１がランダムに選定された、所定の照射点分布パターンのパターン光Ｌ１が照射され得る。

（２）ＬＯＴにおける計測原理について
本発明における、蛍光の多点同時計測の計測システム１について詳細に説明する前に、先ずは、基本的原理を理解するために、励起光Ｌ３を組織表面６ａに照射することで組織６の膜電位分布を求める計測原理の概要について、簡単に説明する。ここでは、図４に示すように、組織６の所定位置を計測座標系の原点を０（図示せず）とし、励起光Ｌ３が照射された照射点Ｐ１の励起光照射位置をｒ_ｓとして説明する。励起光照射位置ｒ_ｓに照射された励起光Ｌ３は、組織６内を拡散し、例えば位置ｒにある色素Ｐ３を励起することで当該位置ｒにて蛍光Ｌ２が発生する。蛍光Ｌ２は、組織６内を拡散伝播し、組織表面６ａ上の蛍光計測位置ｒ_ｄから外部に出射する。このようなとき組織表面６ａ上から計測される静止状態の蛍光Ｆ_ｘ,ｍは次のような式で記述することができる。なお、蛍光Ｆ_ｘ,ｍにおいて、ｘは励起光に依存する関数であることを示し、ｍは蛍光に依存する関数であることを示し、ｔは時間を示す。

ここで、Ｈ_ｘは、組織６内の励起光の強度分布（励起光強度分布とも呼ぶ）を示しており（ｘは励起光に依存する関数であることを示す）、より具体的には励起光波長における組織６の吸収係数、散乱係数、異方性パラメータ、屈折率の関数である。Ｅ_ｍは、組織６内の蛍光の強度分布（蛍光強度分布とも呼ぶ）を示しており（ｍは蛍光に依存する関数であることを示す）、より具体的にはＨ_ｘと同様に蛍光波長における組織６の吸収係数、散乱係数、異方性パラメータ、屈折率の関数である。数２におけるε_ｘは、色素の吸光係数を示しており、励起光の波長の関数である。ｃは色素濃度を示す。η_ｍは、色素の蛍光における量子効率を示しており、蛍光波長による関数である。時間依存の関数ｐ_ｘ,ｍは膜電位の変化に依存して変化するものとなっている。なお、時間依存の関数ｐ_ｘ,ｍにおいて、ｘは励起光に依存する関数であることを示し、ｍは蛍光に依存する関数であることを示す。

ここでは、組織６の色素濃度ｃは時間的に変化せず、色素の励起光および蛍光波長スペクトルは膜電位分布Ｖに応じて変化する、と仮定する。さらに、先行研究の知見から、吸光係数ε_ｘの変化に比べて量子効率η_ｍの変化の方が非常に大きいものとし、次式のように近似を行う。

さらに、ｐ_ｘ,ｍ（ｒ，ｔ）∝Ｖ（ｒ，ｔ）であるため、ｐ_ｘ,ｍ（ｒ，ｔ）を次式のように近似する。ｗはスケーリング定数を示す。

ここで、蛍光のスペクトルシフトは実験的に計測された蛍光から組織６の光学特性変化を与えるには十分でないほどの変化しかないと仮定する。これによりスペクトルシフトによる光学定数の変化を無視している。これらの仮定を置いたうえで、励起光強度分布Ｈ_ｘと蛍光強度分布Ｅ_ｍは時間的に変化しないものとする。これは組織６内部の吸収や拡散は、心臓の拍動中でも蛍光のスペクトルシフトによってもほぼ変化しないためである。ここでは心筋の収縮については考慮しない。計測される組織表面６ａ上の蛍光の変化は次式から得られる。

ここで、励起光強度分布Ｈ_ｘと蛍光強度分布Ｅ_ｍは、モンテカルロシミュレーションによってモデル化することができる。そこで、励起光強度分布Ｈ_ｘと蛍光強度分布Ｅ_ｍは、予めモンテカルロシミュレーションによってモデル化しておく。モンテカルロシミュレーションで励起光強度分布Ｈ_ｘと蛍光強度分布Ｅ_ｍを予めモデル化するためには下記のような考え方に基づく。

例えば心筋組織を計測対象とする場合には、当該心筋組織の光学定数を、励起光波長および蛍光波長で予め調べる必要がある。ここで、光学定数とは、励起光のときの組織６の拡散係数、吸収係数、異方性パラメータの設定、空気の屈折率、組織６の屈折率、蛍光のときの組織６の拡散係数、吸収係数、異方性パラメータの設定、空気の屈折率、組織６の屈折率、励起光波長、蛍光波長等である。そして、得られた光学定数を元に、図５Ａに示すように、励起光の波長で励起光照射位置ｒ_ｓに入射したときの励起光強度分布Ｈ_ｘを求める。これによって位置ｒにおける励起光の強度が得られる。

位置ｒで励起された色素が発する蛍光は、蛍光計測位置ｒ_ｄで計測されるが、ここでは蛍光計測位置ｒ_ｄに照射した蛍光波長が位置ｒに達し、そのときの強度をみる。これはフォトンが位置ｒから発生し、蛍光計測位置ｒ_ｄに達する確率と、蛍光計測位置ｒ_ｄに照射した蛍光波長が位置ｒに達する確率と、がほぼ同じで近似できるためである。そして、これらから得られた励起光強度分布Ｈ_ｘと蛍光強度分布Ｅ_ｍを、各ボクセルで掛け合わせることにより、計測感度行列Ｊ_ｘ,ｍとして統合することができる。計測感度行列Ｊ_ｘ,ｍのｘは励起光に依存する関数であることを示し、ｍは蛍光に依存する関数であることを示す。計測感度行列Ｊ_ｘ,ｍは、励起光照射位置ｒ_ｓと蛍光計測位置ｒ_ｄで関連付けられる行列となっており、図５Ｂに示すようなものとなる。上述した数５は以下のように表される。ｔｏは基準時間を示す。

ここで、蛍光分布画像から計測される静止状態の蛍光Ｆ_ｘ,ｍと、求めたい膜電位分布Ｖは、線型的な関係式にあることから圧縮センシングの導入が可能となる。

（３）上述の計測原理を利用した本発明の計測システムについて
計測システム１では、照射点分布パターンをランダムに変えながらパターン光Ｌ１を組織表面６ａに照射してゆき、励起光を多点照射する照射点分布パターン毎に蛍光分布画像を取得している。そこで、励起光を同時に多点照射する照射点分布パターン毎に蛍光分布画像を取得する計測システム１に、上述した計測原理を適用した例について以下説明する。例えば、蛍光分布画像を取得する回数（計測回数とも呼ぶ）がｐ回であるとき、ｋ番目のときの組織表面６ａ上におけるパターン光の励起光強度Ｐ_ｋは下記の式で表すことができる。

ここで（ｕ，ｖ）は、計測対象となる組織表面６ａ上での座標系であり（図１）、ｈ_ｅｘは対物レンズ１３の励起光のＰｏｉｎｔ Ｓｐｒｅａｄ Ｆｕｎｃｔｉｏｎ（ＰＳＦ）を示す。ｕ^ｇ _ｋはｋ番目の照射点分布パターンの幾何学イメージであり、ｇは幾何座標を示し、次のように表される。（ξ，η）は、マイクロミラー１６がアレイ状に配置されたＤＭＤ１２において、マイクロミラー１６の配置面上での座標系を示す。

ここでｕ_ｋ（α，β）は、ＤＭＤ１２のピクセル位置（α，β）におけるＤＭＤパターンを示しており、それぞれの位置にあるマイクロミラー１６がオン状態またはオフ状態を取る。△_Ｄは、マイクロミラー１６を正方形としたときのマイクロミラー１６の一辺のサイズを示す。励起光照射の後に時間ｔにおける組織表面６ａ上での蛍光は次のように表される。

Ｏ（ｕ，ｖ，ｔ）は、計測対象となる組織６の反応強度（励起光強度に対する蛍光強度）を示し、Ｈ（ｘ，ｙ，ｚ）は、照射点分布パターンによって生じる組織６中の励起光強度の３次元分布（励起光強度分布）を示す。Ｅ（ｘ−ｕ，ｙ−ｖ，ｚ）は、蛍光強度の３次元分布（蛍光強度分布）を示す。（ｘ，ｙ，ｚ）は、組織６内での３次元座標系を示す。ε（ｘ，ｙ，ｚ）は色素の吸光係数を示し、ｃ（ｘ，ｙ，ｚ）は色素濃度を示し、τ（ｘ，ｙ，ｚ，ｔ）は色素の量子効率を示す。ｋ番目の照射点分布パターンによって生じた組織表面６ａ上の蛍光分布を撮像機器３で撮像した際、当該撮像機器３で取得した蛍光分布画像内での蛍光Ｉ_ｋは、以下のように表される。

（ｎ，ｍ）は、撮像機器３のセンサアレイ上の連続的な値の座標系である。ｈ_ｅｍは蛍光波長でのイメージングレンズ１５のＰＳＦを示す。Ｅ^ｋ _ｇ（ｎ，ｍ，ｔ）は、ｋ番目の照射点分布パターンにおいて撮像機器３のセンサアレイ上での蛍光強度を示しており（ｇは幾何座標を示す）、下記のように表される。

Ｍ_ｄはイメージングレンズ１５の拡大率を示す。ｋ番目の蛍光分布画像における蛍光輝度値ｙ_ｋは、撮像機器３のディテクタ位置（ａ，ｂ）においては以下のようになる。ディテクタ位置（ａ，ｂ）とは、センサアレイに配置されたディテクタの位置（ピクセルの位置）を示す。

ρ_ａｂはディテクタ位置（ａ，ｂ）におけるシステムノイズを示す。また、△_ｄｅｔは、ディテクタ（ピクセル）を正方形としたときのディテクタの一辺のサイズを示す。例えば、１番目の照射点分布パターンからｐ番目の照射点分布パターンまでの蛍光輝度値ｙ_１，…，ｙ_ｐをまとめたものを蛍光輝度行列Ｙとし、上述した数１２を簡便に表すと、△Ｙ＝Ｊ△Ｖ＋ρと表すことができる。△Ｙ＝Ｊ△Ｖ＋ρのうちρはノイズ行列である。

ここで、計測感度行列Ｊは、上記の「（２）ＬＯＴにおける計測原理について」にて説明したように、励起光強度分布Ｈと蛍光強度分布Ｅの各ボクセルを掛け合わせることで得られる。計測感度行列Ｊは、撮像機器３の各ディテクタでそれぞれ計測される蛍光の空間感度を表す。膜電位分布Ｖは、計測対象となる組織６の反応強度の支配的な要因となる。従って、上述した数１０において、膜電位分布Ｖは次のように表される。ｓはスケーリング定数を示す。

（４）圧縮センシングの適用について
圧縮センシングの目的はサブナイキスト周波数で少ない計測値（蛍光）から求めたい係数の値（膜電位分布Ｖ）を決定することである。先ず、図６に示すように、一般的な圧縮センシングを示すｙ＝Φｘの線形方程式を考える。ここではｙをＭ×１の計測ベクトルとし、ΦはＭ×Ｎの基底行列とする。ｘがＮ×１の係数ベクトルである。この係数ベクトルはスパースであり、高々Ｋ個（<<Ｎ）の要素が非ゼロであるとする。少ない計測値から係数ベクトルを再構成するためｌ_１ノルム最小化を考える。ｙ＝Φｘの線形方程式の中で、解をｘとすると、ｌ_１ノルム最小化によって、非常に少ない計測ベクトルｙから、解であるｘを精度よく求めることができる。

図６に示すように、一般的な圧縮センシングを示すｙ＝Φｘの線形方程式は、行列をスパース行列にするための空間変換行列Ψを用いることで、ｙ＝ΦΨαの線形方程式で表すことができる。αはスパース行列であり、「ΦΨ」を新たな基底行列として考えることで、スパース行列αを含んだ線形方程式として計測理論モデルを扱うことができる。

計測システム１では、このような圧縮センシングの理論を利用して、膜電位分布Ｖを求めるが、膜電位分布Ｖは非ゼロ要素数が多い疎な行列（スパース行列）とは言えない。よって、この膜電位分布Ｖをスパースにするために、上述した圧縮センシングの理論の空間変換行列Ψの作用を考える。上述した△Ｙ＝Ｊ△Ｖ＋ρの線形方程式（Ｙは蛍光輝度行列、Ｊは計測感度行列、Ｖは膜電位分布）を、空間変換行列Ψおよびスパース行列αを用いて簡単に表すと、次のように表すことができる。

図７は、Ｙ＝ＪΨαの線形方程式について、圧縮センシングの理論を適用したときのデータ構造を示した概略図である。計測画像となる蛍光分布画像は、ＤＭＤ１２で照射点分布パターンを変えた回数取得する。各蛍光分布画像Ｉ_ｍ１，…，Ｉ_ｍｐのそれぞれの総ピクセル数をＮとし、ＤＭＤ１２で作製された照射点分布パターンの数をＰ個（パターン１，…，パターンｐ）とすると、全Ｐ個の照射点分布パターン毎に計測した蛍光輝度値の全サイズはＮ×Ｐとなり、これが蛍光輝度行列Ｙのデータ構造となる。

ＤＭＤ１２で作られる照射点分布パターンは圧縮センシングを適用できるようにランダム性を持ったパターンとなっている。計測感度行列Ｊは、照射点分布パターン内のそれぞれの励起点に対して、撮像機器３の各ディテクタ用にＭ個のモンテカルロシミュレーションによって作成する。最終的に求めたい膜電位分布Ｖのサイズを、例えばＸ×Ｙ×Ｚとすると、各計測感度行列ＪのサイズもＸ×Ｙ×Ｚとなる。計測感度行列Ｊについては、Ｐ個の照射点分布パターン毎に、それぞれ撮像機器３のディテクタ数となるピクセル数Ｎだけ計測感度行列ｊを作成する必要がある。

例えば照射点分布パターンがパターン１のときには（図７）、パターン１により照射された複数点の励起光に対して、１ピクセル目での計測感度行列ｊ_１、２ピクセル目での計測感度行列ｊ_２、３ピクセル目での計測感度行列ｊ_３がそれぞれ求められ、最後のＮピクセル目まで計測感度行列ｊ_Ｎが求められる。結果として、照射点分布パターンがパターン１〜ｐまでの計測感度行列Ｊは、データ構造が（Ｎ×Ｐ）×（Ｘ×Ｙ×Ｚ）となる。

ここで、圧縮センシングの理論の下、膜電位分布Ｖをスパースの行列にするためには、上述したように空間変換行列Ψを用いる必要がある。この空間変換行列ΨのサイズはＸ×Ｙ×Ｚである。スパース行列αのサイズもＸ×Ｙ×Ｚである。求めたい膜電位分布Ｖは、空間変換行列Ψ×スパース行列αと表すことができる。蛍光分布画像Ｉ_ｍ１，…，Ｉ_ｍｐからそれぞれ得られる蛍光輝度値ｙ_１，…，ｙ_ｐをまとめた蛍光輝度行列Ｙは、計測感度行列Ｊ×空間変換行列Ψ×スパース行列αで表すことができる。

本発明による計測システム１では、この関係式を用いて、最初に、計測感度行列Ｊと空間変換行列Ψを新たな基底行列として扱い、蛍光分布画像Ｉ_ｍ１，…，Ｉ_ｍｐにより求めた蛍光輝度値ｙ_１，…，ｙ_ｐ（蛍光輝度行列Ｙ）から、スパース行列αをｌ_１ノルム最小化法によって再構成する。得られたスパース行列αを空間変換行列Ψの逆変換によって、求めたい膜電位分布Ｖに変換する。これにより、計測画像である蛍光分布画像Ｉ_ｍ１，…，Ｉ_ｍｐから、圧縮センシングを用いてスパースではない膜電位分布行列を得ることができる。

ここで、圧縮センシングにおいては基底行列のインコヒーレンス性が重要となるが、本発明による手法の基底行列は、計測感度行列Ｊと空間変換行列Ψを掛け合わせたものとなる。インコヒーレンス性を変化させる要因は、ＤＭＤ１２で生成される照射点分布パターンの変化のみによる。よってＤＭＤ１２の照射点分布パターンをランダム性のあるパターンにすることで、インコヒーレンス性を高めることができる。

ところで、二分法的に言えば、図８に示すように、例えば心筋組織における膜電位は電気的に高い状態と低い状態になる。心筋組織は、高い状態になった後は不応期に入り、外部刺激に対して反応しない状態をとる。ゆえに、本発明の計測システム１の計測対象となる膜電位分布Ｖは、組織６中で電気的に高い部位が塊として存在し、空間的に低周波の物であると考えられる。

離散コサイン変換（Discrete Cosine Transformation：ＤＣＴ）や、離散ウェーブレット変換(Discrete Wavelet Transformation：ＤＷＴ)のような一般的な周波数変換は、逆変換もでき、完全再構成も保証されており、スパースに表現できる有用な方法である。このような計測対象のスパース性は計測回数（照射点分布パターン毎に得られる蛍光分布画像の取得回数）の削減につながってくる。

本発明の計測システム１で用いる変換手法としては、離散コサイン変換（ＤＣＴ）と、離散ウェーブレット変換(ＤＷＴ)を用いる。ＤＣＴはコサイン波の重ね合わせで対象を表現する手法となっており、以下のように得られる。

ここで、Χ_ｋは、複数の照射点分布パターンのうちｋ番目（ｋ＝０，…，Ｎ−1）の照射点分布パターンのＤＣＴ係数であり、データ長さをＮ、元の信号をｘで示す。一方で、ＤＷＴは、マザーウェーブレットと呼ばれる波を伸縮・拡大させて対象を表現する手法となっている。本発明では、マザーウェーブレットとして代表的なＨａａｒウェーブレット（以下、単にＨａａｒとも呼ぶ）の、Ｃｏｉｆｌｅｔウェーブレット（以下、単にＣｏｉｆｌｅｔとも呼ぶ）の、Ｄａｕｂｅｃｈｉｅｓウェーブレット（以下、単にＤａｕｂｅｃｈｉｅｓとも呼ぶ）を用いることができる。ＤＣＴおよびＤＷＴは、３次元変換に拡張して対象となる膜電位分布Ｖを変換する。ここでＤＷＴにおいてはマルチレベル分解構造を取ることができる。これは１次元のＤＷＴで考えたとき、ＤＷＴによって元の信号を高周波成分と低周波数成分に分割する。これを1レベル（level１）の分解と呼び以下のように表す。

ここでｙ_ｌｏｗは低周波数成分を示し、ｙ_ｈｉｇｈは高周波数成分を示し、カッコ内のｎは離散変数を示す。ｘ［ｋ］は、ｋ番目の元信号を示し、ｇは低周波数フィルタ、ｈは高周波数フィルタを示す。ＤＷＴでは、ここで得られた低周波数成分について、図９に示すように、高周波数成分（ｈ［ｎ］）と低周波数成分（ｇ［ｎ］）とに再度分解することができる。ここでは、１レベルの低周波成分を再度分解して高周波数成分と低周波数成分を得る分解を、２レベル（level２）の分解と呼ぶ。このような高周波成分と低周波成分の分解は連続して行うことができ、元信号が２^Ｎ個のデータ数を持つとき、最大でＮレベル（levelＮ）の分解を行うことができる。

この連続的な分解を考えると、蛍光分布画像内の蛍光が、ある特定の周波数をもつ場合においては、ある分解レベルで最もスパースになることが予想される。３次元変換で考えた場合には、図１０に示すように、ｘ方向、ｙ方向、ｚ方向にそれぞれ異なる分解レベルを適用することができる。ここでは、計測対象となる組織６がｘ方向、ｙ方向で興奮伝播が異なる特性を持たないことが考えられるため、ｘｙ平面においては異なる分解レベルを適用しても、ほぼ変化はないものとしてｘ方向、ｙ方向に対しては同様の分解レベルを適用することとした。すなわち、ｘｙ平面とｚ方向とで異なる分解レベルを適用することとした。

図１０に示すように、３次元データをｘ方向、ｙ方向、ｚ方向にそれぞれ１レベルの分解レベルとしたとき、レベル（ｘ，ｙ，ｚ）＝（１，１，１）と表し、ｘｙ平面とｚ方向の分解レベルとして１．１と表す。分解レベルが１．１の３次元データのうち、低周波数成分である「ＬＬＬ」をｚ方向にさらに分解したときは、レベル（ｘ，ｙ，ｚ）＝（１，１，２）＝１．２と表す。分解レベルが１．１の３次元データのうち、低周波数成分である「ＬＬＬ」をｘ方向、ｙ方向にそれぞれ分解したときには、レベル（ｘ，ｙ，ｚ）＝（２，２，１）＝２．１と表す。分解レベルが１．１の３次元データのうち、低周波数成分である「ＬＬＬ」をｘ方向、ｙ方向、ｚ方向にそれぞれ分解したときには、レベル（ｘ，ｙ，ｚ）＝（２，２，２）＝２．２と表す。

ＤＷＴにおいては、フィルタとして使用するマザーウェーブレットをいくつか用意しておき、蛍光分布画像内の蛍光分布や、求める膜電位分布の特性に応じて選択することが望ましい。例えば、１次元波形(128 data points)のステップ関数ではＨａａｒが最もスパースとなり、１次元波形のコサイン波ではＤＣＴが最もスパースになり、１次元波形の三角波ではＤａｕｂｅｃｈｉｅｓ２（Ｄｂ２）が最もスパースになることが確認できている。活動電位波形においては大きな差はなかったが、波形の違いによりＣｏｉｆｌｅｔまたはＨａａｒが最もスパースとなることを確認している。

マザーウェーブレットの選択については計測対象と類似した波形の選択が有効である。点蛍光体のような場合においてはステップ関数とみなすことができるため、Ｈａａｒが有効である。膜電位分布の場合には、計測視野範囲を２０ミリ四方とすれば、その電位分布はステップ関数上になるためＨａａｒが有効である。一方、計測視野範囲を１０［ｃｍ］四方とした場合には、活動電位波形のような信号分布となるため、Ｃｏｉｆｌｅｔのような形状のマザーウェーブレットを選択することが有効である。

（５）演算処理装置について
次に、圧縮センシングの理論を利用して、膜電位分布Ｖを求める演算処理装置４について以下説明する。図１１に示すように、演算処理装置４は、図示しないＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）およびＲＡＭ（Random Access Memory）等からなるマイクロコンピュータ構成の制御部３１と、外部インターフェース３２と、記憶部３３と、画像処理部３４と、キーボードやマウスからなる操作部３５と、再構成演算部３６と、膜電位分布演算部３７と、蛍光輝度計測部３９と、がバス４０を介して相互に接続された構成を有する。画像処理部３４には表示部３８が接続されており、外部インターフェース３２には、光学機器２および撮像機器３が接続されている。

制御部３１は、操作部３５を介して、作業者から各種操作命令が与えられると、ＲＯＭに予め格納している演算処理プログラム等を、操作命令に基づき適宜読み出してＲＡＭに展開することにより、演算処理プログラム等に従って各回路部を制御する。制御部３１は、各種プログラムの実行結果を、画像処理部３４を介して表示部３８に表示させ、例えば、表示部３８を介して、組織６における膜電位分布Ｖを検査作業者に把握させ得る。

制御部３１は、光学機器２のＤＭＤ１２に対して外部インターフェース３２を介して制御信号を送出し、ＤＭＤ１２の各マイクロミラー１６をオンオフ動作させる。光学機器２では、ＤＭＤ１２により生成されるパターン光Ｌ１の照射点分布パターンをランダムに変えながら、組織表面６ａにパターン光Ｌ１を照射する。外部インターフェース３２は、照射点分布パターンがランダムに変わるパターン光Ｌ１が照射されている組織表面６ａを撮像した蛍光分布画像を撮像機器３から受け取り、記憶部３３に記憶する。なお、撮像機器３は、光学機器２と同期しており、光学機器２が異なる照射点分布パターンのパターン光を生成する毎に蛍光分布画像を取得する。

記憶部３３には、スパース行列αを求める際に用いる空間変換行列Ψと、膜電位分布Ｖを求める際に用いる計測感度行列Ｊと、が予め記憶されている。ここで、空間変換行列Ψとしては、上述したように、ＤＣＴ行列（離散コサイン変換行列）や、ＤＷＴ行列（離散ウェーブレット変換行列）が記憶されている。ＤＷＴ行列としては、例えば各分解レベルのＨａａｒ、Ｃｏｉｆｌｅｔ、Ｄａｕｂｅｃｈｉｅｓの各行列が記憶されている。計測感度行列Ｊは、計測対象とする組織６の拡散係数、吸収係数、異方性パラメータの設定等の光学定数を用いて、照射点分布パターン毎にモンテカルロシミュレーションを利用して予め計算して求めたものである。

蛍光輝度計測部３９は、ランダムに変わる照射点分布パターン毎に得られた蛍光分布画像から蛍光輝度値を求め、これら蛍光輝度値からなる蛍光輝度行列Ｙを求める。なお、この実施形態の場合、蛍光輝度計測部３９を演算処理装置４に設けた場合について述べたが、本発明はこれに限らず、蛍光輝度計測部を撮像機器３に設けるようにしてもよい。この場合、撮像機器３は、蛍光輝度計測部により、ランダムに変わる照射点分布パターン毎に得られた蛍光分布画像から蛍光輝度行列Ｙを求め、これを計測結果として演算処理装置４に送出する。

再構成演算部３６は、組織表面６ａに照射したパターン光Ｌ１の各照射点分布パターンに対応する計測感度行列Ｊを、記憶部３３から読み出す。具体的には、例えば、どの位置を照射点位置とするかについて予め決めた、幾つかの照射点分布パターンを用意しておく。そして、用意した各照射点分布パターン毎に、それぞれラスタースキャンにより１個目のディテクタから最後のＮ個目のディテクタまでの各計測感度行列ｊ_１，…，ｊ_Ｎを予め求めておき、各照射点分布パターンに対応する計測感度行列Ｊ（１個目からＮ個目までのディテクタの計測感度行列ｊ_１，…，ｊ_Ｎの集合）を記憶部３３に記憶しておく。

再構成演算部３６は、ＤＭＤ１２で照射されたパターン光Ｌ１が、どの照射点分布パターンであったかを、光学機器２や制御部３１等からの情報により認識し、当該照射点分布パターンに対応した計測感度行列Ｊを記憶部３３から読み出す。

また、再構成演算部３６は、例えば検査作業者により選ばれたＤＣＴ行列を、空間変換行列Ψとして記憶部３３から読み出す。再構成演算部３６は、圧縮センシングの理論の下、計測感度行列Ｊと空間変換行列Ψを掛け合わせ、これら計測感度行列Ｊと空間変換行列Ψとを基底行列として扱い、ｌ_１ノルム最小化法により蛍光輝度行列Ｙからスパース行列αを再構成する。

膜電位分布演算部３７は、再構成演算部３６で求めたスパース行列αを、当該スパース行列αを求める際に用いた空間変換行列Ψ（ＤＣＴ行列）の逆変換により、膜電位分布Ｖを求める。画像処理部３４は、求められた膜電位分布Ｖを３次元情報として表示部３８に表示し、検査作業者に対して組織６の３次元的な膜電位分布状況を呈示し得る。

この際、膜電位分布演算部３７は、例えば、再構成された各ボクセルを、電気的に２値化した後、所定の閾値を基にＨｉｇｈまたはＬｏｗのどちらの状態にする。その後、膜電位分布演算部３７は、モルフォロジー演算処理を実行し、膜電位分布Ｖを３次元画像として表示部３８に表示する。モルフォロジー演算処理では、周りの全てのピクセルがＬｏｗ（またはＨｉｇｈ）であるが中心１ピクセルだけＨｉｇｈ（またはＬｏｗ）の状態のとき、このＨｉｇｈ（Ｌｏｗ）のピクセルをＬｏｗ（Ｈｉｇｈ）に塗りつぶす処理を行う。これにより、膜電位分布演算部３７は、周囲と電位状態が異なる中心の１ピクセルは明らか誤差と見なし、誤差の少ない３次元情報を生成できる。

（６）膜電位分布計測手順について
次に、上述した組織６の膜電位分布計測手順について、図１２に示すフローチャートを用いて以下簡単に説明する。図１２に示すように、開始ステップＲＴ１から膜電位分布計測手順を開始し、次のステップＳＰ１に移る。ステップ１において、計測対象とする組織６の光学定数を決定するとともに、組織表面６ａに照射するパターン光の照射点分布パターンを決定し、これら光学定数および照射点分布パターンの各情報を演算処理装置４の記憶部３３に記憶し、次のステップ２に移る。ここで、組織６の光学定数の決定としては、励起光のときにおける組織６の拡散係数、吸収係数、異方性パラメータの設定、空気の屈折率、組織６の屈折率、蛍光のときにおける組織６の拡散係数、吸収係数、異方性パラメータの設定、空気の屈折率、組織６の屈折率、励起光波長、蛍光波長等を決定する。照射点分布パターンの決定としては、励起光の照射点数と、照射点分布パターンが異なるパターン光の照射回数（計測回数）とを決定する。

ステップＳＰ２において、演算処理装置４の制御部３１は、ステップＳＰ１で決定した光学定数を用いてモンテカルロシミュレーションにより計測感度行列Ｊを計算し、得られた計測感度行列Ｊを記憶部３３に記憶して、次のサブルーチンＳＲＴ１に移る。サブルーチンＳＲＴ１において、計測システム１は、画像取得処理を実行して、各照射点分布パターン毎に蛍光分布画像を取得してゆき、これら蛍光分布画像から蛍光輝度行列Ｙを求め、次のステップＳＰ３に移る。

ステップＳＰ３において、演算処理装置４の再構成演算部３６は、スパース行列αを求める際に用いる空間変換行列Ψを選択し、次のステップＳＰ４に移る。ここで、空間変換行列Ψの選択は、記憶部３３に予め記憶されたＤＣＴ行列やＤＷＴ行列の中から、例えば操作部３５を介して検査作業者により選択させてもよい。その他として、計測対象となる組織６の種類に合わせて最適な種類の空間変換行列Ψを予め対応付けておき、組織６の種類が再構成演算部３６に設定されることで、当該組織６に対応付けられた空間変換行列Ψを自動的に選択させるようにしてもよい。

ステップＳＰ４において、演算処理装置４の再構成演算部３６は、計測感度行列Ｊと空間変換行列Ψとを基底行列として扱い、後述するサブルーチンＳＲＴ１で得た蛍光輝度行列Ｙから、ｌ_１ノルム最小化法によりスパース行列αを再構成し、次のステップＳＰ５に移る。

ステップＳＰ５において、演算処理装置４の膜電位分布演算部３７は、ステップＳＰ４で用いた空間変換行列Ψによる逆変換によって、スパース行列αを膜電位分布Ｖに変換し、次のステップＳＰ６に移り、膜電位分布計測手順を終了する。

ここで、画像取得処理を実行するサブルーチンＳＲＴ１は、図１３に示すように、開始ステップＲＴ２から開始し、次のステップＳＰ１１に移る。ステップＳＰ１１において、光学機器２は、ＤＭＤ１２のマイクロミラー１６をオンオフ動作し、照射点分布パターンがランダムに変わるパターン光を生成して次のステップＳＰ１２に移る。

ステップＳＰ１２において、光学機器２は、ＤＭＤ１２で生成された照射点分布パターンのパターン光を、計測対象である組織６に照射し、次のステップＳＰ１３に移る。ステップＳＰ１３において、撮像機器３は、パターン光が照射された組織表面６ａにおける蛍光分布を撮像してゆき、蛍光分布画像を取得して次のステップＳＰ１４に移る。

ステップＳＰ１４において、演算処理装置４の蛍光輝度計測部３９は、ランダムに変わる照射点分布パターン毎にそれぞれ得られた蛍光分布画像を、撮像機器３から受け取る。ステップＳＰ１４において、蛍光輝度計測部３９は、各蛍光分布画像の蛍光輝度値からなる蛍光輝度行列Ｙを求め、次のステップＳＰ１５に移り、画像取得処理を終了する。

（７）作用および効果
以上の構成において、計測システム１では、照射点分布パターンをランダムに変えながらパターン光Ｌ１を組織表面６ａに照射してゆく光学機器２と、照射点分布パターン毎に組織表面６ａにおける蛍光分布を撮像して組織表面６ａの蛍光分布画像を取得する撮像機器３と、蛍光分布画像に基づいて組織６の膜電位分布Ｖを求める演算処理装置４とを設けるようにした。演算処理装置４では、励起光および蛍光の強度分布関係をモンテカルロシミュレーションから求めた計測感度行列Ｊと、計測感度行列Ｊをスパース行列αにするための空間変換行列Ψと、を基底行列として扱い、蛍光分布画像から求めた蛍光輝度行列Ｙからスパース行列αを再構成するようにした。演算処理装置４では、再構成演算部３６で再構成したスパース行列αを、空間変換行列Ψの逆変換によって膜電位分布Ｖに変換することで、膜電位分布Ｖを求めることができる。

これにより、計測システム１では、パターン光Ｌ１を組織表面６ａに照射することで、蛍光を多点で同時計測できる蛍光分布画像を得ることができ、１点１点蛍光を計測していた従来の計測システムに比べて、膜電位分布Ｖを高速に計測し得る。また、計測システム１では、非ゼロ要素の数が少ないスパース行列αを再構成し、逆変換によってスパースでない膜電位分布Ｖを求めるようにしたことで、圧縮センシングの理論の下、少ない枚数の蛍光分布画像からでも膜電位分布Ｖを求めることができる。よって、計測システム１では、蛍光分布画像の枚数を減らせる分、処理負担が軽減し、膜電位分布を高速に計測し得る。

なお、上述した実施形態においては、照射点分布パターンをランダムに変えながらパターン光Ｌ１を組織表面６ａに照射してゆき、照射点分布パターン毎に蛍光分布画像を取得し、照射点分布パターン毎に得られた複数の蛍光分布画像から蛍光輝度行列Ｙを求めるようにした。しかしながら、本発明はこれに限らず、１つの照射点分布パターンのパターン光Ｌ１を組織表面６ａに照射し、１つの照射点分布パターンで得られた１枚の蛍光分布画像から、当該蛍光分布画像の蛍光分布を示す蛍光輝度値を求めてもよい。この場合、再構成演算部３６は、計測感度行列Ｊおよび空間変換行列Ψを基底行列として扱い、蛍光分布画像から求めた蛍光輝度値からスパース行列αを再構成する。

また、上述した実施形態においては、計測感度行列Ｊを求めるシミュレーションとして、モンテカルロシミュレーションを適用した場合について述べたが、本発明はこれに限らず、放射輸送方程式の他、拡散近似法を用いた解析的な手法を適用しても良い。

（８）シミュレーション実験
次に、図１に示した計測システム１を用いて膜電位分布を求めるシミュレーション実験を行った。計測対象モデルはウサギの心臓の活動電位モデルを用い、図１４Ａ〜図１４Ｆに示すような６つの代表的な膜電位分状態を対象とした。電位は−８０［ｍＶ］から２０［ｍＶ］に変化するとし、活動電位モデルのサイズは２０［ｍｍ］×２０［ｍｍ］×５［ｍｍ］（３２×３２×８ボクセル（合計８１９２ボクセル））とし、各ボクセルサイズを６２５［μｍ］とした。電気的な興奮伝播の速度を３０［ｃｍ／ｓ］とした。

図１４Ａは、上部に高電位領域ＥＲ_Ｈを設け、下部に低電位領域ＥＲ_Ｌを設けて垂直方向に上から下に興奮伝播している状態を表した活動電位モデルを示す。図１４Ｂは、下部に高電位領域ＥＲ_Ｈを設け、上部に低電位領域ＥＲ_Ｌを設けて垂直方向に下から上に興奮伝播している状態を表した活動電位モデルを示す。図１４Ｃは、右側に高電位領域ＥＲ_Ｈを設け、左側に低電位領域ＥＲ_Ｌを設けて垂直な波面で興奮伝播している状態を表した活動電位モデルを示す。

図１４Ｄは、右側に高電位領域ＥＲ_Ｈを設け、左側に低電位領域ＥＲ_Ｌを設けて−４５度に傾いた波面で興奮伝播している状態を表した活動電位モデルを示す。図１４Ｅは、右側に高電位領域ＥＲ_Ｈを設け、左側に低電位領域ＥＲ_Ｌを設けて＋４５度に傾いた波面で興奮伝播している状態を表した活動電位モデルを示す。図１４Ｆは、中心に球状の高電位領域ＥＲ_Ｈを設け、その周辺を取り囲むように低電位領域ＥＲ_Ｌを設けて球面波で興奮伝播している状態を表した活動電位モデルを示す。組織６の膜電位分布は、これら６つの単純な活動電位モデルで表現できるものと考えた。

先ず最初に、これら６つの活動電位モデルが周波数変換において、どれほどのスパース性を持つかについて評価を行った。空間変換行列Ψを用いた変換手法としては、ＤＣＴとＤＷＴとを用いた。評価手法としては、全８１９２ボクセル内の非ゼロ要素数をカウントすることでスパース性を評価した。また、係数のエネルギーからエネルギー評価を行った。ＤＷＴを用いたときの非ゼロ要素数を、図１５Ａに示す。

図１５Ａでは、全ボクセル８１９２（３２×３２×８ボクセル）中の非ゼロ要素数を、各膜電位分布モデル（図１４Ａ〜図１４Ｆ）においてカウントし、平均を四捨五入した値を示している。図１５中のカッコ内は各膜電位分布モデルのばらつきを示す標準偏差である。図１５Ａ中の「Ｌｅｖｅｌ」は、図９および図１０にて説明した分解レベルを示す。ＤＣＴ行列を使用した場合には８１４３（２５）という結果となった。

図１５Ａから、ＤＷＴの変換手法によって非ゼロ要素数は大きく変化することが確認できた。また、分解レベルによっても、非ゼロ要素数が変化することが確認できた。Ｈａａｒの分解レベル１．３の変換が、最もスパース性が高く、１５３９（１３２０）という結果であった。全体としてはＤａｕｂｅｃｈｉｅｓやＣｏｉｆｌｅｔではスパース性がＨａａｒに比べて低かった。

３次元の膜電位分布モデルでは、興奮伝播波面を境界として高電位領域ＥＲ_Ｈと低電位領域ＥＲ_Ｌが存在する。Ｈａａｒでは高電位領域ＥＲ_Ｈや低電位領域ＥＲ_Ｌのフィルタの結果は値を持つか０かのような結果となる。その結果としてＤａｕｂｅｃｈｉｅｓやＣｏｉｆｌｅｔのような複雑な波形を持ったフィルタに比べてＨａａｒはスパースになると考えられる。

Ｈａａｒの変換において、ｘｙ平面のどの分解レベルにおいても、ｚ方向において３レベルのときに最もスパースになった。これはシミュレーション実験で用いた膜電位分布モデルがｚ方向に対して空間周波数が低いものが支配的であったためであると考えられる。シミュレーション実験で用意した膜電位分布モデルは、ｚ方向に対して低い周波数を有するものであったが、いくつかの電位分布モデルについてはｘｙ平面で高い周波数を持つものであった。その点を考慮すると、結果的にｘｙ平面において分解レベル１、ｚ方向で分解レベル３のＨａａｒのときに最もスパース性の高い結果となったことに妥当性がある。

エネルギーに関しては変換後の係数のうち、大きいものから５０個の係数を選択し、そのエネルギーが全体の何パーセントに当たるかを評価した。その結果、図１５Ｂに示すような結果が得られた。ＤＣＴでは９０［％］（１４）という結果であった。エネルギーについても各変換によって大きく変化し、分解レベルによっても変化することを確認した。結果として、Ｈａａｒの分解レベル３．３が最も高いエネルギーであった。

図９に示したように、ＤＷＴでは、高周波係数は連続的に分解されず、その後の分解でも残存する。そのために特定の分解レベルにおいて係数の数が最も小さくなることが分かった。結果として、ｘｙ平面において分解レベル３、ｚ方向において分解レベル３のＨａａｒのときに最も高いエネルギーとなったと考えられる。

次に、図１４Ａ〜図１４Ｆに示した各膜電位分布モデルを計測対象として用い、パターン光における励起光の照射点数を変えたときに、再構成の精度がどのように変わるかについて確認した。ここでは、空間変換行列を用いた変換として、高スパース性のＨａａｒ １．３と、高エネルギーのＨａａｒ ３．３と、ＤＣＴと、を用いた。一般的に用いられている膜電位感受性色素Ｄｉ−４−ＡＮＥＰＰＳの励起光波長および蛍光波長に近いものにするため、励起光波長は５５０［ｎｍ］、蛍光波長は６００［ｎｍ］とした。計測感度行列Ｊはモンテカルロシミュレーションを用いて計算した。モンテカルロシミュレーションにて使用する、拡散係数、吸収係数、異方性パラメータの設定等の光学定数は、図１６に示す。図１６に示す光学定数は、ラットの筋組織の光学定数である。

パターン光の照射点分布パターンは、ランダム行列パターン（３２×３２ pixel）とした。パターン光の照射点直径は１［ｍｍ］とし、パターン光内の照射点数のみ固定し、ランダムに照射点分布パターンを変化させた。１回の再構成で使用する照射点分布パターンは、５パターン、１０パターン、２０パターンとした。すなわち、計測回数５回、１０回、２０回とした。ここでは、広く膜電位計測で使用されている膜電位感受性色素Ｄｉ−４−ＡＮＥＰＰＳを用いるが、この蛍光変化率は経験的に９［％］と定めた。Ｐｏｉｎｔ Ｓｐｒｅａｄ Ｆｕｎｃｔｉｏｎ（ＰＳＦ）については理想的であるとして、光学機器２内にある全てのレンズの倍率は１とした。センサアレイについては３２×３２ pixel（ディテクタが３２×３２に配置）で組織表面上の関心領域（Region of interest：ＲＯＩ）を１倍で投影する状況を考えた。

パターン光内における照射点数を１０点から９０点まで１０点ずつの間隔で変化させた照射点分布パターンのパターン光と、当該照射点数を１００点から１０００点まで１００点ずつの間隔で変化させた照射点分布パターンのパターン光とを用いてそれぞれ評価を行った。図１７には、ランダムに励起光を１００点照射したパターン光と、ランダムに励起光を５００点照射したパターン光と、ランダムに励起光を８００点照射したパターン光と、についてそれぞれ示したものである。再構成に用いる蛍光分布画像の画像枚数（すなわち、計測回数）の変化が再構成の精度に与える影響について検討するため、照射点数はそのままで蛍光分布画像の画像枚数を５枚、１０枚、２０枚と変化させて再構成を行った。

再構成において、重要となるのは膜電位分布の形状である。すなわち、組織の脱分極と再分極の界面の動きを知ることができればよい。再構成されたボクセルは電気的にＨｉｇｈ（脱分極）、またはＬｏｗ（再分極）の状態にあると言える。膜電位分布形状の評価を行うために、再構成された各ボクセルは電気的に２値化してＨｉｇｈまたはＬｏｗのどちらかの状態にした。膜電位分布形状の評価では、２値化処理とモルフォロジー演算処理とを実行した。２値化処理では、図１８に示すように、再構成されたボクセルの最大輝度値と最小輝度値の平均を閾値（ここでは０．５）とし、各ボクセルがＨｉｇｈかＬｏｗかを判断した。

モルフォロジー演算処理では、周りの全てのピクセルがＬｏｗであるが中心１ピクセルだけＨｉｇｈの状態のとき、このＨｉｇｈのピクセルをＬｏｗに塗りつぶす処理を行った。また、この逆として、周り全てがＨｉｇｈであるが中心１ピクセルだけＬｏｗの状態のとき、このＬｏｗのピクセルをＨｉｇｈに塗りつぶす処理を行った。このようなモルフォロジー演算処理を行った理由としては、電気生理学的に１ボクセルの領域（６２５［μｍ］立方）のみ異なる状態にあることは考え難いことから、明らか誤差と見なして取り除くためである。

図１４Ａ〜図１４Ｆに示した６つの膜電位分布モデルの再構成誤差ボクセル数の平均（Ｅｒｒｏｒ ｖｏｘｅｌｓ）を縦軸とし、励起光の照射点数（Ｎｕｍｂｅｒ ｏｆ ｉｒｒａｄｉａｔｉｏｎ ｐｏｉｎｔｓ）を横軸としてグラフを作成したところ、図１９Ａ、図１９Ｂおよび図１９Ｃに示すような結果が得られた。再構成誤差ボクセル数は、図１４Ａ〜図１４Ｆに示した６つの膜電位分布モデルを２値化して正解データを得、再構成された各ボクセルと、正解データとを比較して、正解データと異なっているボクセルを再構成誤差ボクセルとしてカウントした。

図１９Ａ、図１９Ｂおよび図１９Ｃでは、再構成に用いる蛍光分布画像の画像枚数が５枚のときに得られた結果を「５」と表記し、１０枚のときに得られた結果を「１０」と表記し、２０枚のときに得られた結果を「２０」と表記した。計測回数（再構成に用いる蛍光分布画像の画像枚数）が増えるに従って、再構成誤差ボクセル数が減少する傾向を示したが、蛍光分布画像の枚数が少ない５枚のときでも、再構成誤差ボクセル数は十分に小さいことが確認できた。また、励起光の照射点数が増えるに従って、再構成誤差ボクセル数が減少する傾向を示した。

次に、図１４Ａ〜図１４Ｆに示した各膜電位分布モデルについて、計測回数１０回（すなわち、再構成に用いる蛍光分布画像１０枚）のときの再構成誤差ボクセル数を、ＤＣＴ、Ｈａａｒ １．３、およびＨａａｒ ３．３の変換毎にまとめたところ、図２０Ａ〜図２０Ｆに示すような結果が得られた。各変換の再構成誤差ボクセル数は、膜電位分布モデルによって大きく変化することが確認できた。また、励起光の照射点数においても、再構成誤差ボクセル数が最も小さくなる照射点数が膜電位分布モデルによって異なることが分かった。

次に、図１４Ａ〜図１４Ｆに示した各膜電位分布モデルについて、ＤＣＴを用いた再構成結果について調べたところ、図２１Ａ〜図２１Ｆに示すような結果が得られた。計測輝度値に応じてランダムに変化するショットノイズを蛍光分布画像に付加した。ショットノイズは１［％］とした。図２１Ａは、図１４Ａの膜電位分布モデルに照射するパターン光の照射点数を９０点とし、蛍光分布画像を得る計測回数を１回としたときの再構成結果を示す。図２１Ｂは、図１４Ｂの膜電位分布モデルに照射するパターン光の照射点数を５０点とし、蛍光分布画像を得る計測回数を１回としたときの再構成結果を示す。

図２１Ｃは、図１４Ｃの膜電位分布モデルに照射するパターン光の照射点数を１０００点とし、蛍光分布画像を得る計測回数を５回としたときの再構成結果を示す。図２１Ｄは、図１４Ｄの膜電位分布モデルに照射するパターン光の照射点数を１０００点とし、蛍光分布画像を得る計測回数を５回としたときの再構成結果を示す。図２１Ｅは、図１４Ｅの膜電位分布モデルに照射するパターン光の照射点数を９００点とし、蛍光分布画像を得る計測回数を５回としたときの再構成結果を示す。図２１Ｆは、図１４Ｆの膜電位分布モデルに照射するパターン光の照射点数を２０点とし、蛍光分布画像を得る計測回数を１回としたときの再構成結果を示す。

図２１Ａ〜図２１Ｆから、５枚の蛍光分布画像からでも誤差２５［％］以下の再構成結果が得られ、少ない枚数の蛍光分布画像からでも膜電位分布Ｖを求められることが確認できた。また、再構成結果が得られるまでの計測速度についても、蛍光分布画像の枚数を減らせた分、処理負担を軽減し得、膜電位分布を高速に計測できた。

１ 計測システム
２ 光学機器
３ 撮像機器
４ 演算処理装置
３３ 記憶部
３６ 再構成演算部
３７ 膜電位分布演算部
３９ 蛍光輝度計測部

Claims (6)

1. 光源から照射された励起光を基に、所定の照射点分布パターンのパターン光を生成し、前記パターン光を組織表面に照射する光学機器と、
   前記組織表面における蛍光分布を撮像して該組織表面の蛍光分布画像を取得する撮像機器と、
   前記蛍光分布画像に基づいて組織の膜電位分布Ｖを求める演算処理装置と、を備え、
   前記演算処理装置は、
   前記励起光および蛍光の強度分布関係をシミュレーションから求めた計測感度行列Ｊと、前記計測感度行列Ｊをスパース行列αにするための空間変換行列Ψと、を記憶する記憶部と、
   前記計測感度行列Ｊおよび前記空間変換行列Ψを基底行列として扱い、前記蛍光分布画像から求めた蛍光輝度値から前記スパース行列αを再構成する再構成演算部と、
   前記再構成演算部で再構成した前記スパース行列αを、前記空間変換行列Ψの逆変換によって前記膜電位分布Ｖに変換する膜電位分布演算部と、を備える
   ことを特徴とする計測システム。
2. 前記光学機器は、
   前記照射点分布パターンをランダムに変えながら前記パターン光を前記組織表面に照射してゆき、
   前記撮像機器は、
   前記照射点分布パターン毎に前記蛍光分布画像を取得してゆき、
   前記演算処理装置は、
   前記照射点分布パターン毎に得られた各前記蛍光分布画像の前記蛍光輝度値からなる蛍光輝度行列Ｙを求める
   ことを特徴とする請求項１に記載の計測システム。
3. 前記光学機器には、
   複数のマイクロミラーがアレイ状に配置され、各前記マイクロミラーがそれぞれ独立して動き、前記照射点分布パターンをランダムに変えるデジタルマイクロデバイスが設けられている
   ことを特徴とする請求項１または２に記載の計測システム。
4. 前記記憶部には、前記空間変換行列Ψとして、離散コサイン変換行列、離散ウェーブレット変換行列、のうち少なくとも１つ以上が記憶されている
   ことを特徴とする請求項１〜３のいずれか１項に記載の計測システム。
5. 前記記憶部には、前記離散ウェーブレット変換行列として、Ｈａａｒウェーブレット、Ｄａｕｂｅｃｈｉｅｓウェーブレット、Ｃｏｉｆｌｅｔウェーブレットの各行列のうち少なくともいずれか１つ以上が記憶されている
   ことを特徴とする請求項４に記載の計測システム。
6. 前記再構成演算部は、
   ｌ_１ノルム最小化法により前記スパース行列αを再構成する
   ことを特徴とする請求項１〜５のいずれか１項に記載の計測システム。

Images