(57)【特許請求の範囲】
【請求項１】
物体(1)のマルチスペクトル光音響イメージングのためのデバイスであって、
a)2つ以上の異なる照射波長(λ)の電磁放射線を前記物体(1)に照射するように構成された照射ユニット(10)であって、前記電磁放射線が、時間変化する強度を有する、照射ユニット(10)と、
b)前記異なる照射波長(λ)の前記電磁放射線を前記物体(1)に照射する際に前記物体(1)内に生成される音響波を検出するように構成された検出ユニット(20)と、
c)前記異なる照射波長(λ)の各々で前記物体(1)内に生成され、検出された前記音響波に基づいて、前記物体(1)の画像(P(r,λ))を再構成し、
前記物体(1)内の少なくとも1つの電磁放射線吸収体(HbO2、HHb)の濃度に関連する少なくとも1つの第1の濃度値(c'HbO2、c'HHb)の空間分布を決定するように構成された処理ユニット(30)と
を備え、
前記少なくとも1つの第1の濃度値(c'HbO2、c'HHb)の前記空間分布の決定は、
前記異なる照射波長(λ)での再構成された前記画像(P(r,λ))と、
前記物体(1)内の前記少なくとも1つの電磁放射線吸収体(HbO2、HHb)の少なくとも1つの波長依存減衰係数(ε1(λj)、ε2(λj))と、
少なくとも2つのモデルスペクトル(ΦM(λ)、Φi(λ))の線形結合と
に基づき、
前記モデルスペクトル(ΦM(λ)、Φi(λ))は、前記物体(1)内の放射線フルエンス(Φ(r,λ))の波長依存基底関数、または放射線フルエンス(Φ(r,λ))の関数(Φ'(r,λ))の波長依存基底関数を表し、
前記放射線フルエンス(Φ(r,λ))の関数(Φ'(r,λ))は、正規化された放射線フルエンス(Φ(r,λ)/||Φ(r)||)を含む、デバイス。
【請求項２】
前記処理ユニット(30)が、
前記物体(1)内の少なくとも2つの電磁放射線吸収体(HbO2、HHb)の濃度に関連する少なくとも2つの第1の濃度値(c'HbO2、c'HHb)の空間分布を決定し、
前記少なくとも2つの第1の濃度値(c'HbO2、c'HHb)の決定された前記空間分布から第2の濃度値(sO2(r))の空間分布を得るように構成され、
前記第2の濃度値(sO2(r))は、相対濃度値を含み、
前記少なくとも2つの第1の濃度値(c'HbO2、c'HHb)の前記空間分布の決定は、
前記異なる照射波長(λ)での再構成された前記画像(P(r,λ))と、
前記物体(1)内の前記少なくとも2つの電磁放射線吸収体(HbO2、HHb)の前記波長依存減衰係数(ε1(λj)、ε2(λj))と、
前記少なくとも2つのモデルスペクトル(ΦM(λ)、Φi(λ))の線形結合と
に基づく、請求項1に記載のデバイス。
【請求項３】
前記少なくとも2つのモデルスペクトル(ΦM(λ)、Φi(λ))が、前記物体(1)内にあると仮定された前記少なくとも1つの電磁放射線吸収体(HbO2、HHb)のために決定されている最小化した数のモデルスペクトル(ΦM(λ)、Φi(λ))に対応する、請求項1または2に記載のデバイス。
【請求項４】
前記少なくとも2つのモデルスペクトル(ΦM(λ)、Φi(λ))が、
媒体内の異なる深さと、前記媒体内の前記少なくとも1つの電磁放射線吸収体(HbO2、HHb)の異なる濃度とのためのスペクトルパターンのセット(ΦZ,OX(λ))を決定することと、
前記スペクトルパターンのセット(ΦZ,OX(λ))に統計的手順を適用することと
によって決定されている、請求項1から3のいずれか一項に記載のデバイス。
【請求項５】
前記スペクトルパターンのセット(ΦZ,OX(λ))が、
異なる波長(λ)のための拡散方程式または放射伝達方程式の解析解(ΦZ,OX(λ)=Φ0exp(-μeffZ))または数値解、および/または、
異なる波長(λ)のための光伝播のモンテカルロシミュレーション、および/または、
異なる波長(λ)のための実験測定値のセット
から得られている、請求項4に記載のデバイス。
【請求項６】
前記スペクトルパターンのセット(ΦZ,OX(λ))に適用される前記統計的手順が、
主成分分析(PCA)、および/または、
カーネル主成分分析(カーネルPCA)、および/または、
線形判別分析(LDA)
を含む、請求項4または5に記載のデバイス。
【請求項７】
前記モデルスペクトル(ΦM(λ)、Φi(λ))の線形結合が、平均モデルスペクトル(ΦM(λ))と、第1の数(p)のさらなるモデルスペクトル(Φi(λ))の線形結合との加算に対応し、前記さらなるモデルスペクトル(Φi(λ))の各々が、モデルパラメータ(mi)を乗算され、ΦM(λ)+Σpi=1miΦi(λ)となる、請求項1から6のいずれか一項に記載のデバイス。
【請求項８】
前記照射ユニット(10)が、第2の数(w)の異なる照射波長(λ)の電磁放射線を前記物体(1)に照射するようにさらに構成されており、
前記第2の数(w)は、前記さらなるモデルスペクトル(Φi(λ))の前記第1の数(p)と、決定すべき前記第1の濃度値(c'HbO2、c'HHb)の数(nc)との和以上であり、w≧p+ncである、請求項7に記載のデバイス。
【請求項９】
前記物体(1)内の2つの異なる電磁放射線吸収体の濃度に関連する2つの第1の濃度値(c'1(r)、c'2(r))の空間分布の決定が、連立非線形方程式の逆算を含み、
前記連立非線形方程式の逆算は、画像領域内の分散された異なる位置(r)の選択されたセットに対して同時に実行される、請求項1から8のいずれか一項に記載のデバイス。
【請求項１０】
前記物体(1)内の2つの異なる電磁放射線吸収体の濃度に関連する2つの第1の濃度値(c'1(r)、c'2(r))の空間分布の決定が、
P(r,λj)=(ΦM(λj)+Σpi=1mirΦi(λj))(c'1(r)ε1(λj)+c'2(r)ε2(λj))、j=1...w
によって与えられる連立非線形方程式の逆算を含み、
ここで、P(r,λj)が、位置rおよび照射波長λjでの前記物体(1)の再構成された前記画像の強度を表し、ΦM(λj)+Σpi=1mirΦi(λj)が、照射波長λjでのi=1...pとした前記モデルスペクトル(ΦM(λ)、Φi(λ))の線形結合を表し、mirが、位置rでのモデルパラメータの値を表し、c'1(r)ε1(λj)+c'2(r)ε2(λj)が、照射波長λjでの前記2つの電磁放射線吸収体の電磁放射線減衰係数ε1(λj)およびε2(λj)の線形結合を表し、c'1(r)およびc'2(r)が、位置rでの前記2つの第1の濃度値を表す、請求項1から9のいずれか一項に記載のデバイス。
【請求項１１】
物体(1)のマルチスペクトル光音響イメージングのための方法であって、
a)2つ以上の異なる照射波長(λ)の電磁放射線を前記物体(1)に照射するステップであって、前記電磁放射線が、時間変化する強度を有する、ステップと、
b)前記異なる照射波長(λ)の前記電磁放射線を前記物体(1)に照射する際に前記物体(1)内に生成される音響波を検出するステップと、
c)前記異なる照射波長(λ)の各々で前記物体(1)内に生成され、検出された前記音響波に基づいて、前記物体(1)の画像(P(r,λ))を再構成し、前記物体(1)内の少なくとも1つの電磁放射線吸収体(HbO2、HHb)の濃度に関連する少なくとも1つの第1の濃度値(c'HbO2、c'HHb)の空間分布を決定するステップと
を有し、
前記少なくとも1つの第1の濃度値(c'HbO2、c'HHb)の前記空間分布の決定は、
前記異なる照射波長(λ)での再構成された前記画像(P(r,λ))と、
前記物体(1)内の前記少なくとも1つの電磁放射線吸収体(HbO2、HHb)の少なくとも1つの波長依存減衰係数(ε1(λj)、ε2(λj))と、
少なくとも2つのモデルスペクトル(ΦM(λ)、Φi(λ))の線形結合と
に基づき、
前記モデルスペクトル(ΦM(λ)、Φi(λ))は、前記物体(1)内の放射線フルエンス(Φ(r,λ))の波長依存基底関数、または放射線フルエンス(Φ(r,λ))の関数(Φ'(r,λ))の波長依存基底関数を表し、
前記放射線フルエンス(Φ(r,λ))の関数(Φ'(r,λ))は、正規化された放射線フルエンス(Φ(r,λ)/||Φ(r)||)を含む、方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、マルチスペクトル光音響イメージングのためのデバイスおよび方法に関し、具体的には、独立請求項に記載の、物体内の1つまたは複数の発色団の濃度に関する定量的な画像を決定するためのデバイスおよび方法に関する。
【背景技術】
【０００２】
光音響イメージングは、光音響効果に基づいており、光音響効果に従って、超音波が、物体、たとえば生体組織による電磁放射の吸収と、その後の物体の熱弾性膨張とに起因して発生される。励起放射線、たとえば、レーザ光または高周波は、短いパルス持続時間を有するパルス放射線、または、変化する振幅または周波数を有する連続放射線であり得る。
【０００３】
光音響における重要な進歩は、それらの個別の吸収スペクトルを分離することによる異なる部分のイメージングである。これは、典型的には、撮像される物体を複数の波長で照明し、収集されたデータにおいてスペクトルアンミキシング操作を実行することによって達成される。そのような部分のイメージングは、様々な用途にとって重要である。
【０００４】
たとえば生物学的および医学的用途において、組織血液の酸素化は、組織生存性、代謝、低酸素症、および神経細胞活性化に関する価値のある知識を推論することができる重要な指標である。高い空間分解能および時間分解能で深部組織血液酸素化の定量的な情報を提供することができるイメージングモダリティは少ない。いくつかの励起波長を利用することによって、高い空間分解能および時間分解能のいわゆるマルチスペクトル光音響画像(分光光音響画像とも称される)が得られ得、これは、血液酸素化を推論するために使用され得る組織の光吸収に関連する。しかしながら、マルチスペクトル光音響トモグラフィ(multispectral optoacoustic tomography，MSOT)による深部組織中の酸素化ヘモグロビンおよび脱酸素化ヘモグロビンなどの発色団の空間分布の正確な定量化は、特に3mmを超える深さで課題を残している。
【０００５】
光音響画像は、吸収エネルギー密度H(r,λ)(ここで、rは、位置を示し、λは、励起源の波長を示す)に比例し、吸収エネルギー密度は、組織の光吸収係数μα(r,λ)および光フルエンスΦ(r,λ)の空間的な分布に依存し、すなわち、H(r,λ)=μα(r,λ)Φ(r,λ)である。μα(r,λ)が既知であると仮定すると、組織吸収体の濃度は、線形スペクトルアンミキシング法を使用して、各位置rにおいて推定され得る。したがって、価値のある生理学的情報(たとえば、組織内の血液酸素化)および分子情報(たとえば、分子プローブの空間的定位および濃度推定)が抽出され得る。しかしながら、得られた情報は、未知の光フルエンスΦ(r,λ)に結合されるので、μα(r,λ)の分布は、生の光音響画像から明示的に知られない。吸収エネルギー密度H(r,λ)を光フルエンスΦ(r,λ)および光吸収係数μα(r,λ)に分解することは、MSOT問題における大きな課題である。
【０００６】
定量的光音響イメージング(定量的光音響トモグラフィまたはqPATとも称される)の問題を解決するための可能な定式化は、モデルベースの光学特性推定手法である。この定式化の下では、光フルエンスΦ(r,λ;μα(r,λ),μs(r,λ))は、光伝播モデルを使用して、画像化された組織の光吸収係数μα(r,λ)および散乱係数μs(r,λ)の空間分布に関連付けられる。単一または複数の波長における吸収エネルギー密度H(r,λ)を仮定して、汎関数||H(r,λ)-Φ(r,λ;μα(r,λ),μs(r,λ))μα(r,λ)||を最小化することによって、光フルエンスおよび組織の光学特性(吸収係数および散乱係数)を同時に抽出しようとするいくつかの非線形モデルベースの反転方法が公開されている。光学特性推定方法は、信号または画像領域において動作することができ、それらは、(光散乱が、一意性を達成するために既知である必要がある場合)1つの励起波長または複数の励起波長のもとで動作することができ、かつそれらは、音響伝播のモデルと組み合わされてもよい。それらの共通の特徴は、それらが、組織の吸収係数の定量値を得るために汎関数||H(r,λ)-Φ(r,λ;μα(r,λ),μs(r,λ))μα(r,λ)||を最小化しようとすることである。
【０００７】
光学特性推定手法の主な欠点は、それが典型的には、照射された組織容積全体における吸収エネルギー密度H(r,λ)の正確な推定に依存することである。典型的には、光音響画像P(r,λ)は、吸収エネルギー密度H(r,λ)と、空間依存スケーリング係数C(r)とに比例し、すなわち、P(r,λ)=C(r)H(r,λ)であり、ここで、C(r)は、システムの空間インパルス応答と(一般的には未知の)空間的に変化するグリュナイゼン係数とに依存する。したがって、光音響画像からのH(r,λ)の推定は、(システムの空間インパルス応答または超音波減衰などの)典型的には断層再構成に含まれないスケーリング効果、または異なる組織のタイプのグリュナイゼン係数を推定する際の不確実性に起因するエラーを起こしやすい可能性がある。さらに、(限定されたビューシステムの場合、または超音波散乱事象下で共通する)照射された容積の一部で、吸収エネルギー密度が未知であるか大幅に誤っている場合、反転プロセスは、損なわれる可能性がある。最後に、光学特性推定手法は、数千の未知のパラメータによる大規模な反転問題を規定し、これは、高い計算の複雑さをもたらす。これらの理由のために、生物学的組織などの複雑な構造の実験的な組織の測定の場合におけるそのような方法論の適用は、これまでのところ限定されている。
【発明の概要】
【発明が解決しようとする課題】
【０００８】
本発明の目的は、物体内の1つまたは複数の発色団の濃度に関連する少なくとも1つの値の空間分布の改善された定量化を可能にするマルチスペクトル光音響イメージングのためのデバイスおよび方法を提供することである。
【課題を解決するための手段】
【０００９】
この目的は、独立請求項に記載のデバイスおよび方法によって達成される。
【００１０】
本発明による物体のマルチスペクトル光音響イメージングのためのデバイスは、2つ以上の異なる照射波長の電磁放射線を物体に照射するように構成された照射ユニットであって、電磁放射線が時間変化する強度を有する、照射ユニットと、異なる照射波長の電磁放射線を物体に照射する際に物体内に生成される音響波を検出するように構成された検出ユニットと、照射波長の各々において物体内に生成され、検出された音響波に基づいて、物体の画像を再構成し、物体内の少なくとも1つの電磁放射線吸収体の濃度に関連する少なくとも1つの第1の濃度値の空間分布を決定するように構成された処理ユニットとを備え、少なくとも1つの第1の濃度値の空間分布の決定は、異なる照射波長において再構成された画像と、物体内の少なくとも1つの電磁放射線吸収体の少なくとも1つの波長依存減衰係数と、少なくとも2つのモデルスペクトルの線形結合とに基づき、モデルスペクトルは、物体内の放射線フルエンスの波長依存基底関数、または放射線フルエンスの関数の波長依存基底関数を表す。好ましくは、放射線フルエンスは、正規化された放射線フルエンスであり、および/または、モデルスペクトルは、放射線フルエンスの正規化された波長依存性を表す。
【００１１】
本発明による物体のマルチスペクトル光音響イメージングのための方法は、以下のステップ、すなわち、2つ以上の異なる照射波長の電磁放射線を物体に照射するステップであって、電磁放射線が時間変化する強度を有する、ステップと、異なる照射波長の電磁放射線を物体に照射する際に物体内に生成される音響波を検出するステップと、照射波長の各々において物体内に生成され、検出された音響波に基づいて、物体の画像を再構成し、物体内の少なくとも1つの電磁放射線吸収体の濃度に関連する少なくとも1つの第1の濃度値の空間分布を決定するステップとを含み、少なくとも1つの第1の濃度値の空間分布の決定は、異なる照射波長での再構成された画像と、物体内の少なくとも1つの電磁放射線吸収体の少なくとも1つの波長依存減衰係数と、少なくとも2つのモデルスペクトルの線形結合とに基づき、モデルスペクトルは、物体内の放射線フルエンスの波長依存基底関数、または放射線フルエンスの関数の波長依存基底関数を表す。好ましくは、放射線フルエンスは、正規化された放射線フルエンスであり、および/または、モデルスペクトルは、放射線フルエンスの正規化された波長依存性を表す。
【００１２】
本発明は、i)2つ以上の異なる照射波長の過渡電磁放射線を物体に照射する際に物体から放射される音響波から再構成された2つ以上の光音響画像と、ii)物体内の好ましくは正規化された波長依存放射線フルエンスの基底関数に対応する少なくとも2つの波長依存モデルスペクトルの線形結合とを考慮して、少なくとも1つの放射線吸収体の定量的な画像、具体的には相対的な濃度の空間分布、を決定する手法に基づく。本発明の意味の範囲内では、基底関数は、好ましくは、特定の関数空間のための特定の基底の要素に対応し、(ベクトル空間内のすべてのベクトルが基底ベクトルの線形結合として表され得るという事実と同様に)この要素によって、基底関数の線形結合は、特定の関数空間内のすべての連続関数を表すことができる。したがって、本発明の文脈では、「モデルスペクトル」および「基底関数」という用語は、フルエンス固有スペクトルまたは固有スペクトルとも称される。
【００１３】
本発明の手法によれば、物体内の推定された吸収エネルギー密度の画像から放射線フルエンスならびに吸収係数を抽出するために、物体内の光伝播のモデルが使用されるモデルベースの光学特性推定手法を使用する代わりに、波長依存モデルスペクトルのセットの線形結合が使用され、線形結合は、モデルスペクトルのセットの各モデルスペクトルに波長非依存パラメータを乗算し、結果を加算することによって得られる。好ましくは、モデルスペクトルのセットのモデルスペクトルは、あらかじめ、たとえば、少なくとも1つの放射線吸収体を含む代表媒体のシミュレーションおよび/または代表媒体での測定によって決定され、少なくとも1つの放射線吸収体に固有である。
【００１４】
本発明により、具体的には光学特性推定手法を超える以下の利点が達成される。(i)はるかにより小さい次元数を有する小規模な最適化問題が定義され得る。(ii)ピクセル単位での適用が可能になり、これは、照射された容積全体の完全かつ正確な画像再構成が絶対的な必須条件ではないことを意味する。(iii)適用は、画像全体の選択された部分において実行され得、それによって、正確に再構成された領域における収束は、画像の異なる部分内の破損した可能性のある情報に影響されない。(iv)本発明による手法は、吸収エネルギー密度H(r,λ)の絶対値ではなく、吸収エネルギー密度の相対値(すなわち、それによって||H(r)||が位置rにおけるすべての波長にわたる吸収エネルギー密度のノルムに対応する吸収エネルギー密度の波長依存性H(r,λ)/||H(r)||)に依存し、完全なシステムのキャリブレーション、または、組織のグリュナイゼン係数、すなわち、モデルベースの光学特性推定手法のパフォーマンスに影響を与えることが示されているパラメータを推定することなどのさらなる複雑化を回避する。
【００１５】
結果として、本発明は、物体内、具体的にはまた、3mmを超える組織の深さにおける、1つまたは複数の発色団の濃度に関連する少なくとも1つの値の空間分布の正確かつ信頼性の高い定量化を可能にする。
【００１６】
一般に、本発明は、光音響測定を使用して、光を減衰させる媒体の内部の1つまたは複数の発色団の分布を定量化するためのデバイスおよび方法に関する。本発明は、組織酸素化、具体的には、血液酸素飽和度sO2を定量化するのに特に適しているが、有利には、少なくとも1つの光吸収体を含む任意の媒体内部の任意の光吸収体分布の相対濃度を分離および/または定量化する際にも使用され得る。
【００１７】
好ましい実施態様によれば、処理ユニットは、物体内の少なくとも2つの電磁放射線吸収体、たとえば、HbO2およびHHbの濃度に関連する少なくとも2つの第1の濃度値の空間分布を決定するように構成され、少なくとも2つの第1の濃度値の空間分布の決定は、異なる照射波長での再構成された画像と、物体内に存在する少なくとも2つの発色団、たとえば、HbO2およびHHbの少なくとも2つの波長依存減衰係数と、少なくとも2つのモデルスペクトルの線形結合とに基づき、また、処理ユニットは、少なくとも2つの第1の濃度値の決定された空間分布から第2の濃度値の空間分布、具体的には、sO2などの相対濃度値を得るように構成される。
【００１８】
さらに別の好ましい実施態様によれば、少なくとも2つのモデルスペクトルは、物体内にあると仮定された少なくとも1つの電磁放射線吸収体、たとえば、HbO2およびHHbのために決定されている最小化した数のモデルスペクトルに対応する。好ましくは、モデルスペクトルの最小化した数は、3と8との間であり、より詳細には、3と5との間である。
【００１９】
さらに、少なくとも2つの波長依存モデルスペクトルは、媒体内の異なる深さと、媒体内の少なくとも1つの電磁放射線吸収体の異なる濃度とのためのスペクトルパターンのセットを決定し、スペクトルパターンのセットに統計的手順を適用することによって決定されていることが好ましい。
【００２０】
好ましくは、スペクトルパターンのセットは、異なる波長のための拡散方程式または放射伝達方程式の解析解または数値解から、および/または異なる波長のための光伝播のモンテカルロシミュレーションから、および/または異なる波長のための実験測定値のセットから得られている。
【００２１】
好ましくは、スペクトルパターンのセットに適用される統計的手順は、主成分分析(principal component analysis，PCA)および/またはカーネル主成分分析(カーネルPCA)および/または線形判別分析(linear discriminant analysis，LDA)を含む。
【００２２】
別の好ましい実施態様によれば、モデルスペクトルの線形結合は、平均モデルスペクトルΦM(λ)と、第1の数pのさらなるモデルスペクトルΦi(λ)の線形結合との加算に対応し、ここで、さらなるモデルスペクトルΦi(λ)の各々には、モデルパラメータmiが乗算され、ΦM(λ)+Σpi=1miΦi(λ)となる。
【００２３】
好ましくは、照射ユニットは、第2の数wの異なる照射波長の電磁放射線を物体に照射するように構成されており、ここで、第2の数wは、さらなるモデルスペクトルの第1の数pと、決定すべき第1の濃度値の数ncとの和以上であり、w≧p+ncである。
【００２４】
さらに、物体内の2つの異なる電磁放射線吸収体の濃度に関連する2つの第1の濃度値c'1(r)およびc'2(r)の空間分布の決定は、連立非線形方程式の逆算を含むことが好ましく、ここで、連立非線形方程式の逆算は、画像領域における分散された異なる位置rの選択されたセットに対して同時に実行される。
【００２５】
さらに好ましくは、位置rでの物体内の2つの異なる電磁放射線吸収体の濃度に関連する2つの第1の濃度値c'1(r)およびc'2(r)の空間分布の決定は、
P(r,λj)=(ΦM(λj)+Σpi=1mirΦi(λj))(c'1(r)ε1(λj)+c'2(r)ε2(λj))、j=1...w
によって与えられる連立非線形方程式の逆算を含み、ここで、P(r,λj)は、照射波長λjおよび位置rにおける物体の再構成された画像の強度を表し、ΦM(λj)+Σpi=1mirΦi(λj)は、照射波長λjおよび位置rにおけるi=1...pとした波長依存モデルスペクトルΦM(λj)およびΦi(λj)の線形結合を表し、mirは、位置rでのモデルパラメータの値を表し、c'1(r)ε1(λj)+c'2(r)ε2(λj)は、照射波長λjにおける2つの電磁放射線吸収体の電磁放射線減衰係数ε1(λj)およびε2(λj)の線形結合を表し、c'1(r)およびc'2(r)は、位置rでの2つの第1の濃度値を表す。
【００２６】
さらに別の好ましい実施態様によれば、連立非線形方程式の逆算は、n=1...Ngridとして、画像内のグリッドに対応する数Ngridの異なる位置rnに対して同時に実行され、隣接する位置の間の1つまたは複数のモデルパラメータの変動は、制限される。
【００２７】
好ましくは、連立非線形方程式は、画像内のグリッドに対応する複数の位置における再構成された画像強度とモデルとの間の差(すなわち、Σj||P(r,λj)-(ΦM(λj)+Σpi=1mirΦi(λj))(c'1(r)ε1(λj)+c'2(r)ε2(λj))||nrm1pwr1、ここで、nrm1は、適切なノルムを表し、pwr1は、適切な指数を表す)の最小化と、隣接するグリッド点の間のモデルパラメータの変動の同時の最小化とによって解決される。好ましくは、これは、たとえば、画像内のグリッドに対応する複数の位置における再構成された画像強度とモデルとの間の差の線形結合である費用関数fgrid(すなわち、ΣnΣj||P(rn,λj)-(ΦM(λj)+Σpi=1mir(n)Φi(λj))(c'1(rn)ε1(λj)+c'2(rn)ε2(λj))||nrm1pwr1)の最小化と、画像内の隣接する位置に対応するモデルパラメータの変動(すなわち、||mir(n1)-mir(n2)||nrm2pwr2、ここで、rn1およびrn2は、グリッド内の2つの隣接する位置に対応し、nrm2は、適切なノルムを表し、pwr2は、適切な指数を表す)の同時の最小化によって達成される。好ましくは、最小化関数fgridは、以下の形式、fgrid=ΣnΣj||P(rn,λj)-(ΦM(λj)+Σpi=1mir(n)Φi(λj))(c'1(rn)ε1(λj)+c'2(rn)ε2(λj))||nrm1pwr1+Σn,mλn,m||mir(n)-mir(m)||nrm2pwr2をとることになり、ここで、λn,mは、適切な重み係数である。
【００２８】
さらに、連立非線形方程式の逆算において、モデルパラメータmirの値が、限定された領域内にあるように制限されること(すなわち、limmini<mir<limmaxi、ここで、値limminiおよびlimmaxiは、シミュレーションまたは実験測定値によって定義される)、および/または、未知の値mirといくつかの事前定義された値
【数１】
0006669521.tif 000002
との間の差も、費用関数fgridの最小化と共に最小化されることが好ましい。好ましくは、これは、fgridと、モデルパラメータmirと事前定義された値
【数２】
0006669521.tif 000003
(事前情報(priors)とも称される)との間の重み付き差との線形結合を構成する第2の費用関数ggridの最小化によって実行され、すなわち、
【数３】
0006669521.tif 000004
ここで、nrm3は、適切なノルムを表し、pwr3は、適切な指数を表し、wirは、適切な重み係数を表す。好ましくは、事前情報
【数４】
0006669521.tif 000005
は、シミュレーションまたは実験測定値によって得られる。
【００２９】
本発明の上記および他の要素、特徴、特性、および利点は、添付図面を参照する好ましい実施形態の以下の詳細な説明からより明らかになるであろう。
【図面の簡単な説明】
【００３０】
【図１】マルチスペクトル光音響イメージングのためのデバイスの一例を示す図である。
【図２】波長依存光フルエンスの固有スペクトルモデルの一例を示す図であり、(a)〜(d)利用される波長λでのモデルスペクトルΦM(λ)、Φ1(λ)、Φ2(λ)、およびΦ3(λ)、(e)異なるモデル次元のためのトレーニングデータセットにおけるモデルの平均二乗誤差(MSE)、(f)〜(h)深さ(y軸)および組織酸素化(x軸)の関数としてのトレーニングデータセットにおけるパラメータm1、m2、およびm3の値である。
【図３】点のグリッド上の同時反転の場合のフルエンスパラメータm1、m2、およびm3の空間特性における制約の組み込みの一例を示す図であり、(a)点のグリッド上の有向グラフは、深さとともにm2モデルパラメータの値に制約を強制し、(b)点のグリッド上の非有向グラフは、隣接する点の間のフルエンスパラメータの大きな変動にペナルティを課す。
【図４】点のグリッド上の同時反転を用いる固有スペクトル法と、非否定的に制約された線形アンミキシング手法とを使用する、すべてのグリッド点に関するsO2推定誤差の一例を示す図である。
【図５】(a)担癌ヌードマウスのMSOT画像(1波長を提示)であり、画像内の点は、反転のために選択された初期グリッドに対応し、(b)反転のために使用されたピクセルにおけるSNRを最大化するための初期グリッドの自動微細化であり、(c)反転グリッドの凸包におけるフルエンスパラメータm2のマップであり、(d)フルエンス補正後の組織血液sO2の推定である。
【図６】マルチスペクトル光音響イメージングのための方法の処理およびアルゴリズムの一例における概略図である。
【発明を実施するための形態】
【００３１】
図1は、2つ以上の異なる照射波長λjおよび時間変化する強度の電磁放射線を物体1に照射するように構成された照射ユニット10を備えるマルチスペクトル光音響イメージングのためのデバイスの一例を示す。
【００３２】
たとえば、照射ユニット10は、たとえば、好ましくは、5〜20の励起波長の、690nmと900nmとの間の可視光および/または近赤外光を放射するように構成された少なくとも1つのレーザ光源を備えることができる。好ましくは、放射光は、パルス化、および/または、振幅または周波数変調される。
【００３３】
過渡的な、具体的にはパルス化および/または変調された電磁放射線を物体1に照射することによって、音響波が生成され、音響波は、検出ユニット20によって検出され、検出ユニット20は、1つまたは複数の超音波検出器、たとえば、1つまたは複数の圧電型検出器を備える。好ましくは、検出ユニット20は、超音波検出器要素の2次元または3次元アレイを備え、これによって、物体1内で生成された音響波は、検出され得、空間的に分解され得る。
【００３４】
物体1は、好ましくは、ヒトまたは動物の体の一部のような生体組織であり得る。所望の用途、前臨床または臨床用途の特殊性に応じて、ホルダ(図示せず)は、物体1、たとえば、組織サンプルまたは小動物を、照射ユニット10および検出ユニット20に対して保持し、位置決めするために設けられ得る。代替的にまたは追加的に、照射ユニット10および検出ユニット20は、ハンドヘルドデバイス内に統合され得、ハンドヘルドデバイスは、物体1上に手動で配置され得、および/または物体1に対して移動され得る。
【００３５】
検出器ユニット20において、検出された音響波は、処理ユニット30に転送される対応する検出器信号に変換され、処理ユニット30は、照射波長λjの各々で得られた検出器信号に基づいて、物体1の画像を再構成するように構成されている。
【００３６】
処理ユニット30は、さらに、異なる照射波長λjでの物体1の再構成された画像と、物体1における放射線フルエンスまたは放射線フルエンスの関数、たとえば、正規化された放射線フルエンスとを考慮して、物体1に含まれる少なくとも1つの放射線吸収体の相対濃度の空間分布を決定するように構成されており、ここで、放射線フルエンスまたは放射線フルエンスの関数は、放射線フルエンスの線形モデルに基づいている。線形放射線フルエンスモデルによれば、放射線フルエンスまたはその関数の波長依存基底関数を表す少なくとも2つのモデルスペクトルは、線形結合される。
【００３７】
好ましくは、表示ユニット40、たとえば、コンピュータモニタは、再構成された画像、および/または、物体1に含まれる少なくとも1つの放射線吸収体の相対濃度の決定された空間分布を表示するために設けられる。
【００３８】
以下では、物体1に含まれる少なくとも1つの放射線吸収体の相対濃度の空間分布を決定する好ましい態様が、詳細に説明される。
【００３９】
ピクセル単位の非線形スペクトルアンミキシング問題としての定量化問題の新しい定式化に基づく、マルチスペクトル光音響画像からの血液酸素化の定量的な尺度を求めるためのアルゴリズム的解法は、以下の「問題定式化」のセクションでより詳細に説明される。
【００４０】
好ましくは、波長依存基底関数の単純な線形モデルは、トレーニングデータを使用して組織内の光フルエンスのスペクトル摂動を捕捉するために作成され、これは、以下の「光フルエンスの固有スペクトルモデル」のセクションでより詳細に説明される。
【００４１】
好ましくは、このモデルに基づいて、非線形反転方式は、MSOT画像から定量的血液sO2推定値を得るために定式化され、これは、以下の「反転」のセクションで説明される。提案される反転方式は、さらに、反転精度を向上させるための光フルエンスの構造的特性およびスペクトル特性における事前知識を利用する。
【００４２】
最後に、生体内のマルチスペクトル光音響画像における新規のアルゴリズム的解法の適用は、以下の「組織画像への適用」のセクションで実証される。
【００４３】
問題定式化
光フルエンス、および正確なsO2推定の取得
主な組織吸収体である酸素化ヘモグロビン(HbO2)および脱酸素化ヘモグロビン(HHb)を用いて、組織の光音響画像P(r,λ)は、非線形関係によってHbO2およびHHbの濃度に関連付けられ得、すなわち、
【数５】
0006669521.tif 000006
である。式(1)では、C(r)は、グリュナイゼン係数(Γ(r))およびおそらくは追加の物理的現象(たとえば、超音波減衰)の効果と、システムの効果(たとえば、超音波センサの空間および周波数依存感度)とを含む空間依存スケーリング因数である。εHbO2(λ)およびεHHb(λ)は、ヘモグロビンの波長依存モル吸光係数であり、cHbO2(r)およびcHHb(r)は、位置rでの関連する濃度である。||Φ(r)||2は、位置rでのすべての波長にわたる光フルエンスのノルムであり、Φ(r,λ)/||Φ(r)||2は、特定の位置での光フルエンスの正規化された波長依存性であり、本発明に関連してΦ'(r,λ)とも称される項である。
【００４４】
空間のみに依存する因数C(r)および||Φ(r)||2は、HbO2およびHHbの相対濃度(それぞれ、c'HbO2(r)=Const(r)・cHbO2(r)およびc'HHb(r)=Const(r)・cHHb(r))がわかると、比率として計算されるsO2(式(2))の推定値に影響しない。
【００４５】
【数６】
0006669521.tif 000007
【００４６】
光フルエンスの波長依存性の推定値Φ'(r,λ)は、c'HbO2(r)およびc'HHb(r)の相対値を正確に抽出するために必要とされる。しかしながら、Φ'(r,λ)は、一般に未知の関数である。
【００４７】
本発明の好ましい態様によれば、光フルエンスの波長依存性Φ'(r,λ)は、波長の任意の関数ではなく、組織内のスペクトルパターンの限定されたサブ空間内にのみあると想定される。この態様は、光フルエンスのスペクトルが特定の組織吸収体による累積的な光吸収の結果であるので、光フルエンスのスペクトルは、常に、複雑な関係を介してこれらの吸収体のスペクトルに関連することになるという概念に基づいている。
【００４８】
別の好ましい態様によれば、スペクトルパターンのこのサブ空間は、低次元線形モデルを使用して正確に記述され得ることが想定される。「光フルエンスの固有スペクトルモデル」のセクションでは、固有スペクトルとも称される少数の特徴的組織光フルエンススペクトルの線形結合から成るそのようなモデルの作成が、詳細に説明される。
【００４９】
非線形スペクトルアンミキシング定式化
さらに別の態様によれば、組織内の各任意の位置rでのΦ'(r,λ)を正確に特徴付けるために、平均スペクトルΦM(λ)といくつかの光フルエンス固有スペクトルΦi(λ)との重ねあわせに基づいて線形モデルを定義することができることが想定され、次いで、位置rでのMSOT画像のための非線形混合モデルは、以下のように記述され得る。
【００５０】
【数７】
0006669521.tif 000008
【００５１】
P+2の未知数(すなわち、Pのモデルパラメータm1...Pr、および相対的血液濃度c'HbO2(r)およびc'HHb(r))についての解は、次いで、適切な非線形反転方式と、少なくともP+2の励起波長とを使用して与えられ得る。前述の通り、相対的血液濃度c'HbO2(r)およびc'HHb(r)は、スケーリング因数に関して元のcHbO2(r)およびcHHb(r)に比例するが、この事実は、レシオメトリック(相対)測定であるsO2の計算には影響しない。
【００５２】
前述の非線形スペクトルアンミキシング定式化を達成するために、(ΦM(λ),Φ1(λ),...,ΦP(λ))の形式の光フルエンスのための線形モデルが確立される。モデルは、組織内の任意の位置で、周囲の光学特性およびsO2値と独立に、光フルエンスΦ'(r,λ)の波長依存性を正確に記述することができる。
【００５３】
光フルエンスの固有スペクトルモデル
組織内の波長依存光フルエンスは、いくつかの異なるスペクトルパターンをとることができる。波長ごとの光フルエンスの変動は、ランダムではなく、典型的には、組織の深さ、組織散乱、血液sO2値、および全体的な血液濃度によって決定されるパターンに従う。好ましくは、トレーニングセットとも称される異なる光フルエンススペクトルパターンのセットに主成分分析(PCA)を使用することによって、少数のフルエンス固有スペクトルΦi(λ)から成る線形モデル(ΦM(λ),Φ1(λ),...,ΦP(λ))が決定される。PCAは、線形的にトレーニングセットの光フルエンススペクトルパターンのスペクトル変動を捕捉する際に、最小二乗誤差特性を提供するために使用される。
【００５４】
大きく変動し、ランダムに分散された光学特性および酸素化値を有する完全に不均質の媒体内で生成された光フルエンススペクトルパターンのセットにおいて固有スペクトルモデルをテストすることによって、固有スペクトルモデルが、すべてのそのような光フルエンススペクトルパターンをきわめて小さい平均二乗誤差(<2%)で正確に捕捉することができることが見出され、そのような条件(すなわち、光学特性の変動)とは無関係に、光フルエンスのスペクトルパターンは、固有スペクトルモデルによって正確にモデル化され得る限られたスペクトルサブ空間内にあることが確認された。
【００５５】
以下のサブセクションでは、本明細書で使用される光伝播モデル、ならびに、トレーニングデータから固有スペクトルモデルを作成するための手順の分析的な説明が詳細に説明される。
【００５６】
光伝播モデル
好ましくは、拡散方程式に基づく光伝播モデルが用いられ、ここで、波長依存光フルエンスモデル(ΦM(λ),Φ1(λ),...,ΦP(λ))は、無限の均質な媒体のための拡散方程式の解を使用して作成され、すなわち、
ΦZ=Φ0exp(-μeffZ)    (4)
であり、ここで、Zは、ソースからの距離(組織の深さ)を表し、μeffは、有効減衰係数、すなわち
【数８】
0006669521.tif 000009
であり、μS'は、減少した散乱係数であり、μαは、吸収係数である。
【００５７】
代替的には、光学特性の任意の構造のための拡散近似法の有限要素解法が使用され得る。
【００５８】
固有スペクトルモデル
好ましくは、光フルエンス固有スペクトルモデルを生成するためにトレーニングデータセットを生成するために、拡散方程式(式(4))の分析的解法に基づくシミュレーションが実行される。
【００５９】
たとえば、いくつかの光フルエンススペクトルパターンΦZ,OX(λ)は、10nmのステップサイズで波長λ=700nm〜890nmについて、0.143mmのステップサイズでz=0cmからz=1cmまでの深さについて、5%のステップでox=0%からox=100%まで変化する一定の酸素化における均一なヘモグロビン分布について計算される。μα=0.3cm-1(800nmにおけるヘモグロビンの等吸収点における)およびμS'=10cm-1(すべての波長に対して一定)の一定の光学特性が使用される。そのように計算された70×21の光フルエンス関数ΦZ,OX(λ)は、任意のスケーリングΦ'Z,OX(λ)=ΦZ,OX(λ)/||ΦZ,OX||2を無視するために正規化された。1cmの組織深さは、小動物の光音響イメージングに用途を関連付けるために選択された。
【００６０】
計算されたベクトルΦ'Z,OX(λ)は、以下で、図2(a)〜図2(d)に例示的に示すように、固有スペクトルの線形モデル(平均スペクトルΦM(λ)および3つの固有スペクトルΦ1(λ)、Φ2(λ)、Φ3(λ))を導出するために、PCAの文脈でトレーニングデータとして使用される。
【００６１】
PCAは、最小の平均二乗誤差で、少数の第1の係数におけるトレーニングデータセットにおいて使用されるすべてのスペクトルパターンの支配的な変動を捕捉する。3つの固有スペクトルのモデル寸法は、関連する誤差が小さいより単純なモデルを提供するために選択された(図2(e)参照)。
【００６２】
トレーニングデータセットにモデルのパラメータm1、m2、およびm3を関連付け(図2(f)〜図2(h)をそれぞれ参照)、驚くべきことに、m2の値が組織の深さに大きく関連付けられ得、m1の値が背景の組織酸素化の平均レベルに大きく関連付けられ得ることが見出された。これは、第2の固有スペクトルΦ2(λ)が、深さによる光フルエンスの変更と、周囲の組織の平均光学特性とに主に関連付けられ、第1の固有スペクトルΦ1(λ)が、周囲の組織の平均酸素化に関連する光フルエンスの「スペクトル形状」に関連付けられることを示す。さらに、トレーニングデータセットにおけるこれらのパラメータの最高値および最低値は、組織において期待されるそれらの限界[limmini,limmaxi]の指標を提供することができる。
【００６３】
任意の組織において作成されたシミュレーションベースの光フルエンススペクトルに対する固有スペクトルモデル(ΦM(λ),Φ1(λ),...,ΦP(λ))の精度の検証は、わずかな、たとえば3つの固有スペクトルのみを有する単純な線形モデルが、光学特性、深さ、および血液酸素化と無関係に、複雑な組織構造内のΦ'(r,λ)のスペクトル変動を捕捉することができる本発明の手法を支持する、非常に小さいモデル誤差をもたらした。
【００６４】
同様に、組織内に存在する光フルエンススペクトルパターンを捕捉する固有スペクトルモデルの有効性および精度は、生体内および生体外での小動物からの測定値を用いて実験的に調査され、ここで、組織内の光フルエンスは、十分に特徴付けされたスペクトルを有する基準発色団を深部組織内に挿入することによって測定された。この調査の結果は、生体内で得られたノイズの多い実験データの場合であっても、固有スペクトルモデルが、組織内の光フルエンスのスペクトル変動を現実的に捕捉することができることを示した。
【００６５】
反転
光フルエンスのための固有スペクトルモデルを考えると、MSOTを使用する血液sO2の定量的な推定は、非線形最適化手法を使用して式(3)を解くことに帰着する。反転は、画像内の単一のピクセルに対して独立して実行され得、または、画像領域において適切に分布したピクセルのグリッドに対して同時に実行され得る。第2の場合には、固有スペクトルモデルパラメータの値に対する追加の制約が、光フルエンスの空間的特性およびスペクトル特性に基づいて課され得る。
【００６６】
以下では、3つの異なる反転手法が提示され、すなわち、(1)局所的または大域的な最適化を使用するピクセル単位の反転、(2)点のグリッドに対する制約された最適化であって、光フルエンスの空間特性がさらにモデル化される、および、(3)別の制約された最適化であって、光フルエンスの空間特性とスペクトル特性の両方が、反転の安定性を高めるために考慮される、が提示される。
【００６７】
(1)単一ピクセル反転
単一のピクセルの反転の場合、特定の点rのための最適化問題は、式(6)に従って制約される、m1r、m2r、m3r、c'HbO2(r)、およびc'HHb(r)の値に対する、式(5)における関数の最小化として定式化する。
【００６８】
【数９】
0006669521.tif 000010
【００６９】
シミュレートされたデータセットの場合において式(5)および(6)の非線形最適化問題を解くために制約された局所最適化アルゴリズムを使用して、問題は、非凸状(non-convex)であり、したがって、反転は、小さい割合のピクセルのみに正確に収束することが観察され得る。したがって、最適化問題は、いわゆる不良設定問題である。好ましくは、反転の精度をさらに向上させ、唯一の解を達成するために、さらなる制約が、フルエンスパラメータmirの値に課される。
【００７０】
(2)空間フルエンス特性を組み込む点のグリッドにおける反転
好ましくは、光フルエンスの空間特性は、さらに、最適化問題の不良設定性を克服するために利用される。任意に変化し得る組織吸収に反して、光フルエンスは、拡散光伝播の性質のために空間内で必ずなめらかに変化する。
【００７１】
固有スペクトルモデル反転の好ましい態様によれば、そのような事前情報は、画像領域内の点のグリッドに対して同時の反転を試み、隣接するピクセル間のフルエンスモデルパラメータm1r、m2r、m3rの大きな変動にペナルティを課すことによって組み込まれる。これは、好ましくは、1つまたはすべてのモデルパラメータの変動を抑制するために非線形反転関数に適切なラグランジュ乗数λiを組み込むことによって達成される。ラグランジュ乗数の値は、シミュレートされたデータセットに交差検証を使用して最適化され得る。加えて、第2の固有スペクトルは、深さに強く関連付けられるので、m2モデルパラメータは、(図2(g)に従って)組織深さが増加するにつれて、より低い値をとるように制約され得る。
【００７２】
図3は、点のグリッドに対する同時反転の場合のフルエンスパラメータm1r、m2r、およびm3rの空間特性への制約の組み込みを示し、ここで、点のグリッドにおける有向グラフは、深さとともにm2モデルパラメータの値への制約を強制し(図3(a)参照)、点のグリッドにおける非有向グラフは、隣接する点の間のフルエンスパラメータの大きい変動にペナルティを課す(図3(b)参照)。
【００７３】
合計P×Lの点rp,lのPの円弧およびLの放射状線の円形グリッドを想定し(図3参照)、ベクトルmi=[mir1,1,mir1,2,...,mir1,L,mir2,1,...,mirp,l,...,mirP,L]をすべてのそのような点にわたる光フルエンスパラメータi(i=1...3)の値に対応させると、新しい反転問題は、好ましくは、式(8)で定義された制約の下で、式(7)で定義された関数fgridの最小化として定義される。
【００７４】
【数１０】
0006669521.tif 000011
【００７５】
式(7)では、Wは、想定された点のグリッドに対応する重み付けされた連結性行列である。各行列要素は、一対のグリッド点rp1,l1rp2,l2に対応し、点が直接接続されていない場合、ゼロであり、または、点が接続されている場合、それらの距離(w(rp1,l1rp2,l2)∝1/||rp1,l1-rp2,l2||)に反比例する。
【００７６】
点のグリッドにおける同時反転は、驚くべきことに、反転の安定性を有意に高めることが判明している。
【００７７】
空間-スペクトルフルエンス特性を組み込む点のグリッドにおける反転
反転の安定性をさらに高めるために、同様の条件における光フルエンスの以前の実験的な測定値に基づいて、または、再構成されたマルチスペクトル光音響画像から直接導出され得る代替の事前知識に基づいて、追加の制約がフルエンスパラメータの値に課せられ得る。グリッド反転の文脈においてそのような制約を強制する2つの可能な例が、以下で説明される。
【００７８】
m1に対する大域的制約
たとえば、生体内イメージングの特定の場合には、好ましくは、平均組織酸素化が50%よりも高いことが仮定され得、したがって、モデルパラメータm1は、図2(f)によって示されるように、ゼロ以下の値を取得することになる。この好ましい仮定は、複数の生体内対照実験(動物が100%O2ガスを呼吸しており、画像化が脳内または腹部内実行される)によって検証されており、ここで、測定されたフルエンスパラメータm1は、負であると繰り返し見出された。この条件下で、かつ、すべてのグリッド点において負であるようにm1rの値を制約することによって、優れた収束速度が、高いレベルの重畳ノイズの場合にも達成され得る。仮定は、グリッド点の特定の酸素化値に行われず、平均組織酸素化のみ行われることに留意されたい。
【００７９】
モデルパラメータに対する局所的制約
代替的には、事前のフルエンス測定値が利用可能でない場合、モデルパラメータm1のためのいくつかの制約が、利用可能なデータ自体から抽出され得る。データに第1の線形アンミキシング手法を実行することによって、sO2の第1の推定マップが得られ得る。このsO2マップは、組織深さとともに漸次的に誤りを含むが、第1の近似として役割を果たすことができることに留意されたい。このsO2マップおよび平均組織光学特性(たとえば、800nmにおけるμα=0.2cm-1、およびμS=10cm-1)を使用し、組織内の光フルエンスは、シミュレートされ得、1つまたはすべてのモデルパラメータm1r、m2r、およびm3rの値に対する第1の推定(すなわち、
【数１１】
0006669521.tif 000012
)は、抽出され得、ここで、rは、グリッド点の位置に対応する。以下では、1つまたはすべてのモデルパラメータm1r、m2r、およびm3rの値が初期推定値
【数１２】
0006669521.tif 000013
に近い領域に存在するように制約される違いのみで、式(7)〜(8)の最適化問題が解かれる。これらの領域は、深い組織では元の推定値が真の値から通常は逸脱するので、深さとともに漸次的により大きくされる必要があることに留意されたい。各グリッドピクセルに対応するモデルパラメータのための事前情報を組み込むそのような最適化手法は、また、ベイズ反転の文脈において定式化され得る。そのような方法で問題を制約することによって、反転の安定性は、また、高いランダムノイズの場合でも大幅に向上される。
【００８０】
図4は、シミュレートされたデータセットの場合における64の点のグリッドにおける式(7)および(8)の反転の一例に関係し、本発明の方法を使用したとき(青い曲線)と、単純な線形の負に制約されないアンミキシングを使用したとき(赤い曲線)の、すべてのグリッド点に関する局所的血液sO2の推定値における誤差を示す。比較から明らかなように、本発明の方法は、しばしば近似値として使用される線形アンミキシングと比較して、血液酸素化の推定値におけるきわめて小さい誤差に関連付けられる。
【００８１】
組織画像への適用
以下では、本発明の固有スペクトル方法は、実験MSOTデータの場合に例示的に実証される。同所性乳腺4T1腫瘍を有するヌードマウスは、690nmから900nmまでの8つの励起波長において腫瘍領域において画像化された。
【００８２】
図5(a)は、波長900nmにおいて取得された解剖学的光音響画像を示す。画像上の点は、上記で説明した非線形反転方式の適用のために選択されたピクセルの初期グリッドに対応する。グリッドは、腫瘍領域をカバーするように選択された。ピクセルのグリッドは、その後、反転において低い画像強度(かつ、したがって、低いSNR)のピクセルを使用することを回避するように自動的に調整される。
【００８３】
非線形反転方式を点の微細化された点のグリッド(図5(b)参照)に適用した後、フルエンスパラメータm1r、m2r、およびm3rの推定値は、各グリッド点rについて得られる。以下では、グリッドの凸包内の各ピクセルのためのフルエンスパラメータは、キュービック補間を使用して得られ、したがって、画像サブ領域のためのフルエンスマップがもたらされる。
【００８４】
図5(c)は、フルエンスパラメータm2のためのそのようなマップを提示する。光フルエンスの波長依存性は、Φ'(r,λ)=ΦM(λ)+m1rΦ1(λ)+m2rΦ2(λ)+m3rΦ3(λ)のように、それらのマップ内の各ピクセルについて計算され、ここでΦi(λ)は、i番目のフルエンス固有スペクトルである。
【００８５】
最後に、フルエンス補正されたMSOT画像は、元の画像P(r,λ)を各々の位置rおよび波長λにおける正規化された波長依存光フルエンスΦ'(r,λ)で除算した後、得られ、すなわち、Pcorr(r,λ)=P(r,λ)/Φ'(r,λ)となる。フルエンス補正の後、オキシヘモグロビンおよびデオキシヘモグロビンのスペクトルを用いた非負的に制約された線形スペクトルアンミキシングは、組織血液酸素化の推定値をもたらすことができる。
【００８６】
図5(d)は、固有スペクトル法を使用して計算された腫瘍の領域内の推定組織酸素化を示す。組織酸素化は、グレースケールで表された解剖学的画像に疑似カラーでオーバレイされている。
【００８７】
図6は、マルチスペクトル光音響イメージングのための方法のプロセスおよびアルゴリズムの一例の概要を示す。本例では、酸素化ヘモグロビン(HbO2)および脱酸素化ヘモグロビン(HHb)は、物体1内の主な組織吸収体である。
【００８８】
第1のステップでは、MSOT画像は、照射波長λjの過渡電磁放射線を物体1に照射する際に物体1内で生成された検出された音響波に基づいて、異なる照射波長λjにおいて再構成される(図1参照)。
【００８９】
第2のステップでは、再構成されたMSOT画像は、点選択プロセスを受け、ここで、画像領域における信頼できるピクセルriは、半自動的な方法で、またはユーザの注釈によって選択される。
【００９０】
すべての照射波長λjに対する選択されたピクセルriにおける光音響画像強度P(ri,λj)は、第3のステップで処理され、ここで、フルエンスパラメータmkriおよび濃度値、たとえば、sO2値は、画像強度P(ri,λj)と、物体内のフルエンスの所定のモデルとを使用して、グリッド内の各ピクセルriについて計算され、ここで、モデルは、モデルスペクトルΦk(λ)の線形結合に基づく。これは、以下の式を解く反転アルゴリズムによって実行される。
【００９１】
【数１３】
0006669521.tif 000014
【００９２】
計算されたフルエンスパラメータmkri(i=1〜Np、Npは、選択されたピクセルの数であり、k=1〜K、Kは、モデルスペクトルの数である)は、補間によって、グリッドの凸包に対応する全領域内の光フルエンスの波長依存性を推測するために使用され得、すなわち、mk値は、点のグリッドの凸包内のすべての中間ピクセルについてmkriから補間される。
【００９３】
続いて、光フルエンスの波長依存性は、Φ'(r,λ)=ΦM(λ)+m1rΦ1(λ)+m2rΦ2(λ)+m3rΦ3(λ)+...のように、凸包内の各ピクセルについて計算される。
【００９４】
フルエンス補正されたMSOT画像は、元の再構成されたMSOT画像P(r,λ)を各々の位置rおよび波長λにおける正規化された波長依存光フルエンスΦ'(r,λ)で除算した後、得られ、すなわち、Pcorr(r,λ)=P(r,λ)/Φ'(r,λ)となる。
【００９５】
フルエンス補正の後、オキシヘモグロビンおよびデオキシヘモグロビンのスペクトル内の非負的に制約されたスペクトルアンミキシングは、特に正確なsO2マップの決定をもたらす。
【００９６】
好ましくは、フルエンスのモデル、具体的には、モデルスペクトルΦk(λ)は、具体的にはあらかじめ、モデル作成アルゴリズムによる現在の測定値と独立して確立され、ここで、物体内の光フルエンスの波長依存性の基底関数は、用いられた励起波長λjのために作成される。好ましくは、以前のフルエンスモデル化方法と異なり、光フルエンスの波長依存性のみがモデル化され、実際のフルエンスは、絶対的にはモデル化されない。さらに、モデル化は、統計的処理をトレーニングデータに用いることによって、MSOT画像に適した空間光伝播モデルを介さずに実行される。
【００９７】
総括
本発明は、組織酸素化の定量的マルチスペクトル光音響イメージングの問題のための新規な解決策に関する。本発明の好ましい態様によれば、組織内の波長にともなう光フルエンスのスペクトル変動を正常に捕捉するための単純な線形モデルが作成され、検証された。このモデルを使用し、光フルエンスを補正し、組織血液酸素化の正確な定量的推定値を得るための非線形反転方式が定式化された。MSOTの定量化問題を、モデルベースの光学特性推定手法の代わりに、ピクセル単位の非線形アンミキシング問題として定式化することは、組織の実験の画像におけるより高い柔軟性およびより容易な適用可能性を可能にする。本発明の方法は、完全かつ正確な空間情報を想定せず、画像再構成においてしばしばモデル化されない多数の現象(たとえば、超音波アレイの空間感度領域、超音波減衰および散乱、グリュナイゼン係数)による誤差の傾向がある絶対画像値に依存しない。この点で、本発明の方法は、生体内MSOT画像において容易に適用可能である。さらに、提示される反転問題は、低次元のもの(典型的には、グリッド反転の場合に関して300〜500の未知のパラメータ)であるので、計算コストは、かなり低い。具体的には、Intel Core i7-4820Kプロセッサを使用する64点のグリッドの反転に必要な計算時間は、1分未満であり得る。
【００９８】
本発明の手法の検証に基づいて、むしろ単純な3次元モデルは、周囲の組織の光学特性および生理機能と無関係に、光フルエンスのスペクトル変動を正確に捕捉するために十分であると判明した。
【００９９】
基本的には、ヘモグロビンの影響が近赤外線領域においてはるかにより優勢であると仮定して、組織内に同様に存在する水、脂質、およびメラニンの吸収は、適切な波長選択の下で無視され得る。さらに、散乱の波長依存性は、無視され得る。しかしながら、本発明の好ましい実施形態によれば、そのようなパラメータの影響も、固有スペクトルモデルを最適化するために同様に組み込まれ得る。
【０１００】
前述の単純化している仮定にもかかわらず、マウスの生体内および生体外での実験の光音響画像の場合に適したモデルは、理論と実験結果との間で非常に有望な一致を示した。
【０１０１】
その最も単純な形態(ピクセル単位の反転)での反転問題は、非凸状かつ不良設定であるので、ピクセル単位の反転の場合には、光フルエンスおよび局所的酸素化の異なる組み合わせは、乗算後、ほぼ同一の結果を生成する可能性がある。したがって、非線形反転は、追加の制約が適用されないかぎり、間違った解に収束する可能性がある。
【０１０２】
したがって、利用可能な事前知識を使用して、空間-スペクトルフルエンス特性をより複雑な反転方式に組み込むことが、特に好ましい。このようにして、反転安定性は、スペクトルがランダムノイズによって破損している場合にも達成された。好ましい反転方式は、近似として実験研究において今日典型的に使用されている線形スペクトルアンミキシングのものよりもはるかに低い推定誤差で、組織酸素化の正確な推定値を提供する。
【符号の説明】
【０１０３】
1  物体
10  照射ユニット
20  検出ユニット
30  処理ユニット
40  表示ユニット
【図１】
0006669521.tif 000015
【図２】
0006669521.tif 000016
【図３】
0006669521.tif 000017
【図４】
0006669521.tif 000018
【図５】
0006669521.tif 000019
【図６】
0006669521.tif 000020