JPWO2011052602A1 - 超音波イメージング装置、超音波イメージング方法および超音波イメージング用プログラム - Google Patents

Abstract

エコー信号が微弱なノイズ領域の判別が可能な超音波イメージング装置を提供することにある。受信信号を処理して得た２フレーム以上の画像から基準フレームと比較フレームとを選択する。基準フレームに関心領域を設定し、比較フレームに関心領域よりも広い探索領域を設定し、探索領域内に前記関心領域の移動先の候補である複数の候補領域を設定し、関心領域内と候補領域内の画像特性値の類似度を候補領域ごとに計算して、探索領域全体にわたる類似度の分布を求める。これにより、類似度の分布に基づき、関心領域がノイズ領域かどうか判別することが可能になる。例えば、類似度の分布において、類似度の最小値と類似度の全体の値とを比較する統計量を求める構成とする。統計量と閾値とを比較することにより関心領域の信頼度を判定することが可能になる。

Description

本発明は、超音波により生体のイメージングを行う際に組織境界を明瞭に識別することが可能な超音波撮像方法および超音波撮像装置に関する技術である。

医療画像診断に用いられる超音波撮像装置においては、例えば特許文献１に記載されているように、診断動画像（Bモード画像）の小領域の変化量に基づいて組織の弾性係数分布を推定して、硬さをカラーマップに変換して表示する方法が知られている。しかしながら、例えば腫瘍の辺縁部の場合、周囲の組織に対して、音響インピーダンスも弾性率も大きく異ならないことがあり、このような場合には、診断動画像においても弾性画像においても腫瘍と周囲の組織との境界を把握できない。

そこで、特許文献２に記載の技術では、診断動画像の動きベクトルから直接スカラー場画像を作成することにより、音響インピーダンスも弾性率も周囲と大きく異ならない組織の境界を識別可能にする手法が提案されている。

特開２００４−１３５９２９号公報 特開２００８−７９７９２号公報

従来の特許文献２に記載の技術では、２つの診断画像データについてブロックマッチング処理を行うことにより動きベクトルの推定を行っているが、推定処理の際に画像データのノイズの影響により誤りベクトルが発生する。そのため境界の識別度が劣化するという問題があった。特にエコー信号が微弱になるペネトレーション限界領域ではベクトルの推定精度が大きく劣化する。

本発明の目的は、エコー信号が微弱なノイズ領域の判別が可能な超音波イメージング装置を提供することにある。

上記目的を達成するために、本発明の第１の態様によれば、以下のような超音波イメージング装置が提供される。すなわち、本発明の超音波イメージング装置は、対象に向かって超音波を送信する送信部と、対象から到来する超音波を受信する受信部と、受信部の受信信号を処理して２フレーム以上の画像を生成する処理部とを有する。処理部は、生成した２フレーム以上の画像のうちの１のフレームを基準フレームとし、予め定めた位置、または、操作者から受け付けた位置に関心領域を設定する。別の１のフレームを比較フレームとし、関心領域よりも広い探索領域を予め定めた位置、または、操作者から受け付けた位置に設定し、探索領域内に前記関心領域の移動先の候補である複数の候補領域を設定する。関心領域内と候補領域内の画像特性値の類似度を候補領域ごとに計算して、探索領域全体にわたる類似度の分布を求める。これにより、類似度の分布に基づき、関心領域がノイズ領域かどうか判別することが可能になる。

例えば、処理部は、類似度の分布において類似度の最小値と類似度の全体の値とを比較する統計量を求め、当該統計量により関心領域の信頼度を判定する構成とする。具体的には例えば、処理部は、類似度の最小値、平均値および標準偏差を用いて上述の統計量を求め、求めた統計量と閾値とを比較することにより関心領域の信頼度を判定することが可能である。

例えば、処理部は、比較フレームにおける関心領域に対応する位置と、類似度が最小の候補領域の位置とを結ぶベクトルを生成することができ、信頼度が低いと判定された関心領域についてはベクトルをゼロ、もしくは、所定のベクトルに置き換える。これにより、誤りベクトルを除去等することができ、ベクトルの精度を向上させることができる。

例えば、処理部は、類似度分布について類似度の平均、最小値及び標準偏差を算出し、平均と最小値の差分を標準偏差で除算した分離度を統計量として用いることができる。また、例えば、類似度分布について類似度の平均及び標準偏差を算出し、標準偏差を平均で除算した変動係数を統計量として用いることも可能である。

前述した統計量と比較する閾値は以下のようにして求めることができる。例えば、複数の関心領域を設定し、それぞれについて統計量を求め、求めた統計量の値の頻度を示すヒストグラム分布を求め、ヒストグラム分布における中央値、平均値、または、ヒストグラム分布が複数の山形状を呈する場合には、山の間の谷の最小値の統計量を閾値として用いる。

また、類似度分布を平滑化処理した平滑化後類似度分布を生成し、平滑化処理前の類似度分布から平滑化後類似度分布を差し引いた差分類似度分布を求めることが可能である。これにより、類似度分布からノイズによる類似度の変動を除去できる。

上記平滑化処理は、例えば、類似度分布に、予め定めた大きさのフィルタを設定し、フィルタ内の分布を平滑化する処理を、フィルタを所定量移動させながら繰り返す手法を用いる。フィルタの大きさは、以下のように定めることができる。比較フレームにおける関心領域に対応する位置と、平滑化処理前の類似度分布において類似度が最小の候補領域の位置とを結ぶベクトルを、複数の関心領域について予め生成する。生成したベクトルのうち、最大のベクトル長をフィルタの大きさとする。

上記類似度分布をラプラシアンフィルタでフィルタリング処理して輪郭強調分布を作成し、輪郭強調分布において連続する輪郭線を抽出することにより、対象の境界を求めることも可能である。

処理部を次の構成にすることにより、腫瘍の浸潤度を判定することも可能である。すなわち、処理部は、生体腫瘍境界付近に設定された関心領域について類似度分布を求め、類似度を画像特性値とする類似度分布画像を生成し、類似度分布画像上で関心領域に対応する位置を中心として複数の異なる方向に所定長の1次元領域を設定する第1処理手段と、1次元領域内の類似度の総和を、設定した方向ごとに計算する第２処理手段と、類似度総和が最小となる方向の類似度総和と、それと直交する方向の１次元領域の類似度総和との比率を計算する第３処理手段と、比率に基づいて腫瘍の浸潤度を判定する第４処理手段と、を有する構成とする。

第３処理手段で計算された比率が予め設定した一定値よりも小さい場合には、その着目画素が境界線を構成する点であると判断することにより、境界線を求めることも可能である。

また、本発明の第２の態様によれば、以下のような超音波イメージング方法が提供される。すなわち、対象に向かって超音波を送信し、対象から到来する超音波を受信して得た受信信号を処理して２フレーム以上の画像を生成し、画像から基準フレームと比較フレームとを選択し、基準フレームに関心領域を設定し、比較フレームに前記関心領域よりも広い探索領域を設定し、探索領域内に前記関心領域の移動先の候補である複数の候補領域を設定し、関心領域内と候補領域内の画像特性値の類似度を候補領域ごとに計算して、探索領域全体にわたる類似度の分布を求める超音波イメージング方法である。

また、本発明の第３の態様によれば、以下のような超音波イメージング用プログラムが提供される。すなわち、コンピュータに、２フレーム以上の超音波画像から基準フレームと比較フレームとを選択する第１のステップ、基準フレームに関心領域を設定する第２のステップ、比較フレームに関心領域よりも広い探索領域を設定し、探索領域内に関心領域の移動先の候補である複数の候補領域を設定する第３のステップ、関心領域内と候補領域内の画像特性値の類似度を候補領域ごとに計算して、探索領域全体にわたる類似度の分布を求める第４のステップ、を実行させるための超音波イメージング用プログラムである。

本発明によれば、類似度分布に基づき関心領域がノイズの多い領域かどうかを判断することが可能になる。これにより、誤りベクトルの発生を抑制し、且つペネトレーション限界領域でも高精度なベクトル推定が可能になる。推定した動きベクトル場から変換されるスカラー場画像の精度が向上し、より適切な境界検出が可能となる。

実施形態１の超音波イメージング装置のシステム構成例を示すブロック図。 実施形態１の超音波イメージング装置による画像生成の処理手順を示すフローチャート。 図２のステップ２４のブロックマッチング処理の詳細を示すフローチャート。 図２のステップ２４のブロックマッチング処理を２層構造のファントムで説明する図。 （ａ）実施形態１の超音波イメージング装置で生成したＢモード画像例を示す図、（ｂ）実施形態１の超音波イメージング装置で生成した動きベクトル分布画像例を示す図。 （ａ）は図５（ｂ）の位置(3)にＲＯＩを設定した得たＳＡＤ分布画像例を示す図、（ｂ）は図５（ｂ）の位置(5)にＲＯＩを設定した得たＳＡＤ分布画像例を示す図、（ｃ）は図（ａ）に示したＳＡＤ値のヒストグラム、（ｄ）は図（ｂ）に示したＳＡＤ値のヒストグラム。 実施形態１における分離度の定義をＳＡＤ値のヒストグラム上で示す説明図。 図２のステップ２５の処理を、分離度を用いて行う場合の処理の詳細を示すフローチャート。 図８の処理で求めた分離度の分布を画像化した例を示す図。 図２のステップ２５の処理を、変動係数を用いて行う場合の処理の詳細を示すフローチャート。 図１０の処理で求めた変動係数の分布を画像化した例を示す図。 （ａ）図２のステップ２４で生成した動きベクトル分布画像例を示す図、（ｂ）。図２のステップ２５で誤ベクトルを除去したベクトル分布画像例を示す図。 実施形態２のＳＡＤ分布の演算の際にノイズを除去する処理のフローチャート。 （ａ）実施形態２のノイズ除去前のＳＡＤ分布画像の例を示す図、（ｂ）ＳＡＤ分布画像を平滑化（ＬＰＦ）処理して得たＳＡＤ分布画像の例を示す図、（ｃ）ノイズ除去後のＳＡＤ分布画像の例を示す図。 実施形態３の浸潤度を求める処理を示すフローチャート。 （ａ）〜（ｈ）図１５の処理で用いるＳＡＤ分布の領域選択パターンを示す説明図。 （ａ）は図５（ｂ）の位置(1)にＲＯＩを設定して得たＳＡＤ分布画像を示す図、（ｂ）は図５（ｂ）の位置(2)にＲＯＩを設定して得たＳＡＤ分布画像を示す図、（ｃ）は図（ａ）に示したＳＡＤ値分布にラプラシアンフィルタを適用した画像を示す図、（ｄ）は図（ｂ）に示したＳＡＤ値分布にラプラシアンフィルタを適用した画像を示す図。

本発明の一実施形態の超音波イメージング装置について以下説明する。
（実施形態１）
図１に本実施形態の超音波イメージング装置のシステム構成を示す。本装置は、超音波境界検出機能を備えている。図１のように、本装置は、超音波探触子（プローブ）１、ユーザインタフェース２、送波ビームフォーマ３、制御系４、送受切り替えスイッチ５、受波ビームフォーマ６、包絡線検波部７、スキャンコンバータ８、処理部１０、パラメータ設定部１１、合成部１２および表示部１３を備えて構成される。

超音波素子が一次元に配列された超音波探触子１は、生体に超音波ビーム（超音波パルス）を送信し、生体から反射されたエコー信号（受波信号）を受信する。制御系４の制御下で、送波焦点に合わせた遅延時間をもつ送波信号が送波ビームフォーマ３により出力され、送受切り替えスイッチ５を介して超音波探触子１に送られる。生体内で反射又は散乱されて超音波探触子１に戻った超音波ビームは、超音波探触子１によって電気信号に変換され、送受切り替えスイッチ５を介し受波ビームフォーマ６に受波信号として送られる。

受波ビームフォーマ６は、９０度位相がずれた2つの受波信号をミキシングする複素ビームフォーマであり、制御系４の制御下で受信タイミングに応じて遅延時間を調整するダイナミックフォーカスを行って、実部と虚部のＲＦ信号を出力する。このＲＦ信号は包絡線検波部７によって検波されてからビデオ信号に変換され、スキャンコンバータ８に入力されて画像データ（Ｂモード画像データ）に変換される。以上説明した構成は、周知の超音波イメージング装置の構成と同じである。

本発明の装置では、超音波境界検出処理を処理部１０により実現する。処理部１０は、ＣＰＵ１０ａとメモリ１０ｂとを有し、あらかじめメモリ１０ｂに格納しておいたプログラムをＣＰＵ１０ａが実行することにより、以下の処理を行い被検体組織の境界を検出する。すなわち、スキャンコンバータ8から出力される2フレーム以上の画像データに基づいて、処理部１０においてまず動きベクトル場を作成する。次に作成された動きベクトル場をスカラー場に変換処理する。そして、元の画像データと、それに対応する動きベクトル場またはスカラー場とを合成部１２にて合成処理した後、表示部１３にて表示する。

パラメータ設定部１１では、処理部１０での信号処理のためのパラメータや合成部１２での表示画像の選択設定等を行う。これらのパラメータは、オペレータ（装置操作者）によりユーザインタフェース２から入力される。信号処理のためのパラメータとしては、例えば、所望のフレームｍ上の関心領域の設定や、フレームｍとは異なるフレームｍ+Δ上の探索領域の設定をオペレータから受け付けることが可能である。表示画像の選択設定としては、例えば、元画像とベクトル場画像（またはスカラー画像）とを１画像に合成してディスプレイに表示するか、あるいは２画像以上の動画像を並べて表示するかの選択設定を操作者から受け付けることが可能である。

図２に、本発明の処理部１０及び合成部１２における境界検出処理および画像処理の一例のフローチャートを示す。処理部１０は、まず、スキャンコンバータ８から計測信号を取得し、通常の信号処理を施してＢモード動画像を作成する（ステップ２１，２２）。次にＢモード動画像から所望のフレームとそれとは異なる時間のフレームとの2フレームを抽出する（ステップ２３）。例えば所望フレームと次フレームの２フレームを抽出する。そして２つのフレームから動きベクトル場を算出する（ステップ２４）。動きベクトル場の算出は、ブロックマッチング方法に基づいて実施する。算出した動きベクトル場に対してノイズ除去処理を行い（ステップ２５）、ノイズ除去した動きベクトル場をスカラー場に変換する（ステップ２６）。そしてスカラー場画像と、動きベクトル場画像あるいはＢモード画像とを合成表示して１画像分の処理を終了する（ステップ２７）。ステップ２３において所望フレームとして、時系列で順次異なるフレームを選択して上記ステップ２１〜２７の処理を繰り返し、合成画像を連続表示することにより、合成画像の動画を表示することが可能となる。

図３は、ステップ２４の詳細な処理を示すフローチャートであり、図４は、ステップ２４のブロックマッチング処理を説明する図である。図３および図４を用いて、上記ステップ２４の動きベクトル場算出のためのブロックマッチング処理を具体的に説明する。ここでは、ステップ２３においてｍ番目とｍ＋Δ番目のフレームが選択されているとする。例えばΔ＝1フレームとする。まず、処理部１０は、フレームｍにおいて、所定のピクセル数Ｎの関心領域ＲＯＩ(region of interest：基準ブロック)３１を図４のように設定する（ステップ５１）。ＲＯＩ３１に含まれるピクセルの輝度分布をP_m(i_０、j_０)と表す。i_０、j_０は、ＲＯＩ３１内の当該ピクセルの位置を示す。次に処理部１０は、フレームｍ＋Δ番目において、フレームｍのＲＯＩ３１に対応する位置およびその近傍に所定の大きさの探索領域３２を設定する（ステップ５２）。ここでは、ＲＯＩ３１の設定は、処理部１０がフレームｍの画像全体に渡って順次ＲＯＩ３１を設定し、それを中心とした所定の大きさの探索領域３２を設定する構成について説明するが、処理部１０が予め定められた位置および大きさのＲＯＩ３１と、その近傍に予め定められた大きさの探索領域３２を設定することや、パラメータ設定部１１においてオペレータから受け付けた関心領域（ＲＯＩ）および探索領域に、処理部１０がＲＯＩ３１および探索領域３２を設定することも可能である。

探索領域３２は、ＲＯＩ３１と等しい大きさの複数の移動候補領域３３に区切られている。処理部１０は、ＲＯＩ３１の輝度に対して最も類似度の高い移動候補領域３３を計算により求め、移動先領域として選択する。類似度を表す指標としては、差分絶対値和、二乗平均、相互相関値等を用いることができる。ここでは、一例として差分絶対値和を用いる場合について以下説明する。

探索領域３２内の移動候補領域３３に含まれるピクセルの輝度分布をP_m＋Δ(i、j)と表す。i、jは、移動候補領域３３内の当該ピクセルの位置を示す。処理部１０は、ＲＯＩ３１のピクセルの輝度分布P_m＋Δ(i、j)と移動候補領域３３の輝度分布P_m(i_０、j_０)との差分絶対値和ＳＡＤ(sum of absolute difference)を計算する（ステップ５３）。ここでＳＡＤは次式（１）により定義される。

処理部１０は、探索領域３２のすべての移動候補領域３３について、ＲＯＩ３１とのＳＡＤ値を求め、求められたＳＡＤ分布の中でＳＡＤ値が最小となる移動候補領域３３を移動先領域と判定し、ＲＯＩ３１の位置と最小ＳＡＤ値の移動候補領域３３の位置とを結ぶ動きベクトルを決定する（ステップ５４）。

そして処理部１０は、ＲＯＩ３１をフレームｍの画像全体に渡って移動させながら、上記処理を繰り返すことにより、フレームｍの画像全体について動きベクトルを決定する（ステップ５５）。決定したベクトルを例えば矢印で示す画像を生成することにより、動きベクトル場（動きベクトル分布画像）を得る。

図５（ａ）、（ｂ）は、上記処理で得られたＢモード画像例と、動きベクトル分布画像の例である。図５（ａ）のＢモード画像は、ゲル基材ファントム４１、４２を２層に重ね、上側のファントム４１に超音波プローブを固定して横方向に移動させながら撮像したものである。図５（ａ）のＢモード画像をフレームｍとし、次フレーム（フレームｍ＋Δ、Δ＝１フレーム）とのブロックマッチング処理（ステップ２４）により求めた動きベクトルが図５（ｂ）である。

図５（ｂ）のように、超音波プローブ１を固定した上側ファントム４１に対応する上側１／３程度の領域は相対的に静止状態であり、下側のファントム４２は横移動を示すベクトル場（横向きの矢印）となっている。しかしながら図５（ｂ）の下側ファントム４２の下部１／３程度の領域は、矢印が斜め上方を指し、向きも一定ではなく、動きベクトルが乱れている現象が見られる。この現象は、プローブ１からの距離が遠くなる程、検出感度のＳＮ比（ＳＮＲ）が低下する事に起因しており、ペネトレーション限界を示している。すなわち、プローブ１から遠い低ＳＮＲ領域で誤ベクトルが発生している。

図６を用いて、上記ステップ２５の動きベクトル場のノイズ除去処理について説明する。

図６（ａ）、（ｂ）は、図５（ｂ）の位置(3)及び(5)にＲＯＩ３１を各々設定してステップ２４で求めた次フレームの探索領域３２の各移動候補領域３３のＳＡＤ値を、それぞれの移動候補領域３３の濃度として示した図（ＳＡＤ分布図）である。図６（ａ），（ｂ）は、それぞれ全体が探索領域３２を示し、探索領域３２は、２１×２１個の移動候補領域３３に区分けされている。移動候補領域３３のブロックサイズは３０×３０ピクセル、探索領域３２は５０×５０ピクセルとし、移動候補領域３３を探索領域３２内で１ピクセルずつ移動させた。探索領域３２は、ＲＯＩ３１の位置に対応する移動候補領域３３の位置が、探索領域３２の中心に位置するように設定している。

図６（ａ）のＳＡＤ分布では、探索領域３２の中心から右側真横にずれた位置の移動候補領域３３のＳＡＤ値が最小になっている。よって、図５（ｂ）の位置(3)については、ステップ２４において右横向きのベクトルが決定されることが分かる。図５（ｂ）から確認されるように、位置(3)は２層ファントム４１，４２の境界の少し下側（ファントム４２側）付近であり、右横向きのベクトルが表示されている。

ここで図６（ａ）のＳＡＤ値の空間分布に着目すると、ＳＡＤ値が最小の移動候補領域３３の周囲に、横方向、すなわち２層ファントム４１，４２の境界に沿った方向にＳＡＤ値の小さい領域（ＳＡＤ値の谷）が形成されていることがわかる。この現象は、図６（ｂ）のように、ＲＯＩ３１をフレームｍの全領域に移動させながら全ての動きベクトル場を作成せずとも、図６（ａ）のＳＡＤ分布から境界を直接検出できる事を示唆している。

一方、図６（ｂ）のＳＡＤ値分布は、図５（ｂ）の動きベクトルが乱れた領域内の位置(5)のものであるため、ＳＡＤ値が最小になる移動候補領域３３は、ノイズ値が小さくなるプローブ１に近い上部の広い領域に一様に広がっており、動きベクトルは上向きで乱れやすくなることが分かる。また、ＳＡＤ値が最小の移動候補領域３３の周囲に形成されているはずのＳＡＤ値の小さい領域（ＳＡＤ値の谷）も、全体のノイズ変動に埋もれてしまい認識することができないことがわかる。

本発明では、図６（ａ）のように、ノイズの小さい位置(3)にＲＯＩ３１を設定したＳＡＤ分布画像においては、ＳＡＤ値が最小の移動候補領域３３が周辺領域よりも明らかに小さいＳＡＤ値をとるため、はっきりと認識でき、逆に図６（ｂ）のように、ノイズの多い位置(5)にＲＯＩ３１を設定した探索領域３２のＳＡＤ分布画像において、ＳＡＤ値が最小の移動候補領域３３が、はっきり認識できない、という現象を利用して、ノイズの多いＲＯＩ３１を判別する。ノイズの多いＲＯＩ３１については、対応する動きベクトルを除去する処理を行う。

処理部１０は上述の判別のために、例えば、図２のステップ２５において、所定の位置にＲＯＩ３１を設定した場合の図６（ａ），（ｂ）のような探索領域３２のＳＡＤ分布から、ＳＡＤ値の分布を示すヒストグラムを図６（ｃ）及び（ｄ）のように作成する。図５（ｂ）の位置(3)にＲＯＩ３１を設定した場合のヒストグラムは、図６（ｃ）のようにＳＡＤ最小値が、ヒストグラム分布中の頻度の大きいＳＡＤ値からで充分に分離している。すなわち、最小のＳＡＤ値と、頻度の大きいＳＡＤ値の範囲とが十分に離れている。一方、図５（ｂ）の位置(5)にＲＯＩ３１を設定した場合のヒストグラムは、図６（ｄ）のように、ＳＡＤ値の分布に対して頻度の差があまりなく、ヒストグラム分布がブロードに広がっている。このため、ＳＡＤ最小値が、頻度の大きいＳＡＤ値の範囲に含まれており、最小のＳＡＤ値と頻度の大きいＳＡＤ値の範囲との分離は不充分である。

このような特性から、探索領域３２のＳＡＤ分布のヒストグラムのＳＡＤ最小値と、頻度の大きいＳＡＤ値とを比較することにより、この探索領域３２に対応するＲＯＩ３１の信号の信頼度（ノイズの少なさ）、および、その探索領域３２において決定された動きベクトルの信頼度を判別することができる。これにより、信頼度の低い領域を、低ＳＮＲ領域として判別できるとともに、対応する動きベクトルの信頼度も判別することができる。

本発明では、ＳＡＤ値分布のヒストグラムにおいて、ＳＡＤ最小値と頻度の大きいＳＡＤ値との分離度を判断するために、指標を用いる。まず、分離度パラメータを指標とする場合の処理方法をまず説明する。図７に分離度の定義の概念を示す。分離度は、ヒストグラムの分布平均と最小値との距離に相当する値であり、次式（２）で定義する。

式（２）においては、分布の差異による影響を回避するため、標準偏差で規格化を施している。

図８に、分離度を指標とする場合のノイズの大きな領域を判別し、動きベクトルを除去する処理のフローを示す。この処理は、図２のステップ２５の処理を具体的に示すものであり、図３のステップ５１でフレームｍで設定した全ＲＯＩ３１について行う。処理部１０は、最初に対象とするＲＯＩ３１を決定し（ステップ８１）、決定したＲＯＩ３１に対して図３のステップ５２、５３で設定および演算された探索領域３２の全移動候補領域３３のＳＡＤ値を用いて、ＳＡＤ値の平均、最小値、標準偏差を統計処理に従い計算する（ステップ８２）。そして上式（２）で定義された分離度を求める（ステップ８３）。これを全ＲＯＩ３１について繰り返す（ステップ８４）。ステップ８４で求めた分離度が所定値よりも低いＲＯＩ３１については、ノイズが大きく、動きベクトルの信頼性が低い領域であるので、図３のステップ５４で求めた動きベクトルを０に置き換える（ステップ８５）。これにより、動きベクトル画像から、信頼性の低い領域を判別し、ベクトル（誤ベクトル）を除去することができる。

上記ステップ８５において、分離度が低いかどうか判定するための所定値としては、予め定めた閾値や、ステップ８４で全ＲＯＩ３１について求めた分離度の分布の中央値を用いることができる。もしくは、求めた分離度のヒストグラムを生成し、頻度の山が複数形成される場合、分離度が最も低い側に位置する山と、それよりも分離度が高い側に位置する山との間の谷の位置の分離度を閾値として用いることも可能である。

ステップ８３で全ＲＯＩ３１について求めた分離度を画像にしたものを図９に示す。図９においてフレームｍには３３×５１個のＲＯＩ３１が設定され、各ＲＯＩ３１の分離度が濃度により示されている。図９のように分離度は、フレームｍの下部の低ＳＮＲ領域で低くなっており、分離度が動きベクトル推定の信頼性を反映していることが分かる。

上記図８の処理では、分離度を用いたが、ＳＡＤ最小値と頻度の大きいＳＡＤ値との分離度を判断する指標としては、他の指標を用いることも可能である。例えば、変動係数を使用することが可能である。変動係数は次式で定義され、標準偏差を平均で規格化した統計量であり、分布のばらつきの大きさ（すなわち最小値の分離し難さ）を表す。

図１０に変動係数を指標とする場合のノイズの大きな領域のベクトルを除去する処理のフローを示す。この処理は、図８の処理フローと同様に対象とするＲＯＩ３１を決定し（ステップ８１）、決定したＲＯＩ３１に対して図３のステップ５２、５３で設定および演算された探索領域３２の全移動候補領域３３のＳＡＤ値を用いて、ＳＡＤ値の平均と、標準偏差を統計処理に従い計算する（ステップ１０１）。そして上式（３）で定義された変動係数を求める（ステップ１０２）。これを全ＲＯＩ３１について繰り返す（ステップ８４）。ステップ１０２で求めた変動係数が所定値よりも大きいＲＯＩ３１については、図３のステップ５４で求めた動きベクトルを０に置き換える（ステップ８５）。これにより、ノイズの大きいＲＯＩ３１を判別し、動きベクトル画像から、信頼性の低い領域を判別し、ベクトル（誤ベクトル）を除去することができる。

上記ステップ８５において、変動係数が高いかどうか判定するための所定値としては、予め定めた閾値や、ステップ８４で全ＲＯＩ３１について求めた変動係数の分布の中央値を用いることができる。もしくは、求めた変動係数のヒストグラムを生成し、頻度が二つの山を呈する場合、二つの山の間の谷の最小値を所定値として採用することも有効である。

図１１に、ステップ１０２で全ＲＯＩ３１について求めた変動係数を画像にしたものを示す。図１１において各ＲＯＩ３１の変動係数が濃度により示されている。図１１のように変動係数は、フレームｍの下部の低ＳＮＲ領域で大きくなっており、変動係数が動きベクトル推定の信頼性を反映していることが分かる。

図１２（ａ）、（ｂ）に誤ベクトル除去前と除去後の動きベクトル分布画像の例を示す。図１２（ａ）は、図４（ｂ）と同じ誤ベクトル除去前の動きベクトル場であり、図１２（ｂ）は、図１０の処理によりＳＡＤ分布から変動係数分布を求め、変動係数分布の中央値を閾値として、それより変動係数が大きいＲＯＩは信頼性が低いと判定し、動きベクトルを０（静止状態）に設定したものである。図１２（ａ）、（ｂ）を対比すると、動きベクトルが乱れている下側の領域は、誤ベクトルが明らかに除去され、静止状態に置き換えられている。すなわち、下側の領域は、超音波エコーが正確に得られないペネトレーション領域（すなわち信頼性の低い領域）と判定できている。

上述の図８および図１０の処理により、動きベクトル分布の誤りベクトルが除去されたならば、図２のステップ２６，２７により、動きベクトル分布をスカラー分布に変換し、スカラー分布画像と、動きベクトル分布画像（またはＢモード画像）とを合成表示する。

なお、図８および図１０の処理では、信頼性が低い領域の動きベクトルを除去し静止状態にしたが、本発明は、この処理方法に限られるものではない。例えば、動きベクトルを静止状態に設定する代わりに、同じ領域について前に得た動きベクトルの状態をそのまま保持する処理方法にすることも可能である。

（実施形態２）
実施形態１では、ＳＡＤ分布を求めた後に、ＲＯＩ３１が低ＳＮＲ領域かどうかをＳＡＤ分布により判別し、低ＳＮＲ領域の場合には動きベクトルを除去等するものであったが、本実施形態では、ＳＡＤ分布を演算する際に、ノイズを除去し、ノイズ除去後の検出感度の高いＳＡＤ分布を用いて信頼性の高い動きベクトルを求める。

図１３は、本実施形態２のＳＡＤ分布の演算の際にノイズを除去する処理フローである。図１３の処理は、実施形態１の図３のステップ５１〜５４のＳＡＤ演算処理において、ノイズ除去処理（ステップ１３２，１３３）を追加している。図１４（ａ）、（ｂ）、（ｃ）は、図１３の各処理段階のＳＡＤ分布の例を示している。

図１３のように、処理部１０は、まず実施形態１の図３のステップ５１〜５３と同じ処理を行い、ＳＡＤ値分布を求める。求めたＳＡＤ分布画像の例を図１４（ａ）に示す。ＲＯＩ３１は、図５（ｂ）の位置(4)であり、ファントム４２は、本来は相対的に横方向に移動しているにもかかわらず、ノイズの影響で上側の移動候補領域３３のＳＡＤ値が最小となるため、このまま動きベクトルを決定すると、ベクトルの誤検出が生じる。こうした誤検出を回避するため、処理部１０は、ステップ５１〜５３で得られたＳＡＤ分布画像に対して平滑化処理（ローパスフィルタＬＰＦ(low pass filter)処理）を施す（ステップ１３２）。平滑化処理は、例えば、所定の大きさのフィルタをＳＡＤ分布画像にかけ、フィルタ内のＳＡＤ分布の高周波成分をカットして平滑化する。この処理をフィルタを所定量ずつ移動させながら繰り返す。このようにＳＡＤ値分布画像の平滑化により、ファントム４２の移動によるＳＡＤ値の変化は急峻であるため除去されるのに対し、ノイズによるＡＤ値の変動は緩やであるため抽出することができる。平滑化処理により得られたＳＡＤ値分布画像を図１４（ｂ）に示す。

次に元のステップ５３のＳＡＤ値分布と、ステップ１３２の平滑化処理後のＳＡＤ値分布の差分を求める（ステップ１３３）。これにより、ノイズによるＳＡＤ値の変動が除去された、ファントムの移動による本来のＳＡＤ値分布を得ることができる。得られた分布を図１４（ｃ）に示す。得られたＳＡＤ値分布を用いて、図３のステップ５４を行い、最小のＳＡＤ値の移動候補領域３３を移動先と判定し、動きベクトルを決定する。動きベクトル決定後は、実施形態１の図３のステップ５６の処理により動きベクトル分布画像を生成する。また、実施形態１の図２のステップ２５により、動きベクトル分布の信頼性の低いベクトルを除去等する処理をさらに行うことも可能である。

上述したように本実施形態２の処理により、ノイズによるＳＡＤ値変動が除去されたＳＡＤ値分布を用いて動きベクトルを決定することができるため、動きベクトルの信頼性を向上させることができる。

なお、上記平滑化処理では、ＬＰＦを用いたが、これに限定されず、被検体（ファントム）の移動によるＳＡＤ値の分布の空間周波数が高い（すなわちより複雑な形状）場合には、バンドパスフィルタを適用することも有効である。

ステップ１３２のフィルタ処理において用いるフィルタの一辺の大きさは、次のように定めることができる。すなわち、平滑化処理していないＳＡＤ分布について、予め図３のステップ５４を行って動きベクトル場を作成し、その動きベクトル場の最大ベクトル長を求め、最大ベクトル長をローパスフィルタあるいはバンドパスフィルタの一辺の大きさとして設定する。

（実施形態３）
次に、実施形態３として、実施形態１の図３のステップ５３で求めたＳＡＤ値分布を用いて、直接的に組織の境界を求めるとともに、生体腫瘍の正常組織に対する浸潤度を判定する処理方法について説明する。なお、本実施形態３においては、探索領域３２の移動候補領域３３を単に領域３３と呼ぶ。また、ＲＯＩ３１を注目画素とも呼ぶ。

被検体の組織の境界に沿った領域は、組織の類似度が高いため、Ｂモード画像で近い輝度を示す。このため、ＳＡＤは、被検体の組織の境界に沿った領域３３の方が、境界の直交方向に沿った領域３３よりも小さい値を示す特徴がある。一方、生体腫瘍の浸潤が進むと、境界が不明瞭になるため、境界に沿った領域３３のＳＡＤ値が大きくなる。これを利用して、浸潤度を判定する。

図１５に本実施形態３の処理部１０の処理フローを示す。また、図１６（ａ）〜（ｈ）には、対象方向と、それに対応してＳＡＤ値分布画像上で選択される領域３３の８つのパターンを示す。

まず、処理部１０は、Ｂモード画像の所望フレームｍにおいて調査したい組織の境界位置に注目画素（ＲＯＩ）３１を設定し、フレームｍ＋Δに探索領域３２を設定し、探索領域３２のＳＡＤ値分布を求める（ステップ１５１）。フレームの選択方法およびＳＡＤ値分布の演算方法は、実施形態１の図２のステップ２１〜２３ならびに図３のステップ５１〜５３と同様に行う。

このＳＡＤ値分布において、図１６（ａ）のように探索領域３２の中心を通り、所定の対象方向（水平方向）１５１に沿って位置する領域３３を選択し（ステップ６３）、選択した領域３３のＳＡＤ値の和を求める（ステップ６４）。同様に、対象方向１５１の直交方向（垂直方向）１５２に沿って位置する領域３３を選択し、選択した領域３３のＳＡＤ値の和も求める。

このステップ６３，６４の処理を、図１６（ａ）〜（ｈ）の８つのパターンについてすべて行うまで繰り返す（ステップ６２）。図１６（ｂ）のパターンでは、所定の対象方向（水平方向に対して反時計回りに約３０°傾斜した方向）１５１に沿って位置する領域３３のＳＡＤ値の和を求め、その直交方向１５２に沿って位置する領域３３のＳＡＤ値の和を求める。

図１６（ｃ）〜（ｈ）のパターンでは、水平方向に対して反時計回りに約４５°、約６０°、９０°、約１２０°、約１３５°、約１５０°の対象方向１５１と、その直交方向１５２に沿って位置する領域３３のＳＡＤ値の和をそれぞれ求める。

求めた各対象方向１５１のＳＡＤ値の和が最小となる対象方向１５１を選択する（ステップ６５）。この選択した対象方向１５１の方向が組織の境界の方向である。これにより、動きベクトルを求めずに直接的に境界を検出することができる。

次に、選択した対象方向１５１に直交する方向１５２を選択する（ステップ６６）。選択した対象方向１５１のＳＡＤ値の和と、それに直交する方向１５２のＳＡＤ和の比率（対象方向のＳＡＤ値和／直交方向のＳＡＤ値和）を算出する（ステップ６７）。

浸潤度が低く、境界が明確な場合には、境界方向（選択した対象方向１５１）のＳＡＤ和は小さく、直交方向１５２のＳＡＤ和は大きくなるので、比率としては小さな値が得られる。一方、浸潤度が大きくなり境界が不明瞭になるにつれて、境界方向（選択した対象方向１５１）のＳＡＤ和が増大するため、比率は増大して行く。したがって、比率をパラメータとして浸潤度を評価する事ができる。具体的には例えば、予め定めておいた複数の基準値と比率とを対比して、浸潤の度合いを判定し、判定結果を表示する。

なお、比率が予め設定した一定値よりも小さい場合には、その注目画素（ＲＯＩ３１）が境界線を構成する点であると識別して境界を表示することも可能である。

なお、各方向のＳＡＤ値の和を求める際には、方向依存フィルタを用いることができる。方向依存フィルタは、処理画素のフィルタ範囲（探索領域３２）において1次元方向の濃度変化が最小の方向を判定する機能を有するフィルタである。

図１６（ａ）〜（ｈ）では、図示を容易にするために、５×５の領域３３からなる探索領域３２における対象方向１５１および直交方向１５２の領域選択パターンを示しているが、実際の処理では、探索領域３２の領域３３の数に対応して領域選択パターンを設定する。

（実施形態４）
実施形態４として、探索領域３２のＳＡＤ値分布から、動きベクトルを求めずに、直接的に境界を検出する他の方法を説明する。ここでは、２次微分に対応する強調処理を行うラプラシアンフィルタを適用する。

図１７（ａ）及び（ｂ）に、図５（ｂ）の位置(1)と(2)にＲＯＩ３１を設定し、探索領域３２のＳＡＤ分布を示す。位置(1)は、プローブ１に対して相対的に静止しているファントム４１の内部の位置である。位置(2)は、ファントム４１とそれに対して相対的に横移動しているファントム４２の境界付近に位置する。位置(1)、(2)のＳＡＤ分布を求める処理は、実施形態１の図２のステップ２１〜２３および図３のステップ５１〜５３と同様に行う。得られたＳＡＤ分布画像（図１７（ａ）及び（ｂ））に対して、空間２次微分を行うラプラシアンフィルタを適用すると、輪郭強調効果により、ＳＡＤ値の変動の大きい部分が強調された図１７（ｃ）及び（ｄ）の画像がそれぞれ得られる。

探索領域３２中に境界が存在する場合には、図１７（ｄ）に示すように、境界に沿ったＳＡＤ値の変動の大きな領域が強調されて筋状に抽出された輪郭強調分布が生成される。よって、輪郭強調分布において連続する輪郭線の領域（移動候補領域３３）を2値化処理により抽出表示することにより、ＳＡＤ値のラプラシアン画像から境界を検出することができる。一方、探索領域３２中に境界が存在しない場合には、図１７（ｃ）のように、ＳＡＤ値の変動の大きな領域は、中心の領域（ＲＯＩ３１の位置に対応する移動候補領域３３）のみであり、周囲から孤立しており、連続した輪郭線にはならない。よって、境界は、この位置にはないことがわかる。

上記処理によれば、ＳＡＤ分布にラプラシアン処理を施して直接的に境界を抽出可能となる。よって、実施形態１の動きベクトルを決定して画像化するための図３のステップ５４、５５、ならびに、動きベクトルのノイズ除去やスカラー分布への変換を行い、境界推定する図２のステップ２５〜２６を省略することができるため、処理量を大幅に削減できる。

本発明は、医用超音波診断・治療装置、ならびに、超音波を含めた電磁波一般を用いて歪みやずれを計測する装置全般に適用可能である。

１：超音波探触子（プローブ）、２：ユーザインタフェース、３：送波ビームフォーマ、４：制御系、５：送受切り替えスイッチ、６：受波ビームフォーマ、７：包絡線検波部、８：スキャンコンバータ、１０：処理部、１０ａ：ＣＰＵ，１０ｂ：メモリ、１１：パラメータ設定部、１２：合成部、１３：表示部。

探索領域３２内の移動候補領域３３に含まれるピクセルの輝度分布をP_m＋Δ(i、j)と表す。i、jは、移動候補領域３３内の当該ピクセルの位置を示す。処理部１０は、ＲＯＩ３１のピクセルの輝度分布P _m (i _０ 、j _０ )と移動候補領域３３の輝度分布P _m＋Δ (i、j)との差分絶対値和ＳＡＤ(sum of absolute difference)を計算する（ステップ５３）。ここでＳＡＤは次式（１）により定義される。

処理部１０は、探索領域３２のすべての移動候補領域３３について、ＲＯＩ３１とのＳＡＤ値を求め、求められたＳＡＤ分布の中でＳＡＤ値が最小となる移動候補領域３３を移動先領域と判定し、ＲＯＩ３１の位置と最小ＳＡＤ値の移動候補領域３３の位置とを結ぶ動きベクトルを決定する（ステップ５４）。

図１２（ａ）、（ｂ）に誤ベクトル除去前と除去後の動きベクトル分布画像の例を示す。図１２（ａ）は、図５（ｂ）と同じ誤ベクトル除去前の動きベクトル場であり、図１２（ｂ）は、図１０の処理によりＳＡＤ分布から変動係数分布を求め、変動係数分布の中央値を閾値として、それより変動係数が大きいＲＯＩは信頼性が低いと判定し、動きベクトルを０（静止状態）に設定したものである。図１２（ａ）、（ｂ）を対比すると、動きベクトルが乱れている下側の領域は、誤ベクトルが明らかに除去され、静止状態に置き換えられている。すなわち、下側の領域は、超音波エコーが正確に得られないペネトレーション領域（すなわち信頼性の低い領域）と判定できている。

Claims (14)

1. 対象に向かって超音波を送信する送信部と、前記対象から到来する超音波を受信する受信部と、前記受信部の受信信号を処理して２フレーム以上の画像を生成する処理部とを有し、
   前記処理部は、前記生成した２フレーム以上の画像のうち、１のフレームを基準フレームとし、予め定めた位置、または、操作者から受け付けた位置に関心領域を設定し、別の１のフレームを比較フレームとし、前記関心領域よりも広い探索領域を予め定めた位置、または、操作者から受け付けた位置に設定し、前記探索領域内に前記関心領域の移動先の候補である複数の候補領域を設定し、前記関心領域内と前記候補領域内の画像特性値の類似度を前記候補領域ごとに計算して、前記探索領域全体にわたる前記類似度の分布を求めることを特徴とする超音波イメージング装置。
2. 請求項１に記載の超音波イメージング装置において、前記処理部は、前記類似度の分布において前記類似度の最小値と前記類似度の全体の値とを比較する統計量を求め、当該統計量により前記関心領域の信頼度を判定することを特徴とする超音波イメージング装置。
3. 請求項２に記載の超音波イメージング装置において、前記処理部は、前記類似度の最小値、平均値および標準偏差を用いて前記統計量を求め、求めた前記統計量と閾値とを比較することにより前記関心領域の信頼度を判定することを特徴とする超音波イメージング装置。
4. 請求項３に記載の超音波イメージング装置において、前記処理部は、前記比較フレームにおける前記関心領域に対応する位置と、前記類似度が最小の前記候補領域の位置とを結ぶベクトルを生成し、前記信頼度が低いと判定された前記関心領域については前記ベクトルをゼロ、もしくは、所定のベクトルに置き換えることを特徴とする超音波イメージング装置。
5. 請求項２に記載の超音波イメージング装置において、前記処理部は、前記類似度分布について類似度の平均、最小値及び標準偏差を算出し、前記平均と前記最小値の差分を前記標準偏差で除算した分離度を前記統計量として求めることを特徴とする超音波イメージング装置。
6. 請求項２に記載の超音波イメージング装置において、前記処理部は、前記類似度分布について類似度の平均及び標準偏差を算出し、前記標準偏差を前記平均で除算した変動係数を前記統計量として求めることを特徴とする超音波イメージング装置。
7. 請求項３に記載の超音波イメージング装置において、
   前記処理部は、複数の前記関心領域を設定し、それぞれについて前記統計量を求め、求めた統計量の値の頻度を示すヒストグラム分布を求め、
   前記ヒストグラム分布における中央値、平均値、または、ヒストグラム分布が複数の山形状を呈する場合には、山の間の谷の最小値の統計量を前記閾値として用いることを特徴とする超音波イメージング装置。
8. 請求項１に記載の超音波イメージング装置において、
   前記処理部は、前記類似度分布を平滑化処理した平滑化後類似度分布を生成し、平滑化処理前の前記類似度分布から前記平滑化後類似度分布を差し引いた差分類似度分布を求めることを特徴とする超音波イメージング装置。
9. 請求項８に記載の超音波イメージング装置において、前記平滑化処理は、前記類似度分布に、予め定めた大きさのフィルタを設定し、該フィルタ内の分布を平滑化する処理を、フィルタを所定量移動させながら繰り返すものであり、
   前記処理部は、前記比較フレームにおける前記関心領域に対応する位置と、前記平滑化処理前の類似度分布において前記類似度が最小の前記候補領域の位置とを結ぶベクトルを、複数の前記関心領域について生成し、生成したベクトルのうち、最大のベクトル長を前記フィルタの大きさとすることを特徴とする超音波イメージング装置。
10. 請求項１に記載の超音波イメージング装置において、前記処理部は、前記類似度分布をラプラシアンフィルタでフィルタリング処理して輪郭強調分布を作成し、前記輪郭強調分布において連続する輪郭線を抽出することにより、前記対象の境界を求めることを特徴とする超音波イメージング装置。
11. 請求項１に記載の超音波イメージング装置において、前記処理部は、生体腫瘍境界付近に設定された前記関心領域について前記類似度分布を求め、該類似度を画像特性値とする類似度分布画像を生成し、該類似度分布画像上で前記関心領域に対応する位置を中心として複数の異なる方向に所定長の1次元領域を設定する第1処理手段と、
    前記1次元領域内の前記類似度の総和を、設定した方向ごとに計算する第２処理手段と、
    前記類似度総和が最小となる方向の類似度総和と、それと直交する方向の１次元領域の類似度総和との比率を計算する第３処理手段と、
    前記比率に基づいて腫瘍の浸潤度を判定する第４処理手段と、
    を有することを特徴とする超音波イメージング装置。
12. 請求項１に記載の超音波イメージング装置において、前記処理部は、前記第３処理手段で計算された前記比率が予め設定した一定値よりも小さい場合には、その関心領域が境界線を構成する点であると判断することを特徴とする超音波イメージング装置。
13. 対象に向かって超音波を送信し、前記対象から到来する超音波を受信して得た受信信号を処理して２フレーム以上の画像を生成し、
    前記画像から基準フレームと比較フレームとを選択し、
    前記基準フレームに関心領域を設定し、
    前記比較フレームに前記関心領域よりも広い探索領域を設定し、前記探索領域内に前記関心領域の移動先の候補である複数の候補領域を設定し、
    前記関心領域内と前記候補領域内の画像特性値の類似度を前記候補領域ごとに計算して、前記探索領域全体にわたる前記類似度の分布を求めることを特徴とする超音波イメージング方法。
14. コンピュータに、
    ２フレーム以上の超音波画像から基準フレームと比較フレームとを選択する第１のステップ、
    前記基準フレームに関心領域を設定する第２のステップ、
    前記比較フレームに前記関心領域よりも広い探索領域を設定し、前記探索領域内に前記関心領域の移動先の候補である複数の候補領域を設定する第３のステップ、
    前記関心領域内と前記候補領域内の画像特性値の類似度を前記候補領域ごとに計算して、前記探索領域全体にわたる前記類似度の分布を求める第４のステップ、
    を実行させるための超音波イメージング用プログラム。

Images