JP6924530B2

JP6924530B2 - モデルベース画像再構成方法

Info

Publication number: JP6924530B2
Application number: JP2020511196A
Authority: JP
Inventors: エイドリアン・ベッソン; ディミトリス・ペルディオス; マルセル・アルディティ; ジャン＝フィリップ・ティラン
Original assignee: Ecole Polytechnique Federale de Lausanne EPFL
Current assignee: Ecole Polytechnique Federale de Lausanne EPFL
Priority date: 2017-08-23
Filing date: 2018-08-21
Publication date: 2021-08-25
Anticipated expiration: 2038-08-21
Also published as: KR102326149B1; IL272819A; JP2020534043A; EP3447525A1; IL272819B; KR20200040301A; US20200225335A1; CN111213066B; EP3447525B1; WO2019038296A1; ES2812649T3; US11175392B2; CN111213066A

Description

本発明は、例えば超音波イメージングで使用するのに適するモデルベース画像再構成方法に関する。本発明はまた、方法を実行するためのイメージング装置に関する。

超音波（ＵＳ）イメージング用途など、パルス波イメージング用途は、様々な技術分野において使用され得る。医療用途では、それは、腱、筋肉、関節、血管、内臓などの体内構造を見るために使用される。しかしながら、パルス波イメージング技術は、医療分野以外でも広く使用されている。超音波の１つの商業的使用は、非破壊検査（ＮＤＴ）であり、それは、部品またはシステムの保守性を損なうことなく、材料、構成要素、または組み立て部品の不連続性、または特性の相違を検査、試験、または評価するプロセスである。言い換えると、検査または試験が完了しても、その部品は、引き続き使用することができる。ＮＤＴの例は、オイルパイプにある。この場合、超音波は、パイプを真っすぐに伝搬するように使用され、このようにして、パルスの反射として亀裂または欠陥があるかどうかを検出することができる。パイプ内の波の速度、および亀裂からの反射と波の送信と間の時間間隔の知識があれば、欠陥の位置を推定することができる。パイプ内の亀裂検出のため使用されるのと同じ理論を使用して、航空機構造／翼パネルなどの亀裂を検出すること、および例えば、翼の表面上の氷の蓄積を検出することができる。もちろん、パルス波イメージングのための他の多くの可能なＮＤＴ用途がある。パルス波イメージングは、その安全性、携帯性、およびリアルタイム性のため、人気を博している。

従来のＵＳイメージングは、時間分離された複数の集束送信ビームを使用して、走査線に沿って画像を再構成する。したがって、フレームレートは、集束送信ビームの数によって制限され、１秒当たり数十画像を超えることはできない。このフレームレートは、多くの用途には十分であるが、エラストグラフィにおけるせん断波の伝播など、より複雑な現象の理解には、より高いフレームレートが必要である。

この弱みを克服するための１つの方法は、送信ビームの数を可能な限り減らすことにある。合成開口（ＳＡ）法は、少数のトランスデューサ素子のみを使用して、媒体全体を連続的に超音波照射することにより、この問題を克服しようとする。もう１つの選択肢は、超高速ＵＳ（ＵＦＵＳ）イメージングであり、それは、平面波（ＰＷ）または発散波（ＤＷ）を使用して、視野全体を一度に超音波照射するという着想を利用し、ＵＳシステムが１秒当たり数千フレームに到達することを可能にする。

かかるアプローチの重要な弱みは、画質の劣化である。実際、エネルギーが特定の関心領域に集中する集束送信ビームと比較すると、ＰＷまたはＤＷのエネルギーは、超音波照射された媒体に広がり、結果として信号対雑音比が低くなり、空間分解能が低下する。この問題に対処するための１つの方法は、ＰＷイメージングの状況においては異なる超音波照射角で、ＵＦＵＳのＤＷイメージングの状況においては異なる点源位置で、ＳＡイメージングにおいてはトランスデューサ素子の異なるグループで取得された複数の画像を平均化することにある。ＵＦＵＳにおいて、このプロセスは、コヒーレント複合と呼ばれる。かかる技術の実装は、かなり簡易であるが、複数の超音波照射を必要とし、したがってフレームレートが低下する。

複数の超音波照射を平均化することに代わるものは、主に遅延和（ＤＡＳ）ビームフォーミングを中心に展開する従来の技術よりも効率的な再構成方法を使用することにある。方法の１つの一般的なグループは、ＵＳ画像再構成によって引き起こされる不適切な画像形成問題を解くための反復アルゴリズムの使用に依存している。これらの方法は、問題のフォワードモデルに基づいている。既知の解決策が有する主な問題は、通常、対応する行列表示のストレージ要件に変換される、計算の複雑さにある。提案された、いくつかのモデルは、２次元（２Ｄ）の行列係数に数百ギガバイトのストレージを必要とする。この問題は、従来のアプローチを考慮すると、反復法の魅力を大幅に制限する。

本発明の目的は、パルス波イメージングにおける画像再構成に関連する上記で特定された問題の少なくともいくつかを克服することである。

本発明の第１の態様によれば、請求項１に記載の画像再構成方法が提供される。

提案された方法は、再構成された画像が非常に高品質であり、迅速に再構成できるという利点を有する。言い換えると、従来技術の解決策と比較して、画像のコントラストおよび解像度を大幅に向上させることができる。

本発明の第２の態様によれば、請求項１５に記載の本発明の第１の態様による方法を実行するように構成されたイメージング装置が提供される。

本発明の他の態様は、添付の従属請求項に記載されている。

本発明の他の特徴および利点は、添付の図面を参照して、非限定的な例示的実施形態の以下の説明から明らかになるであろう。

本発明の一実施例による、本発明の教示を理解するのに有用であるイメージングシステムのいくつかの素子を示すブロック図である。本発明の一実施例による、提案されたイメージング方法において使用されるイメージング構成を示す概略図である。本発明の一実施例による、提案されたイメージング方法において使用されるイメージング構成および連続領域反射率超曲面を示す概略図である。本発明の一実施例による、提案されたイメージング方法において使用されるイメージング構成および離散領域反射率超曲面を示す概略図である。本発明の一実施例による、画像再構成方法を要約するフローチャートである。本発明の一実施例による、画像再構成方法をより詳細に例解するフローチャートである。本発明の別の実施例による、画像再構成方法をより詳細に例解するフローチャートである。本発明の一実施例による、画像測定モデルを実装するプロセスを例解するフローチャートである。本発明の一実施例による、測定モデルの随伴作用素を実装するプロセスを例解するフローチャートである。

次に、添付の図面を参照して、本発明の一実施形態を詳細に説明する。本実施形態は、ＵＳイメージングの文脈で説明されているが、本発明の教示は、本環境に限定されない。本発明の教示は、パルス波イメージングを使用することができる他の分野に等しく適用可能である。異なる図面に現れる同一または対応する機能的および構造的素子には、同じ参照番号が割り当てられている。

本実施形態における本発明は、以下の形式の不適切な線形逆問題を解くことに基づいて、２Ｄおよび３Ｄ超音波イメージングの文脈における画像形成方法および装置またはシステムを提案する。
ｍ＝Ｈγ＋ｎ、
ここで、ｍは測定値、γは精査中の物体、またはより具体的にはその反射率関数、Ｈは測定モデル、ｎはノイズである。提案された画像形成方法は、２つの主要な柱、つまり、高速かつ行列を含まない測定モデルＨ、および測定値ｍが与えられた場合に精査中の物体γの推定値を検索することが可能になる画像再構成方法に依存している。物体は、より大きな構造の内部部品または素子である。提案された方法は、図１に示すように、ＵＳプローブ、画像形成モジュール、ならびに後処理および表示モジュールを含む標準ＵＳシステムによって実行することができる。

図１は、本発明の教示を理解するのに有用である、イメージングシステムまたは装置１のいくつかの素子を示す簡略ブロック図である。システム１は、超音波プローブ３を含み、本実施例では、トランスデューサとも呼ばれる圧電素子またはセンサの線形アレイを備えている。かかるプローブは、２Ｄにおけるイメージングの目的に適しており、画像平面の近傍に仰角において音響エネルギーを集中させるための円筒形の集束スキームを含み得る。しかしながら、この種のプローブは、単なる一例である。提案された方法は、異なるプローブ形状（コンベックスプローブ、フェーズドアレイなど）または技術（ポリフッ化ビニリデン（ＰＶＤＦ）共重合体および容量性マイクロマシン超音波トランスデューサ（ＣＭＵＴ）などを含む）に容易に適応することができる。同様に、提案された方法は、ボリュームイメージングを提供するために音響ビームを送信し、超音波照射されたボリュームからエコーを収集するように設計された２Ｄ行列プローブを使用する場合に拡張することができ、それにより、イメージングプロセスは、精査中のボリュームの３Ｄ表現を生成する。

ＵＳプローブ３は、画像形成を実行するように構成されている、画像形成モジュールまたは装置５に接続されている。測定モデル装置７において推定された測定モデルは、画像形成プロセスで使用される。後で詳細に説明する提案された測定モデルは、パルスエコー空間インパルス応答モデルに基づいている。本発明は、測定モデルおよびその随伴の行列を含まない定式化を導入し、作用素の随伴は、連続作用素の場合には正方行列の転置行列の連続拡張と定義され、それは、パルス波イメージングにおける新しい画像形成方法を表している。測定モデルの説明は、２つの主要なステップによって実現し得る。
●パルス波の往復飛行時間のパラメトリック方程式が導出される。モデルは、後で説明するように、指定される超曲面上での積分として再計算することができることが示される。
●次に、提案された積分は、後述のように離散化され得る。

モデルは、例えば、ＰＷ、ＤＷ、およびＳＡイメージングのために導出され得る。それは、無線周波数および同相直交（ＩＱ）データと互換性がある。最後に、提案された測定モデルの随伴の行列を含まない定式化が導出される。それはまた、２Ｄイメージングの場合は１Ｄ多様体として、３Ｄイメージングの場合は２Ｄ多様体として定義される超曲面上での積分として再計算することができ、それらは、パラメータ化されていること、および測定モデルのために提案されたものと同じ離散化スキームを使用し得ることも実証される。

画像形成装置５はまた、画像再構成プロセスまたは方法を実行するための画像再構成装置９を含む。画像再構成方法は、精査中の画像の推定値

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

と書くことができ、
ここで、ｆは、ハイパーパラメータ

を備える画像再構成プロセスの関数である。本実施形態では、以下の２つの画像再構成方法について説明する。
（１）

が明示的な定式化（例えば、Ｈの随伴作用素を使用したバックプロパゲーションなど）を有する分析的アプローチ。
（２）

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

ここで、

は、Ｈγとｍとの間の誤差を測定するデータ不一致項とも呼ばれる汎関数Ｆを含む目的関数を示し、Ｒは、非負の汎関数（正則化）を示し、λは、正則化パラメータを示す。この場合、

は、最小化問題を解くために使用される最適化アルゴリズムに関連するハイパーパラメータであり、例えば、反復回数、アルゴリズムの停止基準などである。

画像再構成後、再構成されたデータは、後処理および表示装置またはモジュール１１に送信される。後処理ステップは、Ｂモードイメージング、カラードプライメージング、ベクトルドプライメージング、エラストグラフィなど様々な用途を対象としている。

ＲＦデータを使用したＢモードイメージングの場合、再構成されたデータに包絡線検波が適用される。例えば、包絡線検波は、ヒルベルト変換、次いで、振幅検出および任意選択的にローパスフィルタリングによって実現することができる。それはまた、信号を２乗して、ローパスフィルタリングすることによっても実現することができる。ＢモードイメージングがＩＱデータを用いて使用されると、信号の振幅が抽出される。包絡線検波ステップに続いて、正規化およびダイナミクス圧縮ステップが行われる。ドプライメージングおよびエラストグラフィの場合、再構成されたＲＦまたはＩＱデータは、後処理を伴わずに直接使用される。

提案された方法を詳細に説明する前に、本記載で使用されるいくつかの表記を最初に説明する。本記載で考慮されるパルスエコー構成では、媒体

からのエコーは、

である反射率（関数）γ（ｒ）によって特徴づけられ、センサ１３によって検出される。図２は、ＰＷイメージングの場合の特定の簡略化された測定構成を示している。以下の表記は、本記載全体を通して使用される。
［数学表記または類似表記］
ベクトルは太字で示されている。

は、ベクトル

の
l_pノルムを示し、

は、実数の空間を示す。
●ヒルベルト空間Ｗの２乗可積分関数の空間は、Ｌ_２（Ｗ）と示される。２つの関数

間の内積は、

と定義される。

超音波関連の表記
●ｍ（ｐ，ｔ）は、ｐに配置されたセンサ１３によって時刻ｔに受信された電気信号である。言い換えると、ｍ（ｐ，ｔ）は、測定格子点（ｐ，ｔ）において計算された測定値の連続的な集合の１つのサンプル（素子生データの１つのサンプルと同等）を定義する。
●γ（ｒ）は、反射率格子点ｒにおいて計算された反射率（画像サンプルと同等）の１つの値を定義する。
●Ωは、対象となる媒体または物体、すなわち、反射率サンプルの場所の連続的な集合を示す。
●ｏ（ｒ，ｐ）は、ｐに位置するセンサとｒに位置する点と間の素子指向性である。
●ａ（ｒ，ｐ）は、ｐに位置する点源とｒに位置する点との間の減衰の値を定義する。
●ｖ_pe（ｔ）は、受信したパルス形状（すなわち、パルスエコー波形）である。
●ｔ_Tx（ｒ，ｑ）は、送信伝搬遅延、すなわち、ｑに位置する送信デバイスまたは素子１３によって送信されるパルス波がｒに位置する点に到達するのにかかる時間を示す。
●ｔ_Rx（ｒ，ｐ）は、受信伝搬遅延、すなわち、ｒに位置する点によって送信されるパルスエコー波がｐに位置するセンサに到達するのにかかる時間を示す。

測定格子関連の表記
●送信素子の集合、

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

ここで、Ｎ_elrは、受信センサの数を指し、受信機とも呼ばれる。
●送信素子内の点の集合：

●受信センサ内の点の集合：

●時間サンプルの集合：

反射率格子関連の表記
●反射率の場所の集合：

ここで、Ｎ_γは、反射率の数を指す。

提案された方法は、第１の格子および第２の格子を定義し、共に３Ｄ空間にある。第１の格子は、反射率格子と呼ばれ、それは、２つの交差する格子線の交点における反射率格子点を含み、それは、本例では、互いに対して実質的に９０度の角度を形成する。反射率格子の全ての反射率格子点は、媒体Ωに属する。第２の格子は、測定格子と呼ばれ、それは、２つの交差する格子線の交点における測定格子点を含み、それは、本例では、互いに対して実質的に９０度の角度を形成する。測定格子点は、１次元では受信センサの場所と一致し、他の次元では各センサのサンプリング時刻と一致する。

次に、本発明による所与の測定格子点での測定モデルの推定、評価、計算または実装について説明する。本例では、測定モデルは、各測定格子点で推定される。本開示で使用される測定モデルは、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によって最初に導入されたパルスエコー空間インパルス応答モデルである。本モデルによれば、位置ｐおよび時間ｔで記録される生データサンプルを説明する以下の方程式を書くことができる。

次に、１つの送信点源素子および１つの受信点源センサの構成に関する測定モデルの連続パラメトリック定式化について説明する。

測定モデルは、

と定義される。

上記の方程式は、以下に示すように表現されることができ、測定モデルを超曲面Γ（ｐ、ｑ、ｔ）に沿った積分と定義する。

ここで

上記の方程式では、超曲面、またはより具体的には反射率超曲面は、方程式ｇ（ｒ、ｐ、ｑ、ｔ）＝０によって定義されるΓ（ｐ、ｑ、ｔ）であり、ここで、対象となる変数は、ｒである。Γ（ｐ、ｑ、ｔ）は、ｒが２Ｄ平面（放物線または楕円）にある場合は円錐曲線を定義し、ｒが３Ｄボリュームにある場合は２次曲面を定義する。同じ形状は、測定超曲面にも適用する。図３は、連続領域における測定構成を概略的に例解し、反射率超曲面も示す。

次に、反射率超曲面Γ（ｐ、ｑ、ｔ）は、積分次元を減らすためにパラメータ化することができる。方程式ｇ（ｒ、ｐ、ｑ、ｔ）＝０によって定義される反射率超曲面Γ（ｐ、ｑ、ｔ）は、反射率パラメトリック方程式の集合によって同等に定義されることができる。数学では、パラメトリック方程式は、数量の集合を１つ以上の独立変数（すなわち、パラメータ）の関数と定義する。パラメトリック方程式は、一般に、面または曲線などの幾何学的対象を構成する点の座標を表すために使用され、その場合、方程式は、物体のパラメトリック表現またはパラメータ化と総称される。反射率パラメトリック方程式の集合は、ここで、以下、

と定義されることができる。

ｒは、パラメータα∈Ａ、に依存し、ここで、

ならば

の場合

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

ここで、

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

と定義される。

１つの送信点源素子および１つの受信点源センサの構成に関連して記載したモデルとの主な相違は、任意形状の有限長センサを形成する全ての点源の寄与の和に対応する、もう１つの積分があるという事実にあることがわかる。

対応するパラメトリック定式化は、

によって与えられる。

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

によって示され、ここで、

である。それは、例えば、Πⁱの中心の点の座標であり得る。

それは、

と書くことができる。結果として得られた測定値は、Πⁱによって定義されたトランスデューサ面に沿って積分することにより取得される。ｍ（ｐ，ｔ）を定義する方程式を上記の積分に統合すると、次の関係、

が成り立つ。

対応するパラメトリック定式化は、

によって与えられる。

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

が導かれる。

この方程式は、任意形状の１つの送信有限長素子および任意形状の１つの受信有限長センサの構成に関連して記載した方程式と非常に類似している。唯一の相違は、複数の送信素子を考慮するために導入された和にある。

対応するパラメトリック定式化は、

によって与えられる。

次に、任意形状の複数の送信有限長素子および任意形状の複数の受信有限長センサの構成に関する連続測定モデルについて説明する。

複数の受信センサの集合を使用する場合、点源場所の集合Πⁱおよび座標ξⁱは、受信センサｉに関連付けられている。したがって、受信測定領域は、

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

が成り立つ。

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

を定義するパラメータαのいくつかの離散値α^Ｐを選択することを意味する。パラメータαの離散値の選択が、次にＮ_α点

によってのみ記載される超曲面Γ（ｐ，ｑ，ｔ）の離散化をもたらすことがわかる。図４は、離散領域における測定構成を概略的に例解し、反射率超曲面も示している。測定モデルのパラメトリック定式化の離散化は、

として表現することができ、
ここで、

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

であり、

であり、
ここで、ｗ（ｑ^ｋ），ｕ（α^ｐ）およびｚ（ｐ^ｎ）は、積分の重み（連続積分の離散化に関連する）である。送信素子の各点ｑ^ｋについて、上記式は、反射率超曲面の点上に推定される、超曲面反射率サンプルまたは値と呼ばれる、反射率の集合

を含むことに気づくことができる。

次に、超曲面反射率は、反射率格子点と実質的に一致するように補間される。超曲面反射率は、反射率格子点上にないため、通常は未知である。したがって、それらは、Ｐ．Ｔｈｅｖｅｎａｚ，Ｔ．ＢｌｕａｎｄＭ．Ｕｎｓｅｒ，"Ｉｎｔｅｒｐｏｌａｔｉｏｎｒｅｖｉｓｉｔｅｄ"，ＩＥＥＥＴｒａｎｓａｃｔｉｏｎｓｏｎＭｅｄｉｃａｌＩｍａｇｉｎｇ，ｖｏｌ．１９，ｐ．７３９−７５８，２０００によって導入される補間カーネルによって、反射率格子上で定義された反射率サンプルから近似する必要がある。したがって、

であるような補間カーネル

を導入する。

上記補間関数は、超曲面反射率

（通常、反射格子点上にない）を反射率サンプルγ^Sに関連付ける。

次に、離散化および補間された超曲面反射率サンプルまたは値の総和が実行される。補間関数は、離散化された測定モデルに挿入され、次の総和、

が実行される。

本例では、先の方程式に記載したように、考慮された測定格子点において測定値推定値とも呼ばれる測定値の推定値を取得するために、最終的に総和の出力がパルス形状に畳み込まれる。上記のプロセスは、全ての測定値に対して実行される。したがって、測定モデルの推定または処理の結果は、それぞれの測定格子点における測定値の推定値である。

次に、各反射率格子点における測定モデルの随伴作用素の推定、評価、計算、または実装が実行される。まず、測定モデルの随伴作用素の連続パラメトリック定式化が取得される。随伴作用素の連続定式化は、機能解析ツールを使用して、測定モデルから、すなわち上記で与えられた連続測定モデルのいずれかから導出することができる。以下の方程式は、任意形状の複数の有限長送信素子および任意形状の複数の有限長受信センサの場合から導出される。実際に、

であるような作用素

を定義する場合、Ｈの随伴作用素

は、

と定義される。

以下の内積、

を考慮し、
ここで、

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

を有し、
ここで、

。

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

と定義する。

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

として表現することができる。

次いで、上記の方程式に対応する測定超曲面は、

である。

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

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

を有し、
ここで、ｗ（ｑ^ｋ）およびＺ（ｐ^ｎ）は、積分の重みであり、

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

を含むことに気付くことができる。

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

であるような、補間カーネル

を導入する。

上記の方程式は、超曲面測定サンプル

を畳み込み測定サンプル

に関連付ける。

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

が実行される。

総和の出力は、考慮された反射率格子点における反射率の推定値である。

次に画像再構成手順を説明する。最初に逆問題を定義する。測定値ｍは、上記の測定モデルの随伴作用素を使用して、反射率の推定値

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

を使用して、前記の測定モデルを使用して、ｍに等しくない測定値の推定値

を生成することができる。したがって、測定値の推定値と反射率の推定値と間の関係がわかる。

画像再構成方法は、測定値ｍの集合が与えられた場合、反射率格子点において評価される反射率の推定値

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

によって定義され、
それは、測定格子の各点で推定される測定モデルに対応する。

さらに、作用素

は、γに対して線形であるため（所与の格子点において測定モデルを評価するときに表わされた式から推定することができる）、
ｍ＝Ｈγであるような行列

（線形作用素Ｈに関連付けられている）が存在し、それは、
逆問題を定義する。

次に、逆問題に関連付けられる最適化問題について説明する。上記の逆問題を解くために、それは、以下の最適化問題、

として表現され、
ここで、Ｊ（γ，ｍ）は、最適化問題に関与する目的関数を示し、Ｆ（Ｈγ，ｍ）は、データ不一致項を考慮する下半連続汎関数であり、測定値の推定値Ｈγと測定値ｍとの間の距離を測定し、Ｒ（γ）は、任意選択の事前分布項を記述する下半連続汎関数であり、それは、反射率に関する特定の統計的挙動など追加情報を考慮し、λ＞０は、正則化パラメータである。

データ不一致項に使用される汎関数のいくつかの例を以下に列挙する。
●ユークリッド距離、

、
●半径

の

球体に関する指標関数、

、
●

ノルム、

。

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

ノルムのｐ乗、

、
●ｐ＜１の場合、それは、γのスパース性の尺度である。
●ｐ∈［１，２］の場合、それは、データが一般化ガウス分布（ＧＧＤ）にどれだけ適合するかについての尺度である。
●所与のモデル

における

ノルムのｐ乗：

●ψは、ウェーブレット変換またはフーリエ変換など一般的な変換であり得る。
●ψは、学習された辞書であり得る。

次に、最適化問題を１つの画像事前分布でどのように解くことができるかについて説明する。ＦおよびＲの特性に応じて、異なる種類のアルゴリズムを使用し得る。

次に、ＦおよびＲが微分可能である場合について説明する。この場合、逆問題は、γに関する目的関数Ｊ（γ，ｍ）の導関数の根を見つけることによって解かれ、

と書くことができる。

問題がＦ（Ｈγ，ｍ）の導関数の計算にＨの計算を含むことがわかる。関数Ｆ（Ｈγ，Ｍ）に応じて、導関数の計算はまた、測定モデルＨの随伴作用素の計算を含み得る。かかる問題の１つの有名な例は、Ｔｉｋｈｏｎｏｖ正則化であり、ここで、

である。この場合、以下の解

を有し、ここで、

は、単位行列、つまり、主対角線上は１、他の場所はゼロである、Ｎ_γ×Ｎ_γ正方行列である。

次に、Ｆのみが微分可能であり、Ｒが凸である場合について説明する。かかる問題を解くために使用される方法の１つの一般的なグループは、射影勾配法と呼ばれる。射影勾配法は、最適化問題の解が固定点方程式を満たすという事実を利用する。１つの一般的なアルゴリズムは、以下の固定点方程式のシステムを含む前方−後方分離であり、本具体例では、１つの方程式に限定されている。

ここで、

は、Ｒに関連付けられた近接作用素を示し、近接作用素は、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によって導入されている。上記の方程式は、方程式の両側に

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

を実行することによって問題を解く。

別の一般的なアルゴリズムは、高速反復縮小閾値アルゴリズム（ＦＩＳＴＡ）であり、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によって導入され、それは、上記の方法の変形である。

前方−後方分離の１つの例は、

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

、
ここで、

は、ハイパーパラメータである。

実際、

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

と定義され、ここで、ソフト閾値作用素は、素子ごとにｓｉｇｎ（ｘ_ｉ）ｍａｘ（｜ｘ_ｉ｜−τ、０）と定義され、ここで、ｓｉｇｎ（ｘ_ｉ）は、ｘ_ｉ＞０の場合は１、ｘ_ｉ＝０の場合は０、それ以外の場合は−１である。

次に、ＦおよびＲが微分可能でない場合について説明する。この場合、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によって導入されたプライマルデュアルフォワードバックワード（ＰＤＦＢ）アルゴリズムなど、いくつかの方法を使用することができる。

次に、最適化問題をＭ個の画像事前分布でどのように解くことができるかについて説明する。この場合、事前分布項は、

と、書くことができると考えられることができる。

この場合、プライマルデュアルフォワードバックワード、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によって導入された同時方向乗数法（ＳＤＭＭ）など、いくつかのアルゴリズムを使用することができる。これらのアルゴリズムの共通の特性は、Ｆ（Ｈγ，ｍ）（またはＦ（Ｈγ，ｍ）の導関数）およびＭ個の画像事前分布項の近接作用素の計算を必要とすることである。Ｆ（Ｈγ，ｍ）の導関数または近接作用素の計算は、関数Ｆ（Ｈγ，ｍ）に応じて、測定モデルおよび測定の随伴作用素の推定を含み得る。

上記の画像再構成方法は、図５〜９のフローチャートによって要約することができる。図５のフローチャートは、方法の概要を示す。ステップ２１では、パルス波の集合が、送信素子の集合によって送信される。本例では、送信素子は、有限のサイズおよび任意の形状を有する。送信されたパルスの集合は、１つのパルスのみ、または連続して送信される２つ以上のパルスを含み得る。ステップ２３では、反射されたエコー波形の集合が、本例では有限のサイズおよび任意の形状を有するセンサの集合によって受信される。ステップ２５では、測定値の集合を反射率の集合に関連付ける逆問題が定義される。ステップ２７では、逆問題の解が、データ不一致項および１つ以上の画像事前分布項を含む目的関数を含む最適化問題として表現される。ステップ２９では、最適化問題が、射影を使用して解かれ、反射率の推定値を取得する。

図６のフローチャートは、汎関数Ｆが微分可能である場合の提案された方法の概要をより詳しく示す。最初の４つのステップ、すなわちステップ３１、３３、３５、および３７は、図６のフローチャートの最初の４つのステップと同じである。ステップ３９では、第１の反射率推定値

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

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

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

の場合、Ｆ（Ｈγ，Ｍ）の導関数は、

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

が更新される。ステップ５０では、収束基準が満たされているか否かが判定される。肯定の場合、プロセスは、終了する。基準が満たされていない場合、プロセスは、ステップ４１に続く。

図７のフローチャートは、データ不一致項がユークリッド距離であり、逆問題が１つの画像事前分布項を含む場合の提案された方法の概要を、より詳細に示す。最初の５つのステップ、すなわちステップ５１、５３、５５、５７、５９は、図６のフローチャートの最初の５つのステップと実質的に同じであるが、ここでは、逆問題が１つの画像事前分布項を含むという相違を有する。ステップ６１では、変数

が、図８のフローチャートに従って、各測定格子点に対して計算される。ステップ６３では、残差ｒ^ｎ＝ｖ^ｎ−ｍが計算される。ステップ６５では、残差

についての測定モデルの随伴作用素が図１０のフローチャートに従って推定される。ステップ６７では、値

が計算される。ステップ６９では、画像事前分布項への射影

が計算される。ステップ７０では、収束基準が満たされているかどうかが判定される。肯定の場合、プロセスは、終了する。基準が満たされてない場合、プロセスは、ステップ６１に続く。

図８のフローチャートは、測定モデルを推定または処理するプロセスの例を説明する。ステップ７１では、反射率パラメトリック方程式の集合が、各測定格子点（任意形状の有限サイズ受信センサを定義する）に対して、測定格子点に対応する任意形状の有限サイズ受信機内の各点源場所に対して、任意形状の各有限サイズ送信機に対して、および任意形状の有限サイズ送信機内の各点源場所に対して、生成される。言い換えると、反射率パラメトリック方程式の集合は、各送信機および各受信機内の点源の全てのペア間の各送信機−受信機波経路に対して生成される。反射率パラメトリック方程式は、反射率超曲面の集合を定義する。ステップ７３では、反射率超曲面が離散化されて、超曲面反射率サンプルの集合が取得される。ステップ７５では、超曲面反射率サンプルが反射率格子上に補間される。ステップ７７では、補間され、および離散化された超曲面反射率サンプルが合計される。ステップ７９では、ステップ７７で取得された値が、測定格子の点における測定値推定値に追加される。ステップ８０では、取得された累積測定値推定値がパルス形状に畳み込まれる。

図９のフローチャートは、測定モデルの随伴作用素を推定する、または処理するプロセスの例を説明する。ステップ８１では、測定値推定値がパルス形状の整合フィルタに畳み込まれる。ステップ８３では、測定パラメトリック方程式の集合が、反射率格子の各点に対して、任意形状の各有限サイズ送信機に対して、任意形状の各有限サイズ受信機に対して、任意形状の有限サイズ送信機内の各点源場所に対して、および任意形状の有限サイズ受信機に内の各点源場所に対して、生成される。言い換えると、測定パラメトリック方程式の集合は、各送信機および各受信機内の点源の全てのペア間の各送信機−受信機経路に対して生成される。測定パラメトリック方程式は、測定超曲面の集合を定義する。ステップ８５では、測定超曲面が離散化されて、超曲面測定サンプルの集合が取得される。ステップ８７において、超曲面測定サンプルが測定格子上に補間される。ステップ８９では、補間され、および離散化された超曲面測定サンプルが合計される。ステップ９０では、ステップ８９で取得された値が、反射率格子の点における反射率推定値に追加される。

送信されたパルス波は、電気信号によって励起される少なくとも１つの電気機械変換装置によって生成されることができる。測定値は、少なくとも１つの相互電気機械変換装置によって生成される電気信号から取得することができる。電気機械変換装置は、実質的に１波長ピッチを有する多素子構造を備え、直線、凸状線または凹状線に沿って整列された１つ以上のリニアアレイプローブ、実質的に半波長ピッチを有する多素子構造を備え、直線、凸状線または凹状線に沿って整列された１つ以上のフェーズドアレイプローブ、実質的に１波長ピッチを有する多素子構造を備え、平面上に整列された１つ以上の行列アレイプローブ、または実質的に半波長ピッチを有する多素子構造を備え、平面上に整列された１つ以上の行列アレイプローブ、に空間的に配置される。

本発明は、図面および前述の記載で詳細に例解および記載されたが、かかる例解および記載は、事例的または例示的であり、限定的ではないとみなされるべきであり、本発明は、開示された実施形態に限定されない。他の実施形態および変形形態は、図面、開示および添付の特許請求の範囲の研究に基づいて、請求された発明を実施する際に当業者によって、理解され、達成されることができる。

例えば、本発明の教示は、送信音響波が電磁放射パルスで置き換えられる、光音響イメージングに適用することができる。この特定の場合、電磁放射は、光の速度に近い速度、すなわち通常、音の速度よりも５桁速い速度で伝搬するため、全ての音響送信伝搬遅延をゼロに設定することで、本記載で導入される形式全体を簡素化することができる。このため、受信時の音響伝搬遅延を考慮する場合、送信伝搬遅延は無視することができる。本発明で考慮される別の変形は、容量性マシン超音波トランスデューサ（ＣＭＵＴ）など、圧電材料とは異なるトランスデューサ技術の使用にある。さらに別の変形は、本発明を時間多重方式に適用することであり、それにより送信および／または受信トランスデューサ素子のサブセットが連続して順次アドレス指定され、結果として得られる信号が後続の再結合のために記録される。上記の用途は全て、パルスインバージョン、振幅変調、振幅変調パルスインバージョン、ハーモニックイメージングなど、当技術分野で知られているマイクロバブルベースのコントラスト固有の非線形イメージング技術に変換することもできる。

本発明はまた、イメージング装置１など、電子データ処理装置のコンピューティング手段上でロードおよび実行される場合、上記で説明した方法のステップを実装するための非一時的媒体に記憶された命令を含むコンピュータプログラム製品を提案する。

特許請求の範囲において、「含む」という語は、他の素子またはステップを除外せず、不定冠詞「ａ」または「ａｎ」は複数を除外しない。異なる特徴が相互に異なる従属請求項に記載されているという単なる事実は、これらの特徴の組み合わせを有利に使用することができないことを示していない。

Claims

少なくとも部分的に画像を再構成すべき物体に対し、送信素子を用いてパルス波を送信するステップと、
反射率格子点を含む反射率格子上に定義された反射率を用いて特徴付けられた前記物体から反射され、測定格子点を含む測定格子上に定義された測定値を定義するエコー波形を、センサを用いて受信するステップと、
測定モデルを用いて前記測定値から反射率推定値を取得するための逆問題を定義するステップと、
前記逆問題を、データ不一致項を有する目的関数を含む最適化問題として表現するステップと、
前記最適化問題を、前記データ不一致項の寄与を含む連立方程式として表現するステップと、
前記連立方程式において所定の収束閾値に到達するまで反復を実行することにより、前記反射率の推定値を取得する取得ステップと、
を含むパルス波画像再構成方法であって、
前記反復のうちの少なくともいくつかは、前記測定格子点のそれぞれに対応する前記測定値の推定、および前記反射率格子点のそれぞれに対応する前記反射率の推定を少なくとも含み、
前記取得ステップは、
超曲面反射率を含む反射率超曲面を定義する、前記送信素子とセンサとの間の経路ごとに、反射率パラメトリック方程式を生成し、前記超曲面反射率を、前記反射率超曲面のそれぞれにおいて積分して、前記測定モデルから、前記測定格子点のそれぞれに対応する前記測定値の推定値を取得するステップと、
超曲面測定値を含む測定超曲面を定義する、前記送信素子とセンサとの間の経路ごとに、測定パラメトリック方程式を生成し、前記超曲面測定値を、前記測定超曲面のそれぞれにおいて積分して、前記測定モデルの随伴作用素から、前記反射率格子点のそれぞれに対応する前記反射率の推定値を取得するステップと、
を含むパルス波画像再構成方法。
前記送信素子のそれぞれは、送信点とみなされ、前記センサのそれぞれは、受信点とみなされる、請求項１に記載のパルス波画像再構成方法。
前記超曲面反射率を積分する前に、
前記反射率超曲面のそれぞれを離散化して前記超曲面反射率を導きだし、前記反射率格子点と略一致するように前記超曲面反射率を補間する請求項１または２に記載のパルス波画像再構成方法。
前記超曲面測定値を積分する前に、前記測定超曲面のそれぞれを離散化して前記超曲面測定値を導き出し、前記測定格子点と実質的に一致するように前記超曲面測定値を補間する請求項１〜３のいずれか１項に記載のパルス波画像再構成方法。
前記データ不一致項は、前記測定格子点のそれぞれについて前記測定値と前記測定値の推定値との間の距離を評価する正の汎関数として表わされる、請求項１〜４のいずれか１項に記載のパルス波画像再構成方法。
前記正の汎関数は、微分可能な正の汎関数である、請求項５に記載のパルス波画像再構成方法。
前記目的関数は、前記データ不一致項と少なくとも１つの画像事前分布項との組み合わせとして表わされ、前記画像事前分布項は、前記反射率に関する事前情報を含む、請求項１〜６のいずれか１項に記載のパルス波画像再構成方法。
前記少なくとも１つの画像事前分布項は、所定のモデルにおいて前記反射率推定値をｐとした場合に、l_pのノルムのｐ乗として表わされ、前記l_pのノルムは、ベクトル

に対して

として定義され、

は、実数の空間を表す、請求項７に記載のパルス波画像再構成方法。
前記パルス波は、パルス音響波または電磁放射波である、請求項１〜８のいずれか１項に記載のパルス波画像再構成方法。
前記パルス波は、偏向平面波または発散波を含む、請求項１〜９のいずれか１項に記載のパルス波画像再構成方法。
前記パルス波は、電気信号によって励起される少なくとも１つの電気機械変換装置によって生成され、前記測定値は、少なくとも１つの相互電気機械変換装置によって生成される電気信号から取得される、請求項１〜１０のいずれか１項に記載のパルス波画像再構成方法。
前記電気機械変換装置は、
略１波長ピッチを有する多素子構造を備え、直線、凸状線または凹状線に沿って整列されたリニアアレイプローブ、
略半波長ピッチを有する多素子構造を備え、直線、凸状線または凹状線に沿って整列されたフェーズドアレイプローブ、
略１波長ピッチを有する多素子構造を備え、平面上に整列されたマトリクスアレイプローブ、または、
略半波長ピッチを有する多素子構造を備え、平面上に整列されたマトリクスアレイプローブ、に空間的に配置される請求項１１に記載のパルス波画像再構成方法。
前記取得ステップにおいて、
前記超曲面測定値を積分する際、または前記超曲面反射率を積分する際には、数値積分を行ない、
前記測定格子点のそれぞれについて、取得された前記測定値の推定値を、先に取得された前記測定値の推定値に加算して、測定値推定値の累積値を取得し、
前記累積値を、送信された前記パルス波のパルス形状と畳み込む、請求項１〜１２のいずれか１項に記載のパルス波画像再構成方法。
物体の画像を再構成する画像処理装置であって、
少なくとも部分的に再構成すべき物体に対し、送信素子を用いてパルス波を送信する送信部と、
反射率格子点を含む反射率格子上に定義された反射率を用いて特徴付けられた前記物体から反射され、測定格子点を含む測定格子上に定義された測定値を定義するエコー波形を、センサを用いて受信する受信部と、
測定モデルを用いて前記測定値から反射率推定値を取得するための逆問題を定義する定義部と、
前記逆問題を、データ不一致項を有する目的関数を含む最適化問題として表現する第１表現部と、
前記最適化問題を、前記データ不一致項の寄与を含む連立方程式として表現する第２表現部と、
前記連立方程式において所定の収束閾値に到達するまで反復を実行することにより、前記反射率の推定値を取得する取得部と、
を備え、
前記反復のうちの少なくともいくつかは、前記測定格子点のそれぞれに対応する前記測定値の推定、および前記反射率格子点のそれぞれに対応する前記反射率の推定を少なくとも含み、
前記取得部は、
超曲面反射率を含む反射率超曲面を定義する、前記送信素子とセンサとの間の経路ごとに、反射率パラメトリック方程式を生成し、前記超曲面反射率を、前記反射率超曲面のそれぞれにおいて積分して、前記測定モデルから、前記測定格子点のそれぞれに対応する前記測定値の推定値を取得する第１取得部と、
超曲面測定値を含む測定超曲面を定義する、前記送信素子とセンサとの間の経路ごとに、測定パラメトリック方程式を生成し、前記超曲面測定値を、前記測定超曲面のそれぞれにおいて積分して、前記測定モデルの随伴作用素から、前記反射率格子点のそれぞれに対応する前記反射率の推定値を取得する第２取得部と、
を備えた画像処理装置。