JP2010259658A - 超音波診断装置 - Google Patents

Abstract

【課題】方向性特徴量を利用したノイズ低減処理を、特に時間方向の座標に対して用い、時間方向のノイズやスペックルを低減し組織部分を強調する技術を提供する。
【解決手段】超音波画像データに対して、前記領域の空間方向、又は、時間方向及び前記領域の空間方向に、ｎ次の多重解像度分解を行い、ｋ次（１≦ｋ≦ｎ）における、低域分解画像データと高域分解画像データとを得る多重解像度分解部１０１と、特徴量を算出する特徴量算出部１０２と、特徴量に基づいてフィルタリング処理を施し処理後低域分解画像データを取得するフィルタリング処理部１０３と、元のデータと同じ画像ピクセルサイズを有する処理後超音波画像データを生成するまで多重解像度合成を繰り返す多重解像度合成部１０４とを備える。
【選択図】図１

Description

本発明は、被検体に対して超音波を送受信し超音波画像を生成する超音波診断装置に関する。さらに詳しくは、生成した超音波画像に対し補正処理を行う超音波診断装置に関する。

超音波診断装置は、超音波プローブに内蔵された振動子から発生される超音波を被検体に送信し、被検体組織からの反射波を振動子により受信して被検体の断層像などの超音波画像を生成し表示する装置である。

ところで、複数の近接する被検体組織からの受信信号は、それぞれの位相が異なるために干渉する。この干渉により、振幅のみを合成する場合とは見え方が異なる画像パターン、すなわちスペックルが生成されてしまう。スペックルは、被検体組織の境界の位置や形状の正確な観測を妨げてしまう。また、ノイズの発生によっても正確な観測が妨げられる。こういったスペックルやノイズは空間方向のみではなく時間方向にも発生する。そのためスペックルやノイズを低減するための各種の処理方法が提案されている。

例えば、時間方向のＩＩＲ（Ｉｎｆｉｎｉｔｅ Ｉｍｐｕｌｓｅ Ｒｅｓｐｏｎｓｅ）フィルタを用いて、時間的に変動するノイズやスペックルを低減する技術が提案されている。これは、時間ｔにおいて取得される画像をＩ（ｔ）とし、フィルタによる処理を行った画像をＪ_ｔとしたときに、Ｊ_ｔ＝（１−ａ）Ｉ_ｔ＋ａＪ_ｔ−１（ただし、ａ≦１）を満たすフィルタを用いて処理を行う技術である。

しかし、この技術では、簡単に表現すれば現在のデータＩ_ｔの（１−ａ）の割合に過去のデータＪ_ｔ−１のａの割合を加算するという方法でノイズを低減するために、横隔膜や心臓などの動きの激しい組織の場合、加え合わせた画像における動きの激しい組織のそれぞれの画像が残像のように残ってしまう。

さらに、走査を行った時に１枚の断層像を作成するために必要なデータの集まりの単位であるフレーム毎に画像を生成していく場合、過去のフレームにおける各画像（例えば、臓器の画像）の位置からその画像の動きベクトルを求め、その動きベクトルを用いて次のフレームにおける画像の位置を予測し、その予測した位置に対しフィルタ処理などを施す技術が提案されている。

この技術を用いれば、時間の経過による画像の位置ずれを含んだ画像処理を行うことができるが、上述のようにノイズ低減処理とは別に動きベクトルを求める処理が必要になるため、処理が重くなってしまう。

また従来、組織部分を明確にするためのフィルタであるスペックル／ノイズ低減フィルタとして、２次元や３次元の空間における方向性特徴量を利用したフィルタがある。方向性特徴量とはエッジの部分（大きさ）、エッジの向きを言う。ここで、エッジとは組織の外側などの輝度が大きく変わりかつノイズ以外の部分を指す。そして、前述の方向性特徴量を利用したフィルタとは、エッジの向き方向にはぼやかすフィルタ（たとえば、並んだ点をとって平均値を求めるフィルタなど）を用い、エッジの向きと直交する方向にはコントラストを強くするフィルタ（たとえば、エッジの外側は輝度を低減させ、エッジの内側は輝度を向上させる。具体的には閾値を超える部分の輝度を上げて、閾値を下回る部分の輝度を抑える。）を用いるフィルタである。エッジの向き方向にぼやかすフィルタを用いることで、エッジの向き方向に濃淡が滑らかな絵になる。また、エッジの向きと直交する方向にコントラストを強くするフィルタを用いることで、エッジが明確になり組織部分が強調される。

そして、上述した方向性特徴量を求めるための方法として多重解像度分解を用いる技術が提案されている。多重解像度分解とは、任意の映像信号を異なる周波数帯域（空間周波数の帯域）を有する複数の映像に分解する方法である。

多重解像度分解の一つとしてＷａｂｌｅｔ変換がある。このＷａｂｌｅｔ変換は、入力映像の２^−１サイズの低周波成分の映像と高周波成分の映像に分解する。

このＷａｂｌｅｔ変換を行った映像を基に、拡散フィルタを用いて空間周波数が低域の画像におけるエッジを検出し、画素毎にエッジの向きを求める。そして、エッジの接線方向にぼかしのフィルタを使用し画素の平滑化を行い、エッジの法線方向にコントラストのフィルタを使用し画素の先鋭化を行う技術（例えば、特許文献１参照。）が提案されている。これにより低周波成分の映像の組織部分を明確にすることが可能となる。

特開２００６−１１６３０７号公報

しかし、特許文献１に記載された多重解像度分解及び各フィルタによる画像処理は、空間方向のノイズ低減処理は行っているが、時間方向を含む画像では行われていない。そのような、空間方向のみの処理を行う方法では組織部分が強調されるだけで時間による動きは対応できていない。しかし、時間方向に対しても不要な信号の混入は避けられず、画像の時間方向に対してノイズやスペックルが発生してしまう。また、空間に対するフィルタリング処理をそのまま適用すると時相遅れが発生してしまい、生成された直後の画像（リアルタイムの画像）を表示することができなくなってしまう。

この発明は、空間方向で行っていた方向性特徴量を利用したノイズ低減処理を、時間方向を含む座標（空間＋時間）、特に時間方向の座標に対して用い、時相遅れの発生しない軽い処理を用いて超音波画像における動きの激しい組織の時間方向のノイズやスペックルを低減し組織部分を強調する技術を提供することを目的としている。

上記目的を達成するために、請求項１に記載の超音波診断装置は、被検体に対し超音波を送受信し、所定時間毎に前記被検体における特定の領域の超音波画像データを生成する画像生成手段と、前記超音波画像データに対して、前記領域の空間方向、又は、時間方向及び前記領域の空間方向に、階層的にｎ次（ｎは２次以上の自然数）の多重解像度分解を行い、ｋ次（１≦ｋ≦ｎ）における、全ての前記方向において空間周波数の低い画像である低域分解画像データと、少なくとも一つの前記方向において空間周波数の高い画像である高域分解画像データとを得る多重解像度分解手段と、算出時点の階層の１つ下の次元の出力データ又はｎ次の前記低域分解画像データからエッジの接線方向を含む当該次元における特徴量を算出する特徴量算出手段と、前記１つ下の次元の出力データ又はｎ次の前記低域分解画像データに対して、当該次元における前記特徴量に基づいて各前記方向毎に強調又は平滑化を行うフィルタリング処理を施し処理後低域分解画像データを取得するフィルタリング処理手段と、ｋ次の前記処理後低域分解画像データとｋ次の前記高域分解画像データとを多重解像度合成して（ｋ−１）次の合成画像を生成し、前記超音波画像データと同じ画像ピクセルサイズを有する処理後超音波画像データを生成するまで前記多重解像度合成を繰り返す多重解像度合成手段と、前記処理後超音波画像データを基に超音波画像を表示部に表示させる表示制御手段と、を備えることを特徴とするものである。

請求項９に記載の超音波診断装置は、被検体に対して超音波を送受信して、所定時間毎に前記被検体における同一領域の超音波画像データを生成する画像生成手段と、複数の前記超音波画像データを基に、時間方向及び空間方向に多重解像度分解して、１次については空間周波数に応じて低域分解画像及び高域分解画像に分離し、以後、（ｋ―１）次に分離された低域分解画像をｋ次の低域分解画像及び該ｋ次の高域分解画像に分離することを、ｋ＝２からｎまで行う多重解像度分解手段と、ｎ次においては、該ｎ次の低域分解画像を基にした一方のデータと及び該ｎ次の高域分解画像を基にした他方のデータとを多重解像度合成した該ｎ次の合成画像を生成し、以後、一方のデータをｋ次の合成画像とし、他方のデータをｋ―１次の高域分解画像として合成画像を生成する前記多重解像度分解の逆変換を、ｋ＝ｎから２まで行う多重分解画像合成手段と、前記ｎ次の場合は該ｎ次の多重解像度分解手段と該ｎ次の多重分解画像合成手段との間にあって該ｎ次の低域分解画像を受けて画像におけるエッジの接線方向を含む該ｎ次の特徴量を算出し、前記ｋ―１次の場合は該ｋ―１次と該ｋ次との多重分解画像合成手段の間にあって該ｋ次の合成画像を受けて画像におけるエッジの接線方向を含む該ｋ−１時の特徴量を算出することをｋ＝ｎから２まで行う特徴量算出手段と、該特徴量算出手段から前記ｎ次から１次までの特徴量を受けて、それぞれについて前記エッジの接線方向の強調又は平滑化を行い、各次における前記一方のデータとして前記多重分解画像合成手段へ送るフィルタ処理手段と、前記多重解像度分解手段から前記ｎ次から１次までの高域分解画像を受けて、対応する各次の特徴量で重み付けして前記他方のデータとして前記多重分解画像合成手段へ送る高域レベル制御部と、を備え、前記多重分解画像合成手段の１次の出力として、元の超音波画像データと同じ画像ピクセルサイズの合成画像を得ることを特徴とするものである。

請求項１及び請求項９に記載の超音波診断装置は、空間方向及び時間方向にフィルタリング処理を行う構成である。これにより、時間方向へのノイズ及びスペックルの変化を効率的に低減することが可能となる。

本発明に係る超音波診断装置のブロック図 本発明に係るスペックル／ノイズ低減処理部の詳細な構成を表す図 時間軸を含めた３次元ウェーブレット変換による多重解像度分解の概念図 時間軸を含めた３次元ウェーブレット変換における各座標軸方向に適用されるフィルタの適用パターンを示す図 （Ａ）空間方向への数値微分を説明するための図、（Ｂ）時間方向への数値微分を説明するための図 ３次の階層における非線形異方性拡散フィルタのフィルタ処理の処理手順を示すフローチャート 第１の実施形態に係る超音波診断装置の全体的な動作のフローチャート 第２の実施形態におけるＢモードデータの蓄積及びその蓄積したＢモードデータを使用した処理を説明するための図 変形例１におけるＢモードデータの蓄積及びその蓄積したＢモードデータを使用した処理を説明するための図

〔第１の実施形態〕
以下、この発明の第１の実施形態に係る超音波診断装置について説明する。本実施形態に係かる超音波診断装置は図１に示すブロック図で表わされる。図１は本実施形態に係る超音波診断装置の機能を表すブロック図である。

図１に示すように、本実施形態に係る超音波診断装置は、超音波プローブ００１、送受信部００２、信号処理部００３、画像生成部００４、スペックル／ノイズ低減処理部１００、表示制御部００５、及び、ユーザインタフェース００６を有する。そして、ユーザインタフェース００６は、モニタなどの表示部００７とキーボードやマウスなどの入力部００８を有する。

超音波プローブ００１は、複数の圧電振動子を有している。超音波プローブ００１は、後述する送受信部００２から入力されたパルス信号を圧電振動子で超音波に変換し、被検体に送信する。さらに、超音波プローブ００１は、被検体で反射した超音波（超音波エコー）を圧電素子で受信し電気信号（以下では、「エコー信号」という。）に変換して送受信部００２へ出力する。ここで、各振動子から被検体に向けた１回の超音波の送受信によって得られる超音波エコーの集合を１フレームという。つまり、１フレームにより１枚の超音波画像が生成できる。そして、単位時間あたり何フレームの走査を行うかをフレームレートという。さらに、本実施形態では、超音波プローブ００１は、被検体の同じ領域（走査位置）の走査を所定のフレームレートで繰り返し行う。

送受信部００２は、パルス信号を発生する。そして、送受信部００２は、チャンネル毎に超音波をビーム状に収束し且つ送信指向性を決定するのに必要な遅延時間を発生したパルス信号に与える。そして、送受信部００２は遅延時間を与えたパルス信号を超音波プローブ００１に出力する。

また、送受信部００２は、超音波プローブ００１からエコー信号の入力を受ける。送受信部００２は、チャンネル毎に受信したエコー信号を増幅し、さらに、受信し構成を決定するのに必要な遅延時間を与え、Ａ／Ｄ変換をおこなった後に、それらのチャンネル毎の信号を基に加算処理を行う。そして、送受信部００２は、加算処理を行った信号を信号処理部００３へ出力する。

信号処理部００３は、Ｂモード処理ユニット及びドプラ処理ユニットを有する。信号処理部００３は、送受信部００２から入力された信号を受信する。

Ｂモード処理ユニットは、送受信部００２から受けた信号に対し、対数圧縮や包絡線検波処理などを施し、信号の強度で輝度を表現するＢモードデータを生成する。

ドプラ処理ユニットは、送受信部００２から受けた信号から速度情報を周波数解析し、ドプラ効果による血流や組織、造影剤エコー成分を抽出し、平均速度、分散、パワー等の血流情報のドプラデータを生成する。

信号処理部００３は、Ｂモード処理ユニットで生成されたＢモードデータ及びドプラ処理ユニットで生成されたドプラデータの各データ（以下では、「生データ」という。）をスペックル／ノイズ低減処理部１００へ出力する。

超音波プローブ００１、送受信部００２、及び信号処理部００３が本発明における「画像生成手段」にあたる。

スペックル／ノイズ低減処理部１００は、信号処理部００３から生データの入力を受ける。そして、スペックル／ノイズ低減処理部１００は、入力された生データに対し、スペックル及びノイズの低減処理を実行する。このスペックル及びノイズの低減処理については後で詳細に説明する。そして、スペックル／ノイズ低減処理部１００は、スペックル及びノイズの低減処理を行った生データを画像生成部００４へ出力する。

画像生成部００４は、ＤＳＣ（Ｄｉｇｉｔａｌ Ｓｃａｎ Ｃｏｎｖｅｒｔｏｒ）を有する。画像生成部００４は、スペックル／ノイズ低減処理部１００からスペックル及びノイズの低減処理が行われた生データの入力を受ける。

ＤＳＣは、入力された生データを超音波スキャンの走査線に対応した座標系からモニタ等に表示するための座標系に変換し、画像データを生成する。

画像生成部００４は、メモリ等の記憶領域を有しており、ＤＳＣで生成された画像データを順次記憶していく。この記憶された画像は、例えば診断後に医師や検査技師などの操作者（以下では、単に「操作者」という。）がユーザインタフェース００６を用いて呼び出し、表示部００７に表示させて参照することが可能となっている。画像生成部００４は、生成した画像データを表示制御部００５へ出力する。

表示制御部００５は、画像生成部００４から画像データの入力を受ける。そして、表示制御部００５は、入力された画像データを基に表示部００７へＢモード画像やドプラ画像といった超音波画像を表示させる。以下で、Ｂモード画像とドプラ画像とを区別しない時には、単に「超音波画像」と呼ぶことにする。

ここで、上述した機能を有する画像生成部００４及び表示制御部００５が本発明における「表示制御手段」にあたる。

統括制御部００９は、各機能部の動作のタイミングの制御や各機能部間のデータの受け渡しの制御を行う。図１における一点鎖線は制御命令の流れを示している。

（スペックル／ノイズ低減処理）
次に、スペックル／ノイズ低減処理部１００により実行されるスペックル及びノイズの低減処理について詳細に説明する。スペックル／ノイズ低減処理部１００は、Ｂモードデータ及びドプラデータのいずれに対しても同様のスペックル及びノイズの低減処理を実行できる。そこで、以下の説明では、Ｂモードデータに対するスペックル及びノイズの低減処理を例に説明する。また、以下の説明では、説明を分かり易くするために２次元の超音波画像のデータを基にスペックル／ノイズ低減処理を説明するが、これは３次元の超音波画像のデータでも可能である。３次元の超音波画像のデータを使用する場合には、時間の次元を入れると４次元となるが、計算上は以下で説明する行列の計算が１行１列増えるだけである。また、以下の説明では、１枚のＢモード画像を構成するデータのことを１フレームのＢモードデータという場合がある。さらに、Ｂモード画像の枚数に対応したＢモードデータの数のことをフレーム数という場合がある。１フレームのＢモードデータが本発明における「１枚の超音波画像データ」にあたる。

図２はスペックル／ノイズ低減処理部１００の詳細な構成を示す図である。図２に示すように、スペックル／ノイズ低減処理部１００は、多重解像度分解から多重解像度合成を行うために、一点鎖線で示されるように複数階層（１次〜３次）からなる多重構造を有している。ここで、本実施形態においては、説明をわかり易くするために多重解像度分解から多重解像度合成の実施の最高次数を３次としている。しかしながらこの最高次数は他の数でもよい。すなわち、多重解像度分解から多重解像度合成の実施は、１次からｎ次（ただし、ｎは２以上の自然数）の範囲で行われればよい。また、本実施形態では多重解像度分解としてウェーブレット変換を用いているが、このウェーブレット変換は多重解像度分解の一例であり、多重解像度分解を実施する方法であれば他の方法と用いてもよく、例えばプラシアン・ピラミッド法やガボール変換などの他の多重解像度分解方法を用いてもよい。なお、本実施形態におけるウェーブレット変換／逆変換は、いわゆる離散ウェーブレット変換／逆変換を指す。以下の説明では図２において次数が多いほど（すなわち、図２において下の位置にあるほど）高い階層という。

スペックル／ノイズ低減処理部１００は、図２に示すように、階層毎に多重解像度分解部１０１（１０１ａ、１０１ｂ、１０１ｃ）、特徴量算出部１０２（１０２ａ、１０２ｂ、１０２ｃ）、フィルタリング処理部１０３（１０３ａ、１０３ｂ、１０３ｃ）、多重解像度合成部１０４（１０４ａ、１０４ｂ、１０４ｃ）、及び高域レベル制御部１０５（１０５ａ、１０５ｂ、１０５ｃ）を備えている。さらに、スペックル／ノイズ低減処理部１００はフレーム蓄積部１０６を備える。

フレーム蓄積部１０６は、ハードディスクやメモリなどの記憶媒体を有している。フレーム蓄積部１０６は、信号処理部００３からＢモードデータの入力を受け、最新のフレームまでの所定のフレーム数（枚数）のＢモードデータを記憶する。Ｂモードデータは時間の系列に合わせて順次生成されているので、フレーム蓄積部１０６は、次々とＢモード画像が蓄積されていく。この各Ｂモードデータはそれぞれ同じ被検体の領域を走査した超音波画像のデータであり、且つ走査された時間が異なるものである。ここでは、フレーム蓄積部１０６は、新しいフレームのＢモードデータが入力されると、一番古いフレームのＢモードデータを削除するとする。ただし、フレーム蓄積部１０６での記憶方法は必要数のフレームのＢモードデータを記憶できる方法であれば他の方法でもよく、例えば、古いフレームのＢモードデータの削除を行わずに全てのフレームのＢモードデータを記憶しておく構成でもよい。このフレーム蓄積部１０６が本発明における「処理前データ蓄積手段」にあたる。

フレーム蓄積部１０６は、図２に示すように、記憶している所定のフレーム数のＢモードデータを多重解像度分解部１０１へ出力する。以下の説明では、所定のフレーム数のＢモードデータの集合を集合Ｂモードデータという。

多重解像度分解部１０１は、フレーム蓄積部１０６から所定のフレーム数のＢモードデータ（集合Ｂモードデータ）の入力を受ける。本実施形態に係る多重解像度分解部１０１は、空間方向及び時間方向について多重解像度分解を実施する。ここで、空間方向とは被検体を走査した領域における方向のことを指し、例えば、２次元断層像であればＸ及びＹの２次元の方向、３Ｄ画像であればＸ、Ｙ、及びＺの３次元の方向となる。ここでの説明では、超音波画像が２次元断層像であるので、多重解像度分解部１０１は、時間軸も含めて３次元の座標系（時空間）に対して多重解像度分解を行う。多重解像度分解部１０１は、入力された集合Ｂモードデータに対して３次元的なウェーブレット変換を行い、集合Ｂモードデータを１種類の低域分解画像データと７種類の高域分解画像データに分解する。この多重解像度分解部１０１が本発明における「多重解像度分解手段」にあたる。

図３は時間軸を含めた３次元ウェーブレット変換による多重解像度分解の概念図である。図４は時間軸を含めた３次元ウェーブレット変換における各座標軸方向に適用されるフィルタの適用パターンを示す図である。多重解像度分解部１０１は、図３と図４のように、分開前の空間（ｘ、ｙ）及び時間（ｔ）の次元を有する集合Ｂモードデータ３００に対し、ｘｙｔ直交座標軸の各座標方向（各次元）毎に１次元の低域（Ｌ）フィルタと高域（Ｈ）フィルタとをそれぞれ適用して、集合Ｂモードデータ３００をウェーブレット変換する。このウェーブレット変換により集合Ｂモードデータ３００は、１種類の低域分解画像データ３０１と、７種類の高域分解画像データ３０２〜３０８とに分解される。低域分解画像データ３０１は、分解前の集合Ｂモードデータ３００が有する周波数成分のうち低周波成分を含む。また、各高域分解画像データは、分解前の集合Ｂモードデータ３００が有する周波数成分のうちの少なくとも１つの座標方向に関する高周波成分を含む。図４は、分解後の低域分解画像データ３０１及び高域分解画像データ３０２〜３０８がそれぞれ、座標軸別に低域フィルタと高域フィルタのいずれのフィルタが適用されているかを示している。分解後の低域分解画像データ３０１及び高域分解画像データ３０２〜３０８はそれぞれ、各座標軸あたりの標本数が分解前の集合Ｂモードデータ３００の各座標軸あたりの標本数の半分に縮小される。

多重解像度分解部１０１は多重解像度分解を最高階層に至るまで１回または複数回繰り返す。これにより、図２に示すように何層階か（本実施形態では３層）に分かれる。各階層は字数で表わされる。複数回繰り返す場合には、一つ下の階層（図２では一つ上の位置に表わされている階層）での低域分解画像データが次の多重解像度分解に対する入力となる。具体的には、多重解像度分解部１０１が最高階層（本実施形態では３次の階層）に属さない場合、すなわち、ここでは多重解像度分解部１０１ａ及び多重解像度分解部１０１ｂの場合、多重解像度分解部１０１ａは生成した低域分解画像データを多重解像度分解部１０１ｂに、多重解像度分解部１０１ｂは生成した低域分解画像データを多重解像度分解部１０１ｃに供給する。また最高階層（３次の階層）に属する多重解像度分解部１０１ｃは生成した低域分解画像データを、同階層（３次の階層）の特徴量算出部１０２ｃに供給する。

また、各階層の多重解像度分解部１０１は、生成した高域分解画像データを同階層に属する高域レベル制御部１０５に供給する。

特徴量算出部１０２は、供給された低域分解画像データに含まれる方向性特徴量であるエッジ情報、本実施形態ではエッジの大きさの情報を算出（具体的には、エッジの接線方向を用いて算出）する。特徴量算出部１０２は、算出したエッジの大きさの情報を、同階層の高域レベル制御部１０５及びフィルタリング処理部１０３に供給する。さらに、特徴量算出部１０２は、フィルタリング処理部１０３に低域分解画像データを出力する。この特徴量算出部１０２が本発明における「特徴量算出手段」にあたる。

フィルタリング処理部１０３は、特徴量算出部１０２から低域分解画像データ及びエッジの大きさの情報の入力を受ける。フィルタリング処理部１０３は、供給された低域分解画像データに３次元的な非線形異方性拡散フィルタを適用する。非線形異方性拡散フィルタは、非線形的に異方性拡散する成分を強調するフィルタである。低域分解画像データ中のエッジ成分は、非線形的に異方性拡散している。したがって、低域分解画像データに非線形異方性拡散フィルタが適用されると、エッジ成分の輝度が増大され、非エッジ成分が低減される。そのため、非線形異方性拡散フィルタの適用により、フィルタリング処理部１０３は、低域分解画像データに含まれるエッジ成分を強調し、非エッジ成分を滑らかにした低域分解画像データを生成する。このフィルタリング処理部１０３が本発明における「フィルタリング処理手段」にあたる。

ここで、特徴量算出部１０２による方向性特徴量であるエッジ情報（本実施形態ではエッジの大きさ）の検出及びフィルタリング処理部１０３による非線形異方性拡散フィルタを用いたフィルタリング処理について具体的に説明する。方向性特徴量の検出及び非線形異方性拡散フィルタを用いたフィルタリング処理は、本実施形態では低域分解画像データの構造テンソル（ｓｔｒｕｃｒｕｒｅ ｔｅｎｓｏｒ）の固有値・固有ベクトルを求めることによって実現する。この構造テンソルの固有値・固有ベクトルは、それぞれ低域分解画像データのエッジの向き、すなわち大きさを表す。特に、第１固有値とそれに対応する固有ベクトルは空間・時間的に面的なまとまりをもつ低域分解画像データの特徴を表すものである。以下では、低域分解画像データＩに対して処理を行うものとし、以下の説明では、低域分解画像データＩを単に「画像Ｉ」という。

画像Ｉの中の画素の時間方向も含めた構造テンソルは次のように定義される。

ここでＩ_ｘは画像Ｉのｘ方向の空間微分、Ｉ_ｙは画像Ｉのｙ方向の空間微分、Ｉ_ｔは画像Ｉのｔ方向（時間方向）の時間微分、Ｇ_ρはガウス関数、演算子「＊」は畳み込みを表す。λ_Ｓ１、λ_Ｓ２、及びλ_Ｓ３は、構造テンソルＳの固有値であり、その絶対値の大小関係はλ_Ｓ１≧λ_Ｓ２≧λ_Ｓ３であるとする。この場合λ_Ｓ１が第１固有値である。空間微分Ｉ_ｘ、Ｉ_ｙ、及びＩ_ｔの算出方法は、上記の算出方法に限定されない。例えば、Ｉ_ｘ、Ｉ_ｙ、及びＩ_ｔを計算する代わりに、ソーベルフィルタ（ｓｏｂｅｌ ｆｉｌｔｅｒ）や、多重解像度分解の高域分解画像データを適用してもよい。

構造テンソルＳの各要素ｓが求まっている場合における固有値λ_Ｓ１、λ_Ｓ２、及びλ_Ｓ３と固有ベクトルＲ及びＲ^Ｔとは、線形代数学でよく知られている方法によって算出することができる。すなわち、（式１）中の３×３行列の固有値λ_Ｓ１、λ_Ｓ２、及びλ_Ｓ３は、例えばカルダノ（Ｃａｒｄａｎｏ）の方法等によって、実質的に３次方程式を解くことにより算出できる。また、構造テンソルＳは、実対象行列であるので、固有値λ_Ｓ１、λ_Ｓ２、及びλ_Ｓ３は実数であり、固有ベクトルＲ及びＲ^Ｔは実ベクトルである。固有ベクトルＲ及びＲ^Ｔは互いに直交する。この構造テンソルＳの固有値λ_Ｓ１、λ_Ｓ２、及びλ_Ｓ３がエッジの大きさを表し、固有ベクトルＲ及びＲ^Ｔがエッジの方向（エッジの接線方向）を表している。

方向性特徴に基づくフィルタリングは、本実施形態では次式に示す非線形拡散方程式（偏微分方程式）で表わされる。

ここで、Ｉは画像Ｉの画素レベル（輝度値）、∇Ｉはその勾配ベクトル（ｇｒａｄｉｅｎｔ ｖｅｃｔｏｒ）、ｗは物理的な拡散方程式における時間の次元に属する。言い換えれば、ｗは処理にかかわる時刻であり、実際の処理上では、この拡散方程式の処理回数を示す。Ｄは拡散テンソル（Ｄｉｆｆｕｓｉｏｎ Ｔｅｎｓｏｒ）であって、次のように表わされる。

ここでＲは次の（式４）のように表わされる。

ここで、ω_１、ω_２、及びω_３は拡散テンソルＤの固有ベクトル、λ_Ｄ１、λ_Ｄ２、及びλ_Ｄ３は拡散テンソルＤの固有値である。そして、Ｒ＝（ω_１，ω_２，ω_３）のω_１がエッジの法線方向を表す。

拡散テンソルＤの固有値λ_Ｄ１は、固有ベクトルω_１が示す方向の拡散の強さを表す。同様に、固有値λ_Ｄ２は、固有ベクトルω_２が示す方向の拡散の強さを表し、固有値λ_Ｄ３は、固有ベクトルω_３が示す方向の拡散の強さを表す。そして、固有値λ_Ｄ１、λ_Ｄ２、及びλ_Ｄ３の値を画素毎に制御することによって、非線形異方性拡散フィルタの拡散の強さが制御される。なお、固有ベクトルＲは、非線形異方性拡散フィルタのフィルタ方向と同義である。すなわち、固有ベクトルＲを適宜設定することにより、所望の非線形異方性拡散フィルタのフィルタ方向が設定される。

ここで、拡散テンソルＤの固有ベクトルω_１、ω_２、及びω_３は、上述した画像Ｉの中の画素の構造テンソルＳのそのままの固有ベクトル又は固有ベクトルの順番を入れ替えたものに等しい。また、拡散テンソルＤの固有値λ_Ｄ１、λ_Ｄ２、及びλ_Ｄ３は、上述の構造テンソルＳの固有値λ_Ｓ１、λ_Ｓ２、及びλ_Ｓ３から計算されるエッジの大きさに依存する。

エッジ情報（ｅｄｇｅ）は、第１固有値である固有値λ_Ｓ１と第３固有値である固有値λ_Ｓ３とを用いて、次の（式５）のように表わされる。すなわち、エッジ情報（ｅｄｇｅ）は、特徴量のうちのエッジの大きさを用いて表わされる値である。

ここで、エッジは「大きさ」と「向き」有する。そして、「大きさ」とは画像の強度の変化がなだらかなときに、「小さい」エッジとなり、画像の強度の変化が急峻なときに、「大きな」エッジとなるパラメータである。

エッジ情報は、エッジの大きさを０から１までの範囲に正規化することにより算出されるパラメータである。エッジ情報の値が１に近いほどエッジ成分に近く、０に近いほど非エッジ成分に近い。パラメータｋはエッジ成分の抽出度合いを示すパラメータである。すなわちパラメータｋを小さくすると、エッジ成分が抽出されやすくなる。

さらに、このエッジ情報を用いて拡散テンソルＤの固有値λ_Ｄ１、λ_Ｄ２、及びλ_Ｄ３が算出される。固有値λ_Ｄ１、λ_Ｄ２、及びλ_Ｄ３はそれぞれ以下の（式６）、（式７）、（式８）で表わされる。

（式６）、（式７）、及び（式８）の中のβ（１−ｅｄｇｅ）の部分は非エッジ成分を表し、β・ｅｄｇｅの部分はエッジ成分を表している。この固有値λ_Ｄ１、λ_Ｄ２、及びλ_Ｄ３は上述したように拡散の強さを表すものである。そして、低減対象であるスペックルやノイズが属する非エッジ成分に対しては、その方向性を持たせない必要がある。すなわち、この非エッジ成分を方向に依らず拡散させなければならない。そこで、β_１＝β_３＝β_５＞０に設定する。一方、強調対象であるエッジ成分に対しては、その方向性をより強調させなければならない。すなわち、このエッジ成分の垂直方向（固有値λ_Ｄ１の固有ベクトルω_１が指す方向）を先鋭化させ、それ以外の方向（固有値λ_Ｄ２の固有ベクトルω_２が指す方向及び固有値λ_Ｄ３の固有ベクトルω_３が指す方向）を拡散させる必要がある。そのため、β_２は０に近似した値が設定され、β_４及びβ_６はβ_２よりも大きい所定値に設定される。

以上により、拡散テンソルＤの固有値λ_Ｄ１、λ_Ｄ２、及びλ_Ｄ３、並びに、固有ベクトルω_１、ω_２、及びω_３が算出することができる。したがって、（式３）における拡散テンソルＤの各要素ｄを算出でき、非線形異方性拡散フィルタ（拡散方程式）である（式２）を得ることができる。

さらに、この求めた非線形異方性拡散フィルタ（式２）自体の計算は、偏微分方程式の数値解析的解法によって行う。すなわち、時刻ｗのある点における画素とその周囲の各点における各画素レベル（画素の輝度値）及び拡散テンソルＤの各要素値から、時刻ｗ＋Δｗにおけるその点の新たな画素レベルを求め、次にｗ＋Δｗを新たなｗとして、同様の計算を所定回数繰り返す。

ここで、式２の数値的計算方法を説明する。そこで、まず時間方向への数値微分の方法について説明する。図５（Ａ）は空間方向への数値微分を説明するための図、図５（Ｂ）は時間方向への数値微分を説明するための図である。

図５（Ａ）（Ｂ）のデータ５０１は説明を容易にするために、１次元の空間の座標と時間方向の座標の２つの座標を有している画像としている。図５（Ａ）において、空間方向の数値微分やフィルタリングは現在の位置のデータ５０２を処理中心とし、近傍の正負両側の所定数のピクセルのデータ５０３を用いて計算する。１回の数値微分を例にとると、正負両側の各近傍１ピクセルを用いれば、１階の数値微分の近似は次の（式９）で表わされる。これは数値微分の３点公式と呼ばれる。

これに対し、図５（Ｂ）において時間方向の数値微分やフィルタリングは、現在のデータ５０４を処理中心とするが、未来の時間を処理に用いることはできない。例えば、処理中心を過去に設定して、計算式を空間方向と同じ係数を使用するとした場合、処理結果は現在のフレームより時間的に遅延したものとなり、時相遅れをなくすことができない。そのため、時間方向の数値微分やフィルタリングは現在のデータ５０４を処理中心とし、現在のデータ５０４及び過去の所定数のピクセルのデータ５０５を用いて計算しなくてはならない。

ここで、関数ｆ（ｘ）が区間［ｘ_ｉ，ｘ_ｉ＋ｈ］で連続のとき、

と表わされる。
そして、Ｒ_ｎ＋１は数値計算で打切り誤差となる項であり、｜ｈ｜が十分に小さければ、Ｒ_ｎ＋１＝０となる。

そして（式１０）より以下の式が求められる。

そして、（式１１）及び（式１２）より以下の式が求められる。

この（式１３）の｛｝内を誤差とみなして省略すると、

と表わされる。

（式１４）で示されるように、時間方向の１階の数値微分の近似は現在のデータと過去２フレーム分のピクセルのデータを用いて求めることができる。この（式１４）を用いることで（式２）の数値的計算方法による計算を行うことができる。

具体的には、（式２）は以下の（式１５）ように表わすことができる。

この（式１５）の各項の偏微分を計算するために、上述した（式９）及び（式１４）の近似値が用いられる。

ここで、以上の非線形異方性拡散フィルタの処理手順を図６を参照して説明する。図６は、３次の階層における非線形異方性拡散フィルタのフィルタ処理の処理手順を示すフローチャートである。ここでは、３次の階層のフィルタリング処理を例に説明するが、１次又は２次におけるフィルタリング処理においても同等の処理が行われる。

ステップＳ００１：多重解像度分解部１０１ｃは、多重解像度分解部１０１ｂからフィルタリング処理を行うデータの入力を受ける。

ステップＳ００２：多重解像度分解部１０１ｃは、入力されたデータに対し各座標軸方向（Ｉ_ｘ、Ｉ_ｙ、Ｉ_ｔ）に微分し低域分解画像データ及び高域分解画像データを生成する。

ステップＳ００３：特徴量算出部１０２ｃは、多重解像度分解部１０１ｃから入力された低域分解画像データの構造テンソルＳの各要素を計算する。この計算には、構造テンソルＳ計算前のガウス関数Ｇ_ρの計算も含まれる。

ステップＳ００４：特徴量算出部１０ｃは、構造テンソルＳからその固有値λ_Ｓ１、λ_Ｓ２、及びλ_Ｓ３、並びに、固有ベクトルω_１、ω_２、及びω_３を算出する。

ステップＳ００５：特徴量算出部１０２ｃは、構造テンソルＳの固有値λ_Ｓ１、λ_Ｓ２、及びλ_Ｓ３からエッジの大きさが、固有ベクトルω_１、ω_２、及びω_３からエッジの向きが求められる。

ステップＳ００６：特徴量算出部１０２ｃは、求めたエッジの大きさを同じ階層の高域レベル制御部１０５ｃへ出力する。

ステップＳ００７：フィルタリング処理部１０３ｃは、特徴量算出部１０２ｃで求められたエッジの大きさを基に拡散の強さを計算し、さらに、求めた拡散の強さと構造テンソルＳの固有ベクトルω_１、ω_２、及びω_３から、非線形異方性拡散フィルタの拡散方程式（式２）の数値解析に用いられる各係数を算出する。

ステップＳ００８：フィルタリング処理部１０３ｃは、拡散方程式の数値解析的計算を行う。

ステップＳ００９：フィルタリング処理部１０３ｃは、拡散方程式の数値解析的計算を所定回数行ったか否かを判断する。所定回数行っていない場合にはステップＳ００８に戻る。所定回数行っている場合にはステップＳ０１０に進む。

ステップＳ０１０：フィルタリング処理部１０３ｃは、非線形異方性拡散フィルタによりフィルタリング処理を行った低域分解画像データを同じ階層の多重解像度合成部１０４ｃに出力する。

高域レベル制御部１０５は、同じ階層の多重解像度分解部１０１から８個の高域分解画像データの入力を受ける。さらに、高域レベル制御部１０５は、特徴量算出部１０２からエッジの大きさの入力を受ける。高域レベル制御部１０５は、高域分解画像データに対しエッジの情報（エッジの大きさ）を用いて処理を行う。本実施形態では、拡散テンソルＤの固有値より計算されたエッジの大きさであり、高域レベル制御部１０５は、エッジの大きさと各高域分解画像データとの画素毎の積をとり、さらにその結果に各高域分解画像データの制御係数をかけている。この様にエッジの大きさとの積をとることにより、高域分解画像データにおけるエッジ部分を強調し、エッジ以外の部分を抑制することができる。

他のエッジの大きさ（エッジ情報）を用いた高域分解画像データの処理の例としては、エッジ以外の領域に各高域画像の制御係数をかける方法がある。この様な処理を行うことにより、高域分解画像データにおいてもエッジ部分を強調することができる。高域レベル制御部１０５は、上述したような処理を行った高域分解画像データを同じ階層の多重解像度合成部１０４へ出力する。この高域レベル制御部１０５が本発明における「高域レベル制御手段」にあたる。

多重解像度合成部１０４は、フィルタリング処理部１０３から１個の低域分解画像データ及び高域レベル制御部１０５から７個の高域分解画像データの入力を受ける。多重解像度合成部１０４は、入力された１個の低域分解画像データ及び７個の高域分解画像データを合成し１つの合成画像を生成する。本実施形態では、多重解像度合成部１０４は、１個の低域分解画像データ及び７個の高域分解画像データを基に逆ウェーブレット変換を行うことにより合成画像を生成する。この生成された合成画像は、各座標軸あたりの標本数がフィルタリング処理部１０３から入力された画像の２倍となる。２次以上の階層における多重解像度合成部１０４（本実施形態では、多重解像度合成部１０４ｂ又１０４ｃ）は、形成した画像を１つ下の階層の特徴量算出部１０２に出力する。また、１次の階層における多重解像度合成部１０４（本実施形態では、多重解像度合成部１０４ａ）は、形成した画像を画像生成部００４へ出力する。この多重解像度合成部１０４が本発明における「多重解像度合成手段」にあたる。

ここで、図７を参照して、本実施形態に係る超音波診断装置の全体的な動作の流れを説明する。図７は本実施形態に係る超音波診断装置の全体的な動作のフローチャートである。ここでも、Ｂモード画像を例に説明する。

ステップＳ１０１：送受信部００２は、超音波プローブ００１を介して被検体に対し超音波を送信し、その被検体で反射した超音波エコーに基づくエコー信号を受信する。

ステップＳ１０２：信号処理部００３は、送受信部００２から入力されたエコー信号に対し、対数圧縮や包絡線検波などの信号処理を行う。

ステップＳ１０３：スペックル／ノイズ低減処理部１００は、信号処理部００３から入力されたＢモードデータをフレーム毎にフレーム蓄積部１０６に記憶する。

ステップＳ１０４：スペックル／ノイズ低減処理部１００は、フレーム蓄積部１０６に記憶された最新のＢモードデータから所定フレーム数前までのＢモードデータを基に、多重解像度分解、特徴量算出、フィルタリング処理、高域レベル制御、及び多重解像度合成を行ない、最新のＢモードデータに対するスペックル／ノイズ低減処理を行う。

ステップＳ１０５：画像生成部００４は、スペックル／ノイズ低減処理部１００から入力されたスペックル／ノイズ低減処理が行われたＢモードデータに対し、座標変換などを行い、Ｂモードの超音波画像を生成する。

ステップＳ１０６：表示制御部００５は、画像生成部００４から入力された超音波画像を表示部００７に表示させる。

ステップＳ１０７：統括制御部００９は、検査が終了したか否かを判断する。検査が終了していない場合にはステップＳ１０１に戻る。検査が終了した場合には超音波診断装置の動作を停止する。

以上の説明では、超音波画像として２次元の画像を例に説明したが、これは３Ｄ超音波画像、すなわち３次元ボリュームデータを使用した場合にも同様に時間方向を含めたスペックル／ノイズ低減処理を行うことが可能である。３Ｄ超音波画像に対するスペックル／ノイズ低減処理を行う場合には、以上で説明した処理に１次元増やせばよく、エッジの大きさの算出やフィルタリング処理における計算において、単に行列が１行１列増やして計算を行うことで実施できる。

さらに、以上の説明では信号処理部００３の後にスペックル／ノイズ低減処理部１００を配置し生データに対しスペックル／ノイズ低減処理を行ったが、スペックル／ノイズ低減処理部１００は、画像生成部００４の後に配置されてもよく、その場合座標変換が行われた超音波画像に対しスペックル／ノイズ低減処理が行われる。その場合、超音波プローブ００１、送受信部００２、信号処理部００３、及び画像生成部００４が本発明における「画像生成手段」にあたり、表示制御部００５が本発明における「表示制御手段」にあたる。

以上で説明したように、本実施形態に係る超音波診断装置は、Ｂモードやドプラモードといった各モードにおける生データの空間及び時間方向に多重解像度分解を行い、方向性特徴量を算出し、その方向性特徴量を利用したフィルタを用いることによって、空間方向及び時間方向のノイズ及びスペックルの変化を低減する構成である。これにより、空間方向だけでなく時間方向にも滑らかで、かつ組織部分が強調された画像を生成することが可能となる。

また、本実施形態にかかる超音波診断装置は、時間方向への数値微分については、現在と過去の複数フレームのデータを用いて、空間方向と異なる数式と係数によって計算を行う構成である。これにより、時相遅れのないフィルタリング処理が可能となる。

また、本実施形態では、高域レベル制御部１０５を用いて高域分解画像データに対しても、処理を行い高域分解画像データのエッジ部分を強調することで、よりエッジが強調された画像を生成している。しかし、高域レベル制御部１０５を用いずに、高域分解画像データに対して処理を行わず、そのままの高域分解画像データを用いて多重解像度合成部１０４において合成画像を生成することも可能である。この場合、高域レベル制御部１０５を用いた場合と比較して組織部分の強調の程度は弱まるが、時間方向への滑らかな画像の生成は可能である。

さらに、多重解像度分解部１０１、特徴量算出部１０２、フィルタリング処理部１０３、多重解像度合成部１０４、及び高域レベル制御部１０５は、各階層（図２における１次〜ｎ次までの各次毎）それぞれに配置して構成してもよいし、１組の多重解像度分解部１０１、特徴量算出部１０２、フィルタリング処理部１０３、多重解像度合成部１０４、及び高域レベル制御部１０５の構成で、各階層（各次）の処理を順次処理する構成であってもよい。

〔第２の実施形態〕
以下、この発明の第２の実施形態に係る超音波診断装置について説明する。本実施形態に係る超音波診断装置はスペックル／ノイズ低減処理を行ったＢモードデータを利用して次のフレームに対するスペックル／ノイズ低減処理を行う構成であることが第１の実施形態と異なるものである。本実施形態に係る超音波診断装置の構成も図１のブロック図で表わされる。以下の説明では、第１の実施形態と同一の符号を付された機能部は特に説明のない限り同じ機能を有するものとする。

図８は、本実施形態におけるＢモードデータの蓄積及びその蓄積したＢモードデータを使用した処理を説明するための図である。図８は説明をわかり易くするために、１次元の空間軸と時間軸との２次元として説明する。

フレーム蓄積部１０６は、時間の経過毎に順番に信号処理部００３から入力されたＢモードデータを記憶する。ここで、最新のＢモードデータをＢモードデータ８０１とする。そして、スペックル／ノイズ低減処理部１００は、フレーム蓄積部１０６に記憶されているＢモードデータ８０１を含む所定フレーム数のＢモードデータの集合である集合Ｂモードデータ３００を用いてＢモードデータ８０１に対するスペックル／ノイズ低減処理を行う。そして、フレーム蓄積部１０６は、処理前のＢモードデータ８０１を削除し、処理が行われたＢモードデータ８０１（以下では、処理後Ｂモードデータ８０１）を最新のフレームのＢモードデータとして記憶する。その後、次のフレームのＢモードデータ（これを、Ｂモードデータ８０２という。）が信号処理部００３から入力されると、フレーム蓄積部１０６に最新のフレームとしてＢモードデータ８０２が記憶される。そして、スペックル／ノイズ低減処理部１００は、Ｂモードデータ８０２に対しスペックル／ノイズ低減処理を行い、フレーム蓄積部１０６は処理後Ｂモードデータ８０２を記憶する。以上のような処理を繰り返していくことにより、フレーム蓄積部１０６は、過去のＢモードデータの処理結果が時間方向に蓄積され、その蓄積された処理後のＢモードデータを用いて最新のＢモードデータのスペックル／ノイズ低減処理を行うことができる。ここで、本実施形態におけるフレーム蓄積部１０６が本発明における「処理後フレーム蓄積手段」にあたる。

このように、本実施形態に係る超音波診断装置は、スペックル／ノイズ低減処理を施した過去のＢモードデータを用いて最新のＢモードデータのスペックル／ノイズ低減処理を行う構成である。

これにより、時間方向のスペックル／ノイズ低減処理において既に処理を施した滑らかなＢモードデータを使用することができるため、より時間方向のスペックル／ノイズの変化を効率的に低減し、時相遅れのない処理を行うことが可能となる。

（変形例１）
次に、本実施形態に係る超音波診断装置の他の例について説明する。本変形例では処理前のＢモードデータを蓄積しておき、スペックル／ノイズ低減処理を行うときに最新の処理前のフレームをフレーム蓄積部１０６に出力する構成である。

図９は、本変形例におけるＢモードデータの蓄積及びその蓄積したＢモードデータを使用した処理を説明するための図である。図９は説明をわかり易くするために、１次元の空間軸と時間軸との２次元として説明する。

本変形例に係る超音波診断装置は、第２の実施形態に係る超音波診断装置に点線で表わされる処理前フレーム蓄積部１０７をさらに加えた構成である。

処理前フレーム蓄積部１０７は、メモリやハードディスクなどの記憶媒体を備えている。そして、処理前フレーム蓄積部１０７は、時間の経過に対応させて順番に信号処理部００３から入力されたＢモードデータを記憶していく。これは全フレームのＢモードデータを記憶してもよいし、予め決められたフレーム数のＢモードデータを記憶する構成でもよい。予め決められたフレーム数のＢモードデータを記憶する構成では、新しい処理前のＢモードデータ９０１が入力されると、一番古い処理前のＢモードデータは削除される。

処理前フレーム蓄積部１０７は、フレーム蓄積部１０６からスペックル／ノイズ低減処理を行う処理前のＢモードデータ９０１の入力要求を受けて、フレーム蓄積部１０６に最新の処理前のＢモードデータ９０１をコピーする。

フレーム蓄積部１０６は、所定数−１フレーム前までのスペックル／ノイズ低減処理を行ったＢモードデータ９０２を記憶しておく。フレーム蓄積部１０６は、処理前フレーム蓄積部１０７からコピーされた処理前のＢモードデータ９０１と記憶している処理後のＢモードデータ９０２を多重解像度分解部１０１へ出力する。そして、フレーム蓄積部１０６は、スペックル／ノイズ低減処理部１００によりスペックル／ノイズ低減処理が行われた処理後のＢモードデータを記憶する。このとき、フレーム蓄積部１０６は、一番古い処理後Ｂモードデータを削除する。

このような構成によっても、より時間方向のスペックル／ノイズの変化を効率的に低減し、時相遅れのない処理を行うことが可能である。

００１ 超音波プローブ
００２ 送受信部
００３ 信号処理部
００４ 画像生成部
００５ 表示制御部
００６ ユーザインタフェース
００７ 表示部
００８ 入力部
００９ 統括制御部
１００ スペックル／ノイズ低減処理部
１０１ 多重解像度分解部
１０２ 特徴量算出部
１０３ フィルタリング処理部
１０４ 多重解像度合成部
１０５ 高域レベル制御部
１０６ フレーム蓄積部
１０７ 処理前フレーム蓄積部

Claims (9)

1. 被検体に対し超音波を送受信し、所定時間毎に前記被検体における特定の領域の超音波画像データを生成する画像生成手段と、
   前記超音波画像データに対して、前記領域の空間方向、又は、時間方向及び前記領域の空間方向に、階層的にｎ次（ｎは２次以上の自然数）の多重解像度分解を行い、ｋ次（１≦ｋ≦ｎ）における、全ての前記方向において空間周波数の低い画像である低域分解画像データと、少なくとも一つの前記方向において空間周波数の高い画像である高域分解画像データとを得る多重解像度分解手段と、
   算出時点の階層の１つ下の次元の出力データ又はｎ次の前記低域分解画像データからエッジの接線方向を含む当該次元における特徴量を算出する特徴量算出手段と、
   前記１つ下の次元の出力データ又はｎ次の前記低域分解画像データに対して、当該次元における前記特徴量に基づいて各前記方向毎に強調又は平滑化を行うフィルタリング処理を施し処理後低域分解画像データを取得するフィルタリング処理手段と、
   ｋ次の前記処理後低域分解画像データとｋ次の前記高域分解画像データとを多重解像度合成して（ｋ−１）次の合成画像を生成し、前記超音波画像データと同じ画像ピクセルサイズを有する処理後超音波画像データを生成するまで前記多重解像度合成を繰り返す多重解像度合成手段と、
   前記処理後超音波画像データを基に超音波画像を表示部に表示させる表示制御手段と、
   を備えることを特徴とする超音波診断装置。
2. 前記超音波画像データを生成順に順次記憶していく処理前データ蓄積手段をさらに備え、
   前記フィルタリング処理手段は、ｔ枚目の超音波画像データにおける前記フィルタリング処理として、前記記憶手段に記憶されている前記ｔ−１枚目から所定枚数前までの超音波画像データにおける同じ位置の各画素の各画素値を用いて拡散方程式を数値解析してｔ枚目の超音波画像データの前記同じ位置の各画素の画素値を算出することを特徴とする請求項１に記載の超音波診断装置。
3. 前記多重解像度合成手段が生成した前記処理後低域分解画像データを生成順に順次記憶していく処理後データ蓄積手段をさらに備え、
   前記フィルタリング処理手段は、ｔ枚目の超音波画像データにおけるフィルタリング処理として、前記記憶手段に記憶されている前記ｔ−１枚目から所定枚数前までの前記処理後低域分解画像データにおける同じ位置の各画素の各画素値を用いて拡散方程式を数値解析してｔ枚目の処理後低域分解画像データの前記同じ位置の各画素の画素値を算出することを特徴とする請求項１に記載の超音波診断装置。
4. 前記フィルタリング処理手段は、現在と過去の所定枚数の前記超音波画像データを用いて前記空間方向に対する前記フィルタリング処理を行うことを特徴とする請求項１乃至請求項３のいずれか一つに記載の超音波診断装置。
5. 前記フィルタリング処理手段は、
   非線形拡散方程式で表わされるフィルタを用いて前記フィルタリング処理を行うものであり、
   該非線形拡散方程式の数値解析を行う場合に、
   空間方向の数値微分として、
   ∂Ｉ_ｉ／∂ｘ＝｛Ｉ_ｉ＋１−Ｉ_ｉ−１｝／２
   時間方向の数値微分として、
   ∂Ｉ_ｔ／∂ｔ＝｛３Ｉ_ｔ−４Ｉ_ｔ−１−Ｉ_ｔ−２｝／２
   （Ｉ_ｉ：ｘがｉの時の画素値、Ｉ_ｔ：時刻ｔにおける画素値）
   を用いることを特徴とする請求項１乃至請求項３のいずれか一つに記載の超音波診断装置。
6. 前記特徴量算出手段は、各画素の構造テンソルを算出し、前記構造テンソルを用いて固有値と固有ベクトルを算出し、特徴量として、前記各固有値に基づいてエッジの大きさを得、前記固有ベクトルに基づいて前記エッジの接線方向を得ることを特徴とする請求項１乃至５のいずれか１つに記載の超音波診断装置。
7. 前記フィルタリング手段は、非線形異方性拡散フィルタであることを特徴とする請求項１乃至請求項６のいずれか一つに記載の超音波診断装置。
8. 前記高域分解画像データに対し前記特徴量に基づいて前記エッジを強調するフィルタリング処理を行い、該フィルタリング処理が行われた高域分解画像データを前記多重解像度合成手段に出力する高域レベル制御手段をさらに備えたことを特徴とする請求項１乃至請求項７のいずれか一つに記載の超音波診断装置。
9. 被検体に対して超音波を送受信して、所定時間毎に前記被検体における同一領域の超音波画像データを生成する画像生成手段と、
   複数の前記超音波画像データを基に、時間方向及び空間方向に多重解像度分解して、１次については空間周波数に応じて低域分解画像及び高域分解画像に分離し、以後、（ｋ―１）次に分離された低域分解画像をｋ次の低域分解画像及び該ｋ次の高域分解画像に分離することを、ｋ＝２からｎまで行う多重解像度分解手段と、
   ｎ次においては、該ｎ次の低域分解画像を基にした一方のデータと及び該ｎ次の高域分解画像を基にした他方のデータとを多重解像度合成した該ｎ次の合成画像を生成し、以後、一方のデータをｋ次の合成画像とし、他方のデータをｋ―１次の高域分解画像として合成画像を生成する前記多重解像度分解の逆変換を、ｋ＝ｎから２まで行う多重分解画像合成手段と、
   前記ｎ次の場合は該ｎ次の多重解像度分解手段と該ｎ次の多重分解画像合成手段との間にあって該ｎ次の低域分解画像を受けて画像におけるエッジの接線方向を含む該ｎ次の特徴量を算出し、前記ｋ―１次の場合は該ｋ―１次と該ｋ次との多重分解画像合成手段の間にあって該ｋ次の合成画像を受けて画像におけるエッジの接線方向を含む該ｋ−１時の特徴量を算出することをｋ＝ｎから２まで行う特徴量算出手段と、
   該特徴量算出手段から前記ｎ次から１次までの特徴量を受けて、それぞれについて前記エッジの接線方向の強調又は平滑化を行い、各次における前記一方のデータとして前記多重分解画像合成手段へ送るフィルタ処理手段と、
   前記多重解像度分解手段から前記ｎ次から１次までの高域分解画像を受けて、対応する各次の特徴量で重み付けして前記他方のデータとして前記多重分解画像合成手段へ送る高域レベル制御部と、を備え、
   前記多重分解画像合成手段の１次の出力として、元の超音波画像データと同じ画像ピクセルサイズの合成画像を得る、
   ことを特徴とする超音波診断装置。