特表2024-529002 | 知財ポータル「IP Force」

知財求人 - 知財ポータルサイト「IP Force」

▶ ユニヴェルシテグルノーブルアルプの特許一覧 ▶ センター　ナショナル　ド　ラ　ルシェルシュ　サイエンティフィークの特許一覧

特表2024-529002マルチスペクトル光音響イメージングを介した生体組織における変動の定量分析

書誌要約請求の範囲詳細な説明課題実施例実施するための形態図面の説明

目に優しい文字サイズ小中大 PDF Top

< >

1
2
3A
3B
4A
4B
5A
5B

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公表特許公報(A)

(11)【公表番号】

(43)【公表日】2024-08-01

(54)【発明の名称】マルチスペクトル光音響イメージングを介した生体組織における変動の定量分析

(51)【国際特許分類】

A61B 8/13 20060101AFI20240725BHJP

【ＦＩ】

A61B8/13

【審査請求】未請求

【予備審査請求】未請求

(21)【出願番号】P 2024506758

(86)(22)【出願日】2022-07-20

(85)【翻訳文提出日】2024-04-01

(86)【国際出願番号】 FR2022051444

(87)【国際公開番号】W WO2023012416

(87)【国際公開日】2023-02-09

(31)【優先権主張番号】2108452

(32)【優先日】2021-08-03

(33)【優先権主張国・地域又は機関】FR

(81)【指定国・地域】

(71)【出願人】

【識別番号】323002277

【氏名又は名称】ユニヴェルシテグルノーブルアルプ

(71)【出願人】

【識別番号】506316557

【氏名又は名称】サントルナショナルドゥラルシェルシュシアンティフィック

(74)【代理人】

【識別番号】100108453

【弁理士】

【氏名又は名称】村山靖彦

(74)【代理人】

【識別番号】100110364

【弁理士】

【氏名又は名称】実広信哉

(74)【代理人】

【識別番号】100133400

【弁理士】

【氏名又は名称】阿部達彦

(72)【発明者】

【氏名】バスティアン・アルナル

(72)【発明者】

【氏名】エマニュエル・ボッシー

(72)【発明者】

【氏名】ギヨーム・ゴドフロワ

【テーマコード（参考）】

4C601

【Ｆターム（参考）】

4C601DE16

4C601EE04

4C601EE09

4C601JB28

4C601JB34

4C601JB48

(57)【要約】

光音響画像を処理するための方法が開示され、この方法は、試料の時系列画像を取得するステップ(110)であって、前記画像がMx個の励振パルス波長において光音響イメージングシステムによって獲得され、波長当たりN個の画像が獲得される、ステップと、N^*Mx個のフィルタリングされた画像を取得するために、N^*Mx個の獲得された画像のすべてに適用される特異値分解を介して、マルチスペクトル時空間フィルタリングを実行するステップ(120)と、各波長に対して、N個のフィルタリングされた画像に基づいてフィルタリングされた分散画像を計算するステップ(140)であって、フィルタリングされた分散画像における座標rの画素が、この波長に対して取得されたフィルタリングされた画像における同じ座標rの画素の値の分布の分散に等しい、ステップと、マルチスペクトル時空間フィルタリングが実行された後の残余電子ノイズの分散を減算することによって、フィルタリングされた分散画像を訂正するステップ(150A)であって、前記ノイズが光音響イメージングシステムの超音波のセンサによって生み出される、ステップとを含む。

【特許請求の範囲】

【請求項1】

光音響画像を処理するための方法であって、
- M_λ個の励振パルス波長に対して光音響イメージングシステムによって獲得される時間的に連続する試料の画像を取得するステップ(110)であって、波長当たりN個の画像が獲得される、ステップと、
- N^*M_λ個のフィルタリングされた画像を取得するために、獲得されたすべてのN^*M_λ個の画像に適用される特異値分解による、マルチスペクトル時空間フィルタリング(120)するステップと、
- 各波長に対して、前記N個のフィルタリングされた画像からフィルタリングされた分散画像を計算するステップ(140)であって、前記フィルタリングされた分散画像における座標rの1つの画素が、その波長に対して取得された前記フィルタリングされた画像における同じ座標rの画素値の分布の分散に等しい、ステップと、
- 各波長に対して、前記光音響イメージングシステムの超音波センサによって生み出される、前記マルチスペクトル時空間フィルタリングの後の残余電子ノイズの分散を減算することによって、前記フィルタリングされた分散画像を訂正するステップ(150A)と
を含む、方法。

【請求項2】

各波長に対して、訂正された変動画像の決定を含み、前記訂正された変動画像の各画素が、考慮される前記波長に対する前記減算によって取得される訂正された分散画像の対応する画素の平方根に等しい、請求項1に記載の方法。

【請求項3】

前記画像上の前記超音波センサによって生み出された前記残余電子ノイズの前記分散を推定するステップをさらに含み、前記画像上の前記超音波センサによって生み出された前記残余電子ノイズの前記分散が、特異値への分解による前記マルチスペクトル時空間フィルタリングにより除去されたノイズの量により訂正された、試料の不在時に獲得された光音響信号において生み出された電子ノイズの分散の関数として推定され、除去されたノイズの前記量が、特異値分解による前記マルチスペクトル時空間フィルタリングにより取り除かれた最低のエネルギー成分に対応する特異値に基づいて推定される、請求項1または2に記載の方法。

【請求項4】

考慮される各波長に対して、前記試料による吸収変動を表す吸収変動画像を取得するために、前記光音響イメージングシステムのレーザーパルスフルエンスの関数により訂正された前記変動画像を正規化するステップ(150B)を含む、請求項1から3のいずれか一項に記載の方法。

【請求項5】

特異値分解による前記マルチスペクトル時空間フィルタリング(120)が、取り除かれるべき最高エネルギーの特異値に対応する成分を選択するステップと、選択された成分を取り除くステップとを含み、前記選択が、インデックス値のセットから、コントラスト対ノイズ比が最大である、保たれるべき第1の成分を特定するインデックスを選ぶことによって行われ、あるインデックスに対して決定される前記コントラスト対ノイズ比が、保たれるべき前記第1の成分を特定するために前記インデックスを適用する特異値への分解によるマルチスペクトル時空間フィルタリング(120)によって計算されるフィルタリングされた分散画像に対して決定される、請求項1から4のいずれか一項に記載の方法。

【請求項6】

前記コントラスト対ノイズ比が、少なくとも1つの波長に対して獲得された前記画像の平均値によるコントラストを除去することによって決定される、請求項5に記載の方法。

【請求項7】

少なくとも2つの対応する波長に対して取得された少なくとも2つの吸収変動画像から、酸素飽和度の画像を計算するステップ(160)を含む、請求項4に従属する場合の請求項4から6のいずれか一項に記載の方法。

【請求項8】

前記酸素飽和度の画像の前記計算(160)が、座標rの各画素に対する、総ヘモグロビン濃度と、酸素飽和度と、前記吸収変動画像の画素rにおける値との関係を表現するモデルに基づいて実行される、請求項7に記載の方法。

【請求項9】

プログラムコード命令を備える少なくとも1つのデータメモリと、少なくとも1つのデータプロセッサとを備える、光音響画像処理デバイスであって、前記データプロセッサが、前記プログラムコード命令が前記データプロセッサによって実行されると、請求項1から8のいずれか一項に記載の方法を前記光音響画像処理デバイスに実行させるように構成される、光音響画像処理デバイス。

【請求項10】

プロセッサによって実行されると、請求項1から8のいずれか一項に記載の方法を実行するコンピュータプログラム命令を含む、コンピュータ可読記憶媒体。

【請求項11】

プロセッサによって実行されると、請求項1から8のいずれか一項に記載の方法を実行するコンピュータプログラム命令を備える、コンピュータプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本文書は、音響分解能光音響イメージングに関し、より正確には、マルチスペクトル光音響イメージングシステムおよび関連するデバイスにより獲得される画像を処理する方法に関する。

【背景技術】

【0002】

光音響イメージング(オプトアコースティックとも呼ばれる)の技法は、典型的には300nmから2000nmのスペクトルの範囲内にある、特に650nmから950nmの区間にあるレーザーパルスである電磁放射を試料(典型的には生体組織)に照射することによる励起によって、試料で生み出される超音波の生成に基づく。一般に、光音響画像は、送られた光が試料に吸収されることにより生成される音響信号の測定から再構築される。光音響的な生成の効率に関する理由で、イメージングされることになる生体組織を照射するために、一般に、短いナノ秒のパルスが医用イメージングに使用される。パルス光の吸収は、軟組織において温度の急激な上昇をもたらし、これは熱弾性効果により圧力の上昇を生み出し、そして、この圧力は、生体組織においてパルス音響波の伝播を引き起こすことにより緩和する。

【0003】

音響分解能光音響イメージングは、超音波センサのアレイが使用されるという点で、光分解能光音響イメージングと異なる。音の速さを知っていれば、音源の大きさと位置を示す画像を、受信された音響信号から再構築することができる。一方、光分解能光音響イメージングでは、試料の表面上の点ごとに各光パルスが合焦され、光合焦ゾーンから放出された音が単一の超音波センサによって測定される。

【0004】

音響分解能光音響イメージングは、生体組織における詳細な光吸収コントラストを提供する医用イメージング技法である。たとえば、血液の光音響イメージングが、血管をイメージングするために使用される。血液は実際に、ヘモグロビンを介して、人体において非常に豊富な吸収要素である。

【0005】

しかしながら、超音波を捉えるために圧電センサアレイを使用する固定された光音響イメージングデバイスは可視性が限られており、センサアレイがイメージングされる媒体を取り囲まないとき、帯域幅と開口数に関するセンサの制約により、ある組織構造が獲得された画像に現れない。したがって、開口数と帯域幅が限られている獲得システムの使用は、限られた可視性というアーティファクトを引き起こし、すなわち、その軸がプローブの軸とほぼそろっている血管は見えず、大きな血管も見えない。より具体的には、90°側で、およびその周りで受信された波が入手可能ではなく、概ね水平方向の血管構造のみが検出され、垂直方向の構造は現れない。

【先行技術文献】

【非特許文献】

【0006】

【非特許文献1】「Photoacoustic fluctuation imaging: theory and application to blood flow imaging」、Vilov S.、Godefroy G.、Arnal B.、およびBossy E.、Optica, 7(11)、1495-1505(2020)

【発明の概要】

【発明が解決しようとする課題】

【0007】

そのような問題を解決するために、断層映像法を使用し、広帯域検出器を使用することが可能である。しかしながら、そのような非共鳴検出器は、超音波イメージングの実行を可能にせず、関連する断層映像法は、すべての臨床上および前臨床の問題には適合しない。

【0008】

可視性の問題を克服するために、最近では深層学習を用いた人工知能に基づく手法も提案されているが、そのような方策には、訓練データの利用可能性ならびに一般化という問題がある。

【0009】

したがって、より適応された生体組織イメージングの方策が必要である。

【課題を解決するための手段】

【0010】

第1の態様によれば、光音響画像を処理する方法は、M_λ個の励振パルス波長に対して光音響イメージングシステムによって獲得される時間的に連続する試料の画像を取得するステップであって、波長当たりN個の画像が獲得される、ステップと、N^*M_λ個のフィルタリングされた画像を取得するために、獲得されたすべてのN^*M_λ個の画像に適用される特異値分解による、マルチスペクトル時空間フィルタリングと、各波長に対する、N個のフィルタリングされた画像からフィルタリングされた分散画像(variance image)の計算であって、フィルタリングされた分散画像における座標rの1つの画素が、前記波長に対して取得されたフィルタリングされた画像における同じ座標rの画素値の分布の分散に等しい、計算と、各波長に対する、光音響イメージングシステムの超音波センサによって生み出される、マルチスペクトル時空間フィルタリングの後の残余電子ノイズの分散を減算することによる、フィルタリングされた分散画像の訂正とを含む。このように取得される訂正された分散画像は、各画像点において、関係する画像点に対応する空間点において吸収されたエネルギーの二乗に比例する。

【0011】

1つまたは複数の実施形態において、方法は、各波長に対して、訂正された変動画像(image of fluctuations)の決定を含み、訂正された変動画像の各画素は、考慮される波長に対する前記減算によって取得される訂正された分散画像の対応する画素の平方根に等しい。このように取得される訂正された変動画像は、各画像点において、関係する画像点に対応する空間の点において吸収されたエネルギーに比例する。

【0012】

1つまたは複数の実施形態では、方法はさらに、画像上の超音波センサによって生み出された残余電子ノイズの分散の推定を含み、画像上の超音波センサによって生み出された残余電子ノイズの分散は、特異値への分解によるマルチスペクトル時空間フィルタリングにより除去されたノイズの量により訂正された、試料の不在時に獲得された光音響信号において生み出された電子ノイズの分散の関数として推定され、除去されたノイズの量は、特異値分解によるマルチスペクトル時空間フィルタリングにより取り除かれた最低のエネルギー成分に対応する特異値に基づいて推定される。

【0013】

1つまたは複数の実施形態では、方法は、考慮される各波長に対して、試料による吸収変動を表す吸収変動画像を取得するために、光音響イメージングシステムのレーザーパルスフルエンスの関数により訂正された変動画像の正規化を含む。

【0014】

1つまたは複数の実施形態では、特異値分解によるマルチスペクトル時空間フィルタリングは、取り除かれるべき最高エネルギーの特異値に対応する成分の選択と、選択された成分を取り除くステップとを含み、この選択は、インデックス値のセットから、コントラスト対ノイズ比が最大である、保たれるべき第1の成分を特定するインデックスを選ぶことによって行われ、あるインデックスに対して決定されるコントラスト対ノイズ比は、そのインデックスを適用する特異値への分解によるマルチスペクトル時空間フィルタリングによって計算されるフィルタリングされた分散画像から計算される。

【0015】

1つまたは複数の実施形態では、コントラスト対ノイズ比は、少なくとも1つの波長に対する獲得された画像の平均値によるコントラストを除去することによって決定される。

【0016】

1つまたは複数の実施形態では、方法は、少なくとも2つの対応する波長に対して取得された少なくとも2つの吸収変動画像からの酸素飽和度画像の計算を含む。この計算は、座標rの各画素に対する、総ヘモグロビン濃度と、酸素飽和レベルと、吸収変動画像の画素rにおける値との関係を表現するモデルに基づいて実行され得る。

【0017】

別の態様によれば、光音響画像処理デバイスは、プログラムコード命令を備える少なくとも1つのデータストレージと、少なくとも1つのデータプロセッサとを備え、データプロセッサは、プログラムコード命令がデータプロセッサによって実行されると、実施形態のいずれかによる光音響画像処理方法を光音響画像処理デバイスに実行させるように構成される。

【0018】

別の態様では、コンピュータ可読データストレージは、プロセッサによって実行されると、実施形態のいずれかによる光音響画像処理方法を実行するコンピュータプログラム命令を備える。

【0019】

別の態様によれば、コンピュータプログラムは、プロセッサによって実行されると、実施形態のいずれかによる光音響画像処理方法を実行するコンピュータプログラム命令を備える。

【0020】

同封の図面を参照して、限定するものとしてではなく例示として与えられた実施形態および例に基づいて行われる以下の詳細な説明から、他の特徴および利点が得られる。

【図面の簡単な説明】

【0021】

【図1】マルチスペクトル光音響イメージングシステムによって獲得される画像を処理するための方法の簡略化されたフローチャートである。

【図2】マルチスペクトル光音響イメージングシステムによって獲得される画像を処理するための方法の態様を示す図である。

【図3A】マルチスペクトル光音響イメージングシステムによって獲得される画像を処理するための方法において使用され得る特異値分解によるマルチスペクトル時空間フィルタリング方法の態様を示す図である。

【図3B】マルチスペクトル光音響イメージングシステムによって獲得される画像を処理するための方法において使用され得る特異値分解によるマルチスペクトル時空間フィルタリング方法の態様を示す図である。

【図4A】特異値分解によるマルチスペクトル時空間フィルタリングの間に使用され得る特異値選択方法の簡略化されたフローチャートである。

【図4B】特異値分解によるマルチスペクトル時空間フィルタリングの間に使用され得る特異値選択方法の態様を示す図である。

【図5A】本説明に従って獲得された画像および従来の方法に従って獲得された画像を処理するための方法によって取得される酸素化画像の2つの例を示す図である。

【図5B】本説明に従って獲得された画像および従来の方法に従って獲得された画像を処理するための方法によって取得される酸素化画像の2つの例を示す図である。

【発明を実施するための形態】

【0022】

レーザー励起によりイメージングされるべき試料(一般には生体組織)において生み出される超音波の獲得のために音響波センサのアレイを使用する、マルチスペクトル光音響イメージングシステムによって獲得される画像を処理するための方法が、より詳しく説明される。

【0023】

この方法は、生体組織において発生するあらゆるタイプの変動の光音響イメージングに適用可能である。血液の酸素化率の定量分析へのこの方法の適用が、非限定的な例として説明される。

【0024】

本文書の枠組みでは、変動画像とは各画像点における時間的変動を表す画像を指し、各画像点(または画素)は考慮される変動の標準偏差を表す。したがって、変動画像は時間的に連続する画像から取得することができ、変動画像における座標rの各画素は、時間的に連続する画像における同じ座標rの画素の統計的分布の標準偏差から計算される。分散画像は各変動画像に対応し、二乗された変動画像を指定する。分散画像の各画像点は、考慮される変動の分散、すなわち、考慮される時間的に連続する画像における同じ座標の画素の統計的分布の分散を表す。文脈に応じて、分散画像または変動画像について論じる。

【0025】

本文書の枠組みでは、画像は2Dまたは3D画像であり得る。簡潔にするために、ピクチャ要素を指すために、「画像画素」または「画像点」という用語が全般に本明細書で使用される。画素または画像点は、3D画像の場合はボクセルに対応する。画像処理方法は、獲得された異なる画像からの所与の試料ボクセルに対して検出される変動のマルチスペクトル定量分析を使用する。したがって、各画素において、生体組織の所与のボクセルを表す画素値の変動を表す変動画像(普通は3D画像)が、各波長に対して生成される。たとえば、血管を描写する各画素において、血流が画像ごとに大きさの変動を引き起こす。より具体的には、血管の場合、赤血球は画像ごとに厳密に同じ空間的な構造を持たず、これは時間とともに局所的な吸収を変える。標準的な光音響イメージングのように、理想的な条件(ノイズおよびスプリアスの変動がなく、ヘマトクリットが一定である)のもとでは、光音響変動画像は、局所的な光吸収μ_λ(r)とレーザーパルスのフルエンスφ_λ(r)の積に比例することを示すことができ(たとえば、Vilov S.、Godefroy G.、Arnal B.、およびBossy E.による、「Photoacoustic fluctuation imaging: theory and application to blood flow imaging」という表題の文書、Optica, 7(11)、1495-1505(2020)を参照)、

【0026】

【数1】

【0027】

ここで、rは、獲得された3次元(3D)画像における空間位置を特定するベクトルであり、Aは、超音波イメージャの点拡がり関数の空間的不変性を考慮できる場合、λまたはrに依存しない定数である。点拡がり関数が空間において変動する場合、rへのAの依存性が考慮されて補正され得る。

【0028】

複数の波長の各波長に対して、変動画像を決定することができる。この変動画像は、生体組織による(たとえば、血流による)時間的変動を表すが、レーザーのパルスごとの変動および超音波センサの電子ノイズにより影響される。

【0029】

実際には、生体組織の変動とは別に、スプリアス変動の他の原因がある。すなわち、レーザーパルスのエネルギーの変動、および特にマルチスペクトル光音響イメージングシステムの超音波センサにより生み出される電子ノイズの変動である。したがって、変動画像が、関心対象の変動、ここでは生体試料による吸収の変動のみを(またはそれを主に)表すように、特定の処理を行うことが必要である。本文書において説明される方法により取得される変動画像は、吸収変動の原因である生体試料の発色団の吸収に実際に比例する。そのような画像は「吸収変動画像」と呼ばれる。

【0030】

したがって、レーザーパルスエネルギーの変動により、変動画像は、平均の光音響画像を減算することにより取得される場合でも、平均の光音響画像に対応する残余項とパルスエネルギーの分散を乗じたものを常に含む。特異値分解(SVD)によるマルチスペクトル時空間フィルタリングは、その固有の時空間シグネチャーにより、残余項を実質的に除去する。一般に、第1の特異値に対応する画像の成分を取り除くことによって、レーザーパルスエネルギーの変動が取り除かれるだけではなく、生体試料において生じ得る(たとえば、生理学的な)動き、画像における平均の吸収要素による寄与も取り除かれ、関心対象の時間的な変動に対応する画像の部分だけが保たれる。さらに、SVDは、同じ波長を有する画像に適用されるときよりも、マルチスペクトル画像のセットに対して実行されるときに、性能が高いように見える。

【0031】

次いで、マルチスペクトル時空間フィルタリングの終わりにおいて取得される時間的に連続するフィルタリングされた画像から、フィルタリングされた変動画像が各波長に対して生成される。所与の波長におけるフィルタリングされた変動画像において、変動画像における座標rの画素は、前記波長に対して取得されるフィルタリングされた画像における同じ座標rの画素値の分布の標準偏差に等しい。

【0032】

関心対象の変動だけを保ち、関心対象の対応する変動画像または関心対象の分散画像を取得するために、各波長に対して取得されたフィルタリングされた変動画像の訂正を実行することも可能である。

【0033】

第1のタイプの訂正は、超音波センサの電子ノイズによる変動を訂正することからなる。そのような訂正により、変動画像がφ_λ(r)μ_λ(r)に比例するようになることが可能になる。ノイズは分散の空間において加法的であり、このステップは、超音波センサにより捉えられた電子ノイズから生じる画像空間における背景ノイズを推定し、分散画像の各画素からそれを減算して、訂正された分散画像を取得することからなっていてもよく、訂正された分散画像の平方根、すなわち訂正された変動画像は、各画素において、フルエンスの分布により重み付けられた、試料の対応するボクセルにおいて考慮される生体組織による(この例では血流による)変動のみを表す。それにより、訂正された変動画像は、各画像点において、関係する画像点に対応する空間点において吸収されるエネルギーに比例する。

【0034】

第2のタイプの訂正は、フルエンス関数による正規化からなり、訂正された変動画像は、超音波センサの電子ノイズによる変動を訂正した後に取得される。

【0035】

所与の波長に対してSVDによってフィルタリングされていない生の画像から取得された生のフィルタリングされていない分散画像は、各画素において、パルスのフルエンスによって重み付けられた、試料の対応するボクセルにおいて引き起こされる変動(この例では血流により引き起こされる変動)を表す分散と、画像空間における超音波センサの電子ノイズの分散との合計に対応する。数学的には、

【0036】

【数2】

【0037】

であることを示すことができ、ここで、
σ²[A_i,λ(r)]
は、iによりインデクシングされた実現値に対して計算された波長λについて獲得された生の画像(SVDによってフィルタリングされていない)から取得された生の分散画像における座標rの画素の値であり、実現値は波長λに対して獲得された時間的に連続する画像に対応し、

【0038】

【数3】

【0039】

は座標rの画素における波長λに対するレーザーの平均二乗フルエンスであり、

【0040】

【数4】

【0041】

は波長λについてのパルスエネルギーに対するレーザーパルスの相対的な変動の分散(「パルスエネルギー変動」、PEFと呼ばれる)であり、
|m_λ(r)|²
は波長λにおいて獲得された画像A_i,λ(r)の平均として計算される平均画像m_λにおける座標rの画素の絶対値の二乗であり、

【0042】

【数5】

【0043】

は、

【0044】

【数6】

【0045】

によって正規化される、関心対象の分散画像における座標rの画素の値であり、座標rの画素に対応する所与のボクセルにおける生体組織による時間的変動を表し、

【0046】

【数7】

【0047】

は画像における超音波センサによる電子ノイズの分散である。

【0048】

PEF項は2つの項に影響することがわかる。すなわち、同じものが分散画像についての係数を偏らせ、平均画像の前因子である。平均画像は生体組織による変動より1桁から2桁大きいので、PEFが小さい場合でも、関心対象の分散画像を復元するために訂正が必要である。

【0049】

フォトダイオードが測定できないレーザービームの空間的変動があるので、レーザーの出力におけるフォトダイオードによるPEFの監視も、不十分な訂正方法であり得る。SVD方法は、完全なフィルタリングされた画像を取得するために、PEFの効果を除去する。SVDによるマルチスペクトル時空間フィルタリングの間の第1の特異値の除去には、PEFを受ける平均の吸収物の寄与を取り除く効果がある。SVDフィルタリングは、分散画像の背景レベルも変える。したがって、SVDによるマルチスペクトル時空間フィルタリングの後、フィルタリングされた分散画像が取得され、これは

【0050】

【数8】

【0051】

などのフィルタリングされた変動画像に対応し、

【0052】

【数9】

【0053】

はSVDフィルタリングの後の残余電子ノイズの分散であり、

【0054】

【数10】

【0055】

は波長λにおけるレーザーパルスエネルギーの相対的な変動の分散(平均値の二乗により正規化される分散)である。

【0056】

レーザーパルスの相対的な変動は、普通は数パーセントであり、無視することができる。実際に
σ_φ,λ～10^-2
であり、

【0057】

【数11】

【0058】

である。

【0059】

よって、

【0060】

【数12】

【0061】

と書くことができる。したがって、

【0062】

【数13】

【0063】

である。

【0064】

それにより、生体組織の吸収変動を表す吸収変動画像σ_λ,flow(r)は、2回の訂正を連続して適用することによって取得され得ることがわかる。

【0065】

第1の訂正は、フィルタリングされた分散画像の各画素において、訂正された分散画像、そして平方根をとることにより対応する訂正された変動画像を取得するために、SVDによるフィルタリングの後に残余電子ノイズの分散を減算することからなる。

【0066】

第2の訂正は、吸収変動画像σ_λ,flow(r)または関心対象の変動画像と呼ばれる、訂正され正規化された変動画像を取得するためにレーザーパルスのフルエンスφ_λ(r)で各画素を割ることによって、φ_λ(r)により訂正された変動画像を正規化することからなる。考慮される各波長に対して、2回の訂正が適用される。

【0067】

それにより各励振波長に対して取得される吸収変動画像σ_λ,flow(r)は、生体組織により引き起こされる吸収変動のみを(または少なくとも主にそれを)表す。生体組織における吸収変動が主に血流による吸収変動である場合、σ_λ,flow(r)は血流による吸収変動を表す。

【0068】

よって、イメージングによる変動という手法が、導入部分において言及された可視性の問題を解決するために、かつコントラストを改善するために、一連の画像を獲得するという代償を払って使用され得る。たとえば、大きい血管内での変動についての情報および垂直方向の血管についての情報が取得され得る。

【0069】

そのような方法は、獲得の間に動かす必要なく、所与の領域をイメージングするために被験者の皮膚に置かれる、手で持たれるまたは機械部品により保持される固定獲得デバイス(「ハンドヘルドデバイス」または「片面デバイス」と呼ばれる)による、画像獲得の場合において特に有用である。そのような獲得デバイスでは、音響波センサは、たとえば可視性が限られている(通過帯域および角度受信スペクトルが限られている)圧電センサである。

【0070】

さらに、血液の酸素化率の定量分析が可能である。実際に、異なる波長に対して取得される異なる吸収変動画像からの、血管の酸素化のイメージングに適用する場合、2Dまたは3D酸素化マップが、フルエンスにより正規化された異なる波長における血流による変動マップを酸素化率と結びつけるモデルを反転させることによって決定され得る。血管の光吸収スペクトルの定量化は、それから達成され得る。より正確には、複数の波長における組織励振光の波を使用して、オキシヘモグロビンおよびデオキシヘモグロビンの相対的または絶対的な量が、波長に従って取得された吸収プロファイルの差から決定され得る。

【0071】

酸素化マップを取得するために、組織によるスペクトルカラーレーション(spectral coloration)がない場合に入射フルエンスにより訂正され正規化される変動画像、または、試料内の空間フルエンスにより訂正され正規化される変動画像である吸収変動画像のいずれかから開始することが可能である。第1の場合には、各点において、吸収変動画像との比例関係があるが、深さ方向の光の減衰は訂正されない(これは、組織によるスペクトルカラーレーションがないとき、各波長に対して同じであると仮定される)。

【0072】

それにより、血管の酸素化率の定量的なイメージングが取得される。吸収変動画像は血流に固有であり、他の吸収物により囲まれている低解像(under-resolved)血管が関わるときに特に、定量的評価において固有性を提供し、スペクトル分解の複雑さを回避する。この方法は、腫瘍のイメージング、脳のイメージング、および血管のイメージングに特に適用可能である。上で言及されたように、変動画像は、従来の画像よりコントラストが高く、従来のイメージングに存在する可視性アーティファクトにより影響されない。したがって、従来の画像よりもはるかに豊かな画像について、酸素化率が取得される。

【0073】

図1は、マルチスペクトル光音響イメージングシステムによって獲得される画像を処理するための方法の簡略化されたフローチャートを示す。

【0074】

ステップ110の間に、M_λ個の波長に対する光音響イメージングシステムによって時間的に連続する画像が獲得され、すなわち、波長当たりN個の画像、全体でN^*M_λ個の画像が獲得される。獲得された画像はA_i,jと表記され、iは1とNの間に含まれる整数であり、これは所与の波長に対する獲得された画像のインデックスを表し、jは1とM_λとの間に含まれる整数であり、関係する波長を特定するインデックスを表す。画像A_i,jの画素はA_i,j(r)と表記され、rは画像の画素の座標である。

【0075】

獲得率は一定であり、たとえば100Hzである。通常は、M_λ=6個の波長および波長当たりN=250個の画像が使用される。より一般的には、Nは10と100,000の間に含まれてもよく、M_λは2と100の間に含まれてもよく、獲得周波数は10Hzと100,000Hzの間に含まれてもよい。

【0076】

獲得は、たとえば、最後の波長λ_Mまで、第1の波長λ₁を用いた画像、次いで第2の波長λ₂を用いた画像などのように、画像ごとに波長を周期的に変え、そして、波長λ₁からλ_Mについて同じ獲得サイクルをN回繰り返すことによって行われ得る。このようにすると、同じ波長における2回の獲得と獲得の間の期間は、波長にかかわらず一定であり同一である。2つの波長のみが使用される場合、第1の波長λ₁は第2の波長λ₂と交互にされる。

【0077】

獲得された画像A_i,jは、以下で定義されるステップ120から150を含む後処理を受ける。

【0078】

ステップ120の間に、特異値分解(SVD)によるマルチスペクトル時空間フィルタリングがすべてのN^*M_λ個の画像に適用され、N^*M_λ個のフィルタリングされた画像がステップの終わりに取得され、

【0079】

【数14】

【0080】

と表記される。

【0081】

SVDによるそのようなマルチスペクトル時空間フィルタリングは、分散画像において、平均の吸収物の成分、特に、その大きさによって関心対象の変動を覆い隠すレーザーのパルスエネルギーの変動および平均画像から生じる項

【0082】

【数15】

【0083】

を実質的に取り除くために使用される。複数の波長に対してSVDフィルタリングを使用することで、SVDが統計的により効率的になり、複数の励振波長に応答した組織の変動を表す成分が得られる。SVDによる時空間フィルタリングが単一の波長に対して獲得された画像について使用されるとすると、SVDフィルタリング範囲の選択は波長ごとに変わり、これは手順を面倒なものにする。

【0084】

ステップ140の間に、各波長に対して別々に、フィルタリングされた変動画像(それぞれ、フィルタリングされた分散画像)が決定される。前記波長に対してステップ120の間に取得されたN個のフィルタリングされた画像における同じ座標rの画素値の分布の標準偏差(それぞれ、分散)が、フィルタリングされた変動画像(それぞれ、フィルタリングされた分散画像)における座標rの画素を取得するために計算される。それにより、各波長jに対して、ステップ140の終わりに、M_λ個の分散画像が取得され、これは

【0085】

【数16】

【0086】

と表記され、またはより簡単に、

【0087】

【数17】

【0088】

と表記される。上記の表記

【0089】

【数18】

【0090】

について、

【0091】

【数19】

【0092】

は、波長jについて実現値i=1からNに対して計算される座標rの画素における値の分布の分散

【0093】

【数20】

【0094】

である。

【0095】

ステップ150の間に、生体試料の吸収変動だけを、または基本的にそれを表す、吸収変動画像を取得するために、各々のフィルタリングされた変動画像の訂正が実行される。それにより、フィルタリングされた変動画像の訂正は、試料による吸収変動を表す吸収変動画像を取得するために、試料による関心対象の変動以外の変動を取り除くことを含む。2つのタイプの訂正が、組み合わせて、または個別に使用され得る。

【0096】

第1のタイプの訂正(ステップ150A)は、フィルタリングされた分散画像における座標rの各画素において、SVDの後の残余電子ノイズの分散を減算することによる、超音波センサの電子ノイズによる変動の訂正からなり、SVDの後の残余電子ノイズの分散は、

【0097】

【数21】

【0098】

により表記され、すなわち

【0099】

【数22】

【0100】

である。

【0101】

獲得センサの残余電子ノイズの分散が推定され、次いで、残余電子ノイズの分散

【0102】

【数23】

【0103】

が、所与の波長に対する訂正された分散画像を取得するために、ステップ140の間に取得されたフィルタリングされた分散画像の各画素において減算される。SVDの後の残余電子ノイズの分散は、SVDフィルタリングの限界bが画像の総数に等しい場合、SVDの前の電子ノイズの分散に対応し得る。

【0104】

対応する訂正された変動画像は、訂正された分散画像の平方根を計算することによって取得され得る。

【0105】

次いで、訂正された変動画像は、考慮される波長のフルエンスによって重み付けられ得る。訂正された変動画像は、各画像点において、関係する画像点に対応する空間点において吸収されたエネルギーに比例する。

【0106】

超音波センサによって生み出される電子ノイズは、試料の不在時に獲得された光音響信号において生み出される電子ノイズの分散の関数として推定されてもよく、この分散は、特異値への分解によるマルチスペクトル時空間フィルタリングにより除去されるノイズの量によって訂正され、ノイズの量は、特異値への分解によるマルチスペクトル時空間フィルタリングにより取り除かれる最低エネルギー成分に対応する特異値に基づいて推定され得る。

【0107】

第2のタイプの訂正(ステップ150B)は、考慮される波長におけるレーザーパルスのフルエンスの関数による正規化からなる。フルエンスによって正規化される分散画像は、二乗されたフルエンス空間関数[Math.37](数24)

【0108】

【数24】

【0109】

による正規化によって各波長に対して取得され得る。

【0110】

各波長j=1からM_λに対して、試料によって(たとえば、血流によって)生み出される吸収変動だけを表す正規化された分散画像

【0111】

【数25】

【0112】

を取得するために、

【0113】

【数26】

【0114】

は、座標rの画素において波長λに対して取得される正規化された分散画像の画素値を表す。

【0115】

フルエンスによる正規化は、各画素の平方根を計算することによって、正規化された分散画像を、そして対応する正規化された変動画像を取得するために、二乗フルエンス関数φ² _λ(r)を使用して、訂正された分散画像に対して実行され得る。代替的に、および等価的に、フルエンスによる正規化は、レーザーパルスのフルエンスφ_λ(r)による正規化の前に、電子ノイズ訂正の後に取得された訂正された分散画像の各画素の平方根を計算することによって、訂正された分散画像に対応する変動画像に対して実行され得る。

【0116】

【数27】

【0117】

そのような正規化に続いて、正規化された変動画像が最終的に取得され、これは、試料による吸収変動を表す吸収変動画像σ_λ,flow(r)である。

【0118】

正規化のために使用されるべきフルエンス関数の決定に関して、異なる方法および近似が可能である。

【0119】

フルエンスの空間関数を使用する代わりに、所与の波長における(空間点に無関係な)平均フルエンスが使用されてもよい。ある波長における平均フルエンスは、たとえば、所与の波長の励振パルスが試料に送られる瞬間にレーザーの出力に置かれたフォトダイオードによって取得される励振パルスのエネルギー測定結果に基づいて決定され得る。フォトダイオードによって生み出される電気信号の強さは、較正係数を使用してフルエンスの推定値へと変換され、フルエンスの推定値は、平均のフルエンスを計算し、そして正規化による訂正を実行するために使用される。

【0120】

フルエンス空間関数は、
φ_λ(r)=φ_λ(0)ξ_λ(r)
と定義され得る。

【0121】

ここで、
φ_λ(0)
は試料の表面において測定されるフルエンスであり、これは較正されたフォトダイオードによって取得することができ、
ξ_λ(r)
はフルエンスの相対的な空間的変動を表す。

【0122】

一般的な場合、フルエンスの空間的分布の推定は、モデリングまたは測定によって取得され得る。

【0123】

いくつかの場合、組織によるスペクトルカラーレーションはないと近似することができる。そのような場合、
ξ_λ(r)=ξ(r)
と書くことができる。

【0124】

そのような場合、酸素化の推定は単に、フォトダイオードのフルエンスによる正規化から行われてもよく、それは、同じ定数
ξ(r)
をすべての波長に使用できるからである。

【0125】

フルエンスが空間的に一様であるとき(ニワトリ杯を仮定)、フルエンスは、フォトダイオードにより生み出される信号から直接推定することもでき、すなわち
φ_λ(r)=φ_λ(0)
である。

【0126】

ステップ160の間に、酸素飽和度の画像α(r)が、少なくとも2つの波長λについて、ステップ150の終わりに取得された少なくとも2つの吸収変動画像σ_λ,flowから決定される。代替として、ステップ160は、ステップ140の終わりに取得された対応する分散画像の各画素の平方根を計算することによって取得される、少なくとも2つの訂正された変動画像σ_λ(ステップ150Aによる訂正を伴うがステップ150Bによる訂正を伴わない)からも実行され得る。

【0127】

ステップ160が、ステップ150の終わりに取得された吸収変動画像σ_λ,flowに基づいて実行される(同じ計算が、ステップ150Bによる訂正を経ていない訂正された変動画像以降で使用可能である)、例示的な事例が説明されている。この目的で、座標rの各画素に対して、吸収変動画像の画素rにおける値
σ_λ,flow(r)
と、総ヘモグロビン濃度および電子受信機システムの感度係数を含むK(r)、ならびに酸素飽和度α(r)という2つのパラメータとの関係を表す、モデルが使用される。

【0128】

このモデルは、たとえば、2つの未知数K(r)およびα(r)を伴う式からなり、変動
σ_λ,flow(r)
はK(r)およびα(r)の関数として与えられ、rは画像画素の座標である。したがって、座標rの各画素に対して、酸素飽和レベルα(r)を決定し、したがって、各画素に対して、酸素飽和についての定量的情報(0%と100%の間)を備える画像αを生成することが可能である。同様に、各画素に対して、「血液量」とも呼ばれる総ヘモグロビン濃度K(r)に比例する定量的情報を備える画像Kを生成することが可能である。

【0129】

一実施形態では、このモデルは、

【0130】

【数28】

【0131】

という式に基づく分析モデルであり、

【0132】

【数29】

【0133】

は波長λにおけるオキシヘモグロビンの吸収係数であり、
μ_λ,Hbb
は波長λにおけるデソキシヘモグロビンの吸収係数である。

【0134】

K(r)は、位置rに依存するが波長λには依存しない前因子である。

【0135】

図2は、各ステップにおいて取得される画像の例に基づくステップ120から160の態様を示す。

【0136】

この例では、入力として取得される、獲得された画像A_ijは、3次元(3D)画像である。したがって、
A_1,1
から

【0137】

【数30】

【0138】

により表記される、N^*M_λ個の獲得された3D画像を得る。

【0139】

ステップ120(SVD)の終わりにおいて、

【0140】

【数31】

【0141】

から

【0142】

【数32】

【0143】

により表記される、N^*M_λ個のフィルタリングされた3D画像を得る。

【0144】

ステップ140の終わり(分散画像の計算)において、波長j=λ₁からλ_Mに対して、

【0145】

【数33】

【0146】

から

【0147】

【数34】

【0148】

により表記される、M_λ個の3D変動画像を得る。ステップ150の終わり(ノイズ訂正および/またはフルエンスによる正規化)において、波長j=λ₁からλ_Mに対して、

【0149】

【数35】

【0150】

から

【0151】

【数36】

【0152】

により表記される、吸収変動を表すM_λ個の3D分散画像を得る。最後に、ステップ160の終わり(酸素飽和の計算)において、酸素飽和レベルの3D画像が取得され、α(r)と表記される。

【0153】

図2は、例示的な事例について、異なる3D画像に対する異なる処理の効果を観察することを可能にする。ステップ120に続くSVDの後にフィルタリングされた対応する画像と獲得された画像を比較することによって、血流の変動に対応する画像の点だけが保たれていることが観察される。ステップ120に続くSVDの後のフィルタリングされた画像を、ステップ140に続いて取得されたフィルタリングされた変動画像と比較することによって、いくつかの構造がフィルタリングされた変動画像に現れることが観察される。ステップ140に続いて取得されたフィルタリングされた変動画像を、ステップ150に続いて取得された吸収変動画像と比較することによって、背景レベルが0に下がることが観察される。最終的に、ステップ160の後に、酸素飽和レベルのボリューム画像が取得されることが観察される。

【0154】

図3Aから図3Bは、模範的な実施形態による、ステップ120のために使用され得るSVDによるマルチスペクトル時空間フィルタリング方法を示す。

【0155】

図3Aに示されるものに従って、時間的に連続するマルチスペクトル画像(1000)が、M_λ個の波長のために光音響イメージングシステムによって獲得され、すなわち、波長当たりN個の、全体でN^*M_λ個の画像が獲得される。入力として取得される獲得された画像A_i,jは、同一のサイズn_X ^*n_Y ^*n_Zを有すると仮定される3D画像である。画素の数は、3D直交空間の第1のX軸に沿ってn_X個であり、第2のY軸に沿ってn_Y個であり、第3のZ軸に沿ってn_Z個である。

【0156】

A(r,t)と表記され、空間および時間を表す、時間的に連続するマルチスペクトル画像から、2次元(2D)Casorati行列(301)がステップ301の間に形成され、tは時間インデックスであり、rは3D空間における画素の位置に対応するスカラーであり、N^*M_λ個の画像が定義される。Casorati行列(301)のサイズはn_R ^*(N^*M_λ)であり、n_R=n_X ^*n_Y ^*n_Zである。Casorati行列(301)は、N^*M_λ個の画像の画素を含む。それにより、時間インデックスの値tは、t=1からN^*M_λまで変わることができるので、A(r,1)=A_1,1(r)は第1の波長に対する第1の時間的に連続する画像であり、A(r,2)=A_1,2(r)は第2の波長に対する第2の時間的に連続する画像であり、以下同様である。表記A(r,1)ではrはスカラーであるが、表記A_i,j(r)ではrは3次元空間における座標(x,y,z)のベクトルであることに留意されたい。スカラーrは、ベクトルの座標r=(x,y,z)および獲得された画像のサイズn_X ^*n_Y ^*n_Zから計算され得る。Casorati行列(301)の使用は、時空間領域における変動の分析を実行することを可能にする。

【0157】

ステップ320の間に、Casorati行列(301)は、特異値分解方法に従って分解される。A=USV^*であり、Uは空間特異ベクトルを備えるサイズn_R ^*n_Rの2D行列であり、Sはサイズn_R ^*(N^*M_λ)の特異値の2D行列であり、その対角係数(diagonal coefficient)S_kは正の実数であり、V^*はVの共役転置行列であり、Vはサイズ(N^*M_λ)^*(N^*M_λ)の空間特異ベクトルの行列である。行列S(302)の例は図3Aに示される。一般的な慣習は、値の降順で値S_kを並べることである。

【0158】

図3Bに示されるものに従って、特異値分解に基づいて、行列A(r,t)は、その加重係数が特異値の2D行列の対角係数S_kである、行列の加重和として表現され得る。

【0159】

【数37】

【0160】

上記の合計において、各行列U_k(r)V_k(t)^*は、特異値S_kに対応する成分を定義する。フィルタリングを実行するために、一部の成分のみが保たれるので、加重和の一部の行列のみが保たれる。値S_kが値の降順で並べられるとき、それは、ステップ330の間に、マルチスペクトル時空間フィルタリングの後に取得された行列が

【0161】

【数38】

【0162】

となるように、インデックスkの最小インデックスaおよび最大インデックスbの値を選択するという問題である。

【0163】

最小インデックスaの値の選択は、そのインデックスが、どの高エネルギーの成分がフィルタリングにより除去されるかを決定することによってフィルタリングの効率と妥当性を決定するという点で、重要である。最大インデックスbの値の選択は、フィルタリングに対する影響がより小さく、それは、その選択が、純粋なノイズおよび場合によってはそれに埋め込まれている情報を含む低エネルギーの成分にも関係し、たとえばa+100に等しくなるように選ばれ、またはa+1とN^*M_λの間に含まれ得るからである。

【0164】

最高エネルギーの特異値に対応する成分を除去した後、ステップ340の間に変換が実行され、これは、

【0165】

【数39】

【0166】

により表記される時間的に連続するフィルタリングされた画像を取得するための、行列A_SVD(r,t)からステップ310の間に実行された変換の逆である。

【0167】

図4Aから図4Bは、ステップ120の間に実行されるSVDによるマルチスペクトル時空間フィルタリングの間に取り除かれるべき成分を選択するための方法の態様を示す。

【0168】

図4Aは、ステップ120の間に実行されるSVDによるマルチスペクトル時空間フィルタリングの間に除去されるべき成分を選択するための方法の簡略化されたフローチャートを示し、特に、最小インデックスa(下限とも呼ばれる)の値の選択は、どの成分がフィルタリングにより除去されるかを特定する。

【0169】

最小インデックスaの選択は、最小インデックスaの値を大きくしすぎることによって、画像において表現される物体の一部を取り除く危険性があるという点で、繊細であり得る。本明細書で提案される方法は、画像において表現される物体の一部を削除する危険性を伴わずに、最小インデックス値aを自動的に選択することを可能にする。

【0170】

この方法は、インデックスaの複数の値に対してSVDによって(ステップ150Aのノイズ訂正も、ステップ150Bのフルエンスによる正規化も伴わずに)フィルタリングされる画像について推定されるコントラスト対ノイズ比(CNR)の推定に基づき、コントラスト対ノイズ比を最大にする最小インデックスaの値が決定される。

【0171】

SVDによってフィルタリングされる画像と、マスキングされた画像との比較に基づいて、それらに現れる構造を抽出するために、コントラスト対ノイズ比の評価が実行される。マスクは、ステップ410および420に基づいて取得される。

【0172】

ステップ410の間に、最小インデックスaの基準値a_rが任意で選択される。たとえば、a_r=0.03×N×M_λ=0.03×10×100=30である。そのような値は、少数の代表的なサンプルについての分析から経験的に選択され得る。

【0173】

次いで、単一の波長における基準分散画像が選択され、これは、最小インデックスaのこの基準値a_rに対して、かつ、たとえば生の信号のSNR(「信号対雑音比」)が最高である波長である、所与の波長λ=λ_jmaxに対して、

【0174】

【数40】

【0175】

と表記される。

【0176】

【数41】

【0177】

の計算は、インデックスj_maxの波長に対して獲得される画像のためにステップ120(SVD)および140(フィルタリングされた分散画像)を実行することを含むが、訂正ステップ150の実行は含まない。

【0178】

ステップ420の間に、バイナリマスクとして使用される2つのバイナリ画像が計算される。

【0179】

M_σは、閾値th_σに基づいて、参照分散画像

【0180】

【数42】

【0181】

に基づいて計算されるバイナリマスクである。

【0182】

M_mは、平均画像
<A>_jmax
に基づいて計算されるバイナリマスクであり、平均画像の各画素は、iについての、したがって、たとえばλ=λ_jmaxとして選ばれる少なくとも1つの所与の波長に対して獲得されたすべての画像についての、画素の平均
A_i,jmax(r)
として計算される。別の閾値th_mが使用される。

【0183】

バイナリマスクM_σを取得するための閾値th_σの決定は、

【0184】

【数43】

【0185】

という式に基づいて実行されてもよく、ここで、

【0186】

【数44】

【0187】

は、たとえばステップ150Aについて上で説明されたように決定される、電子ノイズの時空間分散である。

【0188】

【数45】

【0189】

は、基準分散画像の画素分布について計算される画像内標準偏差である。

【0190】

θは、少数の例について、取得されたマスクM_σが画像

【0191】

【数46】

【0192】

上で観察可能な構造に対応することを確実にすることによって経験的に取得される、たとえば0.35に等しい加重係数である。

【0193】

バイナリマスクM_mを取得するための閾値th_mの決定は、

【0194】

【数47】

【0195】

という式に基づいて行われてもよく、ここで、

【0196】

【数48】

【0197】

は、

【0198】

【数49】

【0199】

と表記される平均画像における画素の平均値であり、

【0200】

【数50】

【0201】

は、平均画像

【0202】

【数51】

【0203】

の画素分布について計算された画像内標準偏差である。

【0204】

係数θは、マスクM_σを取得するために使用される係数と同じである。

【0205】

バイナリマスクは、画像において、画像の背景を、残りの部分、特に画像において表される1つまたは複数の物体と区別するために使用される。既知の方法で、画像および閾値に基づいてバイナリマスクを計算するために、画像の画素値が閾値より大きい(またはそれに等しい)かどうかが決定される。大きい場合、バイナリマスクにおける対応する画素の値は1に等しく、大きくない場合、バイナリマスクにおける対応する画素の値は0に等しい。

【0206】

取得される2つのバイナリマスクがCNRの式において使用され、その計算は、SVDによるマルチスペクトル時空間フィルタリングのために使用されるべき最小インデックスaを決定することを可能にする。このフィルタリングはすべての波長に対して実行され、最小インデックスaはすべての波長に対して同じである。

【0207】

ステップ430の間に、最小インデックスaの複数の値が、たとえば1および100などのあらかじめ定められた値のセットから選択され、分散画像はすべての波長において計算され、最小インデックスaの各値に対して

【0208】

【数52】

【0209】

と表記される。

【0210】

ステップ440の間に、コントラスト対ノイズ比の値CNR(a)は、2つのバイナリマスクを使用することによって、最小インデックスaの各値に対して計算される。

【0211】

ステップ450の間に、コントラスト対ノイズ比CNR(a)が最大である最小インデックスaの値は、下限として選択される。処理の残り(特に、図1に関して説明された方法のステップ150および160)を行うための下限に対応する、以下の分散画像が使用される。

【0212】

【数53】

【0213】

コントラスト対ノイズ比(CNR)は、画像のどの画素が画像の背景の部分であるか、および画像のどの画素がその画像において表される1つまたは複数の物体の部分であるかを区別することを可能にするバイナリマスクに基づいて、様々な方法で計算され得る。一般に、コントラスト対ノイズ比は、1つまたは複数の物体の部分を形成する画素の平均値と、画像の背景の部分を形成する画素の平均値との差として決定されるコントラストの、画像の背景の変動に対する比である。基準としてのコントラスト対ノイズ比の選択は特に、画像において病理を視覚的に検出するために、そのコントラストが画像の背景の変動より大きくなければならず、したがってコントラスト対ノイズ比が1より大きくなければならないという事実に起因する。

【0214】

模範的な実施形態によれば、コントラスト対ノイズ比は、最小インデックスaの値に対して取得される各波長における各分散画像

【0215】

【数54】

【0216】

に対して計算され、波長にわたって次のように平均化され、

【0217】

【数55】

【0218】

ここで、Mは波長の数であり、
∩
は2つのバイナリマスクの共通部分を表す演算子であり、論理和関数に対応する共通部分は、2つのマスクの間で画素ごとに実行され、

【0219】

【数56】

【0220】

は画素の分布について計算される画像内標準偏差演算であり、

【0221】

【数57】

【0222】

はマスクM_σの補集合、すなわち1-M_σであり、

【0223】

【数58】

【0224】

は画像の背景に属す画素に対応する項であり、
<>_r
は画素分布について計算される画像内平均化演算である。

【0225】

そのようなタイプの計算式は、平均画像と分散画像の両方に存在する構造を、コントラスト対ノイズ比の計算において除去するために、項
M_σ∩M_m
によって使用される。光音響イメージング技法によって獲得される画像において、画像の平均値は残りよりも約100倍大きく、これは、平均画像値によるコントラストが最小インデックスaの低い値について減算されなかった場合、非常に高いCNR(a)の値をもたらす。それにより、使用されるマスクは、新しい構造が現れるときにコントラストが上昇することを確実にする。それにより、最小インデックスaが何であってもCNR(a)の値を妥当に比較できるように、前記特定の計算式が使用される。

【0226】

図4Bに示されるCNR(a)変動の例示的な曲線(470)は、コントラスト対ノイズ比CNR(a)が、25から28の間に含まれる最小インデックスaの値に対して最大値を有することを示す。別の例では、最小インデックスaの値は33と35の間に含まれる。

【0227】

ステップ150Aの間に使用されるセンサの電子ノイズの推定のために、複数の方法が使用され得る。

【0228】

再構築された画像

【0229】

【数59】

【0230】

の場合に適用可能な、第1のノイズ推定方法によれば、再構築された画像における電子ノイズは、実際の無線周波数信号の空間ノイズ測定結果から差し引かれ得る。

【0231】

【数60】

【0232】

ここで、
N_trans
は再構築において使用されるトランスデューサの数である。

【0233】

【数61】

【0234】

は、測定された無線周波数信号についての電子ノイズの分散であり、試料の不在時に取得された生の信号から決定される。

【0235】

ここで、係数2は、複素数値を持つ信号が再構築のために使用されるという事実に起因する。

【0236】

しかしながら、

【0237】

【数62】

【0238】

の推定は、SVDによるマルチスペクトル時空間フィルタリングが

【0239】

【数63】

【0240】

に等しい量のノイズを取り除くことを考慮することによって訂正されてもよく、ここで、bは
A_SVD(r,t)
を定義する行列の加重和の最大インデックスである。

【0241】

N_XN_YN_Zは再構築された画像における画素/ボクセルの数であり、S_kはSVDによって取得される特異値である。

【0242】

【数64】

【0243】

という表現に存在する、ノイズを表す特異値は、波長に関して不変であることが観察されている。

【0244】

最後に、訂正された分散画像を取得するために、分散画像からステップ150Aの間に減算されるべきSVDの後の残余電子ノイズは、

【0245】

【数65】

【0246】

である。

【0247】

第2のノイズ推定方法によれば、ノイズが前記法則に従わない場合、たとえば、定数の減算を受ける再構築された画像

【0248】

【数66】

【0249】

のノルムL2を最小にすることによって、代替の訂正を実行することが可能である。

【0250】

【数67】

【0251】

ノルムを最小にする定数
Q²
の値は、

【0252】

【数69】

【0253】

から減算されるべき値

【0254】

【数68】

【0255】

を計算するために使用され得る。

【0256】

ノイズのレベルを推定するための他の方法が使用され得る。

【0257】

本文書において説明される光音響画像処理方法は、ニワトリ胚において獲得されたボリューム画像に適用され、可視の構造における平均画像に基づく従来の光音響分光法を用いた酸素化率に近い値を得た。

【0258】

図5Aから図5Bは、本文書で説明される方法を血管の酸素化のイメージングに適用した場合の、その方法による生体組織のいくつかの要素の可視性の改善を示す。2つの例が示されている。両方の例(図5Aおよび図5B)において、第1の線は、左から右に、3つの投影面YZ(a)、XY(b)、およびXZ(c)に沿って平均光音響画像m_λから取得された酸素化画像に対応する。酸素化画像は可視性が限られている。吸収変動イメージングから取得された酸素化3D画像が第2の線上に示されており、投影は同じ平面YZ(d)、XY(e)、およびXZ(f)に従い、従来の光音響イメージング技法を用いた場合よりも多くの構造を見ることができる。

【0259】

図5Aにおいて、イメージングされるゾーンは絨毛尿膜に対応する。吸収変動イメージングの技法(図5A、画像d、e、f)は、酸素化の測定をより多くの血管に拡張することを可能にし、それは、従来の画像(図5A、画像a、b、c)では、向きまたはサイズにより一部の血管が現れなかったからであった。

【0260】

図5Bにおいて、イメージングされるエリアはニワトリ胚の心臓に対応する。変動技法(図5B、画像d、e、f)は、酸素化の測定を器官全体に、ならびに外に向かう血管に拡張することを可能にするが、従来の画像(図5B、画像a、b、c)には、不連続な画素分布としてわずかな情報しか現れない。

【0261】

媒体が透明である(澄んだ媒体により囲まれる血管新生を伴うニワトリ胚の場合のように)場合、生体組織による複数の波長の使用が原因のスペクトルカラーレーションを訂正することは必要ではない。

【0262】

本明細書で説明された血流酸素化イメージング方法は、わずかなハードウェアの改変で従来の検出器の限界を超える。関心対象の変動を抽出して本明細書で説明されるような後処理を実行するには、十分な画像(通常は少なくとも50個)を獲得することで十分である。

【0263】

本文書において提示される機能、モーター、原理図、流れ図、状態遷移図、および/またはフローチャートは、本発明の原理を実施する説明のための回路の概念図を表す。同様に、すべてのフローチャート、流れ図、状態遷移図、疑似コード、および他のものは、方法の様々な態様を表し、方法は、コンピュータ可読媒体に記憶されているコンピュータプログラム命令によって基本的に実施され、それにより、プロセッサまたはデバイスが明確に表現されているかどうかにかかわらず、プロセッサまたはプロセッサを含むデバイスによって実施され得る。

【0264】

一実施形態によれば、光音響画像処理方法のすべてのステップの1つまたは複数が、ソフトウェアまたはコンピュータプログラムによって実施される。

【0265】

それにより、この説明は、データプロセッサによって実行され得るコンピュータプログラムに関係し、コンピュータプログラムは、本文書において説明される実施形態のいずれかに従った光音響画像処理方法のステップの1つまたは複数またはすべてのデバイスによる実行を制御するためのコンピュータプログラム命令を含む。コンピュータプログラム命令は、たとえば、デバイスのメモリに記憶され、ロードされ、そしてこのデバイスのプロセッサによって実行されることが意図される。

【0266】

コンピュータプログラムは、任意のプログラム言語を使用することができ、ソースコード、オブジェクトコード、または、部分的にコンパイルされた形式もしくは任意の他の望ましい形式などの、ソースコードとオブジェクトコードとの間の中間コードの形式であり得る。

【0267】

デバイスは、1つまたは複数の物理的に別個の機械によって実装されてもよく、全体として、コンピュータのアーキテクチャを有し、これは、データメモリ、プロセッサ、通信バス、このコンピューティングデバイスをネットワークまたは別の機器に接続するためのハードウェアインターフェース、ユーザインターフェースなど、そのようなアーキテクチャのコンポーネントを含む。

【0268】

本説明はさらに、データプロセッサによって読み取り可能な、かつ上で言及されたようなコンピュータプログラムの命令を含む、情報媒体に関する。

【0269】

情報媒体は、プログラムを記憶するのに適した任意のエンティティまたはデバイスであり得る。

【0270】

情報媒体は、信号を記憶するのに適した任意のハードウェア手段、エンティティ、またはデバイスであり得る。たとえば、媒体は、ROMまたはRAMメモリ、たとえばCD ROMディスクもしくは磁気記録手段、コンピュータハードディスク、光記憶媒体、フラッシュメモリデバイス、および/または情報を記憶するための他の有形な機械可読媒体などの、記憶手段を含み得る。

【0271】

「コンピュータ可読媒体」という用語は、限定はされないが、ポータブルまたは固定記憶デバイス、光記憶デバイス、命令および/もしくはデータを記憶し、格納し、または輸送するのに適した様々な他の媒体を含み得る。

【0272】

それは、コンピュータ記憶媒体および/もしくは通信媒体であってもよく、または、より一般的には、ある場所から別の場所へのコンピュータプログラムの移送を容易にする任意の媒体であってもよい。

【0273】

コンピュータ可読媒体の例は、限定はされないが、フラッシュドライブまたは他のフラッシュメモリデバイス(たとえば、メモリキー、メモリスティック、USBキードライブ)、CD-ROMまたは他の光学ストレージ、DVD、磁気ディスクストレージまたは他の磁気記憶デバイス、ソリッドステートメモリ、メモリチップ、RAM、ROM、EEPROM、スマートカード、リレーショナルデータベース管理システム、企業データ管理システムなどを含む。

【0274】

情報媒体は、電磁信号(電気、無線、または光信号)などの搬送波の形態で送信可能な媒体であってもよく、これは、電気もしくは光ケーブル、無線もしくは赤外線リンク、または他の手段のような、適切な有線または非有線送信手段を介して搬送され得る。

【0275】

「プロセッサ」という用語は、たとえば、1つまたは複数の処理ユニットまたはハードウェアに基づく1つまたは複数の処理カーネルを備える、任意のマイクロプロセッサ、マイクロコントローラ、コントローラ、集積回路、または中央処理装置(CPU)を指し得る。その上、「プロセッサ」という用語は、コンピュータプログラム命令を実行するのに適したハードウェアだけを指すものとして解釈されるべきではなく、たとえば、デジタルシグナルプロセッサ(DSP)、ネットワークプロセッサ、特定用途向け集積回路(ASIC)、フィールドプログラマブルゲートアレイ(FPGA)、または、プログラム可能かどうかにかかわらず、特定用途向けであるかどうかにかかわらず、他の回路を指し得る。「プロセッサ」という用語は、本明細書で言及された複数の模範的な実施形態の組合せにも対応し得る。

【0276】

本発明は、プログラムコード命令を備える少なくとも1つのデータメモリと、少なくとも1つのデータプロセッサとを備える、光音響画像処理デバイスに関し、データプロセッサは、プログラムコード命令がデータプロセッサによって実行されると、光音響画像処理デバイスが本明細書で説明される実施形態のいずれかによる光音響画像処理方法のステップの1つまたは複数またはすべてを実行するように、構成される。

【0277】

一実施形態では、光音響画像処理デバイスは、
本文書で説明された実施形態のいずれかによる光音響画像処理方法のステップの1つまたは複数またはすべての実行を制御するように設計されるコンピュータプログラム命令を記憶するための、データ記憶手段、たとえば1つまたは複数のメモリと、
本文書で説明された実施形態のいずれかによる光音響画像処理方法のステップの1つまたは複数またはすべてを実施するためにコンピュータプログラム命令を実行するように構成される、データ処理手段、特にデータプロセッサとを備える。

【0278】

より一般的には、光音響画像処理デバイスは、本文書で説明される実施形態のいずれかによる光音響画像処理方法のステップの1つまたは複数またはすべてを実施するための手段を備える。この手段は、たとえばソフトウェアおよび/またはハードウェア手段を備える。

【符号の説明】

【0279】

470 曲線

【図1】