JP2022503946A - 学習非線形マッピングに基づく画像再構成方法 - Google Patents

学習非線形マッピングに基づく画像再構成方法 Download PDF

Info

Publication number
JP2022503946A
JP2022503946A JP2021517765A JP2021517765A JP2022503946A JP 2022503946 A JP2022503946 A JP 2022503946A JP 2021517765 A JP2021517765 A JP 2021517765A JP 2021517765 A JP2021517765 A JP 2021517765A JP 2022503946 A JP2022503946 A JP 2022503946A
Authority
JP
Japan
Prior art keywords
image
functional
image reconstruction
reconstruction method
regularized
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
JP2021517765A
Other languages
English (en)
Other versions
JP7345208B2 (ja
Inventor
ディミトリス・ペルディオス
エイドリアン・ベッソン
フロリアン・マルティネス
マルセル・アルディティ
ジャン=フィリップ・ティラン
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Ecole Polytechnique Federale de Lausanne EPFL
Original Assignee
Ecole Polytechnique Federale de Lausanne EPFL
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Ecole Polytechnique Federale de Lausanne EPFL filed Critical Ecole Polytechnique Federale de Lausanne EPFL
Publication of JP2022503946A publication Critical patent/JP2022503946A/ja
Application granted granted Critical
Publication of JP7345208B2 publication Critical patent/JP7345208B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/04Analysing solids
    • G01N29/043Analysing solids in the interior, e.g. by shear waves
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/04Analysing solids
    • G01N29/06Visualisation of the interior, e.g. acoustic microscopy
    • G01N29/0654Imaging
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/04Analysing solids
    • G01N29/11Analysing solids by measuring attenuation of acoustic waves
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/22Details, e.g. general constructional or apparatus details
    • G01N29/32Arrangements for suppressing undesired influences, e.g. temperature or pressure variations, compensating for signal noise
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/44Processing the detected response signal, e.g. electronic circuits specially adapted therefor
    • G01N29/4481Neural networks
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/88Sonar systems specially adapted for specific applications
    • G01S15/89Sonar systems specially adapted for specific applications for mapping or imaging
    • G01S15/8906Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques
    • G01S15/8977Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques using special techniques for image reconstruction, e.g. FFT, geometrical transformations, spatial deconvolution, time deconvolution
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/52017Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging
    • G01S7/52046Techniques for image enhancement involving transmitter or receiver
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06NCOMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00Computing arrangements based on biological models
    • G06N3/02Neural networks
    • G06N3/08Learning methods
    • G06N3/084Backpropagation, e.g. using gradient descent
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/005Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/002D [Two Dimensional] image generation
    • G06T11/003Reconstruction from projections, e.g. tomography
    • G06T11/006Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00Indexing codes associated with group G01N29/00
    • G01N2291/02Indexing codes associated with the analysed material
    • G01N2291/023Solids
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00Indexing codes associated with group G01N29/00
    • G01N2291/04Wave modes and trajectories
    • G01N2291/044Internal reflections (echoes), e.g. on walls or defects
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00Indexing codes associated with group G01N29/00
    • G01N2291/10Number of transducers
    • G01N2291/106Number of transducers one or more transducer arrays
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/52017Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging
    • G01S7/52077Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging with means for elimination of unwanted signals, e.g. noise or interference
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00Image generation
    • G06T2211/40Computed tomography
    • G06T2211/424Iterative

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Pathology (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • Immunology (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Artificial Intelligence (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Physics (AREA)
  • Signal Processing (AREA)
  • Molecular Biology (AREA)
  • Software Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Biophysics (AREA)
  • Computing Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Biomedical Technology (AREA)
  • Computational Linguistics (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

Figure 2022503946000001
対象物から測定値yとして1セットの波形を受信する受信ステップ(13)と、対象物の未知画像xと測定値yとを結びつける測定モデルH(x)を定義する定義ステップ(15)と、測定モデルH(x)と測定値yとの間の距離を測定するデータ忠実度汎函数D(x)を決定する決定ステップ(17)と、入力データに関して非線形な、学習非線形マッピングfθを用いて、未知画像に関する事前情報を含む正則化汎函数R(x)を決定する第2決定ステップ(19)と、未知画像の第1画像推定値x^を得るために、データ忠実度汎函数Dと正則化汎函数Rとを含む最適化問題を定義する定義ステップ(21)と、最適化問題を解いて、第1画像推定値x^を算出する算出ステップ(23)と、を備えた画像再構成方法。
【選択図】図2

Description

本発明は、画像再構成方法に関するものである。より具体的には、本発明は、学習非線形マッピングを用いて画像を再構成する。本発明は、例えば、超音波撮像において使用することができる。また、本発明は、対応する画像処理装置、および本方法を実施するためのコンピュータプログラム製品に関するものである。
例えばレーダー、ソナー、地震画像、超音波(US)画像などに用いられるアレイ処理技術は、さまざまな技術分野で広く利用されている。特に医療分野では、腱、筋肉、関節、血管、内臓などの体内構造を見るために使用されているが、超音波イメージングは、医療分野以外でも広く利用されている。超音波イメージングの商業的利用の1つに、非破壊検査(NDT)がある。NDTとは、部品やシステムを破壊することなく、材料、コンポーネント、アセンブリの不連続性や特性の違いを検査、試験、評価するプロセスのことである。言い換えれば、NDTによれば検査や試験が完了しても、その部品を使用することができる。NDTの例としては、石油パイプの超音波テストが挙げられる。この場合、例えばパイプ内において超音波のパルス波を伝播させ、パルス波の反射に基づいて、亀裂や欠陥の有無を検出することができる。管内の音響波の速度と、亀裂から反射してパルス波を送信するまでの時間間隔が分かれば、欠陥の位置を推定することができる。このパイプの亀裂検出と同様の理論を用いて、例えば航空機の構造体/翼のパネルなどの亀裂を検出したり、翼の表面に氷が溜まっているのを検出したりすることができる。
超音波イメージングは、その高い柔軟性、コスト効率、安全性、非侵襲性により、イメージングモダリティとしても広く利用されている。パルスエコー法による超音波イメージングは、従来から、電気信号を音響圧力波に変換するトランスデューサエレメントのアレイを用いて行われる。このようなトランスデューサエレメントのアレイは、対象となる媒体を伝搬する短い音響波を送信するために使用される。音響インピーダンスの局所的な変化から生じる後方散乱エコーは、同じアレイで受信され、いわゆる生のエレメント高周波(RF)データとして検出される。高周波画像は、一般的なDAS(delay-and-sum)アルゴリズム(遅延加算アルゴリズム)を用いて、これらの測定値からリアルタイムに再構成される。この高周波画像に包絡線検波とログ圧縮を施して得られるBモード画像は、表示や保存の目的で一般的に使用される。
Clarke F., (2013) Generalized gradients, "Functional Analysis, Caliculus of Variations and Optimal Control", Graduate Texts in Mathematics, vol. 264. Springer, London.
近年、いわゆる超高速超音波イメージングが研究者の間で注目されている。集束した送信ビームを用いて複数回のパルスエコー測定を行う従来の超音波イメージングに比べて、超高速超音波イメージングでは、平面波(PW)や発散波(DW)などの非集束の波面を送信して、視野全体に対して同時に(一度に)音波を当てる("insonify"する)というアイデアを利用している。これにより、非集束波面とそれに伴う後方散乱エコーの往復時間と画像再構成時間のみの制限で、数kHzという極めて高いフレームレートを実現することができる。このようなフレームレートにより、エラストグラフィ、超高速ベクトルフローイメージング、ファンクショナル超音波イメージングなどといった、新たな超音波イメージング診断法を実現できた。超高速超音波イメージングの主な欠点は、連続した集束ビームを使用した場合に比べてSNR(信号対雑音比)が低下し、画質が劣化することである。実際、エネルギーが特定の領域に集中する集束ビーム送信時と比較して、非集束波面のエネルギーは視野全体に広がってしまう。
画質を全体的に向上させる一般的な方法として、様々な方向から音波を当てることにより再構成された複数の画像を平均化する「コヒーレント合成」と呼ばれるプロセスが知られている。この方法は、ある意味、コンピュータ断層撮影(CT)のように複数の投影画像から画像を再構成するのと似ている。このような技術を実現するのは簡単だが、複数の送受信イベントや画像再構成プロセスが必要になるため、フレームレートの低下、データ転送量の増加、計算量の増加といった代償を払わなければならない。
超高速超音波イメージングにおいて、高品質な画像を再構成するという課題は、超音波イメージングのコミュニティにおいて重要な関心事となっている。最近のアプローチでは、焦点の合わない透過波面の数を最小限に抑えて、高品質な画像を再構成しようとしている。既存のアプローチの中には、少ない送受信イベントで高品質な画像を実現するものもあるが、計算が非常に複雑であるため、実時間での使用は困難である。
すなわち、基地のアレイ処理による画像再構成法は、高品質な画像を低速度で提供したり、低品質な画像を超高速で提供したりするなど、最適とはいえず、妥協が必要である。
本発明の目的は、アレイ処理の分野において、画像の再構成に関連する上記の問題の少なくとも一部を克服することである。
本発明の第1の側面によれば、請求項1に記載された画像再構成方法が提供される。
提案された新しいソリューションは、画像を非常に高速に再構成でき、しかも非常に高い品質を実現できるという利点がある。
そこで本発明では、効果の高いデータ駆動型の事前情報を用いて最適化問題を解くことにより、逆問題として定式化することができる画像再構成法を提案する。このような事前情報は、イメージングタスクの目的のために特別に学習された人工ニューラルネットワーク(ANN)などの非線形マッピングを用いて構築される。データ駆動型事前分布の特性に応じて、勾配降下法、プラグアンドプレイ事前分布、投影勾配降下法などの反復アルゴリズムに使用することができる。ここでは、単一のPWまたはDWから再構成された低品質の画像から、合成開口(SA)測定から再構成されたアーチファクトのない高品質の画像への非線形マッピングを学習することが提案されており、基準またはゴールドスタンダードの方法として考えられている。既存のいくつかのソリューションとは対照的に、本アプローチでは、高品質の高周波画像とスペックルパターンを、そのスペクトル内容を維持したまま再構成することに焦点を当てている。提案されたソリューションは、ニューラルネットワーク(NN)に依存している。ニューラルネットワークは、同様のタスクにおいて非常に効率的で堅牢であることが証明されている。ニューラルネットワークを適切に学習させて低品質画像と高品質画像を対応させることで、高品質画像の一般的な構造に関する非常に具体的な事前情報を得ることができる。さらに、このようなマッピングを実現するために必要な計算量は、リアルタイムイメージングと互換性がある。
本発明の第2の側面によれば、特許請求の範囲に記載された方法を実行するように画像処理装置が提供される。
本発明の他の側面は、本明細書に添付された従属請求項に記載されている。
本発明の一実施例としての画像処理システムの構成を示すブロック図である。 本発明の一実施例としての画像再構成法の処理の流れを示すフローチャートである。 本発明の一実施例によるダイレクトノンリニアマッピングを示すブロック図である。 本発明の一実施例による残差非線形マッピングを示すブロック図である。 本発明の一実施例による「エンコーダ-デコーダ」アーキテクチャを有する残差非線形マッピングを示すブロック図である。 本発明の一実施例による、勾配降下アルゴリズムを用いた画像再構成法の提案を示すフローチャートである。 本発明の一実施例による一般的な投影法ベースのアルゴリズムを用いた画像再構成法の提案を示すフローチャートである。 本発明で使用することができる深層残差畳み込みニューラルネットワークの例を示すブロック図である。
本発明の一実施形態を、添付の図を参照して詳細に説明する。この実施形態では、パルスエコー超音波イメージングについて説明されているが、本発明はこれに限定されるものではない。異なる図面に現れる同一または対応する機能および構造要素には、同じ符号を付している。また、本明細書では、ベクトル変数およびベクトルを出力する関数を太字で表記している。
図1は、本発明の一実施形態に係るイメージングシステム(またはイメージング装置)1の構成要素を示す概略ブロック図である。図1に示されているように、本実施形態のイメージングシステム1は、超音波プローブ3を備えており、送受信モジュール5によって対象媒体中に送信された音波ビームの後方散乱エコーを受信することができる。画像は、画像再構成モジュール7を用いて再構成される。この画像再構成モジュール7は、画像再構成アルゴリズムまたはプロセスを実行する際、学習非線形マッピングまたは学習非線形関数を組み合わせて用いる。学習非線形マッピングは、その入力データ(この例では画像)に関して非線形であり、入力ベクトル空間を出力ベクトル空間にマッピングするように配置されている。後処理・表示モジュール9は、再構成された画像を後処理するように配置されている。この画像再構成法によれば、従来、同様の画質を得るために必要だった送信-受信イベントの数と比較して、送信-受信イベントの数が劇的に減少するので、一般的なアレイ処理アーチファクトのない高品質の画像を再構成することができる。
超音波プローブ3は定形の開口部を有し、そこにトランスデューサとも呼ばれる圧電素子のリニアアレイを有している。このような超音波プローブは、2Dのイメージングに適しており、音響エネルギーを画像平面の近傍に集中させるための円筒フォーカス機能を有していてもよい。しかし、この種のプローブは単なる一例に過ぎない。本実施形態の方法は、様々なプローブ形状(凸型プローブ、フェーズドアレイなど)や、様々な技術(ポリフッ化ビニリデン(PVDF)共重合体や容量性マイクロマシン超音波トランスデューサ(CMUT)など)にも容易に適応できる。同様に、本実施形態の方法は、2Dマトリックスプローブを使用する場合にも拡張することができる。この超音波プローブは、ボリュームイメージングを行うために、音波ビームを送信し、音波を受けた媒体からエコーを収集するように設計されており、イメージングプロセスによって、対象となるボリュームの3D表現が生成される。2Dイメージングでは、平面、凸面、凹面の直線に沿って振動子が配置されたプローブが使用されている。3Dイメージングでは、平面、凸面、凹面に振動子を配置したマトリクスアレイ型のプローブが一般的に使用されている。これらの振動子素子は、電気的な加振を音圧波形に変換したり、音圧波形を電気信号に記録したりすることができる。この特性は、対象となる媒体中での音波ビームの送信と、局所的な音響インピーダンスの変化から生じる対応する後方散乱エコーの受信の両方に利用される。したがって、超音波プローブ3は、空間的に分離された、センサまたはアクチュエータとして機能することができる一組のトランスデューサを備えたアレイ構造を有するアレイ処理構成を形成する。
送受信モジュール5は、送信経路と受信経路とを備え、超音波プローブ3に接続されている。送信経路には、一般的にフロントエンドエレクトロニクスが含まれており、制御されたタイミング、電圧振幅および波形を整形する機能を有する電気的な励起回路を備えてもよい。これは、所望の送信ビームを制御および形成するために使用される。受信経路は、物体(組織など)での減衰を消すためのプログラム可能なゲインコントロール(固定ゲインや時間ゲイン補正など)を備えた電気回路と、アナログ/デジタル変換部、およびデジタル信号処理部を含む。受信した信号または測定値は、画像再構成モジュール7に送られる。
画像再構成モジュール7の目的は、受信した信号から精査中の媒体の画像を再構成することである。この目的のために、画像再構成モジュール7は、fθとして定義される学習非線形マッピングを用いた画像再構成方法、プロセス、またはアルゴリズムを実装する。学習非線形マッピングは、その学習可能なパラメータθが、一般的な超音波撮像アーチファクトを低減または除去するように最適化された人工ニューラルネットワークであってもよい。言い換えれば、非線形マッピングは、学習可能なパラメータを有する機械学習構成であってもよい。これらの学習可能なパラメータの学習は、専用のデータセットを用いた最適化アルゴリズムを用いて行われる。学習が終わると、画像再構成法で学習非線形マッピングが使用され、最小数の送受信イベントで高品質の画像を生成する。
画像再構成法は、次のような形式の逆問題を解くものである。
Figure 2022503946000002
ここで、x∈R、つまり未知画像xは実数RのN乗であり、Nは画像を構成する画素数、すなわちピクセル数またはボクセル数である。また、y∈R、つまり、測定値yは実数RのM乗であり、Mは測定サンプル数である。H∈R→R、つまり、測定モデルHは、物理的な測定プロセスに対応するものである(超音波場合、媒体中の音速によって決まる、波形の到達時間と移動距離の関係を利用することができる)、n∈R、つまり測定ノイズnは、実数RのM乗である。以下の形式の最適化問題を解く反復アルゴリズムを用いて未知画像xまたはその推定値x^を見つけることができる。
Figure 2022503946000003
ここで、D∈R→Rであり、データ忠実度汎函数(データ不一致汎函数としても知られている)Dは、実数RのN乗を正の実数Rり、R∈RN→R(Rは正の実数を表す)は、一般に未知画像xに関する事前情報を符号化する"プライア"と呼ばれる正則化汎函数Rであり、λ∈Rは、データ忠実度汎函数に対する正則化汎函数の重み付けを制御するハイパーパラメータである。汎函数とは、ベクトル空間から実数Rへの(例えば線形)マッピングを指すと理解される。超音波イメージングではあまり魅力的でない、固定ベースまたはフレームにおけるスパース性などの標準的な画像処理プライアを使用する代わりに、学習非線形マッピングfθを使用して、非常に効率的で代表的なデータ駆動型正則化汎函数Rを構築することができる。式(2)を解くために、関数Dと関数Rの両方の数学的特性に応じて、勾配降下法、近位勾配降下法、投影勾配降下法など、多くのよく知られた反復アルゴリズムが文献に掲載されている。このようなアルゴリズムは、正則化汎函数Rの微分可能性から恩恵を受けることがあり、これは非線形マッピングfθが微分可能であることを意味する。このような設定では、学習ニューラルネットワークを非線形マッピングfθとして使用することは、微分可能に構築されるか、またはその微分を正確に近似するための標準的な技術を利用することができるので、有益である。さらに、その微分を効率的に評価することで、アルゴリズムの各反復に関わる計算を高速化することができる。
図2のフローチャートは、一例による提案された画像再構成方法をまとめたものである。ステップ11では、プローブ3は、一組の送信素子を用いて、少なくとも部分的に再構成されるべき媒体中の一組の波形を送信する。しかし、対象物が単に物理的な励起にさらされる可能性があるため、ステップ11は行わなくてもよい。物理的な励起には、超音波(音響波形)、熱、または光音響イメージングにおける電磁放射線パルスなどの他の励起が含まれる。ステップ13では、プローブ3、より具体的にはそのセンサ(複数可)が、一組の受信素子を用いて、音波が当てられた媒体の物理的散乱から生じる一組の(エコー)波形を受信または取得する。言い換えれば、このステップでは、物理的な励起によって刺激された(または暴露された)対象物からの測定値として、一連の波形(または信号)が受信または取得される。ステップ15では、記録された一連の測定値に媒体の未知画像をリンクする測定モデルが定義、決定、または識別される。ステップ17では、定義された測定モデルを用いて、データ忠実度汎函数Dを構築、決定、または識別する。ステップ19では、学習非線形マッピングfθを用いて、データ駆動型の正則化汎函数Rまたはプライアを構築、決定、または識別する。ステップ21では、未知画像の意味のあるまたは高品質の推定値x^を得るために、データ忠実度汎函数Dおよび正則化汎函数Rを含む最適化問題を定義、決定、または識別する。ステップ23では、最適化問題を、例えば、反復アルゴリズムまたはソルバーを用いて、収束基準に達するまで解き、精査対象の媒体の画像推定値x^を得る。アルゴリズムは、10~100回の反復、より具体的には10~50回の反復、または1~20回の反復で構成されていてもよい。ステップ15~21は、画像再構成モジュール7によって実行されてもよい。
画像再構成後、再構成されたデータ(すなわち画像)は、図1にさらに示すように、後処理・表示モジュール9に送られる。後処理ステップは、以下のような様々なアプリケーションをカバーしている。
Bモードイメージング
スキャン変換
カラードップラー画像
ベクトルフローイメージング
エラストグラフィ
高周波データを用いたBモードイメージングでは、再構成されたデータ(すなわち画像)に対してエンベロープ検出が行われる。エンベロープ検出は、例えば、ヒルベルト変換の後にマグニチュード検出を行い、オプションでローパスフィルタリングを行うことで実現できる。また、信号を2乗してローパスフィルタリングすることでも実現できる。IQデータを用いたBモードイメージングの場合、複素信号の大きさが抽出される。エンベロープ検出ステップの後、ゲイン調整とダイナミクス圧縮ステップが行われる。ドップラーやエラストグラフィの場合は、後処理を行わずに、再構成された高周波データまたはIQデータ(つまり画像)を直接使用する。
以上の説明では、本発明の概要を説明したが、いくつかの具体的な実施例を参照しながら、本発明をより詳細に説明する。
超音波イメージング、そしてより一般的なアレイ処理には、いくつかの要因によるアーチファクトがある。一般的な超音波プローブは、トランスデューサエレメントのアレイに分割された有限の開口部を有している。これらのトランスデューサエレメントは、それぞれバンドリミティッド(バンドパス)な電気音響インパルス応答と音響電気インパルス応答を持ち、電気的な加振を音圧に変換し、音圧を電気信号に記録する。このことは、組織内の理想的な点状反射体は、帯域制限された(バンドパスされた)RFパルス波形としてのみ記録できることを意味する。トランスデューサ要素の有限次元の性質は、方向性のある感度を持つ回折物体となる。プローブの開口部、開口部内でのトランスデューサ素子の配置、送信ビームを形成するために適用される遅延と振幅は、アレイの送信および受信の空間インパルス応答(SIR)を支配する要因となる。送信ビームの数と種類、およびそれに対応する遅延加算ビームフォーミング動作は、結果として得られる画像品質の重要な要因となる。超音波プローブ、トランスデューサ素子の配置、形状、物理的特性、および送受信設定は、超音波イメージング構成を定義するパラメータである。このようなイメージング構成は、よく知られている遅延加算法などの画像再構成法と組み合わせて、記録された測定値から画像を生成する。結果として得られる画像の標準的な品質指標は、点-広がり-関数(PSF)で与えられる。また、散乱や減衰などの波動伝播効果も全体の画質に影響を与える。物理的なイメージングプロセス全体には、イメージングアーティファクトの発生が伴いる。これは、特にサイドローブやグレーチングローブのレベル、減衰、残響、多重散乱によって引き起こされ、コントラストの低下、不適切な位置のゴーストイメージ、非物理的な縞模様の画像構造、シャドーイングなどとして視覚化される。
例として、表1に示す2つのイメージング構成を定義する。LQと呼ばれる低品質の構成では、128個のトランスデューサエレメントで構成されたリニアプローブを使用し、5。133MHzを中心に67%の帯域幅を確保している。つまり、アレイの全ての素子が同時に媒体中に音響波を送信することで、焦点の合わない準平面状の波面が得られる。後方散乱されたエコーは、アレイの全ての素子で測定され、そこから画像が遅延加算法によって画像推定の全てのピクセル位置で再構成される。HQと呼ばれる高品質な構成では、同じアパーチャ内に2倍の数のトランスデューサエレメントを持つ同様のプローブを使用する。画像は256回のSA測定値から再構成され、これは所定のアパーチャを用いた超音波イメージングのゴールドスタンダードと考えられる。SA測定値は次のようにして得られる。アレイの1つの素子が、焦点の合っていない波面を媒体中に送信する。アレイの全ての素子で後方散乱エコーが測定され、そこから遅延加算法によって中間高周波画像のピクセルが再構成される。この手順をアレイの各要素について繰り返し、全ての中間高周波画像を平均して最終的な高周波画像を得る。
Figure 2022503946000004
画像再構成の結果を評価すると、低品質構成のPSFは高品質構成よりもはるかにアーチファクトの影響を受けていることがわかる。画像のコントラストに大きな影響を与えるサイドローブとゴーストのアーチファクトは、高品質構成では大幅に減少している。画像の解像度に直接影響を与えるメインローブは、高品質構成ではきつくなる。一方、低品質構成では、ボケやシャドーイングに大きな影響を与えるグレーチングローブが存在している。これは、アレイピッチが波長の半分よりも大きいことに起因している。半波長ピッチの高品質構成では、これらのグレーチングローブは完全に除去されている。全体として、高品質構成で再構成された画像は、低品質構成で再構成された画像よりもはるかに高い品質を持っている。
ここで、W⊂RNを、相当数のアーチファクトが発生した撮像条件から再構成された低品質画像の部分空間とする。Nは、2次元空間のピクセル数や3次元空間のボクセル数など、画像の次元を表す。V⊂RNを高品質画像の部分空間とする。このような高品質な画像は、これらのアーチファクトを低減し、あるいは完全に除去する最適な超音波イメージング構成によって再構成されうる。また、これらのアーチファクトに影響されない別のイメージングモダリティで取得された高品質な画像や、利用可能であれば正確な物理媒体の特性である場合もある。高画質画像または高画質画像推定値(第1画像推定値とも呼ばれる)x^∈Vを、対応する低画質画像または低画質画像推定値(第2画像推定値とも呼ばれる)x~∈Wに変換するマッピングF:V→Wを以下のように定義する。
Figure 2022503946000005
部分空間VとWに応じて、多くの異なるマッピングFが存在する可能性があるが、本質的には全てのマッピングが非正規である。実際、低品質の画像に存在するアーチファクトのレベルは、高品質の画像に存在するアーチファクトのレベルよりも高いので、Fは情報の損失につながることを意味している。言い換えれば、写像Fは高品質画像x^に対してぼかし演算子として働く。したがって、対応する低品質画像x~から高品質画像x^を復元するための逆マッピングを見つけることは、最新の手法では不可能である。
例として、表1に定義されているHQおよびLQの構成を考える。結果として得られる画像の部分空間WとVを正式に定義するために、Y⊂RMLとY⊂RMHを低品質構成と高品質構成で取得した測定値の部分空間と定義し、MとMは対応する測定サンプル数である。この例では、高品質構成が低品質構成の2倍のトランスデューサ素子で構成され、256倍の送受信イベントを必要とするため、MはMよりも実質的に大きいことがわかる。これらの定義を用いると、低品質画像x~と高品質画像x^を次のように表すことができる。
Figure 2022503946000006
Figure 2022503946000007
ここで、D:Y→W、D:Y→Vは、それぞれ低品質構成と高品質構成のイメージング演算子である。したがって、低品質画像の部分空間はW=D(Y)、高品質画像の部分空間はV=D(Y)と表すことができる。上述の例では、画像再構成法として遅延加算法を用いているので、演算子DとDはともに線形演算子である。興味深いのは、高品質構成で得られたSA測定値の集合から画像を再構成するために使用される再構成法は、実際には、同じ画像の複数のビューを平均化するプロセスである。高品質構成で再構成された画像は、低品質構成で得られた測定値から再構成された画像よりもはるかに多くのビューで得られているため、本質的に多くの情報を含んでいる。したがって、式(3)で定義されるマッピングFは、この例では線形であり、例えば主成分分析(PCA)によって表現することができる。しかし、その逆マッピングを最新の方法で見つけることは、失われた情報を回復することになるため、不可能であることは明らかである。
この制限を回避するために、画像再構成タスクを逆問題として定式化するのが一般的である。ここで、未知画像x∈RNと、超音波イメージング構成によって感知された測定値y∈Rを定義すると、ここで、Mはプローブのトランスデューサ要素によって記録された測定サンプルの数を表す。これに対応する逆問題は次のように書くことができる。
Figure 2022503946000008
ここで、H:RN→Rは、物理的な測定プロセスを説明する測定モデルまたはフォワードモデル(xを入力の1つとする)であり、n∈Rは測定ノイズである。ここで重要なのは、測定モデルが非線形である可能性がある。提案された定式化は、RN(実数)の値をとる測定モデルで表現されているが、測定モデルがCN(複素数)の値をとる場合にも容易に拡張することができる。線形逆問題の特殊なケースでは、測定モデルは、H(x)=Hxとなるような行列H∈RM×Nとして表すことができる。この場合、式(6)の問題は次のようになる。
Figure 2022503946000009
式(6)の問題は一般的には非可決であるため、制約のない最適化問題を定式化することで解決することができる。
Figure 2022503946000010
汎函数DとRがともに凸である場合、式(8)の問題は、制約付き最小化問題として等価に定式化することができる
Figure 2022503946000011
ここで、C={x∈RN|R(x)≦τ}は閉凸集合であり、τは正則化パラメータλによって一意に定義される。
H(x)と測定値yの間の距離を測定するデータ忠実度汎関数Dは、しばしば二乗lノルムとして表される。
Figure 2022503946000012
未知画像xの性質に応じて、しばしばプライアと呼ばれる異なる正則化汎函数Rが文献で一般的に使用されている。lノルムの二乗、lノルム、ウェーブレット・スパース、トータル・バリエーション(TV)などが有名なプライアである。このようなプライアは、標準的な画像常に特殊であるため、超音波画像の文脈では最適とは言い難い。学習した辞書と適切なノルムを組み合わせたものも、正則化汎函数としてよく用いられる。しかし、これらは計算負荷が高く、画像のパッチにしか適用できないため、このような技術はあまり魅力的ではない。
本発明は、処理タスクにおいてその有効性が実証されているが、超音波画像の性質が非学習可能なパラメータθを有する学習非線形マッピングfθ:RN→RNを使用して、未知画像xに対する非常に効果的なデータ駆動型プライアを構築することを含む戦略を提案する。図3は、そのような非線形マッピングfθの簡略化したブロック図を示している。学習非線形マッピングは、典型的には、学習可能なパラメータθが、超音波画像アーチファクトを低減または除去するように学習ディープニューラルネットワーク(DNN)とすることができる。
学習非線形マッピングfθを用いて、データ駆動型プライアは正則化汎函数Rで表すことができる。Rの典型的な例は、画像xと残差項[x-fθ(x)]の内積として定義され、次のようになる。
Figure 2022503946000013
ここでxTはxの転置を表す。複素ベクトルx∈CNの場合、転置演算はxHと表記される共役転置(またはエルミート転置)に置き換えられるべきである。
このような正則化汎函数を定義する別の方法は、残差の二乗lノルムであり、これは明示的に次のように表される。
Figure 2022503946000014
式(11)と(12)で与えられた両方の例では、画像xがfθ(x)に近づくとR(x)がゼロに近づくことは明らかである。いくつかのシナリオでは、特にディープニューラルネットワークを考慮する場合、非線形マッピングfθを以下のような残差マッピングとして設計することが望ましい場合がある。
Figure 2022503946000015
ここで、rθ:RN→RNは、xに適用される負のノイズを予測するように学習した非線形マッピングである。
図4は、このような残差非線形マッピングのブロック図である。この場合、式(12)で定義される正則化汎函数は、例えば次のように簡略化することができる。
Figure 2022503946000016
この場合、予測される負のノイズrθ(x)がゼロになると、R(x)もゼロに近づくことは明らかである。
もう一つの一般的なニューラルネットワークアーキテクチャは、「エンコーダ-デコーダ」アーキテクチャであり、学習可能なパラメータθeを持つエンコーダeθeは、まず、入力画像xを、eθe:RN→Rとなるような潜在空間にマッピングする。学習可能なパラメータθを持つデコーダdθdは、dθd:R→RNとなるように、潜在空間を画像空間にマッピングし直す。残差形式では、このようなアーキテクチャは次のように表すことができる。
Figure 2022503946000017
ここで、学習可能なパラメータθには、θeとθの両方が含まれている。図5は、「エンコーダ-デコーダ」アーキテクチャを持つ残差非線形マッピングのブロック図である。このような学習済みの非線形マッピングfθに対して、正則化汎函数を次のように表すことができる。
Figure 2022503946000018
これは、計算量の点で有利になるかもしれない。
このデータ駆動型正則化汎函数Rを、用途に応じて適切に学習された非線形写像fθで表現すると、これを適切な画像再構成法に用いることで、従来の画像再構成法では低品質の画像しか得られなかった超音波イメージングの構成から、高品質の画像を得ることができる。提案された画像再構成法は、通常、式(8)のような最適化問題を反復アルゴリズムによって解く必要がある。式(8)を解くために、関数Dと関数Rの凸性や微分可能性などの数学的特性に応じて、多くのアルゴリズムが展開される可能性があり、そのいくつかを以下のセクションで説明する。データ忠実度汎函数に関しては、簡潔で分かりやすくするために、以下の例では式(10)に記載されているような二乗lノルムで表現されたデータ忠実度汎函数に限定するが、本発明がこれに制限されるものではない。
まず、勾配降下法では、式(8)に関して、コスト汎函数F(x)=D(x)+λR(x)を最小化するx∈RNを求める。汎函数Fがxに対して微分可能であると仮定すると、汎函数Dと汎函数Rの両方がxに対して微分可能であることを意味し、次のような勾配降下法の反復を表現することができる。
Figure 2022503946000019
ここで、γ∈Rは反復k∈N(Nは自然数を表す)ごとに変化する可能性のあるステップサイズである。適切な初期化xとステップサイズのシーケンス(γk∈Nにより、シーケンス(xk∈Nは、コスト汎函数Fの所望のローカルまたはグローバル・ミニマムに収束する可能性がある。勾配降下アルゴリズムは、収束基準が満たされるまで、式(17)で定義される反復を再帰的に適用することで構成される。図6のフローチャートは、勾配降下アルゴリズムを用いた提案された画像再構成プロセスを要約したものである。ステップ30では、プローブ3は、関心のある対象物から生じる一連の波形を受信または取得する。言い換えれば、このステップでは、一連の波形(または信号)が、物理的な励起によって刺激された(または暴露された)関心対象物からの測定値として受信または取得される。ステップ31では、媒体の未知画像と記録された一連の測定値を結びつける測定モデルが定義される。ステップ32では、Hを用いて微分可能なデータ忠実度汎函数Dを構築し、ステップ33では、fθを用いて微分可能な正則化汎函数Rを構築する。ステップ35では、アルゴリズムまたはソリューションをxとして初期化する。ステップ37では、データ忠実度勾配∇D(x)が評価または決定される。ステップ39では、正則化勾配∇R(x)が評価または決定される。ステップ41では、反復ステップサイズγが評価または決定される。ステップ43では、xk+1←x-γ[∇D(x)+λ∇R(x)]を計算することにより、推定値または解が更新される。ステップ45では、所定の収束基準が満たされているか否かが判断される。肯定的な場合、プロセスは終了する。一方、基準が満たされていない場合には、ステップ37に処理を継続する。
コスト関数Fの凸性やリプシッツ勾配などの厳密な仮定と、ステップサイズ列(γk∈N)の適切な選択により、勾配降下法の収束が保証される。このようなコスト汎函数Fに関する厳密な仮定は、満たすことが難しい場合がある。状況によっては、汎函数Dや汎函数Rが任意のx∈RNに対して微分可能でない場合もある。そのような場合には、汎函数Dおよび/または汎函数Rが凸型の非微分可能な関数である場合には、式(17)における勾配∇xD(xk) および/または∇xR(xk)をその副勾配sD(xk)∈∂xD(xk)および/またはsR(xk)∈∂xR(xk)で置き換えることができる。汎函数Dおよび/または汎函数Rが局所的なリプシッツ連続性を有する非凸の非微分可能な汎函数である場合、式(17)の勾配∇xD(xk) および/または∇xR(xk) を、その一般化勾配gD(xk)∈∂xD(xk) および/またはgR(xk)∈∂xR(xk)に置き換えることができ、これは、非特許文献1:Clarke F., (2013) Generalized gradients, "Functional Analysis, Caliculus of Variations and Optimal Control", Graduate Texts in Mathematics, vol. 264. Springer, London.によって定義されている。
データ忠実度汎函数は、線形測定モデルを用いて、よく保存された凸およびリプシッツ勾配関数として構築するのが一般的だが、ニューラルネットワークなどの学習非線形マッピングfθを用いて正則化汎函数Rを構築する場合には、このような特性を容易に検証できない場合がある。この場合、微分可能な活性化関数のみを用いることで、fθの微分可能性を強制することができる。一般的な整流線形単位(ReLU)などの非微分活性化関数を使用し、局所的なリプシッツ連続性を確保することができれば、上で定義した一般化された勾配の恩恵を受けることができる。局所的なリプシッツ連続性が確保できない場合でも、ニューラルネットワークの学習で一般的に行われているように、不連続な位置での勾配∇xR(xk)を近似することができる。また、初期化xがあると、特に有利になる。例えば、よく知られている遅延加算法で得られた画像を使用すると、式(17)で定義される反復に効率的な初期化xを提供することができる。勾配降下法のアルゴリズムは、そのペースを速めたり、収束を保証するためのFに関する仮定を緩和したりするために、多くのバリエーションが開発されている。このような変形は、もちろん容易に適用することができる。
式(17)で表される勾配降下反復は、DとRの両方の勾配を含みる。例として、式(10)で定義される二乗lノルムデータフィデリティを考えると、この勾配は次のように表される。
Figure 2022503946000020
ここで、∇xH(x)は、以下のヤコビ行列の転置行列である。
Figure 2022503946000021
二乗lノルムデータ忠実度汎函数Dの勾配は、以下のように表すことができる。
Figure 2022503946000022
正則化汎函数Rの勾配は、特に興味深い。例えば、式(11)で定義された正則化汎函数Rの勾配∇Rは次のように表される。
Figure 2022503946000023
ここで、∇θ(x)xは、方向xに沿った位置xにおける学習非線形マッピングfθの方向微分であり、∇θ(x)は、次のように定義されるヤコビアン行列の転置行列である。
Figure 2022503946000024
しかし、各反復kでヤコビアン行列Jfθ(x)を構築することは、計算上非常に困難である。式(17)に記載されているように、各反復において、項∇θ(x)xの評価のみが必要であることがわかる。この評価がヤコビアン行列を構築することなく可能であれば、反復を効率的に処理することができる。ニューラルネットワークをfθとして使用する際に注目すべき点は、fθが本質的に微分可能であること、あるいは少なくともその勾配を効率的に近似する手段が容易に入手できることである。この微分可能性は、チェーンルールを用いた強力なバックプロパゲーション法に依存する学習戦略の重要な構成要素である。したがって、計算効率の高いバックプロパゲーションアルゴリズムを利用することで、ヤコビアン行列Jfθ(x)xを明示的に構築することなく、行列-ベクトル積∇θ(x)xを直接評価することができる。この戦略は、方向微分の「行列なし」の評価と見なすことができ、学習非線形マッピングfθにニューラルネットワークを使用する際の明確な利点となる。
2つ目の例として、式(14)で定義される正則化汎函数Rを考えると、その勾配は次のように表される。
Figure 2022503946000025
ここで、∇θ(x)rθ(x)は、方向rθ(x)に沿った位置xにおける学習非線形残差マッピングrθの方向微分である。前述の例と同様に、効率的なバックプロパゲーションアルゴリズムを利用することで、各反復kにおいて行列-ベクター積∇θ(x)rθ(x)を評価することができる。
次に、投影法を用いた方法を検討する。近似勾配法は、DとRがともに凸状汎函数で、Dが微分可能、Rが非微分可能である場合に、式(8)の問題を解くのに使用できる。対応するアルゴリズムのイテレーションは次のように与えられる。
Figure 2022503946000026
ここで、γ∈Rは、反復k∈Nにおける勾配ステップサイズである。式(24)は、まず、測定値への近さを確保するデータ忠実度に関する勾配ステップを実行し、次に、未知画像に関する事前知識を考慮する近位ステップを実行する。式(24)に関与する近位演算子はいくつかのμ∈Rに関して、次のように定義される。
Figure 2022503946000027
ここでz∈RNは、中間画像を参照する。正則化汎函数Rの構築に使用される学習非線形マッピングfθに応じて、問題(25)は閉形式解を有する場合があり、その場合には容易に適用することができる。そのような閉形式の解がない場合は、式(25)を反復的に解くことができるが、これは計算量が多くなる可能性がある。
式(25)で定義され、式(24)に関与する近位演算子をさらに詳しく見ると、この演算は概念的には画像xのノイズ除去であることがわかる。例えば、S. V. Venkata Krishnan、 Ch. A. Bouman and Br. Wohlbergの”Plug-and-Play priors for model based reconstruction”、 in 2013 IEEE Global Conference on Signal and Information Processing、 2013、 pp.945-94に説明されているようなプラグアンドプレイの事前アプローチを利用すると、近位演算子を適切なノイズ除去器で置き換えることができる。したがって、近位演算子を学習非線形マッピングfθで置き換えることで、結果として得られる画質を大幅に改善することができる。式(24)で定義される反復は次のようになる。
Figure 2022503946000028
θがある閉凸集合Cへのプロジェクタである特別な場合、式(26)に記載された反復は、よく知られた投影勾配降下アルゴリズムに対応し、式(9)で定義された制約付き最適化問題を解くために使用することができることに留意することができる。したがって、式(26)で定義された反復は、fθがある閉凸集合Cへのプロジェクタである場合に投影勾配降下アルゴリズムに還元される、汎用の投影ベースのアルゴリズムを形成するために使用することができる。一般的な投影ベースのアルゴリズムは、式(26)で定義された反復を、収束基準が満たされるまで再帰的に適用することで構成される。図7のフローチャートは、ジェネリック・プロジェクション・ベース・アルゴリズムを用いた画像再構成プロセスの提案を示している。ステップ50では、プローブ3は、関心対象物から生じる一連の波形を受信または取得する。言い換えれば、このステップでは、物理的な励起によって刺激された(または暴露された)関心対象物からの測定値として、一組の波形(または信号)が受信または取得される。ステップ51では、記録された一連の測定値に媒体の未知画像をリンクする測定モデルが定義される。ステップ52では、微分可能または副微分可能なデータ忠実度汎函数DがHを用いて構築または定義される。ステップ55では、データ忠実度勾配∇D(x)が評価または決定される。ステップ57では、反復ステップサイズγが評価または決定される。ステップ59において、投影x_(k+1)←fθ(x-γD(x))を適用することにより、推定値または解が更新される。ステップ61では、所定の収束基準が満たされているか否かが判断される。肯定的な場合、処理は終了する。一方、基準が満たされていない場合には、ステップ55で処理を続行する。学習非線形マッピングfθが一般的なマッピングである場合、それはプラグアンドプレイアプローチと同様である。学習非線形マッピングfθが、ある閉凸集合へのプロジェクタである場合、投影ベースのアルゴリズムは、よく知られた投影勾配降下法アルゴリズムと同等である。
式(26)をステップサイズγ=1、初期化x=0でk=0回繰り返した場合の特殊なケースが注目される。式(10)で定義される二乗lノルムデータ忠実度汎函数Dを持つ線形測定ケースを仮定し、その勾配は式(20)で表され、推定値x^は次のように与えられる。
Figure 2022503946000029
Figure 2022503946000030
Tyという演算は、線形測定演算子Hに関する特定の仮定がなされた場合、実際には測定値yに適用されるよく知られた遅延加算演算であることが証明されている。したがって、この非常に特殊なケースでは、高品質の画像推定値x^は、単純な直接マッピングによってx~=HTyとなるような遅延加算法によって再構成された低品質の画像推定値x~から得ることができる。
Figure 2022503946000031
このように、まず低画質画像の推定値x~=HTyを計算し、次にx^=fθ(x~)となるようにマッピングを適用して推定値を計算する。
このように単純化することで、遅延加算法による低画質画像の再構成とfθの1回の評価で済むという大きなメリットがある。したがって、学習非線形マッピングfθは、式(3)で定義されたマッピングFの逆数に近似している。単純化されているとはいえ、この特殊なケースでは、再構成された高品質画像x^が、低品質画像x~に比べて大幅に改善されている。
上述の画像再構成の例(勾配降下法と投影法)はいずれも、反復アルゴリズムの各反復において、グローバルに学習非線形マッピングfθに依存していることに留意できる。しかし、このアプローチを、各fθkが反復アルゴリズムの各反復k∈Nにおいて存在する所与のレベルのアーチファクトに対して特別に設計し、学習非線形マッピング{fθkk∈Nのセットに容易に特化することができる。この場合、画像再構成法は、少なくとも2つの学習非線形マッピング(学習非線形マッピングの多重度)に基づいて行われることになる。
上記では、式(8)のような最適化問題を解決するための一般的な方法を説明した。ここでは、学習非線形マッピングfθを用いて正則化汎函数Rを構築する。以下では、提案された手法を十分に活用するために、このような非線形マッピングをどのように設計し、学習するかについて、いくつかの洞察を示す。
例として、図5に示した残差「エンコーダ-デコーダ」アーキテクチャが考えられる。このようなアーキテクチャは、図8に示されているような深層残差畳み込みニューラルネットワークとして実装されてもよい。入力画像xの符号化は、ニューラルネットワークアーキテクチャの下向きの経路(左上から中心下へ)で行われ、復号化は、中間残差出力rθ(x)まで、ニューラルネットワークアーキテクチャの上向きの経路(中心下から右上へ)で行われる。その後、残差出力が入力画像に加算され、最終出力fθ(x)=x+rθ(x)が形成される。実線の矢印は演算を指定し、破線の矢印は固有のスキップ接続を指定している。破線の矢印は、操作ではなく、量の移動を表している。点線のボックスは、それぞれ分解レベルを表している。この例では、3つの分解レベル、すなわちl=0、1、2が考えられる。実線で描かれた長方形は、NNのそれぞれの深さにある特徴マップ、つまり演算やレイヤーの出力を表している。各ボックスの幅(水平線に沿った寸法)は、特徴マップに含まれるチャンネル番号を示している。最初の演算では、カーネルサイズが(1×1)の畳み込み演算を用いて、入力画像xのチャネル番号をユーザ定義可能なチャネル番号Nに拡張する。残差出力rθ(x)の前の最後の演算では、カーネルサイズが(1×1)の畳み込み演算を用いて、チャンネル番号Nを入力画像のチャンネル番号に戻す縮退が行われる。各分解レベル内では、(3×3)の畳み込み演算に続いて、非線形活性化(図8のnl.)が、特徴マップの次元を変えずに各特徴マップ間で実行される。このような非線形活性化は、例えば、周知の整流線形単位(ReLU)、指数線形単位(ELU)、シグモイド関数などである。下降経路の各分解レベルの間では、チャネル数の拡張を伴う(または伴わない)マルチスケール(3×3)の畳み込み演算を用いてダウンサンプリング演算が行われ、その後、非線形活性化が行われる。このようなマルチスケール演算は、例えば、最大プーリング演算の後にコンボリューション演算を行うものであってもよい。また、チャンネル数の拡大を伴う(または伴わない)ストライドコンボリューションでもよい。上昇経路の各分解レベルの間では、マルチスケール(3×3)の畳み込み演算とそれに続く非線形活性化を用いてアップサンプリング演算が行われる。このようなマルチスケール演算は、例えば、チャンネル数の縮小を伴う(または伴わない)転置型およびストライド型の畳み込み演算であってもよい。
ニューラルネットワークアーキテクチャのフィルタ係数は、全て学習可能なパラメータである。これらのパラメータ(重みとも呼ばれる)は、経験的なリスクを最小化することで学習することができる。
Figure 2022503946000032
ここでLは学習損失関数であり、{x^j}j=1 Lと{x~j}j=1 LはL枚の画像の集合である。集合{x^j}j=1 Lと{x~j}j=1 Lは、高品質部分空間Vと低品質部分空間Wからそれぞれサンプリングして構成するのが一般的である。このようなサンプルは、例えば、実験や数値シミュレーションによって得られる。セット{x^j、x~j }j=1 Lは、しばしばトレーニングセットと呼ばれる。確率的勾配降下法(SGD)、バッチSGDまたはAdamオプティマイザー(D. P. Kingma and J. Ba、 "Adam: A Method for Stochastic Optimization", pp.1-15, 2014. arXiv: 1412.6980)などの最適化アルゴリズムは、例えば、式(30)で定義される最小化問題を解くために使用されてもよい。この目的のために、学習可能なパラメータ∇θLに対する損失関数の勾配は、例えば、一般的なチェーンルールを用いたバックプロパゲーションによって計算することができる。典型的な損失関数は、次のように定義される平均二乗誤差(MSE)である。
Figure 2022503946000033
ただし、異なる損失関数、リスク最小化関数、最適化アルゴリズムを使用してもよい。
非線形マッピングfθとしては、O. Ronneberger, Ph. Fischer, and. Brox "U-Net: Convolutional Networks for Biomedical Image Segmentation", pp. 1-8, 2015. arXiv: 1505.04597によってバイオメディカルな画像セグメンテーションの目的で最初に紹介され、最近ではスパースビューCT画像再構成の文脈で使用されている、U-Netと呼ばれるニューラルネットワークアーキテクチャを一例として使用することができる。これは本質的に、図8に示された残差非線形マッピングに非常に似ている。fθのパラメータθは、Adam optimiserを用いて最適化され、一方で、前述のMSEが学習損失として使用される。
学習セットは、例えば特定のシミュレータを使用して生成された約1000~100000または10000~30000の画像で構成されていてもよい。この目的のために合成ファントムを使用してもよく、様々なサイズとエコー源性の介在物をランダムに生成してもよい。介在物の大きさとエコー源性の確率分布は、対象物に存在する分布を忠実に再現するようにモデル化することができる。数値的なファントムを設計し、画像全体に完全に拡散したスペックルパターンを得るために使用することができる。実験データではなくシミュレーションデータを使用することの主な利点は以下の通りである。異なったエコー源性を持つアーチファクトや介在物の分布を完全にコントロールできる。特定のアーチファクトが学習セットにほとんど含まれていない場合、ニューラルネットワークはそれを除去することにほとんど焦点を当てないため、これはニューラルネットワークの学習に役立つ。また、実験データと比較して、トレーニングセットに多様性を持たせることも可能である。また、シミュレーション環境が完全に制御されているため、組織の動きによるボケを防ぐことができる。もう一つの利点は、既存の超音波イメージング構成ではないものをシミュレートできることである。例えば、表1に示した低品質構成は現実的なプローブに基づいているが、高品質構成で使用されるプローブは存在しないため、与えられたアパーチャに対して理想的な画像を再構成することができる。シミュレーションされたトレーニングセットを使用することの主な欠点は、実験データに適用したときにその結果が反映されるかどうかが自明ではないという点にある。
学習した後は、例えば、式(29)で示される直接マッピングアプローチを使用することができる。このようにして、学習非線形マッピングfθを低品質の画像推定値x~に適用するだけで、高品質の画像推定値x^を得ることができる。
要約すると、本発明は、超高速超音波イメージングにおいて、高品質の超音波イメージング構成から画像を再構成することは、代償を伴い、実際には実現できないかもしれないという事実に動機づけられている。実際、高品質の超高速画像は、一般に、順次再構成された低品質の複数の画像を平均化または複合化することによって得られる。この戦略の例として、よく知られているコヒーレントPW合成がある。ここでは、ステアされたPWが順次送信され、各ステアされたPWについて再構成された全ての画像を平均することで最終的な画像が得られる。送受信イベントが連続しているため、非常に高いフレームレートが要求されるシナリオでは、音響波の往復の飛行時間がすでに制限要因となっている。また、3Dイメージングや低消費電力の環境など、多くの画像を合成して再構成するための計算要件も、状況によっては実行できないことがある。さらに、物体(人間の組織など)の動きによって、異なる時間間隔で再構成された連続画像を平均化する際に、結果として得られる画像の実質的なぼやけが誘発される可能性がある。したがって、本発明の目的の一つは、単一または非常に少数のPW測定で高品質の画像を再構成することである。
本発明は、図面および前述の説明において詳細に図示および説明されているが、このような図示および説明は、例示的または模範的なものであり、制限的なものではないと考えられ、本発明は開示された実施形態に限定されるものではない。他の実施形態および変形は、図面、開示および添付の請求項の検討に基づいて、請求された発明を実施する際に、当業者によって理解され、達成され得るものである。例えば、本発明の教示は、送信された音響波が電磁放射線パルスに置き換えられた光音響イメージングに適用することができる。
特許請求の範囲では、「備えた」という言葉は他の要素やステップを除外するものではなく、単数表示は複数を除外するものではない。異なる特徴が相互に異なる従属請求項に記載されているという事実だけで、これらの特徴の組み合わせが有利に使用できないことを示すものではない。

Claims (15)

  1. 対象物から測定値yとして1セットの波形を受信する受信ステップと、
    前記対象物の未知画像xと前記測定値yとを結びつける測定モデルH(x)を定義する定義ステップと、
    前記測定モデルH(x)と前記測定値yとの間の距離を測定するデータ忠実度汎函数D(x)を決定する決定ステップと、
    入力データに関して非線形な、学習非線形マッピングfθを用いて、前記未知画像に関する事前情報を含む正則化汎函数R(x)を決定する第2決定ステップと、
    前記未知画像の第1画像推定値x^を得るために、前記データ忠実度汎函数Dと正則化汎函数Rとを含む最適化問題を定義する定義ステップと、
    前記最適化問題を解いて、前記第1画像推定値x^を算出する算出ステップと、
    を備えた画像再構成方法。
  2. 前記データ忠実度汎函数D(x)および正則化汎函数R(x)を定義する前に、トレーニングデータセットを用いて前記非線形マッピングをトレーニングするステップをさらに含む請求項1に記載の画像再構成方法。
  3. 前記トレーニングデータセットは、実験または数値シミュレーションによって取得された画像を含む請求項2に記載の画像再構成方法。
  4. 前記学習非線形マッピングは、学習済みの人工ニューラルネットワークである請求項1~3のいずれか1項に記載の画像再構成方法。
  5. 前記最適化問題は、x^=argminx∈RN{D(x)+λR(x)}で表され、λ∈Rは、正則化汎函数Rの重み付けとなる正則化パラメータである請求項1~4のいずれか1項に記載の画像再構成方法。
  6. 前記データ忠実度汎函数D(x)および前記正則化汎函数R(x)の少なくともいずれか一方が、凸性、非凸性、微分可能性および非微分可能性のうちの少なくとも1つの特性を有する請求項1~5のいずれか1項に記載の画像再構成方法。
  7. 反復アルゴリズムを実行することにより前記最適化問題を解く請求項1~6のいずれか1項に記載の画像再構成方法。
  8. 前記反復アルゴリズムは、勾配降下法またはジェネリック・プロジェクション・ベース法である請求項7に記載の画像再構成方法。
  9. 前記反復アルゴリズムの反復ステップサイズを1に設定する反復ステップサイズ設定ステップと、
    初期パラメータをゼロに設定する初期パラメータ設定ステップと、
    前記反復アルゴリズムとして1回の反復のみを盛り込むステップと、
    前記対象物の第1画像推定値x^を得るために、前記学習非線形マッピングを第2画像推定値x~に直接適用するステップと、
    をさらに含み、
    前記第2画像推定値x~は、前記第1画像推定値x^よりも画像アーチファクトの点で低品質である請求項7または8に記載の画像再構成方法。
  10. 前記反復アルゴリズムの1つの反復ステップから他の反復ステップへと前記学習非線形マッピングが変化する請求項7~9のいずれか1項に記載の画像再構成方法。
  11. 前記対象物に前記1セットの波形を送信するステップをさらに含む請求項7または8に記載の画像再構成方法。
  12. 送信される前記波形の数は、1以上5以下、1以上3以下、または1である請求項7または8に記載の画像再構成方法。
  13. 前記画像再構成法は超音波画像再構成法であること請求項1~12のいずれか1項に記載の画像再構成方法。
  14. 前記正則化汎函数R(x)は、以下のいずれかである請求項1~13のいずれか1項に記載の画像再構成方法。
    R(x)=1/2x[x-fθ(x)];
    R(x)=1/2||x-fθ(x)||
    R(x)=1/2||rθ(x)|| 、(ここで、rθ:R→Rは、xに適用される負のノイズを予測するように学習された学習非線形マッピング)
    R(x)=1/2||eθe(x)|| (ここでeθeは学習可能なパラメータθを持つエンコーダ)
  15. 対象物からの測定値yとして、一連の波形を受信する受信部と、
    対象物の未知画像xと測定値を結びつける測定モデルH(x)を定義する第1定義部と、
    前記測定モデルH(x)とyとの間の距離を測定してデータ忠実度汎函数D(x)を決定する決定部と、
    入力データに関して非線形な、学習非線形マッピングfθを用いて、前記未知画像xに関する事前情報を含む正則化汎函数R(x)を決定する第2決定部と、
    前記未知画像の第1画像推定値x^を得るために、前記データ忠実度汎函数D(x)と正則化汎函数R(x)とを含む最適化問題を定義する第2定義部と、
    前記最適化問題を解き、前記第1画像推定値x^を算出する算出部と、
    を備えた画像再構成装置。
JP2021517765A 2018-10-08 2019-08-29 学習非線形マッピングに基づく画像再構成方法 Active JP7345208B2 (ja)

Applications Claiming Priority (3)

Application Number Priority Date Filing Date Title
EP18199166.2 2018-10-08
EP18199166.2A EP3637099B1 (en) 2018-10-08 2018-10-08 Image reconstruction method based on a trained non-linear mapping
PCT/EP2019/073152 WO2020074181A1 (en) 2018-10-08 2019-08-29 Image reconstruction method based on a trained non-linear mapping

Publications (2)

Publication Number Publication Date
JP2022503946A true JP2022503946A (ja) 2022-01-12
JP7345208B2 JP7345208B2 (ja) 2023-09-15

Family

ID=63794402

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2021517765A Active JP7345208B2 (ja) 2018-10-08 2019-08-29 学習非線形マッピングに基づく画像再構成方法

Country Status (7)

Country Link
US (1) US11846608B2 (ja)
EP (1) EP3637099B1 (ja)
JP (1) JP7345208B2 (ja)
KR (1) KR20210069655A (ja)
ES (1) ES2886155T3 (ja)
IL (1) IL281888B2 (ja)
WO (1) WO2020074181A1 (ja)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
DE102019204287A1 (de) 2019-03-27 2020-10-01 Siemens Healthcare Gmbh Vorrichtung und Verfahren zur Kontrolle von Akquisitionsparametern bei der Durchführung einer medizinischen Röntgenuntersuchung
CN111783616B (zh) * 2020-06-28 2024-03-26 北京瓦特曼科技有限公司 一种基于数据驱动自学习的无损检测方法
CN111861961B (zh) * 2020-07-25 2023-09-22 安徽理工大学 单幅图像超分辨率的多尺度残差融合模型及其复原方法
CN112674794B (zh) * 2020-12-21 2023-02-10 苏州二向箔科技有限公司 一种结合深度学习与吉洪诺夫正则化反演的超声ct声速重建方法
US11933765B2 (en) * 2021-02-05 2024-03-19 Evident Canada, Inc. Ultrasound inspection techniques for detecting a flaw in a test object
CN113014358B (zh) * 2021-03-05 2022-06-28 苏州至善视界信息技术有限公司 一种应用于反向散射系统的自适应解码方法及解码系统
CN115018943B (zh) * 2022-05-31 2024-02-20 合肥工业大学 一种基于无训练的深度级联网络的电磁逆散射成像方法
EP4285835A1 (en) * 2022-06-03 2023-12-06 Helmholtz Zentrum München - Deutsches Forschungszentrum für Gesundheit und Umwelt (GmbH) Methods and system for optoacoustic and/or ultrasonic imaging, reconstructing optoacoustic and/or ultrasonic images and training an artificial neural network provided therefor
CN115236208B (zh) * 2022-06-27 2023-04-07 哈尔滨工业大学 一种基于信息增强与变步长稀疏表达的钢轨健康监测方法
CN115100069A (zh) * 2022-07-12 2022-09-23 清华大学 超声图像重建方法、系统、设备及介质
CN116612206B (zh) * 2023-07-19 2023-09-29 中国海洋大学 一种利用卷积神经网络减少ct扫描时间的方法及系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2018099881A1 (en) * 2016-11-30 2018-06-07 Koninklijke Philips N.V. Bone and hard plaque segmentation in spectral ct

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8427650B2 (en) * 2008-12-02 2013-04-23 Opteryx, Llc Reconstruction of nonlinear wave propagation
WO2012173095A1 (ja) * 2011-06-13 2012-12-20 株式会社東芝 磁気共鳴イメージング装置及びその制御装置
US9664765B2 (en) * 2011-06-30 2017-05-30 Hitachi, Ltd. Magnetic resonance imaging apparatus and gradient magnetic field waveform estimation method
WO2013009944A1 (en) * 2011-07-12 2013-01-17 Colorado School Of Mines Wave-equation migration velocity analysis using image warping
CN107369189A (zh) * 2017-07-21 2017-11-21 成都信息工程大学 基于特征损失的医学图像超分辨率重建方法
EP3447525B1 (en) * 2017-08-23 2020-05-06 Ecole Polytechnique Federale De Lausanne (Epfl) Model-based image reconstruction method

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2018099881A1 (en) * 2016-11-30 2018-06-07 Koninklijke Philips N.V. Bone and hard plaque segmentation in spectral ct

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
EUNHEE KANG ET AL.: "A deep convolutional neural network using directional wavelets for low-dose X-ray CT reconstruction", MEDICAL PHYSICS, vol. Volume 44, Issue 10, JPN6023004322, 13 October 2017 (2017-10-13), pages 360 - 375, ISSN: 0005104780 *

Also Published As

Publication number Publication date
US20210341436A1 (en) 2021-11-04
ES2886155T3 (es) 2021-12-16
EP3637099A1 (en) 2020-04-15
WO2020074181A1 (en) 2020-04-16
JP7345208B2 (ja) 2023-09-15
CN112771374A (zh) 2021-05-07
IL281888A (en) 2021-05-31
KR20210069655A (ko) 2021-06-11
EP3637099B1 (en) 2021-06-23
US11846608B2 (en) 2023-12-19
IL281888B2 (en) 2023-11-01
IL281888B1 (en) 2023-07-01

Similar Documents

Publication Publication Date Title
JP7345208B2 (ja) 学習非線形マッピングに基づく画像再構成方法
Luijten et al. Adaptive ultrasound beamforming using deep learning
Zhang et al. Ultrasound image reconstruction from plane wave radio-frequency data by self-supervised deep neural network
US11693113B2 (en) Quantitative ultrasound imaging based on seismic full waveform inversion
Besson et al. Ultrafast ultrasound imaging as an inverse problem: Matrix-free sparse image reconstruction
Laroche et al. An inverse approach for ultrasonic imaging from full matrix capture data: Application to resolution enhancement in NDT
KR20210042907A (ko) 초음파를 이용하여 이종 매체를 비침습적으로 특성화하기 위한 방법 및 시스템
CN109875606B (zh) 基于先验反射成像的超声ct声速成像的方法
CN112912758A (zh) 用于对超声信号进行自适应波束形成的方法和系统
JP6924530B2 (ja) モデルベース画像再構成方法
Schretter et al. Ultrasound imaging from sparse RF samples using system point spread functions
Ramkumar et al. Compressed sensing approach for reducing the number of receive elements in synthetic transmit aperture imaging
Shpigler et al. Detection of overlapping ultrasonic echoes with deep neural networks
Pilikos et al. Fast ultrasonic imaging using end-to-end deep learning
Szasz Advanced beamforming techniques in ultrasound imaging and the associated inverse problems
CN112771374B (zh) 基于训练的非线性映射的图像重建方法
Guenther et al. Robust finite impulse response beamforming applied to medical ultrasound
CN112041699B (zh) 重建系统和方法
Hourani et al. Joint deconvolution of fundamental and harmonic ultrasound images
Bangalore Narasimha Prasad Ultrasound Computed Tomography using Deep Learning
Carlson et al. High resolution image reconstruction from full-matrix capture data using minimum mean square error deconvolution of the spatio-temporal system transfer function
Hourani Fundamental and Harmonic Ultrasound Image Joint Restoration
Ben Daya Compensated Row-Column Ultrasound Imaging System
Rodríguez et al. The signal projection method: A general approach to active imaging with arrays
Han et al. Numerical estimation of femoral neck cortical bone thickness based on time domain topological energy and sparse signal approximation

Legal Events

Date Code Title Description
A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A821

Effective date: 20210330

A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20220712

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20230207

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20230417

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20230711

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20230718

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20230801

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20230829

R150 Certificate of patent or registration of utility model

Ref document number: 7345208

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150