特許第6924530号(P6924530)IP Force 特許公報掲載プロジェクト 2022.1.31 β版

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

▶ エコール・ポリテクニーク・フェデラル・デ・ローザンヌ(イーピーエフエル)の特許一覧

<>
  • 特許6924530-モデルベース画像再構成方法 図000121
  • 特許6924530-モデルベース画像再構成方法 図000122
  • 特許6924530-モデルベース画像再構成方法 図000123
  • 特許6924530-モデルベース画像再構成方法 図000124
  • 特許6924530-モデルベース画像再構成方法 図000125
  • 特許6924530-モデルベース画像再構成方法 図000126
  • 特許6924530-モデルベース画像再構成方法 図000127
  • 特許6924530-モデルベース画像再構成方法 図000128
  • 特許6924530-モデルベース画像再構成方法 図000129
< >
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】6924530
(24)【登録日】2021年8月4日
(45)【発行日】2021年8月25日
(54)【発明の名称】モデルベース画像再構成方法
(51)【国際特許分類】
   A61B 8/14 20060101AFI20210812BHJP
【FI】
   A61B8/14ZDM
【請求項の数】14
【全頁数】34
(21)【出願番号】特願2020-511196(P2020-511196)
(86)(22)【出願日】2018年8月21日
(65)【公表番号】特表2020-534043(P2020-534043A)
(43)【公表日】2020年11月26日
(86)【国際出願番号】EP2018072579
(87)【国際公開番号】WO2019038296
(87)【国際公開日】20190228
【審査請求日】2021年2月26日
(31)【優先権主張番号】17187412.6
(32)【優先日】2017年8月23日
(33)【優先権主張国】EP
【早期審査対象出願】
(73)【特許権者】
【識別番号】519424607
【氏名又は名称】エコール・ポリテクニーク・フェデラル・デ・ローザンヌ(イーピーエフエル)
【氏名又は名称原語表記】ECOLE POLYTECHNIQUE FEDERALE DE LAUSANNE(EPFL)
(74)【代理人】
【識別番号】100134430
【弁理士】
【氏名又は名称】加藤 卓士
(74)【代理人】
【識別番号】100133639
【弁理士】
【氏名又は名称】矢野 卓哉
(72)【発明者】
【氏名】エイドリアン・ベッソン
(72)【発明者】
【氏名】ディミトリス・ペルディオス
(72)【発明者】
【氏名】マルセル・アルディティ
(72)【発明者】
【氏名】ジャン=フィリップ・ティラン
【審査官】 伊知地 和之
(56)【参考文献】
【文献】 特開2017−104476(JP,A)
【文献】 特開2005−342140(JP,A)
【文献】 Guillaume David et al.,"Time domain compressive beam forming of ultrasound signals",The Journal of the Acoustical Society of America,2015年05月,Vol.137, No.5,pp.2773-2784
【文献】 Ece Ozkan et al.,"Inverse Problem of Ultrasound Beamforming With Sparsity Constraints and Regularization",IEEE TRANSACTIONS ON ULTRASONICS, FERROELECTRONICS, AND FREQUENCY CONTROL,2018年03月,Vol.65, No.3,pp.356-365
【文献】 Solivan A. Valente et al.,"An Assessment of Iterative Reconstruction Methods for Sparse Ultrasound Imaging",Sensors,2017年03月08日,Vol.17, No.533,pp.1-33
(58)【調査した分野】(Int.Cl.,DB名)
A61B 8/00 − 8/15
G01N 29/00 − 29/52
PubMed
(57)【特許請求の範囲】
【請求項1】
少なくとも部分的に画像を再構成すべき物体に対し、送信素子を用いてパルス波を送信するステップと、
反射率格子点を含む反射率格子上に定義された反射率を用いて特徴付けられた前記物体から反射され、測定格子点を含む測定格子上に定義された測定値を定義するエコー波形を、センサを用いて受信するステップと、
測定モデルを用いて前記測定値から反射率推定値を取得するための逆問題を定義するステップと、
前記逆問題を、データ不一致項を有する目的関数を含む最適化問題として表現するステップと、
前記最適化問題を、前記データ不一致項の寄与を含む連立方程式として表現するステップと、
前記連立方程式において所定の収束閾値に到達するまで反復を実行することにより、前記反射率の推定値を取得する取得ステップと、
を含むパルス波画像再構成方法であって、
前記反復のうちの少なくともいくつかは、前記測定格子点のそれぞれに対応する前記測定値の推定、および前記反射率格子点のそれぞれに対応する前記反射率の推定を少なくとも含み、
前記取得ステップは、
超曲面反射率を含む反射率超曲面を定義する、前記送信素子とセンサとの間の経路ごとに、反射率パラメトリック方程式を生成し、前記超曲面反射率を、前記反射率超曲面のそれぞれにおいて積分して、前記測定モデルから、前記測定格子点のそれぞれに対応する前記測定値の推定値を取得するステップと、
超曲面測定値を含む測定超曲面を定義する、前記送信素子とセンサとの間の経路ごとに、測定パラメトリック方程式を生成し、前記超曲面測定値を、前記測定超曲面のそれぞれにおいて積分して、前記測定モデルの随伴作用素から、前記反射率格子点のそれぞれに対応する前記反射率の推定値を取得するステップと、
を含むパルス波画像再構成方法。
【請求項2】
前記送信素子のそれぞれは、送信点とみなされ、前記センサのそれぞれは、受信点とみなされる、請求項1に記載のパルス波画像再構成方法。
【請求項3】
前記超曲面反射率を積分する前に、
前記反射率超曲面のそれぞれを離散化して前記超曲面反射率を導きだし、前記反射率格子点と略一致するように前記超曲面反射率を補間する請求項1または2に記載のパルス波画像再構成方法。
【請求項4】
前記超曲面測定値を積分する前に、前記測定超曲面のそれぞれを離散化して前記超曲面測定値を導き出し、前記測定格子点と実質的に一致するように前記超曲面測定値を補間する請求項1〜3のいずれか1項に記載のパルス波画像再構成方法。
【請求項5】
前記データ不一致項は、前記測定格子点のそれぞれについて前記測定値と前記測定値の推定値との間の距離を評価する正の汎関数として表わされる、請求項1〜4のいずれか1項に記載のパルス波画像再構成方法。
【請求項6】
前記正の汎関数は、微分可能な正の汎関数である、請求項5に記載のパルス波画像再構成方法。
【請求項7】
前記目的関数は、前記データ不一致項と少なくとも1つの画像事前分布項との組み合わせとして表わされ、前記画像事前分布項は、前記反射率に関する事前情報を含む、請求項1〜6のいずれか1項に記載のパルス波画像再構成方法。
【請求項8】
前記少なくとも1つの画像事前分布項は、所定のモデルにおいて前記反射率推定値をpとした場合に、lpのノルムのp乗として表わされ、前記lpのノルムは、ベクトル
【数1】
に対して
【数2】
として定義され、
【数3】
は、実数の空間を表す、請求項7に記載のパルス波画像再構成方法。
【請求項9】
前記パルス波は、パルス音響波または電磁放射波である、請求項1〜8のいずれか1項に記載のパルス波画像再構成方法。
【請求項10】
前記パルス波は、偏向平面波または発散波を含む、請求項1〜9のいずれか1項に記載のパルス波画像再構成方法。
【請求項11】
前記パルス波は、電気信号によって励起される少なくとも1つの電気機械変換装置によって生成され、前記測定値は、少なくとも1つの相互電気機械変換装置によって生成される電気信号から取得される、請求項1〜10のいずれか1項に記載のパルス波画像再構成方法。
【請求項12】
前記電気機械変換装置は、
略1波長ピッチを有する多素子構造を備え、直線、凸状線または凹状線に沿って整列されたリニアアレイプローブ、
略半波長ピッチを有する多素子構造を備え、直線、凸状線または凹状線に沿って整列されたフェーズドアレイプローブ、
略1波長ピッチを有する多素子構造を備え、平面上に整列されたマトリクスアレイプローブ、または、
略半波長ピッチを有する多素子構造を備え、平面上に整列されたマトリクスアレイプローブ、に空間的に配置される請求項11に記載のパルス波画像再構成方法。
【請求項13】
前記取得ステップにおいて、
前記超曲面測定値を積分する際、または前記超曲面反射率を積分する際には、数値積分を行ない、
前記測定格子点のそれぞれについて、取得された前記測定値の推定値を、先に取得された前記測定値の推定値に加算して、測定値推定値の累積値を取得し、
前記累積値を、送信された前記パルス波のパルス形状と畳み込む、請求項1〜1のいずれか1項に記載のパルス波画像再構成方法。
【請求項14】
物体の画像を再構成する画像処理装置であって、
少なくとも部分的に再構成すべき物体に対し、送信素子を用いてパルス波を送信する送信部と、
反射率格子点を含む反射率格子上に定義された反射率を用いて特徴付けられた前記物体から反射され、測定格子点を含む測定格子上に定義された測定値を定義するエコー波形を、センサを用いて受信する受信部と、
測定モデルを用いて前記測定値から反射率推定値を取得するための逆問題を定義する定義部と、
前記逆問題を、データ不一致項を有する目的関数を含む最適化問題として表現する第1表現部と、
前記最適化問題を、前記データ不一致項の寄与を含む連立方程式として表現する第2表現部と、
前記連立方程式において所定の収束閾値に到達するまで反復を実行することにより、前記反射率の推定値を取得する取得部と、
を備え、
前記反復のうちの少なくともいくつかは、前記測定格子点のそれぞれに対応する前記測定値の推定、および前記反射率格子点のそれぞれに対応する前記反射率の推定を少なくとも含み、
前記取得部は、
超曲面反射率を含む反射率超曲面を定義する、前記送信素子とセンサとの間の経路ごとに、反射率パラメトリック方程式を生成し、前記超曲面反射率を、前記反射率超曲面のそれぞれにおいて積分して、前記測定モデルから、前記測定格子点のそれぞれに対応する前記測定値の推定値を取得する第1取得部と、
超曲面測定値を含む測定超曲面を定義する、前記送信素子とセンサとの間の経路ごとに、測定パラメトリック方程式を生成し、前記超曲面測定値を、前記測定超曲面のそれぞれにおいて積分して、前記測定モデルの随伴作用素から、前記反射率格子点のそれぞれに対応する前記反射率の推定値を取得する第2取得部と、
を備えた画像処理装置。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、例えば超音波イメージングで使用するのに適するモデルベース画像再構成方法に関する。本発明はまた、方法を実行するためのイメージング装置に関する。
【背景技術】
【0002】
超音波(US)イメージング用途など、パルス波イメージング用途は、様々な技術分野において使用され得る。医療用途では、それは、腱、筋肉、関節、血管、内臓などの体内構造を見るために使用される。しかしながら、パルス波イメージング技術は、医療分野以外でも広く使用されている。超音波の1つの商業的使用は、非破壊検査(NDT)であり、それは、部品またはシステムの保守性を損なうことなく、材料、構成要素、または組み立て部品の不連続性、または特性の相違を検査、試験、または評価するプロセスである。言い換えると、検査または試験が完了しても、その部品は、引き続き使用することができる。NDTの例は、オイルパイプにある。この場合、超音波は、パイプを真っすぐに伝搬するように使用され、このようにして、パルスの反射として亀裂または欠陥があるかどうかを検出することができる。パイプ内の波の速度、および亀裂からの反射と波の送信と間の時間間隔の知識があれば、欠陥の位置を推定することができる。パイプ内の亀裂検出のため使用されるのと同じ理論を使用して、航空機構造/翼パネルなどの亀裂を検出すること、および例えば、翼の表面上の氷の蓄積を検出することができる。もちろん、パルス波イメージングのための他の多くの可能なNDT用途がある。パルス波イメージングは、その安全性、携帯性、およびリアルタイム性のため、人気を博している。
【0003】
従来のUSイメージングは、時間分離された複数の集束送信ビームを使用して、走査線に沿って画像を再構成する。したがって、フレームレートは、集束送信ビームの数によって制限され、1秒当たり数十画像を超えることはできない。このフレームレートは、多くの用途には十分であるが、エラストグラフィにおけるせん断波の伝播など、より複雑な現象の理解には、より高いフレームレートが必要である。
【0004】
この弱みを克服するための1つの方法は、送信ビームの数を可能な限り減らすことにある。合成開口(SA)法は、少数のトランスデューサ素子のみを使用して、媒体全体を連続的に超音波照射することにより、この問題を克服しようとする。もう1つの選択肢は、超高速US(UFUS)イメージングであり、それは、平面波(PW)または発散波(DW)を使用して、視野全体を一度に超音波照射するという着想を利用し、USシステムが1秒当たり数千フレームに到達することを可能にする。
【0005】
かかるアプローチの重要な弱みは、画質の劣化である。実際、エネルギーが特定の関心領域に集中する集束送信ビームと比較すると、PWまたはDWのエネルギーは、超音波照射された媒体に広がり、結果として信号対雑音比が低くなり、空間分解能が低下する。この問題に対処するための1つの方法は、PWイメージングの状況においては異なる超音波照射角で、UFUSのDWイメージングの状況においては異なる点源位置で、SAイメージングにおいてはトランスデューサ素子の異なるグループで取得された複数の画像を平均化することにある。UFUSにおいて、このプロセスは、コヒーレント複合と呼ばれる。かかる技術の実装は、かなり簡易であるが、複数の超音波照射を必要とし、したがってフレームレートが低下する。
【0006】
複数の超音波照射を平均化することに代わるものは、主に遅延和(DAS)ビームフォーミングを中心に展開する従来の技術よりも効率的な再構成方法を使用することにある。方法の1つの一般的なグループは、US画像再構成によって引き起こされる不適切な画像形成問題を解くための反復アルゴリズムの使用に依存している。これらの方法は、問題のフォワードモデルに基づいている。既知の解決策が有する主な問題は、通常、対応する行列表示のストレージ要件に変換される、計算の複雑さにある。提案された、いくつかのモデルは、2次元(2D)の行列係数に数百ギガバイトのストレージを必要とする。この問題は、従来のアプローチを考慮すると、反復法の魅力を大幅に制限する。
【発明の概要】
【0007】
本発明の目的は、パルス波イメージングにおける画像再構成に関連する上記で特定された問題の少なくともいくつかを克服することである。
【0008】
本発明の第1の態様によれば、請求項1に記載の画像再構成方法が提供される。
【0009】
提案された方法は、再構成された画像が非常に高品質であり、迅速に再構成できるという利点を有する。言い換えると、従来技術の解決策と比較して、画像のコントラストおよび解像度を大幅に向上させることができる。
【0010】
本発明の第2の態様によれば、請求項15に記載の本発明の第1の態様による方法を実行するように構成されたイメージング装置が提供される。
【0011】
本発明の他の態様は、添付の従属請求項に記載されている。
【図面の簡単な説明】
【0012】
本発明の他の特徴および利点は、添付の図面を参照して、非限定的な例示的実施形態の以下の説明から明らかになるであろう。
【0013】
図1】本発明の一実施例による、本発明の教示を理解するのに有用であるイメージングシステムのいくつかの素子を示すブロック図である。
図2】本発明の一実施例による、提案されたイメージング方法において使用されるイメージング構成を示す概略図である。
図3】本発明の一実施例による、提案されたイメージング方法において使用されるイメージング構成および連続領域反射率超曲面を示す概略図である。
図4】本発明の一実施例による、提案されたイメージング方法において使用されるイメージング構成および離散領域反射率超曲面を示す概略図である。
図5】本発明の一実施例による、画像再構成方法を要約するフローチャートである。
図6】本発明の一実施例による、画像再構成方法をより詳細に例解するフローチャートである。
図7】本発明の別の実施例による、画像再構成方法をより詳細に例解するフローチャートである。
図8】本発明の一実施例による、画像測定モデルを実装するプロセスを例解するフローチャートである。
図9】本発明の一実施例による、測定モデルの随伴作用素を実装するプロセスを例解するフローチャートである。
【発明を実施するための形態】
【0014】
次に、添付の図面を参照して、本発明の一実施形態を詳細に説明する。本実施形態は、USイメージングの文脈で説明されているが、本発明の教示は、本環境に限定されない。本発明の教示は、パルス波イメージングを使用することができる他の分野に等しく適用可能である。異なる図面に現れる同一または対応する機能的および構造的素子には、同じ参照番号が割り当てられている。
【0015】
本実施形態における本発明は、以下の形式の不適切な線形逆問題を解くことに基づいて、2Dおよび3D超音波イメージングの文脈における画像形成方法および装置またはシステムを提案する。
m=Hγ+n、
ここで、mは測定値、γは精査中の物体、またはより具体的にはその反射率関数、Hは測定モデル、nはノイズである。提案された画像形成方法は、2つの主要な柱、つまり、高速かつ行列を含まない測定モデルH、および測定値mが与えられた場合に精査中の物体γの推定値を検索することが可能になる画像再構成方法に依存している。物体は、より大きな構造の内部部品または素子である。提案された方法は、図1に示すように、USプローブ、画像形成モジュール、ならびに後処理および表示モジュールを含む標準USシステムによって実行することができる。
【0016】
図1は、本発明の教示を理解するのに有用である、イメージングシステムまたは装置1のいくつかの素子を示す簡略ブロック図である。システム1は、超音波プローブ3を含み、本実施例では、トランスデューサとも呼ばれる圧電素子またはセンサの線形アレイを備えている。かかるプローブは、2Dにおけるイメージングの目的に適しており、画像平面の近傍に仰角において音響エネルギーを集中させるための円筒形の集束スキームを含み得る。しかしながら、この種のプローブは、単なる一例である。提案された方法は、異なるプローブ形状(コンベックスプローブ、フェーズドアレイなど)または技術(ポリフッ化ビニリデン(PVDF)共重合体および容量性マイクロマシン超音波トランスデューサ(CMUT)などを含む)に容易に適応することができる。同様に、提案された方法は、ボリュームイメージングを提供するために音響ビームを送信し、超音波照射されたボリュームからエコーを収集するように設計された2D行列プローブを使用する場合に拡張することができ、それにより、イメージングプロセスは、精査中のボリュームの3D表現を生成する。
【0017】
USプローブ3は、画像形成を実行するように構成されている、画像形成モジュールまたは装置5に接続されている。測定モデル装置7において推定された測定モデルは、画像形成プロセスで使用される。後で詳細に説明する提案された測定モデルは、パルスエコー空間インパルス応答モデルに基づいている。本発明は、測定モデルおよびその随伴の行列を含まない定式化を導入し、作用素の随伴は、連続作用素の場合には正方行列の転置行列の連続拡張と定義され、それは、パルス波イメージングにおける新しい画像形成方法を表している。測定モデルの説明は、2つの主要なステップによって実現し得る。
●パルス波の往復飛行時間のパラメトリック方程式が導出される。モデルは、後で説明するように、指定される超曲面上での積分として再計算することができることが示される。
●次に、提案された積分は、後述のように離散化され得る。
【0018】
モデルは、例えば、PW、DW、およびSAイメージングのために導出され得る。それは、無線周波数および同相直交(IQ)データと互換性がある。最後に、提案された測定モデルの随伴の行列を含まない定式化が導出される。それはまた、2Dイメージングの場合は1D多様体として、3Dイメージングの場合は2D多様体として定義される超曲面上での積分として再計算することができ、それらは、パラメータ化されていること、および測定モデルのために提案されたものと同じ離散化スキームを使用し得ることも実証される。
【0019】
画像形成装置5はまた、画像再構成プロセスまたは方法を実行するための画像再構成装置9を含む。画像再構成方法は、精査中の画像の推定値
【数1】

を検索することを目指している。この推定値は、
【数2】

と書くことができ、
ここで、fは、ハイパーパラメータ
【数3】

を備える画像再構成プロセスの関数である。本実施形態では、以下の2つの画像再構成方法について説明する。
(1)
【数4】

が明示的な定式化(例えば、Hの随伴作用素を使用したバックプロパゲーションなど)を有する分析的アプローチ。
(2)
【数5】

が以下の問題を解く反復アルゴリズムである、段階的アプローチ。
【数6】


ここで、
【数7】

は、Hγとmとの間の誤差を測定するデータ不一致項とも呼ばれる汎関数Fを含む目的関数を示し、Rは、非負の汎関数(正則化)を示し、λは、正則化パラメータを示す。この場合、
【数8】

は、最小化問題を解くために使用される最適化アルゴリズムに関連するハイパーパラメータであり、例えば、反復回数、アルゴリズムの停止基準などである。
【0020】
画像再構成後、再構成されたデータは、後処理および表示装置またはモジュール11に送信される。後処理ステップは、Bモードイメージング、カラードプライメージング、ベクトルドプライメージング、エラストグラフィなど様々な用途を対象としている。
【0021】
RFデータを使用したBモードイメージングの場合、再構成されたデータに包絡線検波が適用される。例えば、包絡線検波は、ヒルベルト変換、次いで、振幅検出および任意選択的にローパスフィルタリングによって実現することができる。それはまた、信号を2乗して、ローパスフィルタリングすることによっても実現することができる。BモードイメージングがIQデータを用いて使用されると、信号の振幅が抽出される。包絡線検波ステップに続いて、正規化およびダイナミクス圧縮ステップが行われる。ドプライメージングおよびエラストグラフィの場合、再構成されたRFまたはIQデータは、後処理を伴わずに直接使用される。
【0022】
提案された方法を詳細に説明する前に、本記載で使用されるいくつかの表記を最初に説明する。本記載で考慮されるパルスエコー構成では、媒体
【数9】

からのエコーは、
【数10】

である反射率(関数)γ(r)によって特徴づけられ、センサ13によって検出される。図2は、PWイメージングの場合の特定の簡略化された測定構成を示している。以下の表記は、本記載全体を通して使用される。
[数学表記または類似表記]
ベクトルは太字で示されている。
【数11】

は、ベクトル
【数12】


lpノルムを示し、
【数13】

は、実数の空間を示す。
●ヒルベルト空間Wの2乗可積分関数の空間は、L(W)と示される。2つの関数
【数14】

間の内積は、
【数15】

と定義される。
【0023】
超音波関連の表記
●m(p,t)は、pに配置されたセンサ13によって時刻tに受信された電気信号である。言い換えると、m(p,t)は、測定格子点(p,t)において計算された測定値の連続的な集合の1つのサンプル(素子生データの1つのサンプルと同等)を定義する。
●γ(r)は、反射率格子点rにおいて計算された反射率(画像サンプルと同等)の1つの値を定義する。
●Ωは、対象となる媒体または物体、すなわち、反射率サンプルの場所の連続的な集合を示す。
●o(r,p)は、pに位置するセンサとrに位置する点と間の素子指向性である。
●a(r,p)は、pに位置する点源とrに位置する点との間の減衰の値を定義する。
●vpe(t)は、受信したパルス形状(すなわち、パルスエコー波形)である。
●tTx(r,q)は、送信伝搬遅延、すなわち、qに位置する送信デバイスまたは素子13によって送信されるパルス波がrに位置する点に到達するのにかかる時間を示す。
●tRx(r,p)は、受信伝搬遅延、すなわち、rに位置する点によって送信されるパルスエコー波がpに位置するセンサに到達するのにかかる時間を示す。
【0024】
測定格子関連の表記
●送信素子の集合、
【数16】

ここで、Neltは、送信素子の数を指し、送信機とも呼ばれる。
●受信センサの集合、
【数17】

ここで、Nelrは、受信センサの数を指し、受信機とも呼ばれる。
●送信素子内の点の集合:
【数18】


●受信センサ内の点の集合:
【数19】


●時間サンプルの集合:
【数20】

【0025】
反射率格子関連の表記
●反射率の場所の集合:
【数21】

ここで、Nγは、反射率の数を指す。
【0026】
提案された方法は、第1の格子および第2の格子を定義し、共に3D空間にある。第1の格子は、反射率格子と呼ばれ、それは、2つの交差する格子線の交点における反射率格子点を含み、それは、本例では、互いに対して実質的に90度の角度を形成する。反射率格子の全ての反射率格子点は、媒体Ωに属する。第2の格子は、測定格子と呼ばれ、それは、2つの交差する格子線の交点における測定格子点を含み、それは、本例では、互いに対して実質的に90度の角度を形成する。測定格子点は、1次元では受信センサの場所と一致し、他の次元では各センサのサンプリング時刻と一致する。
【0027】
次に、本発明による所与の測定格子点での測定モデルの推定、評価、計算または実装について説明する。本例では、測定モデルは、各測定格子点で推定される。本開示で使用される測定モデルは、G.E.Tupholme,"Generation of acoustic pulses by baffled plane pistons",Mathematika,vol.16,p.209,1969、およびP.R.Stepanishen,"The Time-Dependent Force and Radiation Impedance on a Piston in a Rigid Infinite Planar Baffle",J.Acoust.Soc.Am.,vol.49,p.76,1971によって最初に導入されたパルスエコー空間インパルス応答モデルである。本モデルによれば、位置pおよび時間tで記録される生データサンプルを説明する以下の方程式を書くことができる。
【0028】
次に、1つの送信点源素子および1つの受信点源センサの構成に関する測定モデルの連続パラメトリック定式化について説明する。
【0029】
測定モデルは、
【数22】

と定義される。
【0030】
上記の方程式は、以下に示すように表現されることができ、測定モデルを超曲面Γ(p、q、t)に沿った積分と定義する。
【数23】


ここで
【数24】
【0031】
上記の方程式では、超曲面、またはより具体的には反射率超曲面は、方程式g(r、p、q、t)=0によって定義されるΓ(p、q、t)であり、ここで、対象となる変数は、rである。Γ(p、q、t)は、rが2D平面(放物線または楕円)にある場合は円錐曲線を定義し、rが3Dボリュームにある場合は2次曲面を定義する。同じ形状は、測定超曲面にも適用する。図3は、連続領域における測定構成を概略的に例解し、反射率超曲面も示す。
【0032】
次に、反射率超曲面Γ(p、q、t)は、積分次元を減らすためにパラメータ化することができる。方程式g(r、p、q、t)=0によって定義される反射率超曲面Γ(p、q、t)は、反射率パラメトリック方程式の集合によって同等に定義されることができる。数学では、パラメトリック方程式は、数量の集合を1つ以上の独立変数(すなわち、パラメータ)の関数と定義する。パラメトリック方程式は、一般に、面または曲線などの幾何学的対象を構成する点の座標を表すために使用され、その場合、方程式は、物体のパラメトリック表現またはパラメータ化と総称される。反射率パラメトリック方程式の集合は、ここで、以下、
【数25】

と定義されることができる。
【0033】
rは、パラメータα∈A、に依存し、ここで、
【数26】

ならば
【数27】

【数28】

の場合
【数29】

であるため、同値の右側の方程式は、rの反射率パラメトリック方程式を定義する。測定値またはサンプルごとに反射率パラメトリック方程式の1つの集合がある場合は留意する必要がある。さらに、反射率パラメトリック方程式の集合を描画することにより、超曲面を取得することができる。反射率パラメトリック方程式の集合の上記定式化を備えているため、以下ように測定モデルHの連続パラメトリック定式化を導出することができる。
【数30】

ここで、
【数31】

【0034】
次に、任意形状の1つの送信有限長素子および1つの受信点源センサの構成に関する連続測定モデルについて説明する。この場合、送信素子は、有限の長さおよび任意の形状を有する。それは、点源場所の連続的な集合Χによって特徴付けられる。測定モデルは、
【数32】

と定義される。
【0035】
1つの送信点源素子および1つの受信点源センサの構成に関連して記載したモデルとの主な相違は、任意形状の有限長センサを形成する全ての点源の寄与の和に対応する、もう1つの積分があるという事実にあることがわかる。
【0036】
対応するパラメトリック定式化は、
【数33】

によって与えられる。
【0037】
次に、任意形状の1つの送信有限長素子および1つの受信有限長センサの構成に関する連続測定モデルについて説明する。この場合、受信センサは、有限の長さおよび任意の形状を有する。それは、点源場所の集合Πiによって定義される。結果として得られた信号の座標は、
【数34】

によって示され、ここで、
【数35】

である。それは、例えば、Πiの中心の点の座標であり得る。
【0038】
それは、
【数36】

と書くことができる。結果として得られた測定値は、Πiによって定義されたトランスデューサ面に沿って積分することにより取得される。m(p,t)を定義する方程式を上記の積分に統合すると、次の関係、
【数37】

が成り立つ。
【0039】
対応するパラメトリック定式化は、
【数38】

によって与えられる。
【0040】
次に、任意形状の複数の送信有限長素子および任意形状の1つの受信有限長センサの構成に関する連続測定モデルについて説明する。この場合、点源場所の集合Χは、送信素子jに関連付けられており、以下の式、
【数39】

が導かれる。
【0041】
この方程式は、任意形状の1つの送信有限長素子および任意形状の1つの受信有限長センサの構成に関連して記載した方程式と非常に類似している。唯一の相違は、複数の送信素子を考慮するために導入された和にある。
【0042】
対応するパラメトリック定式化は、
【数40】

によって与えられる。
【0043】
次に、任意形状の複数の送信有限長素子および任意形状の複数の受信有限長センサの構成に関する連続測定モデルについて説明する。
【0044】
複数の受信センサの集合を使用する場合、点源場所の集合Πiおよび座標ξiは、受信センサiに関連付けられている。したがって、受信測定領域は、
【数41】

と定義することができ、
それは、離散座標ξiについてのみ定義されている。本記載の残りの部分では、以下
【数42】

が成り立つ。
【0045】
次に、測定モデルの連続パラメトリック定式化を離散化する。このステップは、任意形状の複数の送信有限長素子および任意形状の1つの受信有限長センサの構成に関連して記載した測定モデルのパラメトリック定式化を離散化することにある。しかしながら、代わりに上記のモデルのいずれかを離散化することができる。離散化とは、以下のパラメータの集合
【数43】

を定義するパラメータαのいくつかの離散値αを選択することを意味する。パラメータαの離散値の選択が、次にNα
【数44】

によってのみ記載される超曲面Γ(p,q,t)の離散化をもたらすことがわかる。図4は、離散領域における測定構成を概略的に例解し、反射率超曲面も示している。測定モデルのパラメトリック定式化の離散化は、
【数45】

として表現することができ、
ここで、
【数46】

は、測定格子の各時間サンプルにおいて評価されるパルス形状であり、
【数47】

であり、
【数48】

であり、
ここで、w(q),u(α)およびz(p)は、積分の重み(連続積分の離散化に関連する)である。送信素子の各点qについて、上記式は、反射率超曲面の点上に推定される、超曲面反射率サンプルまたは値と呼ばれる、反射率の集合
【数49】

を含むことに気づくことができる。
【0046】
次に、超曲面反射率は、反射率格子点と実質的に一致するように補間される。超曲面反射率は、反射率格子点上にないため、通常は未知である。したがって、それらは、P.Thevenaz,T.Blu and M.Unser,"Interpolation revisited",IEEE Transactions on Medical Imaging,vol.19,p.739−758,2000によって導入される補間カーネルによって、反射率格子上で定義された反射率サンプルから近似する必要がある。したがって、
【数50】

であるような補間カーネル
【数51】

を導入する。
【0047】
上記補間関数は、超曲面反射率
【数52】

(通常、反射格子点上にない)を反射率サンプルγSに関連付ける。
【0048】
次に、離散化および補間された超曲面反射率サンプルまたは値の総和が実行される。補間関数は、離散化された測定モデルに挿入され、次の総和、
【数53】

が実行される。
【0049】
本例では、先の方程式に記載したように、考慮された測定格子点において測定値推定値とも呼ばれる測定値の推定値を取得するために、最終的に総和の出力がパルス形状に畳み込まれる。上記のプロセスは、全ての測定値に対して実行される。したがって、測定モデルの推定または処理の結果は、それぞれの測定格子点における測定値の推定値である。
【0050】
次に、各反射率格子点における測定モデルの随伴作用素の推定、評価、計算、または実装が実行される。まず、測定モデルの随伴作用素の連続パラメトリック定式化が取得される。随伴作用素の連続定式化は、機能解析ツールを使用して、測定モデルから、すなわち上記で与えられた連続測定モデルのいずれかから導出することができる。以下の方程式は、任意形状の複数の有限長送信素子および任意形状の複数の有限長受信センサの場合から導出される。実際に、
【数54】

であるような作用素
【数55】

を定義する場合、Hの随伴作用素
【数56】

は、
【数57】

と定義される。
【0051】
以下の内積、
【数58】

を考慮し、
ここで、
【数59】

は、パルスエコー波形の整合フィルタである。次に、
【数60】

を有し、
ここで、
【数61】

【0052】
測定モデルの随伴作用素を備えているため、測定値の推定値が与えられた場合、反射率格子点rにおける反射率の推定値を、
【数62】

と定義する。
【0053】
この場合は、測定モデルよりも簡易である。実際、パラメトリック定式化は、直接的であり、測定パラメトリック方程式は、簡易な飛行時間方程式、
【数63】

として表現することができる。
【0054】
次いで、上記の方程式に対応する測定超曲面は、
【数64】

である。
【0055】
次に、測定モデルの随伴作用素のパラメータ定式化を離散化する。測定モデルの離散化の場合と同じプロセスが適用される。この場合、離散化プロセスは、上に直接適用されるため、より直接的である。離散化されたパラメータの集合は、
【数65】

である。したがって、以下の方程式、
【数66】

を有し、
ここで、w(q)およびZ(p)は、積分の重みであり、
【数67】

は、各時間サンプルで評価したパルス形状の整合フィルタである。上記の方程式は、超曲面測定サンプルまたは値と呼ばれる、測定超曲面上にある、畳み込み測定サンプルの集合
【数68】

を含むことに気付くことができる。
【0056】
次に、超曲面測定サンプルの補間は、測定格子点と実質的に一致するように実行される。超曲面測定サンプルの値は、測定格子点上にないため、通常は未知である。そのため、補間カーネルによって、測定格子上で定義された測定サンプルから近似する必要がある。
【数69】

であるような、補間カーネル
【数70】

を導入する。
【0057】
上記の方程式は、超曲面測定サンプル
【数71】

を畳み込み測定サンプル
【数72】

に関連付ける。
【0058】
この後、離散化および補間された超曲面測定サンプルが合計される。補間方程式は、離散化された測定モデルの随伴作用素に挿入され、以下の総和、
【数73】

が実行される。
【0059】
総和の出力は、考慮された反射率格子点における反射率の推定値である。
【0060】
次に画像再構成手順を説明する。最初に逆問題を定義する。測定値mは、上記の測定モデルの随伴作用素を使用して、反射率の推定値
【数74】

を生成するために使用することができる。同様に、反射率の推定値
【数75】

を使用して、前記の測定モデルを使用して、mに等しくない測定値の推定値
【数76】

を生成することができる。したがって、測定値の推定値と反射率の推定値と間の関係がわかる。
【0061】
画像再構成方法は、測定値mの集合が与えられた場合、反射率格子点において評価される反射率の推定値
【数77】

を検索することを目的とする。測定値と未知の反射率との関係は、以下の線形逆問題、
【数78】

によって定義され、
それは、測定格子の各点で推定される測定モデルに対応する。
【0062】
さらに、作用素
【数79】

は、γに対して線形であるため(所与の格子点において測定モデルを評価するときに表わされた式から推定することができる)、
m=Hγであるような行列
【数80】

(線形作用素Hに関連付けられている)が存在し、それは、
逆問題を定義する。
【0063】
次に、逆問題に関連付けられる最適化問題について説明する。上記の逆問題を解くために、それは、以下の最適化問題、
【数81】

として表現され、
ここで、J(γ,m)は、最適化問題に関与する目的関数を示し、F(Hγ,m)は、データ不一致項を考慮する下半連続汎関数であり、測定値の推定値Hγと測定値mとの間の距離を測定し、R(γ)は、任意選択の事前分布項を記述する下半連続汎関数であり、それは、反射率に関する特定の統計的挙動など追加情報を考慮し、λ>0は、正則化パラメータである。
【0064】
データ不一致項に使用される汎関数のいくつかの例を以下に列挙する。
●ユークリッド距離、
【数82】


●半径
【数83】


【数84】

球体に関する指標関数、
【数85】



【数86】

ノルム、
【数87】

【0065】
事前項に使用される汎関数のいくつかの例を以下に列挙する。

【数88】

ノルムのp乗、
【数89】


●p<1の場合、それは、γのスパース性の尺度である。
●p∈[1,2]の場合、それは、データが一般化ガウス分布(GGD)にどれだけ適合するかについての尺度である。
●所与のモデル
【数90】

における
【数91】

ノルムのp乗:
【数92】

●ψは、ウェーブレット変換またはフーリエ変換など一般的な変換であり得る。
●ψは、学習された辞書であり得る。
【0066】
次に、最適化問題を1つの画像事前分布でどのように解くことができるかについて説明する。FおよびRの特性に応じて、異なる種類のアルゴリズムを使用し得る。
【0067】
次に、FおよびRが微分可能である場合について説明する。この場合、逆問題は、γに関する目的関数J(γ,m)の導関数の根を見つけることによって解かれ、
【数93】

と書くことができる。
【0068】
問題がF(Hγ,m)の導関数の計算にHの計算を含むことがわかる。関数F(Hγ,M)に応じて、導関数の計算はまた、測定モデルHの随伴作用素の計算を含み得る。かかる問題の1つの有名な例は、Tikhonov正則化であり、ここで、
【数94】

である。この場合、以下の解
【数95】

を有し、ここで、
【数96】

は、単位行列、つまり、主対角線上は1、他の場所はゼロである、Nγ×Nγ正方行列である。
【0069】
次に、Fのみが微分可能であり、Rが凸である場合について説明する。かかる問題を解くために使用される方法の1つの一般的なグループは、射影勾配法と呼ばれる。射影勾配法は、最適化問題の解が固定点方程式を満たすという事実を利用する。1つの一般的なアルゴリズムは、以下の固定点方程式のシステムを含む前方−後方分離であり、本具体例では、1つの方程式に限定されている。
【数97】


ここで、
【数98】

は、Rに関連付けられた近接作用素を示し、近接作用素は、P.Combettes and J.-C.Pesquet,"Proximal Splitting Methods in Signal Processing", Fixed-Point Algorithms for Inverse Problems in Science and Engineering,p.185-212,2011によって導入されている。上記の方程式は、方程式の両側に
【数99】

が現れるため、固定点方程式として示される。したがって、収束基準に達するまで、以下の反復、
【数100】

を実行することによって問題を解く。
【0070】
別の一般的なアルゴリズムは、高速反復縮小閾値アルゴリズム(FISTA)であり、A.Beck and M.Teboulle,"A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems",SIAM Journal on Imaging Science,vol.2,p.183-202,2009によって導入され、それは、上記の方法の変形である。
【0071】
前方−後方分離の1つの例は、
【数101】

の場合である。この場合、以下の反復が実行される。
【数102】


ここで、
【数103】

は、ハイパーパラメータである。
【0072】
実際、
【数104】

ノルムに関連付けられた近接作用素は、
【数105】

と定義され、ここで、ソフト閾値作用素は、素子ごとにsign(x)max(|x|−τ、0)と定義され、ここで、sign(x)は、x>0の場合は1、x=0の場合は0、それ以外の場合は−1である。
【0073】
次に、FおよびRが微分可能でない場合について説明する。この場合、S.Boyd,N.Parikh,E.Chu,B.Peleato and J.Eckstein,"Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers",Foundations and Trends in Machine Learning,vol.3,p.1-122,2011によって導入された交互方向乗数法(ADMM)、またはP.Combettes,L.Condat,J.-C.Pesquet and BC Vu,"A forward-backward view of some primal-dual optimization methods in image recovery",Proceedings of the 2014 IEEE International Conference on Image Processing,p.4141-4145,2014によって導入されたプライマルデュアルフォワードバックワード(PDFB)アルゴリズムなど、いくつかの方法を使用することができる。
【0074】
次に、最適化問題をM個の画像事前分布でどのように解くことができるかについて説明する。この場合、事前分布項は、
【数106】

と、書くことができると考えられることができる。
【0075】
この場合、プライマルデュアルフォワードバックワード、N.Pustelnik,C.Chaux and J.-C.Pesquet,"Parallel Proximal Algorithm for Image Restoration Using Hybrid Regularization",IEEE Transactions on Image Processing,vol.20,p.2450-2462,2011によって導入された並列近接アルゴリズム(PPXA)、または、S.Setzer,G.Steidl and T.Teuber,"Deblurring Poissonian Images by Split Bregman Techniques",Journal of Visual Communication and Image Representation,vol.21,p.193-199,2010によって導入された同時方向乗数法(SDMM)など、いくつかのアルゴリズムを使用することができる。これらのアルゴリズムの共通の特性は、F(Hγ,m)(またはF(Hγ,m)の導関数)およびM個の画像事前分布項の近接作用素の計算を必要とすることである。F(Hγ,m)の導関数または近接作用素の計算は、関数F(Hγ,m)に応じて、測定モデルおよび測定の随伴作用素の推定を含み得る。
【0076】
上記の画像再構成方法は、図5〜9のフローチャートによって要約することができる。図5のフローチャートは、方法の概要を示す。ステップ21では、パルス波の集合が、送信素子の集合によって送信される。本例では、送信素子は、有限のサイズおよび任意の形状を有する。送信されたパルスの集合は、1つのパルスのみ、または連続して送信される2つ以上のパルスを含み得る。ステップ23では、反射されたエコー波形の集合が、本例では有限のサイズおよび任意の形状を有するセンサの集合によって受信される。ステップ25では、測定値の集合を反射率の集合に関連付ける逆問題が定義される。ステップ27では、逆問題の解が、データ不一致項および1つ以上の画像事前分布項を含む目的関数を含む最適化問題として表現される。ステップ29では、最適化問題が、射影を使用して解かれ、反射率の推定値を取得する。
【0077】
図6のフローチャートは、汎関数Fが微分可能である場合の提案された方法の概要をより詳しく示す。最初の4つのステップ、すなわちステップ31、33、35、および37は、図6のフローチャートの最初の4つのステップと同じである。ステップ39では、第1の反射率推定値
【数107】

が、初期化として測定値mについての測定モデルの随伴作用素を推定することにより計算される。ステップ41では、変数
【数108】

が、図8のフローチャートに従って計算または推定される。ステップ43では、導関数
【数109】

が、計算される。特定の実装例では、ステップ45において、随伴作用素を使用してデータ不一致項の導関数を計算する。例えば、
【数110】

の場合、F(Hγ,M)の導関数は、
【数111】

として表現され、測定モデルの随伴作用素の計算を含む。ステップ47では、画像事前分布項の近接作用素が計算される。ステップ49では、
【数112】

が更新される。ステップ50では、収束基準が満たされているか否かが判定される。肯定の場合、プロセスは、終了する。基準が満たされていない場合、プロセスは、ステップ41に続く。
【0078】
図7のフローチャートは、データ不一致項がユークリッド距離であり、逆問題が1つの画像事前分布項を含む場合の提案された方法の概要を、より詳細に示す。最初の5つのステップ、すなわちステップ51、53、55、57、59は、図6のフローチャートの最初の5つのステップと実質的に同じであるが、ここでは、逆問題が1つの画像事前分布項を含むという相違を有する。ステップ61では、変数
【数113】

が、図8のフローチャートに従って、各測定格子点に対して計算される。ステップ63では、残差r=v−mが計算される。ステップ65では、残差
【数114】

についての測定モデルの随伴作用素が図10のフローチャートに従って推定される。ステップ67では、値
【数115】

が計算される。ステップ69では、画像事前分布項への射影
【数116】

が計算される。ステップ70では、収束基準が満たされているかどうかが判定される。肯定の場合、プロセスは、終了する。基準が満たされてない場合、プロセスは、ステップ61に続く。
【0079】
図8のフローチャートは、測定モデルを推定または処理するプロセスの例を説明する。ステップ71では、反射率パラメトリック方程式の集合が、各測定格子点(任意形状の有限サイズ受信センサを定義する)に対して、測定格子点に対応する任意形状の有限サイズ受信機内の各点源場所に対して、任意形状の各有限サイズ送信機に対して、および任意形状の有限サイズ送信機内の各点源場所に対して、生成される。言い換えると、反射率パラメトリック方程式の集合は、各送信機および各受信機内の点源の全てのペア間の各送信機−受信機波経路に対して生成される。反射率パラメトリック方程式は、反射率超曲面の集合を定義する。ステップ73では、反射率超曲面が離散化されて、超曲面反射率サンプルの集合が取得される。ステップ75では、超曲面反射率サンプルが反射率格子上に補間される。ステップ77では、補間され、および離散化された超曲面反射率サンプルが合計される。ステップ79では、ステップ77で取得された値が、測定格子の点における測定値推定値に追加される。ステップ80では、取得された累積測定値推定値がパルス形状に畳み込まれる。
【0080】
図9のフローチャートは、測定モデルの随伴作用素を推定する、または処理するプロセスの例を説明する。ステップ81では、測定値推定値がパルス形状の整合フィルタに畳み込まれる。ステップ83では、測定パラメトリック方程式の集合が、反射率格子の各点に対して、任意形状の各有限サイズ送信機に対して、任意形状の各有限サイズ受信機に対して、任意形状の有限サイズ送信機内の各点源場所に対して、および任意形状の有限サイズ受信機に内の各点源場所に対して、生成される。言い換えると、測定パラメトリック方程式の集合は、各送信機および各受信機内の点源の全てのペア間の各送信機−受信機経路に対して生成される。測定パラメトリック方程式は、測定超曲面の集合を定義する。ステップ85では、測定超曲面が離散化されて、超曲面測定サンプルの集合が取得される。ステップ87において、超曲面測定サンプルが測定格子上に補間される。ステップ89では、補間され、および離散化された超曲面測定サンプルが合計される。ステップ90では、ステップ89で取得された値が、反射率格子の点における反射率推定値に追加される。
【0081】
送信されたパルス波は、電気信号によって励起される少なくとも1つの電気機械変換装置によって生成されることができる。測定値は、少なくとも1つの相互電気機械変換装置によって生成される電気信号から取得することができる。電気機械変換装置は、実質的に1波長ピッチを有する多素子構造を備え、直線、凸状線または凹状線に沿って整列された1つ以上のリニアアレイプローブ、実質的に半波長ピッチを有する多素子構造を備え、直線、凸状線または凹状線に沿って整列された1つ以上のフェーズドアレイプローブ、実質的に1波長ピッチを有する多素子構造を備え、平面上に整列された1つ以上の行列アレイプローブ、または実質的に半波長ピッチを有する多素子構造を備え、平面上に整列された1つ以上の行列アレイプローブ、に空間的に配置される。
【0082】
本発明は、図面および前述の記載で詳細に例解および記載されたが、かかる例解および記載は、事例的または例示的であり、限定的ではないとみなされるべきであり、本発明は、開示された実施形態に限定されない。他の実施形態および変形形態は、図面、開示および添付の特許請求の範囲の研究に基づいて、請求された発明を実施する際に当業者によって、理解され、達成されることができる。
【0083】
例えば、本発明の教示は、送信音響波が電磁放射パルスで置き換えられる、光音響イメージングに適用することができる。この特定の場合、電磁放射は、光の速度に近い速度、すなわち通常、音の速度よりも5桁速い速度で伝搬するため、全ての音響送信伝搬遅延をゼロに設定することで、本記載で導入される形式全体を簡素化することができる。このため、受信時の音響伝搬遅延を考慮する場合、送信伝搬遅延は無視することができる。本発明で考慮される別の変形は、容量性マシン超音波トランスデューサ(CMUT)など、圧電材料とは異なるトランスデューサ技術の使用にある。さらに別の変形は、本発明を時間多重方式に適用することであり、それにより送信および/または受信トランスデューサ素子のサブセットが連続して順次アドレス指定され、結果として得られる信号が後続の再結合のために記録される。上記の用途は全て、パルスインバージョン、振幅変調、振幅変調パルスインバージョン、ハーモニックイメージングなど、当技術分野で知られているマイクロバブルベースのコントラスト固有の非線形イメージング技術に変換することもできる。
【0084】
本発明はまた、イメージング装置1など、電子データ処理装置のコンピューティング手段上でロードおよび実行される場合、上記で説明した方法のステップを実装するための非一時的媒体に記憶された命令を含むコンピュータプログラム製品を提案する。
【0085】
特許請求の範囲において、「含む」という語は、他の素子またはステップを除外せず、不定冠詞「a」または「an」は複数を除外しない。異なる特徴が相互に異なる従属請求項に記載されているという単なる事実は、これらの特徴の組み合わせを有利に使用することができないことを示していない。
図1
図2
図3
図4
図5
図6
図7
図8
図9