(57)【特許請求の範囲】
【請求項１】
撮像被写体内へ伝達された電磁エネルギーに対する機械的波動応答を表す熱音響信号を提供する工程と、
該熱音響信号に基づいて被写体の局部的エネルギー吸収を表現するエネルギーデポジション画像を再構成する工程と、
該エネルギーデポジション画像を被写体における局部的吸収係数の分布を表現する定量的吸収画像及び少なくとも一の別な画像成分に分解する工程と、
を含んでなることを特徴とする被写体の熱音響撮像方法。
【請求項２】
請求項１に記載の撮像方法であって、
前記熱音響信号が光パルスの伝達に対する機械的波動応答を表すことを特徴とする撮像方法。
【請求項３】
請求項２に記載の撮像方法であって、
熱音響信号が光パルスの伝達に対する機械的波動応答を表し、該光パルスは少なくとも一の照射波長及び少なくとも一の偏光状態を有してなることを特徴とする撮像方法。
【請求項４】
請求項３に記載の撮像方法であって、
前記エネルギーデポジション画像は、前記少なくとも一の照射波長及び／又は前記少なくとも一の光偏光において取得された熱音響信号に基づいて再構成され、
複数のエネルギーデポジション画像が、異なる照射波長及び／又は異なる光偏光において取得された熱音響信号に基づいて個別に再構成され、又は
エネルギーデポジション画像は、異なる照射波長及び／又は異なる光偏光において取得された熱音響信号の組み合わせに基づいて再構成され、該組み合わせは、信号の線形又は非線形組み合わせを含んでなることを特徴とする撮像方法。
【請求項５】
請求項１〜４のいずれか一に記載の撮像方法であって、さらに
アナログ又はデジタル信号を処理する操作により熱音響信号を前処理する工程を含んでなる撮像方法。
【請求項６】
請求項１〜５のいずれか一に記載の撮像方法であって、前記エネルギーデポジション画像を再構成する工程は、
熱音響信号伝播の順モデル行列を提供する工程と、
反転されたモデル行列の手法を用いた熱音響信号を行列ベースで音響反転する工程と、
を含んでなることを特徴とする撮像方法。
【請求項７】
請求項６に記載の撮像方法であって、
前記反転されたモデル行列の手法は、疑似逆又は特異値分解、最小２乗（ＬＳＱＲ）法、及び代数的再構成法（ＡＲＴ）のうち少なくとも一を含んでなることを特徴とする撮像方法。
【請求項８】
請求項７に記載の撮像方法であって、
モデル行列は、熱音響信号を収集する検出装置の幾何学的形態に依存し、及び／又は
モデル行列は、被写体の音響的非均質性、撮像被写体内の伝達された電磁エネルギーの幾何学的形態、熱音響信号を収集するために使用された検出装置の周波数応答、分解工程の結果のうち少なくとも一に依存して調整されることを特徴とする撮像方法。
【請求項９】
請求項１〜８のいずれか一に記載の撮像方法であって、エネルギーデポジション画像を再構成する工程は、
反復伝播モデルをベースにしたシミュレーションを含んでなることを特徴とする撮像方法。
【請求項１０】
請求項１〜９のいずれか一に記載の撮像方法であって、エネルギーデポジション画像を分解する工程は、
エネルギーデポジション画像の対数表現を提供する工程と、
対数表現の加数として、定量的吸収画像及び少なくとも一の別な画像成分を計算する工程と、
を含んでなることを特徴とする撮像方法。
【請求項１１】
請求項１〜１０のいずれか一に記載の撮像方法であって、エネルギーデポジション画像を分解する工程は、
補完ライブラリーにおける拡張分として定量的吸収画像及び少なくとも一の別な画像成分を提供する工程と、
を含んでなることを特徴とする撮像方法。
【請求項１２】
請求項１１に記載の撮像方法であって、
補完ライブラリーにおける膨張係数は最小化アルゴリズムを用いて得られることを特徴とする撮像方法。
【請求項１３】
請求項１〜１２のいずれか一に記載の撮像方法であって、
前記少なくとも一の別な画像成分は、励起光フルエンスの分布、光学的散乱係数の分布、音速の分布、局部的偏光状態を表現する画像、電場又は磁場の分布、比誘電率の分布、熱弾性膨張係数の分布、及び熱容量の分布のうち少なくとも一を含んでなることを特徴とする撮像方法。
【請求項１４】
請求項１〜１３のいずれか一に記載の撮像方法であって、
前記少なくとも一の別な画像成分は、被写体の局部状態を表現する画像を含んでなることを特徴とする撮像方法。
【請求項１５】
請求項１４に記載の撮像方法であって、
前記局部状態は、被写体の生理学的状態を含んでなることを特徴とする撮像方法。
【請求項１６】
請求項１〜１５のいずれか一に記載の撮像方法であって、
少なくとも一の別な画像成分は、被写体内の組織バイオマーカーの空間的分布を表現する画像を含んでなることを特徴とする撮像方法。
【請求項１７】
請求項１６に記載の撮像方法であって、
前記バイオマーカーは、外因的に投与された分子剤、蛍光タンパク質、及びクロモジェニックアッセイのうち一又は複数を含んでなることを特徴とする撮像方法。
【請求項１８】
撮像被写体へ伝達された電磁エネルギーに対する機械的波動応答を表す熱音響信号を提供する工程と、
熱音響信号に基づいて被写体内の局部的エネルギーデポジションを表現するエネルギーデポジション画像を再構成する工程であって、反転されたモデル行列を用いた熱音響信号の行列ベースの音響反転を含む再構成工程と、
を含んでなることを特徴とする熱音響撮像方法。
【請求項１９】
請求項１８に記載の撮像方法であって、
モデル行列は、熱音響信号を収集する幾何学的形態に依存し、及び／又は
モデル行列は、被写体の音響的非均質性、撮像被写体に対する光パルス照射に関する幾何学的プロファイル、熱音響信号を収集するために用いた検出装置の周波数応答、及び分解工程の結果のうち少なくとも一に依存して調整されることを特徴とする撮像方法。
【請求項２０】
請求項１〜１９のいずれか一に記載の撮像方法であって、さらに
撮像被写体に電磁エネルギーを照射する工程と、
該照射に応答して生成される熱音響信号を記録する工程と、
を含んでなることを特徴とする撮像方法。
【請求項２１】
被写体を熱音響撮像するための撮像装置であって、
撮像被写体内へ伝達される電磁エネルギーに対する機械的波動応答を表す熱音響信号を収集するように構成された信号収集装置と、
収集された熱音響信号に基づいて被写体内の局部的エネルギーデポジションを表現するエネルギーデポジション画像を再構成するようになされたエネルギーデポジション画像を再構成する装置と、
エネルギーデポジション画像を、被写体内の局部的吸収係数の分布を表現する定量的吸収画像内と少なくとも一の別な画像成分に分解するよう構成された吸収画像分解装置と、
を備えてなることを特徴とする撮像装置。
【請求項２２】
請求項２１に記載の撮像装置であって、前記信号収集装置は、
電磁エネルギー源と、
少なくとも一の音響検出要素を含む音響信号記録装置と、
を備えてなることを特徴とする撮像装置。
【請求項２３】
請求項２１又は２２に記載の撮像装置であって、信号収集装置は,
撮像被写体を液状媒体から分離するメンブレンを備え、少なくとも一の音響検出要素が配置されてなることを特徴とする撮像装置。
【請求項２４】
請求項２３に記載の撮像装置であって、
前記メンブレンは、前記液状媒体に音響的にマッチした光学的透明材料で作製されることを特徴とする撮像装置。
【請求項２５】
請求項２１〜２４のいずれか一に記載の撮像装置であって、エネルギーデポジション画像の再構成装置は、
反転されたモデル行列を用いて熱音響信号に行列ベースの音響反転を受けさせるように構成された反転回路を備えてなることを特徴とする撮像装置。
【請求項２６】
請求項２５に記載の撮像装置であって、前記反転回路は、
モデル行列及び／又はその反転を記憶する行列記憶装置、及び／又は
被写体の音響的非均質性、撮像被写体に対する光パルス照射に関する幾何学的プロファイル、熱音響信号を収集するために用いた検出装置の周波数応答、及び分解装置の出力のうち少なくとも一に依存してモデル行列を調整するよう構成された調整装置と、
を備えてなることを特徴とする撮像装置。
【請求項２７】
請求項２６に記載の撮像装置であって、反転回路は、
反復伝播モデルをベースにしたシミュレーションを実行するようになされた反復回路を備えてなることを特徴とする撮像装置。
【請求項２８】
請求項２１〜２７のいずれか一に記載の撮像装置であって、前記吸収画像分解装置は、
エネルギーデポジション画像の対数表現を提供し、
該対数表現の加数として、定量的吸収画像及び少なくとも一の別な画像成分を計算するように構成されてなることを特徴とする撮像装置。
【請求項２９】
請求項２１〜２８のいずれか一に記載の撮像装置であって、
吸収画像分解装置は、補完ライブラリーにおける拡張分として定量的吸収画像及び少なくとも一の別な画像成分を提供することを特徴とする撮像装置。
【請求項３０】
被写体を光音響撮像するための撮像装置であって、
熱音響信号に基づいて被写体内における局部的エネルギーデポジションを表現するエネルギーデポジション画像を再構成するようなされたエネルギーデポジション画像再構成装置を備え、該エネルギーデポジション画像再構成装置は、反転されたモデル行列を用いて熱音響信号に行列ベースの音響反転を受けさせるようになされた反転回路を備えてなることを特徴とする撮像装置。
【請求項３１】
請求項３０に記載の撮像装置であって、前記反転回路は、
反転されたモデル行列を記憶する行列記憶装置、及び／又は
被写体の音響的非均質性、撮像被写体に対する光パルス照射、熱音響信号を収集するために用いた検出装置の周波数応答、及び分解装置の出力のうち少なくとも一に依存してモデル行列を調整するように構成された調整装置を備えてなることを特徴とする撮像装置。
【請求項３２】
請求項２１〜３１のいずれか一に記載の撮像装置であって、前記電磁エネルギー源は、
撮像被写体に電磁パルスを照射するように構成された照射装置を備えてなることを特徴とする撮像装置。
【請求項３３】
請求項２１〜３２のいずれか一に記載の撮像装置であって、前記音響信号記録装置は、
検出器アレイ及び／又は単一の検出器を備えてなることを特徴とする撮像装置。
【請求項３４】
請求項２１〜３３のいずれか一に記載の撮像装置であって、
音響信号記録装置及び被写体は互いに対して可動であることを特徴とする撮像装置。
【請求項３５】
請求項１〜２０の少なくとも一に係る方法を実施するプログラムを格納したコンピュータで読み取り可能な記録媒体に含まれるコンピュータプログラム。
【請求項３６】
請求項１〜２０の少なくとも一に係る方法を実施するためのプログラム命令を格納したコンピュータで読み取り可能な記憶媒体を備える装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、被写体の熱音響画像（thermoacoustic images）を生成する方法に関する。この方法により撮像された被写体内に伝達された電磁エネルギーに対する、機械的波動応答から得られた熱音響信号に基づいて、被写体内における局部的な吸収マップを表現する定量的画像が再構成される。さらに本発明は、上記の方法を実施するよう構成された熱音響撮像装置に関する。本発明の用途として、例えば医用撮像又は材料検査に利用できる。特に本発明は、非侵襲による小動物の撮像や臨床分野での撮像において、バイオマーカーの濃縮体を正確に抽出する用途で利用できる。例えば、血管樹、腫瘍の血管新生、血液酸化作用等の特性の決定に利用したり、染料、吸光性ナノ粒子、蛍光色素等、多様なターゲット造影剤を用いた分子レベルの撮像の研究にも利用できる。
【背景技術】
【０００２】
熱音響撮像は、組織内の局部的エネルギーデポジション（local energy deposition）を高解像度でマッピングするために、急速に発達している非侵襲的なモダリティ（modality）である。特に光音響的トモグラフィ（optoacoustic (or photoacoustic) tomography）は、組織における局部的な光吸収係数のマップを再構成し、さらに再構成後には内的、外的に投与された特定のバイオマーカー及びプローブの濃縮体と関連付けることができる。撮像は、被写体又は関心領域（region of interest：ROI）に高出力の短い電磁パルス、特にレーザーパルスを照射し、内部において瞬間的な温度上昇及び熱膨張を発生させることにより行われる。この結果、被写体の境界へ向かって伝播する高周波の音場が形成され、この後、音場を記録できる。このようにして、被写体周囲のトモグラフィ情報を収集し、光音響反転アルゴリズムを用いることにより、被写体内における局部的エネルギーデポジションを表現した画像を再構成できる。
【０００３】
光音響学における画像形成は複数の方法で実行できるが、これらはすべて、光音響応答のＲａｗ測定値から画像を形成するためにいくつかの数学的反転を必要とする。再構成画像は、撮像された組織における基礎的な吸収特性に関する定量的マップを表現すべきものであり、続いてこれを多様な組織の発色団及びバイオマーカーの分布に関連付けできるのが理想的である。しかしながら、大抵の場合、検出された音響信号は、撮像被写体における局部的エネルギーデポジション全体を表現するものになりがちである。したがって、最初の光音響画像は、関連情報を抽出するために更なる解析を必要とする。一だけ明らかに困難なことは、可視及び近赤外スペクトルにおいて、共通の生物（学的）組織が非常に非均質なものであることである。これはまた、散乱及び減衰により、撮像被写体内において励起光を非常に不均一な分布にしてしまう。その上、現在使用されている光音響反転モデルにおける不正確さが画像の定量化をさらに阻害している。
【０００４】
さらに詳述すれば、現在の光音響反転の技術は、音場で付与されたエネルギー画像を生成するが、実験的及び物理的伝播の要因を正確に説明するための適切な方法がないため、定量化の不正確さ、空間解像度の低さ、又は柔軟性の無さという欠点がある。反転及び画像形成のために現存するアルゴリズムを２つのグループに分けることができるが、一方は光音響的波動方程式の閉型解析解に基づいたものであり、他方は数値計算に基づいたものである。
【０００５】
上記第１グループのアルゴリズムは逆投影アルゴリズム（back-projection algorithms）であり、光音響画像の再構成は、線形のラドン変換に類似した測定信号の積分として求められる。これらの再構成法は通常、正確な光音響方程式へのいくつもの近似値に基づいているので、再構成において実質的なアーチファクトを作り出してしまう。その上、逆投影アルゴリズムは、計測器に依存する複数の要因を反転内に組み入れることができない。例えば、超音波検出器の周波数応答、又は逆投影アルゴリズムにおけるその有限サイズを直接考慮できないこともある。
【０００６】
上記第２グループのアルゴリズムはモデルベース反転であり、例えば有限要素法を用いて周波数又は時間のドメインのいずれかにおいて、音響伝播問題の解を数値的及び反復的に求める。モデルベース熱音響反転の体系は、以前にＨ．ジアン他、「ザ・ジャーナル・オブ・ジ・オプティカル・ソサエティ・オブ・アメリカ」、 Ａ２３巻、２００６年、８７８〜８８８頁）（非特許文献１）及びＧ．パルトーフ他、「ザ・ジャーナル・オブ・ジ・アクースティカル・ソサエティ・オブ・アメリカ」、１１２（４）巻、２００２年、１５３６〜４４頁）（非特許文献２）が試みている。しかしながら、これらの方法に伴う演算の複雑さにより、求め得る解が限定され、実用的な実施が阻害されている。
【０００７】
殆どの現実的な撮像シナリオにおいて、ＲＯＩ全体を均一に照射することは実際的に不可能であるので、最初に形成された画像は、必要とされる定量化された吸収係数の値よりも、組織におけるエネルギーデポジションと吸収を組み合わせたマップを表現することになる。換言すれば、熱音響画像は、被写体内における光学的吸収と光フルエンス間の積で構成することになる。したがって、被写体深部のターゲットは、類似した光学的吸収を有していても照射された表面近くに位置するターゲットよりも弱いように見えることがある。特に容積事象の定量的再構成は、光音響反転アルゴリズム、画像正規化法、及び対応するシステムの導入のいずれにおいても限度があったため、従来不可能であった。これらの不正確さが原因となり、従来の（定性的）光音響撮像の適応可能分野は限定されていた。結果的に、光音響撮像に関してこれまで発表されたシステムは、画質及び画像定量化の性能を損ない、深部の撮像又は容積的な撮像が悪化する
【０００８】
さらに多重スペクトルの光音響画像、すなわち、いくつかの励起波長に対して得られた画像を解析するとき、光フルエンスのスペクトル変化が光音響画像を占有し、関心ターゲットの吸収スペクトルを曖昧にすることもある。これが多重スペクトルの光音響トモグラフィ（ＭＳＯＴ；Multi-Spectral Optoacoustic Tomography）の利用を制限している。これに関して、例えば国際特許出願番号第ＰＣＴ／ＥＰ２００８／００６１４２号明細書（特許文献１；なお本願の優先日時点では未公開）には、スペクトルに依存する吸収を利用した多様なターゲットの濃縮体をマッピングすることに関して記載されている。
【０００９】
また、吸収係数の抽出及び定量化の改良に関する方法が、これまでに提案されてきた。いくつかの手法は、被写体内における光分布を探索して修正するため、光伝播を左右する拡散方程式の解を求めることに基づいている（Ｂ．Ｔ．コックス他、「アプライド・オプティックス」、４５巻、１８６６〜１８７５頁、２００６年（非特許文献３）及びＢ．Ｔ．コックス他、「プロシーディングズ・オブ・ＳＰＩＥ」、６４３７巻、６４３７１Ｔ−１−１０、２００７年（非特許文献４））。仮説を立てた光分布の解が得られると、これを用いて光音響画像を正規化し、吸収係数を抽出する。更なる工程として、抽出された吸収は、被写体内の光分布を計算し直すために用いることができ、そして光音響画像を正規化し直すために用いる。このプロセスは、収斂に達するまで反復方式で繰り返すことができる。しかしながら、上記の方法は、組織の光学的特性に関する経験的な仮定と、他の実験的なパラメータとに依存している。したがって、これらの方法は、ロバストな定量化の正確さに関する能力を限定する、収斂の不安定性という欠点がある（ジェッツフェルナー他、「アプライド・フィジックス・レターズ」、９５（１）巻、２００９年（非特許文献５））。
【００１０】
光音響画像正規化法の主な欠点は、これらの方法が光拡散方程式のような仮説的な光伝播方程式に基づいてモデル化された光フルエンス分布に依存していることが挙げられる。しかしながら、これらの方程式は、媒体の光学的特性に関して事前に正確な知識を必要とする。また、これらの特性は一般に殆ど知られておらず、特に散乱のみならず吸収も生じる。多くの場合は、推定値や構造は、媒体の光学的パラメータで表現できる。しかしながら現実には、付与された光学的特性における僅かな誤差であっても、再構成時には大きな誤差となり得ることが判明している。さらにまた、反復的な自己修復手法を利用する場合であっても、モデル化の不正確さにより、正確な値に収斂することは殆ど不可能である。したがって、今日までインビボのデータ、すなわち光分布が空間的に不均質である媒体に関しては、これらの手法では信頼に足る性能は確立されていない。
【００１１】
さらにこれらの制限は、光音響技術に限ったことではない。似たような問題点は、例えば被写体に対してエネルギーを伝達するために、被写体にレーザーパルスを照射する方法に代えて、無線周波数のパルスを使用した他の熱音響撮像方法においても一般的に発生している。
【先行技術文献】
【特許文献】
【００１２】
【特許文献１】国際特許出願番号第ＰＣＴ／ＥＰ２００８／００６１４２号明細書
【非特許文献】
【００１３】
【非特許文献１】Ｈ．ジアン他、「ザ・ジャーナル・オブ・ジ・オプティカル・ソサエティ・オブ・アメリカ」、 Ａ２３巻、２００６年、８７８〜８８８頁（H. Jiang et al., "J. Opt. Soc. Am." vol. A 23, 2006, p. 878-888）
【非特許文献２】Ｇ．パルトーフ他、「ザ・ジャーナル・オブ・ジ・アコースティカル・ソサエティ・オブ・アメリカ」、１１２巻（４）、２００２年、１５３６〜４４頁（G. Paltauf et al., "The Journal of the Acoustical Society of America", vol. 112(4), 2002, p. 1536-44）
【非特許文献３】Ｂ．Ｔ．コックス他、「アプライド・オプティックス」、４５巻、１８６６〜１８７５頁、２００６年（B. T. Cox et al., "Applied Optics" vol. 45, p. 1866-1875, 2006）
【非特許文献４】Ｂ．Ｔ．コックス他、「プロシーディングズ・オブ・ＳＰＩＥ」、６４３７巻、６４３７１Ｔ−１−１０、２００７年（B. T. Cox et al, "Proc SPIE" 6437, 64371T-1-10, 2007）
【非特許文献５】ジェッツフェルナー他、「アプライド・フィジックス・レターズ」、 ９５（１）巻、２００９年（Jetzfellner et al., "Appl. Phys. Lett." 95(1), 2009）
【非特許文献６】Ｄ．ラザンスキー及びＶ．ヌツィアクリストス、「メディカル・フィジックス」、３４巻、２００７年、４２９３〜４３０１頁（D. Razansky, V. Ntziachristos, "Med.. Phys." Vol. 34, 2007, pp. 4293-4301.）
【非特許文献７】Ｓ．Ｓ．チェン他、「ＳＩＡＭレビュー」、４３巻、２００１年、１２９〜１５９頁（S. S. Chen et al.,"SIAM Review", Vol. 43, 2001, pp. 129-159.）
【発明の概要】
【発明が解決しようとする課題】
【００１４】
本発明の目的は、従来の定性的再構成法における欠点及び限界を回避できる熱音響撮像のための改良された方法及び好ましい幾何学的な実施形態を提供することにある。特に本発明の目的は、表面領域及び表面下の領域だけでなく、体積的に被写体全体において、正確さと再現性を向上させた定量的撮像ができる正確な熱音響撮像方法を提供することにある。さらにまた、本発明の目的は、従来の方法における欠点及び制限を解消可能な熱音響撮像装置を実現することにある。
【課題を解決するための手段及び発明の効果】
【００１５】
上記の目的は、独立請求項の特長をそれぞれ備えた熱音響撮像方法及び／又は熱音響撮像装置を用いて達成される。本発明の好ましい実施形態は、従属請求項に記載されている。
【００１６】
本発明の第１の態様によれば、被写体の熱音響撮像方法が記載されており、エネルギーデポジション画像が分解されて、被写体の局部的吸収係数の幾何学形態的な分布を表現する定量的吸収画像及び少なくとも一の別な画像成分となる。まず、撮像被写体内へ伝達された電磁エネルギー（励起エネルギー）に対する機械的波動応答を表す熱音響信号が提供される。熱音響信号の提供は、例えば、撮像装置で信号を収集又は記録する工程や、データ源から信号を伝達する工程を含む。次に、この熱音響信号に基づき、被写体内における局部的エネルギー吸収を表現する少なくとも一のエネルギーデポジション画像が再構成される。この再構成は、従来の解析法又は数値法（例えば反復伝播モデルベースシミュレーション（iterative propagation model-based simulation））、好ましくは、以下に概説する本発明の第２の態様に係る、行列をベースにした光音響反転を含んでもよい。最後に、エネルギーデポジション画像は分解されて、定量吸収画像及び少なくとも一の別な画像成分となる。
【００１７】
用いられた励起エネルギーの波長に応じて、局部的な吸収係数分布（吸収画像）は、例えば光学的吸収分布又は導電性分布を含む。前記少なくとも一の別な画像成分は、典型例として、被写体における励起エネルギーの分布を含む。これに代わって又はこれに加えて、少なくとも一の別な画像成分は、以下の一又は複数に関する画像を含む。すなわち、光学的散乱係数の分布、音速の分布、局部的偏光状態を表現する画像、電場又は磁場の分布、比誘電率の分布、熱弾性膨張係数の分布、熱容量の分布、及び特定バイオマーカーの分布である。
【００１８】
本発明の第１の態様における本質的な部分を表す画像分解は、例えば現実の全身撮像シナリオで発生するような、非常に不均質な照射条件下の熱音響トモグラフィ再構成に関する新規な定量化法を提供する。好適なことに本発明の画像分解法は、媒体内におけるフルエンスの非均一性を説明するために正確な光子伝播モデルの選択に依存しない。その代わり、この方法は、スパース（sparse）表現の分野からの技法を用いて画像を主要な成分に分解する。本発明の分解法は、媒体内における吸収係数及びエネルギー密度（フルエンス）、特に光学的エネルギー又は光子の密度に関する異なる空間的特徴に依存する。本発明者らは、フルエンスが平滑的で全体的であり、他方、吸収係数は、医用撮像又は物質科学における撮像のような大抵の撮像用途において、より局在しており、また往々にして鋭利である空間的過渡現象を有するという知見を得た。スパース表現をベースにした分解を用いることにより、これらの異なる特徴は、撮像被写体内において吸収係数及びエネルギーデポジション密度のいずれも熱音響画像から抽出するのに利用される。以前の方法とは対照的に、このような画像正規化法は、理論的な光輸送方程式（light transport equations）の解に基づいておらず、ロバストな性能を提供する。これは、照射の幾何学形態、被写体の光学的特性、及びその他未知の、又は緩慢に定義された実験パラメータに関して明示的な知識を必要としないからである。したがって、本発明の方法は、組織バイオマーカーに関する多重投影照射システムや多重スペクトル光音響トモグラフィ研究等、多様な複雑さを呈するトモグラフィの体系において実用的な導入に適している。
【００１９】
本発明の第２の態様によれば、被写体の熱音響撮像方法が提供される。この方法では、被写体内における局部的エネルギーデポジションを表現するエネルギーデポジション画像を再構成するために、被写体へ伝達される電磁エネルギーに対する機械的波動応答を表す熱音響信号が用いられる。この再構成工程は、順モデル行列を用いた熱音響信号における行列をベースにした光音響反転を含む。
【００２０】
本発明の第２の態様における本質的な部分を表す行列ベースの光音響反転は、定量的熱音響画像の再構成のために新規の半解析的なモデルベース反転法を提供する。逆投影ライクなアルゴリズムとは対照的に、本発明の方法は、逆問題に対する解析解に基づくものではない。その代わり、順問題を行列関係に再定式化し、その後に反転される。好適なことに、熱音響信号の生成及び伝播が線形の順問題を呈するので、この種の最適化問題は一つだけの最小値を有する。さらに、モデル行列はスパースであるので、効率的な反転アルゴリズムが使える。理想的なことに、反転は、アーチファクトの無い定量化された熱音響画像の再構成をもたらすことができる。他の利点として、同じ行列反転のアルゴリズムを使いながら、多くの追加的な線形効果をモデル行列に追加できる。したがって、順モデルは非常に柔軟性があり、例えば音響結合、音響分散、及び他の伝播態様のみならず、検出器の周波数応答及び空間的応答を考慮することにより、多様な線形の減衰効果及び分散効果をモデル化するのに利用できる。
【００２１】
本発明がもたらす半解析解（semi-analytical solution）は、区分平面音響源関数（piecewise planar acoustic-source functions）に関して正確である。したがって、この解は音響源関数の平面補間に関しても正確である。このことにより、順モデルに関連した正確さ及び演算速度が大いに改善される。この解を用いて、モデル行列を構築でき、測定された音響信号を熱音響画像と連結する。行列関係を反転させることにより、熱音響画像が回復される。行列要素は設定上の幾何学形態に依存するだけなので、実験が行われる前に計算できる。また行列関係の反転は、行列の擬似逆を計算することによって事前に実行できる。実験が行われると、反転は、事前に計算された擬似逆を、測定されたデータを含むベクトルと乗算することに帰着する。これは、原則的にリアルタイムで行える演算である。数値的及び実験的な評価が示すところによれば、本発明の再構成アルゴリズムは、逆投影に関する再構成アーチファクトの欠点がなく、撮像被写体のレーザーエネルギーデポジションに関しても正確かつ高解像度のマップを提供できる。
【００２２】
行列モデルベース手法（matrix model-based approach）は、逆投影アルゴリズムに比べ、明らかに複数の利点を有する。まず、近似した逆投影の定式化に関する画像アーチファクトを排除できる。すなわち、マイナスの吸収値が生成されることがなく、再構成された画像は、被写体内における真の光減衰及びエネルギーデポジションに対応する。ここで逆投影アルゴリズムによれば、大きなマイナスのオーバシュートを生成することにより、エッジ部及び高速な画像の変化を誤って強調するので、「見栄えの良い」ハイコントラストな画像を生成できることは明白である。しかしながら、その近似の定式化により、組織における実際のレーザーエネルギーデポジション及び強調した光学的吸収値に関して正しい定量的画像を再生することはできない。このような特性は、定量的撮像用途、すなわち分子的撮像の研究において特に重要であり、正しい吸収マップを得ることは非常に重要である。第２に、本発明の方法は、反転手順を変更することなく、より総括的な音響伝播モデルに対して前向きな解法の一般化を許容する。最後に重要なこととして、モデルベース反転は、シームレスに任意の検出幾何学形態に適用できる。従来の解析的な反転定式化（analytical inversion formulation）とは異なり、検出器の位置は特定の幾何学形態に限定されることはない。
【００２３】
本発明の第３の態様によれば、被写体を熱音響撮像する撮像装置は、信号収集装置と、再構成装置と、分解装置を備える。特にこの撮像装置は、前記第１及び第２の態様における方法の一を実施するよう構成されている。信号収集装置は、撮像被写体内へ伝達された電磁エネルギーに応答して熱音響信号を収集するよう構成された装置である。この信号収集装置は、熱音響信号を処理するようになされた再構成装置に接続されている。この再構成装置は、収集された熱音響信号に基づいて、被写体内の局部的エネルギーデポジションを表現するエネルギーデポジション画像を再構成するようになされている。本発明の上記第２の態様を実施するために、再構成装置は、反転されたモデル行列を使って熱音響信号が行列ベースの光音響反転を受けるようになされた反転回路を備えることが好ましい。さらに再構成装置は、分解装置に接続されており、この分解装置はエネルギーデポジション画像を分解して、被写体における局部的な吸収係数の分布を表現する定量的吸収画像と少なくとも一の別な画像成分にするよう構成されている。
【００２４】
本発明の第４の態様によれば、被写体を熱音響撮像するための撮像装置は、熱音響信号に基づいて、被写体内における局部的エネルギーデポジションを表現するエネルギーデポジション画像を再構成するようになされた再構成装置を備える。この再構成装置は、反転されたモデル行列を用いて熱音響信号が行列ベースの光音響反転を受けるようになされた反転回路を備える。
【００２５】
上述した本発明の態様に関して、本発明は、測定された音響信号から被写体内の吸収されたエネルギーデポジションのマッピングを達成し、関数ライブラリーにおけるマッピングをスパースに表現することによりこのマッピングから吸収係数を抽出するシステム及び方法を明らかにする。それぞれの態様は本発明の独立した主題を表しているが、最大限の撮像性能を達成するために、これらの態様を組み合わせることも好ましい。
【００２６】
好ましくは、伝達された電磁エネルギーは、光スペクトルにおいて光子エネルギーを備える。すなわち、被写体は光を照射される。この照射光は、少なくとも１ｎｍ、好ましくは少なくとも４００ｎｍ、より好ましくは少なくとも６５０ｎｍ、又は５０００ｎｍ以下、特に好ましくは１２００ｎｍ以下、一層好ましくは８５０ｎｍ以下の範囲の、少なくとも一の特徴的な波長を有している。このような実施形態のために、信号収集装置は、好ましくは可視又は近赤外のスペクトル内の光パルスを発するようになされたレーザー源装置を備えている。これらの光パルスは、好適には１μｓ以下、特に好ましくは５０ｎｓ以下の持続時間を有する。
【００２７】
一般に照射光は、単一の照射波長及び／又は単一の光偏光を備えている。したがって、エネルギーデポジション画像は、この単一の照射波長及び／又は単一の光偏光において獲得された熱音響信号に基づいて再構成できる。
【００２８】
あるいは、照射光が複数の異なる励起波長及び／又は複数の異なる光偏光を備えておれば、本発明の定量的撮像に関する情報内容を増加できることが好適である。この実施形態では、複数のエネルギーデポジション画像は、異なる照射波長及び／又は異なる光偏光／場偏向（light/field polarizations）において個別に得られた熱音響信号に基づいて個別に再構成できる。他の考え方として、エネルギーデポジション画像は、異なる照射波長及び／又は異なる光偏光／場偏向において得られた熱音響信号の組み合わせに基づいて再構成できる。熱音響信号の組み合わせは、線形又は非線形信号の重畳、例えば減算、加算、積、率等を含む。
【００２９】
あるいは、伝達された電磁エネルギーは、好ましくは１μｓ以下のパルス持続時間、最も好ましくは５０ｎｓ以下のパルス持続時間を有する、無線周波数又はマイクロ波のスペクトル領域におけるパルスを備えていてもよい。したがって、信号収集装置は、撮像被写体に無線周波数を照射するようになされた放射源を備えていてもよい。
【００３０】
本発明の他の実施形態によれば、熱音響信号は、画像の再構成前に前処理を受けることができる。好適なことに、アナログ信号又はデジタル信号の処理作業による前処理、例えばフィルタリング、積分、微分は、画像の再構成結果を改善し得る。
【００３１】
熱音響信号の行列をベースにした光音響反転を用いた本発明の画像再構成工程に関して、モデル行列の要素は、熱音響信号生成及び伝播に影響する線形効果に依存して選択できる。本発明の好ましい変形例によれば、モデル行列は、熱音響信号を収集する機械的波検出装置の幾何学形態に依存して構築される。これに代わって又はこれに加えて、モデル行列は、被写体の音響的不均質性、撮像被写体における伝達された電磁エネルギーの幾何学的プロファイル、熱音響信号を収集するために用いられた音波信号記録装置の周波数応答、及び分解工程の結果のうち、少なくとも一に依存して構築される。好適なことに、モデル行列要素は殆ど被写体の吸収特性とは独立である。
【００３２】
本発明の特に好ましい実施形態によれば、画像分解工程のために、再構成されたエネルギーデポジション画像の対数表現が提供される。したがって、吸収画像及び少なくとも一の別な画像成分の寄与は、別個の加数として得られる。続くエネルギーデポジション画像の分解工程は、エネルギーデポジション画像の対数表現の加数として、定量的吸収画像及び少なくとも一の別な画像成分の計算を含む。これらの工程を実行するため、本発明の撮像装置は、定量的吸収画像及び少なくとも一の別な画像成分を計算する対数回路及びプロセッサーを含む分解装置を備えている。
【００３３】
各々の加数は、補完ライブラリーにおける拡張分（expansion）として提供されるのが好ましい。定量的吸収画像は、第１のライブラリー（関数の集合体、例えば関数基準（function basis））における拡張分として表され、他方、少なくとも一の別な画像成分は、少なくとも一の別なライブラリーにおける拡張分として表される。ライブラリーは、第１のライブラリーを利用するときに必要とされるよりもかなり多数のライブラリー関数を、少なくとも一の別なライブラリーにおける定量的吸収画像の表現が必要とする場合に、補完ライブラリーとなる（逆もまた真である。すなわち、第１のライブラリーにおける少なくとも一の別な画像成分の表現は、少なくとも一の別なライブラリーを利用するときに必要とされるよりもかなり多数のライブラリー関数を必要とする。）好適なことに、補完ライブラリーにおける表現の膨張係数は、最小化アルゴリズムを用いて得られる。
【００３４】
本発明の格別な利点は、分解工程の結果として吸収画像から多様なタイプの画像成分が分離できることである。好ましくは、少なくとも一の別な画像成分は、励起光フルエンス分布画像及び／又は被写体における光学的散乱係数の分布を備える。この場合、被写体におけるエネルギーデポジションの散乱によるアーチファクトを回避できる。あるいは、少なくとも一の別な画像成分は、被写体における、局部的偏光状態（例えば、被写体の複屈折又は二色性の状態）、電場又は磁場の分布、音速の分布、比誘電率の分布、熱弾性膨張係数の分布、及び／又は熱容量の分布を含む。
【００３５】
少なくとも一の別な画像成分が被写体の局部状態、特に生物学的組織の生理学的状態を表現する画像を含む場合、医用撮像にとって格別な利点が得られる。生理学的状態は、例えば、疾患進行、酸化作用の状態、及び／又は代謝状態を表し得る。さらに別の考え方として、少なくとも一の別な画像成分は、組織バイオマーカー（例えば被写体における、外因的に投与された分子剤、蛍光タンパク質、機能化された粒子、及び／又はクロモジェニックアッセイ）の空間的分布を含み得る。
【００３６】
本発明の他の好ましい実施形態によれば、熱音響信号の提供は、信号収集（測定工程）を含む。撮像すべき被写体は、電磁エネルギー、例えばレーザーパルスを照射され、この照射に応答して生成される熱音響信号が記録される。好適なことに、信号収集は、従来の熱音響撮像、特に光音響撮像においてしかるべく知られている通りに実行できる。この目的のため、信号収集装置は、電磁エネルギー源と、少なくとも一の音響検出要素を含む検出器アレイ及び／又は単一の検出器を有する音響信号記録装置とを備える。好ましくは、撮像被写体を例えば水又は生理溶液のような液状媒体から分離する、音響的にマッチしたメンブレンが提供され、そして少なくとも一の音響検出要素が配置される。好適なことに、このメンブレンは、被写体を液状媒体の中に浸漬することなく撮像できるようにする。
【００３７】
好適なことに、本発明の撮像装置に、本発明の前記第１又は第２の態様に係る方法を実施するように構成されたハードウェア回路を設けられる。一例として、熱音響信号を行列ベースの光音響反転するための反転回路は、モデル行列及び／又はこれら行列のそれぞれ（又は両方）の反転バージョン又は切捨てバージョンを記憶し、これらの行列のうち最も強力な要素を収容する少なくとも一の行列記憶モジュールと、撮像条件（例えば、被写体の音響的非均質性、撮像被写体に対する光パルスの照射、熱音響信号及び分解装置の出力を収集するために用いられた音響信号記録装置の周波数応答及び／又は空間的応答のうち少なくとも一）に依存してモデル行列を調整するようになされた調整装置とを備え得る。さらに再構成装置は、反復伝播モデルベースシミュレーションを実施するようになされた反復回路を備えてもよい。
【００３８】
本発明の撮像装置に関する他の好ましい実施形態によれば、音響信号記録装置及び被写体は互いに対して可動である。熱音響信号を収集する工程、エネルギーデポジション画像を再構成する工程、及び画像を分解する工程は、信号記録装置及び被写体が互いに対して変動する幾何学的配向及び／又は配置で反復できる。
【００３９】
本発明の分解工程は、多重スペクトル光音響トモグラフィ（ＭＳＯＴ）法においても第一の解析工程として利用できる。ＭＳＯＴ法の実施に関する詳細は、国際特許出願番号第ＰＣＴ／ＥＰ２００８／００６１４２号明細書（特許文献１）を引用し、本明細書の開示に援用する。好ましい実施形態の一において、分解は、例えば和により正規化された異なる励起波長に対して得られた画像間の差に対して行われる。これにより、スペクトル分解中の光フルエンスのスペクトル依存に関して起こり得るマイナス効果がなくなる。
【００４０】
提起した本方法を効果的に実施するために、十分な量の情報が存在することを確実にするには、組織に対して最適に照射し組織からデータを収集するための適切な幾何学形態が必要とされる。例えば、音響信号を測定するのに用いる時間的及び空間的解像度は、再構成された画像における所望のピクセル／ボクセル解像度に対応することが好ましい。さらに、撮像対象になっている被写体のいかなる部分も、極力多数の音響測定のために極力均一に照射されることが好ましい。これが好ましいのは、高質な全身（体積的又は断面積的）画像を得るときに本方法の利点が明白になるからである。しかしながら、今日ではさらに普遍的に可能にされているように、相対的に均一な励起エネルギーデポジションを有する表面限定の撮像には必須でなくてもよい。したがって、正確な実施に求められる枢要な要素を説明するために、最適に本方法を実施できるシステムも提起されている。
【００４１】
従来の技術とは対照的に、本発明の方法は、例えば生物学的組織における吸収係数を定性的よりも定量的にマッピングすることを可能にしている。したがって、この方法は、インビボでの全身撮像の用途、特に小動物の撮像及び臨床的撮像における深部組織を扱う用途には必須のものである。後者の場合、撮像被写体における現実の光学的特性及び他の実験的パラメータが貧弱に決められているか、知られていないか、実験過程において変動さえもするようなとき、体積的な撮像、すなわち、（光輸送の平均自由経路長（mean free path length：ＭＦＰＬ）に比して）大きな試料の撮像に関連する様相に対して効果的に取り組む技術も提供することが必要である。
【００４２】
本発明の他の側面は、前記第１又は第２の態様に係る方法を実施するプログラムコードと共に、コンピュータで読み取り可能な記録媒体に記録されたコンピュータプログラムと、本発明の方法を実施するためのプログラム命令を収納しているコンピュータで読み取り可能な記録媒体を備えた装置である。
【００４３】
本発明に関する更なる詳細及び利点について、以下添付図面を参照しながら説明する。
【図面の簡単な説明】
【００４４】
【図１】本発明の実施形態に係る撮像方法の主要な工程を示す概観図である。
【図２】本発明の撮像装置に関する実施形態を示す概略図である。
【図３】本発明の撮像装置に関する実施形態を示す概略図である。
【図４】本発明の撮像装置に含まれる信号処理装置の主要な要素を示す概略図である。
【図５】ｘ−ｙ軸方向の再構成平面を並べて表示するのに使用されるグリッドを示すグラフ図である。
【図６】光音響画像における異なる要素と、補完ライブラリーにおけるこれらのスパース表現とを示すグラフ図である。
【図７】光音響画像における異なる要素と、補完ライブラリーにおけるこれらのスパース表現とを示すグラフ図である。
【図８】本発明の撮像装置に含まれるエネルギーデポジション再構成装置に関する他の変形例を示す概略図である。
【図９】本発明の撮像装置に含まれるエネルギーデポジション再構成装置に関する他の変形例を示す概略図である。
【図１０】本発明の撮像装置に含まれる吸収画像分解装置の主要な要素を示す概略図である。
【発明を実施するための形態】
【００４５】
本発明の好ましい実施形態について、以下では光音響トモグラフィの例を参照しながら説明する。なお本発明は、これらの実施形態に限定されるものではなく、光音響相互作用又は熱音響相互作用を用いた他の方法、例えば光音響トモグラフィ（ＰＡＴ：Photo-Acoustic Tomography）、多重スペクトル光音響トモグラフィ（ＭＳＯＴ：multispectral optoacoustic tomography）、 熱音響トモグラフィ（ＴＡＴ：thermoacoustic tomography）、又は近接無線周波数トモグラフィ（ＮＲＴ：Near-Field Radiofrequency Tomography）に対しても、同様に適用できる。
１．概観
【００４６】
光音響撮像を使って被写体、例えば人体又は動物体あるいはその一部について、吸収係数のマッピングは、概ね３つの工程に分けて説明できる。これらの工程に関して概略的な説明を図１に示す。第１の工程Ｓ１は、データ取得を含んでおり、高出力の短パルスレーザーを用いて撮像被写体を照射すること、及び形成された音響波を、被写体の境界近傍に位置する検出器で測定することを要する。次の２つの工程Ｓ２及びＳ３は、本明細書に開示する本発明の中核部分であり、被写体内における光学的吸収係数に関する２次元又は３次元マップを得るために測定データを処理する必要がある。第１の信号収集工程の直後に第２の工程Ｓ２及び第３の工程Ｓ３を実行することは、本発明における任意選択の特長である。またこれに代えて、本発明の画像再構成を、画像データ記憶装置から、及び／又はコンピュータネットワークを介して得られるデータに対して適用することもできる。
【００４７】
第２の工程Ｓ２において、モデルベース音響反転法は、音響波伝播を後戻りさせ、パルス化された光吸収の結果として被写体内に生成された光音響源を見つけ出すために用いられる。測定された音響場と光音響源との間の関係は線形であるので、数値計算を促進するためにこの関係を行列形状へ離散化することが好ましい。この行列は、音響伝播と検出を必要とするすべての線形効果を含んでもよい。例えば、音響検出器の周波数応答、音速の不均質性、又は音響減衰及び音響分散は、この行列関係において説明できる。行列が計算されると、行列関係を反転させ、被写体内における光音響源の分布（画像）を得るために多様な反転の体系を利用し得る。
【００４８】
光音響源の振幅は被写体内に吸収されたエネルギー量と比例している。したがって、光音響源をマッピングしても直接的に吸収係数をもたらすわけではなく、吸収係数及び光フルエンス間の積を出すのである。この結果、光音響源の画像は、よりよく照光されている被写体の領域の方へ偏倚される。第３の工程Ｓ３は、このような偏倚をなくして、光音響源の画像から吸収係数の画像を得ることを含む。この工程Ｓ３は、身体内の吸収係数及び光フルエンスという異なる空間特性を利用する本発明の分解法に基づいており、これらの空間特性は、これら２つの関数に対して異なるスパース表現を規定している。光拡散ゆえに、フルエンスは緩やかに変動する関数としてみなされてもよく、他方、吸収係数は鋭利な過渡現象を有してもよい。多重スペクトル測定において、第３の工程は、異なる励起波長に対して得られた光音響画像に関する正規化された差に適用されてもよい。
【００４９】
好ましい実施形態において、本方法は、２つのライブラリーの一体化に対してロガリズムの最もスパースな表現を計算することによって光音響画像を分解する。最もスパースな表現は、利用可能ないくつかのアルゴリズムの一を使って求めることができる。この表現において、吸収係数のライブラリーに含まれている全要素を用いて、定量化された吸収係数画像を効果的に抽出できる。
【００５０】
図２は、被写体１を熱音響撮像する撮像装置１００の実施形態を、概略的に示している。この撮像装置１００は、信号収集装置１０と、画像処理装置２０と、必要に応じて他の制御装置３０とを備える。これらの要素２０、３０は、共通のハードウェアユニット内に実装してもよい。信号収集装置１０は、少なくとも一のレーザー源１１のような、５０ｎｓ以下の継続時間を有するパルスを発する電磁エネルギー源と、少なくとも一の音響検出要素を有する音響信号記録装置１２（検出装置）とを備える。レーザー源１１は、多重スペクトルの情報が被写体１から、例えば組織バイオマーカーの明確なスペクトルシグネチャ（spectral signatures）を配置するために収集できるよう、異なる波長に同調できる。
【００５１】
信号収集装置１０に関する細部は、「メディカル・フィジックス」、３４巻、２００７年、４２９３〜４３０１頁（非特許文献６）においてＤ．ラザンスキー及びＶ．ヌツィアクリストスによって説明されているように実施できる。特にレーザー源１１は、可視スペクトル（４００〜７００ｎｍ）又は近赤外スペクトル（７００〜１４００ｎｍ）において、繰返し周波数が１０〜１００Ｈｚ範囲内であって、継続時間が１０ｎｓ未満のパルスを提供する、同調可能な光学パラメータ式発振器（ＯＰＯ：optical parametric oscillator）のレーザーを備えていてもよい。一例として、照射波長を７５０ｎｍに設定し、出力エネルギーを２５ｍＪ／パルス以下に保ち、出力した近赤外線ビームを試料平面においてビームの平均直径が約２ｃｍであるように拡張している。
【００５２】
被写体１におけるＲＯＩに対して極力均一に照射するため、レーザービームを拡張して、被写体の直径全体を対象にして２つ以上のビームに分岐し、異なる方向から被写体１に照射できるようにする。さらに、被写体１中への光の伝達は、例えば複数のアームを有する光ファイバ又はファイバ束で、又はミラー、ビームスプリッタ、拡散器等の光学的構体を介して自由ビームを送り込むこと等により実施できる。
【００５３】
音響信号記録装置１２は、例えばゲルのように音を伝達する物質又は液体を使って被写体の表面と音響的に接触した状態で配置される。被写体１は、図３に示すように液体から部分的に又は完全に分離できる。音響信号記録装置１２は、例えば単一の圧電形の点検出器を備える。３６０°の投影に亘って信号を収集するために、被写体１は、特定の角度インクリメント（例えば２°又は３°）にてステージ上で回転できる。他方で検出器は、回転の中心から、例えば４０ｍｍの距離に配置される。音響信号記録装置１２は、円筒状の走査形態で３次元の再構成のためにデータを取得するよう、上下軸に沿って所定数の点に亘って直線移動できる（Ｄ．ラザンスキー及びＶ．ヌツィアクリストスの上記文献を参照）。検出器の表面は、異なるＲＯＩ（線、面等）から光音響信号を優先的に収集できるよう、異なる合焦形状（例えば円筒状、球状）に択一的に似せることができる。典型的な光音響信号が超広帯域の性質を有するために、検出器は理想的に、例えば１〜１０ＭＨｚの範囲、より好ましくは２００ｋＨｚ〜３０ＭＨｚの範囲、最も好ましくは５０ｋＨｚ〜１００ＭＨｚの範囲という広い有効（−６ｄｂ）帯域幅を有するべきである。あるいは音響信号記録装置１２は、複数の投影角度から同時に情報を収集するために、複数の音響検出要素を有する検出器アレイを備える。光音響応答の検出は、光学的手段、例えば押圧の変化に敏感な干渉装置（例えばファブリーペローフィルム及び他の多層構造物）を用いることによって択一的に行える。
【００５４】
一例として、検出器の周波数応答の半値全幅（ＦＷＨＭ）は３．５ＭＨｚとできる。これは、約２００μｍの回折限界の空間的解像度に相当する。検出器の周波数応答は、再構成された画像の空間的解像度を制限するだけでなく、検出された信号の形状も変形させ得るので、余分な画像アーチファクトをもたらしてしまう。したがって、再構成プロセス中、周波数応答を考慮に入れる（下方を参照）。
【００５５】
音響信号記録装置１２で記録される、時間分解された熱音響信号（典型的には電圧信号）は、デジタル化でき、（ノイズ抑圧のために）平均化でき、定量的画像再構成を受け持つ画像処理装置２０の中央演算処理装置（ＣＰＵ）へ送信できる。
【００５６】
図３に示すように、レーザー源１１で照射された被写体１は、透明なメンブレン１３の上面に配置される。このメンブレン１３は、音の伝播も甚だしく変化させないような光学的に透明な材料で作製されており、例えばポリ塩化ビニル系プラスチゾル（ＰＶＣＰ）で作製された音響（学）的にマッチしたメンブレン、又はポリエチレンの薄膜である。このようにして、レーザー放射及び誘発された光音響応答のいずれも効率的な伝達が促進される。またメンブレン１３は、確実に、水と直接接触することなく撮像するため被写体１を首尾良く配置できるようにする。メンブレン１３は、水槽１４と接触した状態で配置され、この中に検出器アレイ１２が浸漬される。この検出器アレイ１２は、６４個の音響検出要素１５を有する湾曲したアレイである。このアレイは、被写体１を取り囲む例えば１８０°の弧を形成する。図３の例では、すべての音響検出要素１５が検出平面において円筒状に合焦され、２次元の検出平面（図面の平面）から信号を同時に収集できる。したがって、この装置はリアルタイムに２次元データを収集できる。
【００５７】
画像処理装置２０の主要な要素は、図４に示されている。収集されたＲａｗ信号を信号インターフェイス部２１において前処理した後、熱音響信号は、エネルギーデポジション画像を定量的に再構成する再構成装置２２へ伝達される（下記第２項を参照）。次にエネルギーデポジション画像は、エネルギーデポジション画像を定量的吸収画像及び少なくとも一の別な画像成分に分解する分解装置２３内で処理される（下記第３項、第４項を参照）。最後に、定量的吸収画像のデータは、インターフェイス２４を経由して出力され、その後、例えば他のデータ処理（例えばデータの記憶、画像データの記録、画像の表示、又はプリントアウトのような画像の出力等）が続く。
２．定量的熱音響トモグラフィのための半解析的モデルベース光音響反転
【００５８】
光音響トモグラフィの順問題は、結局、レーザー照射の瞬間的加熱及び弛緩効果によって生成される音場をモデル化することに帰着する。熱的な閉じ込めを条件として、照射された被写体の各部における温度上昇は、近傍領域における温度上昇に影響されない。すなわち、熱伝導は無視できる。このような条件は、通常、持続時間が１μｓ以下であるレーザーパルスに関して満たされ、この条件により、被写体内で生成された音響源は、吸収された光エネルギーに比例することが保証される。この条件下で、及び音響損失を無視して、音場に関する伝播方程式は、数１より求められる。
［式１］
【００５９】
【数１】
0005566456.tif 000002
【００６０】
式中、ｃ及びρは、組織内における音速及びその密度であり、Γは、グリュナイゼン・パラメータであり、ｈは、組織内に吸収された単位体積当たり及び単位時間当たりのエネルギー量である。Ｈは、その空間成分ｈ及び時間成分ｈtを乗算した積として表せる。本解析において、音響的に均質な媒体における音場の伝播は仮定できる。このような仮定のもとで、上記の方程式は、式２の形で表される。
［式２］
【００６１】
【数２】
0005566456.tif 000003
【００６２】
殆どの実用的な場合、光パルスの持続時間は、デルタ関数、すなわちｈt（ｔ）＝δ（ｔ）で十分近似できるほど短いものである。この場合はすなわち、レーザーパルスによる加熱が瞬間的なものであると仮定し、また媒体が音響的に均質であると仮定して、ポアソン型の積分で解析解が求められる。座標ｘ及び時間ｔにおいて測定された音圧は、式３を具備する。
［式３］
【００６３】
【数３】
0005566456.tif 000004
【００６４】
ここでＲ＝│ｘ−ｘ’｜であり、νは、媒体内における音速であり、ｈ（ｘ）は、光音響画像である。測定を３Ｄで行うと、球体に対して式３の積分が行われる。また測定を２Ｄで行うと、円に対して式３の積分が行われる。以下の説明では、２Ｄの場合を説明しているが、３Ｄの場合に対しても同じ概念で一般化できる。特に、異なる平面において発生する光音響信号を収集することにより、解を容易に３Ｄに一般化できる。
【００６５】
２次元（２Ｄ）の幾何学的形態においては、すべての発生源が平面内にあり、円に対して積分する。所与のセンサ位置ｒ＝（ｘ0，ｙ0）に関し、式３における積分は、式４として明示的に書き直せる。
［式４］
【００６６】
【数４】
0005566456.tif 000005
【００６７】
２Ｄ及び３Ｄの幾何学的形態のいずれにおいても、式４の計算には、積分計算すべき表面とグリッドが一致していないことに起因する数学的な難しさが発生する。積分の導関数を計算しても、数的な難しさを悪化させるだけである。典型的には、式４を正確に計算するためには、グリッドの解像度はｈ（ｒ）を正確に表すために他の場合に必要とされるものよりも、相当高いものでなければならない。これは結果的に、順問題の解を計算する際の非効率さとなってしまう。
【００６８】
本発明によれば、ｈ（ｘ）を散在的に表すようにグリッドが定義されており、グリッドの座標はｘiで記される。式３において積分計算するためには、グリッド上にはない座標に対して閉型の近似値ｈ（ｘ）を使用する。そのような近似値に関して唯一必要とされることとして、任意所与の座標ｘ’に関して、値ｈ（ｘ’）は、メッシュ状座標における値ｈ（ｘ）の線形重畳でなくてはならないことである。
［式５］
【００６９】
【数５】
0005566456.tif 000006
【００７０】
メッシュポイントがあらかじめ決められているので、関数ｆi（ｘn，．．．，ｘm）の値は、測定データとは独立である。式５を式３に代入すれば、測定された音圧及び光音響画像間の線形接続が得られる。
［式６］
【００７１】
【数６】
0005566456.tif 000007
【００７２】
式中、ｇi（ｘ，ｔ，ｘk，．．．，ｘl）の値は、関数ｆi（ｘn，．．．，ｘm，ｘ’）を被積分関数として、式３における積分を解くことによって算出される。
【００７３】
適切な補間式は、次のようにして見い出せる。これらの補間式は、式３における積分に対する解析解を求めるようなものでなくてはならない。さらに、式３における導関数の演算子のため、補間式は、微分可能又は区分的微分可能でなければならない。好ましくは、得られる関数ｈr（ｒ）が区分的平面である線形補間を用いる。
【００７４】
補間は、図５に例示されたグリッド点上に頂点を有する直角三角形にｘ−ｙ再構成平面を並べて表示することによって行われる。各座標（ｘi，ｙi）に関して、ｚ軸上の上昇値としてｈr（ｒi）が割り当てられ、すなわち、ｚi＝ｈr（ｒi）である。したがって、各三角形は、その頂点（ｘl，ｙl，ｚl）、（ｘn，ｙn，ｚn）、及び（ｘm，ｙm，ｚm）である３つの座標セットで表せる。したがって、各三角形内におけるｈr（ｒ）の補間値は、式７を介して計算された平面上昇の値とされる。
［式７］
【００７５】
【数７】
0005566456.tif 000008
【００７６】
ここで、以下の条件が満たされる必要がある。
［式８］
【００７７】
【数８】
0005566456.tif 000009
【００７８】
演算｜・｜は、行列の行列式である。式７を式４に代入して、式９が得られる。
［式９］
【００７９】
【数９】
0005566456.tif 000010
【００８０】
次に、式５における積分は、各三角形内において解析的に解を得ることができ、ｚl、ｚn、及びｚmの線形関数が得られる。続いて、ｚl、ｚn、及びｚmの係数を用いてｇi（ｒ，ｔ）を計算する。
【００８１】
図６に示す例では、各グリッド点が３つの座標を割り当てられている。各直角三角形は、画像が線形補間されている単位セルを表す。単位セルの例は点線で示されており、ＯＡＴの値は、その頂点（ｘl，ｙl，ｚl）、（ｘn，ｙn，ｚn）、及び（ｘm，ｙm，ｚm）である３Ｄ座標によってのみ決定される。破線で示す曲線は、式９の積分が計算される特定の弧を表し、グリッドの太い黒点は、積分計算する際に必要であり式６における非ゼロの要素に対応するグリッド点である。
【００８２】
次に、音響信号が得られるときの空間及び時間座標｛（ｘj，ｔj）｝が定義される。ｐj＝ｐ（ｘj，ｔj）及びｈi＝ｈ（ｘi）を表示して、以下の行列方程式が得られる。
［式１０］
【００８３】
【数１０】
0005566456.tif 000011
【００８４】
式中、ｐは、要素がｐjであるベクトルであり、ｈは、要素がｈiであるベクトルである。
【００８５】
行列Ｍにおける要素は、画像のグリッド、すなわち測定値の空間及び時間座標によって定義され、また用いられている近似の種類に基づいて定義される。特に行列Ｍにおける要素は、ｇi（ｒ，ｔ）の対応値である。したがって、所定の光音響撮像システムに関して、行列Ｍは、撮像被写体に依存しないので、実験を行う前にあらかじめ計算できる。
【００８６】
音響検出器の周波数応答がおおよそ判明している場合、この周波数応答を順モデル内に組み込むことができる。周波数応答は、離散化されてベクトル内に保存される。そして対角行列Ｍsが構築され、その対角線要素はベクトル内の要素に対応する。Ｆは、１次元の離散フーリエ変換の行列として表示される。すなわち、Ｆｐは、ｐの離散フーリエ変換である。したがって、測定された信号ｐsのスペクトルは、ｐs＝ＭsＦｐから求められ、式１１を満たす。
［式１１］
【００８７】
【数１１】
0005566456.tif 000012
【００８８】
また光音響反転法は、システム内に現れる付加的な効果も備える。すなわち撮像被写体が音響的に非均質であって、被写体内での音速変動が判明しておれば、これらを行列方程式に組み込むこともできる。これは、式３の代わりに式１２に関して行列Ｍを計算し直すことによって達成できる。
［式１２］
【００８９】
【数１２】
0005566456.tif 000013
【００９０】
ここでＳ（ｔ）は、各点が検出器まで同じ飛行時間ｔを有する曲線又は表面である。
【００９１】
行列関係に組み込める他の効果は、異なる検出器位置に対する照射プロファイルの変位に由来し得る。このような変位は、被写体が光源及び検出器に対して回転される場合に出現する。この場合、式１２の積分を式１３のように一般化する必要がある。
［式１３］
【００９２】
【数１３】
0005566456.tif 000014
【００９３】
ここでＵ（ｘ’，θ）は、角度θにおける投影に関する被写体内のフルエンスへの推定である。
【００９４】
光音響画像ｈを回復させるために、式１１が反転される。この反転を行うにはいくつかの方法が存在する。最も一般的な方法は、最小二乗誤差の手法を用いることである。この方法で、以下の二乗誤差を最小化するベクトルｈ’が見つかる。
［式１４］
【００９５】
【数１４】
0005566456.tif 000015
【００９６】
最小化は、勾配降下、共役勾配、ＬＳＱＲ等の標準的な演算方法を用いて実行できる。あるいは、行列ＭsＦＭｈの疑似逆を計算することで最小化を行ってもよい。疑似逆を用いる利点は、これがシステムの特徴に依存するだけであって、測定された被写体には依存しないことである。この場合、再構成は、ｐsに所定の行列を掛け合わせることに帰着する。モデル行列はスパースであり、すなわち、その要素の殆どがゼロである。このような特性は、反転過程の改善や行列の記憶に利用できる。
【００９７】
測定データが正確な再構成を果たすのに不十分なものであれば、例えばいくつかの投影が得られないような場合でも、適応的手法を用いれば式１１を反転できる。第一の工程として、光音響画像がスパースであると思われるライブラリーを見つけ出す必要がある。このようなライブラリーが見つかれば、グリッド上の光音響画像は、ｈ＝Ｑαとして表せる。式中、αは、ライブラリー要素に関する係数を含むベクトルであり、Ｑは、ライブラリー要素を定義する行列である。ｈのスパース表現に関して、α内の要素の殆どは、ゼロに近いものでなくてはならない。反転は、式１２に関して少ない誤差を返す最もスパースな信号を見つけ出すことによって実行される。
［式１５］
【００９８】
【数１５】
0005566456.tif 000016
【００９９】
ここで、εは許容限度内の再構成誤差である。なお式１５の最適化問題の解を求めるための方法については、下記の第３項で説明する。
【０１００】
部分的な情報だけを入手できるときに再構成の問題に対処する他の方法は、式１４に正規化項を加えること、又は特異値分解を行うことである。この分解における最後の特異値は、再構成画像におけるアーチファクトに対応し得るのであるが、これらの特異値又はこれに対応する固有の画像に関する何らかの関数によって切り捨てるか減衰させてもよい。
【０１０１】
本発明の行列反転法を試験するため、マサチューセッツ州ネーティックのマスワークス社（Mathworks Inc.）製マトラボ（Matlab）で実装された、２．６７ＧＨｚのインテル（登録商標）コア２クアッド（商品名）プロセッサＣＰＵを用いて数値シミュレーションを行った。一例として、平坦な周波数応答をする検出器を想定して、数ｃｍの大きさの画像を再構成した。ここで反転に用いた投影数は８０、各軸に対するグリッド点の数は８０であった。また音場に用いた時間的解像度は、対応する空間的解像度の２倍であった。行列を構築し反転させるためのランタイムは約９分であった。逆行列が計算されると、約０．２秒で再構成画像が得られた。行列の大きさは２０２４０×６４００、メモリ消費量は約１ＧＢであった。生成された音響信号の正確な解析表現が存在するＯＡＴ画像を再構成して、画像の正確な再構成を得た。この結果を、従来のフィルタ処理済み逆投影アルゴリズムを用いたときの結果と比べたところ、強いアーチファクトが逆投影に現れたが、本発明によるモデルをベースにした再構成では発生しなかった。実験テストにおいて、モデル行列を構築し反転させるためのランタイムは、所与の実験設定に関して約１時間であった。また逆行列を計算した後、各実験に関する画像の回復には約１．３秒を要したのみであった。ここで逆行列のサイズは８２８１×２４１２０であり、メモリ消費量は約３ＧＢであった。
【０１０２】
さらに、上記の方法に準拠して再構成された光音響画像ｈは、得るべきエネルギーデポジション画像としても取り扱える。好ましくは、エネルギーデポジション画像は、以下の分解を受ける。
３．スパース表現分解法
【０１０３】
光音響画像を少なくとも２つの成分に分解するために、スパース表現法が用いられる。スパース表現の基本概念は、自然信号を適切な関数基準又はライブラリーにおいてスパースに表現できるということである。すなわち、これらの自然信号は、少数の初等関数の和として表せる。ノイズ除去、圧縮、及び圧縮付き検出においては、スパース表現法が有用であることが判明している（例えばＳ．Ｓ．チェン他、「ＳＩＡＭレビュー」、４３巻、２００１年、１２９〜１５９頁（非特許文献７）を参照）。
【０１０４】
以下では分解法について、被写体における吸収係数と更なる画像成分の局部的光フルエンスに依拠した画像を例示的に参照しながら説明する。例えば、光学的散乱係数の分布又は特定バイオマーカーの分布のような、他の又は更なる画像成分を分離するために分解することは、アナログ的な方法で実行できる。
【０１０５】
光音響画像ｈ（ｘ，ｙ）は、吸収係数μa（ｘ，ｙ）と局部的光フルエンスＵ（ｘ，ｙ）との積に比例しており、すなわち、ｈ（ｘ，ｙ）∝μa（ｘ，ｙ）Ｕ（ｘ，ｙ）である。スパース表現法を用いるために、まず信号を積の公式から和の公式に変換する必要がある。これは、光音響画像の対数をとることで達成される。
［式１６］
【０１０６】
【数１６】
0005566456.tif 000017
【０１０７】
関数ｌｏｇ［μa（ｘ，ｙ）］及びｌｏｇ［Ｕ（ｘ，ｙ）］は、これら両方をｌｏｇ［ｈ（ｘ，ｙ）］から抽出できるようにする非常に際立った特徴を有する。関数ｌｏｇ［Ｕ（ｘ，ｙ）］は全体え的であり、ゆっくりと変化するが、関数ｌｏｇ［μa（ｘ，ｙ）］は局所的であり、速い変化を有し得る。関数ｌｏｇ［μa（ｘ，ｙ）］の局所的な変化に比して、関数ｌｏｇ［Ｕ（ｘ，ｙ）］の全体的な変化は小さい。
【０１０８】
モデル被写体に関する（ａ）光音響画像、（ｂ）吸収係数、及び（ｃ）フルエンスの３Ｄメッシュを例示的に図６に示す（印刷上の事由によりグレー値で表示）。図６は、吸収係数及びフルエンス両方の空間的特性が光音響画像において識別可能であることを示している。特にｌｏｇ［μa（ｘ，ｙ）］及びｌｏｇ［Ｕ（ｘ，ｙ）］両方の特長は、関数ｌｏｇ［ｈ（ｘ，ｙ）］において識別可能である。
【０１０９】
本発明によるｌｏｇ［Ｕ（ｘ，ｙ）］の分解は、以下の補完的条件を満たせる２つのライブラリー｛φn｝n及び｛ψn｝nを見つけ出せるという仮定に基づいている。
１）関数ｌｏｇ［μa（ｘ，ｙ）］は、｛φn｝nにおいてスパースに表現できるが、｛ψn｝nにおいてはスパースに表現できない。
２）関数ｌｏｇ［Ｕ（ｘ，ｙ）］は、｛ψn｝nにおいてスパースに表現できるが、｛φn｝nにおいてはスパースに表現できない。
【０１１０】
以下において、これらの条件をスパース条件（又は互換性条件）と称することにする。スパース条件が具備されると、共同ライブラリー｛φn，ψm｝n,mにおいて光音響画像をスパースに表現できる。
［式１７］
【０１１１】
【数１７】
0005566456.tif 000018
【０１１２】
係数ｃn及びｄmとして、式１８を用いてｌｏｇ［μa（ｘ，ｙ）］及びｌｏｇ［Ｕ（ｘ，ｙ）］を計算できる。
［式１８］
【０１１３】
【数１８】
0005566456.tif 000019
【０１１４】
光音響画像を分解するために、スパース条件を満たす２つのライブラリー（例えば関数基準）が選択される。
【０１１５】
ｌｏｇ［Ｕ（ｘ，ｙ）］に関しては、２次元離散フーリエ基準が好ましい。フーリエ基準は、ｌｏｇ［Ｕ（ｘ，ｙ）］のように平滑で全体的な関数のスパース表現（sparse representation）には非常に有用であるものの、ｌｏｇ［μa（ｘ，ｙ）］のように局所的な細部を含む関数をスパース表現するためには利用できない。
【０１１６】
ｌｏｇ［μa（ｘ，ｙ）］のような一般画像をスパースに表現できるライブラリーを見つけるという命題は、過去１０年に亘って信号処理技術分野で広く研究されてきた。一つの可能性としては、ｌｏｇ［μa（ｘ，ｙ）］を表現するために２次元の離散ハールウェーブレット（Haar wavelet）基準の利用が挙げられる。ハールウェーブレットは、局所的であり、区分的に一定な関数であって、自然画像をスパースに表現するのに利用できる最も簡単なライブラリーを構成する。これらのウェーブレットは、不連続性を有する画像をスパースに表現することができ、このため特に光音響画像がファントムの場合にはｌｏｇ［μa（ｘ，ｙ）］を表現するのに適していることが判明している。またハールウェーブレットは、生物学的組織をインビボで撮像する際に得られる光音響画像に利用することもできる。しかしながら、このような画像にはもっと複雑なディクショナリが好ましい。特に自然画像をスパースに表現するための最適なライブラリーは、撮像条件及び／又は解析される画像の種類に依拠して選択できる。ライブラリーは、種々の明確な関数及び関数基準を備えていてもよい。並進変異型又は不変異型のウェーブレット（例えばハールウェーブレット（Haar wavelets）、ドベシイウェーブレット（Daubechies wavelets）、コイフレット（Coiflets）等）、カーブレット（curvelets）、コントゥアレット（contourlets）、エッジレット（edgelets）、学習過程を通して習得した関数、ユーザ定義の関数等が挙げられる。
【０１１７】
さらに、この問題の性質に関するアプリオリな情報がライブラリーの選択を改善することもできる。例えば、画像が血管を含んでいると思われるときには、吸収係数のライブラリーにカーブレットを追加してもよい。これは、カーブレットが１次元の曲線をスパースに表現できるからである。さらに、画像の特性に従って手動的又は自動的に関数を定義できる。良好に定義された被写体が画像に出現したら、（手動的又は自動的に）境界を検出でき、境界内の定数及び境界外のゼロを許容する関数を指定できる。またこれらの関数は、吸収係数のライブラリーに加えることもできる。光フルエンスを表すライブラリーに要素を加えることもできる。ゼロ次変形ベッセル関数の逆関数は、散乱性の媒体における光の２Ｄ伝播を説明するものであって、
【０１１８】
【数１９】
0005566456.tif 000020
を加えることができる。異なる減衰率及び最小限のフルエンス位置それぞれに対応して、ｋ又は
【０１１９】
【数２０】
0005566456.tif 000021
のいくつかの値を選択してもよい。
【０１２０】
分解は、局部的特性及び全体的特性間で識別することに基づいているので、吸収係数のバックグラウンド値は抽出できない。したがって、ｌｏｇ［μa（ｘ，ｙ）］の再構成は付加定数まで正確であり、μa（ｘ，ｙ）の再構成はスケーリングまで正確である。実験で得られた吸収係数を拡大縮小するために、実験的な設定を較正する。これを行う一の方法は、被写体の表面への照射を推定することによるものである。これを用いて、表面上のμa（ｘ，ｙ）を見つけ出し、倍率を回復できる。他の方法としては、システムを較正するために、既知の吸収を有するファントムを用いることが考えられよう。
【０１２１】
フーリエ基準及びハール基準の組み合わせがスパース条件を満たすことを示していることを例示的に図７に示す。この例では、関数が図６の関数に合致するようにスケーリングされている。図７ａ及び図７ｂは、２００個及び１００個の要素だけでそれぞれ計算されたｌｏｇ［μa（ｘ，ｙ）］及びｌｏｇ［Ｕ（ｘ，ｙ）］（グリッド上で離散化され表現されている）を示している。この図は、両方の関数がそれぞれの基準においてスパースに表現できることをはっきりと示している。図７ｃ及び図７ｄにおいて、２つの関数は、対向した基準で表現され、要素数を２５０個に増加させている。この図は、増加した要素数にも関わらず、関数が両方の基準において良好には表現されていなかったことをはっきりと示している。
【０１２２】
スパース条件を満たす２つのライブラリーを選択すると、分解の問題は、式１７を満足させる係数ｃn及びｄmの最小数を探すことになる。この離散型バージョンにおいて、分解の問題は、行列形式で書き表せる。
［式１９］
【０１２３】
【数２１】
0005566456.tif 000022
【０１２４】
上式中、αは、要素がｃn及びｄmである列ベクトルであり、α＝［ｃ1，．．．，ｃn，ｄl，．．．，ｄM］である。ｘは、要素としてｌｏｇ（ｈk）を有する列ベクトルであり、Ａは、列がライブラリー関数φn,k及びψn,kに対応するＫ×（Ｎ＋Ｍ）行列である。ゼロノルム（zero norm）｜｜・｜｜0は、所定のベクトルにおける非ゼロ要素の数である。実際のケースにおいて、式１７及び式１８は単なる近似値であるので、式１９を式２０と置換する必要がある。
［式２０］
【０１２５】
【数２２】
0005566456.tif 000023
【０１２６】
ここで、εは許容限度内の再構成誤差である。
【０１２７】
式２０の最適化は非凸問題であり、一般に非多項の複雑性（non-polynomial complexity）を用いて解を求め得るに過ぎない。式２０の解を近似化するための２つの注目すべき方法は、基準追求型アルゴリズム及び照合追求型アルゴリズムである。基準追求型は式２０におけるゼロノルムをｌ1ノルムと置換することに基づいている。ｌ1ノルムにおいては凸問題であり、線形プログラミングを使ってゼロノルム問題よりも効率的に解を求められる。表現が十分なまでにスパースであれば、両方のノルムで得られた結果は、同一のものであることが示されてきた。照合追求型（ＭＰ）は、式１９及び式２０の解を求めるための反復アルゴリズムである。ＭＰにおいて、元のライブラリーの要素から画像を反復的に構築する。各工程において、ライブラリーからの一の要素が画像の表現に追加される。信号が十分な正確さで表現されるか、反復回数が超過してしまうと、反復は停止する。直交する照合追求型（ＯＭＰ）と呼ばれている、ＭＰの変形型も利用され得る。
【０１２８】
２つのＭＰアルゴリズム間の主な違いは、選択されたライブラリー要素に関する係数を計算するときの方法にある。標準的なＭＰにおいては、各ライブラリー要素の係数は、一回だけ、すなわち、要素が信号表現に追加されるときだけ選択される。したがって、反復工程のいずれにおいても誤差は、他の反復において修正できない。ＯＭＰにおいては、すべての要素に関する最適な係数は、反復ごとに計算し直されるので、前方の反復における誤差に関する補償が可能となる。
【０１２９】
式２０の解を求めた後、解ベクトルαI（Ｉはベクトルの大きさを示す）が得られ、そこでは、殆どの要素がゼロに等しい（又はゼロに非常に近い）。そして、αIにおいて非ゼロの要素に対応するＡの行のみを含んでいる換算行列ＡIが得られる。
【０１３０】
行列ＡIから、２つの新たな行列Ｂ及びＣが構築される。行列Ｂ及びＣは、吸収係数基準｛φn｝n及びフルエンス基準｛ψn｝nそれぞれに対応したＡIの列で構成される。したがって、２つのベクトルβ及びγは、αIから構築される。ｂ及びｃそれぞれによって表示されるｌｏｇ［μa（ｘ，ｙ）］及びｌｏｇ［Ｕ（ｘ，ｙ）］のベクトル表現は、式２１を用いて計算される。
［式２１］
【０１３１】
【数２３】
0005566456.tif 000024
【０１３２】
上記の方法に従って再構成された局部的吸収係数の画像μは、従来の画像処理技法又はデータ処理技法を用いてさらに処理できる。
【０１３３】
図７におけるような実施例で、またより複雑な形状をなす画像で行った数値シミュレーションは、モデルのデータがノイズのデータと重ね合わされても、被写体の吸収係数画像を回復することができることを証明した。さらにまた、この方法は、典型的な小動物の全身撮像形態及び対応する組織特性及び他の実験パラメータに似せた、組織を模したファントムを用いた実験において首尾良くテストされた。
４．正規化及びＭＳＯＴへの適用
【０１３４】
異なる励起波長に対する光音響画像を考えると、分解は、画像自体に適用され得るだけでなく、異なる画像同士及びその組み合わせ同士の間の正規化された差及び／又は比率にも適用され得る。例えば、２つの光音響画像ｈ1及びｈ2が２つの異なる照射波長において得られたとすれば、これらの差をｄｈで、これらの和をｈで表示できる。ｄｈ＜＜ｈであると仮定すると、以下の関係が得られる。
［式２２］
【０１３５】
【数２４】
0005566456.tif 000025
【０１３６】
式２２は、光音響画像の差が平均フルエンスと共にフルエンスの変化にも依存することを示している。画像の差を平均的な画像で正規化すれば、式２３となる。
［式２３］
【０１３７】
【数２５】
0005566456.tif 000026
【０１３８】
正規化された示差の画像は、２つの別個な画像の和、すなわち吸収係数同士の和及び光フルエンス同士の和である。この結果、単一の光音響画像の分解を可能にした同じ引数が本明細書でも適用される。ｄｈ／ｈに対して分解アルゴリズムを用いることにより、ｄμ／μ及びｄＵ／Ｕの両方を抽出できる。μ及びＵは、元の光音響的画像を分解することから得られるので、これらを用いて最終的にｄμ及びｄＵを得ることができる。式２３は、別の対数形式の式２４で書き表せる。
［式２４］
【０１３９】
【数２６】
0005566456.tif 000027
【０１４０】
両方の形式に関して同じ解析を利用できる。
【０１４１】
吸収係数の変化を正確に抽出できることは、スペクトルフィッティングを行うことを可能にし、これを使って、スペクトル的に変化する吸収係数を有する多様な物質の濃度を計算できる。連続的な波長λ1、λ2、・・・、λ3に関して得られた光音響画像ｈ1、ｈ2、・・・、ｈnの集合があるとして、異なる吸収係数ｄμ1＝μ1−μ2，．．．，ｄμn-1＝μn-1−μnは、上述のプロセスを用いて計算できる。集合｛ｄμ1，ｄμ2，．．．，ｄμn-1｝は、いかなる他の測定値の波長においても吸収係数と共に、すべての測定波長｛μ1，μ2，．．．，μn｝に関する（対する）吸収係数の集合を一意的に決定する。２つの集合のそれぞれ、又はそれらの処理版のうち任意のものは、これらの吸収スペクトルに従って異なる物質の濃度を決定するスペクトルフィッティングのアルゴリズムのための入力値として使用できる。スペクトルフィッティングは、以下の２乗誤差を最小化することにより得られる。
［式２５］
【０１４２】
【数２７】
0005566456.tif 000028
【０１４３】
式中、ηk'は物質ｋ’であり、Ｓk'（λn'）は、波長λn'に対するこれの吸収スペクトルである。スペクトルフィッティングは、主成分分析、独立成分分析、又は他種の発生源分別のような、従来技術において知られる他の方法で実行できる。画像のスパースをフィッティングにおいて利用してもよい。
５．システムの好ましい特長
【０１４４】
本発明の光音響反転法を行う画像再構成装置２２に関する更なる形態の詳細を図８で概略的に示している（第２項を参照）。この画像再構成装置２２は、圧力信号を信号収集装置１０から受信する第一入力装置２２１（図１を参照）と、伝播行列を調整するための測定条件を受信する少なくとも一の第二入力装置２２２、２２３と、行列調整回路２２４と、光音響反転回路２２５と、画像データ出力部２２６とを備える。光音響反転回路２２５は、例えば式１４の最小化に基づいて、式１１に関する上記反転を実行するよう構成されている。行列調整回路２２５を画像再構成装置２２に統合することは、厳密には必要とされない。伝播行列は別個に事前計算でき、さらに調整することなく光音響反転回路２２４へ入力できる。
【０１４５】
上述したように、エネルギーデポジション画像の再構成は、従来の解析法又は数値法、例えば反復伝播モデルベースのシミュレーションを含んでもよい。これに応じて適用した画像再構成装置２２の実施形態を図９に概略的に示す。この画像再構成装置２２は、圧力信号を信号収集装置１０から受信する入力装置２２１（図１を参照）と、開始モデル調整回路２２７と、伝播モデルシミュレーション回路２２８と、画像データ出力部２２６とを備える。この伝播モデルシミュレーション回路２２８は、圧力信号に関する反復伝播モデルベースのシミュレーションを実行するように構成されている。
【０１４６】
図１０によれば、本発明の分解法（第３項を参照）を実行する吸収画像分解装置２３は、エネルギーデポジション画像データ入力部２３１と、対数回路２３２と、ライブラリー係数最適化回路２３３と、第一出力部２３６及び第二出力部２３７とを備えている。これらの出力部２３６、２３７は、吸収係数画像及び例えばフルエンス画像それぞれを含む分解された画像成分を伝達する。ライブラリー係数最適化回路２３３は、局所ライブラリー選択回路２３４及び全体ライブラリー選択回路２３５に接続される。
【０１４７】
上述の詳細説明、図面、及び特許請求の範囲における本発明の特長は、単独のままであっても組み合わせた状態であっても、多様な実施形態において本発明を実現するために重要なものであり得る。
【符号の説明】
【０１４８】
１…被写体
１０…信号収集装置
１１…レーザー源
１２…音響信号記録装置、検出装置、検出器アレイ
１３…メンブレン
１４…水槽
１５…音響検出要素
２０…画像処理装置
２１…信号インターフェイス部
２２…再構成装置、画像再構成装置
２３…分解装置、吸収画像分解装置
２４…インターフェイス
３０…他の制御装置
１００…撮像装置
２２１…第一入力装置
２２２…第二入力装置
２２３…第二入力装置
２２４…行列調整回路
２２５…光音響反転回路
２２６…画像データ出力部
２２７…開始モデル調整回路
２２８…伝播モデルシミュレーション回路
２３１…エネルギーデポジション画像データ入力部
２３２…対数回路
２３３…ライブラリー係数最適化回路
２３４…局所ライブラリー選択回路
２３５…全体ライブラリー選択回路
２３６…第一出力部
２３７…第二出力部
【図１】
0005566456.tif 000029
【図２】
0005566456.tif 000030
【図３】
0005566456.tif 000031
【図４】
0005566456.tif 000032
【図５】
0005566456.tif 000033
【図８】
0005566456.tif 000034
【図９】
0005566456.tif 000035
【図１０】
0005566456.tif 000036
【図６】
0005566456.tif 000037
【図７】
0005566456.tif 000038