WO2009081599A1

WO2009081599A1 - 画像処理装置、画像処理プログラム、記憶媒体及び超音波診断装置

Info

Publication number: WO2009081599A1
Application number: PCT/JP2008/055361
Authority: WO
Inventors: Yasunari Yokota; Fumio Nogata; Yoko Kawamura
Original assignee: Gifu University NUC
Current assignee: Gifu University NUC
Priority date: 2007-12-20
Filing date: 2008-03-24
Publication date: 2009-07-02
Anticipated expiration: 2010-06-20
Also published as: US8187188B2; EP2226013A4; US20110112402A1; AU2008341769A1; JP2009148395A; EP2226013A1; JP4879872B2; AU2008341769B2

Abstract

　画像処理装置１１のコンピュータ１２は、連続したフレームの超音波Ｂモード画像を取得して、所定のフレームの超音波Ｂモード画像中に含まれる頸動脈短軸断面において頸動脈壁を推定する。コンピュータ１２はまた、推定された頸動脈壁及びその周辺組織を含んだ超音波Ｂモード画像をテンプレート画像として用い、テンプレート画像を変形させることにより作成される変形後のテンプレート画像と前記取得した各フレームの超音波Ｂモード画像との誤差が最小になるように、頸動脈壁の径の大きさを推定することにより、頸動脈壁の時間変化を取得する。

Description

画像処理装置、画像処理プログラム、記憶媒体及び超音波診断装置

　本発明は、画像処理装置、画像処理プログラム、記憶媒体及び超音波診断装置に関し、特に、頸動脈の超音波画像のための画像処理装置及び画像処理プログラム、並びにそれらに関連する記憶媒体及び超音波診断装置に関する。

　動脈硬化症は、生活習慣の不摂生が危険因子となる生活習慣病の一つであり、自覚症状がほとんどないまま若年層から進行し、放っておくと脳梗塞、心筋梗塞、狭心症などに発展する可能性がある。動脈硬化症は治癒が困難であるため、早期に動脈硬化の兆候をとらえて患者に生活習慣の改善を促し、病気の進行を遅らせることが重要である。

　動脈硬化症の検査には様々な診断方法があり、なかでも超音波エコー検査は、非侵襲に血管を直接観察できることから、現在の臨床診断において必要不可欠な検査として定着しつつある。動脈硬化症診断のための超音波エコー検査では、血管内部に付着したプラークと呼ばれる脂肪分の有無や内膜中膜複合体厚（intima-media complex thickness（ＩＭＴ））を評価する。しかし、動脈硬化初期段階ではプラークもなく、ＩＭＴも薄くて計測が困難であるため、こうした形態学的観察では、初期の兆候をとらえることは難しい。

　これに対し、動脈硬化は動脈が硬化して弾力性を失うことであるから、心拍動に伴う動脈の動きを評価すればよいという考え方がある。プラーク、ＩＭＴなどの形態学的計測とは異なり、動脈の動きの定量的計測には、計算機の援用が不可欠である。

　特許文献１では、心拍動に伴う頸動脈径の変動を測定することにより、動脈硬化の指標を得られることが記載されている。
　また、本件出願人は特願２００７－２３９１９８号において、超音波Ｂモード長軸断面像を用いて、心拍動に伴う頸動脈径の変動を推定する方法を提案している。この方法は、頸動脈はその長軸（ｙ軸）上の位置によらず同質であって心拍動に応じて等しい変動を示すという前提条件のもとで行われる。この前提条件は、多くの場合満たされるが、局所的に頸動脈に動脈硬化が生じているようなケースでは満たされない。従って、そのようなケースでは、本件出願人が先に提案している方法を適用することはできない。なお、頸動脈の長軸とは頸動脈の延長方向に延びる軸をいい、長軸に直交する軸は頸動脈の短軸という。

　局所的に頸動脈に動脈硬化が生じているようなケースでは、超音波Ｂモード画像である頸動脈の短軸断面像を使用して、心拍動に伴う頸動脈の変動を推定することが考えられる。しかし、超音波Ｂモード画像の１画素のサイズはおおよそ０．０７～０．１１ｍｍであり、十分に高い解像度とはいえない。そのため、１画素のサイズよりも細かい精度で頸動脈の変動を捉える技術が不可欠である。

　特願２００７－２３９１９８号において本件出願人が提案している方法では、頸動脈長軸断面像の空間的形状に関する情報を利用して推定精度を改善している。しかし、頸動脈の短軸断面像を使用する場合には、以下の理由によりその空間的形状に関する情報を利用して推定精度を改善させることは難しい。１つ目の理由は、頸動脈の短軸断面形状はおおよそ円形ではあるが個人差が大きいため、統一した形状では表現しにくいことであり、２つ目の理由は、超音波Ｂモード画像中の頸動脈の占める面積が長軸断面像の場合に比べて短軸断面像の場合には小さいため、利用可能な情報が少ないことである。

　非特許文献１では、Ｂモード画像上で頸動脈壁上の１点を手動で指定したのち、指定した頸動脈壁上の点の動きを、ＲＦ信号を用いて１０μｍの精度で推定及び追跡する方法が提案されており、この方法は実用化されている。

　非特許文献１で提案されている方法のようにＲＦ信号を用いる場合には、サンプリングレートを高くして位相差追跡法を用いることにより推定精度を改善することができる。しかし、推定精度を改善するためのこのような手法は、Ｂモード画像を用いて頸動脈の動きを推定及び追跡する方法には適用できない。また、頸動脈の動きを追跡する場合、精度の観点からはＲＦ信号を用いることが有利であるが、ＲＦ信号を出力可能な超音波診断装置は一部の機種に限られていて一般的とはいえず、汎用性に欠ける。
特許第３８８２０８４号原田烈光,「エコートラッキング法による動脈硬化指標の計測」，超音波検査技術，Vol.28(3):pp.29-35,2003. ジーエッチ・グランルンド、「フーリエ　プレプロセシング　フォー　ハンド　プリント　キャラクタ　レコグニッシヨン(Fourier preprocessing for hand character recognition)」、IEEE Trans. on Computers, Vol. C-21, pp.195-201,1972. クリス・エル・デ・コルテ他、「キャラクタリゼーション　オブ　プラーク　コンポーネンツ　ウイズ　イントラバスキュラ　ウルトラサウンド　エラストグラフィ　イン　ヒューマン　フェモラル　アンド　コロナリイ　アーテェリーズ　イン　ヴィトロ(Characterization of Plaque Components With Intravascular Ultrasound Elastgraphy in Human Femoral and Coronary Arteries In Vitro)」、ラジ・エー・バルデュシング，クリス・エル・デ・コルテ他，「コンパリソン　オブ　ファイナイト　エレメンツ　モデル　エラストグラムス　アイヴィユーエス　エラストグラムス　アクワイヤド　フロム　ファントム　アンド　アーテェリーズ(Comparison of Finite Elements Model Elastograms and IVUS Elastograms acquired from Phantoms and Arteries)」、IEEE Symposium, pp.1873-1875,2002 林香苗、「日本人並びに日本産医学実験動物の解剖及び生理学係数」、第２版（訂正復刻版）、文京書院、1965,pp9

　従って、本発明の目的は、局所的に頸動脈に動脈硬化が生じている場合にも、超音波Ｂモード画像を使用して十分な精度で心拍動に伴う頸動脈の変動を追跡することができる画像処理装置及び画像処理プログラムを提供することにある。又、本発明の他の目的は、そのような画像処理プログラムを記憶した記憶媒体及びそのような画像処理装置に接続して使用される超音波診断装置を提供することにある。

　上記の目的を達成するために、本発明の第１の態様では、連続したフレームの超音波Ｂモード画像を取得して、所定のフレームの超音波Ｂモード画像中に含まれる頸動脈短軸断面において頸動脈壁を推定する頸動脈壁推定手段と、推定された頸動脈壁及びその周辺組織を含んだ超音波Ｂモード画像をテンプレート画像として用い、テンプレート画像を変形させることにより作成される変形後のテンプレート画像と前記取得した各フレームの超音波Ｂモード画像との誤差が最小になるように、頸動脈壁の径の大きさを推定することにより、頸動脈壁の時間変化を取得する心拍変動追跡手段とを備える画像処理装置が提供される。

　好ましくは、頸動脈壁推定手段は、頸動脈である可能性の高さを表す評価関数が最大になるように、頸動脈の位置及び形状を決める頸動脈パラメータの値を決定することにより、頸動脈壁を推定する。

　好ましくは、心拍変動追跡手段は、フレームごとに超音波Ｂモード画像上の頸動脈の中心にずれがある場合には、そのずれも推定する。
　本発明の第２の態様では、コンピュータを、連続したフレームの超音波Ｂモード画像を取得して、所定のフレームの超音波Ｂモード画像中に含まれる頸動脈短軸断面において頸動脈壁を推定する頸動脈壁推定手段、及び推定された頸動脈壁及びその周辺組織を含んだ超音波Ｂモード画像をテンプレート画像として用い、テンプレート画像を変形させることにより作成される変形後のテンプレート画像と前記取得した各フレームの超音波Ｂモード画像との誤差が最小になるように、頸動脈壁の径の大きさを推定することにより、頸動脈壁の時間変化を取得する心拍変動追跡手段として機能させる画像処理プログラムが提供される。

　好ましくは、心拍変動追跡手段は、フレームごとに超音波Ｂモード画像上の頸動脈の中心にずれがある場合には、そのずれも推定する。
　本発明の第３の態様では、上記第２の態様に係る画像処理プログラムを記憶した記憶媒体が提供される。

　本発明の第４の態様では、上記第１の態様に係る画像処理装置に接続される超音波診断装置が提供される。超音波診断装置は、探触子による超音波の送受信により得られるエコー信号を基に超音波Ｂモード画像を生成し、生成した超音波Ｂモード画像を画像処理装置に出力する。

本発明の一実施形態の画像処理装置の概略図。図１の画像処理装置の概略電気ブロック図。頸動脈短軸断面を含む超音波Ｂモード画像を示す図。二層円筒モデルの説明図。頸動脈短軸断面での座標系を示す図。距離ｒと変位量d(r; P,r₁,r₂,E₁,E₂)の関係を示すグラフ。基準状態からの血圧の変化量Ｐが－10mmHg、0mmHg、20mmHg、40mmHgであるときの距離ｒと変位量d(r; P,r₁,r₂,E₁,E₂)の関係を示すグラフ。重み関数ｗの説明図。心拍変動追跡法の概念図。図１の画像処理装置のコンピュータが実行する画像処理プログラムのフローチャート。コンピュータのＣＰＵにより推定された頸動脈壁内面の輪郭線Ｑの説明図。変化率ｄ_１／ｒ_１と２乗誤差Ｊの関係を示すグラフ。図１３（ａ）～図１３（ｆ）は２０歳の健常男性の頸動脈の心拍変動を推定した結果を示すグラフ、図１３（ｇ）～図１３（ｌ）は６０歳の健常男性の頸動脈の心拍変動を推定した結果を示すグラフ。

　以下、本発明の一実施形態を図１～図１３（ｌ）を参照して説明する。
　図１に示すように、本実施形態の画像処理装置１１は、頸動脈壁推定手段（頸動脈壁推定部）及び心拍変動追跡手段（心拍変動追跡部）として機能するコンピュータ１２、ディスプレイ１３、プリンタ１４、キーボード１５、及び超音波診断装置１６を備えている。

　超音波診断装置１６は、探触子２１、表示部２２、及び図示しない超音波発生源を備えている。超音波診断装置１６の探触子２１は、生体内にパルス波を送波して頸動脈２３からの反射波（エコー信号）を受波する先端面を有する。探触子２１の先端面は長四角形状に形成されている。超音波診断装置１６は、前記エコー信号に基づいてＢモードで取得される、頸動脈短軸断面を含むエコー動画像（すなわち、複数の超音波Ｂモード画像）を生成する。超音波診断装置１６の表示部２２には、被検者の頸動脈２３の拡張及び収縮などの変形が動画像で表示される。図５は、頸動脈短軸断面での座標系を示している。

　超音波Ｂモード画像は、超音波診断装置１６のビデオ出力から出力され、図示しないＤＶコンバータにより所定の解像度（ｍ×ｎ画素）及び所定の毎秒フレーム数でビデオキャプチャされた後、ＩＥＥＥ１３９４信号としてコンピュータ１２に入力される。本実施形態では、所定の毎秒フレーム数として３０フレーム／秒を選定しているが、この値は限定されるものではなく、心拍動に伴う頸動脈の変動を明確に観測できる値であればよい。

　ディスプレイ１３で表示される超音波画像領域の画素数は、水平方向に関してはＮｘ、垂直方向に関してはＮｙである。本実施形態では、Ｎｘ＝５４０（pixel）及びＮｙ＝４２０（pixel）を選定しているが、この数値は限定されるものではない。

　超音波診断装置１６により取得されたエコー動画像は、コンピュータ１２に入力された後、コンピュータ１２が備える記憶装置４４（図２参照）に格納される。記憶装置４４は、例えばハードディスク又は半導体記憶装置からなり、各種の情報の読み出し及び書き込みが可能である。

　図２に示すように、コンピュータ１２は、中央処理装置（ＣＰＵ）４１、ＲＯＭ４２、及びＲＡＭ４３を備え、記憶媒体としてのＲＯＭ４２に格納された画像処理プログラムを実行する。ＲＡＭ４３は、前記プログラムを実行する際の作業用メモリである。

　（実施形態の作用）
　次に、上記のように構成された画像処理装置１１において、ＣＰＵ４１が実行する画像処理プログラムの処理を説明するが、プログラム処理の説明の前にまず、「頸動脈及び周辺組織の応力及び歪み特性」、「二層円筒モデル」及び「頸動脈壁の心拍変動追跡法の概略」について説明する。

　（頸動脈及び周辺組織の応力及び歪み特性）
　頸動脈は心拍動により拡張及び収縮を繰り返す。頸動脈の周辺にある脂肪などの周辺組織は、頸動脈の拡張及び収縮に伴って伸縮する。一般にこうした物質の応力と歪みの関係は材料力学の見地から考察される。

　図３は、頸動脈短軸断面を含む超音波Ｂモード画像（以下、単に画像ともいう）を示す。図３中のほぼ中央に円形をした頸動脈５０が示されている。超音波Ｂモード画像においては、図３中でハッチングされている頸動脈内腔５０ａは黒く描かれ、頸動脈壁５０ｂは白く描かれる。頸動脈は周囲を脂肪などの周辺組織によって囲まれている。

　（二層円筒モデル）
　頸動脈及び周辺組織を、図４に示すような一様な二層円筒（材料力学的には二層厚肉円筒ともいう）でモデル化することができる。頸動脈壁の内半径をｒ_１で表し、外半径をｒ_２で表す。周辺組織は、頸動脈の中心からｒ_２だけ離れた位置から無限遠まで存在するものと仮定する。頸動脈壁の弾性係数をＥ_１で表し、周辺組織の弾性係数をＥ_２で表す。頸動脈壁のポアソン比ｖ_１及び周辺組織のポアソン比ｖ_２はいずれも０．５とみなす。頸動脈壁の内半径がｒ_１、外半径がｒ_２で表される特定の基準状態から頸動脈の内圧（血圧）に変化があったときの頸動脈及び周辺組織の径方向に関する変位量ｄ（ｒ）は、以下の二層円筒モデルの理論式（１）で表すことができる。式（１）中、Ｐは基準状態からの内圧の変化量を表す。

　頸動脈の中心からｒの距離にある組織は、内圧の変化に応じて、頸動脈の中心からｒ＋ｄ（ｒ）の距離の位置に移動する。内圧の変化量Ｐが正の場合には変位量ｄは正であって、つまり拡張を表し、内圧の変化量Ｐが負の場合には変位量ｄは負であって、すなわち収縮を表す。

　式（１）にｒ＝ｒ_１を代入すれば、頸動脈壁内面の変位量ｄ_１は以下の式（２）で表される。

　ここで、生理学的見地から、頸動脈の弾性係数Ｅ_１と周辺組織Ｅ_２の弾性係数は、それほど大きく違わないことが分かっている（非特許文献３，４参照）。また、頸動脈壁の厚みＷ（Ｗ＝ｒ_２－ｒ_１）は、２０歳代及び３０歳代で約0.73mm、４０歳代で約0.79mm、５０歳代で約0.98mm、６０歳代で約1.00mm、７０歳代で約1.04mm程度であることが分かっており（非特許文献５参照）、これは頸動脈の直径（6-12mm）や周辺組織の厚みに比べて十分に薄いといえる。以上のことを考慮すれば、上式（１）は以下の式（３）で近似することができる。

　上式（３）は、心収縮、すなわち血圧の上昇により頸動脈及び周辺組織が拡張するとき、その拡張の度合い、つまり変位量は、頸動脈の中心からの距離ｒに反比例することを意味している。一例として、Ｅ_１＝200kPa、Ｅ_２＝100kPa、ｒ_１＝4.0mm、ｒ_２＝4.8mm、Ｐ＝40mmHg（5.3kPa）のときの頸動脈の中心からの距離ｒと変位量ｄ（r; P, r₁, r₂, E₁, E₂）の関係を図６に示す。また、基準状態からの内圧の変化量Ｐだけを－10mmHg、0mmHg、20mmHg、40mmHgと変化させたときの距離ｒと変位量ｄ（r; P, r₁, r₂, E₁, E₂）の関係を図７に示す。図６及び図７に示すように、距離ｒと変位量ｄとの関係は、内圧の変化量Ｐに関わりなく、反比例で表される。

　（頸動脈壁の心拍変動追跡法の概略）
　上式（３）は変位量ｄが距離ｒに反比例することを表しているため、図６又は図７に示すグラフ中の１点、すなわち距離ｒの値とそれに対応する変位量ｄの値、例えば頸動脈壁の内半径ｒ_１の値と頸動脈壁内面の変位量ｄ_１の値が与えられれば、周辺組織も含めて頸動脈の中心からｒの距離にある組織の変位量は完全に決まる。

　頸動脈壁の内半径ｒ_１が分かる超音波Ｂモード画像Img１（テンプレート画像）があったとき、その画像を式（３）に従って変形させれば、その他の任意の血圧時における超音波Ｂモード画像Img２（変形後のテンプレート画像）を模擬的に作成できる。その際、拡張或いは収縮の度合いを制御するパラメータはｄ_１のみである。テンプレート画像から変形後のテンプレート画像を作成する際に使用されるパラメータｄ_１の値を、変形後のテンプレート画像と実際に計測して取得される画像Img２との誤差が最小になるように決定すれば、その決定されたｄ_１の値は、画像Img１の撮影時から画像Img２の撮影時までの時間における頸動脈壁内面の変位量ｄ_１の推定値となる。すなわち、画像Img１では頸動脈の中心からｒ_１の距離にあった頸動脈壁内面の位置が、画像Img２では頸動脈の中心からｒ_１＋ｄ_１の距離にまで変位したと推定することができる。換言すれば、画像Img１ではｒ_１であった頸動脈壁の内半径の大きさが、画像Img２ではｒ_１＋ｄ_１に変化したと推定することができる。実際に計測して取得されたすべての超音波Ｂモード画像について同様のことを行えば、頸動脈壁の内径の時間変化を追跡することが可能である。

　この方法のポイントは、頸動脈壁の心拍変動を追跡するにあたり、頸動脈の周辺組織が心拍動により変位することを利用している点である。超音波Ｂモード画像中の頸動脈壁の占める面積に比べて周辺組織が占める面積ははるかに大きいため、周辺組織が心拍動により変位することを利用することにより、頸動脈壁の心拍変動の追跡精度を大幅に改善することができる。この追跡精度は、超音波Ｂモード画像中の周辺組織の占める面積が大きいほど、より一層の改善が可能である。

　（画像処理プログラムの処理）
　次に、画像処理装置１１のＣＰＵ４１が実行する画像処理プログラムの処理を説明する。ＣＰＵ４１は、図１０のステップＳ１０～Ｓ４０の処理を行うことにより、頸動脈の心拍変動追跡を行う。心拍変動追跡法の概念図を図９に示す。

　ステップＳ１０では、ＣＰＵ４１は、頸動脈短軸断面を含む連続したフレームの超音波Ｂモード画像を計測及び取得する。ＣＰＵ４１により取得された超音波Ｂモード画像は、Ｉ（ｘ，ｙ，ｔ），ｘ＝１，．．．，Ｎｘ，ｙ＝１，．．．，Ｎｙ，ｔ＝１，．．．，Ｎｔと表される。ただし、ｘ，ｙ，ｔはそれぞれ、水平方向の画素番号、垂直方向の画素番号、フレーム番号を表し、各超音波Ｂモード画像上の左上に位置する画素に対してｘ＝１，ｙ＝１が与えられる。また、Ｎｘ，Ｎｙ，Ｎｔはそれぞれ、水平方向の画素総数、垂直方向の画素総数、フレーム総数を表す。１枚の超音波Ｂモード画像に着目する場合は、その画像は、Ｉ（ｘ，ｙ）、ｘ＝１，．．．，Ｎｘ，ｙ＝１，．．．，Ｎｙと表される。

　ステップＳ２０では、ＣＰＵ４１は、連続したフレームの中から適当な１つのフレーム（所定フレーム）を選択し、選択したフレームの超音波Ｂモード画像Img1について、頸動脈壁内面の位置及び形状、より具体的には頸動脈壁の内半径ｒ_１（θ）を推定する（図９参照）。内半径ｒ_１（θ）は、方位θ（０≦θ≦２π）を変数にもつ頸動脈壁の内半径ｒ_１に関する関数を表す。ステップＳ２０における処理の詳細は後述する。

　ステップＳ３０では、ＣＰＵ４１は、所定フレームの超音波Ｂモード画像を式（３）に従って変形させる。このとき、変形後の画像とその他のフレームの超音波Ｂモード画像との誤差が最小になるようにして、変位量ｄ_１の推定値を得る。全てのフレームの超音波Ｂモード画像について同様のことを行うことにより、変位量ｄ_１の時間変化ｄ_１（ｔ）は得られる。ステップＳ３０における処理の詳細は後述する。

　ステップＳ４０では、頸動脈壁の内半径の時間変化をｒ_１（θ）＋ｄ_１（ｔ）としてディスプレイ１３やプリンタ１４などの出力装置に出力する。
　（ステップＳ２０における処理の詳細）
　前記ステップＳ２０においては、ＣＰＵ４１は、連続したフレームの中から適当な一つのフレームを選択し、選択したフレームの超音波Ｂモード画像Ｉmg1中に含まれる頸動脈短軸断面において頸動脈壁の内半径ｒ_１(θ)を推定する。いずれのフレームを選択するかは、特に限定されるものではなく、無作為に選択されてもよい。あるいは、最初のフレーム、頸動脈が最も鮮明に映っているフレーム、頸動脈が収縮期にあるフレーム、あるいは頸動脈が拡張期にあるフレームが選択されてもよい。

　頸動脈壁内面の輪郭線は、真円もしくは多少歪んだ円である。頸動脈壁内面の輪郭線はは、閉曲線を表現するためのフーリェ記述子でモデル化される。このフーリェ記述子については非特許文献２に記載されている。

　フーリェ記述子を利用して頸動脈壁内面の輪郭線を表現するにあたっては、画像上の座標ｘ，ｙが複素数ｘ＋ｉｙで表される。ただし、ｉは虚数単位を表す。フーリェ記述子を利用することによって、複素平面上の閉曲線ｆ（ｓ）は、媒介変数ｓを用いて次式（４）で表現される。

　上式（４）中、ｃ_ｋは閉曲線の位置や形状を決めるパラメータであるフーリェ記述子を表し、Ｋは任意の整数を表す。整数Ｋの値が大きいほど、複雑な閉曲線を表現することができるようになるが、頸動脈壁内面の輪郭線の場合、Ｋ＝５程度でよい。ただし、整数Ｋの値はこれに限定されるものではない。以下、頸動脈壁内面の形状を規定するフーリェ記述子をまとめて頸動脈パラメータp＝{c_k|k=0,±1,±2,....,±K}と呼ぶことにする。

　頸動脈内腔は、閉曲線ｆ（ｓ）で囲まれる領域であり、これらの領域に位置する画素（ｘ，ｙ）の集合をＴで表す。閉曲線ｆ（ｓ）からの距離が頸動脈壁の厚みＷ以下であって閉曲線ｆ（ｓ）の外側に位置する画素（ｘ，ｙ）の集合をＢで表す。すなわち、集合Ｂの画素（ｘ，ｙ）は、頸動脈壁に対応する領域に位置する。集合Ｂは次式（５）で与えられる。

　画像中、頸動脈壁は強いエコーのために白く描かれ、頸動脈内腔はエコーが弱いため黒く描かれる。すなわち、頸動脈壁を表す集合Ｂの画素では画素値が大きく、頸動脈壁内腔を表す集合Ｔの画素では画素値は小さくなる。そこで、頸動脈の位置及び形状を決める頸動脈パラメータｐに対し、頸動脈である可能性の高さを表す、次式（６）で表される評価関数Ｅ（ｐ）を導入する。

　ただし、上式（６）において、α、βはいずれも正値の定数であり、|・|は、集合の要素数を表す。上式（６）中、右辺第１項は、頸動脈壁と想定される位置で画素値が大きいほど、つまり白いほど大きな値をとり、右辺第２項は、頸動脈内腔と想定される位置で画素値が小さいほど、つまり黒いほど大きな値をとる。上式（６）の右辺第３項において、ｇ（Ａ；Ａmin,Ａmax）は次式（７）で表される関数である。

　上式（７）において、Ａは、頸動脈内腔の面積、すなわち集合Ｔの画素の数を表し、Ａmin及びＡmaxはそれぞれ、頸動脈内腔の面積の考え得る最小値及び最大値を表す。上式（６）の右辺第３項は、頸動脈内腔の面積の考え得る範囲を外れた場合のペナルティを表す。Ａmin及びＡmaxの値は、これまでに収集された画像のサンプルに基づいて予め設定される。例えばＡmin＝２０(mm^２)，Ａmax＝１１０(mm^２)と設定することができるが、Ａmin及びＡmaxの値はこれに限定されるものではない。上式（６）中のβの値としては、コンピュータにおいて丸め誤差が生じない程度に十分な大きな正数が選択される。上式（６）中のαの値は、頸動脈壁と頸動脈内腔でのエコーの強さの比や、超音波診断装置１６の画像出力などの仕様に依存するが、予備試験によりα＝２．０程度を選択すればよいことが明らかになっている。

　上式（６）は、重み関数ｗ（ｘ，ｙ；ｐ）を以下の式（８）で定義した場合には、以下の式（９）で書き換えることができる。

　この重み関数ｗ（ｘ，ｙ；ｐ）の例を図８に示す。超音波Ｂモード画像Ｉ（ｘ，ｙ）における頸動脈の位置及び形状を決める頸動脈パラメータｐの値は、上式（６）若しくは上式（９）で与えられる評価関数Ｅ（ｐ）が最大になるように決定される。そのようなパラメータの値は解析的には求めることができないが、最急降下法などの非線形最適化法を利用して求めることができる。頸動脈に類似した断面円形の組織は超音波Ｂモード画像中に存在しないため、局所解に陥ることは少ない。

　ＣＰＵ４１は、得られた頸動脈壁内面の輪郭線、すなわち閉曲線ｆ（ｓ）（ただし、０≦ｓ≦１）から頸動脈壁の内半径ｒ_１（θ）を以下のようにして求める。まず、頸動脈の中心座標は、０次フーリェ記述子c₀を用いて（Ｒｅ（c₀），Ｉｍ（c₀））で与えられる。ただし、Ｒｅ（・）及びＩｍ（・）はそれぞれ、「・」の実部及び虚部を表す。頸動脈の中心から見た輪郭線は、ｆ（ｓ）－c₀（ただし、０≦ｓ≦１）で表される。そのため、ある媒介変数ｓの値での動径ｒ（ｓ）及び偏角θ（ｓ）はそれぞれ、以下の式（１０）及び式（１１）で表される。

　媒介変数ｓを外し、ｒ（ｓ）をθ（ｓ）の関数と見たものが頸動脈壁内面の形状の極座標表示ｒ_１（θ），０≦θ≦２πとなる。ＣＰＵ４１により推定された頸動脈壁内面の輪郭線Ｑ、すなわち内半径ｒ_１（θ）の一例を図１１に示す。

　（ステップＳ３０の詳細）
　上記のように推定された頸動脈壁内面及び周辺組織を適当に含む領域がＣＰＵ４１によって所定のフレームの超音波Ｂモード画像Img1から切り出される。すなわち、前記ステップＳ２０で推定された頸動脈の中心を原点として超音波Ｂモード画像Img1を極座標表示したものがテンプレート画像template(r,θ)として準備される。

　又、ＣＰＵ４１は、任意のフレーム（ｔ番目のフレーム）の超音波Ｂモード画像について、前記ステップＳ２０で推定された頸動脈の中心を原点として極座標表示したものを画像Ｉ（ｒ，θ，ｔ）として準備する。

　原理的には、任意のフレームの画像Ｉ（ｒ，θ，ｔ）は、テンプレート画像を適当なｄ_１の値を用いて変形、すなわち拡張或いは収縮させることにより、次式（１２）で近似することができる。

　上式（１２）による近似の２乗誤差Jは以下の式（１３）で表される。ＣＰＵ４１は、２乗誤差Jが最小となるように、以下の式（１４）に従ってｄ_１の値を決定する。

こうして決定されたｄ_１の値は、時刻ｔにおける頸動脈壁内面の変位量を表す（図９参照）。
　ところで、ｒは２乗誤差を評価するために利用される領域の半径であり、ｒの値を大きくすると、頸動脈形状の心拍変動の推定に利用できる領域が増える結果、推定精度が高くなることが期待される。しかし、むやみにｒの値を大きくすると、周辺組織の非一様性などの影響を受けて推定精度が逆に悪化するおそれがある。この点、ｒの値は、頸動脈壁の内半径ｒ_１（θ）の２～３倍程度であることが好ましい。

　ここで一例として、あるフレームの超音波Ｂモード画像における変化率ｄ_１／ｒ_１と２乗誤差Ｊの関係を図１２に示す。ただし、ｒ_１は、次式（１５）で表される頸動脈壁の内半径の平均である。

　図１２に示す一例では、変化率ｄ_１／ｒ_１＝0.08で２乗誤差Ｊが最小になる。従って、頸動脈壁内面の変位量ｄ_１の推定値として0.08ｒ_１が得られる。
　テンプレート画像として使用されるフレームの画像を除く取得したすべてのフレームの画像について同様の処理がＣＰＵ４１により行われることにより、変位量の時間変化ｄ_１（ｔ）が得られる。頸動脈壁の内半径の時間変化は、ｒ_１（θ）＋ｄ_１（ｔ）として与えられる。

　（頸動脈の中心が変動している場合）
　ＣＰＵ４１は、フレームごとに超音波Ｂモード画像上の頸動脈の中心にずれがある場合には、ｄ_１だけでなく、中心座標のずれ(x₀,y₀)も合わせて推定する。この場合、画像Ｉ（ｒ，θ，ｔ）の中心座標を(x₀,y₀)だけシフトさせた画像Ｉ’（ｒ，θ，ｔ）が用意される。また、上記の式（１３）及び式（１４）の代わりに以下の式（１６）及び式（１７）が用いられる。

　テンプレート画像template(r,θ)を伸縮させることによりテンプレート画像template(r＋ｄ（ｒ；ｒ_１（θ），ｄ_１）,θ)を作成する際、標本点（画素）の位置もまた伸縮を受けるため、変形後のテンプレート画像の標本点の位置は元のテンプレート画像のそれとは異なる。このままでは、計測画像Ｉ（ｒ，θ，ｔ）との誤差を評価することができないため、標本点を一致させるためにリ・サンプリングを行う必要がある。具体的には、２次元線形補間を利用してもよいし、あるいは、モーフィングなどで利用される様々な補間（内挿）技術、例えば２次元キュービック補間や２次元スプライン補間を利用してもよい。

　　（推定結果）
　２０歳の健常男性及び６０歳の健常男性からそれぞれ複数の頸動脈超音波Ｂモード画像を取得した。取得した画像を利用して、両男性の頸動脈の心拍変動を推定した結果、より具体的には頸動脈の内半径の時間変化ｒ_１（θ）＋ｄ_１（ｔ）（ただし、θ＝０（度）、０≦ｔ≦４（秒））を推定した結果を図１３（ａ）～図１３（ｌ）に示す。図１３（ａ）～図１３（ｆ）は２０歳健常男性の結果を示し、図１３（ｇ）～図１３（ｌ）は６０歳健常男性の結果を示す。図１３（ａ）～図１３（ｌ）において、縦軸は頸動脈の内半径（mm）を表し、横軸は時間（秒）を表す。図１３（ａ）～図１３（ｆ）では、同順で、頸動脈の心拍変動の推定を行っている位置が遠位端（脳に近い端）から近位端（心臓に近い端）に向かって５mmづつ移動している。同様に、図１３（ｇ）～図１３（ｌ）では、同順で、頸動脈の心拍変動の推定を行っている位置が遠位端（脳に近い端）から近位端（心臓に近い端）に向かって５mmづつ移動している。なお、図１３（ｇ）～図１３（ｌ）の縦軸のスケールは、図１３（ａ）～図１３（ｌ）の縦軸のスケールとは異なっている。

　図１３（ａ）～図１３（ｆ）と図１３（ｇ）～図１３（ｌ）を比較して明らかなように、２０歳男性の頸動脈の心拍変動が大きくはっきりとしているのに対し、６０歳男性の頸動脈の心拍変動は小さく、周期性も乏しい。又、２０歳男性の場合は、頸動脈のどの断面でも同じような心拍変動を示しているが、６０歳男性の場合は、図１３（ｈ）及び図１３（ｋ）に結果が示されている断面では、他の断面に比べて心拍変動が小さく、局所的に動脈硬化が進んでいることが予想される。

　本実施形態によれば以下の利点が得られる。
　（１）　本実施形態の画像処理装置１１は、コンピュータ１２が頸動脈壁推定手段として機能することにより、連続したフレームの超音波Ｂモード画像を取得して、所定のフレームの超音波Ｂモード画像中に含まれる頸動脈短軸断面において頸動脈壁を推定する。画像処理装置１１はまた、コンピュータ１２が心拍変動追跡手段として機能することにより、推定された頸動脈壁及びその周辺組織を含んだ超音波Ｂモード画像をテンプレート画像として用い、テンプレート画像を変形させることにより作成される変形後のテンプレート画像と取得した各フレームの超音波Ｂモード画像との誤差が最小になるように、頸動脈壁の径の大きさを推定することにより、頸動脈壁の時間変化を取得する。

　この結果、本実施形態によれば、局所的に頸動脈に動脈硬化が生じている場合にも、超音波Ｂモード画像を使用して十分な精度で頸動脈の心拍変動を追跡することができる。特に、超音波Ｂモード画像の１画素のサイズはおおよそ０．０７～０．１１ｍｍであり、十分に高い解像度とはいえないにも拘わらず、頸動脈だけでなくその周辺組織も考慮に加えていることにより、十分な精度で頸動脈の心拍変動を追跡することが可能であり、その意義は大きい。さらに、超音波Ｂモード画像は一般的な超音波診断装置から取得することができるので、画像処理装置１１は汎用性が高い。

　（２）　本実施形態の画像処理装置１１は、頸動脈である可能性の高さを表す評価関数が最大になるように、頸動脈の位置及び形状を決める頸動脈パラメータの値を決定することにより、頸動脈壁を推定する。この場合、頸動脈壁の推定が容易である。

　（３）　本実施形態の画像処理装置１１は、フレームごとに超音波Ｂモード画像上の頸動脈の中心にずれがある場合には、そのずれを推定する。この場合、フレームごとに超音波Ｂモード画像上の頸動脈の中心にずれがある場合でも、頸動脈壁の推定が容易である。

　（４）　本実施形態の画像処理プログラムは、連続したフレームの超音波Ｂモード画像を取得して、所定のフレームの超音波Ｂモード画像中に含まれる頸動脈短軸断面において頸動脈壁を推定するようにコンピュータ１２を機能させる。画像処理プログラムはまた、推定された頸動脈壁及びその周辺組織を含んだ超音波Ｂモード画像をテンプレート画像として用い、テンプレート画像を変形させることにより作成される変形後のテンプレート画像と取得した各フレームの超音波Ｂモード画像との誤差が最小になるように、頸動脈壁の径の大きさを推定することにより、頸動脈壁の時間変化を取得するようにコンピュータ１２を機能させる。この結果、局所的に頸動脈に動脈硬化が生じている場合にも、超音波Ｂモード画像を使用して十分な精度で心拍動に伴う頸動脈の変動を追跡することができる画像処理装置としてコンピュータを機能させることができる。

　（５）　本実施形態の画像処理プログラムは、頸動脈である可能性の高さを表す評価関数が最大になるように、頸動脈の位置及び形状を決める頸動脈パラメータの値を決定することにより、頸動脈壁を推定するようにコンピュータ１２を機能させる。この場合、頸動脈壁の推定が容易である。

　（６）　本実施形態の画像処理プログラムでは、フレームごとに超音波Ｂモード画像上の頸動脈の中心にずれがある場合には、そのずれを推定するようにコンピュータ１２を機能させる。この場合、フレームごとに超音波Ｂモード画像上の頸動脈の中心にずれがある場合でも、頸動脈の心拍変動の追跡が容易である。

　（７）　本実施形態のＲＯＭ４２は、上記の画像処理プログラムを記憶した記憶媒体として機能する。コンピュータ１２のＣＰＵ４１が画像処理プログラムを実行することにより、局所的に頸動脈に動脈硬化が生じている場合にも、超音波Ｂモード画像を使用して十分な精度で頸動脈の心拍変動を追跡することができる。

　（８）　本実施形態の超音波診断装置１６は、画像処理装置１１に接続されているため、上記（１）乃至（３）の効果を容易に実現することができる。
　前記実施形態を次のように変更して構成することもできる。

　画像処理プログラムを記憶する記憶媒体として、ＲＯＭ４２の代わりに記憶装置４４を用いてもよい。或いは、ＣＤやＤＶＤ、ＵＳＢメモリ等の記憶媒体に画像処理プログラムを格納し、格納された画像処理プログラムを当該記憶媒体からドライバ装置等を用いて読込みして実行するようにしてもよい。

　動脈硬化に関する解析の結果は、ディスプレイ１３またはプリンタ１４のいずれか一方のみから出力されるように構成してもよい。

Claims

　連続したフレームの超音波Ｂモード画像を取得して、所定のフレームの超音波Ｂモード画像中に含まれる頸動脈短軸断面において頸動脈壁を推定する頸動脈壁推定手段と、
　推定された頸動脈壁及びその周辺組織を含んだ超音波Ｂモード画像をテンプレート画像として用い、テンプレート画像を変形させることにより作成される変形後のテンプレート画像と前記取得した各フレームの超音波Ｂモード画像との誤差が最小になるように、頸動脈壁の径の大きさを推定することにより、頸動脈壁の時間変化を取得する心拍変動追跡手段と
を備えることを特徴とする画像処理装置。
　前記頸動脈壁推定手段は、頸動脈である可能性の高さを表す評価関数が最大になるように、頸動脈の位置及び形状を決める頸動脈パラメータの値を決定することにより、頸動脈壁を推定することを特徴とする請求項１に記載の画像処理装置。
　前記心拍変動追跡手段は、フレームごとに超音波Ｂモード画像上の頸動脈の中心にずれがある場合には、そのずれも推定することを特徴とする請求項１又は請求項２に記載の画像処理装置。
　コンピュータを、
　連続したフレームの超音波Ｂモード画像を取得して、所定のフレームの超音波Ｂモード画像中に含まれる頸動脈短軸断面において頸動脈壁を推定する頸動脈壁推定手段、及び
　推定された頸動脈壁及びその周辺組織を含んだ超音波Ｂモード画像をテンプレート画像として用い、テンプレート画像を変形させることにより作成される変形後のテンプレート画像と前記取得した各フレームの超音波Ｂモード画像との誤差が最小になるように、頸動脈壁の径の大きさを推定することにより、頸動脈壁の時間変化を取得する心拍変動追跡手段
として機能させることを特徴とする画像処理プログラム。
　前記頸動脈壁推定手段は、頸動脈である可能性の高さを表す評価関数が最大になるように、頸動脈の位置及び形状を決める頸動脈パラメータの値を決定することにより、頸動脈壁を推定することを特徴とする請求項４に記載の画像処理プログラム。
　前記心拍変動追跡手段は、
　フレームごとに超音波Ｂモード画像上の頸動脈の中心にずれがある場合には、そのずれも推定することを特徴とする請求項４又は請求項５に記載の画像処理プログラム。
　請求項４乃至請求項６のいずれか１項に記載の画像処理プログラムを記憶した記憶媒体。
　請求項１乃至請求項３のいずれか１項に記載の画像処理装置に接続され、探触子による超音波の送受信により得られるエコー信号を基に超音波Ｂモード画像を生成し、生成した超音波Ｂモード画像を画像処理装置に出力することを特徴とする超音波診断装置。