以下、図面を参照して、実施形態に係る超音波診断装置、医用画像処理装置及び医用画像処理プログラムを説明する。なお、実施形態は、以下の実施形態に限られるものではない。また、一つの実施形態に記載した内容は、原則として他の実施形態にも同様に適用される。
(第1の実施形態)
まず、第1の実施形態に係る超音波診断装置の構成について説明する。図1は、第1の実施形態に係る超音波診断装置の構成例を示すブロック図である。図1に例示するように、第1の実施形態に係る超音波診断装置は、超音波プローブ1と、モニタ2と、入力装置3と、装置本体10とを有する。
超音波プローブ1は、超音波の送受信を行なうために、装置本体10に接続される。超音波プローブ1は、例えば、複数の圧電振動子を有する。これら複数の圧電振動子は、装置本体10が有する送受信回路11から供給される駆動信号(駆動パルス)に基づき超音波を発生する。また、超音波プローブ1が有する複数の圧電振動子は、被検体Pからの反射波(エコー)を受信し、受信した反射波を電気信号(反射波信号)に変換し、反射波信号を装置本体10に送信する。また、超音波プローブ1は、圧電振動子に設けられる整合層と、圧電振動子から後方への超音波の伝播を防止するバッキング材等を有する。なお、超音波プローブ1は、装置本体10と着脱自在に接続される。
超音波プローブ1から被検体Pに超音波が送信されると、送信された超音波は、被検体Pの体内組織における音響インピーダンスの不連続面で次々と反射され、反射波として超音波プローブ1が有する複数の圧電振動子にて受信される。受信される反射波の振幅は、超音波が反射される不連続面における音響インピーダンスの差に依存する。なお、送信された超音波パルスが、移動している血流や心臓壁等の表面で反射された場合の反射波は、ドプラ効果により、移動体の超音波送信方向に対する速度成分に依存して、周波数偏移を受ける。
なお、第1の実施形態は、超音波プローブ1が、被検体Pを2次元で走査する1Dアレイプローブであっても、被検体Pを3次元で走査するメカニカル4Dプローブや2Dアレイプローブであっても適用可能である。
入力装置3は、マウス、キーボード、ボタン、パネルスイッチ、タッチコマンドスクリーン、フットスイッチ、トラックボール、ジョイスティック等により実現される。入力装置3は、超音波診断装置の操作者からの各種設定要求を受け付け、受け付けた各種設定要求を装置本体10に転送する。
モニタ2は、超音波診断装置の操作者が入力装置3を用いて各種設定要求を入力するためのGUI(Graphical User Interface)を表示したり、装置本体10において生成された超音波画像データが示す超音波画像等を表示したりする。モニタ2は、液晶モニタやCRT(Cathode Ray Tube)モニタ等によって実現される。
装置本体10は、超音波プローブ1から送信される反射波信号に基づいて超音波画像データを生成する装置である。図1に示す装置本体10は、2次元の反射波信号に基づいて2次元の超音波画像データを生成可能であり、3次元の反射波信号に基づいて3次元の超音波画像データを生成可能な装置である。ただし、第1の実施形態では、装置本体10が、2次元データ専用の装置である場合であっても適用可能である。
装置本体10は、図1に例示するように、送受信回路11と、バッファ12と、Bモード処理回路13と、ドプラ処理回路14と、画像生成回路15と、画像メモリ16と、内部記憶回路17と、処理回路18とを有する。
送受信回路11は、後述する処理回路18の指示に基づいて、超音波プローブ1が行なう超音波走査を制御する。なお、超音波走査とは、例えば、超音波送受信を指す。送受信回路11は、パルス発生器、送信遅延回路、パルサ等を有し、超音波プローブ1に駆動信号を供給する。パルス発生器は、所定の繰り返し周波数(PRF:Pulse Repetition Frequency)で送信超音波を形成するためのレートパルスを繰り返し発生する。また、送信遅延回路は、超音波プローブ1から発生される超音波をビーム状に集束し、かつ送信指向性を決定するために必要な圧電振動子ごとの遅延時間を、パルス発生器が発生する各レートパルスに対し与える。また、パルサは、レートパルスに基づくタイミングで、超音波プローブ1に駆動信号を印加する。すなわち、送信遅延回路は、各レートパルスに対し与える遅延時間を変化させることで、圧電振動子面から送信される超音波の送信方向を任意に調整する。
なお、送受信回路11は、処理回路18の指示に基づいて、所定のスキャンシーケンスを実行するために、送信周波数、送信駆動電圧等を瞬時に変更可能な機能を有している。特に、送信駆動電圧の変更は、瞬間にその値を切り替え可能なリニアアンプ型の発信回路、又は、複数の電源ユニットを電気的に切り替える機構によって実現される。
例えば、送受信回路11は、処理回路18の制御により、フレーム間のデータ列をドプラデータ列として収集する超音波走査を超音波プローブ1に実行させる(例えば、特許第3724846号,特開2014-42823号公報等を参照)。例えば、送受信回路11は、処理回路18の制御により、第1走査範囲内の移動体の運動に関する情報を取得する第1超音波走査を超音波プローブ1に実行させ、第2走査範囲内の組織形状の情報を取得する第2超音波走査として当該第2走査範囲を分割した複数の分割範囲それぞれの超音波走査を、第1超音波走査の間に時分割で超音波プローブ1に実行させる。なお、第1超音波走査は、1フレームに同一走査線を1回走査する超音波走査であってもよい。
また、送受信回路11は、アンプ回路、A/D(Analog/Digital)変換器、受信遅延回路、加算器、直交検波回路等を有し、超音波プローブ1から送信される反射波信号に対して各種処理を行って反射波データを生成する。そして、送受信回路11は、生成した反射波データをバッファ12に格納する。アンプ回路は、反射波信号をチャンネル毎に増幅してゲイン補正処理を行う。A/D変換器は、ゲイン補正された反射波信号をA/D変換する。受信遅延回路は、デジタルデータに受信指向性を決定するのに必要な受信遅延時間を与える。加算器は、受信遅延回路により受信遅延時間が与えられた反射波信号の加算処理を行う。加算器の加算処理により、反射波信号の受信指向性に応じた方向からの反射成分が強調される。なお、各素子の反射波信号ごとに受信遅延による位相調整を行い、加算する処理のことを整相加算処理或いはビームフォーミング処理とも言う。
そして、直交検波回路は、加算器の出力信号をベースバンド帯域の同相信号(I信号、I:In-phase)と直交信号(Q信号、Q:Quadrature-phase)とに変換する。そして、直交検波回路は、I信号及びQ信号(以下、IQ信号と記載する)を反射波データとして、バッファ12に格納する。なお、直交検波回路は、加算器の出力信号を、解析信号に変換した上で、バッファ12に格納しても良い。IQ信号や、解析信号は、位相情報が含まれる信号(受信信号)となる。以下では、送受信回路11が出力する反射波データを、受信信号と記載する場合がある。
送受信回路11は、被検体Pを2次元走査する場合、超音波プローブ1から2次元の超音波ビームを送信させる。そして、送受信回路11は、超音波プローブ1から送信される2次元の反射波信号から2次元の反射波データを生成する。また、送受信回路11は、被検体Pを3次元走査する場合、超音波プローブ1から3次元の超音波ビームを送信させる。そして、送受信回路11は、超音波プローブ1から送信される3次元の反射波信号から3次元の反射波データを生成する。
バッファ12は、送受信回路11が生成した反射波データを一時的に記憶するメモリである。具体的には、バッファ12は、数フレーム分の反射波データ、又は、数ボリューム分の反射波データを記憶する。例えば、バッファ12は、FIFO(First-In/First-Out)メモリであり、送受信回路11による制御を受けて、所定フレーム分の反射波データを記憶する。そして、例えば、バッファ12は、新たに1フレーム分の反射波データが送受信回路11により生成された場合、送受信回路11による制御を受けて、生成された時間が最も古い1フレーム分の反射波データを破棄して、新たに生成された1フレーム分の反射波データを記憶する。例えば、バッファ12は、RAM(Random Access Memory)、フラッシュメモリ等の半導体メモリ素子によって実現される。
Bモード処理回路13及びドプラ処理回路14は、送受信回路11が反射波信号から生成した反射波データに対して、各種の信号処理を行なう信号処理部である。Bモード処理回路13及びドプラ処理回路14は、例えば、プロセッサにより実現される。Bモード処理回路13は、バッファ12から反射波データを読み出し、読み出した反射波データに対して、対数増幅、包絡線検波処理、対数圧縮などを行なって、複数のサンプル点それぞれの信号強度が輝度の明るさで表現されるデータ(Bモードデータ)を生成する。
なお、Bモード処理回路13は、フィルタ処理により、検波周波数を変化させることで、映像化する周波数帯域を変えることができる。このBモード処理回路13の機能を用いることにより、第1の実施形態に係る超音波診断装置は、コントラストハーモニックイメージング(CHI:Contrast Harmonic Imaging)等のハーモニックイメージングを実行可能である。すなわち、Bモード処理回路13は、造影剤が注入された被検体Pの反射波データから、造影剤(微小気泡、バブル)を反射源とするハーモニック成分の反射波データ(高調波データ又は分周波データ)と、被検体P内の組織を反射源とする基本波成分の反射波データ(基本波データ)とを分離する。Bモード処理回路13は、ハーモニック成分の反射波データから、造影画像データを生成するためのBモードデータを生成することができる。
また、このBモード処理回路13のフィルタ処理機能を用いることにより、第1の実施形態に係る超音波診断装置は、ティッシュハーモニックイメージング(THI:Tissue Harmonic Imaging)を実行可能である。すなわち、Bモード処理回路13は、被検体Pの反射波データから、ハーモニック成分の反射波データである高調波データ又は分周波データを分離することができる。そして、Bモード処理回路13は、ハーモニック成分の反射波データから、ノイズ成分を除去した組織画像データを生成するためのBモードデータを生成することができる。
また、CHIやTHIのハーモニックイメージングを行なう際、Bモード処理回路13は、上述したフィルタ処理を用いた方法とは異なる方法により、ハーモニック成分を抽出することができる。ハーモニックイメージングでは、振幅変調(AM:Amplitude Modulation)法や位相変調(PM:Phase Modulation)法、AM法及びPM法を組み合わせたAMPM法と呼ばれる映像法が行なわれる。AM法、PM法及びAMPM法では、同一の走査線に対して振幅や位相が異なる超音波送信を複数回行なう。これにより、送受信回路11は、各走査線で複数の反射波データを生成し、生成した反射波データを出力する。そして、Bモード処理回路13は、各走査線の複数の反射波データを、変調法に応じた加減算処理することで、ハーモニック成分を抽出する。そして、Bモード処理回路13は、ハーモニック成分の反射波データに対して包絡線検波処理等を行なって、Bモードデータを生成する。
例えば、PM法が行なわれる場合、送受信回路11は、処理回路18が設定したスキャンシーケンスにより、例えば(-1,1)のように、位相極性を反転させた同一振幅の超音波を、各走査線で2回送信させる。そして、送受信回路11は、「-1」の送信による反射波データと、「1」の送信による反射波データとを生成し、Bモード処理回路13は、これら2つの反射波データを加算する。これにより、基本波成分が除去され、2次高調波成分が主に残存した反射波データが生成される。そして、Bモード処理回路13は、この反射波データに対して包絡線検波処理等を行なって、THIのBモードデータやCHIのBモードデータを生成する。
又は、例えば、THIでは、反射波データに含まれる2次高調波成分と差音成分とを用いて映像化を行なう方法が実用化されている。差音成分を用いた映像化法では、例えば、中心周波数が「f1」の第1基本波と、中心周波数が「f1」より大きい「f2」の第2基本波とを合成した合成波形の送信超音波を、超音波プローブ1から送信させる。この合成波形は、2次高調波成分と同一の極性を持つ差音成分が発生するように、互いの位相が調整された第1基本波の波形と第2基本波の波形とを合成した波形である。送受信回路11は、合成波形の送信超音波を、位相を反転させながら、例えば、2回送信させる。かかる場合、例えば、Bモード処理回路13は、2つの反射波データを加算することで、基本波成分が除去され、差音成分及び2次高調波成分が主に残存したハーモニック成分を抽出した後、包絡線検波処理等を行なう。
ドプラ処理回路14は、バッファ12から反射波データを読み出し、読み出した反射波データを周波数解析することで、走査範囲内にある移動体のドプラ効果に基づく運動情報を推定し、推定した運動情報を示すデータ(ドプラデータ)を生成する。例えば、ドプラ処理回路14は、移動体の運動情報として、平均速度、平均分散値、平均パワー値等を複数のサンプル点それぞれで推定し、推定した運動情報を示すドプラデータを生成する。ここで、移動体とは、例えば、血流や、心壁等の組織、造影剤である。本実施形態に係るドプラ処理回路14は、血流の運動情報(血流情報)として、血流の平均速度、血流速度の分散値、血流信号のパワー値等を、複数のサンプル点それぞれで推定し、推定した血流情報を示すドプラデータを生成する。
Bモード処理回路13及びドプラ処理回路14は、2次元の反射波データ及び3次元の反射波データの両方について処理可能である。すなわち、Bモード処理回路13は、2次元の反射波データから2次元のBモードデータを生成し、3次元の反射波データから3次元のBモードデータを生成する。また、ドプラ処理回路14は、2次元の反射波データから2次元のドプラデータを生成し、3次元の反射波データから3次元のドプラデータを生成する。
上記のドプラ処理回路14の機能を用いて、本実施形態に係る超音波診断装置は、カラーフローマッピング(CFM:Color Flow Mapping)法とも呼ばれるカラードプラ法を実行可能である。CFM法では、超音波の送受信が複数の走査線上で複数回行なわれる。そして、CFM法では、同一位置のデータ列に対してMTI(Moving Target Indicator)フィルタを掛けることで、静止している組織、又は、動きの遅い組織に由来する信号(クラッタ信号)を抑制して、血流に由来する信号を抽出する。そして、CFM法では、この血流信号から血流の速度、血流の分散、血流のパワー等の血流情報を推定する。後述する画像生成回路15は、推定結果の分布を、例えば、2次元でカラー表示した超音波画像データ(カラードプラ画像データ)を生成する。そして、モニタ2は、カラードプラ画像データが示すカラードプラ画像を表示する。なお、血流情報の推定結果の分布をカラー表示した超音波画像データは、血流画像データとも称される。
MTIフィルタとしては、通常、バタワース型のIIR(Infinite Impulse Response)フィルタや、多項式回帰フィルタ(Polynomial Regression Filter)等、係数が固定されたフィルタが用いられる。一方、本実施形態に係るドプラ処理回路14は、MTIフィルタとして、入力信号に応じて係数を変化させる適応型のMTIフィルタを用いる。具体的には、本実施形態に係るドプラ処理回路14は、適応型のMTIフィルタとして、「Eigenvector Regression Filter」と呼ばれているフィルタを用いる。以下、固有ベクトルを用いた適応型MTIフィルタである「Eigenvector Regression Filter」を、「固有ベクトル型MTIフィルタ」と記載する。
固有ベクトル型MTIフィルタは、相関行列から固有ベクトルを計算し、計算した固有ベクトルから、クラッタ成分抑制処理に用いる係数を計算する。この方法は、主成分分析や、カルーネン・レーベ変換(Karhunen-Loeve transform)、固有空間法で使われている手法を応用したものである。
固有ベクトル型MTIフィルタを用いる第1の実施形態に係るドプラ処理回路14は、図1に例示するように、相関行列計算機能141と、計算機能142と、MTIフィルタ処理機能143と、推定機能144と、ブロック分割機能145と、領域合成機能146とを有する。ここで、例えば、図1に示すドプラ処理回路14の構成要素である相関行列計算機能141と、計算機能142と、MTIフィルタ処理機能143と、推定機能144と、ブロック分割機能145と、領域合成機能146とが実行する各処理機能は、コンピュータによって実行可能なプログラム(医用画像処理プログラム)の形態で内部記憶回路17内に記録されている。ドプラ処理回路14は、例えば、プロセッサにより実現される。ドプラ処理回路14は、内部記憶回路17から各プログラムを読み出し、実行することで読み出した各プログラムに対応する機能を実現する。換言すると、各プログラムを読み出した状態のドプラ処理回路14は、図1のドプラ処理回路14内に示された各機能を有することとなる。
相関行列計算機能141は、例えば、同一位置(同一サンプル点)の連続した反射波データのデータ列から、走査範囲の相関行列を計算する。相関行列計算機能141は、例えば、計算部の一例である。
計算機能142は、例えば、相関行列の固有値及び当該固有値に対応する固有ベクトルを計算する。計算機能142は、例えば、各固有値の大きさに基づいて各固有ベクトルを並べた行列のランクを低減した行列を、クラッタ成分を抑制するフィルタ行列として計算する。計算機能142は、例えば、決定部の一例である。また、フィルタ行列は、固有ベクトル型MTIフィルタに与えられるフィルタ係数の一例である。
MTIフィルタ処理機能143は、フィルタ行列を用いて、同一位置(同一サンプル点)の連続した反射波データのデータ列から、クラッタ成分が抑制され、血流に由来する血流信号が抽出されたデータ列を出力する。
推定機能144は、MTIフィルタ処理機能143が出力したデータを用いた自己相関演算等の演算を行なって、血流情報を推定し、推定した血流情報を示すドプラデータを出力する。
ブロック分割機能145は、複数の走査線で形成される走査範囲を複数の領域に分割する。領域合成機能146は、分割した領域のデータを合成する。なお、第1の実施形態に係るドプラ処理回路14が行なう具体的な処理については、後に詳述する。
画像生成回路15は、Bモード処理回路13及びドプラ処理回路14が生成したデータから超音波画像データを生成する。画像生成回路15は、Bモード処理回路13が生成した2次元のBモードデータから反射波の強度を輝度で表した2次元Bモード画像データを生成する。また、画像生成回路15は、ドプラ処理回路14が生成した2次元のドプラデータから血流情報が映像化された2次元ドプラ画像データを生成する。2次元ドプラ画像データは、速度画像データ、分散画像データ、パワー画像データ、又は、これらを組み合わせた画像データである。画像生成回路15は、ドプラ画像データとして、血流情報がカラーで表示されるカラードプラ画像データを生成したり、1つの血流情報がグレースケールで表示されるドプラ画像データを生成したりする。
ここで、画像生成回路15は、一般的には、超音波走査の走査線信号列を、テレビ等に代表されるビデオフォーマットの走査線信号列に変換(スキャンコンバート)し、表示用の超音波画像データを生成する。具体的には、画像生成回路15は、超音波プローブ1による超音波の走査形態に応じて座標変換を行なうことで、表示用の超音波画像データを生成する。また、画像生成回路15は、スキャンコンバート以外に、種々の画像処理として、例えば、スキャンコンバート後の複数の画像フレームを用いて、輝度の平均値画像を再生成する画像処理(平滑化処理)や、画像内で微分フィルタを用いる画像処理(エッジ強調処理)等を行なう。また、画像生成回路15は、超音波画像データに、種々のパラメータの文字情報、目盛り、ボディーマーク等を合成する。
すなわち、Bモードデータ及びドプラデータは、スキャンコンバート処理前の超音波画像データであり、画像生成回路15が生成するデータは、スキャンコンバート処理後の表示用の超音波画像データである。なお、Bモードデータ及びドプラデータは、生データ(Raw Data)とも呼ばれる。画像生成回路15は、スキャンコンバート処理前の2次元超音波画像データから、表示用の2次元超音波画像データを生成する。
更に、画像生成回路15は、Bモード処理回路13が生成した3次元のBモードデータに対して座標変換を行なうことで、3次元Bモード画像データを生成する。また、画像生成回路15は、ドプラ処理回路14が生成した3次元のドプラデータに対して座標変換を行なうことで、3次元ドプラ画像データを生成する。画像生成回路15は、「3次元のBモード画像データ及び3次元ドプラ画像データ」を「3次元超音波画像データ(ボリュームデータ)」として生成する。
更に、画像生成回路15は、ボリュームデータをモニタ2にて表示するための各種の2次元画像データを生成するために、ボリュームデータに対してレンダリング処理を行なう。画像生成回路15が行なうレンダリング処理としては、例えば、断面再構成法(MPR:Multi Planer Reconstruction)を行なってボリュームデータからMPR画像データを生成する処理がある。また、画像生成回路15が行なうレンダリング処理としては、例えば、3次元の情報を反映した2次元画像データを生成するボリュームレンダリング(VR:Volume Rendering)処理がある。なお、画像生成回路15は、画像生成部の一例である。
画像メモリ16は、画像生成回路15が生成した表示用の画像データを記憶するメモリである。また、画像メモリ16は、Bモード処理回路13やドプラ処理回路14が生成したデータを記憶することも可能である。画像メモリ16が記憶するBモードデータやドプラデータは、例えば、診断の後に操作者が呼び出すことが可能となっており、画像生成回路15を経由して表示用の超音波画像データとなる。また、画像メモリ16は、送受信回路11が出力した反射波データを記憶することも可能である。例えば、画像メモリ16は、RAM、フラッシュメモリ等の半導体メモリ素子、ハードディスク又は光ディスクによって実現される。
内部記憶回路17は、超音波送受信、画像処理及び表示処理を行なうための制御プログラムや、診断情報(例えば、患者ID、医師の所見等)や、診断プロトコルや各種ボディーマーク等の各種データを記憶する。また、内部記憶回路17は、必要に応じて、画像メモリ16が記憶する画像データの保管等にも使用される。また、内部記憶回路17が記憶するデータは、図示しないインターフェースを経由して、外部装置へ転送することができる。また、内部記憶回路17は、外部装置から図示しないインターフェースを経由して転送されたデータを記憶することも可能である。例えば、内部記憶回路17は、フラッシュメモリ等の半導体メモリ素子、ハードディスク又は光ディスクによって実現される。
処理回路18は、超音波診断装置の処理全体を制御する。具体的には、処理回路18は、入力装置3を介して操作者から入力された各種設定要求や、内部記憶回路17から読込んだ各種制御プログラム及び各種データに基づき、送受信回路11、Bモード処理回路13、ドプラ処理回路14及び画像生成回路15の処理を制御する。例えば、処理回路18は、送受信回路11を介して超音波プローブ1を制御することで、超音波走査の制御を行なう。通常、CFM法では、血流像データであるカラードプラ画像データとともに、組織像データであるBモード画像データを表示する。かかる表示を行なうため、処理回路18は、第1走査範囲内の血流情報を取得する第1超音波走査を超音波プローブ1に実行させる。第1超音波走査は、例えば、ドプラモードでカラードプラ画像データを収集するための超音波走査である。また、第1超音波走査は、例えば、1フレームに同一走査線を1回走査する超音波走査である。
また、処理回路18は、第1超音波走査とともに、第2走査範囲内の組織形状の情報を取得する第2超音波走査を超音波プローブ1に実行させる。第2超音波走査は、例えば、BモードでBモード画像データを収集するための超音波走査である。
処理回路18は、送受信回路11を介して超音波プローブ1を制御することで、第1超音波走査及び第2超音波走査を実行させる。なお、第1走査範囲と第2走査範囲は、同じ範囲であっても、第1走査範囲が第2走査範囲より小さい範囲であっても、第2走査範囲が第1走査範囲より小さい範囲であっても良い。
また、処理回路18は、画像メモリ16や内部記憶回路17が記憶する表示用の超音波画像データが示す超音波画像をモニタ2にて表示するように制御する。なお、装置本体10に内蔵される送受信回路11等は、集積回路などのハードウエアで構成されることもあるが、ソフトウエア的にモジュール化されたプログラムである場合もある。処理回路18は、例えば、プロセッサにより実現される。処理回路18は、制御部の一例である。
以上、第1の実施形態に係る超音波診断装置の全体構成について説明した。ここで、一般的な固有ベクトル型MTIフィルタを用いる超音波診断装置を、第1の比較例として説明する。第1の比較例では、固有ベクトル型MTIフィルタは、振幅が比較的大きく、所定の範囲内で一様な動きをする組織成分(クラッタ成分)を顕著に抑制するという特徴を有する。
ここで、第1の比較例において、超音波走査の対象が肝臓であり、強反射体であるとともに周期的に移動する心臓壁や、強反射体であるとともに心臓の周期的な動きの影響を受けて周期的に移動する横隔膜等がサイドローブの反射源として画像範囲外に存在する場合について説明する。この場合には、肝臓の近傍に存在する心臓壁や横隔膜等がサイドローブを反射する。反射されたサイドローブの振幅は、比較的大きくない。また、反射されたサイドローブの位相は位置によって変化するので、位相を含めた動き(又は位相変位)は、一様ではない。また、反射されたサイドローブの振幅及び位相は、血流の振幅及び位相に近い。このため、固有ベクトル型MTIフィルタによるクラッタ成分の抑制が不十分となり、サイドローブが血流信号として固有ベクトル型MTIフィルタにより抽出されることがある。このような現象は、定常的に発生するものではなく、所定の心時相、例えば、心臓壁が最も急峻に変化する収縮期で発生する。この結果、第1の比較例では、抑制しきれないクラッタ成分が、特定の心時相に対応する超音波画像データの特定の領域に周期的に含まれてしまうことがある。そのため、モニタに表示される動画像では、定常的には血流が明瞭に描出されているが、1心拍に1回の特定の時相において、短い期間だけ大きなクラッタが特定の領域に発生する。
図2A及び図2Bは、第1の比較例に係る超音波診断装置により生成された血流画像データが示す血流画像の一例を示す図である。なお、図2A及び図2Bに示す血流画像101,102は、被検体の肝臓の血流情報が描出された画像である。血流画像101,102は、パワー画像データから生成されたパワー画像である。第1の比較例に係る超音波診断装置は、画像化対象範囲(走査範囲)を複数の領域に分割し、領域ごとに相関行列を計算し、領域毎に異なるMTIフィルタ行列を掛けることにより、血流画像101を示す血流画像データを生成する。また、第1の比較例に係る超音波診断装置は、同様の方法により、血流画像102を示す血流画像データを生成する。具体的には、例えば、これらの血流画像データは、例えば、32パケット長のデータ列を領域毎にアンサンブル平均して得られた相関行列の固有値及び固有ベクトルから、固有ベクトル型MTIフィルタのフィルタ係数(MTIフィルタ係数)を計算することにより生成される。ここで、血流画像101に対応する時相(第1の時相)と、血流画像102に対応する時相(第2の時相)とは異なる。
第1の時相の血流画像101では、画像全体的にクラッタが比較的少ない。しかしながら、第2の時相の血流画像102では、画像範囲外にサイドローブの反射源として存在する心臓壁や横隔膜等が大きく移動しているため、図2Bにおいて血流画像102の右側の部分にクラッタ102aが現れている。
図3A及び図3Cは、第1の比較例に係る超音波診断装置により生成された血流画像データが示す血流画像の他の例を示す図である。なお、図3A及び図3Cに示す血流画像103,104も、被検体の肝臓の血流情報が描出された画像であり、パワー画像データから生成されたパワー画像である。第1の比較例に係る超音波診断装置は、図2Aに示す血流画像101を示す血流画像データを生成する方法と同様の方法で、図3Aに示す血流画像103を示す血流画像データ及び図3Cに示す血流画像104を示す血流画像データを生成する。血流画像103,104は、コンベックスプローブにより収集された反射波信号に基づく血流画像データから生成された画像であり、座標変換されずにモニタに表示された画像である。
ここで、上述したように、血流画像103,104を生成する際に、相関行列が複数の領域のそれぞれにおいて計算されるが、図3A及び図3Cに示す5つの符号111~115のそれぞれは、複数の領域のうちの5つの領域のそれぞれを指す。
血流画像103に対応する時相(第3の時相)と、血流画像104に対応する時相(第4の時相)とは異なる。第3の時相では画像範囲外にサイドローブの反射源として存在する心臓壁や横隔膜等の動きが比較的小さいため、血流画像103全体的にクラッタが比較的少ない。しかしながら、第4の時相では心臓壁や横隔膜等が大きく移動しているため、図3Cにおいて血流画像104の右側の部分に多くのクラッタ104aが含まれている。
図3Bは、図3Aに示す血流画像を示す血流画像データを生成する際に計算された各領域における複数の固有値の大きさの一例を示すグラフである。図3Dは、図3Cに示す血流画像を示す血流画像データを生成する際に計算された各領域における複数の固有値の大きさの一例を示すグラフである。
図3B及び図3Dに示すグラフにおいて、横軸は、各領域における全固有値を、固有値の大きさについての降順で並べた場合の順位(順番)を示す。また、縦軸は、全ての領域における全ての相関行列をアンサンブル平均することにより得られた相関行列の固有値のうち最も大きい固有値を0dBとした場合の各領域の各順位の固有値のデジベル単位の大きさを示す。
図3B及び図3Dには、5つの領域111~115のそれぞれについてのグラフが示されている。具体的には、図3B及び図3Dにおいて、グラフ111aは、領域111についてのグラフであり、グラフ112aは、領域112についてのグラフであり、グラフ113aは、領域113についてのグラフであり、グラフ114aは、領域114についてのグラフであり、グラフ115aは、領域115についてのグラフである。また、図3B及び図3Dにおいて、グラフ116aは、全ての領域における全ての相関行列をアンサンブル平均することにより得られた相関行列の全固有値についてのグラフである。
また、図3B及び図3Dには、各領域において、32個の固有値が計算された場合のグラフ111a~116aが示されている。すなわち、グラフ111a~116aのそれぞれは、順位が第1位の固有値(順番が1番目の固有値)から順位が第32位の固有値(順番が32番目の固有値)までの各固有値の大きさを示す。
ここで、第24位の固有値に注目して説明する。図3Bでは、領域112に対応するグラフ112aが示す第24位の固有値の大きさが、他の領域111,113~115に対応するグラフ111a,113a~115aが示す第24位の固有値よりも大きい。そこで、図3Aに示す血流画像103を確認すると、領域112は、血流情報が描出された領域であることが分かる。
一方、図3Dでは、領域115に対応するグラフ115aが示す第24位の固有値の大きさが、他の領域111~114に対応するグラフ111a~114aが示す第24位の固有値よりも大きい。そこで、図3Cに示す血流画像104を確認すると、領域115には、クラッタが含まれていることが分かる。
図4は、各領域における第24位の固有値の大きさの時間変化の一例を示すグラフである。図4には、5つの領域111~115のそれぞれにおける第24位の固有値の大きさの5秒間の時間変化を示す5つのグラフ111b~115bのそれぞれが示されている。
具体的には、図4において、グラフ111bは、領域111における第24位の固有値の大きさの時間変化を示すグラフであり、グラフ112bは、領域112における第24位の固有値の大きさの時間変化を示すグラフである。また、グラフ113bは、領域113における第24位の固有値の大きさの時間変化を示すグラフであり、グラフ114bは、領域114における第24位の固有値の大きさの時間変化を示すグラフであり、グラフ115bは、領域115における第24位の固有値の大きさの時間変化を示すグラフである。また、図4において、グラフ116bは、全ての領域における全ての相関行列をアンサンブル平均することにより得られた相関行列の複数の固有値のうち、第24位の固有値についての大きさの時間変化を示すグラフである。
グラフにおいて、横軸は、時間(秒(sec))を示す。横軸では、現在からの過去の時間が負の値で示されている。例えば、時間「0(sec)」は、現在の時間(現時点)を示す。また、例えば、時間「-0.5(sec)」は、現在から0.5秒前の過去の時間(時点)を示す。縦軸は、全ての領域における全ての相関行列をアンサンブル平均することにより得られた相関行列の複数の固有値のうち最も大きい固有値を0dBとした場合の各領域の第24位の固有値のデジベル単位の大きさを示す。
例えば、図3Aに示す血流画像103に対応する第3の時相は、図4において時間「-3.3(sec)」に対応する。また、図3Cに示す血流画像104に対応する第4時相は、図4において、時間「-3.07(sec)」に対応する。すなわち、血流画像103を示す血流画像データは、時間「-3.3(sec)」に収集された反射波信号に基づくデータであり、血流画像104を示す血流画像データは、時間「-3.07(sec)」に収集された反射波信号に基づくデータである。
領域115に注目して説明する。血流画像103の領域115には、クラッタがほとんど含まれていないが、血流画像104の領域115には、大きなクラッタが含まれている。この様子が、領域115に対応する図4に示すグラフ115bに表れている。例えば、固有値が比較的大きい時相(例えば、時間「-3.07(sec)」)と固有値が比較的小さい時相(例えば、時間「-3.3(sec)」)があるグラフ115bは、固有値が比較的大きい時相においてクラッタが発生し、固有値が比較的小さい時相においてクラッタがほとんど発生しないことを表す。固有値が比較的大きい時相は、例えば、心臓壁の動きが比較的大きいため、クラッタが発生する時相と考えられる。また、固有値が比較的小さい時相は、例えば、心臓壁の動きが比較的小さいため、クラッタがほとんど発生しない時相と考えられる。
また、グラフ115bは、他のグラフ111b~114b,116bと比較して、固有値の変動が比較的大きいグラフである。また、グラフ111b~116bは、被検体Pの心周期に同期して周期的に変動している。したがって、グラフ115bでは、固有値が比較的大きい時相と比較的小さい時相とが周期的に表れる。よって、グラフ115bに対応する領域115には、周期的にクラッタが含まれる。
ここで、固有値の順位と固有値の大きさ(例えば、デジベル値)とから血流か否かを判定する方法(上記特許文献1及び上記非特許文献2参照)がある。しかしながら、図3A~図3D及び図4を参照して説明した以上のことから、一つの時相の固有値の順位と大きさに基づいて血流かクラッタかを正確に分別することは困難である。これは、生体での超音波の減衰や周囲の血流の多寡、組織の反射強度によって、固有値の順位及び大きさは様々に変化するからである。特に、サイドローブを起因とするクラッタは、固有値からでは血流かクラッタかの分別が困難である。
また、超音波診断装置が、心臓壁の移動速度と分散をTDI(Tissue Doppler Imaging)法により計算し、心臓壁の移動速度と分散からMTIフィルタのカットオフ周波数を適応的に変化させる方法を実行することも考えられる。このような方法を実行する超音波診断装置を、第2の比較例に係る超音波診断装置として説明する。しかしながら、第2の比較例に係る超音波診断装置において、TDI法で観測されるのは比較的振幅の大きい組織の動きである。一方で、心臓壁や横隔膜等により反射されたサイドローブの振幅は、比較的大きくない。このため、第2の比較例に係る超音波診断装置が、上述したカットオフ周波数を適応的に変化させる方法を実行する場合であっても、クラッタ成分の抑制が不十分となることがある。
そこで、第1の実施形態に係る超音波診断装置は、心臓壁や横隔膜等の強反射体にサイドローブが反射されることによって発生するクラッタ成分を抑制することが可能なMTIフィルタを生成することができるように、以下に説明する処理を実行する。第1の実施形態に係る超音波診断装置は、1心拍以上の時間幅で、各中ブロック(領域)の所定の順位の固有値の大きさの時間変化を測定し、その時間幅での固有値の最小値と現在の固有値の大きさとの差に基づいて、固有ベクトル型MTIフィルタの特性を決定するランクカット数を決定する。以下の説明では、固有ベクトル型MTIフィルタを、単に、MTIフィルタと表記する。
図5は、第1の実施形態に係る超音波診断装置による処理手順を示すフローチャートである。
図5に示すステップS1からステップS3は、ブロック分割機能145に対応するステップである。ステップS1からステップS3は、ドプラ処理回路14が内部記憶回路17からブロック分割機能145に対応する所定のプログラムを呼び出し実行することにより、ブロック分割機能145が実現されるステップである。ステップS1では、ブロック分割機能145は、処理対象となるスキャンフレーム数分の反射波データが収集されたか否かを判定する(ステップS1)。ここで、ブロック分割機能145は、スキャンフレーム数分の反射波データが収集されていないと判定した場合(ステップS1;No)、再びステップS1の判定処理を実行する。
一方、ブロック分割機能145は、スキャンフレーム数分の反射波データが収集されたと判定した場合(ステップS1;Yes)、ステップS2に移行する。
なお、第1の実施形態では、ステップS2以降の処理は、バッファ12に、新たに1フレーム分の反射波データが格納される度に実行される。また、第1の実施形態では、1フレーム毎に、リアルタイムで、ステップS2以降の各処理が実行される。
ステップS2では、ブロック分割機能145は、全空間地点の全パケットデータをメモリにコピーする。ここで全空間地点の数をN、パケットサイズをLとする。例えば、ブロック分割機能145は、バッファ12内の入力データをドプラ処理回路14内に設けられたメモリにコピーする。ブロック分割機能145は、呼び出すパケット毎にこのコピーする操作を実行する。
ステップS3では、ブロック分割機能145は、画像化対象範囲を複数の領域に分割する。なお、画像化対象範囲は、例えば、超音波の走査範囲である。図6は、第1の実施形態を説明するための図である。図6では、全画像化対象範囲のうち一部の領域を図示している。また、図6では、同一の領域について、(A)から(I)の9種類の分割パターンで領域を分割した場合を図示している。
例えば、ブロック分割機能145は、図6中の小さい升目を最小の分割単位とするように画像化対象範囲を複数の領域に分割する。なお、この最小の分割単位のことを小ブロックと記載する。図6の(A)から(I)において、小ブロックの分割パターンは同一となる。一例をあげると、ブロック分割機能145は、各小ブロックには8ラスタ、8サンプルが含まれるように分割する。なお、ブロック分割機能145は、自由度(≧パケット長L)が確保できるだけの空間サンプル数を有するように空間領域を小ブロックに分割する。
そして、ブロック分割機能145は、複数の小ブロックを用いて中ブロックを構成する。例えば、ブロック分割機能145は、図6中の太線で囲った縦横それぞれ4個ずつの小ブロックを集めて中ブロックとする。すなわち、中ブロックは、4×4の小ブロックから構成される。一例をあげると、ブロック分割機能145は、図6の(A)から(I)に示すように、中ブロックを構成する。
図6の(A)の黒い太枠で囲った中ブロックを例に取って説明する。(A)に示す中ブロックは、隣接する(B)に示す中ブロックから(I)に示す中ブロックまでの8つの中ブロックと一部が重なっている。例えば、(A)の黒塗りで示した小ブロックは(B),(D),(F)に示す中ブロックにも重複して含まれている。
言い換えると、ブロック分割機能145は、各中ブロックの一部が他の中ブロックの一部と重なり合うように、画像化対象範囲を複数の中ブロックに分割する。
ステップS4からステップS6は、相関行列計算機能141に対応するステップである。ステップS4からステップS6は、ドプラ処理回路14が内部記憶回路17から相関行列計算機能141に対応する所定のプログラムを呼び出し実行することにより、相関行列計算機能141が実現されるステップである。ステップS4では、相関行列計算機能141は、各小ブロック内で相関行列R1を計算する。ある地点iでのパケットデータ列ベクトルをxiとすると、相関行列R1は以下の式(1)で表わされる。ここでiはある地点の位置を表す。例えば、2次元スキャンの場合にはx,zの位置を、3次元スキャンの場合にはx,y,zの位置を1つの添字iで表わしている。Nは計算する地点の数である。Hは複素共役転置行列(エルミート転置行列)を示す。パケット長をLとすると相関行列R1はL×Lの行列となる。
なお、相関行列計算機能141は、1フレームに同一走査線を1回走査する上述した第1超音波走査により複数フレームにわたって収集されたパケットデータ列ベクトルxiから、相関行列R1を計算してもよい。ここで、このようなパケットデータ列ベクトルxiは、同一位置の反射波データのデータ列である。
ステップS5では、相関行列計算機能141は、すべての小ブロックの相関行列を加算平均して全体の相関行列Rallを計算する。相関行列Rallは、Mallをすべての小ブロックの数とし、Mall個の相関行列R1のうちm番目の相関行列R1をR1,mとすると、以下の式(2)で表わされる。
ここで、相関行列計算機能141は、複数フレームにわたって収集された上述のパケットデータ列ベクトルxiに基づく相関行列R1,mから、相関行列Rallを計算してもよい。すなわち、相関行列計算機能141は、複数フレームにわたって収集されたパケットデータ列ベクトルxiから、相関行列Rallを計算してもよい。
ステップS6では、相関行列計算機能141は、中ブロックの相関行列R2を計算する。例えば、相関行列計算機能141は、1つの中ブロックに含まれる小ブロックの数をM2とし、M2個の相関行列R1のうちm番目の相関行列R1をR1,mとすると、以下の式(3)を用いて中ブロックの相関行列R2を計算する。
ここで、相関行列計算機能141は、複数フレームにわたって収集された上述のパケットデータ列ベクトルxiに基づく相関行列R1,mから、相関行列R2を計算してもよい。すなわち、相関行列計算機能141は、複数フレームにわたって収集されたパケットデータ列ベクトルxiから、相関行列R2を計算してもよい。この場合、相関行列計算機能141は、複数の走査線により形成される画像化対象範囲での超音波走査により収集された同一位置の反射波データのデータ列から、画像化対象範囲を分割した複数の中ブロックのそれぞれの相関行列R2を計算する。ここで、中ブロックは、画像化対象範囲を分割した領域の一例である。
相関行列計算機能141は、4×4の小ブロックから構成される各中ブロックの相関行列R2を計算する。より具体的には、相関行列計算機能141は、小ブロックで計算したM2個の相関行列R1を加算することにより得られた相関行列を、M2で除算する。すなわち、相関行列計算機能141は、中ブロック内の小ブロックの相関行列R1を加算平均して中ブロックの相関行列R2を計算する。
ステップS7からステップS13は、計算機能142に対応するステップである。ステップS7からステップS13は、ドプラ処理回路14が内部記憶回路17から計算機能142に対応する所定のプログラムを呼び出し実行することにより、計算機能142が実現されるステップである。ステップS7では、計算機能142は、全体の相関行列Rallから最大固有値λ1を計算する。より具体的には、計算機能142は、以下の式(4)を用いて、相関行列Rallを固有値分解することにより、L個の固有値λ(λ1,λ2,・・・λL)と、各固有値λに対応する各固有ベクトルを計算する。
上記式(4)において、Vは固有ベクトルを列ベクトルに持つ行列であり、Dは固有値を対角要素に持つ対角行列である。L個の固有値λと、各固有値λに対応する各固有ベクトルは降順に並んでいるとする。行列Vは以下の式(5)により表わされ、行列Dは以下の式(6)により表わされる。なお、計算機能142は、L個の固有値λの中で最大固有値をλ1とする。
ステップS8では、計算機能142は、中ブロックの相関行列R2から固有値及び固有ベクトルを計算する。ここで、計算機能142は、ステップS7で全体の相関行列Rallから固有値及び固有ベクトルを計算した方法と同様の方法で、中ブロックの相関行列R2から固有値及び固有ベクトルを計算する。
ステップS9では、計算機能142は、固有値を規格化する。例えば、計算機能142は、中ブロック毎に相関行列R2の固有値λを全体の相関行列Rallの全パワー(相関行列Rallのすべての固有値の和)で除算することにより、相関行列R2の固有値λを規格化する。これによりゲインの影響を除去することができる。なお、計算機能142は、全パワーで規格化する代わりに相関行列Rallの最大固有値λ1で除算することにより、相関行列R2の固有値λを規格化してもよい。このように、画像全体の相関行列Rallの最大固有値λ1で規格化することで、ゲイン依存や領域分割による生体組織差による影響を除去することができる。規格化された固有値λがステップS10以降の処理で用いられる。
ステップS10では、計算機能142は、中ブロック毎に、現在から、現在から1.2秒前までの所定時間においてステップS8で計算された相関行列R2の第S位の複数の固有値のうち、第S位の固有値の最小値を特定する。ここで、「S」は、1以上L以下の自然数である。また、「L」は、パケット長である。なお、1.2秒は、例えば、被検体Pの心拍数が50bpm(beats per minute)である場合の1心拍の時間である。なお、所定時間は、1.2秒に限られず、1心拍の時間以上であればよい。ただし、ここでいう所定時間は、長くなるほど、手に保持される超音波プローブ1の動きがノイズとして超音波画像に含まれてやすくなる。このため、所定時間は、1心拍程度の時間が好ましい。
また、ここでいう第S位の固有値とは、例えば、相関行列R2の全固有値を、固有値の大きさについての降順で並べた場合のS番目の固有値である。なお、以下、超音波診断装置が、血流の存在を示す固有値として、第1位から第32位までの固有値のうち、第24位の固有値を用いる場合を例に挙げて説明する。すなわち、「S」が「24」であり、「L」が「32」である場合を例に挙げて説明する。しかしながら、超音波診断装置は、血流の存在を示す固有値であれば、第24位以外の順位の固有値を用いてもよい。例えばこの値はパケット長Lや生体の部位やユーザからの入力によって決定される。
図7は、第1の実施形態に係る計算機能が実行する処理の一例を説明するための図である。図7には、現在(時間「0(sec)」)から、現在から5秒前の時間(「-5.0(sec)」)までの5秒間における第S位の固有値の時間変化を示すグラフ21が示されている。グラフ21において、横軸は、時間(秒(sec))を示す。横軸では、現在からの過去の時間が負の値で示されている。縦軸は、相関行列Rallの複数の固有値のうち最も大きい固有値を0dBとした場合の第S位の固有値のデジベル単位の大きさを示す。また、図7には、時間「0(sec)」~「-1.2(sec)」の時間範囲、時間「-1.2(sec)」~「-2.4(sec)」の時間範囲、及び、時間「-2.4(sec)」~「-3.6(sec)」の時間範囲の各時間範囲における第S位の固有値の最小値を示す線21aが示されている。
例えば、ステップS10では、計算機能142は、図7のグラフ21及び線21aが示すように、現在から、現在から1.2秒前までの所定時間において第S位のデジベル値の固有値の最小値「-71.0(dB)」を特定する。このようにして、ステップS10では、計算機能142は、中ブロック毎に、過去1心拍程度の時間における第S位のデジベル値の固有値の最小値を特定する。上述したように、計算機能142は、中ブロック毎に、相関行列R2の複数の固有値のうちの1つの固有値について、現在から、現在から所定時間(例えば、1.2秒)前までの所定時間の間での最小値を特定する。
そして、ステップS11では、計算機能142は、中ブロック毎に、最終的なランクカット数kfを決定する。なお、ランクカット数とは、例えば、後述するMTIフィルタ行列Wを計算する際に、低減されるランク数のことを指す。ステップS11の処理の一例について説明する。ステップS11では、まず、計算機能142は、中ブロック毎に、現在の第S位のデシベル値の固有値とステップS10で特定した第S位のデシベル値の固有値の最小値との差を計算する。以下の説明では、かかる差を、「固有値の差」と表記する場合がある。具体的には、計算機能142は、現在の第S位のデシベル値の固有値から、ステップS10で特定した第S位のデシベル値の固有値の最小値を減算することにより、デシベル値での固有値の差を計算する。図8は、第1の実施形態に係る計算機能が実行する処理の一例を説明するための図である。図8には、グラフ22が示されている。グラフ22は、図7に示すグラフ21が示す各時間におけるデシベル値の固有値の大きさから、上述した各時間範囲における第S位のデシベル値の固有値の最小値を減算することにより得られるグラフである。
例えば、計算機能142は、図7のグラフ21が示す現在の第S位のデジベル値の固有値「-68.2(dB)」から、ステップS10で特定した第S位のデジベル値の固有値の最小値「-71.0(dB)」を減算する。これにより、図8に示すように、計算機能142は、デジベル値の固有値「-68.2(dB)」とデジベル値の固有値の最小値「-71.0(dB)」とのデジベル値での差「2.8(dB)」を計算する。
そして、ステップS11では、計算機能142は、所定の固有値の大きさ(例えば、3(dB))で、固有値の差を離散化する。例えば、計算機能142は、固有値の差を所定の固有値の大きさで除することにより得られる商(数値)を求める。計算機能142は、かかる商を、最終的なランクカット数kfを決定する際に加味する。
図9は、第1の実施形態に係る計算機能が実行する処理の一例を説明するための図である。図9には、グラフ23が示されている。グラフ23は、図8のグラフ22が示す各時間における上述した固有値の差を、所定のデジベル値の固有値の大きさ「3(dB)」で離散化することにより得られたグラフである。なお、所定のデジベル値の固有値の大きさは、「3(dB)」に限られない。グラフ23では、横軸は、時間(秒(sec))を示す。横軸では、現在からの過去の時間が負の値で示されている。縦軸は、最終的なランクカット数kfを決定する際に加味される商(数値)を示す。
例えば、計算機能142は、図8のグラフ22が示す現在に対応する固有値の差「2.8(dB)」を、所定のデジベル値の固有値の大きさ「3(dB)」で離散化することにより、図9のグラフ23が示すように、数値「0」を計算する。すなわち、計算機能142は、「2.8(dB)」を「3(dB)」で除することにより得られる商「0」を計算する。
また、ステップS11では、計算機能142は、固有値から仮のランクカット数cを決定する。例えば、計算機能142は、中ブロック毎の仮のランクカット数cを、特許文献1(特開2014-158698号公報)に記載されている方法と同様の方法により決定する。例えば、計算機能142は、中ブロック毎に、相関行列R2の固有値の大きさに応じて、仮のランクカット数cを決定する。また、仮のランクカット数cは、予め設定された値、又は、操作者により指定された値であってもよい。
そして、計算機能142は、仮のランクカット数cに、上述した商(数値)を加算して、最終的なランクカット数kfを決定する。例えば、計算機能142は、仮のランクカット数cに、図9の例において計算された商(0)を加算して、最終的なランクカット数kf(c+0)を決定する。上述したように、計算機能142は、所定時間(例えば、上述した1.2秒)の間の相関行列R2の固有値に基づいて、MTIフィルタの特性を決定するランクカット数kfを決定する。より具体的には、計算機能142は、デシベル値の固有値の差に基づいて、ランクカット数kfを決定する。
ステップS12では、計算機能142は、ランクカット数kf(kf≦L)から、ランク数kの固有ベクトル行列Vk(L×k行列)を求める。ここで、固有ベクトル行列Vkは、以下の式(7)で表わされる。また、例えば、kf=kである。
ステップS13において、計算機能142は、以下の式(8)により、固有ベクトル行列VkからMTIフィルタ行列Wを計算する。なお、式(8)において、IはL×Lの単位行列である。
上述したように、計算機能142は、所定時間(例えば、上述した1.2秒)の間の相関行列R2の固有値に基づいて、MTIフィルタに与えられるフィルタ係数としてMTIフィルタ行列Wを決定する。また、計算機能142は、所定時間の間の相関行列R2の固有値に基づいて、MTIフィルタの特性を決定するランクカット数kfを決定し、決定したランクカット数kfに基づいてMTIフィルタ行列Wを決定する。また、計算機能142は、所定時間の間の固有値の差に基づいて、MTIフィルタ行列Wを決定する。
また、上述したように、計算機能142は、複数の中ブロックのそれぞれの相関行列R2の固有値に基づいて、複数の中ブロックのそれぞれに対応するMTIフィルタ行列Wを決定する。また、計算機能142は、被検体の1心拍分の時間の間の相関行列R2の固有値に基づいて、MTIフィルタ行列Wを決定する。
ステップS14は、MTIフィルタ処理機能143に対応するステップである。ステップS14は、ドプラ処理回路14が内部記憶回路17からMTIフィルタ処理機能143に対応する所定のプログラムを呼び出し実行することにより、MTIフィルタ処理機能143が実現されるステップである。ステップS14では、MTIフィルタ処理機能143は、中ブロック毎に、以下に示す式(9)を用いて各点のパケット列ベクトルデータxiにMTIフィルタを施す。すなわち、MTIフィルタ処理機能143は、位置iのパケット列ベクトルデータxiを入力データとし、入力データとMTIフィルタ行列Wとから、位置iの出力データであるパケット列ベクトルデータyiを、以下に示す式(9)により中ブロック毎に計算する。パケット列ベクトルデータyiの要素の個数は、「L」である。
例えば、MTIフィルタ処理機能143は、図6の(A)に示す4×4の小ブロックから構成される中ブロックのデータに、この図6の(A)に示す中ブロックのデータから計算したMTIフィルタを施す。
ステップS15及びステップS16は、推定機能144に対応するステップである。ステップS15及びステップS16は、ドプラ処理回路14が内部記憶回路17から推定機能144に対応する所定のプログラムを呼び出し実行することにより、推定機能144が実現されるステップである。ステップS15及びステップS16では、推定機能144は、データ列から統計的演算を行って移動体情報を推定する。例えば、推定機能144は、複数の走査線で形成される走査範囲を複数の領域に分割して得られた各領域から計算された統計的性質を用いて各領域における移動体情報を推定する。
より具体的には、推定機能144は、ステップS15では、ステップS14でMTIフィルタが施された中ブロックのデータから、以下の式(10)を用いて、位置iにおけるパワーPiを推定する。ここで、推定機能144は、対数圧縮前の値としてパワーPiを推定する。なお、式(10)において、jは列ベクトルの要素番号を示す指標である。また、yi,jは、パケット列ベクトルデータyiのj番目の要素を示す。
また、推定機能144は、ステップS16では、ステップS14でMTIフィルタが施された中ブロックのデータから、以下の式(11)を用いて、位置iの速度Viを推定する。なお、式(11)において、「angle」は複素数の偏角を-πからπの範囲で出力する関数である。また、式(11)において、アスタリスク「*」は、複素共役を示す。
なお、推定機能144は、ステップS15とステップS16の処理順序を入れ替えてもよく、また、同時に実行してもよい。
ステップS17は、領域合成機能146に対応するステップである。ステップS17は、ドプラ処理回路14が内部記憶回路17から領域合成機能146に対応する所定のプログラムを呼び出し実行することにより、領域合成機能146が実現されるステップである。ステップS17では、領域合成機能146は、重み付け加算する。例えば、領域合成機能146は、パワーの領域合成と速度の領域合成とを行う。例えば、領域合成機能146は、重み係数を用いたバイリニア法にて位置に応じた重み係数を乗算して同じ位置のデータを加算することで画素を補間する。
ステップS18は、画像生成回路15により実現されるステップである。ステップS18では、画像生成回路15は、移動体情報から血流画像データ(カラードプラ画像データ)を生成する。例えば、画像生成回路15は、ステップS15において推定したパワーPiを対数圧縮して血流画像データ(パワー画像データ)を生成する。また、画像生成回路15は、ステップS16において推定した速度Viに基づく血流画像データ(速度画像データ)を生成する。
ステップS19は、処理回路18により実現されるステップである。ステップS19では、処理回路18は、モニタ2に、血流画像データが示す血流画像を表示させ、処理を終了する。
図10A及び図10Bは、第1の実施形態に係る超音波診断装置が図5に示す処理を実行することにより生成した血流画像データが示す血流画像の一例を示す図である。
図10A及び図10Bに示す血流画像30,40は、被検体Pの肝臓の血流情報が描出された画像である。血流画像30,40は、パワー画像データから生成されたパワー画像である。例えば、血流画像30,40は、コンベックスプローブである超音波プローブ1により収集された反射波信号に基づく血流画像データから生成された画像であり、座標変換されずにモニタ2に表示された画像である。また、血流画像30,40の範囲外には、サイドローブの反射源として心臓壁や横隔膜等が存在する。
ここで、血流画像30,40を生成する際に、複数の中ブロックのそれぞれにおいて相関行列R2が計算されるが、図10A及び図10Bに示す5つの符号31~35のそれぞれは、複数の中ブロックのうちの5つの中ブロックのそれぞれを指す。
図10Cは、各中ブロックにおける第S位の固有値の大きさの時間変化の一例を示すグラフである。図10Cには、5つの領域31~35のそれぞれにおける第S位の固有値の大きさの5秒間の時間変化を示す5つのグラフ31b~35bのそれぞれが示されている。
具体的には、図10Cにおいて、グラフ31bは、領域31における第S位の固有値の大きさの時間変化を示すグラフであり、グラフ32bは、領域32における第S位の固有値の大きさの時間変化を示すグラフである。また、グラフ33bは、領域33における第S位の固有値の大きさの時間変化を示すグラフであり、グラフ34bは、領域34における第S位の固有値の大きさの時間変化を示すグラフであり、グラフ35bは、領域35における第S位の固有値の大きさの時間変化を示すグラフである。また、図10Cにおいて、グラフ36bは、全ての中ブロックにおける全ての相関行列R2をアンサンブル平均することにより得られた相関行列の複数の固有値のうち、第S位の固有値についての大きさの時間変化を示すグラフである。
グラフにおいて、横軸は、時間(秒(sec))を示す。横軸では、現在からの過去の時間が負の値で示されている。例えば、時間「0(sec)」は、現在の時間(時点)を示す。また、例えば、時間「-0.5(sec)」は、現在から0.5秒前の過去の時間(時点)を示す。縦軸は、全ての中ブロックにおける全ての相関行列R2をアンサンブル平均することにより得られた相関行列の複数の固有値のうち最も大きい固有値を0dBとした場合の各中ブロックの第S位の固有値のデジベル単位の大きさを示す。
例えば、図10Aに示す血流画像30に対応する第5の時相は、図10Cにおいて時間「-3.3(sec)」に対応する。また、図10Bに示す血流画像40に対応する第6時相は、図4において、時間「-3.07(sec)」に対応する。すなわち、血流画像30を示す血流画像データは、時間「-3.3(sec)」に生成されたデータであり、血流画像40を示す血流画像データは、時間「-3.07(sec)」に生成されたデータである。
グラフ31b~36bは、被検体Pの心周期に同期して周期的に変動している。また、グラフ35bでは、固有値が比較的大きい第6の時相と比較的小さい第5の時相とが周期的に表れる。グラフ35bに示すように、固有値が比較的大きい第6の時相(例えば、時間「-3.07(sec)」)では、例えば、心臓壁の動きが比較的大きいと考えられる。しかしながら、図10Bに示すように、第6の時相の血流画像40の領域35には、クラッタがほとんど含まれていない。
図11Aは、第1の比較例において生成された血流画像データが示す血流画像を示す図である。図11Bは、第1の実施形態において生成された血流画像データが示す血流画像を示す図である。
より具体的には、図11Aは、図3Cに示す固有値が比較的大きい時相(第4の時相)の血流画像104に対して座標変換が行われてモニタに表示された血流画像120を示す図である。図11Bは、図100Bに示す固有値が比較的大きい時相(第6の時相)の血流画像40に対して座標変換が行われてモニタ2に表示された血流画像50を示す図である。
図11Aに示すように、第1の比較例では、血流画像120にクラッタ120aが含まれている。しかしながら、図11Bに示すように、第1の実施形態では、血流画像50にクラッタが含まれていないか、又は、血流画像50に含まれるクラッタが大幅に低減される。このため、第1の実施形態によれば、心拍動のように周期的に移動するとともに比較的高速に移動する強反射体にサイドローブが反射されることによって発生するクラッタ成分を抑制することができる。また、第1の実施形態によれば、このようなクラッタ成分を抑制することが可能なMTIフィルタを生成することができる。
なお、上述した第1の実施形態では、超音波診断装置が、肝臓を超音波走査の対象とする場合について説明した。しかしながら、超音波診断装置は、肝臓以外の心臓の周期的な拍動が影響する他の臓器を超音波走査の対象としてもよい。例えば、超音波診断装置は、周期的に動く頸動脈の近傍に存在する甲状腺を超音波走査の対象としてもよい。
また、上述した第1の実施形態では、ブロック分割機能145は、図6に示すように画像化対象範囲を小ブロックに分割する。これにより、相関行列計算機能141は、小ブロックの相関行列の加算で中ブロックの相関行列を計算できる。また、中ブロックは小ブロック単位で重なっているので相関行列の計算に無駄がない。更に、小ブロック単位で空間データとパケットデータをメモリに格納しておくことでMTIフィルタ行列演算をCPUで実行する場合にキャッシュ効率が良い。
また、上述した第1の実施形態では、ブロック分割機能145は、中ブロックの一部が他の中ブロックの一部と重複するように中ブロックを構成する場合について説明したが、実施形態はこれに限定されるものではない。例えば、ブロック分割機能145は、中ブロックの一部が他の中ブロックの一部とは重複しないように中ブロックを構成してもよい。
また、上述した実施形態では、推定機能144は、複数の中ブロックに分割された画像化対象範囲の各中ブロックで作成されたMTIフィルタWを施したデータから移動体情報を抽出するものとして説明したが、実施形態はこれに限定されるものではない。例えば、推定機能144は、画像化対象範囲全体で作成されたMTIフィルタを施したデータから移動体情報を抽出してもよい。すなわち、相関行列計算機能141は、複数の走査線により形成される画像化対象範囲全体での超音波走査により収集された同一位置の反射波データのデータ列から、画像化対象範囲全体の相関行列を計算してもよい。また、計算機能142は、所定時間の間の相関行列の固有値に基づいて、MTIフィルタに与えられるフィルタ係数を決定してもよい。
また、上述した実施形態では、計算機能142は、中ブロックの相関行列R2からランクカット数kfを決定するものとして説明したが、実施形態はこれに限定されるものではない。例えば、計算機能142は、画像化対象範囲全体の相関行列Rallからランクカット数kfを決定してもよい。
(第1の実施形態の変形例)
上述した第1の実施形態では、ドプラ処理回路14において相関行列を固有値分解して固有値及び固有ベクトルを計算し、固有値及び固有ベクトルに基づいてMTIフィルタ行列Wを決定する場合について説明した。しかしながら、ドプラ処理回路14において相関行列を特異値分解して特異値及び右特異ベクトルを計算し、特異値及び右特異ベクトルに基づいてMTIフィルタ行列Wを決定してもよい。そこで、このような変形例を第1の実施形態の変形例として説明する。
第1の実施形態の変形例では、統計的性質として行列XHを特異値分解する場合について説明する。例えば、上記式(1)は、全空間地点の数をN、パケットサイズをLとする場合、行列Xの列ベクトルをxiとしたL×Nの行列Xを用いて以下の式(12)で表わされる。
ドプラ処理回路14の計算機能142は、以下の式(13)を用いて行列XHを特異値分解する。
なお、式(13)における行列Qの列ベクトルは、右特異ベクトルである。
式(13)を式(12)に代入すると、Pはユニタリー行列であるから、相関行列R1は以下の式(14)で表わされる。
式(4)と式(14)を比較すると、行列Vは、以下の式(15)で表わされ、行列Dは、以下の式(16)で表わされる。
式(16)において、ΛはN×Lの行列であり、エルミート行列R1の固有値は正である。このため、行列Λは、以下の式(17)で表わされる。
なお、行列Λの対角要素(λn(n=1,2,・・・L))の値は、特異値である。
そして、計算機能142は、所定時間(例えば、1.2秒)の間の特異値及び右特異ベクトルに基づいて、MTIフィルタ行列Wを決定する。
なお、かかる場合、MTIフィルタ処理機能143は、以下に示す式(18)を用いて各点のパケット列ベクトルデータXにMTIフィルタWを施す。
このように、第1の実施形態の変形例では、ドプラ処理回路14は、行列XHを特異値分解することによって、MTIフィルタ行列Wを決定する。
上記説明において用いた「プロセッサ」という文言は、例えば、CPU(Central Processing Unit)、GPU(Graphics Processing Unit)、或いは、特定用途向け集積回路(Application Specific Integrated Circuit:ASIC)、プログラマブル論理デバイス(例えば、単純プログラマブル論理デバイス(Simple Programmable Logic Device:SPLD)、複合プログラマブル論理デバイス(Complex Programmable Logic Device:CPLD)、及びフィールドプログラマブルゲートアレイ(Field Programmable Gate Array:FPGA))等の回路を意味する。プロセッサは記憶回路に保存されたプログラムを読み出し実行することで機能を実現する。なお、記憶回路にプログラムを保存する代わりに、プロセッサの回路内にプログラムを直接組み込むよう構成しても構わない。この場合、プロセッサは回路内に組み込まれたプログラムを読み出し実行することで機能を実現する。なお、本実施形態の各プロセッサは、プロセッサごとに単一の回路として構成される場合に限らず、複数の独立した回路を組み合わせて1つのプロセッサとして構成し、その機能を実現するようにしてもよい。さらに、図1における複数の構成要素を1つのプロセッサへ統合してその機能を実現するようにしてもよい。
上記の実施形態の説明において、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部又は一部を、各種の負荷や使用状況等に応じて、任意の単位で機能的又は物理的に分散・統合して構成することができる。さらに、各装置にて行なわれる各処理機能は、その全部または任意の一部が、CPUおよび当該CPUにて解析実行されるプログラムにて実現され、或いは、ワイヤードロジックによるハードウエアとして実現され得る。
また、上記の実施形態で説明した医用画像処理方法は、予め用意された医用画像処理プログラムをパーソナルコンピュータやワークステーション等のコンピュータで実行することによって実現することができる。パーソナルコンピュータやワークステーション等のコンピュータは、医用画像処理装置の一例である。医用画像処理プログラムは、インターネット等のネットワークを介して配布することができる。また、この医用画像処理プログラムは、ハードディスク、フレキシブルディスク(FD)、CD-ROM、MO、DVD等のコンピュータで読み取り可能な記録媒体に記録され、コンピュータによって記録媒体から読み出されることによって実行することもできる。
以上説明した少なくとも一つの実施形態によれば、周期的に移動するとともに比較的高速に移動する強反射体にサイドローブが反射されることによって発生するクラッタ成分を抑制することが可能なMTIフィルタを生成することができる。
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。